From 015c90b67d2537e60376032233ec13ebefd70746 Mon Sep 17 00:00:00 2001
From: Dariusz Murakowski <murakdar@mit.edu>
Date: Mon, 15 Jun 2015 18:20:28 -0400
Subject: [PATCH] Remove two-site Hamiltonian.

---
 hamiltonian.cpp | 120 ----------------------------
 hamiltonian.h   |  10 ---
 population.cpp  |  37 ---------
 population.h    |   1 -
 wf.cpp          | 241 ++++++++++++++++++++------------------------------------
 wf.h            |  10 ---
 6 files changed, 87 insertions(+), 332 deletions(-)

diff --git a/hamiltonian.cpp b/hamiltonian.cpp
index 70951c8..71378f7 100644
--- a/hamiltonian.cpp
+++ b/hamiltonian.cpp
@@ -352,7 +352,6 @@ double EpitopeHamiltonian::get_energy(const std::set<unsigned int> &mutated_site
 }
 
 // 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)
 {
@@ -360,7 +359,6 @@ inline bool contains(const std::vector<T> &vec, const T &value)
 }
 
 // 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)
 {
@@ -384,121 +382,3 @@ void EpitopeHamiltonian::generate_allsites_set()
 }
 
 
-
-/***********************************************
- * 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/hamiltonian.h b/hamiltonian.h
index fbdd6ad..b694db5 100644
--- a/hamiltonian.h
+++ b/hamiltonian.h
@@ -75,14 +75,4 @@ public:
 
 };
 
-
-class TwoSiteHamiltonian : public Hamiltonian {
-
-public:
-    
-    TwoSiteHamiltonian(std::string &FILENAME);
-    TwoSiteHamiltonian(double h1, double h2, double J12);
-    
-};
-
 #endif
diff --git a/population.cpp b/population.cpp
index afc1059..6ea1166 100644
--- a/population.cpp
+++ b/population.cpp
@@ -346,40 +346,3 @@ unsigned int Population::escape_variant_all(Hamiltonian &H, std::set<unsigned in
     
 }
 
-
-
-
-// Output population to file
-// two-site two-allele
-
-void Population::write_two_site_population(FILE *output, const Hamiltonian &H, unsigned int generation) {
-
-    fprintf(output,"%d",generation);
-
-    //static const unsigned int v00_arr[] = {};
-    static const unsigned int v10_arr[] = {0};
-    static const unsigned int v01_arr[] = {1};
-    static const unsigned int v11_arr[] = {0,1};
-
-    //std::set<unsigned int> v00_set (v00_arr, v00_arr+sizeof(v00_arr)/sizeof(v00_arr[0]));
-    std::set<unsigned int> v10_set (v10_arr, v10_arr+sizeof(v10_arr)/sizeof(v10_arr[0]));
-    std::set<unsigned int> v01_set (v01_arr, v01_arr+sizeof(v01_arr)/sizeof(v01_arr[0]));
-    std::set<unsigned int> v11_set (v11_arr, v11_arr+sizeof(v11_arr)/sizeof(v11_arr[0]));
-
-    Virus v00_vir (H, mu);//, v00_set);
-    Virus v10_vir (H, mu, v10_set);
-    Virus v01_vir (H, mu, v01_set);
-    Virus v11_vir (H, mu, v11_set);
-
-    fprintf(output,"\t%d",pop[v00_vir]);
-    fprintf(output,"\t%d",pop[v10_vir]);
-    fprintf(output,"\t%d",pop[v01_vir]);
-    fprintf(output,"\t%d",pop[v11_vir]);
-    fprintf(output,"\n");
-
-    fflush(output);
-
-}
-
-
-
diff --git a/population.h b/population.h
index 358f07b..bf9c70a 100644
--- a/population.h
+++ b/population.h
@@ -28,7 +28,6 @@ public:
     void next_generation(const Hamiltonian &H, gsl_rng* r, bool useRelative, bool useVerbose);
     void write_population(std::string filepath);
     void write_population(FILE *output, unsigned int generation);
-    void write_two_site_population(FILE *output, const Hamiltonian &H, unsigned int generation);
     unsigned int compute_population_size();
     void print_population_size();
     void mutate_population(const Hamiltonian &H, gsl_rng* r);
diff --git a/wf.cpp b/wf.cpp
index 8867c13..bd72472 100644
--- a/wf.cpp
+++ b/wf.cpp
@@ -48,163 +48,127 @@ void run(RunParameters &r, unsigned seed) {
     
     r.setFiles();
     FILE *popout=fopen(r.trajectoryOutfile.c_str(),"w");        // <outfile>.dat
-    FILE *supout=r.useTwoSite ? NULL : fopen(r.supplementaryOutfile.c_str(),"w");     // <outfile>.sum
+    FILE *supout=fopen(r.supplementaryOutfile.c_str(),"w");     // <outfile>.sum
     if (popout == NULL) { perror("File error"); exit(1); }
-    if (r.useTwoSite && supout == NULL) { perror("File error"); exit(1); }
+    if (supout == NULL) { perror("File error"); exit(1); }
     
     if (r.importState) importState(r);      // <infile>.st
     if (r.useEpitope)  importEpitope(r);    // <infile>.ep
     fflush(stdout);
     
 
-    if(r.useTwoSite) {
-
-            fprintf(popout,"gen\tV00\tV10\tV01\tV11\n");
-
-            for (unsigned int n=0;n<r.num_runs;n++) {
-
-                //fprintf(popout,"gen\tV00\tV10\tV01\tV11\n");
-            
-                TwoSiteHamiltonian H(r.h1,r.h2,r.J12);
-                H.set_temp(r.bh, r.bJ);
-                Population P(H, r.n, r.mu, r.initPop, r.initFrac);
-                
-                P.write_two_site_population(popout,H,0);
-
-                unsigned int i;
-                for (i=0; i<r.g; ++i) {
-                
-                    //printf("%d\t%d\n",n,i);
-                    P.next_generation(H, rnd, r.useRelative, r.useVerbose);
-
-                    if ((i+1) % r.write_mod == 0) {
-                        P.write_two_site_population(popout,H,i+1);
-                    }
-
-                    //if (r.runUntilEscape && P.escaped(H)) break;
-
-                }
-                
+    // Run (w/ targeted epitope)
+    
+    if (r.useEpitope) {
+    
+        for (unsigned int n=0; n<r.num_runs; n++) {
+        
+            EpitopeHamiltonian H(r.couplingsInfile);
+
+            double penalty = 0.0;
+            if (r.penaltyType == RunParameters::PenaltyEACH)
+                penalty = r.penalty;
+            else if (r.penaltyType == RunParameters::PenaltyTOTAL)
+                penalty = r.penalty / (double)r.numEpitopes;
+            for (unsigned ep=0; ep<r.numEpitopes; ++ep) {
+                H.set_epitope(r.eWT[ep], r.eMut[ep], penalty);
             }
 
-    }
+            H.set_temp(r.bh, r.bJ);
 
+            Population P(H, r.n, r.mu, r.initPop, r.initFrac);
 
-    else {
-        
-        // Run (w/ targeted epitope)
-        
-        if (r.useEpitope) {
-        
-            for (unsigned int n=0; n<r.num_runs; n++) {
-            
-                EpitopeHamiltonian H(r.couplingsInfile);
-
-                double penalty = 0.0;
-                if (r.penaltyType == RunParameters::PenaltyEACH)
-                    penalty = r.penalty;
-                else if (r.penaltyType == RunParameters::PenaltyTOTAL)
-                    penalty = r.penalty / (double)r.numEpitopes;
-                for (unsigned ep=0; ep<r.numEpitopes; ++ep) {
-                    H.set_epitope(r.eWT[ep], r.eMut[ep], penalty);
+            // print epitopes
+            /*
+            std::cout << "-----------\n";
+            for (unsigned ep=0; ep<H.penalty.size(); ++ep) {
+                for (unsigned i=0; i<H.epitopeWT[ep].size(); ++i) {
+                    std::cout << H.epitopeWT[ep][i] << " ";
                 }
-
-                H.set_temp(r.bh, r.bJ);
-
-                Population P(H, r.n, r.mu, r.initPop, r.initFrac);
-
-                // print epitopes
-                /*
-                std::cout << "-----------\n";
-                for (unsigned ep=0; ep<H.penalty.size(); ++ep) {
-                    for (unsigned i=0; i<H.epitopeWT[ep].size(); ++i) {
-                        std::cout << H.epitopeWT[ep][i] << " ";
-                    }
-                    std::cout << "| ";
-                    for (unsigned i=0; i<H.epitopeMut[ep].size(); ++i) {
-                        std::cout << H.epitopeMut[ep][i] << " ";
-                    }
-                    std::cout << "\n";
+                std::cout << "| ";
+                for (unsigned i=0; i<H.epitopeMut[ep].size(); ++i) {
+                    std::cout << H.epitopeMut[ep][i] << " ";
                 }
-                fflush(stdout);
-                */
-                
-                unsigned int i;
-                for (i=0; i<r.g; ++i) {
+                std::cout << "\n";
+            }
+            fflush(stdout);
+            */
             
