From ddf19f780db8b19a6cad513579c4306f5b1a0166 Mon Sep 17 00:00:00 2001 From: Dariusz Murakowski Date: Fri, 12 Jun 2015 19:20:19 -0400 Subject: [PATCH] :retab (replace \t with normal whitespace). --- hamiltonian.cpp | 328 ++++++++++++++++++++++++++++---------------------------- population.cpp | 196 ++++++++++++++++----------------- population.h | 8 +- virus.cpp | 44 ++++---- virus.h | 8 +- wf.cpp | 30 +++--- 6 files changed, 307 insertions(+), 307 deletions(-) diff --git a/hamiltonian.cpp b/hamiltonian.cpp index eff464c..70951c8 100644 --- a/hamiltonian.cpp +++ b/hamiltonian.cpp @@ -10,44 +10,44 @@ //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"); // .j + FILE* input; + input = fopen(FILENAME.c_str(),"r"); // .j if (input == NULL) { perror((std::string("ERROR: ") + FILENAME).c_str()); exit(1); } - - typedef std::vector Key; - std::map JIndex; + + typedef std::vector Key; + std::map 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; @@ -55,45 +55,45 @@ Hamiltonian::Hamiltonian(std::string &FILENAME) { } //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;ibh=1.0; this->bJ=1.0; - FILE* input; - input = fopen(FILENAME.c_str(),"r"); // .j + FILE* input; + input = fopen(FILENAME.c_str(),"r"); // .j if (input == NULL) { perror((std::string("ERROR: ") + FILENAME).c_str()); exit(1); } - - typedef std::vector Key; - std::map JIndex; + + typedef std::vector Key; + std::map 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 &mutated_sites) cons double efield = 0.0; double ecoupling = 0.0; - std::set::iterator iter1; - std::set::iterator iter2; + std::set::iterator iter1; + std::set::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)); @@ -318,18 +318,18 @@ double EpitopeHamiltonian::get_energy(const std::set &mutated_site double efield = 0.0; double ecoupling = 0.0; - std::set::iterator iter1; - std::set::iterator iter2; + std::set::iterator iter1; + std::set::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); @@ -412,41 +412,41 @@ TwoSiteHamiltonian::TwoSiteHamiltonian(double h1, double h2, double J12) { } 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 Key; - std::map JIndex; + + typedef std::vector Key; + std::map 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; @@ -454,45 +454,45 @@ TwoSiteHamiltonian::TwoSiteHamiltonian(std::string &FILENAME) { } //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;iN = 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 > &initPop, const std::vector &initFrac) { - + for (unsigned i=0;iN = 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; @@ -50,32 +50,32 @@ Population::Population(const Hamiltonian &H, unsigned int N, double mu, const st 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) { @@ -88,48 +88,48 @@ void Population::next_generation(const Hamiltonian &H, gsl_rng* r, bool useRelat 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; isecond; + ran = gsl_rng_uniform_int(r,pop_size+1); + virus_map::iterator iter = prev_gen.begin(); + selector = iter->second; - while (selectorsecond; - + 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); @@ -140,29 +140,29 @@ void Population::next_generation(const Hamiltonian &H, gsl_rng* r, bool useRelat 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 ms; + std::set 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::iterator ms_iter=ms.begin(); ms_iter!=ms.end(); ++ms_iter) { + for (std::set::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(); } @@ -171,18 +171,18 @@ void Population::write_population(std::string filename) { void Population::write_population(FILE *output, unsigned int generation) { - //fprintf(output,"%d\n",generation); + //fprintf(output,"%d\n",generation); std::set 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::iterator ms_iter=ms.begin(); ms_iter!=ms.end(); ++ms_iter) fprintf(output,"\t%d",*ms_iter); + //for (std::set::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"); @@ -191,15 +191,15 @@ void Population::write_population(FILE *output, unsigned int generation) { 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); } @@ -208,11 +208,11 @@ void Population::write_population(FILE *output, unsigned int generation) { 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; } @@ -221,7 +221,7 @@ unsigned int Population::compute_population_size() { 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; } @@ -229,42 +229,42 @@ void Population::print_population_size() { //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; ifirst.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; } @@ -273,7 +273,7 @@ unsigned int Population::compute_num_escaped(Hamiltonian &H) { 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; } @@ -283,9 +283,9 @@ bool Population::escaped(Hamiltonian &H) { unsigned int Population::escape_variant(Hamiltonian &H, std::set &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) { @@ -306,11 +306,11 @@ unsigned int Population::escape_variant(Hamiltonian &H, std::set & 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; } @@ -319,7 +319,7 @@ unsigned int Population::compute_num_escaped_all(Hamiltonian &H) { 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; } @@ -329,9 +329,9 @@ bool Population::escaped_all(Hamiltonian &H) { unsigned int Population::escape_variant_all(Hamiltonian &H, std::set &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) { @@ -354,7 +354,7 @@ unsigned int Population::escape_variant_all(Hamiltonian &H, std::set virus_map; 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 > &initPop, const std::vector &initFrac); @@ -40,7 +40,7 @@ public: unsigned int compute_num_escaped_all(Hamiltonian &H); bool escaped_all(Hamiltonian &H); unsigned int escape_variant_all(Hamiltonian &H, std::set &mutant); - + }; #endif diff --git a/virus.cpp b/virus.cpp index f523ca4..4dba473 100644 --- a/virus.cpp +++ b/virus.cpp @@ -11,11 +11,11 @@ //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; } @@ -23,21 +23,21 @@ Virus::Virus(const Hamiltonian &H, double mu) { //Construct a virus and compute its energy. Virus::Virus(const Hamiltonian &H, double mu, const std::set &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; } @@ -46,10 +46,10 @@ void Virus::print_parameters() { //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 mutated_sites; - + Virus(const Hamiltonian &H, double mu); Virus(const Hamiltonian &H, double mu, const std::set &mutated_sites); @@ -27,12 +27,12 @@ public: double survival(double Eavg) const; friend class Population; - -//private: + +private: double mu; unsigned int L; - + void update_energy(const Hamiltonian &H); }; diff --git a/wf.cpp b/wf.cpp index 854be90..8867c13 100644 --- a/wf.cpp +++ b/wf.cpp @@ -40,7 +40,7 @@ void run(RunParameters &r, unsigned seed) { // 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()); @@ -218,11 +218,11 @@ void importState(RunParameters &r) { 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 @@ -230,27 +230,27 @@ void importState(RunParameters &r) { std::set 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 @@ -437,5 +437,5 @@ int main(int argc, char *argv[]) { return 0; -} - +} + -- 2.7.4