Clarify dependencies. Copy hamiltonian.{cpp,h} to ham_ss.{cpp,h}. Copy contents of...

diff --git a/Makefile b/Makefile
index 7557a2f..a96a8b0 100644 (file)
--- 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 pop_ss.cpp reaction.cpp
+SRCS_SS = ss.cpp ham_ss.cpp pop_ss.cpp reaction.cpp
 EXECNAME_SS = ss
 
 

diff --git a/ham_ss.cpp b/ham_ss.cpp
new file mode 100644 (file)
index 0000000..5a0d41d
--- /dev/null
+++ b/ham_ss.cpp
@@ -0,0 +1,504 @@
+#include <iostream>
+#include <string>
+#include <fstream>
+#include <cstdio>
+#include <cstdlib>
+#include <map>
+#include <vector>
+
+#include "ham_ss.h"
+
+
+//Constructor for loading couplings from John Barton's Ising Inversion 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
+    if (input == NULL) { perror((std::string("ERROR: ") + FILENAME).c_str()); exit(1); }
+       
+       typedef std::vector<int> Key;
+       std::map<Key,double> JIndex;
+    char o;
+       int site=0;
+    
+       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) {
+                    
+                                       inputKey[1]=site;
+                    JIndex[inputKey]=neighborJ;
+                                       
+                               }
+                               else {
+                
+                        inputKey[1]=neighborSite;
+                        JIndex[inputKey]=neighborJ;
+                        
+                }
+                //std::cout << site << " " << neighborSite << " " << neighborJ << "\n";
+                
+                       }
+            
+               }
+        
+       }
+               
+       fclose(input);
+       Key inputKey(2,0);
+       J.resize(size);
+    
+       for(int i=0;i<size;i++) J[i].resize(size,0.0);
+    Graph.resize(size);
+    
+       for (int i=0;i<size;i++) {
+               
+        J[i].resize(size,0.0);
+        
+               for(int j=i;j<size;j++){
+        
+                       inputKey[0]=i;
+                       inputKey[1]=j;
+            
+                       if (JIndex[inputKey]!=0.0) {
+            
+                               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);
+                    
+                               }
+                
+                       }
+            
+               }
+        
+       }
+
+    //for(int i=0; i<size; i++)
+        //for(int j=0; j<size; j++)
+            //std::cout << i << " " << j << " " << J[i][j] << "\n";
+            //printf("%d %d %lf\n",i,j,J[i][j]);
+    
+}
+
+
+//Constructor for loading couplings from John Barton's Ising Inversion code
+//Code for this function derived from John Barton's code       
+
+EpitopeHamiltonian::EpitopeHamiltonian(std::string &FILENAME) {
+       
+    this->bh=1.0;
+    this->bJ=1.0;
+    
+       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;
+    char o;
+       int site=0;
+    
+       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) {
+                    
+                                       inputKey[1]=site;
+                    JIndex[inputKey]=neighborJ;
+                                       
+                               }
+                               else {
+                
+                        inputKey[1]=neighborSite;
+                        JIndex[inputKey]=neighborJ;
+                        
+                }
+                
+                       }
+            
+               }
+        
+       }
+               
+       fclose(input);
+       Key inputKey(2,0);
+       J.resize(size);
+    
+       for(int i=0;i<size;i++) J[i].resize(size,0.0);
+    Graph.resize(size);
+    
+       for (int i=0;i<size;i++) {
+               
+        J[i].resize(size,0.0);
+        
+               for(int j=i;j<size;j++){
+        
+                       inputKey[0]=i;
+                       inputKey[1]=j;
+            
+                       if (JIndex[inputKey]!=0.0) {
+            
+                               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);
+                    
+                               }
+                
+                       }
+            
+               }
+        
+       }
+    
+}
+
+
+
+// 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;
+    
+       for (iter1=mutated_sites.begin();iter1!=mutated_sites.end();++iter1) {
+    
+               efield -= J[*iter1][*iter1];
+               
+        iter2 = iter1;
+        ++iter2;
+               
+        for (;iter2!=mutated_sites.end();iter2++) ecoupling -= J[*iter1][*iter2];
+                
+    }
+    
+    return ((efield * bh) + (ecoupling * bJ));
+    
+}
+
+
+// Set the targeted epitope
+
+void EpitopeHamiltonian::set_epitope(const std::vector<unsigned int> &eWT, const std::vector<unsigned int> &eMut, double d) {
+
+    penalty.push_back(d);
+    epitopeWT.push_back(eWT);
+    epitopeMut.push_back(eMut);
+
+}
+
+
+// Did virus escape from *any* immune pressure anywhere?
+bool EpitopeHamiltonian::escaped(const Virus &v) {
+    return EpitopeHamiltonian::escaped(v.mutated_sites);
+}
+
+// Did sequence escape from *any* immune pressure anywhere?
+bool EpitopeHamiltonian::escaped(const std::set<unsigned int> &mutated_sites) const {
+    for (unsigned ep=0; ep<penalty.size(); ++ep) {
+        if (escaped(mutated_sites, epitopeWT[ep], epitopeMut[ep])) {
+            return true;
+        }
+    }
+    return false;
+}
+
+
+// Return true if the virus has escaped from immune pressure at particular epitope
+
+bool EpitopeHamiltonian::escaped(const Virus &v, const std::vector<unsigned int> &eWT, const std::vector<unsigned int> &eMut) {
+
+    return EpitopeHamiltonian::escaped(v.mutated_sites, eWT, eMut);
+
+}
+
+
+// Return true if the sequence has escaped from immune pressure at particular epitope
+
+bool EpitopeHamiltonian::escaped(const std::set<unsigned int> &mutated_sites, const std::vector<unsigned int> &eWT, const std::vector<unsigned int> &eMut) const {
+
+    bool escape=false;
+    
+    for (unsigned i=0;i<eWT.size();i++) {
+    
+        if (mutated_sites.count(eWT[i])==1) { escape=true; break; }
+        
+    }
+    
+    if (!escape) {
+    
+        for (unsigned i=0;i<eMut.size();i++) {
+        
+            if (mutated_sites.count(eMut[i])==0) { escape=true; break; }
+            
+        }
+        
+    }
+    
+    return escape;
+
+}
+
+
+// Did virus escape from *all* immune pressure everywhere?
+bool EpitopeHamiltonian::escaped_all(const Virus &v) {
+    return EpitopeHamiltonian::escaped_all(v.mutated_sites);
+}
+
+// Did sequence escape from *all* immune pressure everywhere?
+bool EpitopeHamiltonian::escaped_all(const std::set<unsigned int> &mutated_sites) const {
+    for (unsigned ep=0; ep<penalty.size(); ++ep) {
+        if (!escaped(mutated_sites, epitopeWT[ep], epitopeMut[ep])) {
+            return false;
+        }
+    }
+    return true;
+}
+
+
+
+// Return the energy for a given set of mutated sites
+// default: do include contribution from epitopes
+
+double EpitopeHamiltonian::get_energy(const std::set<unsigned int> &mutated_sites) const
+{
+    return get_energy<true>(mutated_sites);
+}
+
+// tried to have a forward declaration of explicit specializations,
+// which "must occur before instantiation",
+// but then linker complained about undefined reference
+double ensure_get_energy_false_instantiated(const std::set<unsigned int> &x) { EpitopeHamiltonian A((std::string&)""); return A.get_energy<false>(x); }
+
+// Return the energy for a given set of mutated sites
+
+template <bool include_epitope>
+double EpitopeHamiltonian::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;
+    
+       for (iter1=mutated_sites.begin();iter1!=mutated_sites.end();++iter1) {
+    
+               efield -= J[*iter1][*iter1];
+               
+        iter2 = iter1;
+        ++iter2;
+               
+        for (;iter2!=mutated_sites.end();iter2++) ecoupling -= J[*iter1][*iter2];
+                
+    }
+    
+    double energy = (efield * bh) + (ecoupling * bJ);
+
+    if (include_epitope) {
+
+        // Check for mismatch between targeted epitope and current sequence
+        
+        for (unsigned ep=0; ep<penalty.size(); ++ep) {
+            if (!escaped(mutated_sites, epitopeWT[ep], epitopeMut[ep])) {
+                energy += penalty[ep] * bh;
+            }
+            // alternatively, if(escaped): energy -= penalty
+        }
+
+    }
+
+    return energy;
+    
+}
+
+// helper function to test whether vector contains a value
+// based on http://stackoverflow.com/questions/6194797/what-is-a-c-container-with-a-contains-operation
+template <class T>
+inline bool contains(const std::vector<T> &vec, const T &value)
+{
+        return std::find(vec.begin(), vec.end(), value) != vec.end();
+}
+
+// helper function to test whether set contains a value
+// based on http://stackoverflow.com/questions/6194797/what-is-a-c-container-with-a-contains-operation
+template <class T>
+inline bool contains(const std::set<T> &container, const T &value)
+{
+        return container.find(value) != container.end();
+}
+
+void EpitopeHamiltonian::generate_allsites_set()
+{
+    this->epitopeAllSites = MutatedSiteSequence();
+    for (unsigned ep=0; ep<penalty.size(); ++ep) {
+        epitopeAllSites.insert(epitopeWT[ep].begin(), epitopeWT[ep].end());
+        epitopeAllSites.insert(epitopeMut[ep].begin(), epitopeMut[ep].end());
+    }
+
+    this->notEpitopeAllSites = MutatedSiteSequence();
+    for (unsigned i = 0; i < static_cast<unsigned int>(this->size); ++i) {
+        if (!contains(epitopeAllSites,i)) {
+            notEpitopeAllSites.insert(i);
+        }
+    }
+}
+
+
+
+/***********************************************
+ * my simple 2-site 2-allele system
+ ***********************************************/
+
+// Return the energy for a given set of mutated sites
+
+TwoSiteHamiltonian::TwoSiteHamiltonian(double h1, double h2, double J12) {
+
+    this->bh=1.0;
+    this->bJ=1.0;
+
+    this->size = 2;
+
+    J.resize(size);
+
+    for(int i=0; i<size; ++i) {
+        J[i].resize(size);
+    }
+
+    J[0][0] = h1;
+    J[1][1] = h2;
+    J[0][1] = J12;
+    J[1][0] = J[0][1];
+
+}
+
+TwoSiteHamiltonian::TwoSiteHamiltonian(std::string &FILENAME) {
+       
+    this->bh=1.0;
+    this->bJ=1.0;
+    
+       FILE* input;
+       input = fopen(FILENAME.c_str(),"r");
+    if (input == NULL) { perror((std::string("ERROR: ") + FILENAME).c_str()); exit(1); }
+       
+       typedef std::vector<int> Key;
+       std::map<Key,double> JIndex;
+    char o;
+       int site=0;
+    
+       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) {
+                    
+                                       inputKey[1]=site;
+                    JIndex[inputKey]=neighborJ;
+                                       
+                               }
+                               else {
+                
+                        inputKey[1]=neighborSite;
+                        JIndex[inputKey]=neighborJ;
+                        
+                }
+                //std::cout << site << " " << neighborSite << " " << neighborJ << "\n";
+                
+                       }
+            
+               }
+        
+       }
+               
+       fclose(input);
+       Key inputKey(2,0);
+       J.resize(size);
+    
+       for(int i=0;i<size;i++) J[i].resize(size,0.0);
+    Graph.resize(size);
+    
+       for (int i=0;i<size;i++) {
+               
+        J[i].resize(size,0.0);
+        
+               for(int j=i;j<size;j++){
+        
+                       inputKey[0]=i;
+                       inputKey[1]=j;
+            
+                       if (JIndex[inputKey]!=0.0) {
+            
+                               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);
+                    
+                               }
+                
+                       }
+            
+               }
+        
+       }
+
+    //for(int i=0; i<size; i++)
+        //for(int j=0; j<size; j++)
+            //std::cout << i << " " << j << " " << J[i][j] << "\n";
+            //printf("%d %d %lf\n",i,j,J[i][j]);
+    
+}
+
+