-                    P.next_generation(H, rnd, r.useRelative, r.useVerbose);
-
-                    if ((i+1) % r.write_mod == 0) {
-                        P.write_population(popout,i+1);
-                    }
-                
-                    if (r.runUntilEscape && P.escaped(H)) break;
+            unsigned int i;
+            for (i=0; i<r.g; ++i) {
+        
+                P.next_generation(H, rnd, r.useRelative, r.useVerbose);
 
-                    if (r.runUntilEscape_all && P.escaped_all(H)) break;
-                
+                if ((i+1) % r.write_mod == 0) {
+                    P.write_population(popout,i+1);
                 }
-                
-                std::set<unsigned int> esc_var;
-                unsigned int esc_num = 0;
-                if (!r.runUntilEscape && !r.runUntilEscape_all) esc_num = P.escape_variant(H, esc_var);  // default behavior for backwards compatibility
-                if (r.runUntilEscape_all) esc_num = P.escape_variant_all(H, esc_var);
-                if (r.runUntilEscape) esc_num = P.escape_variant(H, esc_var);
-                // should behave better if both are true
-                // but more likely to escape one than all, so that one is output
-                
-                fprintf(supout,"%d\t%d",i,esc_num);
+            
+                if (r.runUntilEscape && P.escaped(H)) break;
 
-                //for (std::set<unsigned int>::iterator iter=esc_var.begin(); iter!=esc_var.end(); ++iter) fprintf(supout,"\t%d",*iter);
+                if (r.runUntilEscape_all && P.escaped_all(H)) break;
+            
+            }
+            
+            std::set<unsigned int> esc_var;
+            unsigned int esc_num = 0;
+            if (!r.runUntilEscape && !r.runUntilEscape_all) esc_num = P.escape_variant(H, esc_var);  // default behavior for backwards compatibility
+            if (r.runUntilEscape_all) esc_num = P.escape_variant_all(H, esc_var);
+            if (r.runUntilEscape) esc_num = P.escape_variant(H, esc_var);
+            // should behave better if both are true
+            // but more likely to escape one than all, so that one is output
+            
+            fprintf(supout,"%d\t%d",i,esc_num);
 
