//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) {
+ 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++)
//std::cout << i << " " << j << " " << J[i][j] << "\n";
// Return the energy for a given set of mutated sites
double Hamiltonian::get_energy(const std::set<unsigned int> &mutated_sites) const {
-
+
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));
Hamiltonian(std::string &FILENAME);
Hamiltonian() { }
virtual ~Hamiltonian() { }
-
+
virtual double get_energy(const std::set<unsigned int> &mutated_sites) const;
void set_temp(double x) { bh=x; bJ=x; }
//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;
}
//Hamiltonian that determines the energy
void Virus::update_energy(const Hamiltonian &H) {
-
energy=H.get_energy(mutated_sites);
-
}
bool operator<(const Virus& lhs, const Virus& rhs) {
-
return lhs.mutated_sites < rhs.mutated_sites;
-
}
class Virus {
-
public:
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 mu;
unsigned int L;
-
+
void update_energy(const Hamiltonian &H);
};
+// needed for STL map
bool operator<(const Virus& lhs, const Virus& rhs);
std::vector<EpitopeRecognizer> epitopes;
std::vector<double> affinity;
size_t num_ep;
-
+
double recognized(const Virus &v) const;
double recognized(const MutatedSiteSequence &mutated_sites) const;
};
-
+
#endif // POP_SS_H
#include "reaction.h"
-extern gsl_rng *rnd; //pointer to random number generator state
+extern gsl_rng *rnd; //pointer to random number generator state
//static long factorial[] = {1, 1, 2, 6, 24, 120, 720};
if (rand <= a_sum * rate_constant)
break;
}
-
+
// pick how many sites to mutate in the new virus
// adapted from WF Virus::mutate()
//double mu = iter->first.mu;
#include "ss.h"
#include "reaction.h"
-gsl_rng *rnd; //pointer to random number generator state
+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
- 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());
break;
}
}
-
+
// TODO: better handle case with total propensity == 0.0
if (rxn == end)
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
*/
int main(int argc, char *argv[]) {
-
+
RunParameters_SS r;
unsigned seed = sim_random_seed();
-
+
// Process command line input
-
+
for (int i=1;i<argc;i++) {
-
+
if (strcmp(argv[i],"-d")==0) { if (++i==argc) break; else r.directory=argv[i]; }
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],"-epN0")==0) { if (++i==argc) break; else { r.init_CTL_numN.back() = (long)strtodouble(argv[i]); } }
else if (strcmp(argv[i],"-epM0")==0) { if (++i==argc) break; else { r.init_CTL_numM.back() = (long)strtodouble(argv[i]); } }
else if (strcmp(argv[i],"-epT0")==0) { if (++i==argc) break; else { r.init_CTL_numT.back() = (long)strtodouble(argv[i]); } }
-
+
else if (strcmp(argv[i],"-n")==0) { if (++i==argc) break; else r.n=(unsigned int)strtodouble(argv[i]); }
else if (strcmp(argv[i],"-mu")==0) { if (++i==argc) break; else r.mu=strtodouble(argv[i]); }
else if (strcmp(argv[i],"-b")==0) { if (++i==argc) break; else { r.bh=strtodouble(argv[i]); r.bJ=r.bh; } }
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],"-seed")==0) { if (++i==argc) break; else seed=(unsigned)strtodouble(argv[i]); }
else if (strcmp(argv[i],"-rate_h")==0) { if (++i==argc) break; else r.rate_h=strtodouble(argv[i]); }
else if (strcmp(argv[i],"-h")==0 || strcmp(argv[i],"--help")==0) { usage(); return 0; }
-
- else printf("Unrecognized command! '%s'\n",argv[i]);
-
+
+ else printf("Unrecognized argument! '%s'\n",argv[i]);
+
}
-
+
std::cout << "seed = " << seed << "\n";
run(r,seed);
-
+
return 0;
-
-}
+
+}
// Converts generic to string
template <class T>
-
inline std::string tostring (const T& t) {
std::stringstream ss;
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 useEpitope; // Include selective pressure on an epitope
bool useVerbose; // If true, print extra information while program is running
-
+
std::vector<std::set<unsigned int> > initPop; // Initial population sequences
std::vector<double> initFrac; // Initial population fraction
-
+
std::string couplingsInfile;
std::string trajectoryOutfile;
std::string supplementaryOutfile;
double rate_h;
+ // constructor sets default parameters
RunParameters_SS() {
directory="";
-
+
n = 100000;
num_runs = 1000;
mu = 6.0e-5;
bh = 1.0;
bJ = 1.0;
-
+
importState=false;
useVerbose=false;
}
+
void setFiles() {
-
+
if (strcmp(directory.c_str(),"")!=0) {
-
+
couplingsInfile=directory+"/"+infile+".j";
trajectoryOutfile=directory+"/"+outfile+".dat";
supplementaryOutfile=directory+"/"+outfile+".sum";
for (size_t i = 0; i != epitopefiles.size(); ++i) {
epitopefiles[i] = directory+"/"+epitopefiles[i]+".ep";
}
-
+
}
-
+
else {
-
+
couplingsInfile=infile+".j";
trajectoryOutfile=outfile+".dat";
supplementaryOutfile=outfile+".sum";
for (size_t i = 0; i != epitopefiles.size(); ++i) {
epitopefiles[i] = epitopefiles[i]+".ep";
}
-
+
}
-
+
//printf("%s\n%s\n%s\n",couplingsInfile.c_str(),trajectoryOutfile.c_str(),stateInfile.c_str());
-
+
}
+
+
~RunParameters_SS() {}
-
+
};
projects / VirEvoDyn.git / commitdiff