main.cpp.gilps

/*
 *  main.cpp
 *  Population_Dynamics
 *
 *  Created by Greg Hart on May 2015
 *  Copyright 2015 ALF. All rights reserved.
 *
 */
/*
 * IN
 * resIdx.dat 		- ordered list of residue types - as integer codes - present at each site from most to least probable
 * h.dat 		- h parameters
 * J.dat 		- J parameters
 * init_seq.dat		- (optional) a starting sequence, if file is not found then the WT is used.
 *
 * OUT
 * P1_model_traj.dat 	- running trajectory of P1 values of the population
 * P2_model_traj.dat 	- running trajectory of P2 values of the population
 * P1_model.dat 	- terminal P1 values computed by the population path
 * P2_model.dat	 	- terminal P2 values computed by the population path
 * MC_seqs.dat 		- sequences sampled from population trajectory
 * pop_stats.dat	- data on the population at different time points
 * Tcell_traj.dat	- Number of Tcells and targets at different time points
 */
 
#include "main.h"
#include "functions.h"

#include <string>
#include <iostream>
#include <fstream>
#include <sstream>
#include <iomanip>
#include <stdexcept>
#include <vector>
#include <limits>
#include <algorithm>
#include <math.h>
#include <string.h>
#include <unistd.h>
#include <stdio.h>
#include <stdlib.h>
//#include <stdint.h>

#ifdef _OPENMP
#include <omp.h>
#endif

#include <trng/yarn2.hpp>
#include <trng/uniform01_dist.hpp>
#include <trng/poisson_dist.hpp>
#include <boost/numeric/odeint.hpp>

#define PARTPARALLEL 3

using std::cout;
using std::cin;
using std::cerr;
using std::endl;
using std::vector;
using std::sort;
using namespace boost::numeric::odeint;


typedef std::vector<double> state_type;