diff --git a/ham_ss.h b/ham_ss.h
new file mode 100644 (file)
index 0000000..851ac6f
--- /dev/null
+++ b/ham_ss.h
@@ -0,0 +1,88 @@
+#ifndef HAM_SS_H
+#define HAM_SS_H
+
+#include <string>
+#include <vector>
+#include <set>
+
+#include "pop_ss.h"
+
+
+typedef std::vector<std::vector<int> > AdjacencyList;
+typedef std::vector<std::vector<double> > Coupling;
+
+typedef std::set<unsigned int> MutatedSiteSequence;
+
+
+class Virus;
+
+//template <class T> bool contains(const std::vector<T> &vec, const T &value);
+//template <class T> bool contains(const std::set<T> &container, const T &value);
+
+
+class Hamiltonian {
+
+public:
+    
+    int size;
+    double bh;      // fields "inverse temperature" multiplier
+    double bJ;      // couplings "inverse temperature" multiplier
+    Coupling J;
+    AdjacencyList Graph;
+    
+    Hamiltonian(std::string &FILENAME);
+    Hamiltonian() { }
+    virtual ~Hamiltonian() { }
+
+    virtual bool escaped(const Virus &v) { return false; }
+    virtual bool escaped(const std::set<unsigned int> &mutated_sites) const { return false; }
+    virtual bool escaped_all(const Virus &v) { return false; }
+    virtual bool escaped_all(const std::set<unsigned int> &mutated_sites) const { return false; }
+    virtual double get_energy(const std::set<unsigned int> &mutated_sites) const;
+    
+    void set_temp(double x)           { bh=x; bJ=x; }
+    void set_temp(double x, double y) { bh=x; bJ=y; }
+
+};
+
+
+class EpitopeHamiltonian : public Hamiltonian {
+
+public:
+    
+    std::vector<double> penalty;
+    std::vector<std::vector<unsigned int> > epitopeWT;
+    std::vector<std::vector<unsigned int> > epitopeMut;
+
+    MutatedSiteSequence epitopeAllSites;    // all site indices belonging to an epitope (WT or mut ant)
+    MutatedSiteSequence notEpitopeAllSites; // all sites not involved in an epitope
+    
+    EpitopeHamiltonian(std::string &FILENAME);
+    
+    void set_epitope(const std::vector<unsigned int> &eWT, const std::vector<unsigned int> &eMut, double d);
+    void generate_allsites_set();
+    
+    bool escaped(const Virus &v);
+    bool escaped(const std::set<unsigned int> &mutated_sites) const;
+    bool escaped(const Virus &v, const std::vector<unsigned int> &eWT, const std::vector<unsigned int> &eMut);
+    bool escaped(const std::set<unsigned int> &mutated_sites, const std::vector<unsigned int> &eWT, const std::vector<unsigned int> &eMut) const;
+
+    bool escaped_all(const Virus &v);
+    bool escaped_all(const std::set<unsigned int> &mutated_sites) const;
+
+    double get_energy(const std::set<unsigned int> &mutated_sites) const;
+    template<bool include_epitope> double get_energy(const std::set<unsigned int> &mutated_sites) const;
+
+};
+
+
+class TwoSiteHamiltonian : public Hamiltonian {
+
+public:
+    
+    TwoSiteHamiltonian(std::string &FILENAME);
+    TwoSiteHamiltonian(double h1, double h2, double J12);
+    
+};
+
+#endif  // HAM_SS_H

