From: Dariusz Murakowski <murakdar@mit.edu>
Date: Wed, 8 Apr 2015 20:48:00 +0000 (-0400)
Subject: Start of Gillespie test.
X-Git-Url: http://src.murakowski.org/?a=commitdiff_plain;h=37ad5b49f42764af968eeadd2bcabdd12e1587e3;p=VirEvoDyn.git

Start of Gillespie test.
---

diff --git a/Makefile b/Makefile
index d2adedf..7557a2f 100644
--- a/Makefile
+++ b/Makefile
@@ -5,7 +5,7 @@ EXECNAME_WF = wf
 SRCS_MC = mc.cpp hamiltonian.cpp population.cpp virus.cpp
 EXECNAME_MC = mc
 
-SRCS_SS = ss.cpp hamiltonian.cpp population.cpp virus.cpp
+SRCS_SS = ss.cpp hamiltonian.cpp population.cpp virus.cpp pop_ss.cpp reaction.cpp
 EXECNAME_SS = ss
 
 
diff --git a/pop_ss.cpp b/pop_ss.cpp
new file mode 100644
index 0000000..dc480d8
--- /dev/null
+++ b/pop_ss.cpp
@@ -0,0 +1,387 @@
+#include <vector>
+#include <iostream>
+#include <string>
+#include <fstream>
+#include <gsl/gsl_rng.h>
+#include <gsl/gsl_randist.h>
+#include <cmath>
+
+#include "pop_ss.h"
+
+
+//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
new file mode 100644
index 0000000..22a266a
--- /dev/null
+++ b/pop_ss.h
@@ -0,0 +1,74 @@
+#ifndef POP_SS_H
+#define POP_SS_H
+
+#include <vector>
+#include <map>
+#include <string>
+#include <gsl/gsl_rng.h>
+
+#include "virus.h"
+
+
+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
+
+class EpitopeRecognizer {
+public:
+    std::vector<unsigned int> epitopeWT;
+    std::vector<unsigned int> epitopeMut;
+};
+
+class Species {
+public:
+    long count;
+    std::string name;
+};
+
+class SimpleSpecies : public Species {
+};
+
+class VirusSpecies : public Species {
+public:
+    virus_map pop;
+};
+
+class CTLSpecies : public Species {
+public:
+    std::vector<EpitopeRecognizer> epitopes;
+};
+    
+
+#endif  // POP_SS_H
diff --git a/reaction.cpp b/reaction.cpp
new file mode 100644
index 0000000..077b33c
--- /dev/null
+++ b/reaction.cpp
@@ -0,0 +1,64 @@
+
+#include "reaction.h"
+
+//static long factorial[] = {1, 1, 2, 6, 24, 120, 720};
+
+double ElementaryReaction::recalc()
+{
+    //double denom = 1.0;
+    this->propensity = this->rate_constant;
+
+    std::vector<SimpleSpecies*>::iterator reac = this->reactants.begin(),
+                                          reac_end = this->reactants.end();
+    std::vector<int>::iterator reacN = this->reactant_stoich.begin(),
+                               reacN_end = this->reactant_stoich.end();
+
+    for (; (reac != reac_end) && (reacN != reacN_end); ++reac, ++reacN) {
+        long n;
+        for (n = 0; n < *reacN; ++n) {
+            this->propensity *= (*reac)->count - n;
+        }
+        //denom *= factorial[n];
+    }
+    //this->propensity /= denom;
+
+    return this->propensity;
+}
+
+void ElementaryReaction::fire()
+{
+    std::vector<SimpleSpecies*>::iterator reac = this->reactants.begin(),
+                                          reac_end = this->reactants.end();
+    std::vector<int>::iterator reacN = this->reactant_stoich.begin(),
+                               reacN_end = this->reactant_stoich.end();
+
+    for (; (reac != reac_end) && (reacN != reacN_end); ++reac, ++reacN) {
+        //if (reac->count - reacN >= 0) {
+            (*reac)->count = (*reac)->count - *reacN;
+        //}
+    }
+
+    std::vector<SimpleSpecies*>::iterator prod = this->products.begin(),
+                                          prod_end = this->products.end();
+    std::vector<int>::iterator prodN = this->product_stoich.begin(),
+                               prodN_end = this->product_stoich.end();
+
+
+    for (; (prod != prod_end) && (prodN != prodN_end); ++prod, ++prodN) {
+        (*prod)->count = (*prod)->count + *prodN;
+    }
+}
+
+
+double VirusReaction::recalc()
+{
+    // sum stuff
+    return this->propensity;
+}
+
+void VirusReaction::fire()
+{
+    // update copy numbers
+
+}
+
diff --git a/reaction.h b/reaction.h
new file mode 100644
index 0000000..90039ae
--- /dev/null
+++ b/reaction.h
@@ -0,0 +1,45 @@
+#ifndef REACTION_H
+#define REACTION_H
+
+#include <vector>
+
+#include "pop_ss.h"
+
+class Reaction {
+public:
+    double rate_constant;  // must manually ensure it includes division by statistical degeneracy factor (in case of reaction order >1 for any species)
+    double propensity;
+    virtual double recalc() = 0;
+    virtual void fire() = 0;
+};
+
+
+class ElementaryReaction : public Reaction {
+public:
+    std::vector<SimpleSpecies*> reactants;
+    std::vector<SimpleSpecies*> products;
+    std::vector<int> reactant_stoich;   // reaction order
+    std::vector<int> product_stoich;
+    double recalc();
+    void fire();
+};
+
+
+class VirusReaction : public Reaction {
+public:
+    VirusSpecies* V;
+    double recalc();
+    void fire();
+};
+
+
+class KillingReaction : public Reaction {
+public:
+    VirusSpecies* V;
+    CTLSpecies* T;
+    double recalc();
+    void fire();
+};
+
+#endif  // REACTION_H
+
diff --git a/ss.cpp b/ss.cpp
index 862b2c8..c64ee93 100644
--- a/ss.cpp
+++ b/ss.cpp
@@ -4,13 +4,15 @@
 #include <cstdio>
 #include <cstdlib>
 #include <cmath>
+#include <numeric>  // for partial_sum
 
 #include <gsl/gsl_rng.h>
 
 #include "hamiltonian.h"
-#include "population.h"
+#include "pop_ss.h"
 #include "virus.h"
 #include "ss.h"
+#include "reaction.h"
 
 
 // generate a random number generator seed from /dev/urandom
@@ -45,6 +47,150 @@ void run(RunParameters_SS &r, unsigned seed) {
     //srand((unsigned)time(0));
     //gsl_rng_set(rnd,rand());
     gsl_rng_set(rnd,seed);
+
+#if 1
+    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;
+
+    // initialize
+    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;
+
+        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;
+        }
+
+        dt = -log(gsl_rng_uniform(rnd)) / total_propensity;
+        t += dt;
+
+
+        // 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
@@ -201,11 +347,12 @@ void run(RunParameters_SS &r, unsigned seed) {
 
     }
     
-    gsl_rng_free(rnd);  //Free up memory from random number generator
 
     fclose(popout);     // close file handles
     if(!r.useTwoSite) fclose(supout);
+#endif
 
+    gsl_rng_free(rnd);  //Free up memory from random number generator
 }