+#include <iostream>
+#include <cmath>
+#include <gsl/gsl_rng.h>
+#include <gsl/gsl_randist.h>
#include "reaction.h"
+extern gsl_rng *rnd; //pointer to random number generator state
+
//static long factorial[] = {1, 1, 2, 6, 24, 120, 720};
double ElementaryReaction::recalc()
}
}
+double RC_fun(double E)
+{
+ // sigmoid function
+ //double z = exp(-E);
+ //return z/(1.0+z);
+ return (E>0.0) ? 1.0/(1.0+exp(-E)) : exp(E)/(1.0+exp(E));
+}
double VirusReaction::recalc()
{
- // sum stuff
+ V->count = 0;
+ this->propensity = this->rate_constant;
+
+ double partial_sum = 0.0;
+ virus_map::iterator iter = this->V->pop.begin(),
+ end = this->V->pop.end();
+ for(; iter != end; ++iter) {
+ V->count += iter->second;
+ partial_sum += RC_fun(iter->first.energy) * iter->second;
+ }
+
+ this->propensity *= partial_sum;
return this->propensity;
}
void VirusReaction::fire()
{
+ // pick which virus reproduces
+ double rand = gsl_rng_uniform(rnd) * this->propensity;
+ double a_sum = 0.0;
+ virus_map::iterator iter = V->pop.begin(),
+ end = V->pop.begin();
+ for(; iter != end; ++iter) {
+ a_sum += RC_fun(iter->first.energy) * iter->second;
+ if (rand <= a_sum)
+ break;
+ }
+
+ // pick how many sites to mutate in the new virus
+ // adapted from WF Virus::mutate()
+ //double mu = 6.0e-5; // mut prob upon reproduction
+ double mu = 1.0e-1;
+ Virus v(H,mu,iter->first.mutated_sites);
+ unsigned int n = gsl_ran_binomial(rnd,mu,H.size);
+ std::set<unsigned int> sites_to_mutate;
+ while (sites_to_mutate.size() < n) {
+ sites_to_mutate.insert(gsl_rng_uniform_int(rnd,H.size));
+ }
+ for (std::set<unsigned int>::iterator it = sites_to_mutate.begin(),
+ end = sites_to_mutate.end();
+ it != end; ++it) {
+ if (v.mutated_sites.count(*it)==0) v.mutated_sites.insert(*it);
+ else v.mutated_sites.erase(*it);
+ }
+ /*
+ // this version is no good when mu is high or length is too small:
+ // it can propose duplicates, under-estimating the true extent
+ for (unsigned i = 0; i < n; ++i) {
+ unsigned int x = gsl_rng_uniform_int(rnd,H.size);
+ if (v.mutated_sites.count(x)==0) v.mutated_sites.insert(x);
+ else v.mutated_sites.erase(x);
+ }
+ */
+ v.update_energy(H);
+
// update copy numbers
+ if (V->pop.count(v)==0) V->pop[v] = 1; // new variant
+ else V->pop[v] += 1; // add to existing pool
+
+}
+
+double VirusDeathReaction::recalc()
+{
+ V->count = 0;
+ propensity = rate_constant;
+
+ double partial_sum = 0.0;
+ virus_map::iterator iter = this->V->pop.begin(),
+ end = this->V->pop.end();
+ for(; iter != end; ++iter) {
+ V->count += iter->second;
+ partial_sum += iter->second;
+ }
+
+ propensity *= partial_sum;
+ return propensity;
+}
+
+void VirusDeathReaction::fire()
+{
+ // pick which virus dies
+ double rand = gsl_rng_uniform(rnd) * this->propensity;
+ double a_sum = 0.0;
+ virus_map::iterator iter = V->pop.begin(),
+ end = V->pop.begin();
+ for (; iter != end; ++iter) {
+ a_sum += iter->second;
+ if (rand <= a_sum)
+ break;
+ }
+
+ // update copy number
+ V->pop[iter->first] -= 1;
+
+ // delete zero-population strains
+ if (V->pop[iter->first] == 0)
+ V->pop.erase(iter);
+}
+
+
+double KillingReaction::recalc()
+{
+ return this->propensity;
+}
+void KillingReaction::fire()
+{
+ // update copy numbers
}
#include "ss.h"
#include "reaction.h"
+gsl_rng *rnd; //pointer to random number generator state
// generate a random number generator seed from /dev/urandom
// borrowed from SSC src/runtime/ssc_rts.c
// 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
+ rnd = gsl_rng_alloc(gsl_rng_taus2); //pointer to random number generator state
//srand((unsigned)time(0));
//gsl_rng_set(rnd,rand());
gsl_rng_set(rnd,seed);
#if 1
+ std::string couplingsFile = "neutral_2site.j";
+ Hamiltonian H(couplingsFile);
+ //double mu = 6.0e-5;
+ double mu = 1.0e-1;
+
+ std::vector<Species*> species;
+ VirusSpecies s1; s1.count = 100;
+ s1.pop[Virus(H,mu)] = s1.count;
+ species.push_back(&s1);
+
+ std::vector<Reaction*> reactions;
+ // V -> V + V'
+ VirusReaction* r1 = new VirusReaction();
+ r1->rate_constant = 3.0;
+ r1->H = H;
+ r1->V = &s1;
+ reactions.push_back(r1);
+ // V -> 0
+ VirusDeathReaction* r2 = new VirusDeathReaction();
+ r2->rate_constant = 1.5;
+ r2->H = H;
+ r2->V = &s1;
+ reactions.push_back(r2);
+
+/*
+ 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 = 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::set<unsigned int> 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)
+ 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);
+ 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";
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();
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) {
total_propensity += (*rxn)->propensity;
}
- dt = -log(gsl_rng_uniform(rnd)) / total_propensity;
- t += dt;
-
// print copy numbers
std::cout << t << '\t';
projects / VirEvoDyn.git / commitdiff