-                // print sequence as CSV (comma-separated)
-                fprintf(supout,"\t");
-                std::set<unsigned int>::iterator ms_iter=esc_var.begin();
-                if (ms_iter != esc_var.end()) {
-                    fprintf(supout,"%d",*ms_iter);
-                    ++ms_iter;
-                }
-                for (; ms_iter!=esc_var.end(); ++ms_iter)
-                    fprintf(supout,",%d",*ms_iter);
+            //for (std::set<unsigned int>::iterator iter=esc_var.begin(); iter!=esc_var.end(); ++iter) fprintf(supout,"\t%d",*iter);
 
-                fprintf(supout,"\n");
-                
-                fflush(supout);
-            
+            // print sequence as CSV (comma-separated)
+            fprintf(supout,"\t");
+            std::set<unsigned int>::iterator ms_iter=esc_var.begin();
+            if (ms_iter != esc_var.end()) {
+                fprintf(supout,"%d",*ms_iter);
+                ++ms_iter;
             }
+            for (; ms_iter!=esc_var.end(); ++ms_iter)
+                fprintf(supout,",%d",*ms_iter);
+
+            fprintf(supout,"\n");
             
-        }
-        
+            fflush(supout);
         
-        // Run (w/out targeted epitope)
+        }
         
-        else {
+    }
+    
+    
+    // Run (w/out targeted epitope)
+    
+    else {
 
-            for (unsigned int n=0;n<r.num_runs;n++) {
+        for (unsigned int n=0;n<r.num_runs;n++) {
+        
+            Hamiltonian H(r.couplingsInfile);
+            H.set_temp(r.bh, r.bJ);
+            Population P(H, r.n, r.mu, r.initPop, r.initFrac);
             
-                Hamiltonian H(r.couplingsInfile);
-                H.set_temp(r.bh, r.bJ);
-                Population P(H, r.n, r.mu, r.initPop, r.initFrac);
-                
-                unsigned int i;
-                for (i=0; i<r.g; ++i) {
-                
-                    P.next_generation(H, rnd, r.useRelative, r.useVerbose);
-                    P.write_population(popout,i+1);
+            unsigned int i;
+            for (i=0; i<r.g; ++i) {
+            
+                P.next_generation(H, rnd, r.useRelative, r.useVerbose);
+                P.write_population(popout,i+1);
 
-                    if (r.runUntilEscape && P.escaped(H)) break;
+                if (r.runUntilEscape && P.escaped(H)) break;
 
-                }
-                
             }
             
         }
-
+        
     }
-    
+
+
     gsl_rng_free(rnd);  //Free up memory from random number generator
 
     fclose(popout);     // close file handles
-    if(!r.useTwoSite) fclose(supout);
+    fclose(supout);
 
 }
 
