return bind;
}
-//Constructor that assembles a population of N viruses of wild type
-
-#if 0
-Population::Population(const Hamiltonian &H, unsigned int N, double mu) {
-
- Virus V(H, mu);
- pop[V] = N;
- this->N = N;
- this->mu = mu;
-
- p = 1.0-pow((1.0-V.mu),H.size);
- Eavg = V.energy;
-
-}
-
-
-//Constructor that assembles a population of N viruses given starting population fractions
-
-Population::Population(const Hamiltonian &H, unsigned int N, double mu, const std::vector<std::set<unsigned int> > &initPop, const std::vector<double> &initFrac) {
-
- for (unsigned i=0;i<initPop.size();i++) {
-
- Virus V(H, mu, initPop[i]);
- pop[V] = (unsigned int) N * initFrac[i];
- Eavg += V.energy * pop[V];
-
- }
-
- this->N = N;
- this->mu = mu;
-
- p = 1.0-pow((1.0-mu),H.size);
- Eavg /= (double) N;
-
-
- // if not given any particular thing
- // then initialize to N wild type
- if (initPop.size() == 0) {
- Virus V(H, mu);
- pop[V] = N;
- Eavg = V.energy;
- }
-
-}
-
-
-//Step population forward a generation.
-
-void Population::next_generation(const Hamiltonian &H, gsl_rng* r, bool useRelative, bool useVerbose) {
-
- mutate_population(H,r);
-
- unsigned int total_deaths = 0;
- unsigned int num_survive = 0;
- double new_Eavg = 0.0;
-
- virus_map::iterator iter=pop.begin();
-
- for(;iter!=pop.end();++iter){
-
- if (useRelative) num_survive = gsl_ran_binomial(r,iter->first.survival(Eavg),iter->second);
- else num_survive = gsl_ran_binomial(r,iter->first.survival(),iter->second);
-
- total_deaths += iter->second - num_survive;
- iter->second = num_survive;
-
- new_Eavg += iter->first.energy * num_survive;
-
- // Report survival
-
- if (useVerbose) {
-
- std::cout << "survival probability " << ((useRelative) ? iter->first.survival(Eavg) : iter->first.survival())
- << " number survived " << num_survive << " total number "
- << iter->second ; //<< std::endl;
- std::set<unsigned int> ms = iter->first.mutated_sites;
- for (std::set<unsigned int>::iterator v = ms.begin(); v!=ms.end(); ++v) fprintf(stdout,"\t%d",*v);
- fprintf(stdout,"\n");
-
- }
-
- }
-
- if (useVerbose) { std::cout << "checkpoint, " << "total_deaths = " << total_deaths << std::endl; }
-
- //delete zero population strains
-
- iter=pop.begin();
-
- while (iter!=pop.end()) {
-
- if (iter->second==0) pop.erase(iter++);
- else ++iter;
-
- }
-
- unsigned int ran;
- unsigned int selector;
- unsigned int pop_size=compute_population_size();
-
- if (useVerbose) print_population_size();
-
- virus_map prev_gen = pop;
-
- for (unsigned int i=0; i<total_deaths; ++i) {
-
- ran = gsl_rng_uniform_int(r,pop_size+1);
- virus_map::iterator iter = prev_gen.begin();
- selector = iter->second;
-
- while (selector<ran) {
-
- ++iter;
- selector += iter->second;
-
- }
-
- pop[iter->first]++; // segfaults HERE when pop_size == 0 (i.e. none survive)
-
- new_Eavg += iter->first.energy;
-
- }
-
- Eavg = new_Eavg/((double) N);
-
-}
-
-
-// Output population to file
-
-void Population::write_population(std::string filename) {
-
- std::ofstream fout;
- fout.open(filename.c_str());
-
- std::set<unsigned int> ms;
-
- for (virus_map::iterator iter=pop.begin(); iter!=pop.end(); ++iter) {
-
- fout << iter->second << '\t';
- fout << iter->first.energy << '\t';
-
- ms = iter->first.mutated_sites;
-
- for (std::set<unsigned int>::iterator ms_iter=ms.begin(); ms_iter!=ms.end(); ++ms_iter) {
-
- fout << *ms_iter << '\t';
-
- }
-
- fout << std::endl;
-
- }
-
- fout.close();
-
-}
-
-
-// Output population to file
-
-void Population::write_population(FILE *output, unsigned int generation) {
-
- //fprintf(output,"%d\n",generation);
-
- std::set<unsigned int> ms;
-
- for (virus_map::iterator iter=pop.begin(); iter!=pop.end(); ++iter) {
-
- //fprintf(output,"%d\t%.6e",iter->second,iter->first.energy);
- fprintf(output,"%d\t%d\t%.6e",generation,iter->second,iter->first.energy);
-
- ms = iter->first.mutated_sites;
-
- //for (std::set<unsigned int>::iterator ms_iter=ms.begin(); ms_iter!=ms.end(); ++ms_iter) fprintf(output,"\t%d",*ms_iter);
-
- // print sequence as CSV (comma-separated)
- fprintf(output,"\t");
- std::set<unsigned int>::iterator ms_iter=ms.begin();
- if (ms_iter != ms.end()) {
- fprintf(output,"%d",*ms_iter);
- ++ms_iter;
- }
- for (; ms_iter!=ms.end(); ++ms_iter)
- fprintf(output,",%d",*ms_iter);
-
- fprintf(output,"\n");
-
- }
-
- // don't flush, potentially slow with large population?
- //fflush(output);
-
-}
-
-
-// Compute the total population size
-
-unsigned int Population::compute_population_size() {
-
- unsigned int i=0;
-
- for (virus_map::iterator iter=pop.begin(); iter!=pop.end(); ++iter) i+=iter->second;
-
- return i;
-
-}
-
-
-// Print population size to standard out
-
-void Population::print_population_size() {
-
- std::cout << "The total population size is " << compute_population_size() << std::endl;
-
-}
-
-
-//Mutate every virus in the population
-
-void Population::mutate_population(const Hamiltonian &H, gsl_rng* r) {
-
- virus_map prev_generation = pop;
- unsigned int num_mutate;
-
- for(virus_map::iterator iter=prev_generation.begin(); iter!=prev_generation.end(); ++iter) {
-
- num_mutate=gsl_ran_binomial(r,p,iter->second);
-
- pop[iter->first]=pop[iter->first]-num_mutate;
-
- if (pop[iter->first]<=0) pop.erase(iter->first);
-
- for (unsigned int i=0; i<num_mutate; ++i) {
-
- Virus V(H,mu,iter->first.mutated_sites);
- V.mutate(H,r);
-
- if (pop.count(V)==0) pop[V] = 1;
- else pop[V] += 1;
-
- }
-
- }
-
-}
-
-
-// Compute the number of escaped viruses in the population
-
-unsigned int Population::compute_num_escaped(Hamiltonian &H) {
-
- unsigned int i=0;
-
- for (virus_map::iterator iter=pop.begin(); iter!=pop.end(); ++iter) { if (H.escaped(iter->first)) i+=iter->second; }
-
- return i;
-
-}
-
-
-// Check whether most of the population has escaped from immune pressure
-
-bool Population::escaped(Hamiltonian &H) {
-
- if (2*compute_num_escaped(H)>N) return true;
- else return false;
-
-}
-
-
-// Compute the number of escaped viruses in the population
-
-unsigned int Population::escape_variant(Hamiltonian &H, std::set<unsigned int> &mutant) {
-
- unsigned int i=0;
-
- for (virus_map::iterator iter=pop.begin(); iter!=pop.end(); ++iter) {
-
- if (H.escaped(iter->first) && (iter->second)>i) {
-
- i=iter->second;
- mutant=iter->first.mutated_sites;
-
- }
-
- }
-
- return i;
-
-}
-
-
-
-// Compute the number of escaped viruses in the population
-
-unsigned int Population::compute_num_escaped_all(Hamiltonian &H) {
-
- unsigned int i=0;
-
- for (virus_map::iterator iter=pop.begin(); iter!=pop.end(); ++iter) { if (H.escaped_all(iter->first)) i+=iter->second; }
-
- return i;
-
-}
-
-
-// Check whether most of the population has escaped from *all* immune pressure
-
-bool Population::escaped_all(Hamiltonian &H) {
-
- if (2*compute_num_escaped_all(H)>N) return true;
- else return false;
-
-}
-
-
-// Compute the number of escaped viruses in the population
-
-unsigned int Population::escape_variant_all(Hamiltonian &H, std::set<unsigned int> &mutant) {
-
- unsigned int i=0;
-
- for (virus_map::iterator iter=pop.begin(); iter!=pop.end(); ++iter) {
-
- if (H.escaped_all(iter->first) && (iter->second)>i) {
-
- i=iter->second;
- mutant=iter->first.mutated_sites;
-
- }
-
- }
-
- return i;
-
-}
-
-
-
-
-// Output population to file
-// two-site two-allele
-
-void Population::write_two_site_population(FILE *output, const Hamiltonian &H, unsigned int generation) {
-
- fprintf(output,"%d",generation);
-
- //static const unsigned int v00_arr[] = {};
- static const unsigned int v10_arr[] = {0};
- static const unsigned int v01_arr[] = {1};
- static const unsigned int v11_arr[] = {0,1};
-
- //std::set<unsigned int> v00_set (v00_arr, v00_arr+sizeof(v00_arr)/sizeof(v00_arr[0]));
- std::set<unsigned int> v10_set (v10_arr, v10_arr+sizeof(v10_arr)/sizeof(v10_arr[0]));
- std::set<unsigned int> v01_set (v01_arr, v01_arr+sizeof(v01_arr)/sizeof(v01_arr[0]));
- std::set<unsigned int> v11_set (v11_arr, v11_arr+sizeof(v11_arr)/sizeof(v11_arr[0]));
-
- Virus v00_vir (H, mu);//, v00_set);
- Virus v10_vir (H, mu, v10_set);
- Virus v01_vir (H, mu, v01_set);
- Virus v11_vir (H, mu, v11_set);
-
- fprintf(output,"\t%d",pop[v00_vir]);
- fprintf(output,"\t%d",pop[v10_vir]);
- fprintf(output,"\t%d",pop[v01_vir]);
- fprintf(output,"\t%d",pop[v11_vir]);
- fprintf(output,"\n");
-
- fflush(output);
-
-}
-#endif
-
-
-
std::string couplingsFile = "neutral_2site.j";
Hamiltonian H(couplingsFile);
- H.set_temp(0.01);
+ H.set_temp(r.bh, r.bJ);
+ //H.set_temp(0.01);
//double mu = 6.0e-5;
//double mu = 1.0e-1;
double sample_interval = 1.0;
simulate(reactions, species, t_end, max_steps, sample_interval, print_spec);
-
-#if 0
-
-/*
- double total_propensity = 0.0;
- for (std::vector<Reaction*>::iterator rxn = reactions.begin(),
- end = reactions.end();
- rxn != end; ++rxn) {
- (*rxn)->recalc();
- total_propensity += (*rxn)->propensity;
- }
- std::cout << total_propensity << std::endl;
-
- std::cout << ((VirusSpecies*)species.back())->pop.begin()->second << std::endl;
- std::cout << ((VirusReaction*)reactions.back())->V->pop.begin()->second << std::endl;
-
- reactions.back()->fire();
-
- std::cout << ((VirusSpecies*)species.back())->pop.begin()->second << std::endl;
- std::cout << ((VirusReaction*)reactions.back())->V->pop.begin()->second << std::endl;
-*/
-
- double t = 0.0;
- double t_end = 1.0;
- double t_remain = t_end;
- double dt = 0.0;
-
- double total_propensity = 0.0;
-
- Rxn_parray::iterator rxn, end;
-
- for (rxn = reactions.begin(), end = reactions.end(); rxn != end; ++rxn) {
- (*rxn)->recalc();
- total_propensity += (*rxn)->propensity;
- }
-
-
- while (1) {
-
- double rand = (1.0-gsl_rng_uniform_pos(rnd)) * total_propensity; // uniform on (0,R]
- double a_sum = 0.0;
-
- // time to reaction
- dt = -log(gsl_rng_uniform(rnd)) / total_propensity;
- t += dt;
- t_remain -= dt;
-
- // stop if time ran out
- if (t_remain < 0.0)
- break;
-
- // determine next reaction
- for (rxn = reactions.begin(), end = reactions.end(); rxn != end; ++rxn) {
- a_sum += (*rxn)->propensity;
- if (rand <= a_sum) {
- (*rxn)->fire();
- break;
- }
- }
-
- // TODO: handle case with total propensity == 0.0
- if (rxn == end)
- break;
-
- // recalculate rates
- total_propensity = 0.0;
- for (rxn = reactions.begin(), end = reactions.end(); rxn != end; ++rxn) {
- (*rxn)->recalc();
- total_propensity += (*rxn)->propensity;
- }
-
- // print copy numbers
- std::cout << t << '\t';
- for (std::vector<Species*>::iterator spec = species.begin(),
- end = species.end();
- spec != end; ++spec) {
- std::cout << (*spec)->count << '\t';
- }
-
- MutatedSiteSequence ms;
- FILE* output = stdout;
- for (virus_map::iterator iter = s1.pop.begin(),
- end = s1.pop.end();
- iter != end; ++iter) {
- std::cout << iter->second << '(';
- ms = iter->first.mutated_sites;
- // print sequence as CSV (comma-separated)
- MutatedSiteSequence::iterator ms_iter=ms.begin();
- if (ms_iter != ms.end()) {
- fprintf(output,"%d",*ms_iter);
- ++ms_iter;
- }
- for (; ms_iter!=ms.end(); ++ms_iter)
- fprintf(output,",%d",*ms_iter);
- std::cout << ')';
- std::cout << '\t';
- }
- std::cout << '\n';
-
- }
-
-
-#endif
-
-#if 0
- std::vector<Species*> species;
- SimpleSpecies s1; s1.count = 100; s1.name = "A";
- SimpleSpecies s2; s2.count = 50; s2.name = "B";
- species.push_back(&s1);
- species.push_back(&s2);
-
- std::vector<Reaction*> reactions;
- // A -> B
- ElementaryReaction* r1 = new ElementaryReaction();
- r1->rate_constant = 2.0;
- r1->reactants.push_back(&s1);
- r1->products.push_back(&s2);
- r1->reactant_stoich.push_back(1);
- r1->product_stoich.push_back(1);
- reactions.push_back(r1);
- // B -> A
- ElementaryReaction* r2 = new ElementaryReaction();
- r2->rate_constant = 1.0;
- r2->products.push_back(&s1);
- r2->reactants.push_back(&s2);
- r2->reactant_stoich.push_back(1);
- r2->product_stoich.push_back(1);
- reactions.push_back(r2);
-
- double t = 0.0;
- double t_end = 100.0;
- double dt = 0.0;
-
- double total_propensity = 0.0;
- for (std::vector<Reaction*>::iterator rxn = reactions.begin(),
- end = reactions.end();
- rxn != end; ++rxn) {
- (*rxn)->recalc();
- total_propensity += (*rxn)->propensity;
- }
-
-
- while (t < t_end) {
-
- // determine next reaction
-
- double rand = gsl_rng_uniform(rnd) * total_propensity; // uniform on [0,R)
- double a_sum = 0.0;
-
- dt = -log(gsl_rng_uniform(rnd)) / total_propensity;
- t += dt;
-
- for (std::vector<Reaction*>::iterator rxn = reactions.begin(),
- end = reactions.end();
- rxn != end; ++rxn) {
- a_sum += (*rxn)->propensity;
- if (rand <= a_sum) {
- (*rxn)->fire();
- break;
- }
- }
-
- // TODO: handle case with total propensity == 0.0
-
- // recalculate rates
- total_propensity = 0.0;
- for (std::vector<Reaction*>::iterator rxn = reactions.begin(),
- end = reactions.end();
- rxn != end; ++rxn) {
- (*rxn)->recalc();
- total_propensity += (*rxn)->propensity;
- }
-
-
- // print copy numbers
- std::cout << t << '\t';
- for (std::vector<Species*>::iterator spec = species.begin(),
- end = species.end();
- spec != end; ++spec) {
- std::cout << (*spec)->count << '\t';
- }
- std::cout << '\n';
-
- }
-#endif
-
-#if 0
- long NA = 100;
- long NB = 50;
-
- double kAB = 2.0;
- double kBA = 1.0;
-
- double t = 0.0;
- double t_end = 100.0;
-
- double total_propensity = 0.0;
- double rand = 0.0;
- double dt = 0.0;
- std::vector<double> prop(2);
- std::vector<double> partial_sum_prop(2);
- prop[0] = kAB*NA;
- prop[1] = kBA*NB;
- std::partial_sum(prop.begin(),prop.end(),partial_sum_prop.begin());
- //for (std::vector<double>::iterator i = partial_sum_prop.begin(); i != partial_sum_prop.end(); ++i) {
- // std::cout << (*i) << '\n';
- //}
- while (t < t_end) {
- //total_propensity = kAB*NA + kBA*NB;
- total_propensity = partial_sum_prop.back();
- rand = gsl_rng_uniform(rnd) * total_propensity; // uniform on [0,R)
- std::vector<double>::iterator i = partial_sum_prop.begin(),
- end = partial_sum_prop.end();
- for (; i != end; ++i) {
- if (*i > rand) {
- break;
- //std::cout << *i << '\n';
- }
- }
- switch(i - partial_sum_prop.begin()) {
- case 0:
- if (NA != 0) {
- NA -= 1;
- NB += 1;
- prop[0] = kAB*NA;
- }
- break;
- case 1:
- if (NB != 0) {
- NA += 1;
- NB -= 1;
- prop[1] = kBA*NB;
- }
- break;
- }
- std::partial_sum(prop.begin(),prop.end(),partial_sum_prop.begin());
- dt = -log(gsl_rng_uniform(rnd)) / total_propensity;
- t += dt;
- std::cout << t << '\t' << NA << '\t' << NB
- << '\n';
- }
-#endif
-
-
-
-#if 0
-
- r.setFiles();
- FILE *popout=fopen(r.trajectoryOutfile.c_str(),"w"); // <outfile>.dat
- FILE *supout=r.useTwoSite ? NULL : fopen(r.supplementaryOutfile.c_str(),"w"); // <outfile>.sum
- if (popout == NULL) { perror("File error"); exit(1); }
- if (r.useTwoSite && supout == NULL) { perror("File error"); exit(1); }
-
- if (r.importState) importState(r); // <infile>.st
- if (r.useEpitope) importEpitope(r); // <infile>.ep
- fflush(stdout);
-
-
- if(r.useTwoSite) {
-
- fprintf(popout,"gen\tV00\tV10\tV01\tV11\n");
-
- for (unsigned int n=0;n<r.num_runs;n++) {
-
- //fprintf(popout,"gen\tV00\tV10\tV01\tV11\n");
-
- TwoSiteHamiltonian H(r.h1,r.h2,r.J12);
- H.set_temp(r.bh, r.bJ);
- Population P(H, r.n, r.mu, r.initPop, r.initFrac);
-
- P.write_two_site_population(popout,H,0);
-
- unsigned int i;
- for (i=0; i<r.g; ++i) {
-
- //printf("%d\t%d\n",n,i);
- P.next_generation(H, rnd, r.useRelative, r.useVerbose);
-
- if ((i+1) % r.write_mod == 0) {
- P.write_two_site_population(popout,H,i+1);
- }
-
- //if (r.runUntilEscape && P.escaped(H)) break;
-
- }
-
- }
-
- }
-
-
- else {
-
- // Run (w/ targeted epitope)
-
- if (r.useEpitope) {
-
- for (unsigned int n=0; n<r.num_runs; n++) {
-
- EpitopeHamiltonian H(r.couplingsInfile);
-
- double penalty = 0.0;
- if (r.penaltyType == RunParameters_SS::PenaltyEACH)
- penalty = r.penalty;
- else if (r.penaltyType == RunParameters_SS::PenaltyTOTAL)
- penalty = r.penalty / (double)r.numEpitopes;
- for (unsigned ep=0; ep<r.numEpitopes; ++ep) {
- H.set_epitope(r.eWT[ep], r.eMut[ep], penalty);
- }
-
- H.set_temp(r.bh, r.bJ);
-
- Population P(H, r.n, r.mu, r.initPop, r.initFrac);
-
- // print epitopes
- /*
- std::cout << "-----------\n";
- for (unsigned ep=0; ep<H.penalty.size(); ++ep) {
- for (unsigned i=0; i<H.epitopeWT[ep].size(); ++i) {
- std::cout << H.epitopeWT[ep][i] << " ";
- }
- std::cout << "| ";
- for (unsigned i=0; i<H.epitopeMut[ep].size(); ++i) {
- std::cout << H.epitopeMut[ep][i] << " ";
- }
- std::cout << "\n";
- }
- fflush(stdout);
- */
-
- unsigned int i;
- for (i=0; i<r.g; ++i) {
-
- P.next_generation(H, rnd, r.useRelative, r.useVerbose);
-
- if ((i+1) % r.write_mod == 0) {
- P.write_population(popout,i+1);
- }
-
- if (r.runUntilEscape && P.escaped(H)) break;
-
- if (r.runUntilEscape_all && P.escaped_all(H)) break;
-
- }
-
- std::set<unsigned int> esc_var;
- unsigned int esc_num = 0;
- if (!r.runUntilEscape && !r.runUntilEscape_all) esc_num = P.escape_variant(H, esc_var); // default behavior for backwards compatibility
- if (r.runUntilEscape_all) esc_num = P.escape_variant_all(H, esc_var);
- if (r.runUntilEscape) esc_num = P.escape_variant(H, esc_var);
- // should behave better if both are true
- // but more likely to escape one than all, so that one is output
-
- fprintf(supout,"%d\t%d",i,esc_num);
-
- //for (std::set<unsigned int>::iterator iter=esc_var.begin(); iter!=esc_var.end(); ++iter) fprintf(supout,"\t%d",*iter);
-
- // print sequence as CSV (comma-separated)
- fprintf(supout,"\t");
- std::set<unsigned int>::iterator ms_iter=esc_var.begin();
- if (ms_iter != esc_var.end()) {
- fprintf(supout,"%d",*ms_iter);
- ++ms_iter;
- }
- for (; ms_iter!=esc_var.end(); ++ms_iter)
- fprintf(supout,",%d",*ms_iter);
-
- fprintf(supout,"\n");
-
- fflush(supout);
-
- }
-
- }
-
-
- // Run (w/out targeted epitope)
-
- else {
-
- for (unsigned int n=0;n<r.num_runs;n++) {
-
- Hamiltonian H(r.couplingsInfile);
- H.set_temp(r.bh, r.bJ);
- Population P(H, r.n, r.mu, r.initFrac);
-
- unsigned int i;
- for (i=0; i<r.g; ++i) {
-
- P.next_generation(H, rnd, r.useRelative, r.useVerbose);
- P.write_population(popout,i+1);
-
- if (r.runUntilEscape && P.escaped(H)) break;
-
- }
-
- }
-
- }
-
- }
-
-
- fclose(popout); // close file handles
- if(!r.useTwoSite) fclose(supout);
-#endif
-
gsl_rng_free(rnd); //Free up memory from random number generator
}
std::cout << '\n';
}
-#if 0
- // print copy numbers
- std::cout << t << '\t';
- for (std::vector<Species*>::iterator spec = species.begin(),
- end = species.end();
- spec != end; ++spec) {
- std::cout << (*spec)->count << '\t';
- }
-
- MutatedSiteSequence ms;
- FILE* output = stdout;
- for (virus_map::iterator iter = s1.pop.begin(),
- end = s1.pop.end();
- iter != end; ++iter) {
- std::cout << iter->second << '(';
- ms = iter->first.mutated_sites;
- // print sequence as CSV (comma-separated)
- MutatedSiteSequence::iterator ms_iter=ms.begin();
- if (ms_iter != ms.end()) {
- fprintf(output,"%d",*ms_iter);
- ++ms_iter;
- }
- for (; ms_iter!=ms.end(); ++ms_iter)
- fprintf(output,",%d",*ms_iter);
- std::cout << ')';
- std::cout << '\t';
- }
- std::cout << '\n';
-#endif
}
" -i (string) input couplings file\n"
" -o (string) output file stem\n"
" -s (string) input state file, containing initial population fraction\n"
-" -e (string) input file containing targeted epitope\n"
+" -ep (string) input file containing targeted epitope\n"
" -en (int) corresponding number of that CTL (ordered!)\n"
" -n (int/double) population size\n"
" -g (int/double) number of generations\n"
" -mu (double) mutation rate\n"
" -b (double) \"inverse temperature\" multiplier\n"
-" -r flag to use relative energy condition for survival probability\n"
" -esc flag to run until escape observed (or max number of generations reached)\n"
" -esc_all flag to run until escape from *all* immune pressure (or max num gen)\n"
" -v flag for verbose output\n"
" -numruns (int) number of trajectories to simulate\n"
-"\n"
-" -2site flag for two-site two-allele system\n"
-" -h1 (double) value of field at site 1\n"
-" -h2 (double) value of field at site 2\n"
-" -J12 (double) value of coupling between sites 1 and 2\n"
-" (note the sign convention on these is P ~ exp(-H),\n"
-" H = sum_i h_i s_i + sum_{i<j} J_{ij} s_i s_j\n"
-" although internal representation is opposite)\n"
-" -write_mod (int) write out state every 'write_mod' generations\n"
-" -penaltyEach (double) bonus (per epitope) for mutation in epitope\n"
-" -penaltyTotal (double) maximum total bonus if all epitopes contain mutations;\n"
-" bonus per epitope = this / number of epitopes\n"
; std::cout.flush();
}
* -i (string) input couplings file
* -o (string) output file stem
* -s (string) input state file, containing initial population fraction and targeted epitope
- * -e (string) input file containing targeted epitope
+ * -ep (string) input file containing targeted epitope
* -n (int/double) population size
* -g (int/double) number of generations
* -mu (double) mutation rate
* -b (double) "inverse temperature" multiplier
- * -r flag to use relative energy condition for survival probability
* -esc flag to run until escape observed (or max number of generations reached)
* -v flag for verbose output
* -numruns (int) number of trajectories to simulate
else if (strcmp(argv[i],"-i")==0) { if (++i==argc) break; else r.infile=argv[i]; }
else if (strcmp(argv[i],"-o")==0) { if (++i==argc) break; else r.outfile=argv[i]; }
else if (strcmp(argv[i],"-s")==0) { if (++i==argc) break; else { r.statefile=argv[i]; r.importState=true; } }
- else if (strcmp(argv[i],"-e")==0) { if (++i==argc) break; else { r.epitopefiles.push_back(argv[i]); r.useEpitope=true; ++r.num_CTL_clones; } }
+ else if (strcmp(argv[i],"-ep")==0) { if (++i==argc) break; else { r.epitopefiles.push_back(argv[i]); r.useEpitope=true; ++r.num_CTL_clones; } }
else if (strcmp(argv[i],"-en")==0) { if (++i==argc) break; else { r.init_CTL_num.push_back((long)strtodouble(argv[i])); r.useEpitope=true; } }
else if (strcmp(argv[i],"-n")==0) { if (++i==argc) break; else r.n=(unsigned int)strtodouble(argv[i]); }
else if (strcmp(argv[i],"-penaltyEach")==0) { if (++i==argc) break; else { r.penalty=strtodouble(argv[i]); r.penaltyType=RunParameters_SS::PenaltyEACH; } }
else if (strcmp(argv[i],"-penaltyTotal")==0) { if (++i==argc) break; else { r.penalty=strtodouble(argv[i]); r.penaltyType=RunParameters_SS::PenaltyTOTAL; } }
- else if (strcmp(argv[i],"-r")==0) { r.useRelative=true; }
else if (strcmp(argv[i],"-esc")==0) { r.runUntilEscape=true; }
else if (strcmp(argv[i],"-esc_all")==0) { r.runUntilEscape_all=true; }
else if (strcmp(argv[i],"-v")==0) { r.useVerbose=true; }
else if (strcmp(argv[i],"-numruns")==0) { if (++i==argc) break; else r.num_runs=(unsigned int)strtodouble(argv[i]); }
- else if (strcmp(argv[i],"-2site")==0) { r.useTwoSite=true; }
- else if (strcmp(argv[i],"-h1")==0) { if (++i==argc) break; else r.h1=-strtodouble(argv[i]); }
- else if (strcmp(argv[i],"-h2")==0) { if (++i==argc) break; else r.h2=-strtodouble(argv[i]); }
- else if (strcmp(argv[i],"-J12")==0) { if (++i==argc) break; else r.J12=-strtodouble(argv[i]); }
-
- else if (strcmp(argv[i],"-p00")==0) { if (++i==argc) break; else add_to_two_site_pop(r,false,false,strtodouble(argv[i])); }
- else if (strcmp(argv[i],"-p10")==0) { if (++i==argc) break; else add_to_two_site_pop(r,true ,false,strtodouble(argv[i])); }
- else if (strcmp(argv[i],"-p01")==0) { if (++i==argc) break; else add_to_two_site_pop(r,false,true ,strtodouble(argv[i])); }
- else if (strcmp(argv[i],"-p11")==0) { if (++i==argc) break; else add_to_two_site_pop(r,true ,true ,strtodouble(argv[i])); }
-
else if (strcmp(argv[i],"-seed")==0) { if (++i==argc) break; else seed=(unsigned)strtodouble(argv[i]); }
else if (strcmp(argv[i],"-write_mod")==0) { if (++i==argc) break; else r.write_mod=(unsigned int)strtodouble(argv[i]); }
projects / VirEvoDyn.git / commitdiff