From: Dariusz Murakowski <murakdar@mit.edu>
Date: Fri, 17 Apr 2015 18:25:35 +0000 (-0400)
Subject: Clean out WF-only stuff.
X-Git-Url: http://src.murakowski.org/?a=commitdiff_plain;h=7dbaf6b791a7b8dcb756014cdacc34472c236af5;p=VirEvoDyn.git

Clean out WF-only stuff.
---

diff --git a/pop_ss.cpp b/pop_ss.cpp
index e0a197b..94fa159 100644
--- a/pop_ss.cpp
+++ b/pop_ss.cpp
@@ -44,379 +44,3 @@ double CTLSpecies::recognized(const MutatedSiteSequence &mutated_sites) const {
     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
-
-
-
diff --git a/pop_ss.h b/pop_ss.h
index 12e21ab..8c264b5 100644
--- a/pop_ss.h
+++ b/pop_ss.h
@@ -11,42 +11,9 @@
 
 typedef std::map<Virus, unsigned int> virus_map;
 
-
-#if 0
-class Population{
-	
-public:
-    
-    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);
-        
-    void next_generation(const Hamiltonian &H, gsl_rng* r, bool useRelative, bool useVerbose);
-    void write_population(std::string filepath);
-    void write_population(FILE *output, unsigned int generation);
-    void write_two_site_population(FILE *output, const Hamiltonian &H, unsigned int generation);
-    unsigned int compute_population_size();
-    void print_population_size();
-    void mutate_population(const Hamiltonian &H, gsl_rng* r);
-    
-    unsigned int compute_num_escaped(Hamiltonian &H);
-    bool escaped(Hamiltonian &H);
-    unsigned int escape_variant(Hamiltonian &H, std::set<unsigned int> &mutant);
-
-    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
-
 typedef std::set<unsigned int> MutatedSiteSequence;
 
+
 class EpitopeRecognizer {
 public:
     std::vector<unsigned int> epitopeWT;  // sites with state=0 are recognized
diff --git a/ss.cpp b/ss.cpp
index 858e4af..085a74e 100644
--- a/ss.cpp
+++ b/ss.cpp
@@ -59,7 +59,8 @@ void run(RunParameters_SS &r, unsigned seed) {
 
     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;
 
@@ -93,413 +94,6 @@ void run(RunParameters_SS &r, unsigned seed) {
     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
 }
 
@@ -580,35 +174,6 @@ void simulate(Rxn_parray &reactions, Species_parray &species, double t_end, long
             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
 
     }
 
@@ -846,29 +411,16 @@ void usage()
 " -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();
 }
 
@@ -879,12 +431,11 @@ void usage()
  *  -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
@@ -905,7 +456,7 @@ int main(int argc, char *argv[]) {
         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]); }
@@ -918,23 +469,12 @@ int main(int argc, char *argv[]) {
         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]); }
diff --git a/ss.h b/ss.h
index 2c26221..87d4f3f 100644
--- a/ss.h
+++ b/ss.h
@@ -114,20 +114,12 @@ public:
     double mu;                  // Mutation rate
     double bh;                  // Field "inverse temperature" multiplier
     double bJ;                  // Coupling "inverse temperature" multiplier
-    int estart;                 // Epitope start position
-    int eend;                   // Epitope end position
     bool runUntilEscape;        // If true, run the simulation until the population escapes
     bool runUntilEscape_all;    // If true, run the simulation until the population escapes *all* immune pressure
     bool importState;           // If true, import state from a state file
     bool useEpitope;            // Include selective pressure on an epitope 
-    bool useRelative;           // If true, use relative energy for survival probability
     bool useVerbose;            // If true, print extra information while program is running
 
-    bool useTwoSite;            // If true, use two-site two-allele model
-    double h1;
-    double h2;
-    double J12;
-
     unsigned int write_mod;     // Write state of the population every __ generations
     
     std::vector<std::set<unsigned int> > initPop;   // Initial population sequences
@@ -135,7 +127,6 @@ public:
     std::vector<std::vector<unsigned int> >           eWT;       // Sites that are consensus (WT) in the targeted epitope
     std::vector<std::vector<unsigned int> >           eMut;      // Sites that are mutant in the targeted epitope
     double                               penalty;   // Energy penalty if sequence contains the targeted epitope
-    unsigned int numEpitopes;   // how many epitopes are targeted
     
     std::string couplingsInfile;
     std::string trajectoryOutfile;
@@ -147,7 +138,6 @@ public:
 
     PenaltyTYPE penaltyType;
 
-    
 
     std::vector<std::string> epitopefiles;
     std::vector<long> init_CTL_num;
@@ -165,28 +155,11 @@ public:
         bh = 1.0;
         bJ = 1.0;
         
-        estart = 0;
-        eend   = 0;
-        
-        penalty = 10.0;
-        penaltyType = PenaltyTOTAL;
-        
-        runUntilEscape=false;
-        runUntilEscape_all=false;
         importState=false;
-        useEpitope=false;
-        useRelative=false;
         useVerbose=false;
 
-        useTwoSite = false;
-        h1  = 0.0;
-        h2  = 0.0;
-        J12 = 0.0;
-
         write_mod = 1;
 
-        numEpitopes = 0;
-
         num_CTL_clones = 0;
         
     }