@@ -312,20 +276,6 @@ void importEpitope(RunParameters &r) {
 }
 
 
-// insert into initial state of population
-
-void add_to_two_site_pop(RunParameters &r, bool s1, bool s2, double frac) {
-
-    r.initFrac.push_back(frac);
-
-    std::set<unsigned int> tempSet;
-    if (s1) tempSet.insert(0);
-    if (s2) tempSet.insert(1);
-    r.initPop.push_back(tempSet);
-
-}
-
-
 void usage()
 {
     std::cout <<
@@ -346,13 +296,6 @@ void usage()
 " -v                  flag for verbose output\n"
 " -numruns (int)      number of trajectories to simulate\n"
 "\n"
-" -2site              flag for two-site two-allele system\n"
-" -h1  (double)       value of field at site 1\n"
-" -h2  (double)       value of field at site 2\n"
-" -J12 (double)       value of coupling between sites 1 and 2\n"
-"                     (note the sign convention on these is P ~ exp(-H),\n"
-"                      H = sum_i h_i s_i + sum_{i<j} J_{ij} s_i s_j\n"
-"                      although internal representation is opposite)\n"
 " -write_mod  (int)   write out state every 'write_mod' generations\n"
 " -penaltyEach  (double)       bonus (per epitope) for mutation in epitope\n"
 " -penaltyTotal (double)       maximum total bonus if all epitopes contain mutations;\n"
@@ -412,16 +355,6 @@ int main(int argc, char *argv[]) {
         
         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],"-2site")==0)   { r.useTwoSite=true;     }
-        else if (strcmp(argv[i],"-h1")==0)  { if (++i==argc) break; else r.h1=-strtodouble(argv[i]);  }
-        else if (strcmp(argv[i],"-h2")==0)  { if (++i==argc) break; else r.h2=-strtodouble(argv[i]);  }
-        else if (strcmp(argv[i],"-J12")==0) { if (++i==argc) break; else r.J12=-strtodouble(argv[i]); }
-
-        else if (strcmp(argv[i],"-p00")==0) { if (++i==argc) break; else add_to_two_site_pop(r,false,false,strtodouble(argv[i])); }
-        else if (strcmp(argv[i],"-p10")==0) { if (++i==argc) break; else add_to_two_site_pop(r,true ,false,strtodouble(argv[i])); }
-        else if (strcmp(argv[i],"-p01")==0) { if (++i==argc) break; else add_to_two_site_pop(r,false,true ,strtodouble(argv[i])); }
-        else if (strcmp(argv[i],"-p11")==0) { if (++i==argc) break; else add_to_two_site_pop(r,true ,true ,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],"-write_mod")==0) { if (++i==argc) break; else r.write_mod=(unsigned int)strtodouble(argv[i]); }
diff --git a/wf.h b/wf.h
index a2192d7..94adaa9 100644
--- a/wf.h
+++ b/wf.h
@@ -95,11 +95,6 @@ public:
     bool useRelative;           // If true, use relative energy for survival probability
     bool useVerbose;            // If true, print extra information while program is running
 
-    bool useTwoSite;            // If true, use two-site two-allele model
-    double h1;
-    double h2;
-    double J12;
-
     unsigned int write_mod;     // Write state of the population every __ generations
     
     std::vector<std::set<unsigned int> > initPop;   // Initial population sequences
@@ -144,11 +139,6 @@ public:
         useRelative=false;
         useVerbose=false;
 
-        useTwoSite = false;
-        h1  = 0.0;
-        h2  = 0.0;
-        J12 = 0.0;
-
         write_mod = 1;
 
         numEpitopes = 0;
-- 
2.7.4