//Constructor for loading couplings from John Barton's Ising Inversion code
-//Code for this function derived from John Barton's code
+//Code for this function derived from John Barton's code
Hamiltonian::Hamiltonian(std::string &FILENAME) {
-
+
this->bh=1.0;
this->bJ=1.0;
- FILE* input;
-
input = fopen(FILENAME.c_str(),"r"); // <paramfile>.j
+ FILE* input;
+ input = fopen(FILENAME.c_str(),"r"); // <paramfile>.j
if (input == NULL) { perror((std::string("ERROR: ") + FILENAME).c_str()); exit(1); }
-
- typedef std::vector<int> Key;
- std::map<Key,double> JIndex;
+
+ typedef std::vector<int> Key;
+ std::map<Key,double> JIndex;
char o;
- int site=0;
+ int site=0;
+
+ fscanf(input,"%d",&size);
- fscanf(input,"%d",&size);
-
- while (fscanf(input,"%d",&site)==1) {
+ while (fscanf(input,"%d",&site)==1) {
- Key inputKey(2,site);
- while (fscanf(input,"%c",&o)==1) {
-
- if (o==';') break;
-
- else {
-
- int neighborSite=0;
- double neighborJ=0;
- fscanf(input,"%d %c %lf",&neighborSite,&o,&neighborJ);
-
- if (neighborSite==site) {
+ Key inputKey(2,site);
+ while (fscanf(input,"%c",&o)==1) {
+
+ if (o==';') break;
+
+ else {
+
+ int neighborSite=0;
+ double neighborJ=0;
+ fscanf(input,"%d %c %lf",&neighborSite,&o,&neighborJ);
+
+ if (neighborSite==site) {
- inputKey[1]=site;
+ inputKey[1]=site;
JIndex[inputKey]=neighborJ;
-
- }
- else {
+
+ }
+ else {
inputKey[1]=neighborSite;
JIndex[inputKey]=neighborJ;
}
//std::cout << site << " " << neighborSite << " " << neighborJ << "\n";
- }
+ }
- }
+ }
+
+ }
- }
-
- fclose(input);
- Key inputKey(2,0);
- J.resize(size);
+ fclose(input);
+ Key inputKey(2,0);
+ J.resize(size);
- for(int i=0;i<size;i++) J[i].resize(size,0.0);
+ for(int i=0;i<size;i++) J[i].resize(size,0.0);
Graph.resize(size);
- for (int i=0;i<size;i++) {
-
+ for (int i=0;i<size;i++) {
+
J[i].resize(size,0.0);
- for(int j=i;j<size;j++){
+ for(int j=i;j<size;j++){
- inputKey[0]=i;
- inputKey[1]=j;
+ inputKey[0]=i;
+ inputKey[1]=j;
- if (JIndex[inputKey]!=0.0) {
+ if (JIndex[inputKey]!=0.0) {
- J[i][j]+=JIndex[inputKey];
- J[j][i]=J[i][j];
+ J[i][j]+=JIndex[inputKey];
+ J[j][i]=J[i][j];
- if (J[i][j]!=0.0) {
-
- Graph[i].push_back(j); //Populate the adjacency list
- if(j!=i) Graph[j].push_back(i);
+ if (J[i][j]!=0.0) {
+
+ Graph[i].push_back(j); //Populate the adjacency list
+ if(j!=i) Graph[j].push_back(i);
- }
+ }
- }
+ }
- }
+ }
- }
+ }
//for(int i=0; i<size; i++)
//for(int j=0; j<size; j++)
//Constructor for loading couplings from John Barton's Ising Inversion code
-//Code for this function derived from John Barton's code
+//Code for this function derived from John Barton's code
EpitopeHamiltonian::EpitopeHamiltonian(std::string &FILENAME) {
-
+
this->bh=1.0;
this->bJ=1.0;
- FILE* input;
-
input = fopen(FILENAME.c_str(),"r"); // <paramfile>.j
+ FILE* input;
+ input = fopen(FILENAME.c_str(),"r"); // <paramfile>.j
if (input == NULL) { perror((std::string("ERROR: ") + FILENAME).c_str()); exit(1); }
-
- typedef std::vector<int> Key;
- std::map<Key,double> JIndex;
+
+ typedef std::vector<int> Key;
+ std::map<Key,double> JIndex;
char o;
- int site=0;
+ int site=0;
- fscanf(input,"%d",&size);
-
- while (fscanf(input,"%d",&site)==1) {
+ fscanf(input,"%d",&size);
+
+ while (fscanf(input,"%d",&site)==1) {
Key inputKey(2,site);
- while (fscanf(input,"%c",&o)==1) {
-
- if (o==';') break;
-
- else {
-
- int neighborSite=0;
- double neighborJ=0;
- fscanf(input,"%d %c %lf",&neighborSite,&o,&neighborJ);
-
- if (neighborSite==site) {
+ while (fscanf(input,"%c",&o)==1) {
+
+ if (o==';') break;
+
+ else {
+
+ int neighborSite=0;
+ double neighborJ=0;
+ fscanf(input,"%d %c %lf",&neighborSite,&o,&neighborJ);
+
+ if (neighborSite==site) {
- inputKey[1]=site;
+ inputKey[1]=site;
JIndex[inputKey]=neighborJ;
-
- }
- else {
+
+ }
+ else {
inputKey[1]=neighborSite;
JIndex[inputKey]=neighborJ;
}
- }
+ }
- }
+ }
- }
-
- fclose(input);
- Key inputKey(2,0);
- J.resize(size);
+ }
+
+ fclose(input);
+ Key inputKey(2,0);
+ J.resize(size);
- for(int i=0;i<size;i++) J[i].resize(size,0.0);
+ for(int i=0;i<size;i++) J[i].resize(size,0.0);
Graph.resize(size);
- for (int i=0;i<size;i++) {
-
+ for (int i=0;i<size;i++) {
+
J[i].resize(size,0.0);
- for(int j=i;j<size;j++){
+ for(int j=i;j<size;j++){
- inputKey[0]=i;
- inputKey[1]=j;
+ inputKey[0]=i;
+ inputKey[1]=j;
- if (JIndex[inputKey]!=0.0) {
+ if (JIndex[inputKey]!=0.0) {
- J[i][j]+=JIndex[inputKey];
- J[j][i]=J[i][j];
+ J[i][j]+=JIndex[inputKey];
+ J[j][i]=J[i][j];
- if (J[i][j]!=0.0) {
-
- Graph[i].push_back(j); //Populate the adjacency list
- if(j!=i) Graph[j].push_back(i);
+ if (J[i][j]!=0.0) {
+
+ Graph[i].push_back(j); //Populate the adjacency list
+ if(j!=i) Graph[j].push_back(i);
- }
+ }
- }
+ }
- }
+ }
- }
+ }
}
double efield = 0.0;
double ecoupling = 0.0;
- std::set<unsigned int>::iterator iter1;
- std::set<unsigned int>::iterator iter2;
+ std::set<unsigned int>::iterator iter1;
+ std::set<unsigned int>::iterator iter2;
- for (iter1=mutated_sites.begin();iter1!=mutated_sites.end();++iter1) {
+ for (iter1=mutated_sites.begin();iter1!=mutated_sites.end();++iter1) {
- efield -= J[*iter1][*iter1];
-
+ efield -= J[*iter1][*iter1];
+
iter2 = iter1;
++iter2;
-
+
for (;iter2!=mutated_sites.end();iter2++) ecoupling -= J[*iter1][*iter2];
-
+
}
return ((efield * bh) + (ecoupling * bJ));
double efield = 0.0;
double ecoupling = 0.0;
- std::set<unsigned int>::iterator iter1;
- std::set<unsigned int>::iterator iter2;
+ std::set<unsigned int>::iterator iter1;
+ std::set<unsigned int>::iterator iter2;
- for (iter1=mutated_sites.begin();iter1!=mutated_sites.end();++iter1) {
+ for (iter1=mutated_sites.begin();iter1!=mutated_sites.end();++iter1) {
- efield -= J[*iter1][*iter1];
-
+ efield -= J[*iter1][*iter1];
+
iter2 = iter1;
++iter2;
-
+
for (;iter2!=mutated_sites.end();iter2++) ecoupling -= J[*iter1][*iter2];
-
+
}
double energy = (efield * bh) + (ecoupling * bJ);
}
TwoSiteHamiltonian::TwoSiteHamiltonian(std::string &FILENAME) {
-
+
this->bh=1.0;
this->bJ=1.0;
- FILE* input;
- input = fopen(FILENAME.c_str(),"r");
+ FILE* input;
+ input = fopen(FILENAME.c_str(),"r");
if (input == NULL) { perror((std::string("ERROR: ") + FILENAME).c_str()); exit(1); }
-
- typedef std::vector<int> Key;
- std::map<Key,double> JIndex;
+
+ typedef std::vector<int> Key;
+ std::map<Key,double> JIndex;
char o;
- int site=0;
+ int site=0;
- fscanf(input,"%d",&size);
-
- while (fscanf(input,"%d",&site)==1) {
+ fscanf(input,"%d",&size);
+
+ while (fscanf(input,"%d",&site)==1) {
Key inputKey(2,site);
- while (fscanf(input,"%c",&o)==1) {
-
- if (o==';') break;
-
- else {
-
- int neighborSite=0;
- double neighborJ=0;
- fscanf(input,"%d %c %lf",&neighborSite,&o,&neighborJ);
-
- if (neighborSite==site) {
+ while (fscanf(input,"%c",&o)==1) {
+
+ if (o==';') break;
+
+ else {
+
+ int neighborSite=0;
+ double neighborJ=0;
+ fscanf(input,"%d %c %lf",&neighborSite,&o,&neighborJ);
+
+ if (neighborSite==site) {
- inputKey[1]=site;
+ inputKey[1]=site;
JIndex[inputKey]=neighborJ;
-
- }
- else {
+
+ }
+ else {
inputKey[1]=neighborSite;
JIndex[inputKey]=neighborJ;
}
//std::cout << site << " " << neighborSite << " " << neighborJ << "\n";
- }
+ }
- }
+ }
+
+ }
- }
-
- fclose(input);
- Key inputKey(2,0);
- J.resize(size);
+ fclose(input);
+ Key inputKey(2,0);
+ J.resize(size);
- for(int i=0;i<size;i++) J[i].resize(size,0.0);
+ for(int i=0;i<size;i++) J[i].resize(size,0.0);
Graph.resize(size);
- for (int i=0;i<size;i++) {
-
+ for (int i=0;i<size;i++) {
+
J[i].resize(size,0.0);
- for(int j=i;j<size;j++){
+ for(int j=i;j<size;j++){
- inputKey[0]=i;
- inputKey[1]=j;
+ inputKey[0]=i;
+ inputKey[1]=j;
- if (JIndex[inputKey]!=0.0) {
+ if (JIndex[inputKey]!=0.0) {
- J[i][j]+=JIndex[inputKey];
- J[j][i]=J[i][j];
+ J[i][j]+=JIndex[inputKey];
+ J[j][i]=J[i][j];
- if (J[i][j]!=0.0) {
-
- Graph[i].push_back(j); //Populate the adjacency list
- if(j!=i) Graph[j].push_back(i);
+ if (J[i][j]!=0.0) {
+
+ Graph[i].push_back(j); //Populate the adjacency list
+ if(j!=i) Graph[j].push_back(i);
- }
+ }
- }
+ }
- }
+ }
- }
+ }
//for(int i=0; i<size; i++)
//for(int j=0; j<size; j++)
//Constructor that assembles a population of N viruses of wild type
Population::Population(const Hamiltonian &H, unsigned int N, double mu) {
-
- Virus V(H, mu);
- pop[V] = N;
- this->N = N;
+
+ Virus V(H, mu);
+ pop[V] = N;
+ this->N = N;
this->mu = mu;
- p = 1.0-pow((1.0-V.mu),H.size);
+ 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]);
Eavg += V.energy * pop[V];
}
-
+
this->N = N;
this->mu = mu;
- p = 1.0-pow((1.0-mu),H.size);
+ p = 1.0-pow((1.0-mu),H.size);
Eavg /= (double) N;
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;
+ 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){
+ 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);
+ 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;
+ iter->second = num_survive;
new_Eavg += iter->first.energy * num_survive;
-
+
// Report survival
if (useVerbose) {
fprintf(stdout,"\n");
}
-
- }
+
+ }
- if (useVerbose) { std::cout << "checkpoint, " << "total_deaths = " << total_deaths << std::endl; }
+ if (useVerbose) { std::cout << "checkpoint, " << "total_deaths = " << total_deaths << std::endl; }
- //delete zero population strains
+ //delete zero population strains
- iter=pop.begin();
+ iter=pop.begin();
- while (iter!=pop.end()) {
+ while (iter!=pop.end()) {
- if (iter->second==0) pop.erase(iter++);
- else ++iter;
+ if (iter->second==0) pop.erase(iter++);
+ else ++iter;
- }
+ }
- unsigned int ran;
- unsigned int selector;
- unsigned int pop_size=compute_population_size();
+ unsigned int ran;
+ unsigned int selector;
+ unsigned int pop_size=compute_population_size();
- if (useVerbose) print_population_size();
+ if (useVerbose) print_population_size();
- virus_map prev_gen = pop;
+ virus_map prev_gen = pop;
- for (unsigned int i=0; i<total_deaths; ++i) {
+ 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;
+ ran = gsl_rng_uniform_int(r,pop_size+1);
+ virus_map::iterator iter = prev_gen.begin();
+ selector = iter->second;
- while (selector<ran) {
-
+ while (selector<ran) {
+
++iter;
- selector += iter->second;
-
+ selector += iter->second;
+
}
- pop[iter->first]++; // segfaults HERE when pop_size == 0 (i.e. none survive)
+ pop[iter->first]++; // segfaults HERE when pop_size == 0 (i.e. none survive)
new_Eavg += iter->first.energy;
- }
+ }
Eavg = new_Eavg/((double) N);
void Population::write_population(std::string filename) {
- std::ofstream fout;
- fout.open(filename.c_str());
+ std::ofstream fout;
+ fout.open(filename.c_str());
- std::set<unsigned int> ms;
+ std::set<unsigned int> ms;
- for (virus_map::iterator iter=pop.begin(); iter!=pop.end(); ++iter) {
-
+ for (virus_map::iterator iter=pop.begin(); iter!=pop.end(); ++iter) {
+
fout << iter->second << '\t';
- fout << iter->first.energy << '\t';
+ fout << iter->first.energy << '\t';
- ms = iter->first.mutated_sites;
+ ms = iter->first.mutated_sites;
- for (std::set<unsigned int>::iterator ms_iter=ms.begin(); ms_iter!=ms.end(); ++ms_iter) {
+ for (std::set<unsigned int>::iterator ms_iter=ms.begin(); ms_iter!=ms.end(); ++ms_iter) {
- fout << *ms_iter << '\t';
+ fout << *ms_iter << '\t';
- }
+ }
- fout << std::endl;
+ fout << std::endl;
- }
+ }
- fout.close();
+ fout.close();
}
void Population::write_population(FILE *output, unsigned int generation) {
- //fprintf(output,"%d\n",generation);
+ //fprintf(output,"%d\n",generation);
std::set<unsigned int> ms;
- for (virus_map::iterator iter=pop.begin(); iter!=pop.end(); ++iter) {
-
+ 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);
+ //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");
fprintf(output,"%d",*ms_iter);
++ms_iter;
}
- for (; ms_iter!=ms.end(); ++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);
+ //fflush(output);
}
unsigned int Population::compute_population_size() {
- unsigned int i=0;
+ unsigned int i=0;
- for (virus_map::iterator iter=pop.begin(); iter!=pop.end(); ++iter) i+=iter->second;
-
- return i;
+ for (virus_map::iterator iter=pop.begin(); iter!=pop.end(); ++iter) i+=iter->second;
+
+ return i;
}
void Population::print_population_size() {
- std::cout << "The total population size is " << compute_population_size() << std::endl;
+ 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;
+
+ virus_map prev_generation = pop;
unsigned int num_mutate;
- for(virus_map::iterator iter=prev_generation.begin(); iter!=prev_generation.end(); ++iter) {
-
+ 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;
+ pop[iter->first]=pop[iter->first]-num_mutate;
- if (pop[iter->first]<=0) pop.erase(iter->first);
+ if (pop[iter->first]<=0) pop.erase(iter->first);
- for (unsigned int i=0; i<num_mutate; ++i) {
+ for (unsigned int i=0; i<num_mutate; ++i) {
- Virus V(H,mu,iter->first.mutated_sites);
- V.mutate(H,r);
-
+ 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;
+ 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;
+ for (virus_map::iterator iter=pop.begin(); iter!=pop.end(); ++iter) { if (H.escaped(iter->first)) i+=iter->second; }
+
+ return i;
}
bool Population::escaped(Hamiltonian &H) {
- if (2*compute_num_escaped(H)>N) return true;
+ if (2*compute_num_escaped(H)>N) return true;
else return false;
}
unsigned int Population::escape_variant(Hamiltonian &H, std::set<unsigned int> &mutant) {
- unsigned int i=0;
+ unsigned int i=0;
- for (virus_map::iterator iter=pop.begin(); iter!=pop.end(); ++iter) {
+ for (virus_map::iterator iter=pop.begin(); iter!=pop.end(); ++iter) {
if (H.escaped(iter->first) && (iter->second)>i) {
unsigned int Population::compute_num_escaped_all(Hamiltonian &H) {
- unsigned int i=0;
+ 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;
+ for (virus_map::iterator iter=pop.begin(); iter!=pop.end(); ++iter) { if (H.escaped_all(iter->first)) i+=iter->second; }
+
+ return i;
}
bool Population::escaped_all(Hamiltonian &H) {
- if (2*compute_num_escaped_all(H)>N) return true;
+ if (2*compute_num_escaped_all(H)>N) return true;
else return false;
}
unsigned int Population::escape_variant_all(Hamiltonian &H, std::set<unsigned int> &mutant) {
- unsigned int i=0;
+ unsigned int i=0;
- for (virus_map::iterator iter=pop.begin(); iter!=pop.end(); ++iter) {
+ for (virus_map::iterator iter=pop.begin(); iter!=pop.end(); ++iter) {
if (H.escaped_all(iter->first) && (iter->second)>i) {
void Population::write_two_site_population(FILE *output, const Hamiltonian &H, unsigned int generation) {
- fprintf(output,"%d",generation);
+ fprintf(output,"%d",generation);
//static const unsigned int v00_arr[] = {};
static const unsigned int v10_arr[] = {0};
fprintf(output,"\t%d",pop[v11_vir]);
fprintf(output,"\n");
- fflush(output);
+ fflush(output);
}
class Population{
-
+
public:
- unsigned int N; // Population size
+ unsigned int N; // Population size
double p; // Probability that a sequence has one or more mutations
double Eavg; // Average energy of the sequence population
double mu; // Mutation rate
virus_map pop;
-
+
Population(const Hamiltonian &H, unsigned int N, double mu);
Population(const Hamiltonian &H, unsigned int N, double mu, const std::vector<std::set<unsigned int> > &initPop, const std::vector<double> &initFrac);
unsigned int compute_num_escaped_all(Hamiltonian &H);
bool escaped_all(Hamiltonian &H);
unsigned int escape_variant_all(Hamiltonian &H, std::set<unsigned int> &mutant);
-
+
};
#endif
//Constuct a virus object of wildtype
Virus::Virus(const Hamiltonian &H, double mu) {
-
- this->mutated_sites.clear();
- this->mu=mu;
- L=H.size;
- energy=0;
+
+ this->mutated_sites.clear();
+ this->mu=mu;
+ L=H.size;
+ energy=0;
}
//Construct a virus and compute its energy.
Virus::Virus(const Hamiltonian &H, double mu, const std::set<unsigned int> &mutated_sites) {
-
- this->mutated_sites=mutated_sites;
- this->mu=mu;
- L=H.size;
- update_energy(H);
+
+ this->mutated_sites=mutated_sites;
+ this->mu=mu;
+ L=H.size;
+ update_energy(H);
-}
+}
//Print key numerical parameters of the object to the terminal for diagnostics.
void Virus::print_parameters() {
-
- std::cout << "mu = " << mu << std::endl;
- std::cout << "energy = " << energy << std::endl;
+
+ std::cout << "mu = " << mu << std::endl;
+ std::cout << "energy = " << energy << std::endl;
}
//to an instance of a gsl_rng random number generator
void Virus::mutate(const Hamiltonian &H, gsl_rng* r) {
-
- unsigned int n=gsl_ran_binomial(r,mu,L);
- while (n<1) n=gsl_ran_binomial(r,mu,L);
+ unsigned int n=gsl_ran_binomial(r,mu,L);
+
+ while (n<1) n=gsl_ran_binomial(r,mu,L);
for (unsigned i=0;i<n;i++) {
double Virus::survival() const {
- double z = exp(-energy);
- return z/(1+z);
+ double z = exp(-energy);
+ return z/(1+z);
}
double Virus::survival(double Eavg) const {
- double z = exp(Eavg-energy);
- return z/(1+z);
+ double z = exp(Eavg-energy);
+ return z/(1+z);
}
bool operator<(const Virus& lhs, const Virus& rhs) {
-
+
return lhs.mutated_sites < rhs.mutated_sites;
}
double energy;
std::set<unsigned int> mutated_sites;
-
+
Virus(const Hamiltonian &H, double mu);
Virus(const Hamiltonian &H, double mu, const std::set<unsigned int> &mutated_sites);
double survival(double Eavg) const;
friend class Population;
-
-//private:
+
+private:
double mu;
unsigned int L;
-
+
void update_energy(const Hamiltonian &H);
};
// Initialize RNG and set initial state, if importing from file
- static gsl_rng *rnd = gsl_rng_alloc(gsl_rng_taus2); //pointer to random number generator state
+ static gsl_rng *rnd = gsl_rng_alloc(gsl_rng_taus2); //pointer to random number generator state
//srand((unsigned)time(0));
//gsl_rng_set(rnd,rand());
char o;
double frac;
- unsigned int site;
+ unsigned int site;
// Import initial state of the population
- while (fscanf(input,"%le",&frac)==1) {
+ while (fscanf(input,"%le",&frac)==1) {
std::cout << frac; // XXX
std::set<unsigned int> tempSet;
- while (fscanf(input,"%c",&o)==1) {
+ while (fscanf(input,"%c",&o)==1) {
std::cout << o; // XXX
-
- if (o=='\n' || o==';') break;
-
- else {
-
- fscanf(input,"%u",&site);
+
+ if (o=='\n' || o==';') break;
+
+ else {
+
+ fscanf(input,"%u",&site);
std::cout << site; // XXX
tempSet.insert(site);
}
-
- }
+
+ }
r.initPop.push_back(tempSet);
if (o==';') break;
- }
+ }
std::cout << "\n"; // XXX
return 0;
-}
-
+}
+
projects / VirEvoDyn.git / commitdiff