Start framework for targeting multiple epitopes.

diff --git a/hamiltonian.cpp b/hamiltonian.cpp
index b76615d..af10101 100644 (file)
--- a/hamiltonian.cpp
+++ b/hamiltonian.cpp
@@ -223,39 +223,55 @@ double Hamiltonian::get_energy(const std::set<unsigned int> &mutated_sites) cons
 
 void EpitopeHamiltonian::set_epitope(const std::vector<unsigned int> &eWT, const std::vector<unsigned int> &eMut, double d) {
 
-    penalty=d;
-    epitopeWT=eWT;
-    epitopeMut=eMut;
+    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
 
-bool EpitopeHamiltonian::escaped(const Virus &v) {
+bool EpitopeHamiltonian::escaped(const Virus &v, const std::vector<unsigned int> &eWT, const std::vector<unsigned int> &eMut) {
 
-    return EpitopeHamiltonian::escaped(v.mutated_sites);
+    return EpitopeHamiltonian::escaped(v.mutated_sites, eWT, eMut);
 
 }
 
 
 // Return true if the sequence has escaped from immune pressure
 
-bool EpitopeHamiltonian::escaped(const std::set<unsigned int> &mutated_sites) const {
+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<epitopeWT.size();i++) {
+    for (unsigned i=0;i<eWT.size();i++) {
     
-        if (mutated_sites.count(epitopeWT[i])==1) { escape=true; break; }
+        if (mutated_sites.count(eWT[i])==1) { escape=true; break; }
         
     }
     
     if (!escape) {
     
-        for (unsigned i=0;i<epitopeMut.size();i++) {
+        for (unsigned i=0;i<eMut.size();i++) {
         
-            if (mutated_sites.count(epitopeMut[i])==0) { escape=true; break; }
+            if (mutated_sites.count(eMut[i])==0) { escape=true; break; }
             
         }
         
@@ -288,8 +304,16 @@ double EpitopeHamiltonian::get_energy(const std::set<unsigned int> &mutated_site
     
     // Check for mismatch between targeted epitope and current sequence
     
-    if (escaped(mutated_sites)) return ((efield * bh) + (ecoupling * bJ));
-    else                        return ((efield * bh) + (ecoupling * bJ) + (penalty * bh));
+    double energy = (efield * bh) + (ecoupling * bJ);
+
+    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;
     
 }
 

diff --git a/hamiltonian.h b/hamiltonian.h
index 197f8ba..32ca26c 100644 (file)
--- a/hamiltonian.h
+++ b/hamiltonian.h
@@ -43,9 +43,9 @@ class EpitopeHamiltonian : public Hamiltonian {
 
 public:
     
-    double penalty;
-    std::vector<unsigned int> epitopeWT;
-    std::vector<unsigned int> epitopeMut;
+    std::vector<double> penalty;
+    std::vector<std::vector<unsigned int> > epitopeWT;
+    std::vector<std::vector<unsigned int> > epitopeMut;
     
     EpitopeHamiltonian(std::string &FILENAME);
     
@@ -53,6 +53,8 @@ public:
     
     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;
     double get_energy(const std::set<unsigned int> &mutated_sites) const;
 
 };

diff --git a/wf.cpp b/wf.cpp
index 71ecb24..bc23b18 100644 (file)
--- a/wf.cpp
+++ b/wf.cpp
@@ -299,7 +299,7 @@ int main(int argc, char *argv[]) {
         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],"-s")==0) { if (++i==argc) break; else { r.statefile=argv[i]; r.importState=true; } }
+        else if (strcmp(argv[i],"-s")==0) { if (++i==argc) break; else { r.statefile=argv[i]; r.importState=true; r.useEpitope=true; } }
         
         else if (strcmp(argv[i],"-n")==0)  { if (++i==argc) break; else r.n=(unsigned int)strtodouble(argv[i]); }
         else if (strcmp(argv[i],"-g")==0)  { if (++i==argc) break; else r.g=(unsigned int)strtodouble(argv[i]); }