From 07536fb80d5ef156f6d4902f33c5bf0aa01b1208 Mon Sep 17 00:00:00 2001
From: Dariusz Murakowski <murakdar@mit.edu>
Date: Sat, 26 Apr 2014 18:56:23 -0400
Subject: [PATCH] Command line parameter for bonus/penalty applied to each
 epitope or total immune pressure.

---
 wf.cpp | 21 ++++++++++++++++++---
 wf.h   |  5 +++++
 2 files changed, 23 insertions(+), 3 deletions(-)

diff --git a/wf.cpp b/wf.cpp
index 929cd18..9994584 100644
--- a/wf.cpp
+++ b/wf.cpp
@@ -96,12 +96,21 @@ void run(RunParameters &r, unsigned seed) {
         
         if (r.useEpitope) {
         
-            for (unsigned int n=0;n<r.num_runs;n++) {
+            for (unsigned int n=0; n<r.num_runs; n++) {
             
                 EpitopeHamiltonian H(r.couplingsInfile);
-                for (unsigned ep=0; ep<r.numEpitopes; ++ep)
-                    H.set_epitope(r.eWT[ep], r.eMut[ep], r.penalty/(double)r.numEpitopes);
+
+                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);
 
                 // print epitopes
@@ -344,6 +353,9 @@ void usage()
 "                      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"
+"                              bonus per epitope = this / number of epitopes\n"
 ;   std::cout.flush();
 }
 
@@ -387,6 +399,9 @@ int main(int argc, char *argv[]) {
         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],"-bh")==0) { if (++i==argc) break; else r.bh=strtodouble(argv[i]);              }
         else if (strcmp(argv[i],"-bJ")==0) { if (++i==argc) break; else r.bJ=strtodouble(argv[i]);              }
+
+        else if (strcmp(argv[i],"-penaltyEach")==0) { if (++i==argc) break; else { r.penalty=strtodouble(argv[i]); r.penaltyType=RunParameters::PenaltyEACH; } }
+        else if (strcmp(argv[i],"-penaltyTotal")==0) { if (++i==argc) break; else { r.penalty=strtodouble(argv[i]); r.penaltyType=RunParameters::PenaltyTOTAL; } }
         
         else if (strcmp(argv[i],"-r")==0)   { r.useRelative=true;    }
         else if (strcmp(argv[i],"-esc")==0) { r.runUntilEscape=true; }
diff --git a/wf.h b/wf.h
index d849c63..a2192d7 100644
--- a/wf.h
+++ b/wf.h
@@ -114,6 +114,10 @@ public:
     std::string supplementaryOutfile;
     std::string stateInfile;
     std::string epitopeInfile;
+
+    enum PenaltyTYPE { PenaltyEACH, PenaltyTOTAL };
+
+    PenaltyTYPE penaltyType;
     
     
     RunParameters() {
@@ -131,6 +135,7 @@ public:
         eend   = 0;
         
         penalty = 10.0;
+        penaltyType = PenaltyTOTAL;
         
         runUntilEscape=false;
         runUntilEscape_all=false;
-- 
2.7.4