From: Dariusz Murakowski <murakdar@mit.edu>
Date: Thu, 16 Jul 2015 03:30:09 +0000 (-0400)
Subject: Optionally track relative contribution to energy from fields and couplings during... 
X-Git-Url: http://src.murakowski.org/?a=commitdiff_plain;h=ed0c1445b46e274362c9a811b37c92b6e4f5e5df;p=VirEvoDyn.git

Optionally track relative contribution to energy from fields and couplings during MC sampling.
---

diff --git a/mainQEE.cpp b/mainQEE.cpp
index 6b67dbe..3cd263b 100755
--- a/mainQEE.cpp
+++ b/mainQEE.cpp
@@ -132,6 +132,8 @@ int main(int argc, char *argv[]) {
         else if (strcmp(argv[i],"-p3")==0)  { r.saveP3=true;                                                                }
         else if (strcmp(argv[i],"-pk")==0)  { r.savePk=true;                                                                }
         else if (strcmp(argv[i],"-v")==0)   { r.useVerbose=true;                                                            }
+
+        else if (strcmp(argv[i],"-trackE")==0)   { r.trackE=true;                                                           }
         
         else printf("Unrecognized command! '%s'\n",argv[i]);
         
diff --git a/monteCarlo.cpp b/monteCarlo.cpp
index 9b34b80..c8b1d86 100644
--- a/monteCarlo.cpp
+++ b/monteCarlo.cpp
@@ -474,7 +474,15 @@ void getSample(const Vector &expJ, const IntVector &neighbor, MT::MersenneTwist
 
 // Get Monte Carlo sample with some sites fixed
 
-void getFixedSample(const Vector &expJ, const IntVector &neighbor, const std::vector<int> &fixedSites, MT::MersenneTwist &mt, int N, int numSteps, int stepSize, std::vector<int> &lattice, std::vector<int> &nonzero, Vector &p, Vector &expH, std::vector<double> &sumExpH) {
+void getFixedSample(const Vector &expJ, const IntVector &neighbor, const std::vector<int> &fixedSites, MT::MersenneTwist &mt, int N, int numSteps, int stepSize, std::vector<int> &lattice, std::vector<int> &nonzero, Vector &p, Vector &expH, std::vector<double> &sumExpH, const Vector &J, bool trackEnergy, FILE *EOut) {
+
+    double Efield = 0.0;
+    double Ecoupling = 0.0;
+    double Etotal = 0.0;
+    if (trackEnergy) {
+        Etotal = calculateEnergy(J, lattice, Efield, Ecoupling);
+        fprintf(EOut,"%.6e\t%.6e\t%.6e\n",Etotal,Efield,Ecoupling);
+    }
 
     for (int n=0;n<numSteps;n++) {
     
@@ -488,6 +496,11 @@ void getFixedSample(const Vector &expJ, const IntVector &neighbor, const std::ve
         }
         
         updateCorrelations(lattice, nonzero, p);
+
+        if (trackEnergy) {
+            Etotal = calculateEnergy(J, lattice, Efield, Ecoupling);
+            fprintf(EOut,"%.6e\t%.6e\t%.6e\n",Etotal,Efield,Ecoupling);
+        }
         
         // Cut off if too many mutations
         int mutCut = (N/2 < 80) ? N/2 : 120;
@@ -807,9 +820,40 @@ void getErrorMCLearn(const Vector &q, const Vector &J, Vector &expJ, double B, i
 }
 
 
+// calculate energy due to fields and couplings
+double calculateEnergy(const Vector &J, const std::vector<int> &config, double &Efield, double &Ecoupling)
+{
+    Efield = 0;
+    Ecoupling = 0;
+    
+    for (unsigned i=0;i<config.size();i++) {
+
+        //assert(config[i] <= J[i].size());
+        
+        if (config[i]!=J[i].size()) {
+        
+            Efield -= J.at(i).at(config[i]);
+            
+            for (unsigned j=i+1;j<config.size();j++) {
+            
+                if (config[j]!=J[j].size())
+                    Ecoupling -= J.at(index(i,j,config.size())).at(sindex(config[i],config[j],J[i].size(),J[j].size()));
+                
+            }
+            
+        }
+        
+    }
+    
+    double E = Efield + Ecoupling;
+
+    return E;
+}
+
+
 // Compute correlations using rejection-free Monte Carlo (getError_RF) with epitope sites fixed
 
-void getEpitopeCorrelations(const Vector &J, const Vector &expJ, const std::vector<std::vector<int> > &epitopeStates, const std::vector<int> &epitopeStartSite, int numSteps, int numRuns, Vector &p, std::vector<int> &latticeStart) {
+void getEpitopeCorrelations(const Vector &J, const Vector &expJ, const std::vector<std::vector<int> > &epitopeStates, const std::vector<int> &epitopeStartSite, int numSteps, int numRuns, Vector &p, std::vector<int> &latticeStart, bool trackEnergy, const std::string &energyOutFileName) {
 
     // Create and initialize Monte Carlo variables
 
@@ -844,6 +888,8 @@ void getEpitopeCorrelations(const Vector &J, const Vector &expJ, const std::vect
     // MC variables
 
     int stepSize = 4*N;
+    //if (trackEnergy)
+    //    stepSize = N;
     
     // Set variables for parallel
     
@@ -855,6 +901,13 @@ void getEpitopeCorrelations(const Vector &J, const Vector &expJ, const std::vect
 
     MT::MersenneTwist mt;
     mt.init_genrand(seed[0]);
+
+    // possibly track energy
+    FILE *EOut = NULL;
+    if (trackEnergy) {
+        EOut = fopen(energyOutFileName.c_str(),"w");
+        if (EOut == NULL) { perror("File error"); exit(1); }
+    }
         
     // Thermalize and get initial error
         
@@ -883,7 +936,7 @@ void getEpitopeCorrelations(const Vector &J, const Vector &expJ, const std::vect
         } }
         
         setEnergies(expJ, neighbor, lattice, nonzero, expH, sumExpH);
-        getFixedSample(expJ, neighbor, fixedSites, mt, N, numSteps, stepSize, lattice, nonzero, p, expH, sumExpH);
+        getFixedSample(expJ, neighbor, fixedSites, mt, N, numSteps, stepSize, lattice, nonzero, p, expH, sumExpH, J, trackEnergy, EOut);
             
     }
         
@@ -891,6 +944,9 @@ void getEpitopeCorrelations(const Vector &J, const Vector &expJ, const std::vect
     
     for (int i=0;i<p.size();i++) { for (int j=0;j<p[i].size();j++) p[i][j] /= (double) (numThreads * numSteps * numRuns); }
 
+    if (trackEnergy)
+        fclose(EOut);
+
 }
 
 
diff --git a/monteCarlo.h b/monteCarlo.h
index 2eb5c8d..1bfce0e 100755
--- a/monteCarlo.h
+++ b/monteCarlo.h
@@ -22,7 +22,7 @@ void dynamicsRF(const Vector &, const IntVector &, int, double, std::vector<int>
 void updateCorrelations(const std::vector<int> &, const std::vector<int> &, Vector &);
 void updateCorrelationsP3(const std::vector<int> &, const std::vector<int> &, Vector &, std::vector<std::vector<std::vector<std::vector<double> > > > &);
 void getSample(const Vector &, const IntVector &, MT::MersenneTwist &, int, int, int, std::vector<int> &, std::vector<int> &, Vector &, Vector &, std::vector<double> &);
-void getFixedSample(const Vector &, const IntVector &, const std::vector<int> &, MT::MersenneTwist &, int, int, int, std::vector<int> &, std::vector<int> &, Vector &, Vector &, std::vector<double> &);
+void getFixedSample(const Vector &, const IntVector &, const std::vector<int> &, MT::MersenneTwist &, int, int, int, std::vector<int> &, std::vector<int> &, Vector &, Vector &, std::vector<double> &, const Vector &, bool, FILE *);
 void getFixedSampleP3Pk(const Vector &, const IntVector &, const std::vector<int> &, MT::MersenneTwist &, int, int, int, std::vector<int> &, std::vector<int> &, Vector &, std::vector<std::vector<std::vector<std::vector<double> > > > &, std::vector<double> &, Vector &, std::vector<double> &);
 void writeFixedSample(FILE *, const Vector &, const IntVector &, const std::vector<int> &, MT::MersenneTwist &, int, int, int, std::vector<int> &, std::vector<int> &, Vector &, Vector &, std::vector<double> &);
 
@@ -33,9 +33,10 @@ void getNeighbors(const Vector &, int, double, IntVector &);
 // Control functions
 void getError(const Vector &, const Vector &, int, double, int, int, double, double, std::vector<double> &);
 void getErrorMCLearn(const Vector &, const Vector &, Vector &, double, int, int, double, double, Vector &, std::vector<double> &, std::vector<int> &);
-void getEpitopeCorrelations(const Vector &, const Vector &, const std::vector<std::vector<int> > &, const std::vector<int> &, int, int, Vector &, std::vector<int> &);
+void getEpitopeCorrelations(const Vector &, const Vector &, const std::vector<std::vector<int> > &, const std::vector<int> &, int, int, Vector &, std::vector<int> &, bool, const std::string &);
 void getEpitopeCorrelationsP3Pk(const Vector &, const Vector &, const std::vector<std::vector<int> > &, const std::vector<int> &, int, int, Vector &, std::vector<std::vector<std::vector<std::vector<double> > > > &, std::vector<double> &, std::vector<int> &);
 void getEpitopeStates(FILE *, const Vector &, const Vector &, const std::vector<std::vector<int> > &, const std::vector<int> &, int, int, Vector &, std::vector<int> &);
 
+double calculateEnergy(const Vector &, const std::vector<int> &, double &, double &);
 
 #endif
diff --git a/qEpitopeEval.cpp b/qEpitopeEval.cpp
index c1472dd..4807def 100755
--- a/qEpitopeEval.cpp
+++ b/qEpitopeEval.cpp
@@ -92,7 +92,7 @@ void runEpitopeEval(RunParameters &r, std::vector<std::vector<int> > &eStates, s
     // Run MC and get error
     
     if (r.saveP3) getEpitopeCorrelationsP3Pk(J, expJ, eStates, eStartSite, r.b, r.runs, p, p3, pk, lattice);
-    else          getEpitopeCorrelations(J, expJ, eStates, eStartSite, r.b, r.runs, p, lattice);
+    else          getEpitopeCorrelations(J, expJ, eStates, eStartSite, r.b, r.runs, p, lattice, r.trackE, r.getEnergyOutFile());
     
     // Write out correlations
     
diff --git a/qEpitopeEval.h b/qEpitopeEval.h
index fcc4c67..a4d3e24 100755
--- a/qEpitopeEval.h
+++ b/qEpitopeEval.h
@@ -34,6 +34,8 @@ public:
     bool saveP2;                    // If true, record pair correlations in additions to one-point correlations
     bool saveP3;                    // If true, record three-point correlations
     bool savePk;                    // If true, record P(k) mutations curve
+
+    bool trackE;                    // if true, save relative contribution to energy from fields and couplings
     
     RunParameters() {
         
@@ -49,6 +51,8 @@ public:
         saveP2=false;
         saveP3=false;
         savePk=false;
+
+        trackE = false;
         
     }
     std::string getInfile()                   { return (directory+"/"+infile+".j");                   }
@@ -59,6 +63,7 @@ public:
     std::string getPkOutfile()                { return (directory+"/"+outfile+".pk");                 }
     std::string getPkOutfile(int i)           { return (directory+"/"+outfile+"-"+tostring(i)+".pk"); }
     std::string getStartInfile()              { return (directory+"/"+startfile+".dat");              }
+    std::string getEnergyOutFile()            { return (directory+"/"+outfile+".E");                 }
     ~RunParameters() {}
     
 };