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

diff --git a/mainQEE.cpp b/mainQEE.cpp
index 6b67dbe..3cd263b 100755 (executable)
--- 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;                                                                }\r
         else if (strcmp(argv[i],"-pk")==0)  { r.savePk=true;                                                                }\r
         else if (strcmp(argv[i],"-v")==0)   { r.useVerbose=true;                                                            }\r
+\r
+        else if (strcmp(argv[i],"-trackE")==0)   { r.trackE=true;                                                           }\r
         \r
         else printf("Unrecognized command! '%s'\n",argv[i]);\r
         \r

diff --git a/monteCarlo.cpp b/monteCarlo.cpp
index 9b34b80..c8b1d86 100644 (file)
--- a/monteCarlo.cpp
+++ b/monteCarlo.cpp
@@ -474,7 +474,15 @@ void getSample(const Vector &expJ, const IntVector &neighbor, MT::MersenneTwist
 \r
 // Get Monte Carlo sample with some sites fixed\r
 \r
-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) {\r
+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) {\r
+\r
+    double Efield = 0.0;\r
+    double Ecoupling = 0.0;\r
+    double Etotal = 0.0;\r
+    if (trackEnergy) {\r
+        Etotal = calculateEnergy(J, lattice, Efield, Ecoupling);\r
+        fprintf(EOut,"%.6e\t%.6e\t%.6e\n",Etotal,Efield,Ecoupling);\r
+    }\r
 \r
     for (int n=0;n<numSteps;n++) {\r
     \r
@@ -488,6 +496,11 @@ void getFixedSample(const Vector &expJ, const IntVector &neighbor, const std::ve
         }\r
         \r
         updateCorrelations(lattice, nonzero, p);\r
+\r
+        if (trackEnergy) {\r
+            Etotal = calculateEnergy(J, lattice, Efield, Ecoupling);\r
+            fprintf(EOut,"%.6e\t%.6e\t%.6e\n",Etotal,Efield,Ecoupling);\r
+        }\r
         \r
         // Cut off if too many mutations\r
         int mutCut = (N/2 < 80) ? N/2 : 120;\r
@@ -807,9 +820,40 @@ void getErrorMCLearn(const Vector &q, const Vector &J, Vector &expJ, double B, i
 }\r
 \r
 \r
+// calculate energy due to fields and couplings\r
+double calculateEnergy(const Vector &J, const std::vector<int> &config, double &Efield, double &Ecoupling)\r
+{\r
+    Efield = 0;\r
+    Ecoupling = 0;\r
+    \r
+    for (unsigned i=0;i<config.size();i++) {\r
+\r
+        //assert(config[i] <= J[i].size());\r
+        \r
+        if (config[i]!=J[i].size()) {\r
+        \r
+            Efield -= J.at(i).at(config[i]);\r
+            \r
+            for (unsigned j=i+1;j<config.size();j++) {\r
+            \r
+                if (config[j]!=J[j].size())\r
+                    Ecoupling -= J.at(index(i,j,config.size())).at(sindex(config[i],config[j],J[i].size(),J[j].size()));\r
+                \r
+            }\r
+            \r
+        }\r
+        \r
+    }\r
+    \r
+    double E = Efield + Ecoupling;\r
+\r
+    return E;\r
+}\r
+\r
+\r
 // Compute correlations using rejection-free Monte Carlo (getError_RF) with epitope sites fixed\r
 \r
-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) {\r
+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) {\r
 \r
     // Create and initialize Monte Carlo variables\r
 \r
@@ -844,6 +888,8 @@ void getEpitopeCorrelations(const Vector &J, const Vector &expJ, const std::vect
     // MC variables\r
 \r
     int stepSize = 4*N;\r
+    //if (trackEnergy)\r
+    //    stepSize = N;\r
     \r
     // Set variables for parallel\r
     \r
@@ -855,6 +901,13 @@ void getEpitopeCorrelations(const Vector &J, const Vector &expJ, const std::vect
 \r
     MT::MersenneTwist mt;\r
     mt.init_genrand(seed[0]);\r
+\r
+    // possibly track energy\r
+    FILE *EOut = NULL;\r
+    if (trackEnergy) {\r
+        EOut = fopen(energyOutFileName.c_str(),"w");\r
+        if (EOut == NULL) { perror("File error"); exit(1); }\r
+    }\r
         \r
     // Thermalize and get initial error\r
         \r
@@ -883,7 +936,7 @@ void getEpitopeCorrelations(const Vector &J, const Vector &expJ, const std::vect
         } }\r
         \r
         setEnergies(expJ, neighbor, lattice, nonzero, expH, sumExpH);\r
-        getFixedSample(expJ, neighbor, fixedSites, mt, N, numSteps, stepSize, lattice, nonzero, p, expH, sumExpH);\r
+        getFixedSample(expJ, neighbor, fixedSites, mt, N, numSteps, stepSize, lattice, nonzero, p, expH, sumExpH, J, trackEnergy, EOut);\r
             \r
     }\r
         \r
@@ -891,6 +944,9 @@ void getEpitopeCorrelations(const Vector &J, const Vector &expJ, const std::vect
     \r
     for (int i=0;i<p.size();i++) { for (int j=0;j<p[i].size();j++) p[i][j] /= (double) (numThreads * numSteps * numRuns); }\r
 \r
+    if (trackEnergy)\r
+        fclose(EOut);\r
+\r
 }\r
 \r
 \r

diff --git a/monteCarlo.h b/monteCarlo.h
index 2eb5c8d..1bfce0e 100755 (executable)
--- 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 &);\r
 void updateCorrelationsP3(const std::vector<int> &, const std::vector<int> &, Vector &, std::vector<std::vector<std::vector<std::vector<double> > > > &);\r
 void getSample(const Vector &, const IntVector &, MT::MersenneTwist &, int, int, int, std::vector<int> &, std::vector<int> &, Vector &, Vector &, std::vector<double> &);\r
-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> &);\r
+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 *);\r
 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> &);\r
 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> &);\r
 \r