diff --git a/pop_ss.cpp b/pop_ss.cpp
index 94fa159..3d625d0 100644 (file)
--- a/pop_ss.cpp
+++ b/pop_ss.cpp
@@ -1,4 +1,5 @@
 #include <vector>
+#include <set>
 #include <iostream>
 #include <string>
 #include <fstream>
@@ -8,6 +9,101 @@
 
 #include "pop_ss.h"
 
+
+//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;
+
+}
+
+
+//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);
+
+}              
+
+
+//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;
+
+}
+
+
+//Mutate the entire virus.  Takes a Hamiltonian object and a pointer
+//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);
+    
+    for (unsigned i=0;i<n;i++) {
+    
+        unsigned int x=gsl_rng_uniform_int(r,H.size);
+        
+        if (mutated_sites.count(x)==0) mutated_sites.insert(x);
+        else                           mutated_sites.erase(x);
+        
+    }
+        
+    update_energy(H);
+    
+}
+
+
+//Computes the odds of survival
+
+double Virus::survival() const {
+
+       double z = exp(-energy);
+       return z/(1+z);
+    
+}
+
+
+//Computes the odds of survival
+
+double Virus::survival(double Eavg) const {
+
+       double z = exp(Eavg-energy);
+       return z/(1+z);
+    
+}
+
+
+//Update the energy of the virus object.  Takes as an argument the
+//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;
+
+}
+
+
 // Did virus escape from immune pressure?
 bool EpitopeRecognizer::escaped(const Virus &v) const {
     return EpitopeRecognizer::escaped(v.mutated_sites);

diff --git a/pop_ss.h b/pop_ss.h
index 7451e4f..02b4483 100644 (file)
--- a/pop_ss.h
+++ b/pop_ss.h
@@ -2,11 +2,45 @@
 #define POP_SS_H
 
 #include <vector>
+#include <set>
 #include <map>
 #include <string>
 #include <gsl/gsl_rng.h>
 
-#include "virus.h"
+#include "ham_ss.h"
+
+
+class Hamiltonian;
+
+
+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);
+    
+    void print_parameters();
+    void mutate(const Hamiltonian &H, gsl_rng* r);
+    double survival() const;
+    double survival(double Eavg) const;
+    
+    friend class Population;
+       
+//private:
+    
+    double mu;
+    unsigned int L;
+       
+    void update_energy(const Hamiltonian &H);
+    
+};
+
+
+bool operator<(const Virus& lhs, const Virus& rhs);
 
 
 typedef std::map<Virus, unsigned int> virus_map;

diff --git a/ss.cpp b/ss.cpp
index ba25e86..739a1d5 100644 (file)
--- a/ss.cpp
+++ b/ss.cpp
@@ -9,9 +9,9 @@
 
 #include <gsl/gsl_rng.h>
 
-#include "hamiltonian.h"
+#include "ham_ss.h"
 #include "pop_ss.h"
-#include "virus.h"
+//#include "virus.h"
 #include "ss.h"
 #include "reaction.h"
 

diff --git a/ss.h b/ss.h
index 0b86e4e..ff80c08 100644 (file)
--- a/ss.h
+++ b/ss.h
@@ -1,5 +1,5 @@
-#ifndef WF_H
-#define WF_H
+#ifndef SS_H
+#define SS_H
 
 #include <vector>
 #include <set>