// main function preforming Fisher-wright dynamics on the viral side and intergrating ODEs for the Tcell side
int main (int argc, char *argv[]) {
	
    cout.setf(std::ios_base::scientific);
    cout.precision(3);	
	
    #ifdef UNIT_TESTING
	
    #endif  // #ifdef UNIT_TESTING

	
    // start timer
    double wall_start = get_wall_time();
    double cpu_start  = get_cpu_time();
	
	
    // inputs with default values
    long n_inputs=15;
	
    long seed = -240164;
    long m=3;
    long N=10000;
    double rate=0.0001;
    long progeny = 10;
    long n_cycles = 20000;
    long burnin = 2000;
    long sample_mod = 10000;
    double T = 1;
    long print_mod = 1;
    long write_mod = 5;
    long n_epitope = 1;
    double n_T = 1000;
    int rep_lim = 9;
    double T_penalty = 1;
    
    double a = 1e-7;
    double ap = 1e-6;
    double b = 1e-4;
    double d = 0.2;//0.5;
    double dp = 3;
    double re = 6;//4;
    double e = 3*1e-4;
    double g = 0.03;
    double hr = 0.03;
    double w = .1;//7.5;

    std::ifstream fin_inputs;
    fin_inputs.open("./inputs.dat");
    if (!fin_inputs) {
	cerr << "Cannot open input file inputs.dat in the current directory; aborting." << endl;
	exit(-1);
    }
    long input_cntr=0;
    while (!fin_inputs.eof()) {
			
	std::string input_name, tmp_str;
	fin_inputs >> input_name;
			
	if (strcmp(input_name.c_str(),"seed")==0) {
	    fin_inputs >> seed;
	    input_cntr++;
	} else if (strcmp(input_name.c_str(),"m")==0) {
	    fin_inputs >> m;
	    input_cntr++;
	} else if (strcmp(input_name.c_str(),"N")==0) {
	    fin_inputs >> N;
	    input_cntr++;
	} else if (strcmp(input_name.c_str(),"rate")==0) {
	    fin_inputs >> rate;
	    input_cntr++;
        } else if (strcmp(input_name.c_str(),"progeny")==0) {
            fin_inputs >> progeny;
            input_cntr++;
	} else if (strcmp(input_name.c_str(),"n_cycles")==0) {
	    fin_inputs >> n_cycles;
	    input_cntr++;
	} else if (strcmp(input_name.c_str(),"burnin")==0) {
	    fin_inputs >> burnin;
	    input_cntr++;
	} else if (strcmp(input_name.c_str(),"sample_mod")==0) {
	    fin_inputs >> sample_mod;
	    input_cntr++;
	} else if (strcmp(input_name.c_str(),"T")==0) {
	    fin_inputs >> T;
	    input_cntr++;
	} else if (strcmp(input_name.c_str(),"print_mod")==0) {
	    fin_inputs >> print_mod;
	    input_cntr++;
	} else if (strcmp(input_name.c_str(),"write_mod")==0) {
	    fin_inputs >> write_mod;
	    input_cntr++;
    } else if (strcmp(input_name.c_str(),"n_epitope")==0) {
            fin_inputs >> n_epitope;
            input_cntr++;
    } else if (strcmp(input_name.c_str(),"n_T")==0) {
        fin_inputs >> n_T;
        input_cntr++;
    } else if (strcmp(input_name.c_str(),"rep_lim")==0) {
        fin_inputs >> rep_lim;
        input_cntr++;
    } else if (strcmp(input_name.c_str(),"T_penalty")==0) {
            fin_inputs >> T_penalty;
            input_cntr++;
	}
	getline(fin_inputs,tmp_str);
    }
	
    if (input_cntr!=n_inputs) {
	cerr << "ERROR - did not read expected number of inputs from inputs.dat; aborting." << endl;
	exit(-1);
    }
	
    fin_inputs.close();
	
    // loading parameters
    //- loading ordered residue indices at each site (resIdx.dat), computing total number of h and J parameters, 
    //  and candidateSites list of sites containing more than one residue type as only sites at which we attempt mutations
    std::vector<long> nRes(m,-1);
    std::vector< std::vector<long> > resIdx(m);
    load_resIdx("./resIdx.dat",nRes,resIdx,m);
    long n_h=0;
    long n_J=0;
    for(long i=0; i<m; i++) {
	n_h += nRes[i];
        for(long j=i+1; j<m; j++) {
		n_J += nRes[i]*nRes[j];
	}
    }
   
    // 6 hours to replicate, 7 day half life, 40 replication complexes per cell.
//    long progeny = 4*7*40*(n_h-m)/(19*m);
    cout << "progeny = " << progeny << endl;
 
    std::vector<long> epi_start(n_epitope,-1);
    std::vector<long> epi_end(n_epitope,-1);
    std::vector< std::vector<double> > chi(n_epitope);
    //setup T cell populations
    std::vector<state_type> Tcells(n_epitope, state_type(rep_lim+2));
    double chiI[n_epitope];
    std::vector< std::vector<int> > place_value(n_epitope);
    std::vector< std::vector<int> > n_WTepitopes(n_epitope);
    if(n_epitope > 0){
        load_epitopes("./epitopes.dat", epi_start, epi_end, chi, n_epitope, resIdx, nRes, Tcells, rep_lim);
    }
    for (long i=0; i<n_epitope; i++){
        (n_WTepitopes[i]).resize((chi[i]).size(), 0);
        // create a way to index all mutations in an epitope
        (place_value[i]).resize(epi_end[i]-epi_start[i], 0);
        int size = 1;
        int count = 0;
        for(int j=epi_start[i]; j<epi_end[i]; j++){
            place_value[i][count] = size;
            size *= nRes[j];
            count++;
        }
    }
    
    cout << "\n";
    cout << "$$ # parameters: $$" << endl;
    cout << "  # h parameters = " << n_h << endl;
    cout << "  # J parameters = " << n_J << endl;
    cout << "  # epitopes = " << n_epitope << endl;
    cout << endl;
	
    std::vector <long> candidateSites(0);
    for(long i=0; i<m; i++) {
	if (nRes[i]>1) {
	    candidateSites.push_back(i);
	}
    }
    long n_candidateSites = candidateSites.size();
	
		
    //- loading h parameters
    std::vector< std::vector<double> > h(m);
    for(long i=0; i<m; i++) {
	(h[i]).resize(nRes[i],(double)0);
    }
    load_X1("./h.dat",nRes,h,m);
	
	
    //- loading J parameters
    // \-> loading only upper triangle 
    // \-> loading residuewise elements in column major order (i.e., over i residues first)
    std::vector< std::vector< std::vector< std::vector<double> > > > J(m);
    for(long i=0; i<m; i++) {
	(J[i]).resize(m);
        for(long j=i+1; j<m; j++) {
	    (J[i][j]).resize(nRes[i]);
            for(long p=0; p<nRes[i]; p++) {
		(J[i][j][p]).resize(nRes[j],(double)0);
	    }
	}
    }
    load_X2("./J.dat",nRes,J,m);
	
	
    // initializing trajectory files
    FILE* fout_P1_traj;
    {
        std::string fstr_P1_traj="./P1_model_traj.dat";
        fout_P1_traj = fopen(fstr_P1_traj.c_str(),"wb");
        if (fout_P1_traj==NULL) {
	    cerr << "Cannot open output file " << fstr_P1_traj << " in the current directory; aborting." << endl;
	    exit(-1);
        }
        int32_t sizes [2];
        sizes[0] = sizeof(long);
        sizes[1] = sizeof(double);
        fwrite(sizes,4,2,fout_P1_traj);    
    } 
	
    FILE* fout_P2_traj;
    {
        std::string fstr_P2_traj="./P2_model_traj.dat";
        fout_P2_traj = fopen(fstr_P2_traj.c_str(),"wb");
        if (fout_P2_traj==NULL) {
	    cerr << "Cannot open output file " << fstr_P2_traj << " in the current directory; aborting." << endl;
	    exit(-1);
	}
        int32_t sizes [2];
        sizes[0] = sizeof(long);
        sizes[1] = sizeof(double);
        fwrite(sizes,4,2,fout_P2_traj);    
   }

   FILE* fout_MC_seqs;
   {
        std::string fstr_MC_seqs="./MC_seqs.dat";
        fout_MC_seqs = fopen(fstr_MC_seqs.c_str(),"wb");
        if (fout_MC_seqs==NULL) {
	    cerr << "Cannot open output file " << fstr_MC_seqs << " in the current directory; aborting." << endl;
	    exit(-1);
        }
        int32_t sizes [2];
        sizes[0] = sizeof(long);
        sizes[1] = sizeof(int8_t);
        fwrite(sizes,4,2,fout_MC_seqs);
        fwrite(&N,sizeof(long),1,fout_MC_seqs); 
    }
    
    
    FILE* fout_pop_stats;
    {
        std::string fstr_pop_stats="./pop_stats.dat";
        fout_pop_stats = fopen(fstr_pop_stats.c_str(),"wb");
        if (fout_pop_stats==NULL) {
            cerr << "Cannot open output file " << fstr_pop_stats << " in the current directory; aborting." << endl;
            exit(-1);
        }
        int32_t sizes [2];
        sizes[0] = sizeof(long);
        sizes[1] = sizeof(double);
        fwrite(sizes,4,2,fout_pop_stats);    
    }
   
 
    FILE* fout_fitness;
    {
        std::string fstr_fitness="./fitness.dat";
        fout_fitness = fopen(fstr_fitness.c_str(),"wb");
        if (fout_fitness==NULL) {
            cerr << "Cannot open output file " << fstr_fitness << " in the current directory; aborting." << endl;
            exit(-1);
        }
        int32_t sizes [2];
        sizes[0] = sizeof(long);
        sizes[1] = sizeof(double);
        fwrite(sizes,4,2,fout_fitness);    
        fwrite(&N,sizeof(long),1,fout_fitness); 
    }
    

    FILE* fout_Tcell_traj;
    {
        std::string fstr_Tcell_traj="./Tcell_traj.dat";
        fout_Tcell_traj = fopen(fstr_Tcell_traj.c_str(),"wb");
        if (fout_Tcell_traj==NULL) {
            cerr << "Cannot open output file " << fstr_Tcell_traj << " in the current directory; aborting." << endl;
            exit(-1);
        }
        int32_t sizes [3];
        sizes[0] = sizeof(long);
        sizes[1] = sizeof(double);
        sizes[2] = n_epitope;
        fwrite(sizes,4,3,fout_Tcell_traj);    
    }

    // Out put files for epitope seqs
    FILE* fout_epi_traj[n_epitope];
    if(n_epitope > 0){
        std::string fstr = "./epitopes.dat";
        std::ifstream fin;
        fin.open(fstr.c_str());
        if (!fin) {
            std::cerr << "Cannot open input file " << fstr << "; aborting." << std::endl;
            exit(-1);
        }
    
        std::string epi_file;
        for(long i=0; i<n_epitope; i++) {
            if (fin.eof()) {
                std::cerr << "Prematurely encountered eof in " << fstr << "; aborting." << std::endl;
                exit(-1);
            }
        
            std::getline(fin, epi_file);
            std::stringstream readStrStrm(epi_file);
            readStrStrm >> epi_file;
            size_t pos = 0;
            pos = epi_file.find_last_of('/');
            epi_file = epi_file.substr(pos+1);
            std::string token;
            pos = epi_file.find('_');
            epi_file.erase(0, pos + 1);
            pos = epi_file.find('_');
            token = epi_file.substr(0, pos);
            pos = token.find('-');
            {
                fout_epi_traj[i] = fopen(epi_file.c_str(),"wb");
                if (fout_epi_traj==NULL) {
                    cerr << "Cannot open output file " << epi_file << " in the current directory; aborting." << endl;
                    exit(-1);
                }
                int32_t sizes [2];
                sizes[0] = sizeof(long);
                sizes[1] = sizeof(int8_t);
                fwrite(sizes,4,2,fout_epi_traj[i]);    
                fwrite(&N,sizeof(long),1,fout_epi_traj[i]); 
            } 
        }
    
        while(!fin.eof()) {
            std::getline(fin,epi_file);
            if(!epi_file.empty()){
                 std::cerr << "Elements remain in " << fstr << " after populating array; aborting." << std::endl;
                 exit(-1);
            }
        }
    }
	
    // initializing model probability arrays
    std::vector< std::vector<double> > n1(h);
    std::vector< std::vector< std::vector< std::vector<double> > > > n2(J);
    for(long i=0; i<m; i++) {
        for(long p=0; p<nRes[i]; p++) {
	    n1[i][p] = 0;
            for(long j=i+1; j<m; j++) {
                for(long q=0; q<nRes[j]; q++) {
		    n2[i][j][p][q] = 0;
		}
	    }
	}
    }
    
    printf("PARTPARALLEL = %d\n", PARTPARALLEL);
	
    // Population dynamics over empirical fitness landscape
    cout << "Commencing evolutionary dynamics on the Potts model..." << endl;
	
    // initialize population arrays
    std::vector< std::vector<int8_t> > population(N, std::vector<int8_t>(m));
    std::vector< std::vector<int8_t> > temp_pop(progeny*N, std::vector<int8_t>(m));
    double* boltzfac = new double[N];
    double* temp_boltz = new double[N*progeny];
    double* part_sum = new double[N*progeny];
    bool* part_sum_used = new bool[N*progeny];
    double Z = 1.0;
    double Ztemp = 1.0;
    double* boltzfac_eff = new double[N];
    double Z_eff = 1.0;
    double* E = new double[N];
    double* E_eff = new double[N];
    double E_ave = 0.0;
    double E_eff_ave = 0.0;
    int* pop_init = new int[m];
    bool* occupied = new bool[N*progeny];
    int* idx = new int[N*progeny];
	
    // Initialize every sequence in the population to the sequence found in the file or the WT if no file.
    load_seq("init_seq.dat", pop_init, nRes, resIdx, m);
    for(long i=0; i<N; i++){
        for(long j=0; j<m; j++){
            population[i][j]=pop_init[j];
	}
    }
    delete[] pop_init;
	
    // initialize occupied arrays
    for(long i = 0; i < N*progeny; i++){
        occupied[i] = false;
    }

    double temp = 0;
    for(long i=0; i<m; i++) {
        for(long p=0; p<nRes[i]; p++) {
            temp += h[i][p];
        }
    }
    T_penalty = T_penalty*temp/n_h;
    cout << "T penalty: " << T_penalty << endl;

    // Calulate the boltzmann factor for each sequence and the partition function as if only the sequences in the population are possible.
    E[0] = energy(population[0],m,h,J);
    E_eff[0] = E[0];
    E_ave = E[0];
    E_eff_ave = E[0];
    boltzfac[0] = exp(-E[0]/T);
    boltzfac_eff[0] = boltzfac[0];
    for(long i=1; i<N; i++){
        E[i] = E[0];
        E_eff[i] = E[0];
        boltzfac[i]=boltzfac[0];
        boltzfac_eff[i] = boltzfac[0];
    }
    Z = boltzfac[0];
    Ztemp = Z;
    Z_eff = boltzfac_eff[0];
    
    // calculate initial refrencies of amino acids
    for(long k=0; k<N; k++){
        for(long i=0; i<m; i++) {
             n1[i][population[k][i]] += 1.0;
             for(long j=i+1; j<m; j++) {
                n2[i][j][population[k][i]][population[k][j]] += 1.0;
             }
        }
    }
    
    // write initial state to file
    long zero = 0;
    fwrite(&zero, sizeof(long), 1, fout_pop_stats);
    fwrite(&N, sizeof(long), 1, fout_pop_stats);
    fwrite(&Ztemp, sizeof(double), 1, fout_pop_stats);
    fwrite(&Z, sizeof(double), 1, fout_pop_stats);
    fwrite(&Z_eff, sizeof(double), 1, fout_pop_stats);
    fwrite(&E_ave, sizeof(double), 1, fout_pop_stats);
    fwrite(&E_eff_ave, sizeof(double), 1, fout_pop_stats);

    fwrite(&zero, sizeof(long), 1, fout_P1_traj);
    for(long i=0; i<m; i++) {
        for(long p=0; p<nRes[i]; p++) {
            n1[i][p] = n1[i][p]/N;
            fwrite(&n1[i][p], sizeof(double), 1, fout_P1_traj);
            n1[i][p] = n1[i][p]*N;
        }
    }
                 
    fwrite(&zero, sizeof(long), 1, fout_P2_traj);
    for(long i=0; i<m; i++) {
        for(long j=i+1; j<m; j++) {
            for(long q=0; q<nRes[j]; q++) {
                for(long p=0; p<nRes[i]; p++) {
                    n2[i][j][p][q] = n2[i][j][p][q]/N;
                    fwrite(&n2[i][j][p][q], sizeof(double), 1, fout_P2_traj);
                    n2[i][j][p][q] = n2[i][j][p][q]*N;
                }
            }
        }
    }

    // reset amino acid frequencies
    for(long i=0; i<m; i++) {
        for(long p=0; p<nRes[i]; p++) {
            n1[i][p] = 0.0;
            for(long j=i+1; j<m; j++) {
                for(long q=0; q<nRes[j]; q++) {
                    n2[i][j][p][q] = 0.0;
                }
            }
        }
    }

    fwrite(&zero, sizeof(long), 1, fout_Tcell_traj);
    for(long i=0; i<n_epitope; i++){
        double Etotal = 0.0;
        fwrite(&Tcells[i][0], sizeof(double), 1, fout_Tcell_traj);
        for(long j=1; j<=rep_lim; j++){
            Etotal += Tcells[i][j];
        }
        fwrite(&Etotal, sizeof(double), 1, fout_Tcell_traj);
        fwrite(&Tcells[i][rep_lim+1], sizeof(double), 1, fout_Tcell_traj);
        int index = 0;
        int count = 0;
        for(long position=epi_start[i]; position<epi_end[i]; position++){
            index += place_value[i][count]*population[0][position];
            count++;
        }
        chi[i][index] = chi[i][index]*N;
        fwrite(&chi[i][index], sizeof(double), 1, fout_Tcell_traj);
        chi[i][index] = chi[i][index]/N;
    }
                
    fwrite(&zero, sizeof(long), 1, fout_MC_seqs);
    for(long k=0; k<N; k++){
        for(long i=0; i<m; i++) {
            fwrite(&resIdx[i][population[k][i]], sizeof(int8_t), 1, fout_MC_seqs);
        }
    }

    for(long j=0; j<n_epitope; j++){
        fwrite(&zero, sizeof(long), 1, fout_epi_traj[j]);
        for(long k=0; k<N; k++){
            for(long i=epi_start[j]; i<epi_end[j]; i++) {
                fwrite(&resIdx[i][population[k][i]], sizeof(int8_t), 1, fout_epi_traj[j]);
            }
        }
    }
	
    fwrite(&zero, sizeof(long), 1, fout_fitness);
    fwrite(E, sizeof(double), N, fout_fitness);
    fwrite(E_eff, sizeof(double), N, fout_fitness);
    fwrite(boltzfac, sizeof(double), N, fout_fitness);
    fwrite(boltzfac_eff, sizeof(double), N, fout_fitness);
	
	// N.B. For execution acceleration, this subroutine operates under fake residue numbering scheme (fakeres) whereby 
	//      most probable residue at site i is coded as 0, next most probable as 1, ..., least probable as nRes[i]-1
	//      to obviate frequent (slow) lookup of residue types in resIdx
	//
	//		Under this scheme, the fakeres code is identical to its position in resIdx, facilitating rapid lookup of 
	//		corresponding n1,n2,h,J matrix elements
	//
	
    long pos = 0;
    long n_samples = 0;
    double Etemp[rep_lim];
    // \-> safeUnity guards against blue moon segmentation faults for ran2 returning precisely unity due to interaction with floor call (ran2 precisely zero is not problematic)
    static double safeUnity = 1.0 - std::numeric_limits<double>::epsilon();
    double Energy = 0.0;
    double penalty_S;
    int thread_num;
    #pragma omp parallel shared(thread_num)
    {
	if(omp_get_thread_num()==0){
	        thread_num = omp_get_num_threads();
        }
    }
    int seq_S;
    #if PARTPARALLEL == 2 || PARTPARALLEL == 1 // unused in 1,0
    int *seq_arr = new int[thread_num];
    #elif PARTPARALLEL == 3 || PARTPARALLEL == 0
    int *seq_arr = new int[N];
    #endif
    vector<vector <int>> energy_idx; // This matrix of index is for parallelizing
    // the energy calculation in PART 2
    init_energy_index(m, thread_num, energy_idx);
    double *key = new double[N*progeny];
    double *key_l = new double[N*progeny];
		
    // If omp is in use parallelize all operations on the population for speed.
    #pragma omp parallel default(shared) shared(stdout, part_sum_used, seq_arr, energy_idx, Energy, seq_S, penalty_S, nRes, resIdx, h, J, T, population, temp_pop, temp_boltz, m, N, boltzfac, boltzfac_eff, Z, Z_eff, Ztemp, fout_P1_traj, fout_P2_traj, fout_MC_seqs, fout_pop_stats, fout_Tcell_traj, fout_fitness, fout_epi_traj, write_mod, print_mod, sample_mod, n_samples, cout, cerr, n_cycles, n1, n2, burnin, progeny, rate, seed, safeUnity, pos, Tcells, Etemp, n_T, chi, chiI, T_penalty, n_WTepitopes, n_epitope, epi_start, epi_end, rep_lim, occupied, idx, part_sum, place_value, E, E_ave, E_eff, E_eff_ave)
    {
    
    // Parallel safe PRNG
    trng::yarn2 r((long)fabs(seed));        // PRNG
    trng::uniform01_dist<> u;               // uniform distribution for probabilities
    trng::poisson_dist IntGenerator(3*m*rate);  // poisson distribution for number of mutants in a sequence

    #ifdef _OPENMP
    const int size=omp_get_num_threads();    // get total number of processes
    const int rank=omp_get_thread_num();     // get rank of current process
    r.split(size, rank);               // split PRN sequences by leapfrog method - choose sub-stream no. rank out of size streams
    #else
    const int size = 1;
    const int rank = 0;
    #endif
	
    if(rank == 0) {
        printf("OPENMP: %d threads\n", size);
        fflush(stdout);
    }

    // break up the population into a chunk for each processor.
    const long len = N/size;
    const long begin = rank*len;
    const long end = rank + 1 == size ? N : begin+len;
	
    // break up DNA in later computation
    const long len_m = m/size;
    const long begin_m = rank*len_m;
    const long end_m = rank + 1 == size ? m : begin_m+len_m;
    
    // start main loop
    for (long cycle=1; cycle<=n_cycles; cycle++) {
        
	// reset values used in every iteration of the loop
        #pragma omp single
        {
            pos = 0;
            Ztemp = 0.0;
            E_ave = 0.0;
            E_eff_ave = 0.0;
            for(long i=0; i<n_epitope; i++){
               for(long j=0; j<(n_WTepitopes[i]).size(); j++){
                   n_WTepitopes[i][j] = 0;
               }
	    }
        }
        #pragma omp barrier

        int ranInt1 = 0;	// IntGenerator(r)=number of mution for the parent sequence;
        double ranNum1 = 0;     //=u(r)*safeUnity;	// range [0,1)
        long parent = 0;
        long site = 0;
        long res = 0;
        double Z_local = 0;
        double Ztemp_local = 0;
        double Z_local_eff = 0;
        
		
        // Have every sequence in the population produce offspring equal to the value of progeny. Mutating each offspring sequence as it is copied
        for(long k=0; k<progeny; k++){       // loop over the number of progeny
            for(long i=begin; i<end; i++){   // loop over the processor's part of the population
                ranInt1 = IntGenerator(r);   // number of mution for the parent sequence;
                for(long j=0; j<m; j++){     // copy parent sequence
                    temp_pop[i+(N*k)][j] = population[i][j];
                }
                for(long j=0; j<ranInt1; j++){ // introduce mutations
                    ranNum1 = u(r)*safeUnity;
             	    site = (long)floor((double)ranNum1 * (double)m);
                    ranNum1 = u(r)*safeUnity;
                    res = (long)floor((double)ranNum1 * (double)(nRes[site]-1));
                    if (temp_pop[i+(N*k)][site] <= res) {
     	     		    res++;
          	        }
                    temp_pop[i+(N*k)][site] = res;
                }
				
                idx[i+(N*k)] = i+(N*k);
                double penalty = 0;
                for(long j=0; j<n_epitope; j++){ // calculate T-cell susceptibility
                    int index = 0;
                    int count = 0;
                    for(long position=epi_start[j]; position<epi_end[j]; position++){
                        index += place_value[j][count]*temp_pop[i+(N*k)][position];
                        count++;
                    }
                    double num_Ecells = 0;
                    for(long t=1; t<=rep_lim; t++){
                        num_Ecells += Tcells[j][t];
                    }
                    penalty += T_penalty*chi[j][index]*num_Ecells;
                }
                // calculate boltzmann factor of daughter sequence and add it to the partition function
                temp_boltz[i+(N*k)] = exp(-(energy(temp_pop[i+(N*k)],m,h,J))/T-penalty/T);
                Ztemp_local += temp_boltz[i+(N*k)];
                occupied[i+(N*k)] = true;
            }
        }

        // combine partition function from each processor
        #pragma omp atomic
            Ztemp += Ztemp_local;
	
        #pragma omp barrier
		

#if PARTPARALLEL == 0		
        #pragma omp single
        {
            sort(idx, idx+N*progeny, [&temp_boltz](size_t i1, size_t i2){return temp_boltz[i1] < temp_boltz[i2];});
            part_sum[0] = temp_boltz[idx[0]];
            for(long i=1; i<N*progeny; i++){
                part_sum[i] = part_sum[i-1] + temp_boltz[idx[i]];
            }
	
            // pick sequences that survive to the final population based on their fitness
            Z = 0.0;
            Z_eff = 0.0;
            E_ave = 0.0;
            E_eff_ave = 0.0;
            while(pos < N){
                ranNum1 = u(r)*safeUnity*part_sum[N*progeny-1];
                long seq = std::lower_bound(part_sum, part_sum + N*progeny,ranNum1) - part_sum;
	            seq_arr[pos] = seq;

                if(seq >= N*progeny){
                    cout << "Error in find, seq = " << seq << ", ranNum1 = " << ranNum1 << ", Ztemp = " << Ztemp << endl;
                }
                for(long j=0; j<m; j++){
                    population[pos][j]=temp_pop[idx[seq]][j];
                }
                // prime T cells
                double penalty = 0;
                for(long j=0; j<n_epitope; j++){
                    int index = 0;
                    int count = 0;
                    for(long position=epi_start[j]; position<epi_end[j]; position++){
                        index += place_value[j][count]*population[pos][position];
                        count++;
                    }
                    double num_Ecells = 0;
                    for(long t=1; t<=rep_lim; t++){
                        num_Ecells += Tcells[j][t];
                    } 
                    penalty += T_penalty*chi[j][index]*num_Ecells;
                    n_WTepitopes[j][index]+=1.0;
                }
                E[pos] = energy(population[pos],m,h,J);
                boltzfac[pos] = exp(-(E[pos])/T);
                E_eff[pos] = E[pos]+penalty;
                boltzfac_eff[pos] = exp(-E_eff[pos]/T);
                E_ave += E[pos];
                E_eff_ave += E_eff[pos];
                Z += boltzfac[pos];
                Z_eff += boltzfac_eff[pos];
                pos++;
                occupied[idx[seq]] = false;
                if(seq==0){
                    part_sum[seq] = 0;
                    seq++;
                }
                for(long i=seq; i<N*progeny; i++){
                    if(occupied[idx[i]]){
                        part_sum[i] = temp_boltz[idx[i]] + part_sum[i-1];
                    } else {
                        part_sum[i] = part_sum[i-1];
                    }
                }
            }
		
        } // end #pragma omp single
        #pragma omp barrier
#endif

#if PARTPARALLEL == 1
///////////////////////////////
// PART 2 Introducing Mutation
///////////////////////////////

	#pragma omp single
	{
	    // make sure smaller numbers don't get rounded off
	    sort(idx, idx+N*progeny, [&temp_boltz](size_t i1, size_t i2){return temp_boltz[i1] < temp_boltz[i2];});
            part_sum[0] = temp_boltz[idx[0]];
            for(long i=1; i<N*progeny; i++)
            {
                part_sum[i] = part_sum[i-1] + temp_boltz[idx[i]];
            }
            // pick sequences that survive to the final population based on their fitness
            Z = 0.0;
            Z_eff = 0.0;
            E_ave = 0.0;
            E_eff_ave = 0.0;
	}
	
        while(pos < N){
            if(rank==0){ // only main thread start
		        double ranNum1 = u(r)*safeUnity*part_sum[N*progeny-1];
                long seq = std::lower_bound(part_sum, part_sum + N*progeny,ranNum1) - part_sum;
		        seq_S = seq;
		        if(seq >= N*progeny){
		            cout << "Error in find, seq = " << seq << ", ranNum1 = " << ranNum1 << ", Ztemp = " << Ztemp << endl;
		        }
		for(long j=0; j<m; j++){
		    population[pos][j]=temp_pop[idx[seq]][j];
		}
		// prime T cells
		double penalty = 0;
		for(long j=0; j<n_epitope; j++){
		    int index = 0;
		    int count = 0;
		    for(long position=epi_start[j]; position<epi_end[j]; position++){
		        index += place_value[j][count]*population[pos][position];
			count++;
		    }
		    double num_Ecells = 0;
		    for(long t=1; t<=rep_lim; t++){
			num_Ecells += Tcells[j][t];
		    } 
		    penalty += T_penalty*chi[j][index]*num_Ecells;
		    n_WTepitopes[j][index]+=1.0;
		}
		penalty_S = penalty;
	    } // only mean thread end 
			
	    #pragma omp barrier 
						
	    double E_P = 0.0;
	    auto sequence = population[pos];
	    for(auto i=energy_idx[rank].begin(); i!=energy_idx[rank].end(); i++) {
		    E_P += h[*i][sequence[*i]];
		    for (long j=(*i)+1; j<m; j++){
		        E_P += J[*i][j][sequence[*i]][sequence[j]];
		    }
	    }
	    #pragma omp atomic
		    Energy += E_P;
				
	    #pragma omp barrier	
				
	    if(rank==0)	{	// only main thread start
                E[pos] = Energy;
                boltzfac[pos] = exp(-(E[pos])/T);
                E_eff[pos] = E[pos]+penalty_S;
                boltzfac_eff[pos] = exp(-E_eff[pos]/T);
                E_ave += E[pos];
                E_eff_ave += E_eff[pos];
                Z += boltzfac[pos];
                Z_eff += boltzfac_eff[pos];
                pos++;
                occupied[idx[seq_S]] = false;
                Energy = 0.0;
	    }// only main thread end

	    if(size == 1 || rank != 0){
 	        #pragma omp for schedule(static)
		    for(long i=seq_S; i<N*progeny; i++){
		        part_sum[i] -= temp_boltz[idx[seq_S]];
		    } 
	    }
	    #pragma omp barrier
        }// end while

//////////////////////////////
// END Part 2
/////////////////////////////


#endif 

#if PARTPARALLEL == 2
///////////////////////////////
// PART 2 Introducing Mutation
///////////////////////////////

	#pragma omp single
	{
	    // make sure smaller numbers don't get rounded off
	    sort(idx, idx+N*progeny, [&temp_boltz](size_t i1, size_t i2){return temp_boltz[i1] < temp_boltz[i2];});
            part_sum[0] = temp_boltz[idx[0]];
            for(long i=1; i<N*progeny; i++){
                part_sum[i] = part_sum[i-1] + temp_boltz[idx[i]];
            }
	    // pick sequences that survive to the final population based on their fitness
            Z = 0.0;
            Z_eff = 0.0;
            E_ave = 0.0;
            E_eff_ave = 0.0;
	}

	while(pos < N){
	    int pos_L = pos + rank;
	    double ranNum1 = u(r)*safeUnity*part_sum[N*progeny-1];
        long seq = std::lower_bound(part_sum, part_sum + N*progeny,ranNum1) - part_sum;
	    int active_threads = N > pos + (size - 1) ? size : N - pos;
	    seq_arr[rank] = seq;
  	    #pragma omp barrier
	    #pragma omp single
	    {
		int dup_idx;
		while( (dup_idx = check_duplicate(seq_arr, active_threads)) != 0){
		    ranNum1 = u(r)*safeUnity*part_sum[N*progeny-1];
            seq_arr[dup_idx] = std::lower_bound(part_sum, part_sum + N*progeny,ranNum1) - part_sum;
		}
	    } 
	    seq = seq_arr[rank];
	    if(seq >= N*progeny){
		cout << "Error in find, seq = " << seq << ", ranNum1 = " << ranNum1 << ", Ztemp = " << Ztemp << endl;
	    }
	    if(pos_L < N){
		    for(long j=0; j<m; j++){
		        population[pos_L][j]=temp_pop[idx[seq]][j];
		    }
		// prime T cells
		double penalty = 0;
		for(long j=0; j<n_epitope; j++){
		    int index = 0;
		    int count = 0;
		    for(long position=epi_start[j]; position<epi_end[j]; position++){
			index += place_value[j][count]*population[pos_L][position];
			count++;
		    }
		    double num_Ecells = 0;
		    for(long t=1; t<=rep_lim; t++){
			num_Ecells += Tcells[j][t];
		    } 
		    penalty += T_penalty*chi[j][index]*num_Ecells;
		    #pragma omp atomic
		    n_WTepitopes[j][index]+=1.0;
		}
		E[pos_L] = energy(population[pos_L],m,h,J); 
                E_eff[pos_L] = E[pos_L] + penalty;
                boltzfac[pos_L] = exp(-(E[pos_L])/T);
                boltzfac_eff[pos_L] = exp(-E_eff[pos_L]/T);

		#pragma omp critical
		{
            E_ave += E[pos_L];
            E_eff_ave += E_eff[pos_L];
		    Z += boltzfac[pos_L]; 
		    Z_eff += boltzfac_eff[pos_L]; 
		}
			
		occupied[idx[seq]] = false;
	    }
		
 	    for(int j = 0; j < active_threads; j++){
		#pragma omp barrier
		#pragma omp for schedule(static)
		for(long i=seq_arr[j]; i<N*progeny; i++){
		    part_sum[i] -= temp_boltz[idx[seq_arr[j]]];
		}
	    }

	    #pragma omp single
	    {
		pos += size;
	    }
	}

//////////////////////////////
// END Part 2
/////////////////////////////

#endif

#if PARTPARALLEL == 3
///////////////////////////////
// PART 2 Introducing Mutation
///////////////////////////////

	#pragma omp single
	{
        // initialize
            Z = 0.0;
            Z_eff = 0.0;
            E_ave = 0.0;
            E_eff_ave = 0.0;
        seq_S = 0;
	}
		
    // algorithm A, select numbers with largest keys
    #pragma omp for schedule(static)
    for(long i = 0; i < N*progeny; i++)	{
       key[i] = (1.0 / temp_boltz[idx[i]]) * std::log(u(r)); // take log for accuracy
       key_l[i] = key[i];
	}
    #pragma omp barrier
			
    // choose N largest elements
    #pragma omp single
    {
       std::nth_element(key, key + N, key + N*progeny, std::greater<double>());
	
       for(long i = 0; i < N*progeny; i++){
          if(key_l[i] > key[N]){
             seq_arr[seq_S] = i;
             seq_S++;
		   }
	   }
       printf("%d matched bigger than %e\n", seq_S, key[N]);
	}
		
    #pragma omp for schedule(dynamic) reduction(+: E_ave, E_eff_ave, Z, Z_eff)
	for(int pos_L = 0; pos_L < N; pos_L++){
		int seq = seq_arr[pos_L];
		for(long j=0; j<m; j++)	{
			population[pos_L][j]=temp_pop[idx[seq]][j];
		}
		// prime T cells
		double penalty = 0;
		for(long j=0; j<n_epitope; j++){
			int index = 0;
			int count = 0;
			for(long position=epi_start[j]; position<epi_end[j]; position++){
				index += place_value[j][count]*population[pos_L][position];
				count++;
			}
			double num_Ecells = 0;
			for(long t=1; t<=rep_lim; t++){
				num_Ecells += Tcells[j][t];
			} 
			penalty += T_penalty*chi[j][index]*num_Ecells;
			#pragma omp atomic
			    n_WTepitopes[j][index]+=1.0;
		}
		E[pos_L] = energy(population[pos_L],m,h,J); 
		E_eff[pos_L] = E[pos_L] + penalty; 
		boltzfac[pos_L] = exp(-E[pos_L]/T); 
		boltzfac_eff[pos_L] = exp(-E_eff[pos_L]/T); 
		
		E_ave += E[pos_L]; 
		E_eff_ave += E_eff[pos_L];
		Z += boltzfac[pos_L]; 
		Z_eff += boltzfac_eff[pos_L];
        occupied[idx[seq]] = false;
	}  
	
	#pragma omp barrier


//////////////////////////////
// END Part 2
/////////////////////////////
#endif

    #pragma omp for
    for(int i=0; i<n_epitope; i++){
        chiI[i] = 0.0;
        for(int j=0; j<(chi[i]).size(); j++){
            chiI[i] += chi[i][j]*n_WTepitopes[i][j];  
        }      
     }

    #pragma omp single
    {
    double t = 0.0;
    int react_num = 6+3*rep_lim-1;   //number of reactions for a single T-cell line
    std::vector<double> rates(n_epitope*react_num);
    double rate_tot = 0.0;
    int totalT = 0;
    int count = 0;
//cout << "rate_tot at start: " << rate_tot;
    for(long i=0; i<n_epitope; i++){
        // 4 reactions for naive cells: creation, replication, death, and activation
        rates[i*react_num+0] = w;
        rates[i*react_num+1] = Tcells[i][0]*b;
        rates[i*react_num+2] = Tcells[i][0]*e;
        rates[i*react_num+3] = Tcells[i][0]*a*chiI[i];
        totalT += Tcells[i][0];
        // 2 reactions for naive cells: activation and death
        rates[i*react_num+4] = Tcells[i][rep_lim+1]*ap*chiI[i];
        rates[i*react_num+5] = Tcells[i][rep_lim+1]*hr;
        totalT += Tcells[i][rep_lim+1];
        // effector cells becoming memory cells
        for(long j=0; j<rep_lim; j++){
            rates[i*react_num+6+j] = Tcells[i][j+1]*g;
	       // while we are doing this loop addup the total number of effector cells of all types
	       totalT += Tcells[i][j+1];
        }
        // effector cells dying and replicating
        for(long j=0; j<rep_lim-1; j++){
            rates[i*react_num+6+j+rep_lim] = Tcells[i][j+1]*d;
            rates[i*react_num+6+j+2*rep_lim] = Tcells[i][j+1]*re;
        }
        // terminal effector cells dying
        rates[i*react_num+6+2*rep_lim-1] = Tcells[i][rep_lim]*dp;
    }
    // find total rate
    for(long i=0; i<react_num*n_epitope; i++){
        rate_tot += rates[i];
    }

    t += -log(u(r)*safeUnity)/rate_tot;
    // reactions occur until a day has passed
    while(t<=7.0){
    // calculate rates
//        rate_tot = 0.0;
//        totalT = 0;
    
        count++;    
        ranNum1 = u(r)*safeUnity*rate_tot;
        int reaction = 0;
        double temp = 0.0;
        for(long i=0; i<react_num*n_epitope; i++){
            temp += rates[i];
            if(ranNum1<temp){
                reaction = i;
                break;
            }
        }
        int cell_line = reaction/react_num;
        int cell_type = reaction % react_num;
        int gen = 0;
        if(cell_type < 6){
            switch (cell_type){
                case 0:
                    if(totalT < n_T){
                    rates[cell_line*react_num+1] += b;
                    rates[cell_line*react_num+2] += e;
                    rates[cell_line*react_num+3] += chiI[cell_line]*a;
                    rate_tot += b + e + chiI[cell_line]*a;
                    Tcells[cell_line][0]++;
                    totalT++;
                    }
                    break;
                case 1:
                    if(totalT < n_T){
                    rates[cell_line*react_num+1] += b;
                    rates[cell_line*react_num+2] += e;
                    rates[cell_line*react_num+3] += chiI[cell_line]*a;
                    rate_tot += b + e + chiI[cell_line]*a;
                    Tcells[cell_line][0]++;
                    totalT++;
                    }
                    break;
                case 2:
                    rates[cell_line*react_num+1] -= b;
                    rates[cell_line*react_num+2] -= e;
                    rates[cell_line*react_num+3] -= chiI[cell_line]*a;
                    rate_tot -= b + e + chiI[cell_line]*a;
                    Tcells[cell_line][0]--;
                    totalT--;
                    break;
                case 3:
                    // only let the effector cell replicate if the population cap hasn't be hit.
//            	    if(totalT < n_T){
                        rates[cell_line*react_num+1] -= b;
                        rates[cell_line*react_num+2] -= e;
                        rates[cell_line*react_num+3] -= chiI[cell_line]*a;
                        rate_tot -= b + e + chiI[cell_line]*a;
                        rates[cell_line*react_num+6] += g;
                        rates[cell_line*react_num+6+rep_lim] += d;
                        rates[cell_line*react_num+6+2*rep_lim] += re;
                        rate_tot += g + d + re;
                        Tcells[cell_line][0]--;
                        Tcells[cell_line][1]++;
//                        totalT++;
//                    }
                    break;
                case 4:
                    // only let the effector cell replicate if the population cap hasn't be hit.
//            	    if(totalT < n_T){
                        rates[cell_line*react_num+4] -= chiI[cell_line]*ap;
                        rates[cell_line*react_num+5] -= hr;
                        rate_tot -= chiI[cell_line]*ap + hr;
                        rates[cell_line*react_num+6] += g;
                        rates[cell_line*react_num+6+rep_lim] += d;
                        rates[cell_line*react_num+6+2*rep_lim] += re;
                        rate_tot += g + d + re;
                        Tcells[cell_line][rep_lim+1]--;
                        Tcells[cell_line][1]++; 
//                        totalT++;
//                    }
                    break;
                case 5:
                    rates[cell_line*react_num+4] -= chiI[cell_line]*ap;
                    rates[cell_line*react_num+5] -= hr;
                    rate_tot -= chiI[cell_line]*ap + hr;
                    Tcells[cell_line][rep_lim+1]--;
                    totalT--;
                    break;
            }
        } else if(cell_type >= 6 && cell_type < (6+rep_lim)) {
            gen = (cell_type - 6) % rep_lim;
            if(gen < rep_lim-1){
                rates[cell_line*react_num+cell_type] -= g;
                rates[cell_line*react_num+cell_type+rep_lim] -= d;
                rates[cell_line*react_num+cell_type+2*rep_lim] -= re;
                rate_tot -= g + d + re;
                rates[cell_line*react_num+4] += chiI[cell_line]*ap;
                rates[cell_line*react_num+5] += hr;
                rate_tot += chiI[cell_line]*ap + hr;
            } else {
                rates[cell_line*react_num+cell_type] -= g;
                rates[cell_line*react_num+cell_type+rep_lim] -= dp;
                rate_tot -= g + dp;
                rates[cell_line*react_num+4] += chiI[cell_line]*ap;
                rates[cell_line*react_num+5] += hr;
                rate_tot += chiI[cell_line]*ap + hr;
            }
            Tcells[cell_line][gen+1]--;
            Tcells[cell_line][rep_lim+1]++;
//            totalT--;
        } else if(cell_type >= (6+rep_lim) && cell_type < (6+2*rep_lim)) {
            gen = (cell_type - 6) % rep_lim;
            if(gen < rep_lim-1){
                rates[cell_line*react_num+cell_type-rep_lim] -= g;
                rates[cell_line*react_num+cell_type] -= d;
                rates[cell_line*react_num+cell_type+rep_lim] -= re;
                rate_tot -= g + d + re;
            } else {
                rates[cell_line*react_num+cell_type-rep_lim] -= g;
                rates[cell_line*react_num+cell_type] -= dp;
                rate_tot -= g + dp;
            }
            Tcells[cell_line][gen+1]--;
            totalT--;
        } else if(cell_type >= (6+2*rep_lim) && cell_type <(6+3*rep_lim-1)) {
            // only let the effector cell replicate if the population cap hasn't be hit.
            if(totalT < n_T){
                gen = (cell_type - 6) % rep_lim;
                rates[cell_line*react_num+cell_type-2*rep_lim] -= g;
                rates[cell_line*react_num+cell_type-rep_lim] -= d;
                rates[cell_line*react_num+cell_type] -= re;
                rate_tot -= g + d + re;
                if(gen < rep_lim-2){
                    rates[cell_line*react_num+cell_type+1-2*rep_lim] += 2*g;
                    rates[cell_line*react_num+cell_type+1-rep_lim] += 2*d;
                    rates[cell_line*react_num+cell_type+1] += 2*re;
                    rate_tot += 2*(g + d + re);
                } else {
                    rates[cell_line*react_num+cell_type+1-2*rep_lim] += 2*g;
                    rates[cell_line*react_num+cell_type+1-rep_lim] += 2*dp;
                    rate_tot += 2*(g + dp);
                }
                Tcells[cell_line][gen+1]--;
                Tcells[cell_line][gen+2] += 2;
                totalT++;
            }
        } else {
            cout << "ERROR in picking T-cell reaction .... exiting" << endl;
            exit;
//            return -1;
        }
    t += -log(u(r)*safeUnity)/rate_tot;
    }
//cout << "\trate_tot at end: " << rate_tot << endl;
//cout << "Reaction this iteration (" << cycle << "): " << count << endl;
    } //end omp single
        
    #pragma omp barrier
    #pragma omp single
    {
        Z /= N;
        Z_eff /= N;
        Ztemp /= (N*progeny);
        E_ave /= N;
        E_eff_ave /= N;
    }
        
    // sample block
    if (!(cycle % sample_mod) && cycle > burnin) {
        for(long k=0; k<N; k++){
            for(long i=begin_m; i<end_m; i++) {
                n1[i][population[k][i]] += 1.0;
                for(long j=i+1; j<m; j++) {
                    n2[i][j][population[k][i]][population[k][j]] += 1.0;
	            }
	        }
        }

        #pragma omp single
        {
            n_samples += N;
	    }
    }
        
    #pragma omp barrier
    #pragma omp single
    {

	// print block
	if (cycle % print_mod == 0 && cycle > 0) {
	    cout << "  Completed cycle " << std::setw(5) << cycle << " of " << std::setw(5) << n_cycles << "\n";
            fflush(stdout);
	}
		
		
	// write block
        if (cycle % write_mod == 0 && cycle > burnin) {
            fwrite(&cycle, sizeof(long), 1, fout_pop_stats);
            fwrite(&pos, sizeof(long), 1, fout_pop_stats);
            fwrite(&Ztemp, sizeof(double), 1, fout_pop_stats);
            fwrite(&Z, sizeof(double), 1, fout_pop_stats);
            fwrite(&Z_eff, sizeof(double), 1, fout_pop_stats);
            fwrite(&E_ave, sizeof(double), 1, fout_pop_stats);
            fwrite(&E_eff_ave, sizeof(double), 1, fout_pop_stats);

            fwrite(&cycle, sizeof(long), 1, fout_P1_traj);
            for(long i=0; i<m; i++) {
                for(long p=0; p<nRes[i]; p++) {
                     n1[i][p] = n1[i][p]/n_samples;
                     fwrite(&n1[i][p], sizeof(double), 1, fout_P1_traj);
                     n1[i][p] = n1[i][p]*n_samples;
	            }
	        }
			
            fwrite(&cycle, sizeof(long), 1, fout_P2_traj);
            for(long i=0; i<m; i++) {
                for(long j=i+1; j<m; j++) {
                    for(long q=0; q<nRes[j]; q++) {
                        for(long p=0; p<nRes[i]; p++) {
                            n2[i][j][p][q] = n2[i][j][p][q]/n_samples;
                            fwrite(&n2[i][j][p][q], sizeof(double), 1, fout_P2_traj);
                            n2[i][j][p][q] = n2[i][j][p][q]*n_samples;
			}
		    }
		}
	    }

        fwrite(&cycle, sizeof(long), 1, fout_Tcell_traj);
        for(long i=0; i<n_epitope; i++){
             double Etotal = 0;
             fwrite(&Tcells[i][0], sizeof(double), 1, fout_Tcell_traj);
             for(long j=1; j<=rep_lim; j++){
                 Etotal += Tcells[i][j];
             }
             fwrite(&Etotal, sizeof(double), 1, fout_Tcell_traj);
             fwrite(&Tcells[i][rep_lim+1], sizeof(double), 1, fout_Tcell_traj);
             fwrite(&chiI[i], sizeof(double), 1, fout_Tcell_traj);
        }
			
        fwrite(&cycle, sizeof(long), 1, fout_MC_seqs);
	    for(long k=0; k<N; k++){
             for(long i=0; i<m; i++) {
                fwrite(&resIdx[i][population[k][i]], sizeof(int8_t), 1, fout_MC_seqs);
            }
	    }

         for(long j=0; j<n_epitope; j++){
             fwrite(&cycle, sizeof(long), 1, fout_epi_traj[j]);
             for(long k=0; k<N; k++){
                 for(long i=epi_start[j]; i<epi_end[j]; i++) {
                     fwrite(&resIdx[i][population[k][i]], sizeof(int8_t), 1, fout_epi_traj[j]);
                 }
             }
         }
	
         fwrite(&cycle, sizeof(long), 1, fout_fitness);
         fwrite(E, sizeof(double), N, fout_fitness);
         fwrite(E_eff, sizeof(double), N, fout_fitness);
         fwrite(boltzfac, sizeof(double), N, fout_fitness);
         fwrite(boltzfac_eff, sizeof(double), N, fout_fitness);
	}
    } // end #pragma omp single
    #pragma omp barrier
    } // end main loop
    } // end #pragma omp parallel
    cout << "DONE!" << endl << endl;

    delete[] boltzfac;
    delete[] temp_boltz;
    delete[] part_sum;
    delete[] boltzfac_eff;
    delete[] E;
    delete[] E_eff;
    delete[] idx;
    delete[] occupied;
//	delete []seq_arr;
	delete []part_sum_used;
    delete[] key;
    delete[] key_l;
	
    // closing traj files
    fclose(fout_P1_traj);
    fclose(fout_P2_traj);
    fclose(fout_MC_seqs);
    fclose(fout_pop_stats);
    fclose(fout_Tcell_traj);
    fclose(fout_fitness);
    for(int i=0; i<n_epitope; i++){
        fclose(fout_epi_traj[i]);
    }
	
	
    // reporting final probabilities
    FILE* fout_P1;
    {
	std::string fstr_P1="./P1_model.dat";
        fout_P1 = fopen(fstr_P1.c_str(),"wb");
        if (fout_P1==NULL) {
	    cerr << "Cannot open output file " << fstr_P1 << " in the current directory; aborting." << endl;
	    exit(-1);
	}
        int32_t sizes [1];
        sizes[0] = sizeof(double);
        fwrite(sizes,4,1,fout_P1);    
    } 
	

    FILE* fout_P2;
    {
	std::string fstr_P2="./P2_model.dat";
        fout_P2 = fopen(fstr_P2.c_str(),"wb");
        if (fout_P2==NULL) {
	    cerr << "Cannot open output file " << fstr_P2 << " in the current directory; aborting." << endl;
	    exit(-1);
	}
        int32_t sizes [1];
        sizes[0] = sizeof(double);
        fwrite(sizes,4,1,fout_P2);    
    } 

    for(long k=0; k<N; k++){
        for(long i=0; i<m; i++) {
            n1[i][population[k][i]] += 1.0;
            for(long j=i+1; j<m; j++) {
                n2[i][j][population[k][i]][population[k][j]] += 1.0;
	    }
	}
    }
    n_samples += N;
	
    for(long i=0; i<m; i++) {
        for(long p=0; p<nRes[i]; p++) {
              n1[i][p] = n1[i][p]/n_samples;
              fwrite(&n1[i][p], sizeof(double), 1, fout_P1);
	}
    }
		
    for(long i=0; i<m; i++) {
        for(long j=i+1; j<m; j++) {
            for(long q=0; q<nRes[j]; q++) {
                for(long p=0; p<nRes[i]; p++) {
                    n2[i][j][p][q] = n2[i][j][p][q]/n_samples;
                    fwrite(&n2[i][j][p][q], sizeof(double), 1, fout_P2);
                }
            }
        }
    }
	
    fclose(fout_P1);
    fclose(fout_P2);

#if PARTPARALLEL == 0 || PARTPARALLEL == 3
    FILE* seq_file;
    // write unsorted
    seq_file = fopen("seq_unsorted.txt", "w");
    for(long i =0; i < N; i++)
    	fprintf(seq_file, "%d\n", seq_arr[i]);
    fclose(seq_file);

    // write sorted
    sort(seq_arr, seq_arr + N);
    seq_file = fopen("seq_sorted.txt", "w");
    for(long i =0; i < N; i++)
    	fprintf(seq_file, "%d\n", seq_arr[i]);
    fclose(seq_file);
#endif // check sequence
    delete []seq_arr;


	
    // stop timer
    double wall_stop = get_wall_time();
    double cpu_stop  = get_cpu_time();
	
    printf("Wall Time = %e s\n",wall_stop - wall_start);
    printf("CPU Time  = %e s\n",cpu_stop - cpu_start);
	
    return 0;
}