diff --git a/lib/msprime.c b/lib/msprime.c index fc84f3cc3..b5d6c000b 100644 --- a/lib/msprime.c +++ b/lib/msprime.c @@ -4225,6 +4225,26 @@ msp_migration_event(msp_t *self, population_id_t source_pop, population_id_t des return ret; } +/* Migration event in an arbitrary label + */ +static int MSP_WARN_UNUSED +msp_migration_event_in_background( + msp_t *self, population_id_t source_pop, population_id_t dest_pop, label_id_t label) +{ + int ret = 0; + uint32_t j; + avl_node_t *node; + avl_tree_t *source = &self->populations[source_pop].ancestors[label]; + size_t index = ((size_t) source_pop) * self->num_populations + (size_t) dest_pop; + + self->num_migration_events[index]++; + j = (uint32_t) gsl_rng_uniform_int(self->rng, avl_count(source)); + node = avl_at(source, j); + tsk_bug_assert(node != NULL); + ret = msp_move_individual(self, node, source, dest_pop, label); + return ret; +} + static int MSP_WARN_UNUSED msp_reset_memory_state(msp_t *self) { @@ -5803,6 +5823,44 @@ msp_sweep_initialise(msp_t *self, double switch_proba) return ret; } +/* Set up the intial populations for a sweep by moving individuals from + * label 0 to label 1 with the specified probability for multiple demes. + Should be merged with ms_sweep_initialise at some point. + */ +static int +msp_sweep_reverse_initialise(msp_t *self, double *switch_proba) +{ + int ret = 0; + uint32_t j; + avl_node_t *node, *next; + avl_tree_t *pop; + + /* We only support two labels for now */ + if (self->num_labels != 2) { + ret = MSP_ERR_UNSUPPORTED_OPERATION; + goto out; + } + + /* Move ancestors to new labels. */ + for (j = 0; j < self->num_populations; j++) { + tsk_bug_assert(avl_count(&self->populations[j].ancestors[1]) == 0); + pop = &self->populations[j].ancestors[0]; + node = pop->head; + while (node != NULL) { + next = node->next; + if (gsl_rng_uniform(self->rng) < switch_proba[j]) { + ret = msp_move_individual(self, node, pop, (population_id_t) j, 1); + if (ret != 0) { + goto out; + } + } + node = next; + } + } +out: + return ret; +} + /* Finalise the sweep by moving all lineages back to label 0. */ static int @@ -5882,6 +5940,230 @@ msp_sweep_recombination_event( return ret; } +/* Recombination event during a sweep within a specified background with multiple demes. +Should be merged with msp_sweep_recombination_event at some point. + */ +static int +msp_sweep_reverse_recombination_event(msp_t *self, label_id_t label, double sweep_locus, + double *mutant_population_frequency) +{ + int ret = 0; + lineage_t *left_lin, *right_lin; + double r; + + ret = msp_recombination_event(self, label, &left_lin, &right_lin); + if (ret != 0) { + goto out; + } + + tsk_bug_assert(left_lin != NULL); + tsk_bug_assert(right_lin != NULL); + + /* NOTE: we can look at rhs->left when we compare to the sweep site. */ + r = gsl_rng_uniform(self->rng); + if (label == 0) { + if (sweep_locus < right_lin->head->left) { + if (r < mutant_population_frequency[right_lin->population]) { + ret = msp_change_label(self, right_lin, 1); + if (ret != 0) { + goto out; + } + } + } else { + if (r < mutant_population_frequency[left_lin->population]) { + ret = msp_change_label(self, left_lin, 1); + if (ret != 0) { + goto out; + } + } + } + } + if (label == 1) { + if (sweep_locus < right_lin->head->left) { + if (r < 1.0 - mutant_population_frequency[right_lin->population]) { + ret = msp_change_label(self, right_lin, 0); + if (ret != 0) { + goto out; + } + } + } else { + if (r < 1.0 - mutant_population_frequency[left_lin->population]) { + ret = msp_change_label(self, left_lin, 0); + if (ret != 0) { + goto out; + } + } + } + } +out: + return ret; +} + +/* Reads the trajectory of a sweep when fed a file. +The file contains (in order) +1. Number of steps +2. Number of demes +3. Pop size +4. Migration rate +5. State after finishing the backward sim +6. state at the start of the backward sim +7. Time, type, and start and end locations for each event +No. 7 is currently read backward in time as well + */ +static int +genic_selection_read_trajectory(const char *filename, size_t *num_steps_ret, + int *num_demes_ret, int *tot_pop_ret, double *migration_rate_ret, + int **final_mut_pop_ret, int **mut_pop_ret, double **allele_frequency_mut_ret, + double **t_of_forward_ev_ret, int **ev_type_ret, int **start_deme_ret, + int **end_deme_ret) +{ + int ret = 0; + size_t num_steps; + int num_demes; + double *allele_frequency_mut = NULL; + double *t_of_forward_ev = NULL; + size_t j = 0; + int i = 0; + double migration_rate = 0.0; + int tot_pop = 0; + int *mut_pop = NULL, *final_mut_pop = NULL; + int *ev_type = NULL, *start_deme = NULL, *end_deme = NULL; + FILE *file = fopen(filename, "r"); + + fread(&num_steps, sizeof(size_t), 1, file); + fread(&num_demes, sizeof(int), 1, file); + fread(&tot_pop, sizeof(int), 1, file); + fread(&migration_rate, sizeof(double), 1, file); + + mut_pop = (int *) malloc(sizeof(*mut_pop) * (size_t) num_demes); + allele_frequency_mut + = (double *) malloc(sizeof(*allele_frequency_mut) * (size_t) num_demes); + final_mut_pop = (int *) malloc(sizeof(*final_mut_pop) * (size_t) num_demes); + + for (i = 0; i < num_demes; i++) { + fread(&final_mut_pop[i], sizeof(int), 1, file); + } + + for (i = 0; i < num_demes; i++) { + fread(&mut_pop[i], sizeof(int), 1, file); + allele_frequency_mut[i] = 1.0 * mut_pop[i] / tot_pop; + } + + t_of_forward_ev = (double *) malloc(sizeof(*t_of_forward_ev) * num_steps); + ev_type = (int *) malloc(sizeof(*ev_type) * num_steps); + start_deme = (int *) malloc(sizeof(*start_deme) * num_steps); + end_deme = (int *) malloc(sizeof(*end_deme) * num_steps); + for (j = 1; j <= num_steps; j++) { + fread(&t_of_forward_ev[num_steps - j], sizeof(double), 1, file); + fread(&ev_type[num_steps - j], sizeof(int), 1, file); + fread(&start_deme[num_steps - j], sizeof(int), 1, file); + fread(&end_deme[num_steps - j], sizeof(int), 1, file); + } + + fclose(file); + *num_demes_ret = num_demes; + *num_steps_ret = num_steps; + *tot_pop_ret = tot_pop; + *migration_rate_ret = migration_rate; + *final_mut_pop_ret = final_mut_pop; + *mut_pop_ret = mut_pop; + *allele_frequency_mut_ret = allele_frequency_mut; + *t_of_forward_ev_ret = t_of_forward_ev; + *ev_type_ret = ev_type; + *start_deme_ret = start_deme; + *end_deme_ret = end_deme; + return ret; +} + +static int +sweep_forward_event(msp_t *self, double t_of_next_forward_ev, int *ev_type, + size_t *curr_step, int *start_deme, int *end_deme, double **ancestral_bg_num_ind, + int **mut_pop, double **allele_frequency_mut, int tot_pop, int num_demes) +{ + int ret = 0; + int curr_ev_type = 0, start_deme_index = 0; + double p_forward_ev, tmp_rand; + + self->time = t_of_next_forward_ev; + curr_ev_type = ev_type[*curr_step]; + start_deme_index = start_deme[*curr_step]; + + if (start_deme_index == end_deme[*curr_step]) { + if (curr_ev_type == 0) { // mutant replaced with wildtype, mutant coalascence + p_forward_ev + = (float) ancestral_bg_num_ind[start_deme_index][1] + * (ancestral_bg_num_ind[start_deme_index][1] - 1.0) + / ((*mut_pop[start_deme_index]) * (*mut_pop[start_deme_index] - 1.0)); + tmp_rand = gsl_rng_uniform(self->rng); + if (tmp_rand <= p_forward_ev) { + ret = self->common_ancestor_event(self, start_deme_index, 1); + } + (*mut_pop[start_deme_index])--; + *allele_frequency_mut[start_deme_index] + = 1.0 * *mut_pop[start_deme_index] / tot_pop; + (*curr_step)++; + + } else if (curr_ev_type == 1) { // wildtype replaced with mutant, wt coalescence + p_forward_ev = (float) ancestral_bg_num_ind[start_deme_index][0] + * (ancestral_bg_num_ind[start_deme_index][0] - 1.0) + / ((tot_pop - *mut_pop[start_deme_index]) + * ((tot_pop - *mut_pop[start_deme_index]) - 1.0)); + tmp_rand = gsl_rng_uniform(self->rng); + if (tmp_rand <= p_forward_ev) { + ret = self->common_ancestor_event(self, start_deme_index, 0); + } + (*mut_pop[start_deme_index])++; + *allele_frequency_mut[start_deme_index] + = 1.0 * *mut_pop[start_deme_index] / tot_pop; + (*curr_step)++; + } + + } else if (num_demes > 1) { + + if (curr_ev_type == 2) { // wildtype replaced by mutant, mutant migration + // + coalascence + p_forward_ev = 1.0 * ancestral_bg_num_ind[start_deme_index][1] + / *mut_pop[start_deme_index]; + tmp_rand = gsl_rng_uniform(self->rng); + if (tmp_rand <= p_forward_ev) { + ret = msp_migration_event_in_background( + self, start_deme_index, end_deme[*curr_step], 1); + p_forward_ev = 1.0 * ancestral_bg_num_ind[end_deme[*curr_step]][1] + / *mut_pop[end_deme[*curr_step]]; + tmp_rand = gsl_rng_uniform(self->rng); + if (tmp_rand <= p_forward_ev) { + ret = self->common_ancestor_event(self, end_deme[*curr_step], 1); + } + } + (*mut_pop[start_deme[*curr_step]])--; + *allele_frequency_mut[start_deme[*curr_step]] + = 1.0 * *mut_pop[start_deme[*curr_step]] / tot_pop; + (*curr_step)++; + + } else if (curr_ev_type == 3) { // mutant replaced by wildtype, wildtype + // migration + coalascence + p_forward_ev = 1.0 * ancestral_bg_num_ind[start_deme_index][0] + / (tot_pop - *mut_pop[start_deme_index]); + tmp_rand = gsl_rng_uniform(self->rng); + if (tmp_rand <= p_forward_ev) { + ret = msp_migration_event_in_background( + self, start_deme_index, end_deme[*curr_step], 0); + p_forward_ev = 1.0 * ancestral_bg_num_ind[end_deme[*curr_step]][0] + / (tot_pop - *mut_pop[end_deme[*curr_step]]); + tmp_rand = gsl_rng_uniform(self->rng); + if (tmp_rand <= p_forward_ev) { + ret = self->common_ancestor_event(self, end_deme[*curr_step], 0); + } + } + (*mut_pop[start_deme[*curr_step]])++; + *allele_frequency_mut[start_deme[*curr_step]] + = 1.0 * *mut_pop[start_deme[*curr_step]] / tot_pop; + (*curr_step)++; + } + } + return ret; +} + static int msp_run_sweep(msp_t *self) { @@ -6064,6 +6346,304 @@ msp_run_sweep(msp_t *self) return ret; } +/* Runs a backward in time sweep after reading the forward trajectory from a file + */ +static int +msp_run_sweep_reverse(msp_t *self) +{ + int ret = 0; + simulation_model_t *model = &self->model; + size_t curr_step = 0; + size_t num_steps = 0; + double *allele_frequency_mut = NULL, *t_of_forward_ev = NULL; + int *ev_type = NULL, *start_deme = NULL, *end_deme = NULL; + double sweep_locus = model->params.sweep_reverse.position; + const char *filename = model->params.sweep_reverse.filename; + size_t j = 0; + int i = 0; + double recomb_mass = 0.0; + label_id_t label; + double rec_rates[] = { 0.0, 0.0 }; + double **ancestral_bg_num_ind = NULL; + double *p_coal_wild = NULL, *p_coal_mut = NULL, *p_mig_wild_right = NULL, + *p_mig_mut_right = NULL, *p_mig_wild_left = NULL, *p_mig_mut_left = NULL; + double p_rec_wild = 0.0, p_rec_mut = 0.0; + double tmp_rand = 0.0, p_any_ev = 0.0, p_cum_ev = 0.0, t_of_next_forward_ev = 0.0; + double t_start = 0.0, t_end = 0.0, t_of_next_ev = 0.0, t_current = 0.0; + int num_demes = 0; + double migration_rate = 0.0; + int tot_pop = 0; + int *mut_pop = NULL, *final_mut_pop = NULL; + bool no_event_yet = true; + + if (rate_map_get_total_mass(&self->gc_map) != 0.0) { // arrow?+ + /* Could be, we just haven't implemented it */ + ret = MSP_ERR_SWEEPS_GC_NOT_SUPPORTED; + goto out; + } + + ret = genic_selection_read_trajectory(filename, &num_steps, &num_demes, &tot_pop, + &migration_rate, &final_mut_pop, &mut_pop, &allele_frequency_mut, + &t_of_forward_ev, &ev_type, &start_deme, &end_deme); + + tsk_bug_assert(num_demes > 0); + tsk_bug_assert(tot_pop > 0); + + p_coal_wild = (double *) malloc(sizeof(double) * (size_t) num_demes); + p_coal_mut = (double *) malloc(sizeof(double) * (size_t) num_demes); + if (num_demes > 1) { + p_mig_wild_right = (double *) malloc(sizeof(double) * (size_t)(num_demes - 1)); + p_mig_mut_right = (double *) malloc(sizeof(double) * (size_t)(num_demes - 1)); + p_mig_wild_left = (double *) malloc(sizeof(double) * (size_t)(num_demes - 1)); + p_mig_mut_left = (double *) malloc(sizeof(double) * (size_t)(num_demes - 1)); + } + + ancestral_bg_num_ind = (double **) malloc(sizeof(double *) * (size_t) num_demes); + for (i = 0; i < num_demes; i++) { + ancestral_bg_num_ind[i] = (double *) malloc(sizeof(double) * 2); + } + + curr_step = 0; + + ret = msp_sweep_reverse_initialise(self, allele_frequency_mut); + if (ret != 0) { + goto out; + } + + /* Sets the initial time of the sweep. Adds a small increment after the final event. + * Also makes the times of the forward step consistent with msprime + */ + t_start = self->time; + t_current = t_start; + t_of_next_ev = 0.0; + t_end = t_of_forward_ev[0]; + for (j = 0; j < num_steps; j++) { /*move this step to file io*/ + t_of_forward_ev[j] = t_start + t_end - t_of_forward_ev[j] + + 1e-5; /*Add a small jitter to avoid time travel*/ + } + + while (msp_get_num_ancestors(self) > 0 && curr_step < num_steps + && self->time < (t_start + t_end)) { + /* Set num ancestral individuals & rec_rates */ + for (j = 0; j < self->num_labels; j++) { + label = (label_id_t) j; + recomb_mass = self->recomb_mass_index == NULL + ? 0 + : fenwick_get_total(&self->recomb_mass_index[label]); + for (i = 0; i < num_demes; i++) { + ancestral_bg_num_ind[i][j] + = avl_count(&self->populations[i].ancestors[label]); + /* We can get small negative rates by numerical jitter which causes + * problems in later calculations and leads to assertion trips. + * See https://github.com/tskit-dev/msprime/issues/1966 + */ + } + rec_rates[j] = TSK_MAX(0, (recomb_mass)); + } + + /* Set coal and mig rates*/ + + for (i = 0; i < num_demes; i++) { + if (ancestral_bg_num_ind[i][0] > 1) { + tsk_bug_assert(mut_pop[i] != tot_pop); + p_coal_wild[i] + = ((ancestral_bg_num_ind[i][0] * (ancestral_bg_num_ind[i][0] - 1.0)) + * 0.5) + / (tot_pop - mut_pop[i]); + } else { + p_coal_wild[i] = 0; + } + if (ancestral_bg_num_ind[i][1] > 1) { + tsk_bug_assert(mut_pop[i] != 0); + p_coal_mut[i] + = ((ancestral_bg_num_ind[i][1] * (ancestral_bg_num_ind[i][1] - 1.0)) + * 0.5) + / mut_pop[i]; + } else { + p_coal_mut[i] = 0; + } + + if (num_demes > 1) { + if (i > 0) { + p_mig_wild_left[i - 1] = migration_rate * ancestral_bg_num_ind[i][0] + * (1.0 - mut_pop[i - 1] / tot_pop); + p_mig_mut_left[i - 1] = migration_rate * ancestral_bg_num_ind[i][1] + * mut_pop[i - 1] / tot_pop; + } + if (i < num_demes - 1) { + p_mig_wild_right[i] = migration_rate * ancestral_bg_num_ind[i][0] + * (1.0 - mut_pop[i + 1] / tot_pop); + p_mig_mut_right[i] = migration_rate * ancestral_bg_num_ind[i][1] + * mut_pop[i + 1] / tot_pop; + } + } + } + + p_rec_wild = rec_rates[0]; + p_rec_mut = rec_rates[1]; + p_any_ev = 0; + + /* Find rate of any event */ + for (i = 0; i < num_demes; i++) { + p_any_ev += p_coal_wild[i]; + p_any_ev += p_coal_mut[i]; + } + + if (num_demes > 1) { + for (i = 0; i < (num_demes - 1); i++) { + p_any_ev += p_mig_wild_left[i]; + p_any_ev += p_mig_mut_left[i]; + p_any_ev += p_mig_wild_right[i]; + p_any_ev += p_mig_mut_right[i]; + } + } + + p_any_ev += p_rec_wild; + p_any_ev += p_rec_mut; + + t_current = self->time; + t_of_next_ev = gsl_ran_exponential(self->rng, 1 / p_any_ev); + tsk_bug_assert(t_of_next_ev > 0); + + t_of_next_forward_ev = t_of_forward_ev[curr_step]; + + if (t_current + t_of_next_ev >= t_of_next_forward_ev) { + + ret = sweep_forward_event(self, t_of_next_forward_ev, ev_type, &curr_step, + start_deme, end_deme, ancestral_bg_num_ind, &mut_pop, + &allele_frequency_mut, tot_pop, num_demes); + + } else { + self->time += t_of_next_ev; + p_cum_ev = 0.0; + no_event_yet = true; + tmp_rand = gsl_rng_uniform(self->rng); + + for (i = 0; i < num_demes; i++) { + p_cum_ev += p_coal_wild[i]; + if (tmp_rand < p_cum_ev / p_any_ev) { + ret = self->common_ancestor_event(self, i, 0); + no_event_yet = false; + break; + } else { + p_cum_ev += p_coal_mut[i]; + if (tmp_rand < p_cum_ev / p_any_ev) { + ret = self->common_ancestor_event(self, i, 1); + no_event_yet = false; + break; + } + } + } + + if (num_demes > 1) { + if (no_event_yet) { + for (i = 0; i < (num_demes - 1); i++) { + p_cum_ev += p_mig_wild_left[i]; + if (tmp_rand < p_cum_ev / p_any_ev) { + ret = msp_migration_event_in_background(self, i + 1, i, 0); + no_event_yet = false; + break; + } else { + p_cum_ev += p_mig_mut_left[i]; + if (tmp_rand < p_cum_ev / p_any_ev) { + ret = msp_migration_event_in_background( + self, i + 1, i, 1); + no_event_yet = false; + break; + } else { + p_cum_ev += p_mig_wild_right[i]; + if (tmp_rand < p_cum_ev / p_any_ev) { + ret = msp_migration_event_in_background( + self, i, i + 1, 0); + no_event_yet = false; + break; + } else { + p_cum_ev += p_mig_mut_right[i]; + if (tmp_rand < p_cum_ev / p_any_ev) { + ret = msp_migration_event_in_background( + self, i, i + 1, 1); + no_event_yet = false; + break; + } + } + } + } + } + } + } + + if (no_event_yet) { + p_cum_ev += p_rec_wild; + if (tmp_rand < p_cum_ev / p_any_ev) { + ret = msp_sweep_reverse_recombination_event( + self, 0, sweep_locus, allele_frequency_mut); + no_event_yet = false; + } + } + if (no_event_yet) { + p_cum_ev += p_rec_mut; + if (tmp_rand < p_cum_ev / p_any_ev) { + ret = msp_sweep_reverse_recombination_event( + self, 1, sweep_locus, allele_frequency_mut); + no_event_yet = false; + } + } + } + + if (ret != 0) { + goto out; + } + } + + /* TODO we should probably support fixed events here using + * msp_apply_fixed_events() like we do in the coalescent models. + * The point below about computing population sizes should be easily + * worked around using msp_compute_population_size(). */ + + /* Check if any demographic events should have happened during the + * event and raise an error if so. This is to keep computing population + * sizes simple */ + if (self->next_demographic_event != NULL + && self->next_demographic_event->time <= self->time) { + ret = MSP_ERR_EVENTS_DURING_SWEEP; + goto out; + } + /* Rule out sampling events for now, but there's no reason we can't + * support them. */ + if (self->next_sampling_event < self->num_sampling_events + && self->sampling_events[self->next_sampling_event].time <= self->time) { + ret = MSP_ERR_EVENTS_DURING_SWEEP; + goto out; + } + ret = msp_sweep_finalise(self); + if (ret != 0) { + goto out; + } + ret = MSP_EXIT_MODEL_COMPLETE; + +out: + msp_safe_free(p_coal_mut); + msp_safe_free(p_coal_wild); + if (num_demes > 1) { + msp_safe_free(p_mig_wild_left); + msp_safe_free(p_mig_mut_left); + msp_safe_free(p_mig_wild_right); + msp_safe_free(p_mig_mut_right); + } + for (i = 0; i < num_demes; i++) { + msp_safe_free(ancestral_bg_num_ind[i]); + } + msp_safe_free(ancestral_bg_num_ind); + msp_safe_free(mut_pop); + msp_safe_free(allele_frequency_mut); + msp_safe_free(final_mut_pop); + msp_safe_free(t_of_forward_ev); + msp_safe_free(ev_type); + msp_safe_free(start_deme); + msp_safe_free(end_deme); + return ret; +} + /* Runs the simulation backwards in time until either the sample has coalesced, * or specified maximum simulation time has been reached or the specified maximum * number of events has been reached. @@ -6093,6 +6673,8 @@ msp_run(msp_t *self, double max_time, unsigned long max_events) } else if (self->model.type == MSP_MODEL_SWEEP) { /* FIXME making sweep atomic for now as it's non-rentrant */ ret = msp_run_sweep(self); + } else if (self->model.type == MSP_MODEL_SWEEP_REVERSE) { + ret = msp_run_sweep_reverse(self); } else { ret = msp_run_coalescent(self, max_time, max_events); } @@ -6345,6 +6927,9 @@ msp_get_model_name(msp_t *self) case MSP_MODEL_SWEEP: ret = "single-sweep"; break; + case MSP_MODEL_SWEEP_REVERSE: + ret = "single-sweep, reverse only"; + break; default: ret = "BUG: bad model in simulator!"; break; @@ -8171,7 +8756,8 @@ msp_set_simulation_model(msp_t *self, int model) if (model != MSP_MODEL_HUDSON && model != MSP_MODEL_SMC && model != MSP_MODEL_SMC_PRIME && model != MSP_MODEL_SMC_K && model != MSP_MODEL_DIRAC && model != MSP_MODEL_BETA && model != MSP_MODEL_DTWF - && model != MSP_MODEL_WF_PED && model != MSP_MODEL_SWEEP) { + && model != MSP_MODEL_WF_PED && model != MSP_MODEL_SWEEP + && model != MSP_MODEL_SWEEP_REVERSE) { ret = MSP_ERR_BAD_MODEL; goto out; } @@ -8454,3 +9040,29 @@ msp_set_simulation_model_sweep_genic_selection(msp_t *self, double position, out: return ret; } + +int +msp_set_simulation_model_sweep_genic_selection_reverse( + msp_t *self, double position, const char *filename) +{ + int ret = 0; + simulation_model_t *model = &self->model; + double L = self->sequence_length; + + /* Check the inputs to make sure they make sense */ + if (position < 0 || position >= L) { + ret = MSP_ERR_BAD_SWEEP_POSITION; + goto out; + } + + ret = msp_set_simulation_model(self, MSP_MODEL_SWEEP_REVERSE); + if (ret != 0) { + goto out; + } + + model->params.sweep_reverse.position = position; + model->params.sweep_reverse.filename = filename; + +out: + return ret; +} \ No newline at end of file diff --git a/lib/msprime.h b/lib/msprime.h index c352ef539..e84799371 100644 --- a/lib/msprime.h +++ b/lib/msprime.h @@ -43,8 +43,9 @@ #define MSP_MODEL_DIRAC 4 #define MSP_MODEL_DTWF 5 #define MSP_MODEL_SWEEP 6 -#define MSP_MODEL_WF_PED 7 -#define MSP_MODEL_SMC_K 8 +#define MSP_MODEL_SWEEP_REVERSE 7 +#define MSP_MODEL_WF_PED 8 +#define MSP_MODEL_SMC_K 9 /* Exit codes from msp_run to distinguish different reasons for exiting * before coalescence. */ @@ -231,6 +232,11 @@ typedef struct _sweep_t { void (*print_state)(struct _sweep_t *self, FILE *out); } sweep_t; +typedef struct _sweep_reverse_t { + double position; + const char *filename; +} sweep_reverse_t; + typedef struct _simulation_model_t { int type; union { @@ -238,6 +244,7 @@ typedef struct _simulation_model_t { beta_coalescent_t beta_coalescent; dirac_coalescent_t dirac_coalescent; sweep_t sweep; + sweep_reverse_t sweep_reverse; } params; /* If the model allocates memory this function should be non-null. */ void (*free)(struct _simulation_model_t *model); @@ -491,6 +498,8 @@ int msp_set_simulation_model_dirac(msp_t *self, double psi, double c); int msp_set_simulation_model_beta(msp_t *self, double alpha, double truncation_point); int msp_set_simulation_model_sweep_genic_selection(msp_t *self, double position, double start_frequency, double end_frequency, double s, double dt); +int msp_set_simulation_model_sweep_genic_selection_reverse( + msp_t *self, double position, const char *filename); int msp_set_start_time(msp_t *self, double start_time); int msp_set_store_migrations(msp_t *self, bool store_migrations); diff --git a/lib/tests/test_sweeps.c b/lib/tests/test_sweeps.c index 93f107750..b070bfc3c 100644 --- a/lib/tests/test_sweeps.c +++ b/lib/tests/test_sweeps.c @@ -330,7 +330,7 @@ sweep_genic_selection_mimic_msms_single_run(unsigned long int seed) double s = 10000; double recom_rate = 0.0004; double start_frequency = 0.5 / 10000; - double end_frequency = 0.9; + double end_frequency = 1 - 0.5 / 10000; double dt = 1.0 / 400000; msp_t msp; gsl_rng *rng = safe_rng_alloc(); @@ -370,6 +370,61 @@ sweep_genic_selection_mimic_msms_single_run(unsigned long int seed) tsk_table_collection_free(&tables); } +static void +sweep_genic_selection_reverse_single_run(unsigned long int seed) +{ + + /* Try to mimic the msms parameters used in verification.py + "100 300 -t 200 -r 200 500000" + " -SF 0 0.9 -Sp 0.5 -SaA 5000 -SAA 10000 -N 10000" + */ + int ret; + const char *filename + = "/home/aalhadbhatt/Desktop/msprime/notebooks/PyNBmsprime/events.bin"; + uint32_t n = 10; + double num_loci = 500001; + double position = num_loci / 2; + size_t num_demes = 1; + double recom_rate = 0.0004; + msp_t msp; + gsl_rng *rng = safe_rng_alloc(); + tsk_table_collection_t tables; + + // Test over differnt seeds + gsl_rng_set(rng, seed); + + ret = build_sim(&msp, &tables, rng, num_loci, num_demes, NULL, n); + CU_ASSERT_EQUAL(ret, 0); + CU_ASSERT_EQUAL_FATAL(msp_set_recombination_rate(&msp, recom_rate), 0); + ret = msp_set_num_labels(&msp, 2); + CU_ASSERT_EQUAL(ret, 0); + + // To mimic the verfication.py call + msp_set_discrete_genome(&msp, 0); + msp_set_gene_conversion_rate(&msp, 0); + msp_set_gene_conversion_tract_length(&msp, 1); + msp_set_avl_node_block_size(&msp, 65536); + msp_set_node_mapping_block_size(&msp, 65536); + msp_set_segment_block_size(&msp, 65536); + + ret = msp_set_simulation_model_sweep_genic_selection_reverse( + &msp, position, filename); + CU_ASSERT_EQUAL(ret, 0); + + ret = msp_initialise(&msp); + CU_ASSERT_EQUAL(ret, 0); + + ret = msp_run(&msp, DBL_MAX, UINT32_MAX); + CU_ASSERT_EQUAL(ret, MSP_EXIT_MODEL_COMPLETE); + + msp_verify(&msp, 0); + + msp_free(&msp); + gsl_rng_free(rng); + // free(samples); + tsk_table_collection_free(&tables); +} + static void test_sweep_genic_selection_mimic_msms(void) { @@ -379,6 +434,14 @@ test_sweep_genic_selection_mimic_msms(void) } } +static void +test_sweep_genic_selection_reverse(void) +{ + /* To mimic the nrepeats = 300 parameter in msms cmdline arguments*/ + for (int i = 0; i < 300; i++) + sweep_genic_selection_reverse_single_run(i + 1); +} + int main(int argc, char **argv) { @@ -395,6 +458,7 @@ main(int argc, char **argv) test_sweep_genic_selection_time_change }, { "test_sweep_genic_selection_mimic_msms", test_sweep_genic_selection_mimic_msms }, + { "test_sweep_genic_selection_reverse", test_sweep_genic_selection_reverse }, CU_TEST_INFO_NULL, }; diff --git a/msprime/__init__.py b/msprime/__init__.py index 63b18fe4b..cc903a225 100644 --- a/msprime/__init__.py +++ b/msprime/__init__.py @@ -47,6 +47,7 @@ StandardCoalescent, SweepGenicSelection, FixedPedigree, + SweepGenicSelectionReverse, TimeUnitsMismatchWarning, NodeType, ) diff --git a/msprime/_msprimemodule.c b/msprime/_msprimemodule.c index b3de9398e..a3c794e05 100644 --- a/msprime/_msprimemodule.c +++ b/msprime/_msprimemodule.c @@ -151,6 +151,20 @@ get_dict_number(PyObject *dict, const char *key_str) return ret; } +static PyObject * +get_dict_string(PyObject *dict, const char *key_str) +{ + PyObject *ret = NULL; + PyObject *value; + value = get_required_dict_value(dict, key_str); + if (value == NULL) { + goto out; + } + ret = value; +out: + return ret; +} + static int parse_rate_map(PyObject *py_rate_map, size_t *ret_size, PyArrayObject **ret_position, PyArrayObject **ret_rate) @@ -1063,6 +1077,35 @@ Simulator_parse_sweep_genic_selection_model(Simulator *self, PyObject *py_model) return ret; } +static int +Simulator_parse_sweep_genic_selection_reverse_model(Simulator *self, PyObject *py_model) +{ + int ret = -1; + int err; + double position; + const char *filename; + PyObject *value; + value = get_dict_number(py_model, "position"); + if (value == NULL) { + goto out; + } + position = PyFloat_AsDouble(value); + value = get_dict_string(py_model, "filename"); + if (value == NULL) { + goto out; + } + filename = PyUnicode_AsUTF8(value); + err = msp_set_simulation_model_sweep_genic_selection_reverse(self->sim, + position, filename); + if (err != 0) { + handle_input_error("sweep genic selection reverse", err); + goto out; + } + ret = 0; +out: + return ret; +} + static int Simulator_parse_simulation_model(Simulator *self, PyObject *py_model) { @@ -1078,9 +1121,10 @@ Simulator_parse_simulation_model(Simulator *self, PyObject *py_model) PyObject *dirac_s = NULL; PyObject *beta_s = NULL; PyObject *sweep_genic_selection_s = NULL; + PyObject *sweep_genic_selection_reverse_s = NULL; PyObject *value; int is_hudson, is_dtwf, is_smc, is_smc_prime, is_smc_k, is_dirac, is_beta, - is_sweep_genic_selection, is_fixed_pedigree; + is_sweep_genic_selection, is_sweep_genic_selection_reverse, is_fixed_pedigree; double psi, c, alpha, truncation_point, hull_offset; hudson_s = Py_BuildValue("s", "hudson"); @@ -1120,6 +1164,11 @@ Simulator_parse_simulation_model(Simulator *self, PyObject *py_model) goto out; } + sweep_genic_selection_reverse_s = Py_BuildValue("s", "sweep_genic_selection_reverse"); + if (sweep_genic_selection_reverse_s == NULL) { + goto out; + } + py_name = get_required_dict_value(py_model, "name"); if (py_name == NULL) { goto out; @@ -1242,8 +1291,20 @@ Simulator_parse_simulation_model(Simulator *self, PyObject *py_model) } } + is_sweep_genic_selection_reverse = PyObject_RichCompareBool(py_name, + sweep_genic_selection_reverse_s, Py_EQ); + if (is_sweep_genic_selection_reverse == -1) { + goto out; + } + if (is_sweep_genic_selection_reverse) { + ret = Simulator_parse_sweep_genic_selection_reverse_model(self, py_model); + if (ret != 0) { + goto out; + } + } + if (! (is_hudson || is_dtwf || is_smc || is_smc_prime || is_smc_k || is_dirac - || is_beta || is_sweep_genic_selection || is_fixed_pedigree)) { + || is_beta || is_sweep_genic_selection || is_sweep_genic_selection_reverse || is_fixed_pedigree)) { PyErr_SetString(PyExc_ValueError, "Unknown simulation model"); goto out; } @@ -1262,6 +1323,7 @@ Simulator_parse_simulation_model(Simulator *self, PyObject *py_model) Py_XDECREF(beta_s); Py_XDECREF(dirac_s); Py_XDECREF(sweep_genic_selection_s); + Py_XDECREF(sweep_genic_selection_reverse_s); return ret; } @@ -2024,6 +2086,16 @@ Simulator_get_model(Simulator *self, void *closure) Py_DECREF(value); value = NULL; /* TODO fill in the parameters for the different types of trajectories. */ + } else if (model->type == MSP_MODEL_SWEEP_REVERSE) { + value = Py_BuildValue("d", model->params.sweep.position); + if (value == NULL) { + goto out; + } + if (PyDict_SetItemString(d, "locus", value) != 0) { + goto out; + } + Py_DECREF(value); + value = NULL; } else if (model->type == MSP_MODEL_SMC_K) { value = Py_BuildValue("d", model->params.smc_k_coalescent.hull_offset); if (value == NULL) { diff --git a/msprime/ancestry.py b/msprime/ancestry.py index 72b45b5d8..74af53c23 100644 --- a/msprime/ancestry.py +++ b/msprime/ancestry.py @@ -833,13 +833,6 @@ def _parse_sim_ancestry( models = _parse_model_arg(model) is_dtwf = isinstance(models[0], DiscreteTimeWrightFisher) is_pedigree = any(isinstance(model, FixedPedigree) for model in models) - is_smck = any(isinstance(model, SmcKApproxCoalescent) for model in models) - - if is_smck and gene_conversion_rate is not None: - raise ValueError( - "Gene conversion is not supported for the SmcKApproxCoalescent model. " - "Please refer to issue #2399 on GitHub for details." - ) if record_full_arg: if coalescing_segments_only is not None: @@ -1486,7 +1479,9 @@ def _choose_num_labels(self, models): """ num_labels = 1 for model in models: - if isinstance(model, SweepGenicSelection): + if isinstance(model, SweepGenicSelection) or isinstance( + model, SweepGenicSelectionReverse + ): num_labels = 2 return num_labels @@ -2116,3 +2111,50 @@ def __init__( self.end_frequency = end_frequency self.s = s self.dt = dt + + +@dataclasses.dataclass +class SweepGenicSelectionReverse(ParametricAncestryModel): + """ + This is a modification of SweepGenicSelection such that the user feeds the + filename for a binary file containing the forward in time trajectory of the sweep + for a single allele and the simulator computes the trajectory for the entire + genome. + + For more details see the definition of class SweepGenicSelection + + + .. warning:: + Currently models with more than one population and a selective sweep + are not implemented. Population size changes during the sweep + are not yet possible in msprime. + + :param float position: the location of the beneficial allele along the + chromosome. + :param filename: The path of the file that contains the forward trajectory + computed from a Guillespie simulation of a logistic birth and death + process stored in a binary file + Currently the file has the following format: + 1. Number of events (uintp) + 2. Number of demes (int) - currently 0 + 3. Population size (int) + 4. Migration rate (float) - currently 0 + 5. The initial state of the forward sim (1 int per deme) + 6. the final state of the forward sim (1 int per deme) + 7. For each event there is also: + a. time of event (float) + b. type of event (0 for wt birth, 1 for mutant birth - int) + c. deme of birth (int) + d. deme of parent (int) + """ + + name = "sweep_genic_selection_reverse" + + position: float | None + filename: str | None + + # We have to define an __init__ to enforce keyword-only behaviour + def __init__(self, *, duration=None, position=None, filename=None): + self.duration = duration + self.position = position + self.filename = filename diff --git a/notebooks/Forward_con_mig.ipynb b/notebooks/Forward_con_mig.ipynb new file mode 100644 index 000000000..bb3590760 --- /dev/null +++ b/notebooks/Forward_con_mig.ipynb @@ -0,0 +1,349 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "import argparse\n", + "import os\n", + "import math\n", + "import numpy as np\n", + "import random\n", + "import struct\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "class GillespieSelectionSimulation:\n", + " def __init__(self, L, N, s, migration_matrix, tfinal, seed=None, l0=0):\n", + " self.L = L # Total number of demes\n", + " self.N = N # Number of individuals in each deme\n", + " self.s = s # Selection coefficient\n", + " self.migration_matrix = migration_matrix\n", + " self.m = 0 #migration_matrix[0, 1]\n", + " self.tfinal = tfinal\n", + " self.l0 = l0 # Initial location for the mutant\n", + " np.random.seed(seed) # Setting the random seed\n", + " print('Initialization complete. Seed:', seed)\n", + "\n", + " def simulate(self):\n", + " t = 0.0\n", + " mutants = np.zeros(self.L, dtype = np.uint32)\n", + " mutants[self.l0] = 1 #Starting mutation at location l0\n", + " start_state = mutants.copy()\n", + " events = []\n", + " mut_traj = [[t, mutants.copy()]] # Initialize with the starting state\n", + "\n", + " print(\"Starting simulation...\")\n", + " while t < self.tfinal:\n", + "\n", + " if np.sum(mutants) == 0:\n", + " print(f\"All mutants extinct at time {t}. Reintroducing...\")\n", + " mutants = np.zeros(self.L, dtype = np.uint32)\n", + " mutants[self.l0] = 1\n", + " t = 0.0\n", + " events = []\n", + " mut_traj = [[t, mutants.copy()]] \n", + " continue\n", + " \n", + " \n", + " f = mutants / self.N\n", + " rate_mutant_replaces_wt = self.N * f * (1 - f) * (1 + self.s)\n", + " rate_wt_replaces_mutant = self.N * f * (1 - f)\n", + " rate_migration_mutant_to_i_from_j = np.zeros((self.L, self.L))\n", + " rate_migration_wt_to_i_from_j = np.zeros((self.L, self.L))\n", + "\n", + " for i in range(self.L):\n", + " for j in range(self.L):\n", + " if i != j:\n", + " # Calculate migration rates of mutants and wild types from j to i\n", + " rate_migration_mutant_to_i_from_j[i][j] = self.N * f[j] * (1 - f[i]) * (1 + self.s) * self.migration_matrix[i][j]\n", + " rate_migration_wt_to_i_from_j[i][j] = self.N * f[i] * (1 - f[j]) * self.migration_matrix[i][j]\n", + "\n", + " total_rate = np.sum(rate_mutant_replaces_wt + rate_wt_replaces_mutant) + np.sum(rate_migration_mutant_to_i_from_j) + np.sum(rate_migration_wt_to_i_from_j)\n", + "\n", + " if total_rate == 0:\n", + " print(\"No possible reactions. Skipping...\")\n", + " break\n", + "\n", + " tau = np.random.exponential(scale = 1 / total_rate)\n", + " t += tau\n", + "\n", + " combined_rates = np.concatenate((rate_mutant_replaces_wt, rate_wt_replaces_mutant, \n", + " rate_migration_mutant_to_i_from_j.flatten(), rate_migration_wt_to_i_from_j.flatten()))\n", + " reaction_probabilities = combined_rates / total_rate\n", + " cumulative_prob = np.cumsum(reaction_probabilities)\n", + " r = np.random.rand()\n", + " reaction_index = np.where(cumulative_prob > r)[0][0]\n", + "\n", + " if reaction_index < self.L:\n", + " event_type = 0\n", + " updated_deme = reaction_index\n", + " parent_deme = reaction_index\n", + " mutants[updated_deme] += 1\n", + " elif reaction_index < 2 * self.L:\n", + " event_type = 1\n", + " updated_deme = reaction_index - self.L\n", + " parent_deme = updated_deme\n", + " mutants[updated_deme] -= 1\n", + " elif reaction_index < 2 * self.L + self.L**2:\n", + " event_type = 2\n", + " index = reaction_index - 2 * self.L\n", + " i, j = index // self.L, index % self.L\n", + " updated_deme = i\n", + " parent_deme = j\n", + " mutants[i] += 1\n", + " else:\n", + " event_type = 3\n", + " index = reaction_index - (2 * self.L + self.L**2)\n", + " i, j = index // self.L, index % self.L\n", + " updated_deme = i\n", + " parent_deme = j\n", + " mutants[i] -= 1\n", + "\n", + " mut_traj.append([t, mutants.copy()])\n", + " events.append([t, event_type, updated_deme, parent_deme])\n", + " # print(f\"Event at time {t}: {event_type} in deme {updated_deme}, from deme {parent_deme}\")\n", + " end_state = mutants.copy()\n", + " demes = len(mutants)\n", + " return events, mut_traj, start_state, end_state\n", + " def save_events(self, events, start_state, end_state, filename):\n", + " with open(filename, 'wb') as f:\n", + " #Using struct.pack is not very efficient, but it seems to work fine for our sizes and I don't really care. It will cause problems for large deme sizes\n", + " # Write the number of events as size_t. Can be L (8 bytes) or I (4 bytes)\n", + " f.write(struct.pack('L', len(events))) # 'L' is long unsigned integer\n", + " f.write(struct.pack('i', self.L))\n", + " f.write(struct.pack('i', self.N)) # 'i' is int\n", + " f.write(struct.pack('d', self.m)) # 'd' is double\n", + " for i in start_state:\n", + " f.write(struct.pack('i', i))\n", + " for i in end_state:\n", + " f.write(struct.pack('i', i))\n", + " for i in events:\n", + " f.write(struct.pack('d', i[0]))\n", + " f.write(struct.pack('i', i[1]))\n", + " f.write(struct.pack('i', i[2]))\n", + " f.write(struct.pack('i', i[3]))\n", + " def load_events(self, filename):\n", + " with open(filename, 'rb') as f:\n", + " output = []\n", + " num_events = struct.unpack('L', f.read(8))\n", + " demes = struct.unpack('i', f.read(4))\n", + " pop_size = struct.unpack('i', f.read(4))\n", + " m = struct.unpack('d', f.read(8))\n", + " output.append(num_events) # Read and unpack one integers\n", + " output.append(demes)\n", + " output.append(pop_size)\n", + " output.append(m)\n", + " for i in range(demes[0]):\n", + " output.append(struct.unpack('i', f.read(4)))\n", + " output.append(struct.unpack('i', f.read(4)))\n", + " for i in range(num_events[0]):\n", + " output.append(struct.unpack('d', f.read(8)))\n", + " output.append(struct.unpack('i', f.read(4)))\n", + " output.append(struct.unpack('i', f.read(4)))\n", + " output.append(struct.unpack('i', f.read(4)))\n", + " return output" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "def plot_mutant_trajectories(mut_traj, L):\n", + " times = [state[0] for state in mut_traj]\n", + " mutant_counts = [state[1] for state in mut_traj]\n", + " plt.figure(figsize=(10, 6))\n", + " for deme in range(L):\n", + " plt.plot(times, [m[deme] for m in mutant_counts], label=f'Deme {deme+1}')\n", + " plt.xlabel('Time')\n", + " plt.ylabel('Number of Mutants')\n", + " plt.title('Mutant Trajectories Over Time')\n", + " #plt.legend()\n", + " plt.show()\n" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "class MigrationMatrixGenerator:\n", + " def __init__(self, deme_dimension, migration_rate):\n", + " self.deme_dimension = deme_dimension\n", + " self.migration_rate = migration_rate\n", + " \n", + "\n", + " def generate_1Dmatrix(self):\n", + " M = np.zeros((self.deme_dimension, self.deme_dimension))\n", + "\n", + " for i in range(self.deme_dimension):\n", + " if i > 0:\n", + " M[i, i - 1] = self.migration_rate\n", + " if i < self.deme_dimension - 1:\n", + " M[i, i + 1] = self.migration_rate\n", + "\n", + " return M\n", + " \n", + " \n", + " def generate_2Dmatrix(self):\n", + " \"\"\"\n", + " Generates a 2D migration matrix for a square grid of demes.\n", + " Each deme can migrate to its immediate neighbors (up, down, left, right).\n", + " Diagonal elements are set to 0, indicating no self-migration.\n", + " \"\"\"\n", + " z = math.sqrt(self.deme_dimension)\n", + " if not z.is_integer():\n", + " raise ValueError(\"Deme dimension must be a perfect square for a 2D grid.\")\n", + " z = int(z)\n", + "\n", + " M = np.zeros((self.deme_dimension, self.deme_dimension))\n", + "\n", + " for i in range(self.deme_dimension):\n", + " for j in range(self.deme_dimension):\n", + " # Check for right neighbor\n", + " if j == i + 1 and (i + 1) % z != 0:\n", + " M[i, j] = self.migration_rate\n", + " # Check for down neighbor\n", + " elif j == i + z and j < self.deme_dimension:\n", + " M[i, j] = self.migration_rate\n", + "\n", + " # Ensure migration is bidirectional\n", + " M[j, i] = M[i, j]\n", + "\n", + " return M\n", + "\n", + " def generate_3Dmatrix(self):\n", + " \"\"\"\n", + " Generates a 3D migration matrix for a cubic grid of demes.\n", + " Each deme can migrate to its immediate neighbors along the x, y, and z axes.\n", + " Diagonal elements are set to 0, indicating no self-migration.\n", + " \"\"\"\n", + " cube_root = round(self.deme_dimension ** (1/3))\n", + " if cube_root ** 3 != self.deme_dimension:\n", + " raise ValueError(\"Deme dimension must be a cube number for a 3D grid.\")\n", + "\n", + " M = np.zeros((self.deme_dimension, self.deme_dimension))\n", + "\n", + " for i in range(self.deme_dimension):\n", + " x, y, z = np.unravel_index(i, (cube_root, cube_root, cube_root))\n", + "\n", + " # Check neighbors in each direction (left, right, up, down, front, back)\n", + " for dx, dy, dz in [(-1, 0, 0), (1, 0, 0), (0, -1, 0), (0, 1, 0), (0, 0, -1), (0, 0, 1)]:\n", + " nx, ny, nz = x + dx, y + dy, z + dz\n", + " if 0 <= nx < cube_root and 0 <= ny < cube_root and 0 <= nz < cube_root:\n", + " neighbor_index = np.ravel_multi_index((nx, ny, nz), (cube_root, cube_root, cube_root))\n", + " M[i, neighbor_index] = self.migration_rate\n", + "\n", + " return M\n", + "\n", + " def generate_island_model_matrix(self):\n", + " \"\"\"\n", + " Generates a migration matrix based on the island model.\n", + " Each deme has an equal migration rate to every other deme.\n", + " Diagonal elements are set to 0, indicating no self-migration.\n", + " \"\"\"\n", + " M = np.full((self.deme_dimension, self.deme_dimension), self.migration_rate)\n", + " np.fill_diagonal(M, 0) # Set diagonal elements to 0\n", + " return M " + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Initialization complete. Seed: None\n", + "Starting simulation...\n", + "No possible reactions. Skipping...\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2QAAAIjCAYAAABswtioAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjUsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvWftoOwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAbNhJREFUeJzt3Xd8FHXixvFnN2XTC4QkBAKEIkiRqhQBRaKIWLCjnoJy6CmoiKfCqdgL2E5sqKfizztROU9UVBCpinTpvYROEloKCWm7398fISsLARLYZDbJ5/165ZXszHdnn4lrkoeZ+Y7NGGMEAAAAAKh0dqsDAAAAAEBNRSEDAAAAAItQyAAAAADAIhQyAAAAALAIhQwAAAAALEIhAwAAAACLUMgAAAAAwCIUMgAAAACwCIUMAAAAACxCIQMAVHmNGjXSoEGDrI5xRiZMmCCbzaZt27ZZHaVaePrpp2Wz2ayOAQBlRiEDgApW8ge3zWbTb7/9dsJ6Y4wSExNls9l05ZVXntFr/Pjjj3r66afPMmnZ/P7773r66aeVkZFxynGzZ8927/fpPnzdiy++qMmTJ1sdw+vmzZuna6+9VnFxcXI4HGrUqJHuuece7dixw+poHho1alSm99GECROsjgoA5WYzxhirQwBAdTZhwgTdeeedCgoK0p133ql3333XY/3s2bPVq1cvORwOJScna8qUKeV+jWHDhumdd95RZfxIf/XVV/XII48oJSVFjRo1Oum4tLQ0TZ8+3WPZqFGjFBYWpscff9xj+V/+8pezypSfny+73a6AgICz2s7JhIWF6YYbbqiQP/idTqcKCwvlcDgqtZy+9dZbevDBB9W4cWMNGjRIdevW1bp16/Svf/1LUnHJ79atW6XlOZXJkyfr8OHD7sc//vijJk6cqDfeeEMxMTHu5d26dVODBg1UVFSkoKAgK6ICQLn5Wx0AAGqKK664QpMmTdK4cePk7//nj9/PP/9cHTt21P79+y1M531xcXEnFK2XX35ZMTExpyxgLpdLBQUF5fqD2uFwnHFOq+Tk5Cg0NFR+fn7y8/Or1NeeN2+ehg8fru7du2vq1KkKCQlxr7v33nt14YUX6oYbbtCaNWsUHR1dablKvifH69+/v8fj1NRUTZw4Uf379y/1HwWO/f8LAHwdpywCQCW55ZZbdODAAY+jRgUFBfrvf/+rW2+99YTxJaf8zZ4922P5tm3bPE7PGjRokN555x1JKvU0wFdffVXdunVT7dq1FRwcrI4dO+q///3vCa9ns9k0bNgwTZ48Wa1bt5bD4VCrVq00depU95inn35ajzzyiCQpKSnJ/Vpnc/1Tyev+5z//UatWreRwONyvWdbspV1DlpGRoeHDhysxMVEOh0NNmzbVmDFj5HK5PMa5XC69+eabatOmjYKCglSnTh1dfvnlWrJkiTtfTk6OPv30U/f+Hvtay5YtU9++fRUREaGwsDD17t1bCxYs8HiNktNW58yZo/vuu0+xsbGqX7++x7rjv4c//fSTevToodDQUIWHh6tfv35as2aNx5jU1FTdeeedql+/vhwOh+rWratrrrnmtP89nnvuOdlsNn366aceZUySmjRporFjx2rv3r16//333f8dbDabtm/ffsK2Ro0apcDAQB06dMi9bOHChbr88ssVGRmpkJAQXXTRRZo3b57H80qu9Vq7dq1uvfVWRUdHq3v37qfMXRalXUNW8h6bNGmSWrZsqeDgYHXt2lWrVq2SJL3//vtq2rSpgoKCdPHFF5f6/SvLPgHAmaCQAUAladSokbp27aqJEye6l/3000/KzMzUgAEDzni799xzjy699FJJ0meffeb+KPHmm2+qffv2evbZZ/Xiiy/K399fN954o3744YcTtvXbb7/pvvvu04ABAzR27Fjl5eXp+uuv14EDByRJ1113nW655RZJ0htvvOF+rTp16pxxfkmaOXOmHnroId18881688033Uc9ypP9WLm5ubrooov073//W3fccYfGjRunCy+8UKNGjdKIESM8xg4ePNhd3MaMGaORI0cqKCjIXao+++wzORwO9ejRw72/99xzjyRpzZo16tGjh1asWKFHH31UTz75pFJSUnTxxRdr4cKFJ+S67777tHbtWo0ePVojR448af7PPvtM/fr1U1hYmMaMGaMnn3xSa9euVffu3T3KwvXXX69vvvnGfSrsAw88oOzs7FNeA5abm6sZM2aoR48eSkpKKnXMzTffLIfD4T599qabbpLNZtNXX311wtivvvpKl112mftI2syZM9WzZ09lZWXpqaee0osvvqiMjAxdcsklWrRo0QnPv/HGG5Wbm6sXX3xRQ4YMOWnus/Xrr7/q4Ycf1sCBA/X0009r3bp1uvLKK/XOO+9o3Lhxuu+++/TII49o/vz5uuuuuzyeW959AoByMQCACvXJJ58YSWbx4sXm7bffNuHh4SY3N9cYY8yNN95oevXqZYwxpmHDhqZfv37u582aNctIMrNmzfLYXkpKipFkPvnkE/eyoUOHmpP9SC95rRIFBQWmdevW5pJLLvFYLskEBgaazZs3u5etWLHCSDJvvfWWe9krr7xiJJmUlJQyfw9KtGrVylx00UUnvK7dbjdr1qw54+wNGzY0AwcOdD9+7rnnTGhoqNm4caPHuJEjRxo/Pz+zY8cOY4wxM2fONJLMAw88cMJru1wu99ehoaEe2y/Rv39/ExgYaLZs2eJetmfPHhMeHm569uzpXlbyHujevbspKiry2EbJupLvZ3Z2tomKijJDhgzxGJeammoiIyPdyw8dOmQkmVdeeeWEXKeyfPlyI8k8+OCDpxx33nnnmVq1arkfd+3a1XTs2NFjzKJFi4wk83//93/GmOLvWbNmzUyfPn08vn+5ubkmKSnJXHrppe5lTz31lJFkbrnllnLlN+bU78GS7R5LknE4HB7j33//fSPJxMfHm6ysLPfyUaNGeWy7PPsEAGeCI2QAUIluuukmHTlyRFOmTFF2dramTJlS6umK3hQcHOz++tChQ8rMzFSPHj30xx9/nDA2OTlZTZo0cT8+77zzFBERoa1bt1ZoxosuukgtW7Y8YXl5sh9r0qRJ6tGjh6Kjo7V//373R3JyspxOp+bOnStJ+vrrr2Wz2fTUU0+dsI3TTbDhdDr1888/q3///mrcuLF7ed26dXXrrbfqt99+U1ZWlsdzhgwZctrrxaZPn66MjAzdcsstHtn9/PzUuXNnzZo1y/29CQwM1OzZsz1OFzyd7OxsSVJ4ePgpx4WHh3vkv/nmm7V06VJt2bLFvezLL7+Uw+HQNddcI0lavny5Nm3apFtvvVUHDhxwZ8/JyVHv3r01d+7cE04Z/dvf/lbm7Gejd+/eHtebde7cWVLxUcZjvxcly0ve82eyTwBQHlz1CgCVqE6dOkpOTtbnn3+u3NxcOZ1O3XDDDRX6mlOmTNHzzz+v5cuXKz8/3728tMLRoEGDE5ZFR0eX6w/+M3GyU+fKk/1YmzZt0sqVK096KmV6erokacuWLUpISFCtWrXKnXnfvn3Kzc1V8+bNT1h37rnnyuVyaefOnWrVqpV7+cn28/jsknTJJZeUuj4iIkJS8UQmY8aM0cMPP6y4uDh16dJFV155pe644w7Fx8efdPsl5aOkmJ1Mdna2R1G58cYbNWLECH355Zf6xz/+IWOMJk2a5L5+7tjsAwcOPOl2MzMzPSYKKcv3xBuOf29HRkZKkhITE0tdXvKeP5N9AoDyoJABQCW79dZbNWTIEKWmpqpv376KiooqddzJSofT6Szza/3666+6+uqr1bNnT7377ruqW7euAgIC9Mknn+jzzz8/YfzJjt6YCp5O/9gjYSXKm/1YLpdLl156qR599NFS159zzjleyV1epe3n8UqOtnz22WelFqtjZxAcPny4rrrqKk2ePFnTpk3Tk08+qZdeekkzZ85U+/btS91+06ZN5e/vr5UrV540Q35+vjZs2KBOnTq5lyUkJKhHjx766quv9I9//EMLFizQjh07NGbMmBOyv/LKK2rXrl2p2w4LC/N4XJbviTec7L19uvf8mewTAJQHhQwAKtm1116re+65RwsWLNCXX3550nEl/+J+/A2YS5vp7mTl7euvv1ZQUJCmTZvmMTX8J598cgbJT/1a3nY22Zs0aaLDhw8rOTn5tOOmTZumgwcPnvIoWWn7XKdOHYWEhGjDhg0nrFu/fr3sdvsJR1/KouSU0djY2NPmLxn/8MMP6+GHH9amTZvUrl07vfbaa/r3v/9d6vjQ0FD16tVLM2fO1Pbt29WwYcMTxnz11VfKz88/4UblN998s+677z5t2LBBX375pUJCQnTVVVedkD0iIqJM2auC6rhPAHwL15ABQCULCwvTe++9p6efftrjj9njNWzYUH5+fu7rnUocf2NpSe57Nx1f3vz8/GSz2TyOqm3btk2TJ08+4/wney1vO5vsN910k+bPn69p06adsC4jI0NFRUWSiq8fMsbomWeeOWHcsUcFQ0NDS/3eXnbZZfr22289Zj5MS0vT559/ru7du7tP5SuPPn36KCIiQi+++KIKCwtPWL9v3z5JxbMl5uXleaxr0qSJwsPDPU7vLM0TTzwhY4wGDRqkI0eOeKxLSUnRo48+qrp167pnkyxx/fXXy8/PTxMnTtSkSZN05ZVXetw3rGPHjmrSpIleffVVjxs5H5+9KqmO+wTAt3CEDAAscKrrUUpERkbqxhtv1FtvvSWbzaYmTZpoypQp7uufjtWxY0dJ0gMPPKA+ffrIz89PAwYMUL9+/fT666/r8ssv16233qr09HS98847atq06SlPWTuVktd6/PHHNWDAAAUEBOiqq64q9Ya+Z+Nssj/yyCP67rvvdOWVV2rQoEHq2LGjcnJytGrVKv33v//Vtm3bFBMTo169eun222/XuHHjtGnTJl1++eVyuVz69ddf1atXLw0bNsy9z7/88otef/11JSQkKCkpSZ07d9bzzz+v6dOnq3v37rrvvvvk7++v999/X/n5+Ro7duwZ7XdERITee+893X777erQoYMGDBigOnXqaMeOHfrhhx904YUX6u2339bGjRvVu3dv3XTTTWrZsqX8/f31zTffKC0t7bS3UejZs6deffVVjRgxQuedd54GDRqkunXrav369frwww/lcrn0448/nnBdVGxsrHr16qXXX39d2dnZuvnmmz3W2+12/etf/1Lfvn3VqlUr3XnnnapXr552796tWbNmKSIiQt9///0ZfV+sUh33CYCPsXKKRwCoCY6d9v5Ujp/23hhj9u3bZ66//noTEhJioqOjzT333GNWr159wrT3RUVF5v777zd16tQxNpvNY9rvjz76yDRr1sw4HA7TokUL88knn5x0avChQ4eWmuv4Kd+fe+45U69ePWO328s1Bf7Jpr0v7XXLk720jNnZ2WbUqFGmadOmJjAw0MTExJhu3bqZV1991RQUFLjHFRUVmVdeecW0aNHCBAYGmjp16pi+ffuapUuXusesX7/e9OzZ0wQHBxtJHq/1xx9/mD59+piwsDATEhJievXqZX7//XePLKd6Dxw/7X2JWbNmmT59+pjIyEgTFBRkmjRpYgYNGmSWLFlijDFm//79ZujQoaZFixYmNDTUREZGms6dO5uvvvqq1O9laebOnWuuueYaExMTYwICAkyDBg3MkCFDzLZt2076nA8//NBIMuHh4ebIkSOljlm2bJm57rrrTO3atY3D4TANGzY0N910k5kxY4Z7TMl/x3379pU5b4kzmfb++PdYye0jjr9tQMntJiZNmlTufQKAM2EzpoKv1AYAoIIlJiaqT58++te//mV1FAAAyoVryAAAVVphYaEOHDigmJgYq6MAAFBuXEMGAKiypk2bpi+++EJHjhxR7969rY4DAEC5ccoiAKDK6tWrlzZv3qx7771X//jHP6yOAwBAuVHIAAAAAMAiXEMGAAAAABahkAEAAACARZjUw0tcLpf27Nmj8PBw2Ww2q+MAAAAAsIgxRtnZ2UpISJDdfupjYBQyL9mzZ48SExOtjgEAAADAR+zcuVP169c/5RgKmZeEh4dLKv6mR0REWJwGAAAAgFWysrKUmJjo7ginQiHzkpLTFCMiIihkAAAAAMp0KROTegAAAACARShkAAAAAGARChkAAAAAWIRCBgAAAAAWoZABAAAAgEUoZAAAAABgEQoZAAAAAFiEQgYAAAAAFqGQAQAAAIBFKGQAAAAAYBEKGQAAAABYhEIGAAAAABahkAEAAACARShkAAAAAGARSwvZ3LlzddVVVykhIUE2m02TJ0/2WG+M0ejRo1W3bl0FBwcrOTlZmzZt8hhz8OBB3XbbbYqIiFBUVJQGDx6sw4cPe4xZuXKlevTooaCgICUmJmrs2LEnZJk0aZJatGihoKAgtWnTRj/++KPX9xcAAAAAjmVpIcvJyVHbtm31zjvvlLp+7NixGjdunMaPH6+FCxcqNDRUffr0UV5ennvMbbfdpjVr1mj69OmaMmWK5s6dq7vvvtu9PisrS5dddpkaNmyopUuX6pVXXtHTTz+tDz74wD3m999/1y233KLBgwdr2bJl6t+/v/r376/Vq1dX3M4DAAAAqPFsxhhjdQhJstls+uabb9S/f39JxUfHEhIS9PDDD+vvf/+7JCkzM1NxcXGaMGGCBgwYoHXr1qlly5ZavHixOnXqJEmaOnWqrrjiCu3atUsJCQl677339Pjjjys1NVWBgYGSpJEjR2ry5Mlav369JOnmm29WTk6OpkyZ4s7TpUsXtWvXTuPHjy9T/qysLEVGRiozM1MRERHe+rYAAACgmjHGyBjJaYxcR792GSOny8hlitd7fH10jDnmucXbOfpZnuuLl5Ws/3NdyXPc449+7bGtY9b/uR3jsb0/xx39bMwJ29fJXvO4jH9+U/7cl+PzHDv2+Of+2WSKv7DZbOrTKv6U3//KUJ5u4F9JmcotJSVFqampSk5Odi+LjIxU586dNX/+fA0YMEDz589XVFSUu4xJUnJysux2uxYuXKhrr71W8+fPV8+ePd1lTJL69OmjMWPG6NChQ4qOjtb8+fM1YsQIj9fv06fPCadQHis/P1/5+fnux1lZWV7YawAAAHiD02WUW1CkIwVO5RQ4lZNfpOy8Ih3OL1JOfpGy8gqVlpWnIwUuFTidyi90qcDp+vNzkVMFRS7lF7lUUOSSyxQXJJfLuL92uoxMydfmmK9dJxYt99fHFCp4n8Pfrg3P97U6Rrn4bCFLTU2VJMXFxXksj4uLc69LTU1VbGysx3p/f3/VqlXLY0xSUtIJ2yhZFx0drdTU1FO+TmleeuklPfPMM2ewZwAAADiZgiKXMnILdCi3UAdzCtxfH8ot0KGcAh3MKdCh3ALlFbpU6HSp0GVU5HSpyGlU4HQp80ihcguKlFfosnpXzprNJtltNtltxUd+JMl2dLlNtqOfS8bair8+fpmt5Dm2Y5YXD/xz3dHl8hx/bI6yvObJMuq48SXbK3nNY1+n5PnHLrCdZL2tlPUBflVvzkKfLWS+btSoUR5H1bKyspSYmGhhIgAAAOsZY3Q4v0gZuYU6nF+krCOF2n+4QJlHCpWVV6isI4XKzisq5esiZecVKqfA6dU8dpsUGuivUIe/woKOfnb4Kczhr1qhgYoKCVSgn12OAHvxZ3+7HP5+CvQv/jrQ364AP7v87TbZjpYjv+O+th8tPiVf291FynZ0rE78+uh6m11/fn3cNo4tRai+fLaQxccXn/uZlpamunXrupenpaWpXbt27jHp6ekezysqKtLBgwfdz4+Pj1daWprHmJLHpxtTsr40DodDDofjDPYMAACg6jDG6EihU1lHinQot0Bb9h1Wela+9h/O14HDBUrNylNqZp4O5xe5P5yuszsnz26TokICFRUSoOiQwKMfAYoODVSt0EDVCglUUKCfAuw2+fvZFeBnc5emyJAAdwELCfSTw99OsYFP89lClpSUpPj4eM2YMcNdwLKysrRw4ULde++9kqSuXbsqIyNDS5cuVceOHSVJM2fOlMvlUufOnd1jHn/8cRUWFiogIECSNH36dDVv3lzR0dHuMTNmzNDw4cPdrz99+nR17dq1kvYWAADAOoVOl9Kz87V9f442pmVrQ9ph7TqUq92HjmjHwVwVnUHBCgqwK8zh7z4SVSs0UJHBgQoP8ldEcIAigvwVERSgiGB/hQcFKCIoQOFB/ooKKf7abqdEoWawtJAdPnxYmzdvdj9OSUnR8uXLVatWLTVo0EDDhw/X888/r2bNmikpKUlPPvmkEhIS3DMxnnvuubr88ss1ZMgQjR8/XoWFhRo2bJgGDBighIQESdKtt96qZ555RoMHD9Zjjz2m1atX680339Qbb7zhft0HH3xQF110kV577TX169dPX3zxhZYsWeIxNT4AAEBVlVfo1OrdmdqUflj7s4uPbu0/XKA9mUe08+ARHcjJP+1EE352myKDA5QQFaSGtUNVJ8yh2qGBigl3qH50sMKDAhTm8Fd4UPFHSKDP/rs/4FMsnfZ+9uzZ6tWr1wnLBw4cqAkTJsgYo6eeekoffPCBMjIy1L17d7377rs655xz3GMPHjyoYcOG6fvvv5fdbtf111+vcePGKSwszD1m5cqVGjp0qBYvXqyYmBjdf//9euyxxzxec9KkSXriiSe0bds2NWvWTGPHjtUVV1xR5n1h2nsAAOALnC6jlP2H9ceODM3fckDrU7O1OT1bhc5T/8nnb7epfnSwmsaGq3l8mBrWDlV8RJCaxIYpOiRAwQF+nPoHlFF5uoHP3IesqqOQAQAAKxzOL9LCrQf02+b9mrNhn7buzyl1XEyYQ63rRSg23KGYsOKP2AiHGtUOVXxkkGqFBHKaIOAl1eI+ZAAAAChdelaevl+5V18s2qFN6YdLHdOpYbQ6N66l9onRah4frvrRwRzhAnwQhQwAAMDHGWO0Kf2wvl2+Wz+vSVPK/pwTJtq4tXMDdW1cW01jw3ROXLj8ONoFVAkUMgAAAB+TV+jUb5v2KzUrT8t3ZmjGujQdyi30GNO+QZSuOi9BPc+poyZ1Qjn6BVRRFDIAAAAfsW1/jp6dslYz16efsM7hb1fXJrV1Q8f6at8gWvWigi1ICMDbKGQAAAAWMsZoxa5MTV62W18u3qkjhU5JUkSQv+pFh+iCRtHq26auOjSIVqC/3eK0ALyNQgYAAGCBvEKn5m3er3dnb9HS7Yfcyzs0iNIL17ZRi/hwTkMEagAKGQAAQCXafiBH783eosnLdyuv0CVJCvCz6ZIWsbr5/ET1ah5LEQNqEAoZAABABcvJL9LS7Yf0zbLd+m7FHjmPzpBYNzJIl7aM030XN1V8ZJDFKQFYgUIGAABQQdanZundWVs0Y12acgqc7uUXnVNHwy5pqk4NozkaBtRwFDIAAAAv25yereemrNOcjfvcy+pFBatL49oa2K2hzqsfZV04AD6FQgYAAOAleYVOvTt7i96bvVmFTiO7TbqsZbzuvLCRLkiqxdEwACegkAEAAHjBgq0H9I//rdLW/TmSpN4tYjX6qpZqWDvU4mQAfBmFDAAA4CzsOJCrZ6es1S/r0iRJseEOPX11K/VtHc8RMQCnRSEDAAA4Qz+vSdXDX61Qdn6RJOkvXRro0ctbKCIowOJkAKoKChkAAEA5bd13WK9P36gpK/dKklrXi9BrN7ZT8/hwi5MBqGooZAAAAOXw46q9eujL5covcslmk27r3ED/uOJchQTyZxWA8uMnBwAAQBkYY/TenC0aO3WDJKlzUi090a+l2tSPtDgZgKqMQgYAAHAaOw/m6rWfN2jy8j2SpEHdGunJK1vKz86kHQDODoUMAADgJFyu4qNib0zfqCJX8X3FnrqqlQZ2a2R1NADVBIUMAACgFLsO5Wrk16v02+b9kqTzG0XrsctbqFOjWhYnA1CdUMgAAACOYYzRD6v26vFvVivzSKH87TY9dnkL/bVHEvcVA+B1FDIAAICjnC6jUf9bqa+W7JIknVc/Us9c3UrtG0RbnAxAdUUhAwAAOGr8nC3uMvbX7kl69PIWCvS3W5wKQHVGIQMAADWey2X0zqzNeuOXjZKkp65qqTsvTLI4FYCagEIGAABqtOy8Qj34xXLNXJ8uSRrYtaEGMYsigEpCIQMAADVWbkGRBn+6RItSDsrhb9dz/Vvrpk6JVscCUINQyAAAQI2UkVug2/61UGv2ZCkowK6JQ7oweQeASsdVqgAAoMbJK3Tqr58u0Zo9WaodGqh/D+5MGQNgCY6QAQCAGqXI6dL9E5dpyfZDCg/y1+dDuqh5fLjVsQDUUBwhAwAANUZeoVN/n7RC09emKdDfrg/v6EQZA2ApjpABAIAaISO3QHd8vEgrd2XK327TuAHt1KVxbatjAajhKGQAAKDay8ordJex6JAAvXFzO13cPNbqWABAIQMAANVbQZFLgycs1spdmaoVGqgv7u6ic+I4TRGAb+AaMgAAUK09N2WtFm8rnsDjs8EXUMYA+BQKGQAAqLYmLdmpzxZsl80mjRvQXq0SIq2OBAAeKGQAAKBaWrUrU49PXi1JGt77HPVqwTVjAHwPhQwAAFQ7szek69YPF6igyKXeLWJ1/yVNrY4EAKViUg8AAFCt/LZpvwZ/ukROl9H5jaL1+s3tZLfbrI4FAKWikAEAgGojPTtPI75aLqfLqG/reP1zQDs5/P2sjgUAJ8UpiwAAoFooKHLpvn//ofTsfDWNDdOYG86jjAHweRQyAABQ5Rlj9PT3a7Rke/H09h/c3lERQQFWxwKA06KQAQCAKu+LxTv1+cIdstmkNwe0U+M6YVZHAoAyoZABAIAqbfuBHD03Za0k6ZE+zXVJiziLEwFA2VHIAABAlVXkdOn+icuUW+BU56Ra+lvPJlZHAoByoZABAIAq65N527RyV6ZCA/300nVtmN4eQJVDIQMAAFXS9gM5em36BknSU1e14roxAFUShQwAAFQ5xhj945tVyit0qVuT2rqxU32rIwHAGaGQAQCAKufzRTs0b/MBOfztevHaNrLZOFURQNVEIQMAAFXK5vRsvfTjeknFsyo2igm1OBEAnDkKGQAAqDLSs/M06JPFOpxfpE4NozWwWyOrIwHAWaGQAQCAKqGgyKWh//lDuw4dUb2oYL1/e0cF+PGnDICqjZ9iAACgSnjt5w1avO2Qwhz++njQ+aod5rA6EgCcNQoZAADweZvSsvXJvG2SpJevb6Pm8eHWBgIAL6GQAQAAn2aM0d8nrVCB06VezeuoX5u6VkcCAK+hkAEAAJ82Y126VuzKVHCAn8bccB5T3AOoVihkAADAZx0pcOrFH9dJkgZ2a6TY8CCLEwGAd1HIAACAz/p4Xoq27s9RbLhD9/RsbHUcAPA6ChkAAPBJuzOO6K2ZmyRJI/u2UHRooMWJAMD7KGQAAMAnPT9lrfIKXbqgUS1d276e1XEAoEJQyAAAgM9ZsPWAflqdKj+7Tc9c04qJPABUWxQyAADgU4wxenXaBknSTZ0SdW7dCIsTAUDFoZABAACf8v3KvVqy/ZCCA/z0QO+mVscBgApFIQMAAD7jSIFTLx2d5v6+i5uobmSwxYkAoGJRyAAAgM94+rs12puZp3pRwRrCNPcAagAKGQAA8AkrdmboyyU7JUlPX91KQQF+FicCgIpHIQMAAD7hi8XFZezqtgm6tGWcxWkAoHJQyAAAgOVW7MzQl4t3SJIGXJBocRoAqDwUMgAAYKmCIpce/e9KuUzx0bFuTWKsjgQAlYZCBgAALPXu7M3akJatWqGBeuqqllbHAYBKRSEDAACWycor1PtztkoqnsijdpjD4kQAULkoZAAAwDITF+7QkUKnzokL01Xn1bU6DgBUOgoZAACwRHZeocbP2SJJ+muPxrLZbBYnAoDKRyEDAACWeHvWZh3KLVTjOqG6rn09q+MAgCV8upA5nU49+eSTSkpKUnBwsJo0aaLnnntOxhj3GGOMRo8erbp16yo4OFjJycnatGmTx3YOHjyo2267TREREYqKitLgwYN1+PBhjzErV65Ujx49FBQUpMTERI0dO7ZS9hEAgJpo9e5MffxbiiTp8SvOlb+fT/9JAgAVxqd/+o0ZM0bvvfee3n77ba1bt05jxozR2LFj9dZbb7nHjB07VuPGjdP48eO1cOFChYaGqk+fPsrLy3OPue2227RmzRpNnz5dU6ZM0dy5c3X33Xe712dlZemyyy5Tw4YNtXTpUr3yyit6+umn9cEHH1Tq/gIAUFO8MX2jCp1GyefGqfe53AQaQM1lM8cebvIxV155peLi4vTRRx+5l11//fUKDg7Wv//9bxljlJCQoIcfflh///vfJUmZmZmKi4vThAkTNGDAAK1bt04tW7bU4sWL1alTJ0nS1KlTdcUVV2jXrl1KSEjQe++9p8cff1ypqakKDAyUJI0cOVKTJ0/W+vXry5Q1KytLkZGRyszMVEREhJe/EwAAVB+/b96vW/+1UDab9MuIi9SkTpjVkQDAq8rTDXz6CFm3bt00Y8YMbdy4UZK0YsUK/fbbb+rbt68kKSUlRampqUpOTnY/JzIyUp07d9b8+fMlSfPnz1dUVJS7jElScnKy7Ha7Fi5c6B7Ts2dPdxmTpD59+mjDhg06dOhQqdny8/OVlZXl8QEAAE7vnzOKLy24sWN9yhiAGs/f6gCnMnLkSGVlZalFixby8/OT0+nUCy+8oNtuu02SlJqaKkmKi/M81SEuLs69LjU1VbGxsR7r/f39VatWLY8xSUlJJ2yjZF10dPQJ2V566SU988wzXthLAABqjo1p2VqUclD+dpuGJ59jdRwAsJxPHyH76quv9J///Eeff/65/vjjD3366ad69dVX9emnn1odTaNGjVJmZqb7Y+fOnVZHAgDApxlj9NKP6yRJFzevo4SoYIsTAYD1fPoI2SOPPKKRI0dqwIABkqQ2bdpo+/bteumllzRw4EDFx8dLktLS0lS37p83k0xLS1O7du0kSfHx8UpPT/fYblFRkQ4ePOh+fnx8vNLS0jzGlDwuGXM8h8Mhh8Nx9jsJAEANMWtDumZt2KdAf7tG9m1hdRwA8Ak+fYQsNzdXdrtnRD8/P7lcLklSUlKS4uPjNWPGDPf6rKwsLVy4UF27dpUkde3aVRkZGVq6dKl7zMyZM+VyudS5c2f3mLlz56qwsNA9Zvr06WrevHmppysCAIDyMcbozRmbJUl3dmukprHhFicCAN/g04Xsqquu0gsvvKAffvhB27Zt0zfffKPXX39d1157rSTJZrNp+PDhev755/Xdd99p1apVuuOOO5SQkKD+/ftLks4991xdfvnlGjJkiBYtWqR58+Zp2LBhGjBggBISEiRJt956qwIDAzV48GCtWbNGX375pd58802NGDHCql0HAKBaWbErUyt2ZijQz64hPRtbHQcAfIZPn7L41ltv6cknn9R9992n9PR0JSQk6J577tHo0aPdYx599FHl5OTo7rvvVkZGhrp3766pU6cqKCjIPeY///mPhg0bpt69e8tut+v666/XuHHj3OsjIyP1888/a+jQoerYsaNiYmI0evRoj3uVAQCAM/fvBdslSZe2ilNMGKf8A0AJn74PWVXCfcgAAChdQZFLnZ6frqy8In1xdxd1aVzb6kgAUKGqzX3IAABA1Td/6wFl5RUpJixQ5zeqZXUcAPApFDIAAFChPvotRZLUt3Vd+dltFqcBAN9CIQMAABVm/pYDmrtxn/ztNv21R5LVcQDA51DIAABAhfnw162SpAEXJKph7VCL0wCA76GQAQCACpGyP0ezN6RLku66kKNjAFAaChkAAKgQ783eLJeRLmkRq8Z1wqyOAwA+iUIGAAC8bnfGEf3vj92SpKG9mlqcBgB8F4UMAAB43Ydzt6rIZdS1cW11bBhtdRwA8FkUMgAA4FX7svM1cdEOSdKwSzg6BgCnQiEDAABe9a/ftiq/yKV2iVHq1qS21XEAwKdRyAAAgNfsOJCrj4/eCPr+S5rKZuNG0ABwKhQyAADgNR/9tlWFTqPuTWN0SYtYq+MAgM+jkAEAAK/IyC3QF4t3SpLuvbgJR8cAoAwoZAAAwCu+/mO38otcOrduBNeOAUAZUcgAAMBZM8a4Z1a8tXMDjo4BQBlRyAAAwFn7fcsBbU4/rJBAP/Vvl2B1HACoMihkAADgrH10dGbFGzvWV3hQgMVpAKDqoJABAICzsm1/jmauT5fNJt15YZLVcQCgSqGQAQCAs/Lt8j2SpB7N6qhRTKjFaQCgaqGQAQCAM5ZX6NRnC7ZJkq5tz7VjAFBeFDIAAHDGZm/Yp/2HC1Q7NFD92lDIAKC8KGQAAOCMTVpSfCPoGzrVV6A/f1YAQHnxkxMAAJyR9Ow8zd64T5J0U6dEi9MAQNVEIQMAAGfku+V75HQZdWgQpSZ1wqyOAwBVEoUMAACUmzFGk5bskiRd1ZZrxwDgTFHIAABAuU1bk6oNadkKCrDr2vb1rI4DAFUWhQwAAJTbP3/ZJEka1C1JUSGBFqcBgKqLQgYAAMpl5a4MrU/NVqCfXX+7qLHVcQCgSqOQAQCAcvlu+R5JUp/W8RwdA4CzRCEDAABldqTAqe9WFBeyfm3qWpwGAKo+ChkAACiz//6xS+nZ+aobGaSLm9exOg4AVHkUMgAAUCbGGP1nwXZJ0t09GysowM/iRABQ9VHIAABAmSzZfkjrU7MV4GfT1dx7DAC8gkIGAADK5POFOyQVXztWO8xhcRoAqB4oZAAA4LTSsvL0/dHJPAZ3Z6p7APAWChkAADit71fsUZHLqGPDaLWpH2l1HACoNihkAADglIwx+n7lXknSFUx1DwBeRSEDAACntHxnhlbszFCgn13XtGMyDwDwJgoZAAA4pX8vODqZx3l1FcNkHgDgVRQyAABwUin7c/Tt8t2SpDu6NrQ4DQBUPxQyAABwUj+u2qsil1G3JrXVvkG01XEAoNqhkAEAgFIZY9xHx/oymQcAVAgKGQAAKNXsjfu0Me2wHP52XdE63uo4AFAtUcgAAECp/rNguyTp2vb1VJvJPACgQlDIAADACfZmHtHM9emSpCE9G1ucBgCqLwoZAAA4wbgZm+UyUuekWmpSJ8zqOABQbVHIAACAh7xCp75fsUeS9GDvZhanAYDqjUIGAAA8/L5lvw7nFykuwqEujWtbHQcAqjUKGQAA8DB5WfHRsctbxctut1mcBgCqNwoZAABwO5xfpJ/XpkqSru1Q3+I0AFD9UcgAAIDbtNWpyit0qXFMqNrWj7Q6DgBUexQyAADg9s2y3ZKk/u3ryWbjdEUAqGgUMgAAIElKy8rTvC37JUn929WzOA0A1AzlLmQ7d+7Url273I8XLVqk4cOH64MPPvBqMAAAULm+W75HxkidGkarQe0Qq+MAQI1Q7kJ26623atasWZKk1NRUXXrppVq0aJEef/xxPfvss14PCAAAKsf/jp6ueG0Hjo4BQGUpdyFbvXq1LrjgAknSV199pdatW+v333/Xf/7zH02YMMHb+QAAQCXYkJqtdXuzFOhnV782da2OAwA1RrkLWWFhoRwOhyTpl19+0dVXXy1JatGihfbu3evddAAAoFL8vKZ4qvue58QoKiTQ4jQAUHOUu5C1atVK48eP16+//qrp06fr8ssvlyTt2bNHtWvX9npAAABQ8X5emyZJuqRFnMVJAKBmKXchGzNmjN5//31dfPHFuuWWW9S2bVtJ0nfffec+lREAAFQdB3MKtHpPpiSp97mxFqcBgJrFv7xPuPjii7V//35lZWUpOjravfzuu+9WaGioV8MBAICKN2NdmoyRzq0bobiIIKvjAECNUu4jZJdccomys7M9ypgk1apVSzfffLPXggEAgMrx3Yo9kqS+reMtTgIANU+5C9ns2bNVUFBwwvK8vDz9+uuvXgkFAAAqx5ECp5ZsOyRJuqwV148BQGUr8ymLK1eudH+9du1apaamuh87nU5NnTpV9epx3xIAAKqSORvTdaTQqXpRwWoeF251HACoccpcyNq1ayebzSabzaZLLrnkhPXBwcF66623vBoOAABUrOlr0yUVn65os9ksTgMANU+ZC1lKSoqMMWrcuLEWLVqkOnXquNcFBgYqNjZWfn5+FRISAAB4nzFGv27aJ0nq1YLZFQHACmUuZA0bNpQkuVyuCgsDAAAqz6rdmUrPzldIoJ86Now+/RMAAF5X7mnvJWnTpk2aNWuW0tPTTyhoo0eP9kowAABQsb5dXjy74sXN6ygogLNcAMAK5S5kH374oe69917FxMQoPt7zfHObzUYhAwCgCnC5jH5ctVeS1L8dk3IBgFXKXcief/55vfDCC3rssccqIg8AAKgEy3Ye0t7MPIU5/NXznDqnfwIAoEKU+z5khw4d0o033lgRWQAAQCX5YWXx7WuSz43ldEUAsFC5C9mNN96on3/+uSKyAACASuByGf20uvh0xSva1LU4DQDUbOU+ZbFp06Z68skntWDBArVp00YBAQEe6x944AGvhQMAAN7H6YoA4DvKXcg++OADhYWFac6cOZozZ47HOpvNRiEDAMDHTVlZfHTs0pZxnK4IABYr9ymLKSkpJ/3YunWr1wPu3r1bf/nLX1S7dm0FBwerTZs2WrJkiXu9MUajR49W3bp1FRwcrOTkZG3atMljGwcPHtRtt92miIgIRUVFafDgwTp8+LDHmJUrV6pHjx4KCgpSYmKixo4d6/V9AQDAasYYTV1dfP0YpysCgPXKXcgq06FDh3ThhRcqICBAP/30k9auXavXXntN0dF/3rxy7NixGjdunMaPH6+FCxcqNDRUffr0UV5ennvMbbfdpjVr1mj69OmaMmWK5s6dq7vvvtu9PisrS5dddpkaNmyopUuX6pVXXtHTTz+tDz74oFL3FwCAirZiV6b2ZuYpJNBPPZrFWB0HAGo8mzHGlPdJu3bt0nfffacdO3aooKDAY93rr7/utXAjR47UvHnz9Ouvv5a63hijhIQEPfzww/r73/8uScrMzFRcXJwmTJigAQMGaN26dWrZsqUWL16sTp06SZKmTp2qK664Qrt27VJCQoLee+89Pf7440pNTVVgYKD7tSdPnqz169eXKWtWVpYiIyOVmZmpiIgIL+w9AADe99JP6/T+nK3qd15dvXNrB6vjAEC1VJ5uUO4jZDNmzFDz5s313nvv6bXXXtOsWbP0ySef6OOPP9by5cvPNHOpvvvuO3Xq1Ek33nijYmNj1b59e3344Yfu9SkpKUpNTVVycrJ7WWRkpDp37qz58+dLkubPn6+oqCh3GZOk5ORk2e12LVy40D2mZ8+e7jImSX369NGGDRt06NChUrPl5+crKyvL4wMAAF83Z8M+SdLlreItTgIAkM6gkI0aNUp///vftWrVKgUFBenrr7/Wzp07ddFFF3n9/mRbt27Ve++9p2bNmmnatGm699579cADD+jTTz+VJKWmFp8DHxcX5/G8uLg497rU1FTFxsZ6rPf391etWrU8xpS2jWNf43gvvfSSIiMj3R+JiYlnubcAAFSstKw8rU/Nls0mdW1S2+o4AACdQSFbt26d7rjjDknFxebIkSMKCwvTs88+qzFjxng1nMvlUocOHfTiiy+qffv2uvvuuzVkyBCNHz/eq69zJkaNGqXMzEz3x86dO62OBADAKf2yLk2S1C4xSjFhDovTAACkMyhkoaGh7uvG6tatqy1btrjX7d+/33vJjm6/ZcuWHsvOPfdc7dixQ5IUH198ukVaWprHmLS0NPe6+Ph4paene6wvKirSwYMHPcaUto1jX+N4DodDERERHh8AAPiy6WuLf7clnxt3mpEAgMpS7kLWpUsX/fbbb5KkK664Qg8//LBeeOEF3XXXXerSpYtXw1144YXasGGDx7KNGzeqYcOGkqSkpCTFx8drxowZ7vVZWVlauHChunbtKknq2rWrMjIytHTpUveYmTNnyuVyqXPnzu4xc+fOVWFhoXvM9OnT1bx5c48ZHQEAqKoO5xfp980HJEmXtaSQAYCvKHche/31191F5plnnlHv3r315ZdfqlGjRvroo4+8Gu6hhx7SggUL9OKLL2rz5s36/PPP9cEHH2jo0KGSim9EPXz4cD3//PP67rvvtGrVKt1xxx1KSEhQ//79JRUfUbv88ss1ZMgQLVq0SPPmzdOwYcM0YMAAJSQkSJJuvfVWBQYGavDgwVqzZo2+/PJLvfnmmxoxYoRX9wcAAKv8unGfCpwuNawdoqaxYVbHAQAcdUbT3lemKVOmaNSoUdq0aZOSkpI0YsQIDRkyxL3eGKOnnnpKH3zwgTIyMtS9e3e9++67Ouecc9xjDh48qGHDhun777+X3W7X9ddfr3Hjxiks7M9fSCtXrtTQoUO1ePFixcTE6P7779djjz1W5pxMew8A8GUjvlyu/y3brb92T9ITV7Y8/RMAAGesPN2g3IWscePGWrx4sWrX9pydKSMjQx06dNDWrVvLn7gaoJABAHxVkdOlTi/8oozcQn15dxd1bswMiwBQkSr0PmTbtm2T0+k8YXl+fr52795d3s0BAIAKtmT7IWXkFioqJEAdG3JtNAD4Ev+yDvzuu+/cX0+bNk2RkZHux06nUzNmzFCjRo28Gg4AAJy9ktkVL2kRK3+/cv9bLACgApW5kJVMkmGz2TRw4ECPdQEBAWrUqJFee+01r4YDAABnxxjjvv8YsysCgO8pcyFzuVySiqeaL5n4AgAA+LZN6Ye1/UCuAv3t6tGsjtVxAADHKXMhK5GSklIROQAAQAWYsS5dktStSW2FOsr9ax8AUMHK/ZP52WefPeX60aNHn3EYAADgPcYY/bBqjyTponM4OgYAvqjcheybb77xeFxYWKiUlBT5+/urSZMmFDIAAHzE6t1ZWr07S4H+dvU7r67VcQAApSh3IVu2bNkJy7KysjRo0CBde+21XgkFAADO3ueLtksqnswjNjzI4jQAgNJ4Ze7biIgIPfPMM3ryySe9sTkAAHCWMnML9fXS4vuDDuzWyNowAICT8trNSDIzM5WZmemtzQEAgLMwY32aCpwuNYsN0/mNalkdBwBwEuU+ZXHcuHEej40x2rt3rz777DP17dvXa8EAAMCZ+2ZZ8dGxvq3jLU4CADiVcheyN954w+Ox3W5XnTp1NHDgQI0aNcprwQAAwJnZnXFEv27aL0m6oWOixWkAAKfCfcgAAKhmFmw5IElqmxilBrVDLE4DADgVr11DBgAAfMOvm/ZJki5sUtviJACA0ynzEbK77rqrTOM+/vjjMw4DAADOTqHTpTkbiwsZN4MGAN9X5kI2YcIENWzYUO3bt5cxpiIzAQCAM7R0+yEdyi1UdEiAOjSMtjoOAOA0ylzI7r33Xk2cOFEpKSm688479Ze//EW1ajGNLgAAvmT2hj+PjgX4cWUCAPi6Mv+kfuedd7R37149+uij+v7775WYmKibbrpJ06ZN44gZAAA+4rfNRwtZc05XBICqoFz/dOZwOHTLLbdo+vTpWrt2rVq1aqX77rtPjRo10uHDhysqIwAAKIPtB3K0Zk+WJOnCJjEWpwEAlMUZn8tgt9tls9lkjJHT6fRmJgAAcAZmrEuXMVLnpFqKjQiyOg4AoAzKVcjy8/M1ceJEXXrppTrnnHO0atUqvf3229qxY4fCwsIqKiMAACiDktkVe7WItTgJAKCsyjypx3333acvvvhCiYmJuuuuuzRx4kTFxHA6BAAAvqCgyKUFW4tvCN2rOYUMAKqKMhey8ePHq0GDBmrcuLHmzJmjOXPmlDruf//7n9fCAQCAspmzcZ/yi1yqHRqoc+I4awUAqooyF7I77rhDNputIrMAAIAzNGnJTknSdR3q8fsaAKqQct0YGgAA+J7D+UWauT5dknRDx0SL0wAAyoM7RgIAUMX9unGfilxG9aOD1Tw+3Oo4AIByoJABAFDFTV6+W5J0VdsEi5MAAMqLQgYAQBXmchktSjkoSUo+N87iNACA8qKQAQBQhS3bmaFDuYUKc/irTb1Iq+MAAMqpTIWsQ4cOOnTokCTp2WefVW5uboWGAgAAZfP75v2SpIvOqaNAf/6dFQCqmjL95F63bp1ycnIkSc8884wOHz5coaEAAEDZ/Hq0kF2QVMviJACAM1Gmae/btWunO++8U927d5cxRq+++qrCwkq/6eTo0aO9GhAAAJQuM7dQS7YVXz/W+9xYi9MAAM5EmQrZhAkT9NRTT2nKlCmy2Wz66aef5O9/4lNtNhuFDACASvLb5v1yGalZbJjqR4dYHQcAcAbKVMiaN2+uL774QpJkt9s1Y8YMxcbyL3EAAFhpzsbim0FfdE4di5MAAM5UmQrZsVwuV0XkAAAA5WCM0dyNxdeP9aSQAUCVVe5CJklbtmzRP//5T61bt06S1LJlSz344INq0qSJV8MBAIDSbUw7rNSsPAUF2JnQAwCqsHLPjztt2jS1bNlSixYt0nnnnafzzjtPCxcuVKtWrTR9+vSKyAgAAI4zd+M+SVLnpNoKCvCzOA0A4EyV+wjZyJEj9dBDD+nll18+Yfljjz2mSy+91GvhAABA6ZbvzJAkdWtS29ogAICzUu4jZOvWrdPgwYNPWH7XXXdp7dq1XgkFAABOzhjjLmTnxIVbGwYAcFbKXcjq1Kmj5cuXn7B8+fLlzLwIAEAl2Jh2WLszjigowK6uHCEDgCqt3KcsDhkyRHfffbe2bt2qbt26SZLmzZunMWPGaMSIEV4PCAAAPC1KOSBJ6tgwmuvHAKCKK3che/LJJxUeHq7XXntNo0aNkiQlJCTo6aef1gMPPOD1gAAAwNOM9cX3H+vWJMbiJACAs2UzxpgzfXJ2drYkKTyc89ezsrIUGRmpzMxMRUREWB0HAFBN5eQXqf2z01XgdGn6Qz3VjGvIAMDnlKcbnNF9yEpQxAAAqFzT16apwOlSg1ohahobZnUcAMBZKvekHgAAwDpfLdkpSbquQz3ZbDaL0wAAzhaFDACAKuLA4Xwt2Fo8ocd17etbnAYA4A0UMgAAqojpa9PkMlKrhAg1qB1idRwAgBeUq5AVFhaqd+/e2rRpU0XlAQAAJ/HT6lRJUt/W8RYnAQB4S7kKWUBAgFauXFlRWQAAwElkHinU71v2S5Iub13X4jQAAG8p9ymLf/nLX/TRRx9VRBYAAHASM9alqdBp1DQ2jNkVAaAaKfe090VFRfr444/1yy+/qGPHjgoNDfVY//rrr3stHAAAKDZ34z5JUp9WcRYnAQB4U7kL2erVq9WhQwdJ0saNGz3WMf0uAADed6TAqRnr0iVJPZrVsTgNAMCbyl3IZs2aVRE5AADASXy7fLey84uUWCtYFzSqZXUcAIAXnfG095s3b9a0adN05MgRSZIxxmuhAADAnxalHJQkXdO2nux2zkYBgOqk3IXswIED6t27t8455xxdccUV2rt3ryRp8ODBevjhh70eEACAmswYo/lHbwbduTFHxwCguil3IXvooYcUEBCgHTt2KCTkz5tS3nzzzZo6dapXwwEAUNNtO5CrvZl5CvSz63xOVwSAaqfc15D9/PPPmjZtmurXr++xvFmzZtq+fbvXggEAAGn2huLJPNolRikowM/iNAAAbyv3EbKcnByPI2MlDh48KIfD4ZVQAACg2Hcr9kiSLm8db3ESAEBFKHch69Gjh/7v//7P/dhms8nlcmns2LHq1auXV8MBAFCTbUzL1rIdGbLbpCvPq2t1HABABSj3KYtjx45V7969tWTJEhUUFOjRRx/VmjVrdPDgQc2bN68iMgIAUCNNOXp07JIWcYqNCLI4DQCgIpT7CFnr1q21ceNGde/eXddcc41ycnJ03XXXadmyZWrSpElFZAQAoEaafvRm0H05XREAqq1yHyGTpMjISD3++OPezgIAAI7akJqtdXuz5G+36aLmdayOAwCoIGdUyA4dOqSPPvpI69atkyS1bNlSd955p2rVYjpeAAC8YfLy3ZKknufUUUwYk2YBQHVV7lMW586dq0aNGmncuHE6dOiQDh06pHHjxikpKUlz586tiIwAANQ4S7YdlCRd2jLO4iQAgIpU7iNkQ4cO1c0336z33ntPfn7F90NxOp267777NHToUK1atcrrIQEAqElSM/O0ZPshSVKPZjEWpwEAVKRyHyHbvHmzHn74YXcZkyQ/Pz+NGDFCmzdv9mo4AABqou9X7JExUqeG0aoffeK9PwEA1Ue5C1mHDh3c144da926dWrbtq1XQgEAUFMZY/T1H7skSde0r2dxGgBARSvTKYsrV650f/3AAw/owQcf1ObNm9WlSxdJ0oIFC/TOO+/o5ZdfrpiUAADUEMt3Zmh9araCA/x0FTeDBoBqz2aMMacbZLfbZbPZdLqhNptNTqfTa+GqkqysLEVGRiozM1MRERFWxwEAVFFv/rJJb/yyUZe3itf42ztaHQcAcAbK0w3KdIQsJSXFK8EAAMCpTVm5R5LU+9xYi5MAACpDmQpZw4YNKzoHAAA13u6MI9qUflh+dpsuaxVvdRwAQCU4oxtD79mzR7/99pvS09Plcrk81j3wwANeCQYAQE0zd+M+SVLzuHBFBgdYnAYAUBnKXcgmTJige+65R4GBgapdu7ZsNpt7nc1mo5ABAHCG5mwoLmQ9z6ljcRIAQGUp97T3Tz75pEaPHq3MzExt27ZNKSkp7o+tW7dWREa3l19+WTabTcOHD3cvy8vL09ChQ1W7dm2FhYXp+uuvV1pamsfzduzYoX79+ikkJESxsbF65JFHVFRU5DFm9uzZ6tChgxwOh5o2baoJEyZU6L4AAHCsrLxCzTl6hKxfG2ZXBICaotyFLDc3VwMGDJDdXu6nnpXFixfr/fff13nnneex/KGHHtL333+vSZMmac6cOdqzZ4+uu+4693qn06l+/fqpoKBAv//+uz799FNNmDBBo0ePdo9JSUlRv3791KtXLy1fvlzDhw/XX//6V02bNq3S9g8AULP9unG/jhQ61TgmVK3rMVsvANQU5W5VgwcP1qRJkyoiy0kdPnxYt912mz788ENFR0e7l2dmZuqjjz7S66+/rksuuUQdO3bUJ598ot9//10LFiyQJP38889au3at/v3vf6tdu3bq27evnnvuOb3zzjsqKCiQJI0fP15JSUl67bXXdO6552rYsGG64YYb9MYbb1TqfgIAaq45G9MlSZe0iPW4HAAAUL2V+xqyl156SVdeeaWmTp2qNm3aKCDA86Lj119/3WvhSgwdOlT9+vVTcnKynn/+effypUuXqrCwUMnJye5lLVq0UIMGDTR//nx16dJF8+fPV5s2bRQXF+ce06dPH917771as2aN2rdvr/nz53tso2TMsadGHi8/P1/5+fnux1lZWV7YUwBATWSMcZ+ueFFzrh8DgJrkjArZtGnT1Lx5c0k6YVIPb/viiy/0xx9/aPHixSesS01NVWBgoKKiojyWx8XFKTU11T3m2DJWsr5k3anGZGVl6ciRIwoODj7htV966SU988wzZ7xfAACUWJ+arbSsfAUH+On8RrWsjgMAqETlLmSvvfaaPv74Yw0aNKgC4njauXOnHnzwQU2fPl1BQUEV/nrlMWrUKI0YMcL9OCsrS4mJiRYmAgBUVbOPzq7YpXEtBQX4WZwGAFCZyn0NmcPh0IUXXlgRWU6wdOlSpaenq0OHDvL395e/v7/mzJmjcePGyd/fX3FxcSooKFBGRobH89LS0hQfX3xDzfj4+BNmXSx5fLoxERERpR4dk4q/DxERER4fAACciaXbD0mSujfjdEUAqGnKXcgefPBBvfXWWxWR5QS9e/fWqlWrtHz5cvdHp06ddNttt7m/DggI0IwZM9zP2bBhg3bs2KGuXbtKkrp27apVq1YpPT3dPWb69OmKiIhQy5Yt3WOO3UbJmJJtAABQUZwuoyXbD0qS2taPtDgNAKCylfuUxUWLFmnmzJmaMmWKWrVqdcKkHv/73/+8Fi48PFytW7f2WBYaGqratWu7lw8ePFgjRoxQrVq1FBERofvvv19du3ZVly5dJEmXXXaZWrZsqdtvv11jx45VamqqnnjiCQ0dOlQOh0OS9Le//U1vv/22Hn30Ud11112aOXOmvvrqK/3www9e2xcAAEqzcleGMnILFR7kr3aJUVbHAQBUsnIXsqioKI/7fFntjTfekN1u1/XXX6/8/Hz16dNH7777rnu9n5+fpkyZonvvvVddu3ZVaGioBg4cqGeffdY9JikpST/88IMeeughvfnmm6pfv77+9a9/qU+fPlbsEgCgBlm+M0OS1LFhtPz9KvcenwAA69mMMcbqENVBVlaWIiMjlZmZyfVkAIAyG/jxIs3ZuE8j+7bQ3y5qYnUcAIAXlKcb8E9xAABYJK/QqcXbiq8f68mEHgBQI5X7lMWkpKRT3m9s69atZxUIAICaYvraNOUWOFU/Olgt4sOtjgMAsEC5C9nw4cM9HhcWFmrZsmWaOnWqHnnkEW/lAgCg2pu6JlWSdFXbBNntJ//HTgBA9VXuQvbggw+Wuvydd97RkiVLzjoQAAA1QX6RU7PXF9+SpU+reIvTAACs4rVryPr27auvv/7aW5sDAKBa+33LAeUUOBUX4dB59bj/GADUVF4rZP/9739Vq1Ytb20OAIBqbfraNEnSpS3jOF0RAGqwcp+y2L59e49JPYwxSk1N1b59+zzu/wUAAErncpljChmnKwJATVbuQta/f3+Px3a7XXXq1NHFF1+sFi1aeCsXAADV1vJdGdqXna9wh7+6Nq5tdRwAgIXKXcieeuqpisgBAECN8fOa4qNjF7eIVaA/twQFgJqM3wIAAFQiY4x+XLVXUvH1YwCAmq3MR8jsdvspbwgtSTabTUVFRWcdCgCA6mp9arZ2HMyVw9+u5HNjrY4DALBYmQvZN998c9J18+fP17hx4+RyubwSCgCA6uqn1cU3g+7RrI5CAst95QAAoJop82+Ca6655oRlGzZs0MiRI/X999/rtttu07PPPuvVcAAAVDeLUw5Kki5pwdExAMAZXkO2Z88eDRkyRG3atFFRUZGWL1+uTz/9VA0bNvR2PgAAqo1DOQVauuOQJOn8RtEWpwEA+IJyFbLMzEw99thjatq0qdasWaMZM2bo+++/V+vWrSsqHwAA1cbXf+xSQZFLrRIi1DQ2zOo4AAAfUOZTFseOHasxY8YoPj5eEydOLPUURgAAcHLT1hRfP3bz+YmnnSgLAFAz2IwxpiwD7Xa7goODlZycLD8/v5OO+9///ue1cFVJVlaWIiMjlZmZqYiICKvjAAB8TOaRQnV4brqcLqNfH+2lxFohVkcCAFSQ8nSDMh8hu+OOO/jXPAAAztAva9PkdBklxYRSxgAAbmUuZBMmTKjAGAAAVG8/ry0+XfHK8+panAQA4EvOaJZFAABQdrkFRfpt035JUvK5cRanAQD4EgoZAAAV7MdVqcopcKpR7RCdVz/S6jgAAB9CIQMAoIJNWrJTknRDx/pcjw0A8EAhAwCgAu04kKuFKQdls0nXdahvdRwAgI+hkAEAUIH+u7T46Fj3pjFKiAq2OA0AwNdQyAAAqCDGGP1v2W5J0o2dEi1OAwDwRRQyAAAqyLq92dp16IgC/ey6lNkVAQCloJABAFBBpqzcI0m6uHkdBQf6WZwGAOCLKGQAAFSAgiKXJh89XfGqtgkWpwEA+CoKGQAAFeDLxTu0JzNPseEOXdqS0xUBAKWjkAEA4GXGGP17wQ5J0r0XN1FQAKcrAgBKRyEDAMDLVu7K1Ia0bDn87bquPfceAwCcHIUMAAAv+98fuyRJfVrFKzIkwOI0AABfRiEDAMCLXC6jX9alS5KuacdkHgCAU6OQAQDgRQtSDmh3xhGFOfzVrUmM1XEAAD6OQgYAgBf9748/p7rn3mMAgNOhkAEA4CV5hU5NW50qidMVAQBlQyEDAMBLZqxLV3Z+kepFBeuCRrWsjgMAqAIoZAAAeMk3y4pPV7ymXYLsdpvFaQAAVQGFDAAALziYU6DZG4pnV7y2fT2L0wAAqgoKGQAAXvDDyj0qchm1SohQs7hwq+MAAKoIChkAAF7wxeKdkqTrOtS3OAkAoCqhkAEAcJY2p2drzZ4s+dltnK4IACgXChkAAGdpzsb9kqRuTWqrVmigxWkAAFUJhQwAgLO0ZNtBSdKFTWMsTgIAqGooZAAAnIXtB3L0y7o0SdIFSdx7DABQPhQyAADOwgs/rFOh06hHsxh1aBBtdRwAQBVDIQMA4Ayt3JWhn9emKcDPplF9z7U6DgCgCqKQAQBwhr48OtX9FW3qqmVChMVpAABVEYUMAIAzkJNfpG+X75Ek3dQp0eI0AICqikIGAMAZ+H7FHh3OL1Kj2iHq2ri21XEAAFUUhQwAgDPw+aIdkqRbLmggu91mcRoAQFVFIQMAoJxW7srQyl2ZCvSz64aO9a2OAwCowihkAACU02fzt0uS+raJV+0wh8VpAABVGYUMAIBySMvK0+TluyVJd3RtZG0YAECVRyEDAKAcPp6XokKn0QWNaqljQ24EDQA4OxQyAADKKDuvUJ8vKJ7M4+6ejS1OAwCoDihkAACU0cRFO5SdX6SmsWG6pEWs1XEAANUAhQwAgDJwuow+/m2bJGlIjySmugcAeAWFDACAMvht836lZuUpOiRA17SrZ3UcAEA1QSEDAKAMPv19mySp33l1FRTgZ20YAEC1QSEDAOA0dmcc0awN6ZKkQd0aWRsGAFCtUMgAADiNb/7YJWOkLo1rqWlsuNVxAADVCIUMAIBTMMbo6z+KbwR9fYf6FqcBAFQ3FDIAAE7hjx0ZStmfo+AAP/VtU9fqOACAaoZCBgDAKXz9xy5JUt/W8Qpz+FucBgBQ3VDIAAA4iYzcAn1TcrpiR05XBAB4H4UMAICT+GTeNh0pdKpxnVB1aVzb6jgAgGqIQgYAQCmcLqNJS3ZKkob1aio/u83iRACA6ohCBgBAKWasS9OezDxFBgfoCibzAABUEAoZAAClGD9niyTp1s4NFBTgZ3EaAEB1RSEDAOA4a/Zk6o8dGQrws+muC5OsjgMAqMYoZAAAHOfj37ZJki5rGa864Q5rwwAAqjUKGQAAx9iXna/vVhRPdX9blwYWpwEAVHc+XcheeuklnX/++QoPD1dsbKz69++vDRs2eIzJy8vT0KFDVbt2bYWFhen6669XWlqax5gdO3aoX79+CgkJUWxsrB555BEVFRV5jJk9e7Y6dOggh8Ohpk2basKECRW9ewAAH/Tl4h0qdBq1TYxSV6a6BwBUMJ8uZHPmzNHQoUO1YMECTZ8+XYWFhbrsssuUk5PjHvPQQw/p+++/16RJkzRnzhzt2bNH1113nXu90+lUv379VFBQoN9//12ffvqpJkyYoNGjR7vHpKSkqF+/furVq5eWL1+u4cOH669//aumTZtWqfsLALDWvux8fTp/uyTplvMTZbMx1T0AoGLZjDHG6hBltW/fPsXGxmrOnDnq2bOnMjMzVadOHX3++ee64YYbJEnr16/Xueeeq/nz56tLly766aefdOWVV2rPnj2Ki4uTJI0fP16PPfaY9u3bp8DAQD322GP64YcftHr1avdrDRgwQBkZGZo6dWqZsmVlZSkyMlKZmZmKiIjw/s4DACrc09+t0YTft6lxnVD9cH8PBQcyuyIAoPzK0w18+gjZ8TIzMyVJtWrVkiQtXbpUhYWFSk5Odo9p0aKFGjRooPnz50uS5s+frzZt2rjLmCT16dNHWVlZWrNmjXvMsdsoGVOyjdLk5+crKyvL4wMAUHWlZ+fpy8XFN4J++qpWlDEAQKWoMoXM5XJp+PDhuvDCC9W6dWtJUmpqqgIDAxUVFeUxNi4uTqmpqe4xx5axkvUl6041JisrS0eOHCk1z0svvaTIyEj3R2Ji4lnvIwDAOp/+vk1HCp1qWz9SPZrFWB0HAFBDVJlCNnToUK1evVpffPGF1VEkSaNGjVJmZqb7Y+fOnVZHAgCcodyCIn0yb5sk6Z6LmnDtGACg0vhbHaAshg0bpilTpmju3LmqX7++e3l8fLwKCgqUkZHhcZQsLS1N8fHx7jGLFi3y2F7JLIzHjjl+Zsa0tDRFREQoODi41EwOh0MOB/emAYDq4N1ZW5Rb4FTD2iHq2zre6jgAgBrEp4+QGWM0bNgwffPNN5o5c6aSkpI81nfs2FEBAQGaMWOGe9mGDRu0Y8cOde3aVZLUtWtXrVq1Sunp6e4x06dPV0REhFq2bOkec+w2SsaUbAMAUH3ty87XfxYWz6x4T0+OjgEAKpdPHyEbOnSoPv/8c3377bcKDw93X/MVGRmp4OBgRUZGavDgwRoxYoRq1aqliIgI3X///eratau6dOkiSbrsssvUsmVL3X777Ro7dqxSU1P1xBNPaOjQoe4jXH/729/09ttv69FHH9Vdd92lmTNn6quvvtIPP/xg2b4DACqeMUYPfrFMh3IL1TgmVDd1qn/6JwEA4EU+Pe39yf6V8pNPPtGgQYMkFd8Y+uGHH9bEiROVn5+vPn366N1333WfjihJ27dv17333qvZs2crNDRUAwcO1Msvvyx//z/76OzZs/XQQw9p7dq1ql+/vp588kn3a5QF094DQNXz66Z9uv2j4tPa/3dfN3VoEG1xIgBAdVCebuDThawqoZABQNVijNHN7y/Qom0HdeeFjfTUVa2sjgQAqCaq7X3IAADwliXbD2nRtoMK9LPr7p6NrY4DAKihKGQAgBpp/OwtkqRr29dT3cjSZ9QFAKCiUcgAADXOkm0HNWN9uvzsNt1zEUfHAADWoZABAGoUY4zGTt0gSbqxY301rhNmcSIAQE1GIQMA1CizN+4rvnbM364Hk5tZHQcAUMNRyAAANUah06Wnvl0jSRrYtSHXjgEALEchAwDUGF8v3aUdB3MVExaoB3pzdAwAYD0KGQCgRjiUU6AxU9dLkv52UROFBwVYnAgAAAoZAKCGGDttvQ7lFqp5XLgGdmtkdRwAACRRyAAANcDCrQc0cdFOSdLz17ZWgB+//gAAvoHfSACAai2/yKl/fLNKknRdh3o6v1EtixMBAPAnChkAoFobO3WDtuzLUe3QQD3Zr6XVcQAA8EAhAwBUW7+sTdNHv6VIkp7r31rRoYEWJwIAwBOFDABQLRljNG7mJknSXRcm6Yo2dS1OBADAiShkAIBqadaGdK3clalAf7uG9mpidRwAAEpFIQMAVDt5hU49/8M6SdLtXRqqdpjD4kQAAJSOQgYAqFbyCp0a9vkf2np0Io9hvZpaHQkAgJOikAEAqo2CIpfun7hMv6xLV6CfXW/d0p6JPAAAPo1CBgCoNkZ+vVLT16bJZpM+uKOjujWNsToSAACnRCEDAFQL//p1q/63bLck6c0B7XVx81iLEwEAcHoUMgBAlTdn4z73JB73XtxEV7dNsDgRAABlQyEDAFRpLpfRmJ/WS5KubV9Pj/ZpbnEiAADKjkIGAKjSPl+0Q2v3Zikk0E+jr2wpm81mdSQAAMqMQgYAqLIycwv1z182SpJGXHoOMyoCAKocChkAoMp67oe12n+4QI1qh2hgt0ZWxwEAoNwoZACAKmnCvBT9d+kuSdLL15+nAD9+pQEAqh5+ewEAqpzF2w66Z1V8sHczdWlc2+JEAACcGQoZAKBK2X84X3f/3xIVuYyuPK+uhic3szoSAABnjEIGAKgy8gqdevCLZTqUW6gW8eEae8N5zKoIAKjSKGQAgCrjpR/Xad7mA3L42/XqjW0VEuhvdSQAAM4KhQwAUCV8u3y3Pp2/XZL09q0d1LpepMWJAAA4exQyAIDPW5+apUf/u1KSdGvnBrq0ZZzFiQAA8A4KGQDAp2XnFeq+//yh/CKXejSL0XPXtLY6EgAAXkMhAwD4LGOMhn6+TFv35Sg23KHXbmwrPzuTeAAAqg8KGQDAJxljNPrbNZq7cZ+CAuz6aOD5io0IsjoWAABeRSEDAPikj35L0WcLtstmk57v30Zt6jOJBwCg+qGQAQB8zk+r9uqFH9dJkh6/4lzd0LG+xYkAAKgY3MAFAOBT3pu9RWOnrZcxxTMqDu6eZHUkAAAqDIUMAOATnC6jp75brX8v2CFJurlTop69upVsNibxAABUXxQyAIDltu47rBFfrdDynRmSpKG9muiRPi2sDQUAQCWgkAEALLVk20EN/nSJMo8UyuFv1+irWurWCxpYHQsAgEpBIQMAWGb2hnQN/3K5Mo8Uql1ilP55czs1igm1OhYAAJWGQgYAsMQ7szbrlWkbJElt60dq4pAuCg70szgVAACVi0IGAKh0b/6ySW/8slGSdGPH+nr2mtaUMQBAjUQhAwBUmvSsPD0+ebWmr02TJN1/SVONuPQcZlIEANRYFDIAQKXYeTBXAz5YoN0ZRyQV3/B5SM/GFqcCAMBaFDIAQIXblJate//zh3ZnHFGDWiF67y8d1Coh0upYAABYjkIGAKhQ3y7frZFfr9KRQqfiIhz6z187K7FWiNWxAADwCRQyAECFOHA4X+/O3qKPfkuRJPVoFqPXb2qnOuEOi5MBAOA7KGQAAK/KK3Tq3Vmb9cGvW5VX6JIk3XVhkp7od67sdibvAADgWBQyAIDXzN6Qricmr9auQ8UTdzSNDdM/rmihS1rEWZwMAADfRCEDAJy13IIivTptoz6eV3x6YnxEkEZd0UJXt01gSnsAAE6BQgYAOCtzNu7TqK9Xak9mniTplgsS9Xi/lgpz8CsGAIDT4bclAOCM5BU69eKP6/R/87dLkhIig/TsNa2V3JLTEwEAKCsKGQCgXIwxmrYmTa/+vEGb0w9Lkm7r3ECP9ztXIYH8WgEAoDz4zQkAKLPM3EI9PnmVpqzcK0mKCPLXuFva6+LmsRYnAwCgaqKQAQBO60iBU+/M2qz/Lt2l1Kw82WzS3y5qont6NlZUSKDV8QAAqLIoZACAU9q677DumrBY2w7kSpIa1ArRGze3VceGtSxOBgBA1UchAwCUakNqtj5bsE1fLdmlgiKX6oQ7NKRHkm7t3JAZFAEA8BJ+owIA3Fwuo/fnbtVHv23V/sMF7uVdGtfSW7d0UJ1wh4XpAACofihkAAC5XEazNqTr9ekbtWZPlnv5BUm1dH2HerqxY6Lsdm7wDACAt1HIAKAGS8/K089r0/TurM3uGzuHOfz10KXn6JYLEpnGHgCACsZvWgCoQYwx2n4gVzPWp+vfC7YrZX+Oe12Yw183dKyvB3o3U61QZk4EAKAyUMgAoJozxuiPHRn6cvEOzVyf7nFtmM0mnRMbrj6t43V3z8ZM1gEAQCXjNy8AVFO7M45o5ro0fTp/uzanH3YvD/Szq3W9CF3bvp6ubldPkcEBFqYEAKBmo5ABQDVyOL9Iv2/erw9/3arF2w65lwcF2HVFm7q6oUN9dWwULYe/n4UpAQBACQoZAFRhBw7na1HKQa3YlallOw5p6fZDKnIZ9/q2iVG6pHmsBnVrpMgQjoQBAOBrKGQAUIVk5BZoY9phrd2TqR9W7fU4ClaifnSwLmkRq7suTFKjmFALUgIAgLKikAGAj8ordGpDarZW7MrQkm2HtHp3prYeMytiiSZ1QtWhQbTaJkapS+NaahobbkFaAABwJihkAOAjjhQ4tWJXhn5YuVcb0rK1KOVgqePqRQWreXy4WsSH67YuDVUvKriSkwIAAG+hkAFAJXO6jFL2H9bWfTnalH5Y61OztSE1SxvTDpc6vn2DKPVsVkdtEyPVpl6U6oQ7KjkxAACoKBQyAPAip8sot6BIRwqcyi1wKjuvSNsP5ihlX45S9udo877DWrkr86TPrxPuUMcG0WqVEKErzqurxjGhstlslbgHAACgMlHIAKAcCp0uHThcoH3Z+dp/OF+7DuVqy74cbdlXfMRrT+YRGXP67UhS63oRahwTphZ1w3VufIRa1A1XfEQQBQwAgBqEQnacd955R6+88opSU1PVtm1bvfXWW7rgggusjgXAi4wxOlLoVE6+U7kFRdp/uEDpWXnKzitSdn6Rso4UFh/lKnTqSIFLRwqLdCinUAdy8rV1X47HtPInY7NJoYH+CnX4KTE6REkxoUqqE6rGMaFKiglT4zqhCvCzV8LeAgAAX0YhO8aXX36pESNGaPz48ercubP++c9/qk+fPtqwYYNiY2OtjgdUa8YYFbmMnEc/ilxGhU6XDuUU6FBuoQqdrqPrXSpy/jnG6TLKK3QqNStPaVl5ysl3Kr/Iqfwil/ILXcovciqv0KXcgiLlFDiVm1+k3EJnmY9ilcbPblNMWKDqhDsUHxGkJnWKC1bjOmFqWCtEEcEBcvjbOdIFAABOy2bM2fxZUr107txZ559/vt5++21JksvlUmJiou6//36NHDnylM/NyspSZGSkMjMzFRERURlxT2rZjkNKy8pz/8FZ8h/4z8fmuMd/Ov7tcNLnnLBt47mtYzZzstc7bZ7jtnlstONf72R5jt+PU75eKZmOf+4Jr++x3dNvw2NTx3/Pyvjc0vZFpeY49euf7PvvGdEot8DpUX6KXC45XSouRseUJ5fHmJLHLo/nHlu2PB+7VIaDThUiNNBPkcEBSogKVniQv8KCAoo/O/wVFOCnkEA/BQf4KSokQNEhgUqIClLjmDDZ7ZQtAABQuvJ0A46QHVVQUKClS5dq1KhR7mV2u13JycmaP3/+CePz8/OVn5/vfpyVlVUpOcviX7+m6IdVe62OAXhFSREK9LPLz26Tv5+t+LO95LNdgf52xYY7FB8ZpPCg4qNTDn+7HAF+cvjbFRTgp9BAP4U6/BUa6K/gQD+FOvwU5O9HsQIAAJaikB21f/9+OZ1OxcXFeSyPi4vT+vXrTxj/0ksv6ZlnnqmseOXSpE6ozm8ULUmy6egfm56fZHM/tnk8Pt264se2k2yr9PWe27Cd5Dmlv+bxy4/d6Mleryx5dKr9PyHzqb9Px275ZPt1+tc6sRSc6rmlfW+P3UYp0cq0r7ZS9iUowE+BfnbZPUpQ8We/Y0qRn/34onTyMR7j/Gzysx0zxq94nd129DOFCQAAVGMUsjM0atQojRgxwv04KytLiYmJFib604jLmlsdAQAAAEAZUMiOiomJkZ+fn9LS0jyWp6WlKT4+/oTxDodDDgc3ZwUAAABw5phz+ajAwEB17NhRM2bMcC9zuVyaMWOGunbtamEyAAAAANUVR8iOMWLECA0cOFCdOnXSBRdcoH/+85/KycnRnXfeaXU0AAAAANUQhewYN998s/bt26fRo0crNTVV7dq109SpU0+Y6AMAAAAAvIH7kHmJL92HDAAAAIB1ytMNuIYMAAAAACxCIQMAAAAAi1DIAAAAAMAiFDIAAAAAsAiFDAAAAAAsQiEDAAAAAItQyAAAAADAIhQyAAAAALAIhQwAAAAALEIhAwAAAACLUMgAAAAAwCIUMgAAAACwCIUMAAAAACzib3WA6sIYI0nKysqyOAkAAAAAK5V0gpKOcCoUMi/Jzs6WJCUmJlqcBAAAAIAvyM7OVmRk5CnH2ExZahtOy+Vyac+ePQoPD5fNZrM0S1ZWlhITE7Vz505FRERYmgU4Fu9N+Crem/BFvC/hq3hvnp4xRtnZ2UpISJDdfuqrxDhC5iV2u13169e3OoaHiIgI/ieBT+K9CV/FexO+iPclfBXvzVM73ZGxEkzqAQAAAAAWoZABAAAAgEUoZNWQw+HQU089JYfDYXUUwAPvTfgq3pvwRbwv4at4b3oXk3oAAAAAgEU4QgYAAAAAFqGQAQAAAIBFKGQAAAAAYBEKGQAAAABYhELmA9555x01atRIQUFB6ty5sxYtWnTK8ZMmTVKLFi0UFBSkNm3a6Mcff/RYb4zR6NGjVbduXQUHBys5OVmbNm3yGHPw4EHddtttioiIUFRUlAYPHqzDhw97jFm5cqV69OihoKAgJSYmauzYsd7ZYVQZvvje3LZtm2w22wkfCxYs8N6Ow6dZ8b584YUX1K1bN4WEhCgqKqrU19mxY4f69eunkJAQxcbG6pFHHlFRUdFZ7SuqFl99b5b2M/OLL744q31F1VLZ781t27Zp8ODBSkpKUnBwsJo0aaKnnnpKBQUFHtvhb82jDCz1xRdfmMDAQPPxxx+bNWvWmCFDhpioqCiTlpZW6vh58+YZPz8/M3bsWLN27VrzxBNPmICAALNq1Sr3mJdfftlERkaayZMnmxUrVpirr77aJCUlmSNHjrjHXH755aZt27ZmwYIF5tdffzVNmzY1t9xyi3t9ZmamiYuLM7fddptZvXq1mThxogkODjbvv/9+xX0z4FN89b2ZkpJiJJlffvnF7N271/1RUFBQcd8M+Ayr3pejR482r7/+uhkxYoSJjIw84XWKiopM69atTXJyslm2bJn58ccfTUxMjBk1apTXvwfwTb763jTGGEnmk08+8fiZeew2UL1Z8d786aefzKBBg8y0adPMli1bzLfffmtiY2PNww8/7N4Gf2v+iUJmsQsuuMAMHTrU/djpdJqEhATz0ksvlTr+pptuMv369fNY1rlzZ3PPPfcYY4xxuVwmPj7evPLKK+71GRkZxuFwmIkTJxpjjFm7dq2RZBYvXuwe89NPPxmbzWZ2795tjDHm3XffNdHR0SY/P9895rHHHjPNmzc/yz1GVeGr782SQrZs2TKv7CeqFivel8f65JNPSv2j98cffzR2u92kpqa6l7333nsmIiLC4+coqi9ffW8aU1zIvvnmm3LuEaoLq9+bJcaOHWuSkpLcj/lb80+csmihgoICLV26VMnJye5ldrtdycnJmj9/fqnPmT9/vsd4SerTp497fEpKilJTUz3GREZGqnPnzu4x8+fPV1RUlDp16uQek5ycLLvdroULF7rH9OzZU4GBgR6vs2HDBh06dOgs9xy+zpffmyWuvvpqxcbGqnv37vruu+/ObodRJVj1viyL+fPnq02bNoqLi/N4naysLK1Zs6bM20HV5MvvzRJDhw5VTEyMLrjgAn388ccy3Ia2RvCl92ZmZqZq1arl8Tr8rVmMQmah/fv3y+l0evwCl6S4uDilpqaW+pzU1NRTji/5fLoxsbGxHuv9/f1Vq1YtjzGlbePY10D15cvvzbCwML322muaNGmSfvjhB3Xv3l39+/enlNUAVr0vy4KfmTWbL783JenZZ5/VV199penTp+v666/Xfffdp7feeqtc20DV5Cvvzc2bN+utt97SPffcc9rXOfY1agp/qwMAQHnExMRoxIgR7sfnn3++9uzZo1deeUVXX321hckAwDc9+eST7q/bt2+vnJwcvfLKK3rggQcsTIWaYvfu3br88st14403asiQIVbH8UkcIbNQTEyM/Pz8lJaW5rE8LS1N8fHxpT4nPj7+lONLPp9uTHp6usf6oqIiHTx40GNMads49jVQffnye7M0nTt31ubNm8uwZ6jKrHpflgU/M2s2X35vlqZz587atWuX8vPzz2o78H1Wvzf37NmjXr16qVu3bvrggw/K9DrHvkZNQSGzUGBgoDp27KgZM2a4l7lcLs2YMUNdu3Yt9Tldu3b1GC9J06dPd49PSkpSfHy8x5isrCwtXLjQPaZr167KyMjQ0qVL3WNmzpwpl8ulzp07u8fMnTtXhYWFHq/TvHlzRUdHn+Wew9f58nuzNMuXL1fdunXLv6OoUqx6X5ZF165dtWrVKo9/UJg+fboiIiLUsmXLMm8HVZMvvzdLs3z5ckVHR8vhcJzVduD7rHxv7t69WxdffLE6duyoTz75RHa7Z+3gb81jWD2rSE33xRdfGIfDYSZMmGDWrl1r7r77bhMVFeWeqev22283I0eOdI+fN2+e8ff3N6+++qpZt26deeqpp0qdijQqKsp8++23ZuXKleaaa64pdWrx9u3bm4ULF5rffvvNNGvWzGNq8YyMDBMXF2duv/12s3r1avPFF1+YkJCQGjkVaU3lq+/NCRMmmM8//9ysW7fOrFu3zrzwwgvGbrebjz/+uBK+K7CaVe/L7du3m2XLlplnnnnGhIWFmWXLlplly5aZ7OxsY8yf095fdtllZvny5Wbq1KmmTp06THtfg/jqe/O7774zH374oVm1apXZtGmTeffdd01ISIgZPXp0JX1nYDUr3pu7du0yTZs2Nb179za7du3yuOVCCf7W/BOFzAe89dZbpkGDBiYwMNBccMEFZsGCBe51F110kRk4cKDH+K+++sqcc845JjAw0LRq1cr88MMPHutdLpd58sknTVxcnHE4HKZ3795mw4YNHmMOHDhgbrnlFhMWFmYiIiLMnXfe6f7hXWLFihWme/fuxuFwmHr16pmXX37ZuzsOn+eL780JEyaYc88914SEhJiIiAhzwQUXmEmTJnl/5+GzrHhfDhw40Eg64WPWrFnuMdu2bTN9+/Y1wcHBJiYmxjz88MOmsLDQ6/sP3+WL782ffvrJtGvXzoSFhZnQ0FDTtm1bM378eON0OivkewDfVNnvzU8++aTU9+Xxx4L4W7OYzRjmPQUAAAAAK3ANGQAAAABYhEIGAAAAABahkAEAAACARShkAAAAAGARChkAAAAAWIRCBgAAAAAWoZABAAAAgEUoZAAAAABgEQoZAADlNGjQIPXv39/qGACAasDf6gAAAPgSm812yvVPPfWU3nzzTRljKikRAKA6o5ABAHCMvXv3ur/+8ssvNXr0aG3YsMG9LCwsTGFhYVZEAwBUQ5yyCADAMeLj490fkZGRstlsHsvCwsJOOGXx4osv1v3336/hw4crOjpacXFx+vDDD5WTk6M777xT4eHhatq0qX766SeP11q9erX69u2rsLAwxcXF6fbbb9f+/fsreY8BAFaikAEA4AWffvqpYmJitGjRIt1///269957deONN6pbt276448/dNlll+n2229Xbm6uJCkjI0OXXHKJ2rdvryVLlmjq1KlKS0vTTTfdZPGeAAAqE4UMAAAvaNu2rZ544gk1a9ZMo0aNUlBQkGJiYjRkyBA1a9ZMo0eP1oEDB7Ry5UpJ0ttvv6327dvrxRdfVIsWLdS+fXt9/PHHmjVrljZu3Gjx3gAAKgvXkAEA4AXnnXee+2s/Pz/Vrl1bbdq0cS+Li4uTJKWnp0uSVqxYoVmzZpV6PdqWLVt0zjnnVHBiAIAvoJABAOAFAQEBHo9tNpvHspLZG10ulyTp8OHDuuqqqzRmzJgTtlW3bt0KTAoA8CUUMgAALNChQwd9/fXXatSokfz9+XUMADUV15ABAGCBoUOH6uDBg7rlllu0ePFibdmyRdOmTdOdd94pp9NpdTwAQCWhkAEAYIGEhATNmzdPTqdTl112mdq0aaPhw4crKipKdju/ngGgprAZY4zVIQAAAACgJuKf4AAAAADAIhQyAAAAALAIhQwAAAAALEIhAwAAAACLUMgAAAAAwCIUMgAAAACwCIUMAAAAACxCIQMAAAAAi1DIAAAAAMAiFDIAAAAAsAiFDAAAAAAs8v+jxNOb8/2DpQAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[(9999,), (1,), (10000,), (0.0,), (1,), (10000,), (4.771818163941268e-05,), (0,), (0,), (0,), (0.0001269114533750033,), (0,), (0,), (0,), (0.00020479840692238415,), (0,), (0,), (0,), (0.0003868512648608837,), (0,), (0,), (0,), (0.00040125281940051754,), (0,), (0,), (0,), (0.0004396525870125214,), (0,), (0,), (0,), (0.0004409992368832654,), (0,), (0,), (0,), (0.0004516841209499256,), (0,), (0,), (0,), (0.00046635198488061293,), (0,), (0,), (0,), (0.0004817560144229769,), (0,), (0,), (0,), (0.0004936453628786726,), (0,), (0,), (0,), (0.000508589160153316,), (0,), (0,), (0,), (0.0005171767755911625,), (0,), (0,), (0,), (0.0005255528771130577,), (0,), (0,), (0,), (0.0005281068653429891,), (0,), (0,), (0,), (0.0005312812683494216,), (0,), (0,), (0,), (0.0005335997649747014,), (0,), (0,), (0,), (0.0005346801579408522,), (0,), (0,), (0,), (0.000537125171186308,), (0,), (0,), (0,), (0.0005380844238624733,), (0,), (0,), (0,), (0.0005404373053466639,), (0,), (0,), (0,), (0.0005405694505743195,), (0,), (0,), (0,), (0.000561812604564046,), (0,), (0,), (0,), (0.0005633734852782369,), (0,)]\n" + ] + } + ], + "source": [ + "# Example usage\n", + "L = 1\n", + "N = 10000\n", + "s = 10000\n", + "m = 0\n", + "migration_matrix = MigrationMatrixGenerator(L, m).generate_1Dmatrix()\n", + "\n", + "# migration_matrix = np.zeros((L, L)) # Example migration matrix\n", + "\n", + "tfinal = 5000\n", + "simulation = GillespieSelectionSimulation(L, N, s, migration_matrix, tfinal)\n", + "events, mut_traj, start_state, end_state = simulation.simulate()\n", + "plot_mutant_trajectories(mut_traj, L)\n", + "simulation.save_events(events, start_state, end_state, 'events.bin')\n", + "loaded_events = simulation.load_events('events.bin')\n", + "\n", + "print(loaded_events[0:100]) # Optional: Print to verify contents\n", + "\n", + "# print(mut_traj)\n", + "#for event in events[0:10]:\n", + "# print(event) # Print each event\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "msprime_stable", + "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.11.13" + }, + "orig_nbformat": 4 + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/notebooks/events.bin b/notebooks/events.bin new file mode 100644 index 000000000..107739bf5 Binary files /dev/null and b/notebooks/events.bin differ