Delete unused functions for SS.

diff --git a/ham_ss.cpp b/ham_ss.cpp
index 5a0d41d..536d72a 100644 (file)
--- a/ham_ss.cpp
+++ b/ham_ss.cpp
@@ -103,94 +103,6 @@ Hamiltonian::Hamiltonian(std::string &FILENAME) {
 }
 
 
-//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
 
@@ -216,289 +128,3 @@ double Hamiltonian::get_energy(const std::set<unsigned int> &mutated_sites) cons
     
 }
 
-
-// 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
index 851ac6f..0e30271 100644 (file)
--- a/ham_ss.h
+++ b/ham_ss.h
@@ -34,10 +34,6 @@ public:
     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; }
@@ -45,44 +41,4 @@ public:
 
 };
 
-
-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 3d625d0..8e43690 100644 (file)
--- a/pop_ss.cpp
+++ b/pop_ss.cpp
@@ -44,49 +44,6 @@ void Virus::print_parameters() {
 }
 
 
-//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
 

diff --git a/pop_ss.h b/pop_ss.h
index 02b4483..22ec9a5 100644 (file)
--- a/pop_ss.h
+++ b/pop_ss.h
@@ -24,12 +24,7 @@ public:
     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;