@@ -33,9 +33,10 @@ void getNeighbors(const Vector &, int, double, IntVector &);
 // Control functions\r
 void getError(const Vector &, const Vector &, int, double, int, int, double, double, std::vector<double> &);\r
 void getErrorMCLearn(const Vector &, const Vector &, Vector &, double, int, int, double, double, Vector &, std::vector<double> &, std::vector<int> &);\r
-void getEpitopeCorrelations(const Vector &, const Vector &, const std::vector<std::vector<int> > &, const std::vector<int> &, int, int, Vector &, std::vector<int> &);\r
+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 &);\r
 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> &);\r
 void getEpitopeStates(FILE *, const Vector &, const Vector &, const std::vector<std::vector<int> > &, const std::vector<int> &, int, int, Vector &, std::vector<int> &);\r
 \r
+double calculateEnergy(const Vector &, const std::vector<int> &, double &, double &);\r
 \r
 #endif\r

diff --git a/qEpitopeEval.cpp b/qEpitopeEval.cpp
index c1472dd..4807def 100755 (executable)
--- 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\r
     \r
     if (r.saveP3) getEpitopeCorrelationsP3Pk(J, expJ, eStates, eStartSite, r.b, r.runs, p, p3, pk, lattice);\r
-    else          getEpitopeCorrelations(J, expJ, eStates, eStartSite, r.b, r.runs, p, lattice);\r
+    else          getEpitopeCorrelations(J, expJ, eStates, eStartSite, r.b, r.runs, p, lattice, r.trackE, r.getEnergyOutFile());\r
     \r
     // Write out correlations\r
     \r

diff --git a/qEpitopeEval.h b/qEpitopeEval.h
index fcc4c67..a4d3e24 100755 (executable)
--- a/qEpitopeEval.h
+++ b/qEpitopeEval.h
@@ -34,6 +34,8 @@ public:
     bool saveP2;                    // If true, record pair correlations in additions to one-point correlations\r
     bool saveP3;                    // If true, record three-point correlations\r
     bool savePk;                    // If true, record P(k) mutations curve\r
+\r
+    bool trackE;                    // if true, save relative contribution to energy from fields and couplings\r
     \r
     RunParameters() {\r
         \r
@@ -49,6 +51,8 @@ public:
         saveP2=false;\r
         saveP3=false;\r
         savePk=false;\r
+\r
+        trackE = false;\r
         \r
     }\r
     std::string getInfile()                   { return (directory+"/"+infile+".j");                   }\r
@@ -59,6 +63,7 @@ public:
     std::string getPkOutfile()                { return (directory+"/"+outfile+".pk");                 }\r
     std::string getPkOutfile(int i)           { return (directory+"/"+outfile+"-"+tostring(i)+".pk"); }\r
     std::string getStartInfile()              { return (directory+"/"+startfile+".dat");              }\r
+    std::string getEnergyOutFile()            { return (directory+"/"+outfile+".E");                 }\r
     ~RunParameters() {}\r
     \r
 };\r