From: Dariusz Murakowski <murakdar@mit.edu>
Date: Mon, 6 Jul 2015 16:20:15 +0000 (-0400)
Subject: Code for doing Potts Monte Carlo (optionally with fixed epitope), from John Barton... 
X-Git-Url: http://src.murakowski.org/?a=commitdiff_plain;h=5bf53deb3213dd88f38d11f8c7521770b7cd791e;p=VirEvoDyn.git

Code for doing Potts Monte Carlo (optionally with fixed epitope), from John Barton 20150630.
---

diff --git a/MC code.zip b/MC code.zip
new file mode 100644
index 0000000..bfb9668
Binary files /dev/null and b/MC code.zip differ
diff --git a/MC_code.txt b/MC_code.txt
new file mode 100644
index 0000000..b75259f
--- /dev/null
+++ b/MC_code.txt
@@ -0,0 +1,21 @@
+from:   John P Barton <jpbarton@mit.edu>
+to: Dariusz Krzysztof Murakowski <murakdar@mit.edu>
+date:   Tue, Jun 30, 2015 at 6:21 PM
+subject:    MC Code
+
+Hi Dariusz,
+
+Here’s the MC code. The names are a little bit weird since I’ve been using it to evaluate epitope (or “epitope”) strengths. Example usage:
+
+g++ mainQEE.cpp io.cpp monteCarlo.cpp qEpitopeEval.cpp tools.cpp -O2 -march=native -o qee.out
+./qee.out -i (input couplings) -o (place to send the output) -mcb (number of MC samples to take)
+
+mainQEE is basically just used to read in the command line options. qEpitopeEval contains the main routine, and as you might have guessed monteCarlo holds all the MC routines. If you need any help with running it, let me know!
+
+Thanks,
+John
+
+
+Attachments:    
+MC code.zip (26 KB)
+
diff --git a/dataStructures.h b/dataStructures.h
new file mode 100755
index 0000000..3939481
--- /dev/null
+++ b/dataStructures.h
@@ -0,0 +1,87 @@
+#ifndef DATASTRUCTURES_H
+#define DATASTRUCTURES_H
+
+
+#include <vector>
+#include <math.h>
+#include <string>
+#include <sstream>
+#include <cstring>
+#include <iostream>
+#include <bitset>
+#include <stdio.h>
+
+#include "tools.h"  // Numerical tools
+
+
+
+// PROGRAM SETTINGS
+
+// This class holds the parameters needed for running the adaptive cluster algorithm
+
+class RunParameters {
+    
+public:
+    
+    std::string directory;  // Path to the directory where the inut file is located
+                            // Output will also be sent to this directory
+    std::string infile;     // Input file (prefix)
+    std::string ssinfile;   // Input file for secondary structure (default selected clusters)
+    std::string outfile;    // Output file (prefix)
+    
+    int kmin;               // Minimum cluster size before computation is truncated
+    int kmax;               // Maximum cluster size before computation is truncated
+    double sampleB;         // Number of samples used to compute correlations from data
+    double gamma0;          // Reference entropy regularization strength
+    double gamma2;          // Cluster entropy regularization strength
+    double thetaMin;        // The minimum cutoff value (loop starts here)
+    double thetaMax;        // The maximum cutoff value (loop ends here)
+    double thetaStep;       // Size of logarithmic steps when looping over theta
+    int b;                  // Number of data points to take in a single Monte Carlo run
+    int runs;               // Number of times to run Monte Carlo dynamics
+    bool computeGamma;      // If true, compute regularization strength based on number of samples
+    bool useLax;            // If true, use a laxer cluster construction rule
+    bool useRef;            // If false, do not use the reference entropy S0 (i.e. set S0, J0 to zero)
+    bool useSparse;         // If true, use sparse L0-norm regularization, in addition to the default L2
+    bool useSecStruct;      // If true, select first pairs from the secondary structure list if it exists
+    bool useVerbose;        // If true, print extra information while program is running
+    bool useVeryVerbose;    // If true, use very verbose output
+
+    
+    RunParameters() {
+        
+        directory=".";
+        infile="input";
+        ssinfile="input";
+        outfile="output";
+        
+        kmin=0;
+        kmax=0;
+        sampleB=10000;
+        gamma0=1.0e-4;
+        gamma2=0.0;
+        thetaMin=1.0e-10;
+        thetaMax=1.0e+0;
+        thetaStep=1.05;
+        b=40000;
+        runs=1;
+        computeGamma=false;
+        useLax=false;
+        useRef=false;
+        useSparse=false;
+        useSecStruct=false;
+        useVerbose=false;
+        useVeryVerbose=false;
+        
+    }
+    std::string getCorrelationsInfile()   { return (directory+"/"+infile+".p");      }
+    std::string getSecStructInfile()      { return (directory+"/"+ssinfile+".wuss"); }
+    std::string getCouplingsOutfile()     { return (directory+"/"+outfile+".j");     }
+    std::string getCorrelationsOutfile()  { return (directory+"/"+outfile+".p");     }  //#FLAG
+    std::string getSupplementaryOutfile() { return (directory+"/"+outfile+".sce");   }
+    ~RunParameters() {}
+    
+};
+
+
+#endif
diff --git a/io.cpp b/io.cpp
new file mode 100755
index 0000000..53aad26
--- /dev/null
+++ b/io.cpp
@@ -0,0 +1,97 @@
+#include <vector>
+
+#include "io.h"     // Data input and output
+#include "tools.h"  // Numerical tools
+
+
+// Read correlations and number of sites in from a file
+
+void getCorrelations(FILE *input, Vector &p) {
+
+    double c;
+    char o;
+    
+    while (fscanf(input,"%le",&c)==1) {
+    
+        p.push_back(std::vector<double>());
+        (p.back()).push_back(c);
+        
+        while (fscanf(input,"%c",&o)==1) {
+    
+            if (o=='\n' || o=='\r') break;
+            
+            fscanf(input,"%le",&c);
+            (p.back()).push_back(c);
+            
+        }
+        
+    }
+	
+}
+
+
+// Read couplings in from a file
+
+void getCouplings(FILE *input, Vector &p) {
+
+    getCorrelations(input, p);
+	
+}
+
+
+// Read secondary structure from file 
+
+void getSS(FILE *input, IntVector &p) {
+
+    int c;
+    char o;
+    
+    while (fscanf(input,"%d",&c)==1) {
+    
+        p.push_back(std::vector<int>());
+        (p.back()).push_back(c);
+        
+        while (fscanf(input,"%c",&o)==1) {
+    
+            if (o=='\n' || o=='\r') break;
+            
+            fscanf(input,"%d",&c);
+            (p.back()).push_back(c);
+            
+        }
+    
+    }
+
+}
+
+
+// Print final couplings out to file in a form understood by ising.cpp
+
+void printCouplings(FILE *output, const Vector &J) {
+    
+    for (int i=0;i<J.size();i++) {
+    
+        fprintf(output,"%.6e",J[i][0]);
+        
+        for (int j=1;j<J[i].size();j++) fprintf(output,"\t%.6e",J[i][j]);
+        
+        fprintf(output,"\n");
+        
+    }
+    
+    fflush(output);
+	
+}
+
+
+// Print supplementary information to file
+
+void printSupplementaryOutput(FILE *output, double theta, const std::vector<double> &error, double finalS, int maxClusterSize, unsigned long numClusters, unsigned long numSignificantClusters) {
+    
+    fprintf(output,"%le",theta);
+    for (int i=0;i<error.size();i++) fprintf(output,"\t%le",error[i]);
+    fprintf(output,"\t%le\t%d\t%lu\t%lu\n",finalS,maxClusterSize,numClusters,numSignificantClusters);
+    fflush(output);
+	
+}
+
diff --git a/io.h b/io.h
new file mode 100755
index 0000000..5c4bb05
--- /dev/null
+++ b/io.h
@@ -0,0 +1,19 @@
+#ifndef IO_H
+#define IO_H
+
+
+#include <vector>
+#include <math.h>
+#include <stdio.h>
+
+#include "tools.h"  // Numerical tools
+
+
+void getCorrelations(FILE *, Vector &);
+void getCouplings(FILE *, Vector &);
+void getSS(FILE *, IntVector &p);
+void printCouplings(FILE *, const Vector &);
+void printSupplementaryOutput(FILE *, double, const std::vector<double> &, double, int, unsigned long, unsigned long);
+
+
+#endif
diff --git a/mainQEE.cpp b/mainQEE.cpp
new file mode 100755
index 0000000..6b67dbe
--- /dev/null
+++ b/mainQEE.cpp
@@ -0,0 +1,169 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string>
+#include <math.h>
+#include <cstring>
+#include <iostream>
+#include <sstream>
+#include <cmath>
+#include <vector>
+#include <map>
+#include <algorithm>
+
+#include "qEpitopeEval.h"   // Evaluate epitopes
+#include "monteCarlo.h"     // Monte Carlo
+#include "tools.h"          // Numerical tools
+#include "io.h"             // Data I/O
+
+
+/*********************************************************************
+ 
+                   COMMAND LINE INPUT FORMAT
+ 
+ Command line instructions tell the program where to look for input
+ files and where to send output, as well as the setting of various
+ parameters (gamma, theta, etc) and flags (useSparse, etc).
+ Note that numerical parameters may be entered in either scientific
+ (recommended) or standard decimal notation. True/false switches are
+ off by default - if entered in the command line, the corresponding
+ option is set to true. Defaults can be found in dataStructures.h.
+ The conventions are given below:
+ 
+ -(flag name): (type of input expected)
+ 
+ -d: string
+    Default: "." (current directory)
+    Path to the directory where the data file is located, and where
+    output will be written.
+ 
+ -i: string
+    Default: "input"
+    The location of the file containing a set of couplings, the
+    starting values for the Monte Carlo learning algorithm.
+ 
+ -o: string
+    Default: "output"
+    The location of the file where output is to be sent. Each
+    different type of output file will have a different file type, 
+    e.g. .j for couplings.
+    
+ -c: string
+    Default: "input"
+    The location of the file containing the set of correlations to
+    reproduce (i.e. the correlations obtained from the data).
+    
+ -s: string
+    Default: none
+    Starting configuration for MC simulations.
+    
+ -g2: real number
+    Default: 0.0
+    The L2 regularization strength. L2 regularization is enabled by setting the 
+    regularization strength to a nonzero value using this flag, or by using the
+    -ag flag below.
+ 
+ -ag: none
+    Attempt to set the L2 regularization strengths to its optimal value,
+    based on the number of samples (input) in the data.
+    
+ -b: real number
+    Default: 1.0e+4
+    Number of samples used to compute the correlations. Used to 
+    determine the inference error.
+    
+ -mcb: integer
+    Default: 8.0e+5
+    Number of Monte Carlo samples to take to check inference error.
+    
+ -mcr: integer
+    Default: 1
+    Number of independent Monte Carlo runs to perform.
+    
+ -e: real number
+    Default: 1.0
+    Maximum tolerable error threshold. The MC learning algorithm will 
+    continue to run until the error on the one- and two-point correlations 
+    falls below this level.
+    
+ -v: none
+    Enable verbose output.
+ 
+ *********************************************************************/
+
+
+// MAIN PROGRAM
+
+int main(int argc, char *argv[]) {
+    
+	RunParameters r;
+    
+    // Process command line input
+    
+    for (int i=1;i<argc;i++) {
+        
+        // Location of input/output files
+        
+        if      (strcmp(argv[i],"-d")==0)   { if (++i==argc) break; else r.directory=argv[i];                               }
+        else if (strcmp(argv[i],"-i")==0)   { if (++i==argc) break; else r.infile=argv[i];                                  }
+        else if (strcmp(argv[i],"-o")==0)   { if (++i==argc) break; else r.outfile=argv[i];                                 }
+        else if (strcmp(argv[i],"-s")==0)   { if (++i==argc) break; else { r.useStart=true; r.startfile=argv[i]; }          }
+        else if (strcmp(argv[i],"-on")==0)  { if (++i==argc) break; else r.on=strtoint(argv[i]);                            }
+        
+        // Epitope input
+        
+        else if (strcmp(argv[i],"-n")==0)   { r.epitopeStates.push_back(std::vector<std::vector<int> >());
+                                              r.epitopeStartSite.push_back(std::vector<int>());                             }
+        else if (strcmp(argv[i],"-e")==0)   { if (++i==argc) break;
+                                              else if (strcmp(argv[i],"[")==0) {
+                                                std::vector<int> tempvec;
+                                                while (strcmp(argv[++i],"]")!=0) tempvec.push_back(strtoint(argv[i]));
+                                                (r.epitopeStates.back()).push_back(tempvec);
+                                              } }
+        else if (strcmp(argv[i],"-es")==0)  { if (++i==argc) break; else (r.epitopeStartSite.back()).push_back(strtoint(argv[i])); }
+        
+        // Monte Carlo settings
+        
+        else if (strcmp(argv[i],"-mcb")==0) { if (++i==argc) break; else r.b=strtoint(argv[i]);                             }
+        else if (strcmp(argv[i],"-mcr")==0) { if (++i==argc) break; else r.runs=strtoint(argv[i]);                          }
+        
+        // Optional output
+        
+        else if (strcmp(argv[i],"-p2")==0)  { r.saveP2=true;                                                                }
+        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 printf("Unrecognized command! '%s'\n",argv[i]);
+        
+    }
+    
+    if (r.epitopeStates.size()==0) {
+        
+        std::vector<std::vector<int> > eStates;
+        std::vector<int> eStartSite;
+        
+        runEpitopeEval(r, eStates, eStartSite, -1);
+        
+    }
+    
+    else if (r.epitopeStates.size()==1) {
+    
+        std::vector<std::vector<int> > eStates(r.epitopeStates[0]);
+        std::vector<int> eStartSite(r.epitopeStartSite[0]);
+        
+        runEpitopeEval(r, eStates, eStartSite, -1);
+    
+    }
+    
+    else { for (int i=0;i<r.epitopeStates.size();i++) {
+    
+        std::vector<std::vector<int> > eStates(r.epitopeStates[i]);
+        std::vector<int> eStartSite(r.epitopeStartSite[i]);
+    
+        runEpitopeEval(r, eStates, eStartSite, r.on+i);
+        
+    } }
+    
+	return 0;
+
+}
diff --git a/monteCarlo.cpp b/monteCarlo.cpp
new file mode 100644
index 0000000..9b34b80
--- /dev/null
+++ b/monteCarlo.cpp
@@ -0,0 +1,1076 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string>
+#include <math.h>
+#include <cstring>
+#include <iostream>
+#include <sstream>
+#include <cmath>
+#include <vector>
+#include <map>
+#include <algorithm>
+
+#include "tools.h"      // Numerical tools
+#include "monteCarlo.h" // MC declarations
+#include "mtrnd.h"      // Random number generator
+
+
+// Compute epsilonP for two vectors of correlations
+
+double epsilonP(const Vector &q, const Vector &p, int N, double maxPrecision, const Vector &J, double gamma, double alpha) {
+    
+    double err  = 0.0;
+    double Neff = 0.0;
+    
+    for (int i=0;i<N;i++) { for (int a=0;a<q[i].size();a++) {
+    
+        Neff++;
+        double gradError = pow(q[i][a] - p[i][a] - (2 * gamma * alpha * J[i][a]) , 2.0);
+    
+        if (q[i][a]<maxPrecision) err += gradError / (maxPrecision * (1 - maxPrecision));
+        else                      err += gradError / (q[i][a]      * (1 - q[i][a]));
+        
+    }   }
+    
+    return sqrt( err / (Neff * maxPrecision) );
+    
+}
+
+
+// Compute epsilonP2 for two vectors of correlations
+
+double epsilonP2(const Vector &q, const Vector &p, int N, double maxPrecision, const Vector &J, double gamma) {
+    
+    double err   = 0;
+    double NJeff = 0;
+    
+    for (int i=0;i<N;i++) { for (int j=i+1;j<N;j++) {
+    
+        int idx = index(i,j,N);
+        
+        //Get quantities for computing K_ij^ab
+    
+        double inv_qi = 1/((double)J[i].size()+1);
+        double inv_qj = 1/((double)J[j].size()+1);
+    
+        double suma[J[i].size()];
+        double sumb[J[j].size()];
+        double sumall=0;
+        
+        for (int a=0;a<J[i].size();a++) suma[a]=0;
+        for (int b=0;b<J[j].size();b++) sumb[b]=0;
+        
+        for (int a=0;a<J[i].size();a++) { for (int b=0;b<J[j].size();b++) {
+            
+            suma[a] += J[idx][sindex(a,b,J[i].size(),J[j].size())];
+            sumb[b] += J[idx][sindex(a,b,J[i].size(),J[j].size())];
+            
+        } }
+        
+        for (int a=0;a<J[i].size();a++) sumall += suma[a];
+        
+        for (int a=0;a<J[i].size();a++) { for (int b=0;b<J[j].size();b++) {
+        
+            NJeff++;
+
+            int    sab       = sindex(a,b,J[i].size(),J[j].size());
+            //double K         = J[idx][sab];
+            double K         = J[idx][sab] - (suma[a] * inv_qj) - (sumb[b] * inv_qi) + (sumall * inv_qi * inv_qj);
+            double gradError = pow(p[idx][sab] - q[idx][sab] + (2 * gamma * K), 2.0);
+    
+            if (q[idx][sab]<maxPrecision) err += gradError/(maxPrecision * (1 - maxPrecision));
+            else                          err += gradError/( q[idx][sab] * (1 - q[idx][sab]));
+
+        } }
+            
+    } }
+    
+    return sqrt( err / (NJeff * maxPrecision) );
+    
+}
+
+
+// Compute epsilonC for two vectors of correlations
+
+double epsilonC(const Vector &q, const Vector &p, int N, double maxPrecision, const Vector &J, double gamma) {
+    
+    double err   = 0;
+    double NJeff = 0;
+    
+    for (int i=0;i<N;i++) { for (int j=i+1;j<N;j++) {
+    
+        int idx = index(i,j,N);
+        
+        // Get quantities for computing K_ij^ab
+    
+        double inv_qi = 1/((double)J[i].size()+1);
+        double inv_qj = 1/((double)J[j].size()+1);
+    
+        double suma[J[i].size()];
+        double sumb[J[j].size()];
+        double sumall=0;
+        
+        for (int a=0;a<J[i].size();a++) suma[a]=0;
+        for (int b=0;b<J[j].size();b++) sumb[b]=0;
+        
+        for (int a=0;a<J[i].size();a++) { for (int b=0;b<J[j].size();b++) {
+            
+            suma[a] += J[idx][sindex(a,b,J[i].size(),J[j].size())];
+            sumb[b] += J[idx][sindex(a,b,J[i].size(),J[j].size())];
+            
+        } }
+        
+        for (int a=0;a<J[i].size();a++) sumall += suma[a];
+        
+        for (int a=0;a<J[i].size();a++) { for (int b=0;b<J[j].size();b++) {
+        
+            NJeff++;
+
+            int    sab       = sindex(a,b,J[i].size(),J[j].size());
+            //double K         = J[idx][sab];
+            double K         = J[idx][sab] - (suma[a] * inv_qj) - (sumb[b] * inv_qi) + (sumall * inv_qi * inv_qj);
+            double gradError = pow((q[idx][sab]-q[i][a]*q[j][b]) - (p[idx][sab]-p[i][a]*p[j][b]) - (2 * gamma * K) , 2.0);
+            
+            if (q[idx][sab]<maxPrecision) err += gradError / pow(sqrt(maxPrecision * (1 - maxPrecision)) + q[i][a] * sqrt(q[j][b] * (1 - q[j][b])) + q[j][b] * sqrt(q[i][a] * (1 - q[i][a])),2);
+            else                          err += gradError / pow(sqrt(q[idx][sab]  * (1 -  q[idx][sab])) + q[i][a] * sqrt(q[j][b] * (1 - q[j][b])) + q[j][b] * sqrt(q[i][a] * (1 - q[i][a])),2);
+
+        } }
+            
+    } }
+    
+    return sqrt( err / (NJeff * maxPrecision) );
+    
+}
+
+
+// Compute maximum error on the one- and two-point correlations (sparse)
+
+double getMaxError(const Vector &q, const Vector &p, const IntVector &nz, double maxPrecision) {
+
+    int N = sizetolength(q.size());
+    double maxError = 0;
+    
+    // Check all single sites
+        
+    for (int i=0;i<N;i++) { for (int j=0;j<q[i].size();j++) {
+    
+        double gradError = pow(q[i][j] - p[i][j], 2.0);
+        
+        if (q[i][j]<maxPrecision) gradError /= (1 - maxPrecision);
+        else                      gradError /= (maxPrecision * q[i][j] * (1 - q[i][j]));
+            
+        if (gradError>maxError) maxError=gradError;
+        
+    } }
+        
+    // Check nonzero pairs
+        
+    for (int n=0;n<nz.size();n++) {
+    
+        int i=nz[n][0];
+        int j=nz[n][1];
+            
+        int idx=index(i,j,N);
+            
+        for (int k=0;k<q[idx].size();k++) {
+            
+            double gradError = pow(q[idx][k] - p[idx][k], 2.0);
+            
+            if (q[idx][k]<maxPrecision) gradError /= (1 - maxPrecision);
+            else                        gradError /= (maxPrecision * q[idx][k] * (1 - q[idx][k]));
+
+            if (gradError>maxError) maxError=gradError;
+            
+        }
+            
+    }
+    
+    return ( maxError / sqrt(2 * log(N*(N+1)/2) ) );
+
+}
+
+
+// Compute maximum error on the one- and two-point correlations (dense, including regularization)
+
+double getMaxError(const Vector &q, const Vector &p, double maxPrecision, const Vector &J, double gamma, double alpha) {
+
+    int    N        = sizetolength(q.size());
+    double maxError = 0;
+    
+    // Check all single sites
+        
+    for (int i=0;i<N;i++) { for (int j=0;j<q[i].size();j++) {
+    
+        double gradError = pow(q[i][j] - p[i][j] - (2 * gamma * alpha * J[i][j]) , 2.0);
+        
+        if (q[i][j]<maxPrecision) gradError /= (1 - maxPrecision);
+        else                      gradError /= (maxPrecision * q[i][j] * (1 - q[i][j]));
+            
+        if (gradError>maxError) maxError = gradError;
+        
+    } }
+        
+    // Check all pairs
+        
+    for (int i=0;i<N;i++) { for (int j=i+1;j<N;j++) {
+    
+        int idx = index(i,j,N);
+        
+        // Get quantities for computing K_ij^ab
+    
+        double inv_qi = 1/((double)J[i].size()+1);
+        double inv_qj = 1/((double)J[j].size()+1);
+    
+        double suma[J[i].size()];
+        double sumb[J[j].size()];
+        double sumall=0;
+        
+        for (int a=0;a<J[i].size();a++) suma[a]=0;
+        for (int b=0;b<J[j].size();b++) sumb[b]=0;
+        
+        for (int a=0;a<J[i].size();a++) { for (int b=0;b<J[j].size();b++) {
+            
+            suma[a] += J[idx][sindex(a,b,J[i].size(),J[j].size())];
+            sumb[b] += J[idx][sindex(a,b,J[i].size(),J[j].size())];
+            
+        } }
+        
+        for (int a=0;a<J[i].size();a++) sumall += suma[a];
+        
+        for (int a=0;a<J[i].size();a++) { for (int b=0;b<J[j].size();b++) {
+
+            int    sab       = sindex(a,b,J[i].size(),J[j].size());
+            //double K         = J[idx][sab];
+            double K         = J[idx][sab] - (suma[a] * inv_qj) - (sumb[b] * inv_qi) + (sumall * inv_qi * inv_qj);
+            double gradError = pow(p[idx][sab] - q[idx][sab] + (2 * gamma * K), 2.0);
+    
+            if (q[idx][sab]<maxPrecision) gradError /= (1 - maxPrecision);
+            else                          gradError /= (maxPrecision * q[idx][sab] * (1 - q[idx][sab]));
+
+            if (gradError>maxError) maxError=gradError;
+            
+        } }
+            
+    } }
+    
+    return ( maxError / sqrt(2 * log(N*(N+1)/2) ) );
+
+}
+
+
+// Set energies based on an initial configuration
+
+void setEnergies(const Vector &expJ, const IntVector &neighbor, const std::vector<int> &lattice, const std::vector<int> &nonzero, Vector &expH, std::vector<double> &sumExpH) {
+    
+    // Get energy of each state at each site (modulo fields/couplings at other sites)
+    
+    for (int i=0;i<sumExpH.size();i++) { sumExpH[i] = 1; for (int j=0;j<expH[i].size();j++) expH[i][j]=1; }
+
+    for (int i=0;i<expH.size();i++) { for (int j=0;j<expH[i].size()-1;j++) {
+        
+        expH[i][j] = expJ[i][j];
+            
+        for (int k=1;k<neighbor[i].size();k++) { if (nonzero[neighbor[i][k]]) {
+        
+            if (i<neighbor[i][k]) expH[i][j] *= expJ[index(i,neighbor[i][k],lattice.size())][sindex(j,lattice[neighbor[i][k]],expJ[i].size(),expJ[neighbor[i][k]].size())];
+            else                  expH[i][j] *= expJ[index(neighbor[i][k],i,lattice.size())][sindex(lattice[neighbor[i][k]],j,expJ[neighbor[i][k]].size(),expJ[i].size())];
+                
+        } }
+            
+        sumExpH[i] += expH[i][j];
+            
+    } }
+    
+}
+
+
+// Update energies after a successful spin flip
+
+void updateEnergies(const Vector &expJ, const IntVector &neighbor, const std::vector<int> &lattice, const std::vector<int> &nonzero, int x, int prev, Vector &expH, std::vector<double> &sumExpH) {
+
+    // Initial state zero
+    
+    if (prev==expJ[x].size()) {
+    
+        for (int i=1;i<neighbor[x].size();i++) {
+            
+            int n = neighbor[x][i];
+            
+            sumExpH[n] = 1;
+            
+            if (x<n) { for (int j=0;j<expH[n].size()-1;j++) {
+            
+                expH[n][j] *= expJ[index(x,n,lattice.size())][sindex(lattice[x],j,expJ[x].size(),expJ[n].size())];
+                sumExpH[n] += expH[n][j];
+                
+            } }
+                
+            else { for (int j=0;j<expH[n].size()-1;j++) {
+            
+                expH[n][j] *= expJ[index(n,x,lattice.size())][sindex(j,lattice[x],expJ[n].size(),expJ[x].size())];
+                sumExpH[n] += expH[n][j];
+                
+            } }
+            
+        }
+        
+    }
+    
+    // Final state zero
+    
+    else if (lattice[x]==expJ[x].size()) {
+    
+        for (int i=1;i<neighbor[x].size();i++) {
+            
+            int n = neighbor[x][i];
+            
+            sumExpH[n] = 1;
+            
+            if (x<n) { for (int j=0;j<expH[n].size()-1;j++) {
+                
+                expH[n][j] /= expJ[index(x,n,lattice.size())][sindex(prev,j,expJ[x].size(),expJ[n].size())];
+                sumExpH[n] += expH[n][j];
+                    
+            } }
+            
+            else { for (int j=0;j<expH[n].size()-1;j++) {
+                
+                expH[n][j] /= expJ[index(n,x,lattice.size())][sindex(j,prev,expJ[n].size(),expJ[x].size())];
+                sumExpH[n] += expH[n][j];
+                
+            } }
+        
+        }
+        
+    }
+    
+    // Both initial and final state are nonzero
+    
+    else {
+    
+        for (int i=1;i<neighbor[x].size();i++) {
+            
+            int n = neighbor[x][i];
+            
+            sumExpH[n] = 1;
+            
+            if (x<n) { for (int j=0;j<expH[n].size()-1;j++) {
+                
+                expH[n][j] *= expJ[index(x,n,lattice.size())][sindex(lattice[x],j,expJ[x].size(),expJ[n].size())]/expJ[index(x,n,lattice.size())][sindex(prev,j,expJ[x].size(),expJ[n].size())];
+                sumExpH[n] += expH[n][j];
+                
+            } }
+            
+            else { for (int j=0;j<expH[n].size()-1;j++) {
+            
+                expH[n][j] *= expJ[index(n,x,lattice.size())][sindex(j,lattice[x],expJ[n].size(),expJ[x].size())]/expJ[index(n,x,lattice.size())][sindex(j,prev,expJ[n].size(),expJ[x].size())];
+                sumExpH[n] += expH[n][j];
+                
+            } }
+            
+        }
+        
+    }
+
+}
+
+
+// Rejection-free Monte Carlo step (update energies on flip only)
+
+void dynamicsRF(const Vector &expJ, const IntVector &neighbor, int x, double y, std::vector<int> &lattice, std::vector<int> &nonzero, Vector &expH, std::vector<double> &sumExpH) {
+    
+    // Check for spin flip, choose new state
+    
+    if ( y < ((sumExpH[x] - expH[x][lattice[x]]) / sumExpH[x]) ) {
+        
+        double tempSum=0;
+        int prev=lattice[x];
+        
+        for (int i=0;i<expH[x].size();i++) { if (i!=lattice[x]) {
+            
+            tempSum += expH[x][i]/sumExpH[x];
+            
+            if (y<tempSum) {
+                
+                if (lattice[x]==expJ[x].size()) nonzero[x]=1;
+                else if (i==expJ[x].size())     nonzero[x]=0;
+                
+                lattice[x]=i;
+                break;
+                
+            }
+            
+        } }
+        
+        updateEnergies(expJ, neighbor, lattice, nonzero, x, prev, expH, sumExpH);
+        
+    }
+
+}
+
+
+// Update correlations based on the current lattice configuration
+
+void updateCorrelations(const std::vector<int> &lattice, const std::vector<int> &nonzero, Vector &p) {
+
+    for (int i=0;i<lattice.size();i++) { if (nonzero[i]) {
+            
+        p[i][lattice[i]]+=1;
+            
+        int off=offset(i, lattice.size());
+            
+        for (int j=i+1;j<lattice.size();j++) { if (nonzero[j]) {
+        
+            p[off + j][sindex(lattice[i],lattice[j],p[i].size(),p[j].size())]+=1;
+                
+        } }
+            
+    } }
+
+}
+
+
+// Update correlations based on the current lattice configuration
+
+void updateCorrelationsP3(const std::vector<int> &lattice, const std::vector<int> &nonzero, Vector &p, std::vector<std::vector<std::vector<std::vector<double> > > > &p3) {
+
+    for (int i=0;i<lattice.size();i++) { if (nonzero[i]) {
+            
+        p[i][lattice[i]]+=1;
+            
+        int off=offset(i, lattice.size());
+            
+        for (int j=i+1;j<lattice.size();j++) { if (nonzero[j]) {
+        
+            p[off + j][sindex(lattice[i],lattice[j],p[i].size(),p[j].size())]+=1;
+            
+            for (int k=j+1;k<lattice.size();k++) { if (nonzero[k]) {
+            
+                p3[i][j-i][k-j][sindex3(lattice[i],lattice[j],lattice[k],p[i].size(),p[j].size(),p[k].size())]+=1;
+                
+            } }
+                
+        } }
+            
+    } }
+
+}
+
+
+// Get Monte Carlo sample
+
+void getSample(const Vector &expJ, const IntVector &neighbor, 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) {
+
+    for (int n=0;n<numSteps;n++) {
+    
+        //for (int k=0;k<stepSize;k++) dynamicsRF(expJ, neighbor, mt.genrand_int31()%N, mt.genrand_real1(), lattice, nonzero, expH, sumExpH);
+        for (int k=0;k<stepSize;k++) dynamicsRF(expJ, neighbor, k%N, mt.genrand_real1(), lattice, nonzero, expH, sumExpH);
+        updateCorrelations(lattice, nonzero, p);
+        
+    }
+
+}
+
+
+// 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) {
+
+    for (int n=0;n<numSteps;n++) {
+    
+        for (int k=0;k<stepSize;k++) {
+            
+            bool allowed = true;
+            for (int i=0;i<fixedSites.size();i++) { if (fixedSites[i]==k%N) allowed = false; }
+            
+            if (allowed) dynamicsRF(expJ, neighbor, k%N, mt.genrand_real1(), lattice, nonzero, expH, sumExpH);
+            
+        }
+        
+        updateCorrelations(lattice, nonzero, p);
+        
+        // Cut off if too many mutations
+        int mutCut = (N/2 < 80) ? N/2 : 120;
+        int nmut   = 0;
+        for (int i=0;i<N;i++) { if (nonzero[i]) nmut++; }
+        
+        if (nmut>mutCut) {
+        
+            for (int i=0;i<lattice.size();i++) {
+            
+                bool allowed = true;
+                for (int k=0;k<fixedSites.size();k++) { if (fixedSites[k]==i) allowed = false; }
+                
+                if (allowed) { lattice[i] = (int) p[i].size(); nonzero[i] = 0; }
+                
+            }
+                    
+            setEnergies(expJ, neighbor, lattice, nonzero, expH, sumExpH);
+        
+        }
+        // Cut off if too many mutations
+        
+    }
+
+}
+
+
+// Get Monte Carlo sample with some sites fixed
+
+void getFixedSampleP3Pk(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, std::vector<std::vector<std::vector<std::vector<double> > > > &p3, std::vector<double> &pk, Vector &expH, std::vector<double> &sumExpH) {
+
+    for (int n=0;n<numSteps;n++) {
+    
+        for (int k=0;k<stepSize;k++) {
+            
+            bool allowed = true;
+            for (int i=0;i<fixedSites.size();i++) { if (fixedSites[i]==k%N) allowed = false; }
+            
+            if (allowed) dynamicsRF(expJ, neighbor, k%N, mt.genrand_real1(), lattice, nonzero, expH, sumExpH);
+            
+        }
+        
+        updateCorrelationsP3(lattice, nonzero, p, p3);
+        
+        // Cut off if too many mutations
+        int mutCut = (N/2 < 80) ? N/2 : 120;
+        int nmut   = 0;
+        for (int i=0;i<N;i++) { if (nonzero[i]) nmut++; }
+        
+        pk[nmut] += 1;  // COUNT NUMBER OF MUTATIONS
+        
+        if (nmut>mutCut) {
+        
+            for (int i=0;i<lattice.size();i++) {
+            
+                bool allowed = true;
+                for (int k=0;k<fixedSites.size();k++) { if (fixedSites[k]==i) allowed = false; }
+                
+                if (allowed) { lattice[i] = (int) p[i].size(); nonzero[i] = 0; }
+                
+            }
+                    
+            setEnergies(expJ, neighbor, lattice, nonzero, expH, sumExpH);
+        
+        }
+        // Cut off if too many mutations
+        
+    }
+
+}
+
+
+// Write Monte Carlo sample to file with some sites fixed
+
+void writeFixedSample(FILE* output, 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) {
+
+    for (int n=0;n<numSteps;n++) {
+    
+        for (int k=0;k<stepSize;k++) {
+            
+            bool allowed = true;
+            for (int i=0;i<fixedSites.size();i++) { if (fixedSites[i]==k%N) allowed = false; }
+            
+            if (allowed) dynamicsRF(expJ, neighbor, k%N, mt.genrand_real1(), lattice, nonzero, expH, sumExpH);
+            
+        }
+        
+        for (int i=0;i<lattice.size();i++) fprintf(output,"%d ",lattice[i]);
+        //for (int i=0;i<lattice.size();i++) { if (nonzero[i]) fprintf(output,"%d %d\t",i,lattice[i]); }
+        fprintf(output,"\n");
+        
+    }
+
+}
+
+
+// Get effective number of fields and couplings
+
+void getNeff(const Vector &q, int N, int &Neff, int &NJeff) {
+
+    for (int i=0;i<N;i++) {
+    
+        Neff+=(int) q[i].size();
+        for (int j=i+1;j<N;j++) NJeff+=(int) (q[i].size()*q[j].size());
+        
+    }
+
+}
+
+
+// Get the set of interacting "neighbors" for each site
+
+void getNeighbors(const Vector &J, int N, double threshold, IntVector &neighbor) {
+
+    for (int i=0;i<N;i++) neighbor[i].push_back(i);
+    
+    for (int i=0;i<N;i++) {
+        
+        int off=offset(i, N);
+        
+        for (int j=i+1;j<N;j++) { for (int k=0;k<J[off + j].size();k++) {
+            
+            if (fabs(J[off + j][k])>threshold) {
+                
+                neighbor[i].push_back(j);
+                neighbor[j].push_back(i);
+                break;
+               
+            }
+                
+        } }
+        
+    }
+
+}
+
+
+// Compute the inferrence error for a given set of couplings using rejection-free Monte Carlo (getError_RF)
+
+void getError(const Vector &q, const Vector &J, int N, double B, int numSteps, int numRuns, double gamma, double alpha, std::vector<double> &error) {
+
+    // Create and initialize Monte Carlo variables
+    
+    std::vector<double> lastError(2,0);
+    Vector p(q.size(),std::vector<double>());
+    for (int i=0;i<p.size();i++) p[i].resize(q[i].size(),0);
+    
+    std::vector<int> lattice(N,0);
+    std::vector<int> nonzero(N,0);
+    
+    // Couplings and neighbors
+    
+    IntVector neighbor(N);
+    getNeighbors(J, N, 0, neighbor);
+    Vector expJ(J);
+    for (int i=0;i<J.size();i++) { for (int j=0;j<J[i].size();j++) expJ[i][j]=exp(J[i][j]); }
+    
+    // Energies
+    
+    Vector expH;
+    std::vector<double> sumExpH;
+    
+    sumExpH.resize(lattice.size(),0);
+    expH.resize(lattice.size(),std::vector<double>());
+    for (int i=0;i<expH.size();i++) expH[i].resize(expJ[i].size()+1,1);
+    
+    // MC variables
+    
+    double multiplier   = 1;
+    double multMax      = 64;
+    int    stepSize     = 4*N;
+    double maxPrecision = 1/B;
+    bool   runMC        = true;
+    
+    // Set variables for parallel
+    
+    int numThreads=1;
+    
+    srand((unsigned)time(0));
+    unsigned long seed[numThreads];
+    for (int i=0;i<numThreads;i++) seed[i]=rand();
+    
+    for (int i=0;i<lattice.size();i++) { lattice[i]=(int) q[i].size(); nonzero[i]=0; }
+
+    MT::MersenneTwist mt;
+    mt.init_genrand(seed[0]);
+    
+    // Thermalize and get initial error
+    
+    setEnergies(expJ, neighbor, lattice, nonzero, expH, sumExpH);
+    getSample(expJ, neighbor, mt, N, numSteps, stepSize, lattice, nonzero, p, expH, sumExpH);
+    
+    // Normalize correlations and compute initial error
+
+    for (int i=0;i<q.size();i++) { for (int j=0;j<q[i].size();j++) p[i][j] /= (double) (numThreads * numSteps * multiplier); }
+    
+    error[0] = epsilonP(q, p, N, maxPrecision, J, gamma, alpha);
+    error[1] = epsilonC(q, p, N, maxPrecision, J, gamma);
+    error[2] = getMaxError(q, p, maxPrecision, J, gamma, alpha);
+    
+    for (int i=0;i<q.size();i++) { for (int j=0;j<q[i].size();j++) p[i][j] = 0; }
+    
+    // Go to MC
+    
+    while (runMC) {
+                            
+        for (int i=0;i<lattice.size();i++) { lattice[i]=(int) q[i].size(); nonzero[i] = 0; }
+    
+        setEnergies(expJ, neighbor, lattice, nonzero, expH, sumExpH);
+        getSample(expJ, neighbor, mt, N, numSteps*multiplier, stepSize, lattice, nonzero, p, expH, sumExpH);
+    
+        // Normalize correlations and compute error
+
+        for (int i=0;i<q.size();i++) { for (int j=0;j<q[i].size();j++) p[i][j] /= (double) (numThreads * numSteps * multiplier); }
+    
+        lastError[0] = error[0];
+        lastError[1] = error[1];
+    
+        error[0] = epsilonP(q, p, N, maxPrecision, J, gamma, alpha);
+        error[1] = epsilonC(q, p, N, maxPrecision, J, gamma);
+        error[2] = getMaxError(q, p, maxPrecision, J, gamma, alpha);
+    
+        // Check exit condition
+            
+        double deltaP = fabs(error[0] - lastError[0]);
+        double deltaC = fabs(error[1] - lastError[1]);
+  
+        bool pBad  = ( (error[0] > 1.0) && (error[0] - 3 * deltaP > 1) );
+        bool cBad  = ( (error[1] > 1.0) && (error[1] - 3 * deltaC > 1) );
+        bool pGood = ( (error[0] < 1.0) && (error[0] + 1 * deltaP < 1) );
+        bool cGood = ( (error[1] < 1.0) && (error[1] + 1 * deltaC < 1) );
+            
+        if ( pBad || cBad || (pGood && cGood) || multiplier==multMax ) runMC=false;
+        
+        else {
+            
+            for (int i=0;i<q.size();i++) { for (int j=0;j<q[i].size();j++) p[i][j] = 0; }
+            multiplier*=2;
+            //printf("m=%d, e_p=%.4e, e_c=%.4e\n",(int)multiplier,error[0],error[1]);
+                
+        }
+    
+    }
+
+}
+
+
+// Compute the inferrence error for a given set of couplings using rejection-free Monte Carlo (getError_RF)
+
+void getErrorMCLearn(const Vector &q, const Vector &J, Vector &expJ, double B, int numSteps, int numRuns, double gamma, double alpha, Vector &p, std::vector<double> &error, std::vector<int> &latticeStart) {
+
+    // Create and initialize Monte Carlo variables
+
+    for (int i=0;i<p.size();i++) { for (int j=0;j<p[i].size();j++) p[i][j]=0; }
+
+    int N = sizetolength(J.size());
+    
+    std::vector<int> lattice(N,0);
+    std::vector<int> nonzero(N,0);
+    
+    // Neighbors
+    
+    IntVector neighbor(N);
+    getNeighbors(J, N, 1.0e-5, neighbor);
+    
+    // Energies
+    
+    Vector expH;
+    std::vector<double> sumExpH;
+    
+    sumExpH.resize(lattice.size(),0);
+    expH.resize(lattice.size(),std::vector<double>());
+    for (int i=0;i<expH.size();i++) expH[i].resize(expJ[i].size()+1,1);
+    
+    // MC variables
+
+    int    stepSize     = 4*N;
+    int    mutCut       = (N/2 < 80) ? N/2 : 120;
+    double maxPrecision = 1/B;
+    
+    // Set variables for parallel
+    
+    int numThreads=1;
+    
+    srand((unsigned)time(0));
+    unsigned long seed[numThreads];
+    for (int i=0;i<numThreads;i++) seed[i]=rand();
+
+    MT::MersenneTwist mt;
+    mt.init_genrand(seed[0]);
+        
+    // Thermalize and get initial error
+        
+    for (int m=0;m<numRuns;m++) {
+        
+        for (int i=0;i<lattice.size();i++) {
+        
+            lattice[i]=latticeStart[i];
+            if (lattice[i]==q[i].size()) nonzero[i]=0;
+            else                         nonzero[i]=1;
+            
+        }
+        
+        setEnergies(expJ, neighbor, lattice, nonzero, expH, sumExpH);
+        getSample(expJ, neighbor, mt, N, numSteps, stepSize, lattice, nonzero, p, expH, sumExpH);
+            
+    }
+        
+    // Normalize correlations and compute error
+    
+    for (int i=0;i<q.size();i++) { for (int j=0;j<q[i].size();j++) p[i][j] /= (double) (numThreads * numSteps * numRuns); }
+        
+    error[0] = epsilonP( q, p, N, maxPrecision, J, gamma, alpha);
+    error[1] = epsilonP2(q, p, N, maxPrecision, J, gamma);
+    error[2] = getMaxError( q, p, maxPrecision, J, gamma, alpha);
+
+}
+
+
+// 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) {
+
+    // Create and initialize Monte Carlo variables
+
+    for (int i=0;i<p.size();i++) { for (int j=0;j<p[i].size();j++) p[i][j]=0; }
+
+    int N = sizetolength(J.size());
+    
+    std::vector<int> lattice(N,0);
+    std::vector<int> nonzero(N,0);
+    std::vector<int> fixedSites;
+    
+    for (int i=0;i<epitopeStartSite.size();i++) { for (int j=0;j<epitopeStates[i].size();j++) {
+    
+        fixedSites.push_back(epitopeStartSite[i]+j);
+        
+    } }
+    
+    // Neighbors
+    
+    IntVector neighbor(N);
+    getNeighbors(J, N, 1.0e-5, neighbor);
+    
+    // Energies
+    
+    Vector expH;
+    std::vector<double> sumExpH;
+    
+    sumExpH.resize(lattice.size(),0);
+    expH.resize(lattice.size(),std::vector<double>());
+    for (int i=0;i<expH.size();i++) expH[i].resize(expJ[i].size()+1,1);
+    
+    // MC variables
+
+    int stepSize = 4*N;
+    
+    // Set variables for parallel
+    
+    int numThreads=1;
+    
+    srand((unsigned)time(0));
+    unsigned long seed[numThreads];
+    for (int i=0;i<numThreads;i++) seed[i]=rand();
+
+    MT::MersenneTwist mt;
+    mt.init_genrand(seed[0]);
+        
+    // Thermalize and get initial error
+        
+    for (int m=0;m<numRuns;m++) {
+        
+        for (int i=0;i<lattice.size();i++) {
+        
+            lattice[i]=latticeStart[i];
+            
+            if (lattice[i]==J[i].size()) nonzero[i]=0;
+            else                         nonzero[i]=1;
+            
+        }
+        
+        // Fix epitope region
+        int k=0;
+        for (int i=0;i<epitopeStates.size();i++) { for (int j=0;j<epitopeStates[i].size();j++) {
+        
+            lattice[fixedSites[k]]=epitopeStates[i][j];
+            
+            if (lattice[fixedSites[k]]==J[fixedSites[k]].size()) nonzero[fixedSites[k]]=0;
+            else                                                 nonzero[fixedSites[k]]=1;
+            
+            k++;
+            
+        } }
+        
+        setEnergies(expJ, neighbor, lattice, nonzero, expH, sumExpH);
+        getFixedSample(expJ, neighbor, fixedSites, mt, N, numSteps, stepSize, lattice, nonzero, p, expH, sumExpH);
+            
+    }
+        
+    // Normalize correlations and compute error
+    
+    for (int i=0;i<p.size();i++) { for (int j=0;j<p[i].size();j++) p[i][j] /= (double) (numThreads * numSteps * numRuns); }
+
+}
+
+
+// Compute correlations using rejection-free Monte Carlo (getError_RF) with epitope sites fixed
+
+void getEpitopeCorrelationsP3Pk(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<std::vector<std::vector<std::vector<double> > > > &p3, std::vector<double> &pk, std::vector<int> &latticeStart) {
+
+    // Create and initialize Monte Carlo variables
+
+    for (int i=0;i<p.size();i++) { for (int j=0;j<p[i].size();j++) p[i][j]=0; }
+
+    int N = sizetolength(J.size());
+    
+    std::vector<int> lattice(N,0);
+    std::vector<int> nonzero(N,0);
+    std::vector<int> fixedSites;
+    
+    for (int i=0;i<epitopeStartSite.size();i++) { for (int j=0;j<epitopeStates[i].size();j++) {
+    
+        fixedSites.push_back(epitopeStartSite[i]+j);
+        
+    } }
+    
+    // Neighbors
+    
+    IntVector neighbor(N);
+    getNeighbors(J, N, 1.0e-5, neighbor);
+    
+    // Energies
+    
+    Vector expH;
+    std::vector<double> sumExpH;
+    
+    sumExpH.resize(lattice.size(),0);
+    expH.resize(lattice.size(),std::vector<double>());
+    for (int i=0;i<expH.size();i++) expH[i].resize(expJ[i].size()+1,1);
+    
+    // MC variables
+
+    int stepSize = 4*N;
+    
+    // Set variables for parallel
+    
+    int numThreads=1;
+    
+    srand((unsigned)time(0));
+    unsigned long seed[numThreads];
+    for (int i=0;i<numThreads;i++) seed[i]=rand();
+
+    MT::MersenneTwist mt;
+    mt.init_genrand(seed[0]);
+        
+    // Thermalize and get initial error
+        
+    for (int m=0;m<numRuns;m++) {
+        
+        for (int i=0;i<lattice.size();i++) {
+        
+            lattice[i]=latticeStart[i];
+            
+            if (lattice[i]==J[i].size()) nonzero[i]=0;
+            else                         nonzero[i]=1;
+            
+        }
+        
+        // Fix epitope region
+        int k=0;
+        for (int i=0;i<epitopeStates.size();i++) { for (int j=0;j<epitopeStates[i].size();j++) {
+        
+            lattice[fixedSites[k]]=epitopeStates[i][j];
+            
+            if (lattice[fixedSites[k]]==J[fixedSites[k]].size()) nonzero[fixedSites[k]]=0;
+            else                                                 nonzero[fixedSites[k]]=1;
+            
+            k++;
+            
+        } }
+        
+        setEnergies(expJ, neighbor, lattice, nonzero, expH, sumExpH);
+        getFixedSampleP3Pk(expJ, neighbor, fixedSites, mt, N, numSteps, stepSize, lattice, nonzero, p, p3, pk, expH, sumExpH);
+            
+    }
+        
+    // Normalize correlations and compute error
+    
+    double norm = (double) (numThreads * numSteps * numRuns);
+    
+    for (int i=0;i<p.size();i++) { for (int j=0;j<p[i].size();j++) p[i][j] /= norm; }
+    
+    for (int i=0;i<p3.size();i++) { for (int j=0;j<p3[i].size();j++) { for (int k=0;k<p3[i][j].size();k++) {
+    
+        for (int a=0;a<p3[i][j][k].size();a++) p3[i][j][k][a] /= norm;
+        
+    } } }
+    
+    for (int i=0;i<pk.size();i++) pk[i] /= norm;
+
+}
+
+
+// Get MC sample of Potts configurations using rejection-free Monte Carlo with epitope sites fixed
+
+void getEpitopeStates(FILE *output, 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) {
+
+    // Create and initialize Monte Carlo variables
+
+    for (int i=0;i<p.size();i++) { for (int j=0;j<p[i].size();j++) p[i][j]=0; }
+
+    int N = sizetolength(J.size());
+    
+    std::vector<int> lattice(N,0);
+    std::vector<int> nonzero(N,0);
+    std::vector<int> fixedSites;
+    
+    for (int i=0;i<epitopeStartSite.size();i++) { for (int j=0;j<epitopeStates[i].size();j++) {
+    
+        fixedSites.push_back(epitopeStartSite[i]+j);
+        
+    } }
+    
+    // Neighbors
+    
+    IntVector neighbor(N);
+    getNeighbors(J, N, 1.0e-5, neighbor);
+    
+    // Energies
+    
+    Vector expH;
+    std::vector<double> sumExpH;
+    
+    sumExpH.resize(lattice.size(),0);
+    expH.resize(lattice.size(),std::vector<double>());
+    for (int i=0;i<expH.size();i++) expH[i].resize(expJ[i].size()+1,1);
+    
+    // MC variables
+
+    int stepSize = 4*N;
+    
+    // Set variables for parallel
+    
+    int numThreads=1;
+    
+    srand((unsigned)time(0));
+    unsigned long seed[numThreads];
+    for (int i=0;i<numThreads;i++) seed[i]=rand();
+
+    MT::MersenneTwist mt;
+    mt.init_genrand(seed[0]);
+        
+    // Thermalize and get initial error
+        
+    for (int m=0;m<numRuns;m++) {
+        
+        for (int i=0;i<lattice.size();i++) {
+        
+            lattice[i]=latticeStart[i];
+            
+            if (lattice[i]==J[i].size()) nonzero[i]=0;
+            else                         nonzero[i]=1;
+            
+        }
+        
+        // Fix epitope region
+        int k=0;
+        for (int i=0;i<epitopeStates.size();i++) { for (int j=0;j<epitopeStates[i].size();j++) {
+        
+            lattice[fixedSites[k]]=epitopeStates[i][j];
+            
+            if (lattice[fixedSites[k]]==J[fixedSites[k]].size()) nonzero[fixedSites[k]]=0;
+            else                                                 nonzero[fixedSites[k]]=1;
+            
+            k++;
+            
+        } }
+        
+        setEnergies(expJ, neighbor, lattice, nonzero, expH, sumExpH);
+        writeFixedSample(output, expJ, neighbor, fixedSites, mt, N, numSteps, stepSize, lattice, nonzero, p, expH, sumExpH);
+            
+    }
+
+}
+
+
diff --git a/monteCarlo.h b/monteCarlo.h
new file mode 100755
index 0000000..2eb5c8d
--- /dev/null
+++ b/monteCarlo.h
@@ -0,0 +1,41 @@
+#ifndef MONTECARLO_H
+#define MONTECARLO_H
+
+
+#include <string>
+
+#include "tools.h"  // Numerical tools
+#include "mtrnd.h"  // Random number generator
+
+
+// Error measures
+double epsilonP(const Vector &, const Vector &, int, double, const Vector &, double, double);
+double epsilonP2(const Vector &, const Vector &, int, double, const Vector &, double);
+double epsilonC(const Vector &, const Vector &, int, double, const Vector &, double);
+double getMaxError(const Vector &, const Vector &, const IntVector &, double);
+double getMaxError(const Vector &, const Vector &, double, const Vector &, double, double);
+
+// MC functions
+void setEnergies(const Vector &, const IntVector &, const std::vector<int> &, const std::vector<int> &, Vector &, std::vector<double> &);
+void updateEnergies(const Vector &, const IntVector &, const std::vector<int> &, const std::vector<int> &, int, int, Vector &, std::vector<double> &);
+void dynamicsRF(const Vector &, const IntVector &, int, double, std::vector<int> &, std::vector<int> &, Vector &, std::vector<double> &);
+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 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> &);
+
+// Auxiliary
+void getNeff(const Vector &, int, int &, int &);
+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 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> &);
+
+
+#endif
diff --git a/mtrnd.h b/mtrnd.h
new file mode 100644
index 0000000..62892c6
--- /dev/null
+++ b/mtrnd.h
@@ -0,0 +1,188 @@
+/* 
+A C-program for MT19937, with initialization improved 2002/1/26.
+Coded by Takuji Nishimura and Makoto Matsumoto.
+
+Before using, initialize the state by using init_genrand(seed)  
+or init_by_array(init_key, key_length).
+
+Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
+All rights reserved.                          
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions
+are met:
+
+1. Redistributions of source code must retain the above copyright
+notice, this list of conditions and the following disclaimer.
+
+2. Redistributions in binary form must reproduce the above copyright
+notice, this list of conditions and the following disclaimer in the
+documentation and/or other materials provided with the distribution.
+
+3. The names of its contributors may not be used to endorse or promote 
+products derived from this software without specific prior written 
+permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR
+CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+
+Any feedback is very welcome.
+http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
+email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space)
+
+Reference: M. Matsumoto and T. Nishimura, "Mersenne Twister: 
+A 623-Dimensionally Equidistributed Uniform Pseudo-Random Number Generator",
+ACM Transactions on Modeling and Computer Simulation, Vol. 8, No. 1,
+January 1998, pp 3--30. 
+*/
+
+// C++ wrapped up by Changjiang Yang, more important: make it thread-safe.
+#ifndef _MT_RND_H_
+#define _MT_RND_H_
+
+namespace MT {
+/* Period parameters */  
+const int N = 624;
+const int M = 397;
+const unsigned long MATRIX_A = 0x9908b0dfUL;   /* constant vector a */
+const unsigned long UPPER_MASK = 0x80000000UL; /* most significant w-r bits */
+const unsigned long LOWER_MASK = 0x7fffffffUL; /* least significant r bits */
+
+class MersenneTwist {
+    unsigned long mt[N]; /* the array for the state vector  */
+    int mti; /* mti==N+1 means mt[N] is not initialized */
+public:
+
+    MersenneTwist() : mti(N+1) {}
+
+    /* initializes mt[N] with a seed */
+    //   Note:  Initializing TWISTER to the scalar integer state 0 actually
+    //   corresponds to the C call init_genrand(5489).
+    void init_genrand(unsigned long s)
+    {
+        mt[0]= s & 0xffffffffUL;
+        for (mti=1; mti<N; mti++) {
+            mt[mti] = 
+                (1812433253UL * (mt[mti-1] ^ (mt[mti-1] >> 30)) + mti); 
+            /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
+            /* In the previous versions, MSBs of the seed affect   */
+            /* only MSBs of the array mt[].                        */
+            /* 2002/01/09 modified by Makoto Matsumoto             */
+            mt[mti] &= 0xffffffffUL;
+            /* for >32 bit machines */
+        }
+    }
+
+    /* initialize by an array with array-length */
+    /* init_key is the array for initializing keys */
+    /* key_length is its length */
+    /* slight change for C++, 2004/2/26 */
+    void init_by_array(unsigned long init_key[], int key_length)
+    {
+        int i, j, k;
+        init_genrand(19650218UL);
+        i=1; j=0;
+        k = (N>key_length ? N : key_length);
+        for (; k; k--) {
+            mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525UL))
+                + init_key[j] + j; /* non linear */
+            mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
+            i++; j++;
+            if (i>=N) { mt[0] = mt[N-1]; i=1; }
+            if (j>=key_length) j=0;
+        }
+        for (k=N-1; k; k--) {
+            mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1566083941UL))
+                - i; /* non linear */
+            mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
+            i++;
+            if (i>=N) { mt[0] = mt[N-1]; i=1; }
+        }
+
+        mt[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */ 
+    }
+
+    /* generates a random number on [0,0xffffffff]-interval */
+    unsigned long genrand_int32(void)
+    {
+        unsigned long y;
+        static unsigned long mag01[2]={0x0UL, MATRIX_A};
+        /* mag01[x] = x * MATRIX_A  for x=0,1 */
+
+        if (mti >= N) { /* generate N words at one time */
+            int kk;
+
+            if (mti == N+1)   /* if init_genrand() has not been called, */
+                init_genrand(5489UL); /* a default initial seed is used */
+
+            for (kk=0;kk<N-M;kk++) {
+                y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
+                mt[kk] = mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1UL];
+            }
+            for (;kk<N-1;kk++) {
+                y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
+                mt[kk] = mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1UL];
+            }
+            y = (mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK);
+            mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1UL];
+
+            mti = 0;
+        }
+
+        y = mt[mti++];
+
+        /* Tempering */
+        y ^= (y >> 11);
+        y ^= (y << 7) & 0x9d2c5680UL;
+        y ^= (y << 15) & 0xefc60000UL;
+        y ^= (y >> 18);
+
+        return y;
+    }
+
+    /* generates a random number on [0,0x7fffffff]-interval */
+    long genrand_int31(void)
+    {
+        return (long)(genrand_int32()>>1);
+    }
+
+    /* generates a random number on [0,1]-real-interval */
+    double genrand_real1(void)
+    {
+        return genrand_int32()*(1.0/4294967295.0); 
+        /* divided by 2^32-1 */ 
+    }
+
+    /* generates a random number on [0,1)-real-interval */
+    double genrand_real2(void)
+    {
+        return genrand_int32()*(1.0/4294967296.0); 
+        /* divided by 2^32 */
+    }
+
+    /* generates a random number on (0,1)-real-interval */
+    double genrand_real3(void)
+    {
+        return (((double)genrand_int32()) + 0.5)*(1.0/4294967296.0); 
+        /* divided by 2^32 */
+    }
+
+    /* generates a random number on [0,1) with 53-bit resolution*/
+    double genrand_res53(void) 
+    { 
+        unsigned long a=genrand_int32()>>5, b=genrand_int32()>>6; 
+        return(a*67108864.0+b)*(1.0/9007199254740992.0); 
+    }
+};
+}
+#endif
diff --git a/qEpitopeEval.cpp b/qEpitopeEval.cpp
new file mode 100755
index 0000000..b0acecc
--- /dev/null
+++ b/qEpitopeEval.cpp
@@ -0,0 +1,156 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string>
+#include <math.h>
+#include <cstring>
+#include <iostream>
+#include <sstream>
+#include <cmath>
+#include <vector>
+#include <map>
+#include <algorithm>
+
+#include "qEpitopeEval.h"   // Evaluate epitopes
+#include "monteCarlo.h"     // Monte Carlo
+#include "tools.h"          // Numerical tools
+#include "io.h"             // Data I/O
+#include "mtrnd.h"          // Random number generator
+
+
+// Runs the main program
+
+void runEpitopeEval(RunParameters &r, std::vector<std::vector<int> > &eStates, std::vector<int> &eStartSite, int nout) {
+
+    // Define variables
+    
+    Vector p;
+    Vector J, expJ;
+    std::vector<std::vector<std::vector<std::vector<double> > > > p3;  // 3-point correlations
+    std::vector<double> pk;                                            // mutation probability
+    
+    // Retrieve couplings from file
+    
+    FILE *dataIn=fopen(r.getInfile().c_str(),"r");
+        
+    if (dataIn!=NULL) getCouplings(dataIn,J);
+    else { printf("Error reading input from file %s",r.getInfile().c_str()); exit(1); }
+    
+    // Resize correlations and expJ
+    
+    int N = sizetolength(J.size());
+    
+    pk.resize(N+1, 0);
+    for (int i=0;i<J.size();i++) { p.push_back(std::vector<double>(J[i].size(),0)); expJ.push_back(std::vector<double>(J[i].size(),0)); }
+    for (int i=0;i<J.size();i++) { for (int j=0;j<J[i].size();j++) expJ[i][j] = exp(J[i][j]); }
+    
+    if (r.saveP3) {
+        
+        p3.resize(N);
+        
+        for (int i=0;i<(N-2);i++) {
+          
+            p3[i].resize(N-i);
+            
+            for (int j=i+1;j<(N-1);j++) {
+                
+                p3[i][j-i].resize(N-j);
+                
+                for (int k=j+1;k<N;k++) p3[i][j-i][k-j].resize(p[i].size()*p[j].size()*p[k].size(),0);
+            
+            }
+        
+        }
+    
+    }
+    
+	
+    // Prepare Monte Carlo
+    
+    if (r.useVerbose) printf("Got N=%d, len(h[0])=%d\n",N,(int)J[0].size());
+    
+    std::vector<int> lattice(N);
+    
+    if (r.useStart) {
+    
+        FILE *startIn=fopen(r.getStartInfile().c_str(),"r");
+        for (int i=0;i<N;i++) fscanf(startIn,"%d",&lattice[i]);
+    
+    }
+    
+    else { for (int i=0;i<N;i++) lattice[i]=(int) p[i].size(); }
+    
+    // Prepare to simulate
+    
+    srand((unsigned)time(0));
+    
+//    // Run MC and output Potts configurations
+//    
+//    FILE *dataOut=fopen(r.getCorrelationsOutfile(nout).c_str(),"w");
+//    getEpitopeStates(dataOut, J, expJ, eStates, eStartSite, r.b, r.runs, p, lattice);
+//    fclose(dataOut);
+    
+    // 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);
+    
+    // Write out correlations
+    
+    FILE *dataOut=NULL;
+    
+    if (nout>=0) dataOut=fopen(r.getCorrelationsOutfile(nout).c_str(),"w");
+    else         dataOut=fopen(r.getCorrelationsOutfile().c_str(),    "w");
+    
+    for (int i=0;i<N;i++) {
+    
+        fprintf(dataOut,"%.6e",p[i][0]);
+        for (int j=1;j<p[i].size();j++) fprintf(dataOut,"\t%.6e",p[i][j]);
+        fprintf(dataOut,"\n");
+        
+    }
+    
+    if (r.saveP2) { for (int i=N;i<p.size();i++) {
+    
+        fprintf(dataOut,"%.6e",p[i][0]);
+        for (int j=1;j<p[i].size();j++) fprintf(dataOut,"\t%.6e",p[i][j]);
+        fprintf(dataOut,"\n");
+        
+    } }
+    
+    fclose(dataOut);
+    
+    // Write out extended correlations
+    
+    FILE *p3Out=NULL;
+    
+    if (r.saveP3) {
+    
+        if (nout>=0) p3Out=fopen(r.getP3Outfile(nout).c_str(),"w");
+        else         p3Out=fopen(r.getP3Outfile().c_str(),    "w");
+    
+        for (int i=0;i<p3.size();i++) { for (int j=0;j<p3[i].size();j++) { for (int k=0;k<p3[i][j].size();k++) {
+    
+            for (int a=0;a<p3[i][j][k].size();a++) fprintf(p3Out,"%.6e\n",p3[i][j][k][a]);
+        
+        } } }
+        
+    }
+    
+    fclose(p3Out);
+    
+    
+    FILE *pkOut=NULL;
+    
+    if (r.savePk) {
+    
+        if (nout>=0) pkOut=fopen(r.getPkOutfile(nout).c_str(),"w");
+        else         pkOut=fopen(r.getPkOutfile().c_str(),    "w");
+    
+        for (int i=0;i<pk.size();i++) fprintf(pkOut,"%.6e\n",pk[i]);
+        
+    }
+    
+    fclose(pkOut);
+    
+}
+
diff --git a/qEpitopeEval.h b/qEpitopeEval.h
new file mode 100755
index 0000000..fcc4c67
--- /dev/null
+++ b/qEpitopeEval.h
@@ -0,0 +1,71 @@
+#ifndef QLEARNSPARSE_H
+#define QLEARNSPARSE_H
+
+
+#include <iostream>
+#include <string>
+
+#include "tools.h"  // Numerical tools
+
+
+// This class holds the parameters needed for running the ising program
+
+class RunParameters {
+    
+public:
+    
+    std::string directory;          // Path to the directory where the inut file is located
+                                    // Output will also be sent to this directory
+    std::string infile;             // Input file from which couplings are to be read
+    std::string outfile;            // Output file (prefix) where data is to be written
+    
+    double b;                       // Number of data points to take in a single run
+    double runs;                    // Number of times to run dynamics
+    
+    std::vector<std::vector<std::vector<int> > > epitopeStates; // Configuration of sites in the targeted epitope
+    std::vector<std::vector<int> >            epitopeStartSite; // Index of the first site in the epitope
+    
+    int on;                         // Start number for keeping track of all epitopes and epitope combinations
+    
+    bool useStart;                  // Use different starting sequence than all wild-type
+    std::string startfile;          // File containing the starting sequence
+    
+    bool useVerbose;                // If true, print extra information while program is running
+    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
+    
+    RunParameters() {
+        
+        directory=".";
+        infile="input";
+        outfile="output";
+        b=1000000;
+        runs=1;
+        on=0;
+        
+        useStart=false;
+        useVerbose=false;
+        saveP2=false;
+        saveP3=false;
+        savePk=false;
+        
+    }
+    std::string getInfile()                   { return (directory+"/"+infile+".j");                   }
+    std::string getCorrelationsOutfile()      { return (directory+"/"+outfile+".p");                  }
+    std::string getCorrelationsOutfile(int i) { return (directory+"/"+outfile+"-"+tostring(i)+".p");  }
+    std::string getP3Outfile()                { return (directory+"/"+outfile+".p3");                 }
+    std::string getP3Outfile(int i)           { return (directory+"/"+outfile+"-"+tostring(i)+".p3"); }
+    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");              }
+    ~RunParameters() {}
+    
+};
+
+
+//void runEpitopeEval(RunParameters &);
+void runEpitopeEval(RunParameters &, std::vector<std::vector<int> > &, std::vector<int> &, int);
+
+
+#endif
diff --git a/tools.cpp b/tools.cpp
new file mode 100755
index 0000000..68bfacd
--- /dev/null
+++ b/tools.cpp
@@ -0,0 +1,1387 @@
+#include <vector>
+#include <math.h>
+
+#include "tools.h"  // Numerical tools
+
+
+// Calculates gamma from the data for the L0 norm regularization
+
+double computeGamma_L0(const Vector &p, double B) { return 1/B; }
+
+
+// Calculates gamma from the data for the L2 norm regularization
+
+double computeGamma_L2(const Vector &p, double B) { return 1/B; }
+
+
+// Returns the inner product of two arrays of doubles
+
+double innerProduct(double list1[], double list2[], size_t length) {
+
+    double product=0;
+    
+    for (int i=0;i<length;i++) product += list1[i] * list2[i];
+    
+    return product;
+    
+}
+
+
+// Returns the inner product of two vectors of doubles
+
+double innerProduct(const std::vector<double> &list1, const std::vector<double> &list2) {
+
+    double product=0;
+    
+    for (int i=0;i<list1.size();i++) product += list1[i] * list2[i];
+    
+    return product;
+    
+}
+
+
+// Returns the inner product of two Vectors of doubles
+
+double innerProduct(const Vector &list1, const Vector &list2) {
+
+    double product=0;
+    
+    for (int i=0;i<list1.size();i++) {
+    
+        for (int j=0;j<list1[i].size();j++) product += list1[i][j] * list2[i][j];
+        
+    }
+    
+    return product;
+    
+}
+
+
+// Returns the inner product of two Vectors of doubles (sparse)
+
+double innerProduct(const Vector &list1, const Vector &list2, const IntVector &sparseSet) {
+
+    double product=0;
+    
+    for (int i=0;i<sparseSet.size();i++) {
+    
+        for (int j=0;j<sparseSet[i].size();j++) product += list1[i][sparseSet[i][j]] * list2[i][sparseSet[i][j]];
+        
+    }
+    
+    return product;
+    
+}
+
+
+// Computes the L1 norm of an array of doubles
+
+double L1(double list[], size_t length) {
+    
+    double norm=0;
+    
+    for (int i=0;i<length;i++) norm += fabs(list[i]);
+    
+    return norm;
+    
+}
+
+
+// Computes the L1 norm of a vector of doubles
+
+double L1(const std::vector<double> &list) {
+    
+    double norm=0;
+    
+    for (int i=0;i<list.size();i++) norm += fabs(list[i]);
+    
+    return norm;
+    
+}
+
+
+// Computes the L1 norm of a Vector of doubles
+
+double L1(const Vector &list) {
+    
+    double norm=0;
+    
+    for (int i=0;i<list.size();i++) {
+    
+        for (int j=0;j<list[i].size();j++) norm += fabs(list[i][j]);
+        
+    }
+    
+    return norm;
+    
+}
+
+
+// Computes the L2 norm of an array of doubles
+
+double L2(double list[], size_t length) {
+    
+    double norm=0;
+    
+    for (int i=0;i<length;i++) norm += pow(list[i],2.0);
+    
+    return sqrt(norm);
+    
+}
+
+
+// Computes the L2 norm of a vector of doubles
+
+double L2(const std::vector<double> &list) {
+    
+    double norm=0;
+    
+    for (int i=0;i<list.size();i++) norm += pow(list[i],2.0);
+    
+    return sqrt(norm);
+    
+}
+
+
+// Computes the L2 norm of a Vector of doubles
+
+double L2(const Vector &list) {
+    
+    double norm=0;
+    
+    for (int i=0;i<list.size();i++) {
+    
+        for (int j=0;j<list[i].size();j++) norm += pow(list[i][j],2.0);
+        
+    }
+    
+    return sqrt(norm);
+    
+}
+
+
+// Computes the L2 norm of a Vector of doubles (sparse version)
+
+double L2(const Vector &list, const IntVector &sparseSet) {
+    
+    double norm=0;
+    
+    for (int i=0;i<sparseSet.size();i++) {
+    
+        for (int j=0;j<sparseSet[i].size();j++) norm += pow(list[i][sparseSet[i][j]],2.0);
+        
+    }
+    
+    return sqrt(norm);
+    
+}
+
+
+// Computes the LInfinity norm of an array of doubles
+
+double LInfinity(double list[], size_t length) {
+    
+    double norm=0;
+    
+    for (int i=0;i<length;i++) { if (fabs(list[i])>norm) norm = fabs(list[i]); }
+    
+    return norm;
+    
+}
+
+
+// Computes the LInfinity norm of a vector of doubles
+
+double LInfinity(const std::vector<double> &list) {
+    
+    double norm=0;
+    
+    for (int i=0;i<list.size();i++) { if (fabs(list[i])>norm) norm = fabs(list[i]); }
+    
+    return norm;
+    
+}
+
+
+// Computes the LInfinity norm of a vector of doubles
+
+double LInfinity(const Vector &list) {
+    
+    double norm=0;
+    
+    for (int i=0;i<list.size();i++) {
+    
+        for (int j=0;j<list[i].size();j++) {
+        
+            if (fabs(list[i][j])>norm) norm = fabs(list[i][j]);
+            
+        }
+        
+    }
+    
+    return norm;
+    
+}
+
+
+// Computes the LInfinity norm of a vector of doubles (sparse)
+
+double LInfinity(const Vector &list, const IntVector &sparseSet) {
+    
+    double norm=0;
+    
+    for (int i=0;i<sparseSet.size();i++) {
+    
+        for (int j=0;j<sparseSet[i].size();j++) {
+        
+            if (fabs(list[i][sparseSet[i][j]])>norm) norm = fabs(list[i][sparseSet[i][j]]);
+            
+        }
+        
+    }
+    
+    return norm;
+    
+}
+
+
+// Insert an element into vector in order (smallest to largest)
+
+void insertInPlace(std::vector<int> &list, int element) {
+
+    if (list.empty()) list.push_back(element);
+    
+    else {
+    
+        bool insert=true;
+        int  imin=-1;
+        int  imax=(int)list.size();
+    
+        while (imax-imin>1) {
+      
+            int imid = (imin+imax)/2;
+ 
+            if      (list[imid] <  element) imin = imid;
+            else if (list[imid] >  element) imax = imid;
+            else if (list[imid] == element) { insert=false; break; }
+        
+        }
+        
+        if (insert) list.insert(list.begin()+imax,element);
+        
+    }
+
+}
+
+
+// Insert an element into vector in order (smallest to largest)
+
+void insertInPlace(int list[], int length, int element) {
+
+    if (length==1) list[0]=element;
+    
+    else {
+    
+        bool insert=true;
+        int  imin=-1;
+        int  imax=length;
+    
+        while (imax-imin>1) {
+      
+            int imid = (imin+imax)/2;
+ 
+            if      (list[imid] <  element) imin = imid;
+            else if (list[imid] >  element) imax = imid;
+            else if (list[imid] == element) { insert=false; break; }
+        
+        }
+        
+        if (insert) {
+        
+            for (int i=length;i>imax;i--) list[i] = list[i-1];
+            list[imax] = element;
+            
+        }
+        
+    }
+
+}
+
+
+// Erase an element in a vector, maintaining order (smallest to largest)
+
+void eraseInPlace(std::vector<int> &list, int element) {
+
+    bool erase=true;
+    int  imin=-1;
+    int  imax=(int)list.size();
+    int  imid;
+    
+    while (imax-imin>1) {
+      
+        imid = (imin+imax)/2;
+ 
+        if      (list[imid] <  element) imin = imid;
+        else if (list[imid] >  element) imax = imid;
+        else if (list[imid] == element) { erase=true; break; }
+        
+    }
+        
+    if (erase) list.erase(list.begin()+imid);
+
+
+}
+
+
+// Erase an element in a vector, maintaining proper order
+
+void eraseInPlace(int list[], int length, int element) {
+
+    bool erase=true;
+    int  imin=-1;
+    int  imax=length;
+    int  imid;
+    
+    while (imax-imin>1) {
+      
+        imid = (imin+imax)/2;
+ 
+        if      (list[imid] <  element) imin = imid;
+        else if (list[imid] >  element) imax = imid;
+        else if (list[imid] == element) { erase=true; break; }
+        
+    }
+        
+    if (erase) {
+    
+        for (int i=imid;i<length-1;i++) list[i]=list[i+1];
+        
+    }
+
+
+}
+
+
+// Sorts the list of doubles (with corresponding index also arranged properly) using the quickSort algorithm
+// Sort order is FROM HIGH TO LOW
+
+void quickSort(double a[][2], int left, int right) {
+    
+    int i=left;
+    int j=right;
+    
+    double temp0, temp1;
+    double pivot = a[(left + right) / 2][0];
+    
+    // Partition
+    
+    while (i <= j) {
+        
+        while (a[i][0] > pivot) i++;
+        while (a[j][0] < pivot) j--;
+        
+        if (i <= j) {
+            
+            temp0 = a[i][0];
+            temp1 = a[i][1];
+            
+            a[i][0] = a[j][0];
+            a[i][1] = a[j][1];
+            
+            a[j][0] = temp0;
+            a[j][1] = temp1;
+            
+            i++;
+            j--;
+
+        }
+
+    }
+
+    // Recursion
+    
+    if (left < j)  quickSort(a, left, j);
+    if (i < right) quickSort(a, i, right);
+
+}
+
+
+// Sorts the list of doubles (with corresponding index also arranged properly) using the quickSort algorithm
+// Sort order is FROM HIGH TO LOW
+
+void quickSort(double a[][3], int left, int right) {
+    
+    int i=left;
+    int j=right;
+    
+    double temp0, temp1, temp2;
+    double pivot = a[(left + right) / 2][0];
+    
+    // Partition
+    
+    while (i <= j) {
+        
+        while (a[i][0] > pivot) i++;
+        while (a[j][0] < pivot) j--;
+        
+        if (i <= j) {
+            
+            temp0 = a[i][0];
+            temp1 = a[i][1];
+            temp2 = a[i][2];
+            
+            a[i][0] = a[j][0];
+            a[i][1] = a[j][1];
+            a[i][2] = a[j][2];
+            
+            a[j][0] = temp0;
+            a[j][1] = temp1;
+            a[j][2] = temp2;
+            
+            i++;
+            j--;
+
+        }
+
+    }
+
+    // Recursion
+    
+    if (left < j)  quickSort(a, left, j);
+    if (i < right) quickSort(a, i, right);
+
+}
+
+
+// Sorts the vector of integers (with corresponding index also arranged properly) using the quickSort algorithm
+// Sort order is FROM LOW TO HIGH
+
+void quickSort(std::vector<int> &a, int left, int right) {
+    
+    int i=left;
+    int j=right;
+    
+    int temp;
+    int pivot = a[(left + right) / 2];
+    
+    // Partition
+    
+    while (i <= j) {
+        
+        while (a[i] < pivot) i++;
+        while (a[j] > pivot) j--;
+        
+        if (i <= j) {
+            
+            temp = a[i];
+            a[i] = a[j];
+            a[j] = temp;
+            
+            i++;
+            j--;
+            
+        }
+        
+    }
+    
+    // Recursion
+    
+    if (left < j)  quickSort(a, left, j);
+    if (i < right) quickSort(a, i, right);
+    
+}
+
+
+// Check to see if two sets of spins of size k have a superset of size k+1
+
+bool checkSize(const Key &a, const Key &b) {
+    
+    int i=0,j=0,error=0;
+    
+    while (error<3) {
+        
+        if      (i==(int)a.size()) { error += (int)b.size()-j; break; }
+		else if (j==(int)b.size()) { error += (int)a.size()-i; break; }
+        
+		else if (a[i]<b[j])   { i++; error++; continue; }
+		else if (a[i]>b[j])   { j++; error++; continue; }
+        else if (a[i]==b[j])  { i++; j++;     continue; }
+		else break;
+        
+	}
+    
+    if (error!=2) return false;
+    else          return true;
+    
+}
+
+
+// Returns the superset of two sets of spins
+
+void getSuperset(const Key &a, const Key &b, Key &superset) {
+	
+	int i=0,j=0,index=0;
+	
+	while (true) {
+		
+		if      (i==a.size()) { for (int k=j;k<b.size();k++) { superset[index]=b[k]; index++; } break; }
+		else if (j==b.size()) { for (int k=i;k<a.size();k++) { superset[index]=a[k]; index++; } break; }
+        
+		else if (a[i]==b[j]) { superset[index]=a[i]; index++; i++; j++; continue; }
+		else if (a[i]<b[j])  { superset[index]=a[i]; index++; i++;      continue; }
+        else if (a[i]>b[j])  { superset[index]=b[j]; index++;      j++; continue; }
+		else break;
+        
+	}
+	
+}
+
+
+// Take a Hessian matrix in Vector format and "unpack" it into the form of a list of upper triangular
+// matrix entries in linearHess. The shape of the coupling Vector is used to help unpack the Hessian.
+
+void unpack(Vector &hess, std::vector<double> &linearHess, const Vector &J) {
+
+    int count=0;
+    
+    for (int i=0;i<J.size();i++) {
+    
+        for (int j=0;j<J[i].size();j++) {
+        
+            for (int l=j;l<J[i].size();l++) {     linearHess[count] = hess[hindex(i,i,J.size())][hindex(j,l,J[i].size())];             count++; }
+        
+            for (int k=i+1;k<J.size();k++) {
+            
+                for (int l=0;l<J[k].size();l++) { linearHess[count] = hess[hindex(i,k,J.size())][sindex(j,l,J[i].size(),J[k].size())]; count++; }
+            
+            }
+            
+        }
+        
+    }
+
+}
+
+
+// Perform the inverse of the unpack function above.
+
+void pack(Vector &hess, std::vector<double> &linearHess, const Vector &J) {
+
+    int count=0;
+    
+    for (int i=0;i<J.size();i++) {
+    
+        for (int j=0;j<J[i].size();j++) {
+        
+            for (int l=j;l<J[i].size();l++) {     hess[hindex(i,i,J.size())][hindex(j,l,J[i].size())] = linearHess[count];             count++; }
+        
+            for (int k=i+1;k<J.size();k++) {
+            
+                for (int l=0;l<J[k].size();l++) { hess[hindex(i,k,J.size())][sindex(j,l,J[i].size(),J[k].size())] = linearHess[count]; count++; }
+            
+            }
+            
+        }
+        
+    }
+
+}
+
+
+// Take a Hessian matrix in Vector format and "unpack" it into the form of a list of upper triangular
+// matrix entries in linearHess. The shape of the coupling Vector is used to help unpack the Hessian.
+// Sparse version for L1 regularization.
+
+void unpack_sparse(Vector &hess, std::vector<double> &linearHess, const Vector &J, const IntVector &workingSet) {
+
+    int count=0;
+    
+    for (int i=0;i<workingSet.size();i++) {
+    
+        for (int j=0;j<workingSet[i].size();j++) {
+        
+            for (int l=j;l<workingSet[i].size();l++) {     linearHess[count] = hess[hindex(i,i,J.size())][hindex(workingSet[i][j],workingSet[i][l],J[i].size())];             count++; }
+            
+            for (int k=i+1;k<workingSet.size();k++) {
+            
+                for (int l=0;l<workingSet[k].size();l++) { linearHess[count] = hess[hindex(i,k,J.size())][sindex(workingSet[i][j],workingSet[k][l],J[i].size(),J[k].size())]; count++; }
+            
+            }
+        
+        }
+    
+    }
+    
+    linearHess.resize(count);
+
+}
+
+
+// Perform the inverse of the unpack function above.
+// Sparse version for L1 regularization.
+
+void pack_sparse(Vector &hess, std::vector<double> &linearHess, const Vector &J, const IntVector &workingSet) {
+
+    int count=0;
+    
+    for (int i=0;i<workingSet.size();i++) {
+    
+        for (int j=0;j<workingSet[i].size();j++) {
+        
+            for (int l=j;l<workingSet[i].size();l++) {     hess[hindex(i,i,J.size())][hindex(workingSet[i][j],workingSet[i][l],J[i].size())] = linearHess[count];             count++; }
+            
+            for (int k=i+1;k<workingSet.size();k++) {
+            
+                for (int l=0;l<workingSet[k].size();l++) { hess[hindex(i,k,J.size())][sindex(workingSet[i][j],workingSet[k][l],J[i].size(),J[k].size())] = linearHess[count]; count++; }
+            
+            }
+        
+        }
+    
+    }
+
+}
+
+
+// This function performs the modified Cholesky factorization on a positive-definite, but
+// possibly near singular matrix and uses this to invert the matrix. The inversion is regularized
+// to insure that the condition number of the matrix does not exceed a certain maximum value.
+
+void modifiedCholeskyInverse(double m[], size_t length) {
+    
+    double c[length][length];
+    double d[length];
+    double beta =pow(10.0,  2.0);
+    // Note: comp should be = beta * sqrt(delta)
+    double comp =pow(10.0, -3.0);
+    
+    // First compute modified LDLT Cholesky factorization
+    
+    for (int i=0;i<length;i++) {
+        
+        c[i][i] = m[hindex(i,i,length)];
+        
+        for (int j=0;j<i;j++) c[i][i] -= d[j] * pow(c[j][i],2);
+        
+        double max=comp;
+        
+        for (int j=i+1;j<length;j++) {
+            
+            c[i][j] = m[hindex(i,j,length)];
+            
+            for (int k=0;k<i;k++) c[i][j] -= d[k] * c[k][i] * c[k][j];
+            
+            if (fabs(c[i][j])>max) max=fabs(c[i][j]);
+            
+        }
+        
+        max = pow(max/beta,2);
+        
+        if (fabs(c[i][i])>max) d[i] = fabs(c[i][i]);
+        else                   d[i] = max;
+        
+        for (int j=i;j<length;j++) c[i][j] /= d[i];
+            
+    }
+    
+    // Now take the inverse
+    
+    double n[length][length];
+    
+    // First get D^-1 and L^-1
+    
+    for (int i=0;i<length;i++) {
+        
+        n[i][i] = 1;
+        
+        for (int j=i+1;j<length;j++) {
+            
+            n[i][j] = -c[i][j];
+            
+            for (int k=i+1;k<j;k++) n[i][j] -= n[i][k] * c[k][j];
+            
+        }
+        
+    }
+    
+    // Multiply out (L^-1)T D^-1 L^-1 to get the inverse
+    
+    for (int i=0;i<length;i++) {
+        
+        for (int j=i;j<length;j++) {
+            
+            m[hindex(i,j,length)]=0;
+            
+            for (int k=j;k<length;k++) m[hindex(i,j,length)] += n[i][k] * n[j][k] / d[k];
+            
+        }
+        
+    }
+    
+}
+
+
+// This function performs the modified Cholesky factorization on a positive-definite, but
+// possibly near singular matrix and uses this to invert the matrix. The inversion is regularized
+// to insure that the condition number of the matrix does not exceed a certain maximum value.
+
+void modifiedCholeskyInverse(Vector &hess, std::vector<double> &linearHess, int length, const Vector &J) {
+
+    unpack(hess, linearHess, J);
+    modifiedCholeskyInverse(linearHess, length);
+    pack(hess, linearHess, J);
+    
+}
+
+
+// This function performs the modified Cholesky factorization on a positive-definite, but
+// possibly near singular matrix and uses this to invert the matrix. The inversion is regularized
+// to insure that the condition number of the matrix does not exceed a certain maximum value.
+// Sparse version for L1 regularization.
+
+void modifiedCholeskyInverse_sparse(Vector &hess, std::vector<double> &linearHess, std::vector<double> &inverseC, std::vector<double> &inverseN, std::vector<double> &inverseD, const Vector &J, const IntVector &workingSet) {
+
+    unpack_sparse(hess, linearHess, J, workingSet);
+    modifiedCholeskyInverse(linearHess, inverseC, inverseN, inverseD, inverseD.size());
+    pack_sparse(hess, linearHess, J, workingSet);
+    
+}
+
+
+// This function performs the modified Cholesky factorization on a positive-definite, but
+// possibly near singular matrix and uses this to invert the matrix. The inversion is regularized
+// to insure that the condition number of the matrix does not exceed a certain maximum value.
+
+void modifiedCholeskyInverse(std::vector<double> &m, std::vector<double> &c, std::vector<double> &n, std::vector<double> &d, size_t length) {
+    
+    //double delta=pow(10.0,-10.0);
+    double beta =pow(10.0,  2.0);
+    // Note: comp should be = beta * sqrt(delta)
+    double comp =pow(10.0, -3.0);
+    
+    // First compute modified LDLT Cholesky factorization
+    
+    for (int i=0;i<length;i++) {
+    
+        d[i]=0;
+        
+        c[hindex(i,i,length)] = m[hindex(i,i,length)];
+        
+        for (int j=0;j<i;j++) c[hindex(i,i,length)] -= d[j] * pow(c[hindex(j,i,length)],2);
+        
+        double max=comp;
+        
+        for (int j=i+1;j<length;j++) {
+            
+            c[hindex(i,j,length)] = m[hindex(i,j,length)];
+            
+            for (int k=0;k<i;k++) c[hindex(i,j,length)] -= d[k] * c[hindex(k,i,length)] * c[hindex(k,j,length)];
+            
+            if (fabs(c[hindex(i,j,length)])>max) max = fabs(c[hindex(i,j,length)]);
+            
+        }
+        
+        max = pow(max/beta,2);
+        
+        if (fabs(c[hindex(i,i,length)])>max) d[i] = fabs(c[hindex(i,i,length)]);
+        else                                 d[i] = max;
+        
+        for (int j=i;j<length;j++) c[hindex(i,j,length)] /= d[i];
+            
+    }
+    
+    // First get D^-1 and L^-1
+    
+    for (int i=0;i<length;i++) {
+        
+        n[hindex(i,i,length)] = 1;
+        
+        for (int j=i+1;j<length;j++) {
+            
+            n[hindex(i,j,length)] = -c[hindex(i,j,length)];
+            
+            for (int k=i+1;k<j;k++) n[hindex(i,j,length)] -= n[hindex(i,k,length)] * c[hindex(k,j,length)];
+            
+        }
+        
+    }
+    
+    // Multiply out (L^-1)T D^-1 L^-1 to get the inverse
+    
+    for (int i=0;i<length;i++) {
+        
+        for (int j=i;j<length;j++) {
+            
+            m[hindex(i,j,length)]=0;
+            
+            for (int k=j;k<length;k++) m[hindex(i,j,length)] += n[hindex(i,k,length)] * n[hindex(j,k,length)] / d[k];
+            
+        }
+        
+    }
+    
+}
+
+
+// This function performs the modified Cholesky factorization on a positive-definite, but
+// possibly near singular matrix and uses this to invert the matrix. The inversion is regularized
+// to insure that the condition number of the matrix does not exceed a certain maximum value.
+
+void modifiedCholeskyInverse(std::vector<double> &m, int length) {
+    
+    std::vector<std::vector<double> > c;
+    c.resize(length,std::vector<double>(length,0));
+    std::vector<double> d(length,0);
+    
+    double beta =pow(10.0,  2.0);
+    // Note: comp should be = beta * sqrt(delta)
+    double comp =pow(10.0, -3.0);
+    
+    // First compute modified LDLT Cholesky factorization
+    
+    for (int i=0;i<length;i++) {
+        
+        c[i][i] = m[hindex(i,i,length)];
+        
+        for (int j=0;j<i;j++) c[i][i] -= d[j] * pow(c[j][i],2);
+        
+        double max=comp;
+        
+        for (int j=i+1;j<length;j++) {
+            
+            c[i][j] = m[hindex(i,j,length)];
+            
+            for (int k=0;k<i;k++) c[i][j] -= d[k] * c[k][i] * c[k][j];
+            
+            if (fabs(c[i][j])>max) max=fabs(c[i][j]);
+            
+        }
+        
+        max = pow(max/beta,2);
+        
+        if (fabs(c[i][i])>max) d[i] = fabs(c[i][i]);
+        else                   d[i] = max;
+        
+        for (int j=i;j<length;j++) c[i][j] /= d[i];
+            
+    }
+    
+    // Now take the inverse
+    
+    std::vector<std::vector<double> > n;
+    n.resize(length,std::vector<double>(length,0));
+    
+    // First get D^-1 and L^-1
+    
+    for (int i=0;i<length;i++) {
+        
+        n[i][i] = 1;
+        
+        for (int j=i+1;j<length;j++) {
+            
+            n[i][j] = -c[i][j];
+            
+            for (int k=i+1;k<j;k++) n[i][j] -= n[i][k] * c[k][j];
+            
+        }
+        
+    }
+    
+    // Multiply out (L^-1)T D^-1 L^-1 to get the inverse
+    
+    for (int i=0;i<length;i++) {
+        
+        for (int j=i;j<length;j++) {
+            
+            m[hindex(i,j,length)]=0;
+            
+            for (int k=j;k<length;k++) m[hindex(i,j,length)] += n[i][k] * n[j][k] / d[k];
+            
+        }
+        
+    }
+    
+}
+
+
+// Same as above, but also returns the determinant
+
+double modifiedCholeskyInverseDet(double m[], size_t length) {
+    
+    double c[length][length];
+    double d[length];
+    double beta =pow(10.0,  2.0);
+    // Note: comp should be = beta * sqrt(delta)
+    double comp =pow(10.0, -3.0);
+    
+    // First compute modified LDLT Cholesky factorization
+    
+    for (int i=0;i<length;i++) {
+        
+        c[i][i] = m[hindex(i,i,length)];
+        
+        for (int j=0;j<i;j++) c[i][i] -= d[j] * pow(c[j][i],2);
+        
+        double max=comp;
+        
+        for (int j=i+1;j<length;j++) {
+            
+            c[i][j] = m[hindex(i,j,length)];
+            
+            for (int k=0;k<i;k++) c[i][j] -= d[k] * c[k][i] * c[k][j];
+            
+            if (fabs(c[i][j])>max) max=fabs(c[i][j]);
+            
+        }
+        
+        max = pow(max/beta,2);
+        
+        if (fabs(c[i][i])>max) d[i] = fabs(c[i][i]);
+        else                   d[i] = max;
+        
+        for (int j=i;j<length;j++) c[i][j] /= d[i];
+        
+    }
+    
+    // Compute the determinant
+    
+    double det=1;
+    
+    for (int i=0;i<length;i++) det *= d[i];
+    
+    // Now take the inverse
+    
+    double n[length][length];
+    
+    // First get D^-1 and L^-1
+    
+    for (int i=0;i<length;i++) {
+        
+        n[i][i] = 1;
+        
+        for (int j=i+1;j<length;j++) {
+            
+            n[i][j] = -c[i][j];
+            
+            for (int k=i+1;k<j;k++) n[i][j] -= n[i][k] * c[k][j];
+            
+        }
+        
+    }
+    
+    // Multiply out (L^-1)T D^-1 L^-1 to get the inverse
+    
+    for (int i=0;i<length;i++) {
+        
+        for (int j=i;j<length;j++) {
+            
+            m[hindex(i,j,length)]=0;
+            
+            for (int k=j;k<length;k++) m[hindex(i,j,length)] += n[i][k] * n[j][k] / d[k];
+            
+        }
+        
+    }
+    
+    return det;
+    
+}
+
+
+// Returns the determinant of a symmetric matrix
+
+double symmetricDet(double m[], int length) {
+    
+    double c[length][length];
+    double d[length];
+    double beta =pow(10.0,  2.0);
+    // Note: comp should be = beta * sqrt(delta)
+    double comp =pow(10.0, -3.0);
+    
+    // First compute modified LDLT Cholesky factorization
+    
+    for (int i=0;i<length;i++) {
+        
+        c[i][i] = m[hindex(i,i,length)];
+        
+        for (int j=0;j<i;j++) c[i][i] -= d[j] * pow(c[j][i],2);
+        
+        double max=comp;
+        
+        for (int j=i+1;j<length;j++) {
+            
+            c[i][j] = m[hindex(i,j,length)];
+            
+            for (int k=0;k<i;k++) c[i][j] -= d[k] * c[k][i] * c[k][j];
+            
+            if (fabs(c[i][j])>max) max=fabs(c[i][j]);
+            
+        }
+        
+        max = pow(max/beta,2);
+        
+        if (fabs(c[i][i])>max) d[i] = fabs(c[i][i]);
+        else                   d[i] = max;
+        
+        for (int j=i;j<length;j++) c[i][j] /= d[i];
+        
+    }
+    
+    // Compute the determinant
+    
+    double det=1;
+    
+    for (int i=0;i<length;i++) det *= d[i];
+    
+    return det;
+    
+}
+
+
+// Given a vector x of length n, this function computes parameters for a Householder matrix H,
+// such that H = 1 - beta v v^T is orthogonal and H x = L2(x) e_1. See Matrix Computations (Golub) for details.
+
+void householder(double x[], int length, double v[], double &beta) {
+    
+    double sigma = 0;
+    
+    v[0] = 1;
+    
+    for (int i=1;i<length;i++) { sigma += pow(x[i],2); v[i] = x[i]; }
+    
+    if (sigma==0) beta=0;
+    
+    else {
+        
+        double mu = sqrt(pow(x[0],2) + sigma);
+        
+        if (x[0]<=0) v[0] = x[0] - mu;
+        else         v[0] = -sigma / (x[0] + mu);
+        
+        beta = 2 * pow(v[0],2) / (sigma + pow(v[0],2));
+        
+        for (int i=length-1;i>=0;i--)  v[i] /= v[0];
+        
+    }
+    
+}
+
+
+// Given a symmetric matrix A, this algorithm overwrites A with T = Q^T A Q, with T tridiagonal,
+// and Q orthogonal. See Matrix Computations (Golub) for details. This algorithm assumes Q
+// has already been initialized.
+
+void tridiagonalize(double A[], int length, std::vector<std::vector<double> > &Q) {
+    
+    for (int k=0;k<length-2;k++) {
+        
+        // Compute the Householder matrix
+        
+        double v[length-(k+1)], a[length-(k+1)];
+        double beta;
+        
+        for (int i=k+1;i<length;i++) a[i-(k+1)] = A[hindex(k,i,length)];
+        
+        householder(a,length-(k+1),v,beta);
+        
+        // Update A
+        
+        double p[length-(k+1)], w[length-(k+1)];
+        double pv=0;
+        
+        for (int i=k+1;i<length;i++) {
+            
+            p[i-(k+1)]                             = 0;
+            for (int j=k+1;j<i;j++)    p[i-(k+1)] += A[hindex(j,i,length)] * v[j-(k+1)];
+            for (int j=i;j<length;j++) p[i-(k+1)] += A[hindex(i,j,length)] * v[j-(k+1)];
+            p[i-(k+1)]                            *= beta;
+            
+            w[i-(k+1)]=p[i-(k+1)];
+            
+            pv += p[i-(k+1)] * v[i-(k+1)];
+            
+        }
+        
+        for (int i=0;i<length-(k+1);i++) w[i] -= beta * pv * v[i] / 2;
+        
+        A[hindex(k,k+1,length)]=L2(a,length-(k+1));
+        
+        for (int i=k+1;i<length;i++) {
+            
+            for (int j=i;j<length;j++)  A[hindex(i,j,length)] -= v[i-(k+1)] * w[j-(k+1)] + v[j-(k+1)] * w[i-(k+1)];
+            
+        }
+        
+        // Update Q
+        
+        double ww[length];
+        
+        for (int i=0;i<length;i++) {
+            
+            ww[i]                               = 0;
+            for (int j=k+1;j<length;j++) ww[i] += Q[i][j] * v[j-(k+1)];
+            ww[i]                              *= beta;
+            
+            
+        }
+        
+        for (int i=0;i<length;i++) {
+            
+            for (int j=k+1;j<length;j++)  Q[i][j] -= ww[i] * v[j-(k+1)];
+            
+        }
+        
+    }
+    
+    for (int i=0;i<length;i++) {
+        
+        for (int j=i+2;j<length;j++) A[hindex(i,j,length)]=0;
+        
+    }
+    
+}
+
+
+// Given an unreduced symmetric tridiagonal matrix T, this algorithm overwrites T with Z^T T Z,
+// where Z is a product of Givens rotations. See Matrix Computations (Golub) for details.
+
+void QRStep(double T[], int length, std::vector<std::vector<double> > &Z) {
+    
+    double d = (T[hindex(length-2,length-2,length)] - T[hindex(length-1,length-1,length)]) / 2;
+    
+    double mu = T[hindex(length-1,length-1,length)] - pow(T[hindex(length-2,length-1,length)],2)
+                     / (d * (1 + sqrt(pow(d,2) + pow(T[hindex(length-2,length-1,length)],2)) / fabs(d) ) );
+    
+    double x = T[0] - mu;
+    double z = T[hindex(0,1,length)];
+    
+    for (int k=0;k<length-1;k++) {
+        
+        double c, s;
+        givens(x, z, c, s);
+        
+        double G[4][4] = {{1, 0, 0, 0}, {0, c, s, 0}, {0, -s, c, 0}, {0, 0, 0, 1}};
+        double newT[(length*(length+1))/2];
+        
+        // Multiply T by G on the right and G^T on the left
+        
+        // Start with i=k-1 --> m=k-1 
+        
+        if (k>0) {
+        
+            for (int j=k-1;j<k+2;j++) {
+            
+                                          newT[hindex(k-1,j,length)]  = 0;
+                for (int n=k-1;n<k+2;n++) newT[hindex(k-1,j,length)] += T[hindex(k-1,n,length)] * G[n-(k-1)][j-(k-1)];
+                
+            }
+                
+        }
+        
+        // i=k --> m=k, k+1
+        
+        for (int j=k;j<k+3 && j<length;j++) {
+            
+                                                  newT[hindex(k,j,length)]  = 0;
+            for (int n=k;n<k+3 && n<length;n++)   newT[hindex(k,j,length)] +=  c * T[hindex(k,n,length)]   * G[n-(k-1)][j-(k-1)];
+                                                  newT[hindex(k,j,length)] += -s * T[hindex(k,k+1,length)] * G[1][j-(k-1)];
+            for (int n=k+1;n<k+3 && n<length;n++) newT[hindex(k,j,length)] += -s * T[hindex(k+1,n,length)] * G[n-(k-1)][j-(k-1)];
+            
+        }
+        
+        // i=k+1 --> m=k, k+1
+        
+        for (int j=k+1;j<k+3 && j<length;j++) {
+            
+                                                  newT[hindex(k+1,j,length)]  = 0;
+            for (int n=k;n<k+3 && n<length;n++)   newT[hindex(k+1,j,length)] += s * T[hindex(k,n,length)]   * G[n-(k-1)][j-(k-1)];
+                                                  newT[hindex(k+1,j,length)] += c * T[hindex(k,k+1,length)] * G[1][j-(k-1)];
+            for (int n=k+1;n<k+3 && n<length;n++) newT[hindex(k+1,j,length)] += c * T[hindex(k+1,n,length)] * G[n-(k-1)][j-(k-1)];
+            
+        }
+        
+        // Copy values to T
+        
+        if (k>0) { for (int j=k-1;j<k+2;j++)  T[hindex(k-1,j,length)] = newT[hindex(k-1,j,length)]; }
+        for (int j=k;  j<k+3 && j<length;j++) T[hindex(k,j,length)]   = newT[hindex(k,j,length)];
+        for (int j=k+1;j<k+3 && j<length;j++) T[hindex(k+1,j,length)] = newT[hindex(k+1,j,length)];
+        
+        // Update Z
+        
+        for (int i=0;i<length;i++) {
+            
+            double z1 = Z[i][k];
+            double z2 = Z[i][k+1];
+            
+            Z[i][k]   = c * z1 - s * z2;
+            Z[i][k+1] = s * z1 + c * z2;
+            
+        }
+        
+        // Adjust x and z
+        
+        if (k<length-2) {
+            
+            x = T[hindex(k,k+1,length)];
+            z = T[hindex(k,k+2,length)];
+            
+        }
+        
+    }
+    
+}
+
+
+// Given a symmetric matrix A, this algorithm computes its eigenvalues and eigenvectors
+// via the (symmetric) QR algorithm. See Matrix Computations (Golub) for details. This function
+// assumes eigenvectors is preset to the identity matrix before beginning the algorithm.
+
+void symmetricQR(double A[], int length, double eigenvalues[], std::vector<std::vector<double> > &eigenvectors) {
+    
+    double eps = pow(10,-8);    // error tolerance
+    
+    // Tridiagonalize the matrix
+    
+    tridiagonalize(A, length, eigenvectors);
+    
+    // Determine the unreduced subset
+    
+    int q=0;
+    
+    for (int i=length-2;i>=0 && fabs(A[hindex(i,i+1,length)])==0;i--) q++;
+    
+    int p=length-q-2;
+    
+    for (int i=p-1;i>=0 && fabs(A[hindex(i,i+1,length)])!=0;i--) p--;
+    
+    
+    // Loop until convergence
+    
+    //DEBUG
+    int loop=0;
+    
+    while (q<length-1) {
+        
+        //DEBUG
+        loop++;
+        
+        // Apply the QR step to the unreduced part of the matrix between p and length-q-1 (inclusive)
+        
+        int n=length-q-p;
+        
+        double d = (A[hindex(p+n-2,p+n-2,length)] - A[hindex(p+n-1,p+n-1,length)]) / 2;
+        double mu;
+        
+        if (d==0) mu = A[hindex(p+n-1,p+n-1,length)] - fabs(A[hindex(p+n-2,p+n-1,length)]);
+        else      mu = A[hindex(p+n-1,p+n-1,length)] - pow(A[hindex(p+n-2,p+n-1,length)],2)
+                       / (d + sign(d) * sqrt(pow(d,2) + pow(A[hindex(p+n-2,p+n-1,length)],2)));
+        
+        double x = A[hindex(p,p,length)] - mu;
+        double z = A[hindex(p,p+1,length)];
+        
+        for (int k=0;k<n-1;k++) {
+            
+            double c, s;
+            givens(x, z, c, s);
+            
+            double G[4][4] = {{1, 0, 0, 0}, {0, c, s, 0}, {0, -s, c, 0}, {0, 0, 0, 1}};
+            double newA[(n*(n+1))/2];
+            
+            // Multiply A by G on the right and G^T on the left
+            
+            // Start with i=k-1
+            
+            if (k>0) {
+                
+                for (int j=k-1;j<k+2;j++) {
+                    
+                    newA[hindex(k-1,j,n)]                            = 0;
+                    for (int m=k-1;m<k+2;m++) newA[hindex(k-1,j,n)] += A[hindex(p+k-1,p+m,length)] * G[m-(k-1)][j-(k-1)];
+                    
+                }
+                
+            }
+            
+            // i=k
+            
+            for (int j=k;j<k+3 && j<n;j++) {
+                
+                newA[hindex(k,j,n)]                                   = 0;
+                for (int m=k;m<k+3 && m<n;m++)   newA[hindex(k,j,n)] +=  c * A[hindex(p+k,p+m,length)]   * G[m-(k-1)][j-(k-1)];
+                newA[hindex(k,j,n)]                                  += -s * A[hindex(p+k,p+k+1,length)] * G[1][j-(k-1)];
+                for (int m=k+1;m<k+3 && m<n;m++) newA[hindex(k,j,n)] += -s * A[hindex(p+k+1,p+m,length)] * G[m-(k-1)][j-(k-1)];
+                
+            }
+            
+            // i=k+1
+            
+            for (int j=k+1;j<k+3 && j<n;j++) {
+                
+                newA[hindex(k+1,j,n)]                                   = 0;
+                for (int m=k;m<k+3 && m<n;m++)   newA[hindex(k+1,j,n)] += s * A[hindex(p+k,p+m,length)]   * G[m-(k-1)][j-(k-1)];
+                newA[hindex(k+1,j,n)]                                  += c * A[hindex(p+k,p+k+1,length)] * G[1][j-(k-1)];
+                for (int m=k+1;m<k+3 && m<n;m++) newA[hindex(k+1,j,n)] += c * A[hindex(p+k+1,p+m,length)] * G[m-(k-1)][j-(k-1)];
+                
+            }
+            
+            // Copy values to A
+            
+            if (k>0) { for (int j=k-1;j<k+2;j++) A[hindex(p+k-1,p+j,length)] = newA[hindex(k-1,j,n)]; }
+            for (int j=k;  j<k+3 && j<n;j++)     A[hindex(p+k,p+j,length)]   = newA[hindex(k,j,n)];
+            for (int j=k+1;j<k+3 && j<n;j++)     A[hindex(p+k+1,p+j,length)] = newA[hindex(k+1,j,n)];
+            
+            // Update matrix of eigenvectors
+            
+            for (int i=0;i<length;i++) {
+                
+                double z1 = eigenvectors[i][p+k];
+                double z2 = eigenvectors[i][p+k+1];
+                
+                eigenvectors[i][p+k]   = c * z1 - s * z2;
+                eigenvectors[i][p+k+1] = s * z1 + c * z2;
+                
+            }
+            
+            // Adjust x and z
+            
+            if (k<(n-1)) {
+                
+                x = A[hindex(p+k,p+k+1,length)];
+                z = A[hindex(p+k,p+k+2,length)];
+                
+            }
+            
+        }
+        
+        // Determine the remaining unreduced subset
+        
+        q=0;
+        
+        for (int i=length-2;i>=0 && fabs(A[hindex(i,i+1,length)]) <= eps * (fabs(A[hindex(i,i,length)]) + fabs(A[hindex(i,i+1,length)])); i--) {
+        
+            q++;
+            A[hindex(i,i+1,length)]=0;
+            
+        }
+        
+        p=length-q-2;
+        
+        for (int i=p-1;i>=0 && fabs(A[hindex(i,i+1,length)])!=0;i--) {
+            
+            p--;
+            
+            if (fabs(A[hindex(i,i+1,length)]) <= eps * (fabs(A[hindex(i,i,length)]) + fabs(A[hindex(i,i+1,length)]))) {
+                
+                A[hindex(i,i+1,length)]=0;
+                p++;
+                break;
+                
+            }
+            
+        }
+        
+    }
+    
+    for (int i=0;i<length;i++) eigenvalues[i]=A[hindex(i,i,length)];
+    
+}
diff --git a/tools.h b/tools.h
new file mode 100644
index 0000000..67c6b13
--- /dev/null
+++ b/tools.h
@@ -0,0 +1,263 @@
+#ifndef TOOLS_H
+#define TOOLS_H
+
+
+#include <string>
+#include <cstring>
+#include <iostream>
+#include <sstream>
+#include <math.h>
+#include <stdio.h>
+
+
+// Typedefs
+typedef std::vector<unsigned long> Key;
+typedef std::vector<std::vector<double> > Vector;
+typedef std::vector<std::vector<int> > IntVector;
+
+
+double computeGamma_L0(const Vector &, double);
+double computeGamma_L2(const Vector &, double);
+
+double innerProduct(double[], double[], size_t);
+double innerProduct(const std::vector<double> &, const std::vector<double> &);
+double innerProduct(const Vector &, const Vector &);
+double innerProduct(const Vector &, const Vector &, const IntVector &);
+double L1(double[], size_t);
+double L1(const std::vector<double> &);
+double L1(const Vector &);
+double L2(double[], size_t);
+double L2(const std::vector<double> &);
+double L2(const Vector &);
+double L2(const Vector &, const IntVector &);
+double LInfinity(double[], size_t);
+double LInfinity(const std::vector<double> &);
+double LInfinity(const Vector &);
+double LInfinity(const Vector &, const IntVector &);
+
+void insertInPlace(std::vector<int> &, int);
+void insertInPlace(int[], int, int);
+void eraseInPlace(std::vector<int> &, int);
+void eraseInPlace(int[], int, int);
+void quickSort(double[][2], int, int);
+void quickSort(double[][3], int, int);
+void quickSort(std::vector<int> &, int, int);
+bool checkSize(const Key &, const Key &);
+void getSuperset(const Key &, const Key &, Key &);
+
+void unpack(Vector &, std::vector<double> &, const Vector &);
+void pack(Vector &, std::vector<double> &, const Vector &);
+void unpack_sparse(Vector &, std::vector<double> &, const Vector &, const IntVector &);
+void pack_sparse(Vector &, std::vector<double> &, const Vector &, const IntVector &);
+
+void modifiedCholeskyInverse(double[], size_t);
+void modifiedCholeskyInverse(Vector &, std::vector<double> &, int, const Vector &);
+void modifiedCholeskyInverse(std::vector<double> &, std::vector<double> &, std::vector<double> &, std::vector<double> &, size_t);
+void modifiedCholeskyInverse_sparse(Vector &, std::vector<double> &, int, const Vector &, const IntVector &);
+void modifiedCholeskyInverse_sparse(Vector &, std::vector<double> &, std::vector<double> &, std::vector<double> &, std::vector<double> &, const Vector &, const IntVector &);
+void modifiedCholeskyInverse(std::vector<double> &, int);
+double modifiedCholeskyInverseDet(double[], size_t);
+double modifiedCholeskyInverseDet(Vector &, int, const Vector &);
+double symmetricDet(double[], int);
+
+void givens(double, double, double &, double &);
+void householder(double[], int, double[], double &);
+void tridiagonalize(double[], int, std::vector<std::vector<double> > &);
+void QRStep(double[], int, std::vector<std::vector<double> > &);
+void symmetricQR(double[], int, double[], std::vector<std::vector<double> > &);
+
+
+
+
+
+// STRING MANIPULATION
+
+// Converts generic to string
+
+template <class T>
+
+inline std::string tostring (const T& t) {
+    
+    std::stringstream ss;
+    ss << t;
+    return ss.str();
+    
+}
+
+
+// Converts a string to an integer
+
+inline int strtoint(const std::string &s) {
+    
+    std::istringstream i(s);
+    int x;
+    
+    if (!(i >> x)) printf("String to integer conversion failed!");
+    return x;
+    
+}
+
+
+// Converts a string to a double
+
+inline double strtodouble(const std::string &s) {
+    
+    std::istringstream i(s);
+    double x;
+    
+    if (!(i >> x)) printf("String to double conversion failed!");
+    return x;
+
+}
+
+
+// MISC FUNCTIONS
+
+// Given the size of a coupling or correlation vector, return the corresponding number of spins
+
+inline int sizetolength(size_t size) {
+    
+    return (int) ((sqrt(1 + 8 * size) - 1) / 2);
+    
+}
+
+
+// Checks to see whether the ith bit of an integer in a key is occupied
+
+inline int checkBit(unsigned long &key, int i) {
+    
+    return (int) (key & 1<<i);
+    
+}
+
+
+// Flip the ith bit of an integer
+
+inline void flipBit(unsigned long &key, int i) {
+    
+    key ^= 1<<i;
+    
+}
+
+
+// For j>i, the pair {i,j} in the list {{0,1}, {0,2}, ... } is located at offset(i,length)+j
+
+inline int offset(int i, size_t length) {
+    
+    return (int) (length + i * (length - 2) - (i * (i - 1)) / 2 - 1);
+    
+}
+
+
+// For j>i, the pair {i,j} in the list {{0,0}, {0,1}, ... } is located at hindex(i,j,length)
+// This routine is necessary when the diagonal term is also included
+
+inline int hindex(int i, int j, size_t length) {
+    
+    return (int) (i * length - (i * (i + 1)) / 2 + j);
+    
+}
+
+
+// Returns as above, but with a check to make sure that j>i
+
+inline int safeHindex(int i, int j, size_t length) {
+    
+    if (j>i) return (int) (i * length - (i * (i + 1)) / 2 + j);
+    else     return (int) (j * length - (j * (j + 1)) / 2 + i);
+    
+}
+
+
+// Return the location of a pair of states in canonical order, {{0,0}, {0,1}, ..., {0,qj}, {1,0}, ...}
+
+inline int sindex(int i, int j, size_t qi, size_t qj) {
+    
+    return (int) (i * qj + j);
+    
+}
+
+
+// Return the location of a triple of states in canonical order
+
+inline int sindex3(int i, int j, int k, size_t qi, size_t qj, size_t qk) {
+    
+    return (int) (i * qj * qk + j * qk + k);
+    
+}
+
+
+// For j>i, the pair {i,j} in the list {{0},{1},...,{length-1},{0,1}, {0,2}, ... } is located at index(i,j,length)
+
+inline int index(int i, int j, size_t length) {
+    
+    return (int)(length + i * (length - 2) - (i * (i - 1)) / 2 - 1 + j);
+    
+}
+
+
+// For k>j>i, return the index of {i,j,k} in the full list
+
+inline int index(int i, int j, int k, int length) {
+    
+    return ( ( (i * (i + 1) * ( i + 2 )) / 3 
+                  - (length-1) * (i * (i + 2) - 2 * j + 1) 
+                  + ((int) pow(length-1, 2)) * (i - 1)
+                  - (j * (j + 1) - 2 * k)                ) / 2 );
+    
+}
+
+
+// For l>k>j>i, return the index of {i,j,k,l} in the full list
+
+inline int index(int i, int j, int k, int l, int length) {
+    
+    return ( (-i * (i + 1) * (i + 2) * (i + 3) 
+              + 4 * (j * (1 + j) * (j + 2) - 3 * k * (k + 1) + 6 * l - 5 * (length-1)) 
+              + 2 * (i * (11 + i * (9 + 2 * i)) - 6 * j * (j + 2) + 12 * k) * (length-1) 
+              - 6 * (i * (i + 3) - 2 * j) * ((int) pow(length-1, 2))
+              + 4 * (i - 1) * ((int) pow(length-1,3))                                   ) / 24 );
+    
+}
+
+
+// Returns the sign of an integer x (and 0 if x==0)
+
+inline int sign(int x) {
+                
+    static const int lookup_table[] = { -1, 1, 0, 0 };
+    return lookup_table[(x==0) * 2 + (x>0)];
+
+}
+
+
+// Returns the sign of a double x (and 0 if x==0)
+
+inline double sign(double x) {
+    
+    //double eps=pow(10.0,-20.0);
+    
+    static const double lookup_table[] = { -1, 1, 0, 0 };
+    return lookup_table[(x==0) * 2 + (x>0)];
+    
+}
+
+
+// Given components (a, b), this function computes parameters for a Givens rotation matrix G = ((c, s), (-s, c)),
+// such that G^T (a, b) = (r, 0). See Matrix Computations (Golub) for details.
+
+inline void givens(double a, double b, double &c, double &s) {
+ 
+    if (b==0) { c=1; s=0; }
+    
+    else {
+        
+        if (fabs(b) > fabs(a)) { s = 1 / sqrt(1 + pow(a,2) / pow(b,2)); c = -a * s / b; }
+        else                   { c = 1 / sqrt(1 + pow(b,2) / pow(a,2)); s = -b * c / a; }
+        
+    }
+    
+}
+
+
+#endif