From: Don Armstrong <don@donarmstrong.com>
Date: Fri, 7 Feb 2014 23:18:51 +0000 (-0800)
Subject: import initial version of genome
X-Git-Url: https://git.donarmstrong.com/?a=commitdiff_plain;h=d49beaff6c97583411fd8c554217da48b86d0e43;p=genome.git

import initial version of genome
---

d49beaff6c97583411fd8c554217da48b86d0e43
diff --git a/Error.cpp b/Error.cpp
new file mode 100644
index 0000000..5701232
--- /dev/null
+++ b/Error.cpp
@@ -0,0 +1,88 @@
+////////////////////////////////////////////////////////////////////// 
+// Error.cpp 
+//////////////////////////////////////////////////////////////////////////////
+//              COPYRIGHT NOTICE FOR GENOME CODE
+//
+// Copyright (C) 2006 - 2009, Liming Liang and Goncalo Abecasis,
+// 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.
+//
+///////////////////////////////////////////////////////////////////////////////
+
+#include "Error.h"
+
+#include "stdlib.h"
+#include "stdarg.h"
+#include "stdio.h"
+
+// Declare a dummy class to ensure that compilers recognize this as C++ code
+class String;
+
+void error ( const char * msg, ... )
+   {
+   va_list  ap;
+
+   va_start(ap, msg);
+
+   printf("\nFATAL ERROR - \n");
+   vprintf(msg, ap);
+   printf("\n\n");
+
+   va_end(ap);
+
+   exit(EXIT_FAILURE);
+   }
+
+void warning ( const char * msg, ... )
+   {
+   va_list  ap;
+
+   va_start(ap, msg);
+
+   printf("\n\aWARNING - \n");
+   vprintf(msg, ap);
+   printf("\n");
+
+   va_end(ap);
+   }
+
+void numerror ( const char * msg , ... )
+   {
+   va_list  ap;
+
+   va_start(ap, msg);
+
+   printf("\nFATAL NUMERIC ERROR - ");
+   vprintf(msg, ap);
+   printf("\n\n");
+
+   va_end(ap);
+
+   exit(EXIT_FAILURE);
+   }
diff --git a/Error.h b/Error.h
new file mode 100644
index 0000000..7235025
--- /dev/null
+++ b/Error.h
@@ -0,0 +1,55 @@
+////////////////////////////////////////////////////////////////////// 
+// Error.h
+//////////////////////////////////////////////////////////////////////////////
+//              COPYRIGHT NOTICE FOR GENOME CODE
+//
+// Copyright (C) 2006 - 2009, Liming Liang and Goncalo Abecasis,
+// 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.
+//
+///////////////////////////////////////////////////////////////////////////////
+
+
+#ifndef _ERROR_H_
+#define _ERROR_H_
+
+// #ifdef __cplusplus
+// extern "C" {
+// #endif
+
+void error(const char * msg, ...);
+void warning(const char * msg, ...);
+void numerror(const char * msg, ...);
+
+// #ifdef __cplusplus
+//   };
+// #endif
+
+
+#endif
diff --git a/LICENSE.GENOME b/LICENSE.GENOME
new file mode 100644
index 0000000..69e9ac4
--- /dev/null
+++ b/LICENSE.GENOME
@@ -0,0 +1,36 @@
+
+COPYRIGHT NOTICE FOR GENOME
+=====================================
+
+   GENOME coded by Liming Liang and Goncalo Abecasis.   
+
+   Copyright (C) 2006 - 2007, Liming Liang and Goncalo Abecasis,
+   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.
+
diff --git a/LICENSE.twister b/LICENSE.twister
new file mode 100644
index 0000000..2d5b11b
--- /dev/null
+++ b/LICENSE.twister
@@ -0,0 +1,37 @@
+Mersenne twister code is included in the file Random.cpp
+
+COPYRIGHT NOTICE FOR MERSENNE TWISTER
+=====================================
+
+   Mersenne twister coded by Takuji Nishimura and Makoto Matsumoto.   
+
+   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.
+
diff --git a/Main.cpp b/Main.cpp
new file mode 100644
index 0000000..7beabb9
--- /dev/null
+++ b/Main.cpp
@@ -0,0 +1,211 @@
+////////////////////////////////////////////////////////////////////// 
+// Main.cpp 
+//////////////////////////////////////////////////////////////////////////////
+//              COPYRIGHT NOTICE FOR GENOME CODE
+//
+// Copyright (C) 2006 - 2009, Liming Liang and Goncalo Abecasis,
+// 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.
+//
+///////////////////////////////////////////////////////////////////////////////
+
+
+#include <iostream>
+#include <math.h>
+#include <vector> 
+#include <time.h>
+#include "genome.h"
+#include "Random.h"
+#include <stdio.h> 
+#include <iomanip>
+    
+using namespace std;
+
+
+int main(int argc, char ** argv)
+   {
+	string popSize("10000"), rec("0.0001");
+	
+	int numPieces=100, pieceLen=10000, numIndepRegion=1, SNP=-1;
+	int nSubPOP=2; 
+	vector<int> nSubSample(nSubPOP);
+	for(int i=0;i<nSubPOP;i++)nSubSample[i]=10;
+	
+	double mut=1e-8, mig=2.5e-4;
+	double maf=0.0,prop=1.0;
+	
+	vector< vector<bool> > data;
+
+	long seed=-1;
+	
+	int drawtree=0;
+
+	if(argc==1){
+		cout << endl  << "GENOME-0.2: Whole Genome Coalescent Simulator. (2009.2) Liming Liang, Goncalo Abecasis"  << endl << endl;
+		cout << "     Parameters and default values:" << endl;
+		cout << "     -pop     " << "number of subpopulations and size of subsamples [ " <<nSubPOP<<" "; 
+									for(int i=0;i<nSubPOP;i++)cout<<nSubSample[i]<<" ";cout<<"]"<< endl;
+		cout << "     -N       " << "effective size of each subpopulation or the filename for population profile [" <<popSize.c_str()<<"]"<< endl;
+		cout << "     -c       " << "number of independent regions (chromosomes) to simulate [" <<numIndepRegion<<"]"<< endl;
+		cout << "     -pieces  " << "number of fragments per independent region [" <<numPieces<<"]"<< endl;
+		cout << "     -len     " << "length in base of each fragment [" <<pieceLen<<"]"<< endl;
+		cout << "     -s       " << "fixed number of SNPs per independent region (chromosome), -1 = number of SNPs follows Poisson distribution [" <<SNP<<"]"<< endl;
+		cout << "     -rec     " << "recombination rate bewteen consecutive fragments"<<" or the filename for recombination rate distribution ["<<rec.c_str()<<"]"<< endl;
+		cout << "     -mut     " << "mutation rate per generation per base pair [" <<mut<<"]"<< endl;
+		cout << "     -mig     " << "migration rate per generation per individual [" <<mig<<"]"<< endl;
+		cout << "     -seed    " << "random seed, -1 = use time as the seed [" <<seed<<"]"<< endl;
+		cout << "     -tree    " << "1=draw the genealogy trees, 0=do not output [" <<drawtree<<"]"<< endl;
+		cout << "     -maf     " << "Output SNPs with minor allele frequency greater than [" <<maf<<"]"<< endl;
+		cout << "     -prop    " << "To keep this proportion of SNPs with MAF < the value of -maf parameter [" <<prop<<"]"<< endl;
+		exit(0);
+	}
+ 
+
+	int k=0;
+    for(int i=1;i<argc;i++){
+        if(strcmp(argv[i],"-pop")==0){ 
+	        nSubPOP=atoi(argv[++i]); k++;
+	        nSubSample.clear();
+	        nSubSample.resize(nSubPOP);
+	        for(int iter=0;iter<nSubPOP;iter++)nSubSample[iter]=atoi(argv[++i]);
+	    }
+	    else if(strcmp(argv[i],"-N")==0){ popSize.assign(argv[++i]); k++;}
+	    else if(strcmp(argv[i],"-c")==0){ numIndepRegion=atoi(argv[++i]); k++;}
+	    else if(strcmp(argv[i],"-pieces")==0){ numPieces=atoi(argv[++i]); k++;}
+	    else if(strcmp(argv[i],"-len")==0){ pieceLen=atoi(argv[++i]); k++;}
+	    else if(strcmp(argv[i],"-s")==0){ SNP=atoi(argv[++i]); k++;}
+	    else if(strcmp(argv[i],"-rec")==0){ rec.assign(argv[++i]); k++;}
+	    else if(strcmp(argv[i],"-mut")==0){ mut=atof(argv[++i]); k++;}
+	    else if(strcmp(argv[i],"-mig")==0){ mig=atof(argv[++i]); k++;}
+	    else if(strcmp(argv[i],"-seed")==0){ seed=atol(argv[++i]); k++;}
+	    else if(strcmp(argv[i],"-tree")==0){ drawtree=atoi(argv[++i]); k++;}
+	    else if(strcmp(argv[i],"-maf")==0){ maf=atof(argv[++i]); k++;}
+	    else if(strcmp(argv[i],"-prop")==0){ prop=atof(argv[++i]); k++;}
+    }
+    
+    if(k!=(argc-1)/2 && k!=(argc-1-nSubPOP)/2){
+        cout << "Parameters do not match! " << argc<< " " << k << endl; 
+        exit(1);
+    } 
+
+    cout << endl 
+	    << "GENOME-0.2: Whole Genome Coalescent Simulator. (2009.2) Liming Liang, Goncalo Abecasis"  << endl << endl;
+	cout << "     Parameters and effective values:" << endl;
+	cout << "     -pop     " << "number of subpopulations and size of subsamples [ " <<nSubPOP<<" "; 
+									for(int i=0;i<nSubPOP;i++)cout<<nSubSample[i]<<" ";cout<<"]"<< endl;
+		cout << "     -N       " << "effective size of each subpopulation or the filename for population profile [" <<popSize.c_str()<<"]"<< endl;
+		cout << "     -c       " << "number of independent regions (chromosomes) to simulate [" <<numIndepRegion<<"]"<< endl;
+		cout << "     -pieces  " << "number of fragments per independent region [" <<numPieces<<"]"<< endl;
+		cout << "     -len     " << "length in base of each fragment [" <<pieceLen<<"]"<< endl;
+		cout << "     -s       " << "fixed number of SNPs per independent region (chromosome), -1 = number of SNPs follows Poisson distribution [" <<SNP<<"]"<< endl;
+		cout << "     -rec     " << "recombination rate bewteen consecutive fragments or the filename for recombination rate distribution ["<<rec.c_str()<<"]"<< endl;
+		cout << "     -mut     " << "mutation rate per generation per base pair [" <<mut<<"]"<< endl;
+		cout << "     -mig     " << "migration rate per generation per individual [" <<mig<<"]"<< endl;
+		cout << "     -seed    " << "random seed, -1 = use time as the seed [" <<seed<<"]"<< endl;
+		cout << "     -tree    " << "1=draw the genealogy trees, 0=do not output [" <<drawtree<<"]"<< endl;
+		cout << "     -maf     " << "Output SNPs with minor allele frequency greater than [" <<maf<<"]"<< endl;
+		cout << "     -prop    " << "To keep this proportion of SNPs with MAF < the value of -maf parameter [" <<prop<<"]"<< endl<<endl;
+
+	if(atoi(popSize.c_str())>0 && atof(rec.c_str())>0){	
+		printf("     Scaled recombination rate = %.4e\n", 2 * atoi(popSize.c_str()) * (numPieces - 1) * atof(rec.c_str()));
+		printf("     Scaled migration rate = %.4e\n", 2 * atoi(popSize.c_str()) * mig);
+		printf("     Scaled mutation rate = %.4e\n\n", 2 * atoi(popSize.c_str()) * mut * numPieces * pieceLen);
+	}
+	
+	genome(popSize, nSubPOP, nSubSample, numPieces, pieceLen, numIndepRegion, SNP, rec, mut, mig, data, seed,drawtree);
+	
+// output format for haploview
+
+// 	for(int i=0;i<data.size();i++){
+// 		printf("FAMILY1 MEMBER%d ",i/2+1);
+// 		for(int j=0;j<data[i].size();j++){
+// 			cout<<data[i][j]+1<<" ";
+// 		}
+// 		cout<<endl;
+// 	}
+
+
+// 0 = ancestral state, 1 = mutated state
+
+	cout<<"Total number of SNPs = "<<data[0].size()<<endl;
+	cout<<"Samples:"<<endl;
+	
+	if(maf>0.0 && prop < 1.0){ // filter out the SNPs with MAF < maf 
+		
+		cout<<"Rare SNPs filtered according to parameters."<<endl;
+		
+		int n=0;
+		for(int k=0;k<nSubPOP;k++)n+=nSubSample[k];
+		
+		int m=(int)(data[0].size());
+		vector<bool> keep(m,true);
+			
+		float p;
+		for(int i=0;i<m;i++){
+			p=0.0;
+			for(int j=0;j<n;j++){
+				p += data[j][i];	
+			}
+			p /= n;
+			
+			if( (p<maf || (1-p)<maf) && globalRandom.Next()>prop)keep[i]=0;
+		}
+		
+
+		int index=0;
+		for(int k=0;k<nSubPOP;k++){
+			for(int i=0;i<nSubSample[k];i++){
+				cout<<"POP"<<k+1<<": ";
+				for(unsigned int j=0;j<data[index].size();j++){
+					if(keep[j])cout<<data[index][j];
+				}
+				cout<<endl;
+				index++;
+			}	
+		}
+	
+	}
+	else{
+		int index=0;
+		for(int k=0;k<nSubPOP;k++){
+			for(int i=0;i<nSubSample[k];i++){
+				cout<<"POP"<<k+1<<": ";
+				for(unsigned int j=0;j<data[index].size();j++){
+					cout<<data[index][j];
+				}
+				cout<<endl;
+				index++;
+			}	
+		}
+	}
+
+  
+}
+
+
diff --git a/MathConstant.h b/MathConstant.h
new file mode 100644
index 0000000..d471dac
--- /dev/null
+++ b/MathConstant.h
@@ -0,0 +1,73 @@
+////////////////////////////////////////////////////////////////////// 
+// MathConstant.h 
+//////////////////////////////////////////////////////////////////////////////
+//              COPYRIGHT NOTICE FOR GENOME CODE
+//
+// Copyright (C) 2006 - 2009, Liming Liang and Goncalo Abecasis,
+// 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.
+//
+///////////////////////////////////////////////////////////////////////////////
+
+
+
+#ifndef __MATHCONSTANT_H__
+#define __MATHCONSTANT_H__
+
+#ifdef  _MSC_VER
+#define _USE_MATH_DEFINES 
+#endif
+
+#include "math.h"
+#include "stdlib.h"
+
+// Constants for numerical routines
+//
+
+#define TINY    1.0e-30        // A small number
+#define ITMAX   200            // Maximum number of iterations
+#define EPS     3.0e-7         // Relative accuracy
+#define ZEPS    3.0e-10        // Precision around zero
+#define FPMIN   1.0e-30        // Number near the smallest representable number
+#define FPMAX   1.0e+100       // Number near the largest representable number
+#define TOL     1.0e-6         // Zero SVD values below this
+#define GOLD    0.61803399     // Golden ratio
+#define CGOLD   0.38196601     // Complement of golden ratio
+
+inline double square(double a)         { return a * a; }
+inline double sign(double a, double b) { return b >= 0 ? fabs(a) : -fabs(a); }
+inline double min(double a, double b)  { return a < b ? a : b; }
+inline double max(double a, double b)  { return a > b ? a : b; }
+
+inline int square(int a)      { return a * a; }
+inline int sign(int a, int b) { return b >= 0 ? abs(a) : -abs(a); }
+inline int min(int a, int b)  { return a < b ? a : b; }
+inline int max(int a, int b)  { return a > b ? a : b; }
+
+#endif
diff --git a/Random.cpp b/Random.cpp
new file mode 100644
index 0000000..a181e97
--- /dev/null
+++ b/Random.cpp
@@ -0,0 +1,387 @@
+
+//////////////////////////////////////////////////////////////////////////////
+// This file includes code derived from the original Mersenne Twister Code
+// by Makoto Matsumoto and Takuji Nishimura
+// and is subject to their original copyright notice copied below:
+//////////////////////////////////////////////////////////////////////////////
+
+//////////////////////////////////////////////////////////////////////////////
+// COPYRIGHT NOTICE FOR MERSENNE TWISTER CODE
+// 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.
+//
+///////////////////////////////////////////////////////////////////////////////
+
+#include "Random.h"
+#include "MathConstant.h"
+#include "Error.h"
+
+#include <math.h>
+
+//Constants used internally by Mersenne random number generator
+#define MERSENNE_N 624
+#define MERSENNE_M 397
+
+// constant vector a
+#define MATRIX_A 0x9908b0dfUL
+
+// most significant w-r bits
+#define UPPER_MASK 0x80000000UL
+
+// least significant r bits
+#define LOWER_MASK 0x7fffffffUL
+
+
+// Constants used internally by Park-Miller random generator
+#define IA 16807
+#define IM 2147483647
+#define AM (1.0 / IM)
+#define IQ 127773
+#define IR 2836
+#define NTAB 32
+#define NDIV (1+(IM-1)/NTAB)
+#define RNMX (1.0-EPS)
+
+Random::Random(long s)
+   {
+#ifndef __NO_MERSENNE
+   mt = new unsigned long [MERSENNE_N];
+   mti = MERSENNE_N + 1;
+   mersenneMult = 1.0/4294967296.0;
+#else
+   shuffler = new long [NTAB];
+#endif
+   Reset(s);
+   }
+
+Random::~Random()
+   {
+#ifndef __NO_MERSENNE
+   delete [] mt;
+#else
+   delete [] shuffler;
+#endif
+   }
+
+void Random::Reset(long s)
+   {
+   normSaved = 0;
+
+#ifndef __NO_MERSENNE
+   InitMersenne(s);
+#else
+   // 'Continuous' Random Generator
+   if ((seed = s) < 1)
+      seed = s == 0 ? 1 : -s;  // seed == 0 would be disastrous
+
+   for (int j=NTAB+7; j>=0; j--)  // Warm up and set shuffle table
+      {
+      long k = seed / IQ;
+      seed = IA * (seed - k * IQ) - IR * k;
+      if (seed < 0) seed += IM;
+      if (j < NTAB) shuffler[j] = seed;
+      }
+   last=shuffler[0];
+#endif
+   }
+
+// initializes mt[MERSENNE_N] with a seed
+void Random::InitMersenne(unsigned long s)
+   {
+   mt[0]= s & 0xffffffffUL;
+   for (mti = 1; mti < MERSENNE_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;
+      }
+   }
+
+int Random::Binary()
+   {
+   return Next() > 0.5 ? 1 : 0;
+   }
+
+#ifndef __NO_MERSENNE
+
+double Random::Next()
+   {
+   unsigned long y;
+
+   // mag01[x] = x * MATRIX_A for x=0,1
+   static unsigned long mag01[2]={0x0UL, MATRIX_A};
+
+   if (mti >= MERSENNE_N)
+      {
+      /* generate MERSENNE_N words at one time */
+      int kk;
+
+      // If InitMersenne() has not been called, a default initial seed is used
+      if (mti == MERSENNE_N+1)
+         InitMersenne(5489UL);
+
+      for (kk=0; kk < MERSENNE_N-MERSENNE_M; kk++)
+         {
+         y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
+         mt[kk] = mt[kk+MERSENNE_M] ^ (y >> 1) ^ mag01[y & 0x1UL];
+         }
+
+      for (; kk < MERSENNE_N-1; kk++)
+         {
+         y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
+         mt[kk] = mt[kk+(MERSENNE_M - MERSENNE_N)] ^ (y >> 1) ^ mag01[y & 0x1UL];
+         }
+
+      y = (mt[MERSENNE_N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
+      mt[MERSENNE_N-1] = mt[MERSENNE_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 (mersenneMult * ((double) y + 0.5));
+   }
+
+// Generates a random number on [0,0xffffffff]-interval
+
+unsigned long Random::NextInt()
+   {
+   unsigned long y;
+
+   // mag01[x] = x * MATRIX_A for x=0,1
+   static unsigned long mag01[2]={0x0UL, MATRIX_A};
+
+   if (mti >= MERSENNE_N)
+      {
+      /* generate MERSENNE_N words at one time */
+      int kk;
+
+      // If InitMersenne() has not been called, a default initial seed is used
+      if (mti == MERSENNE_N + 1)
+         InitMersenne(5489UL);
+
+      for (kk= 0; kk < MERSENNE_N - MERSENNE_M; kk++)
+         {
+         y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
+         mt[kk] = mt[kk+MERSENNE_M] ^ (y >> 1) ^ mag01[y & 0x1UL];
+         }
+
+      for (; kk< MERSENNE_N-1; kk++)
+         {
+         y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
+         mt[kk] = mt[kk+(MERSENNE_M - MERSENNE_N)] ^ (y >> 1) ^ mag01[y & 0x1UL];
+         }
+
+      y = (mt[MERSENNE_N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
+      mt[MERSENNE_N-1] = mt[MERSENNE_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;
+   }
+
+#else
+
+double Random::Next()
+   {
+   // Compute seed = (IA * seed) % IM without overflows
+   // by Schrage's method
+   long k = seed / IQ;
+   seed = IA * (seed - k * IQ) - IR * k;
+   if (seed < 0) seed += IM;
+
+   // Map to 0..NTAB-1
+   int j = last/NDIV;
+
+   // Output value is shuffler[j], which is in turn replaced by seed
+   last = shuffler[j];
+   shuffler[j] = seed;
+
+   // Map to 0.0 .. 1.0 excluding endpoints
+   double temp = AM * last;
+   if (temp > RNMX) return RNMX;
+   return temp;
+   }
+
+unsigned long Random::NextInt()
+   {
+   // Compute seed = (IA * seed) % IM without overflows
+   // by Schrage's method
+   long k = seed / IQ;
+   seed = IA * (seed - k * IQ) - IR * k;
+   if (seed < 0) seed += IM;
+
+   // Map to 0..NTAB-1
+   int j = last/NDIV;
+
+   // Output value is shuffler[j], which is in turn replaced by seed
+   last = shuffler[j];
+   shuffler[j] = seed;
+
+   return last;
+   }
+
+#endif
+
+double Random::Normal()
+   {
+   double v1, v2, fac, rsq;
+
+   if (!normSaved)  // Do we need new numbers?
+      {
+      do {
+      v1 = 2.0 * Next() - 1.0;  // Pick two coordinates from
+      v2 = 2.0 * Next() - 1.0;  // -1 to +1 and check if they
+      rsq = v1*v1 + v2*v2;  // are in unit circle...
+      } while (rsq >= 1.0 || rsq == 0.0);
+
+      fac = sqrt(-2.0 * log(rsq)/rsq);  // Apply the Box-Muller
+      normStore = v1 * fac;  // transformation and save
+      normSaved = 1;  // one deviate for next time
+      return v2 * fac;
+      }
+   else
+      {
+      normSaved = 0;
+      return normStore;
+      }
+   }
+
+void Random::Choose(int * array, int n, int k)
+   {
+   int choices = 1, others = 0;
+
+   if (k > n / 2)
+      {
+      choices = 0;
+      others = 1;
+      k = n - k;
+      }
+
+   for (int i = 0; i < n; i++)
+      array[i] = others;
+
+   while (k > 0)
+      {
+      int i = NextInt() % n;
+
+      if (array[i] == choices) continue;
+
+      array[i] = choices;
+      k--;
+      }
+   }
+
+void Random::Choose(int * array, float * weights, int n, int k)
+   {
+   int choices = 1, others = 0;
+
+   if (k > n / 2)
+      {
+      choices = 0;
+      others = 1;
+      k = n - k;
+      }
+
+   // First calculate cumulative sums of weights ...
+   float * cumulative = new float [n + 1];
+
+   cumulative[0] = 0;
+   for (int i = 1; i <= n; i++)
+      cumulative[i] = cumulative[i - 1] + weights[i - 1];
+
+   float & sum = cumulative[n], reject = 0.0;
+
+   for (int i = 0; i < n; i++)
+      array[i] = others;
+
+   while (k > 0)
+      {
+      float weight = (float)(Next() * sum);
+
+      int hi = n, lo = 0, i = 0;
+
+      while (hi >= lo)
+         {
+         i = (hi + lo) / 2;
+
+         if (cumulative[i + 1] <= weight)
+            lo = i + 1;
+         else if (cumulative[i] >= weight)
+            hi = i - 1;
+         else break;
+         }
+
+      if (array[i] == choices) continue;
+
+      array[i] = choices;
+      reject += weights[i];
+
+      // After selecting a substantial number of elements, update the cumulative
+      // distribution -- to ensure that at least half of our samples produce a hit
+      if (reject > sum * 0.50)
+         {
+         cumulative[0] = 0;
+         for (int i = 1; i <= n; i++)
+            if (array[i] != choices)
+               cumulative[i] = cumulative[i - 1] + weights[i - 1];
+            else
+               cumulative[i] = cumulative[i - 1];
+
+         reject = 0.0;
+         }
+
+      k--;
+      }
+   }
+
+Random globalRandom;
+
+
diff --git a/Random.h b/Random.h
new file mode 100644
index 0000000..cc6900f
--- /dev/null
+++ b/Random.h
@@ -0,0 +1,115 @@
+
+//////////////////////////////////////////////////////////////////////////////
+// This file includes code derived from the original Mersenne Twister Code
+// by Makoto Matsumoto and Takuji Nishimura
+// and is subject to their original copyright notice copied below:
+//////////////////////////////////////////////////////////////////////////////
+
+//////////////////////////////////////////////////////////////////////////////
+//              COPYRIGHT NOTICE FOR MERSENNE TWISTER CODE
+//
+// 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.
+//
+///////////////////////////////////////////////////////////////////////////////
+
+
+#ifndef __RANDOM_H__
+#define __RANDOM_H__
+
+// Define a quick and dirty generator
+#define RANDMUL 1664525L
+#define RANDADD 1013904223L
+
+#define RAND(seed) ((seed = seed * RANDMUL + RANDADD) & 0xFFFFFFFF)
+
+class Random
+// Implements the Mersenne Twister as default random number generator.
+// Compilation flag __NO_MERSENNE sets default generator to
+// a minimal Park-Miller with Bays-Durham shuffle and added safe guards.
+   {
+   protected:
+      // values for "minimal random values"
+      long  seed;
+      long  last;
+      long  * shuffler;
+
+      // and for normal deviates
+      int      normSaved;
+      double   normStore;
+
+      double mersenneMult;
+
+      // Array for Mersenne state vector
+      unsigned long * mt;
+
+      // Used to signal that Mersenne state vector is not initialized
+      int mti;
+
+
+   public:
+
+      Random(long s = 0x7654321);
+      ~Random();
+
+      // Next bit in series of 0s and 1s
+      int    Binary();     // Next bit in series of 0s and 1s
+
+      // Next value in series, between 0 and 1
+      double Next();
+
+      // Next integer
+      unsigned long NextInt();
+
+      // Random number form N(0,1)
+      double Normal();
+
+      void   Reset(long s);
+      void   InitMersenne(unsigned long s);
+
+      // Random number between 0 and 1
+      operator double()
+         { return Next(); }
+
+      // Random number between arbitrary bounds
+      double Uniform(double lo = 0.0, double hi = 1.0)
+         {
+         return lo + (hi - lo) * Next();
+         }
+
+      void Choose(int * array, int n, int k);
+      void Choose(int * array, float * weights, int n, int k);
+
+   };
+
+extern Random globalRandom;
+
+#endif
+
diff --git a/genome.cpp b/genome.cpp
new file mode 100644
index 0000000..f46c661
--- /dev/null
+++ b/genome.cpp
@@ -0,0 +1,1396 @@
+////////////////////////////////////////////////////////////////////// 
+// genome.cpp
+//////////////////////////////////////////////////////////////////////////////
+//              COPYRIGHT NOTICE FOR GENOME CODE
+//
+// Copyright (C) 2006 - 2009, Liming Liang and Goncalo Abecasis,
+// Version 0.2
+// 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.
+//
+///////////////////////////////////////////////////////////////////////////////
+
+
+#define _CRT_SECURE_NO_DEPRECATE 
+
+#include "Random.h"
+#include "Error.h"
+#include "stochastic.h" 
+#include <iostream>
+#include <fstream>
+#include <math.h>
+#include <vector> 
+#include <time.h>
+#include <algorithm>
+#include <map>
+#include <sstream>
+#include <stdio.h>
+
+using namespace std;
+
+//static bool debug = false;
+static bool fixedPOP = true;
+static bool noChange = true;
+
+
+static int N = 10000, n = 100, L = 1, poolsize = 1024 * 1024 * 64, length=10000, numChr=1, SNP=-1; // length = number of base per fragment
+static int maxN;
+static unsigned int currentProfile;
+static int nPOP=1;
+static vector<int> nSample(nPOP,3);
+static double recombination = 1e-8, migration = 0.0, mutation=1e-8;
+
+static int * genepool, * genetimes, * fragments, * MRCA, * MRCAtimes;
+static int *** children, *** parents, * cblocks, * pblocks;
+static int * parent_First, * parent_Last, * child_First, * child_Last;
+static bool * child_flags, * parent_flags;
+static int nextFree = 0, activeSegments = 0, pblockcounter = 0;
+static int numBLOCKs;
+static int BUFFER=1500;
+static int INCREMENT=1000;
+
+#define BLOCKSIZE   128
+#define SWAP(tmp,a,b)  tmp=a; a=b; b=tmp;
+
+static vector< vector< vector<bool> > > chromosome;
+static int * fragmentLength;
+static int * branchFragment;
+static int * geneIndex;
+static int * branchMutation;
+
+static float * recombVector=NULL;
+
+struct popProfilesStruct{
+	int generation;
+	vector<int> popsize;
+	vector<int> popStart;
+	vector< vector<float> > selectProb;
+	vector< vector<int> > popChange;
+};
+
+static vector<struct popProfilesStruct> popProfile;
+
+void print(){
+	cout<<"\n*****************\nParents:\t";
+	for(int i=0;i<N;i++){
+		if(parent_flags[i]){
+			for(int j=0;j<L;j++) cout<<"["<<((parents[i][j / BLOCKSIZE]==NULL)?999:(parents[i][j / BLOCKSIZE][j % BLOCKSIZE]))<<"]";
+			cout<<"\t";	
+		}
+	}	
+	cout<<endl;
+	
+	cout<<"Children:\t";
+	for(int i=0;i<N;i++){
+		if(child_flags[i]){
+			for(int j=0;j<L;j++) cout<<"["<<((children[i][j / BLOCKSIZE]==NULL)?999:(children[i][j / BLOCKSIZE][j % BLOCKSIZE]))<<"]";
+			cout<<"\t";	
+		}
+	}	
+	cout<<endl;
+	
+	
+	
+// 	cout<<"parent_flag:\t";
+// 	for(int i=0;i<N;i++){
+// 		cout<<parent_flags[i]<<"\t";
+// 	}	
+// 	cout<<endl;
+// 	
+// 	cout<<"children_flag:\t";
+// 	for(int i=0;i<N;i++){
+// 		cout<<child_flags[i]<<"\t";
+// 	}	
+// 	cout<<endl;
+	
+	cout<<"genepool:\t";
+	for(int i=0;i<2*n*L-L;i++)cout<<genepool[i]<<" ";
+	cout<<endl;
+	
+	cout<<"genetimes:\t";
+	for(int i=0;i<2*n*L-L;i++)cout<<genetimes[i]<<" ";
+	cout<<endl;
+	
+	cout<<"MCRA:\t";
+	for(int i=0;i<L;i++)cout<<MRCA[i]<<" ";
+	cout<<endl;
+	
+	cout<<"MCRAtime:\t";
+	for(int i=0;i<L;i++)cout<<MRCAtimes[i]<<" ";
+	cout<<endl;
+	
+	cout<<"fragments:\t";
+	for(int i=0;i<L;i++)cout<<fragments[i]<<" ";
+	cout<<endl;
+	
+	cout<<"NextFree="<<nextFree<<endl;
+}
+
+int *** AllocateIntPtrMatrix(int rows, int columns){
+   int *** result = new int ** [rows];
+
+   for (int i = 0; i < rows; i++)
+      result[i] = new int * [columns];
+  
+   return result;
+}
+   
+void FreeIntPtrMatrix(int *** & matrix, int rows){
+	
+   for (int i = 0; i < rows; i++)
+      delete [] matrix[i];
+
+   delete [] matrix;
+
+   matrix = NULL;
+}   
+   
+void NewGeneration(){
+	
+   int *** ptrmatrix, * vector;
+   bool * boolvector;
+   
+//    for(int i=0;i<maxN;i++){
+// 	if(parent_flags[i])printf("numBlock=%d parent=%d first=%d last=%d\n",numBLOCKs,i,parent_First[i],parent_Last[i]);   
+//    }
+   
+   SWAP(ptrmatrix, children, parents);
+   SWAP(vector, child_First, parent_First);
+   SWAP(vector, child_Last, parent_Last);
+   SWAP(boolvector, child_flags, parent_flags);
+   SWAP(vector, cblocks, pblocks);
+   
+   if(!noChange){
+	   currentProfile++;
+	   nPOP = (int)(popProfile[currentProfile].popsize.size());
+	   N = popProfile[currentProfile].popStart.back()+popProfile[currentProfile].popsize.back();
+	   noChange = true;
+   }
+
+   for (int i = 0; i < maxN; i++)
+      parent_flags[i] = 0;
+      
+//    for (int i = 0; i < maxN; i++)
+//        for (int j = 0; j < numBLOCKs; j++)
+//            parents[i][j] = NULL;
+      
+   pblockcounter = 0;
+
+}
+
+   
+void AllocateMemory(vector< vector<bool> > &return_chromosome){
+   
+   
+//    if(SNP<=0){
+//  	   	cout<<"Reserve size for return_chromosome="<<max(1,(int)(mutation*(double)numChr*(double)L*(double)length*40*(double)N/(double)nPOP))<<endl;
+//    }
+//    else{
+//   		cout<<"Reserve size for return_chromosome="<<SNP*numChr<<endl;
+//    }
+//    cout<<"Reserve size for working chromosome="<<max(1,(int)(mutation*(double)length*40*(double)N/(double)nPOP))<<endl;
+   
+ 
+   return_chromosome.clear();
+   return_chromosome.resize(n);
+//    for(int i=0;i<n;i++){
+// 	   if(SNP<=0){
+//  	   		//return_chromosome[i].reserve(max(1,40*(int)(mutation*(double)numChr*(double)L*(double)length*(double)maxN/(double)nPOP)));
+// 	   }
+//  	   else{
+//  	   		//return_chromosome[i].reserve(SNP*numChr);
+// 	   }
+// 	}
+   
+   chromosome.resize(n);
+   for(int i=0;i<n;i++){
+ 		chromosome[i].resize(L);
+// 		for(int j=0;j<L;j++){
+//  		   if(SNP<=0){
+// 	 		   //chromosome[i][j].reserve(max(1,40*(int)(mutation*(double)length*(double)maxN/(double)nPOP)));
+// 		   }
+// 	 	   else{
+// 	 	   	   //chromosome[i][j].reserve(SNP/L+1);
+// 		   }
+// 		}
+   }
+  
+   genepool = new int [poolsize];
+
+   numBLOCKs= (int)ceil((double)L / (double)BLOCKSIZE);
+ 
+   pblockcounter = 0;
+
+}
+
+
+void AllocateMutationMemory(){
+	
+	fragmentLength = new int [L];
+	
+	branchFragment = new int [poolsize];
+	
+	geneIndex = new int [poolsize];   
+	   
+}
+   
+   
+   
+void FreeMemory(){
+	
+   delete [] genepool;
+   
+   if(recombVector!=NULL)delete [] recombVector;
+   
+}
+   
+  
+void FreeMutationMemory(){
+	
+   delete [] fragmentLength;
+   delete [] geneIndex;
+   delete [] branchMutation;
+   
+	
+}
+
+void FreeCoalescentMemory(){
+	
+   delete [] MRCA;
+   delete [] MRCAtimes;
+
+   delete [] fragments;
+
+   FreeIntPtrMatrix(children,maxN);
+   
+   FreeIntPtrMatrix(parents,maxN);
+
+   delete [] child_flags;
+   delete [] parent_flags;
+
+   delete [] cblocks;
+   delete [] pblocks;
+   
+   delete [] parent_First;
+   delete [] parent_Last;
+
+   delete [] child_First;
+   delete [] child_Last;
+	
+}
+
+void reallocBlock(){
+	
+	pblocks = (int *)realloc(pblocks, (unsigned)(sizeof(int)*(n*(numBLOCKs)*BLOCKSIZE+BUFFER*BLOCKSIZE+INCREMENT*BLOCKSIZE)) );
+	cblocks = (int *)realloc(cblocks, (unsigned)(sizeof(int)*(n*(numBLOCKs)*BLOCKSIZE+BUFFER*BLOCKSIZE+INCREMENT*BLOCKSIZE)) );
+	BUFFER += INCREMENT;	
+	
+}
+
+void setnull(int parent, int position){
+	
+  if (parent_flags[parent] == 0){
+		parent_First[parent]=position / BLOCKSIZE;
+		parent_Last[parent]=position / BLOCKSIZE;	
+		parents[parent][position / BLOCKSIZE] = NULL;  
+  }
+  else if(parent_First[parent]>position / BLOCKSIZE){
+	  for(int i=position / BLOCKSIZE;i<parent_First[parent];i++){
+		  parents[parent][i] = NULL;
+	  }
+	  parent_First[parent]=position / BLOCKSIZE;
+  }
+  else if(parent_Last[parent]<position / BLOCKSIZE){
+	  for(int i=position / BLOCKSIZE;i>parent_Last[parent];i--){
+		  parents[parent][i] = NULL;
+	  }
+	  parent_Last[parent]=position / BLOCKSIZE;
+  }
+  
+}
+
+
+
+void TouchParent(int parent, int position){
+
+  
+  setnull(parent,position);
+           	
+  if (parents[parent][position / BLOCKSIZE] != NULL)return;
+ 
+
+  if (pblockcounter >= n*numBLOCKs+BUFFER){
+	    
+	    int * pblocks_old=pblocks;
+	    int * cblocks_old=cblocks;
+	    
+	    reallocBlock();
+	    
+	    for(int i=0;i<N;i++)
+	      if(parent_flags[i])
+	    	for(int j=parent_First[i];j<=parent_Last[i];j++)
+	    		if(parents[i][j]!=NULL)parents[i][j] += pblocks - pblocks_old;
+	    
+	    for(int i=0;i<N;i++)
+	      if(child_flags[i])
+	    	for(int j=child_First[i];j<=child_Last[i];j++)
+	    		if(children[i][j]!=NULL)children[i][j] += cblocks - cblocks_old;
+	    	
+ 		printf("pblock and cblock enlarged, BUFFER=%d.!\n",BUFFER);
+ 		
+		if(pblocks==NULL || cblocks==NULL)error("Blocks memory reallocation error.\n");
+  }
+
+  parents[parent][position / BLOCKSIZE] = pblocks + BLOCKSIZE * pblockcounter++;
+  
+  for (int i = 0; i < BLOCKSIZE; i++)
+      parents[parent] [position / BLOCKSIZE] [i] = -1;
+      
+  parent_flags[parent] = 1;
+  
+}
+
+void Initialize(){
+	
+   genetimes = new int [poolsize];
+	   
+   MRCA = new int [L];
+   MRCAtimes = new int [L];
+
+   fragments = new int [L];
+
+   children = AllocateIntPtrMatrix(maxN, numBLOCKs);
+   parents = AllocateIntPtrMatrix(maxN, numBLOCKs);
+   
+   parent_First = new int [maxN];
+   parent_Last = new int [maxN];
+
+   child_First = new int [maxN];
+   child_Last = new int [maxN];
+
+   cblocks = (int *)malloc((unsigned)(sizeof(int)*(n*(numBLOCKs)*BLOCKSIZE+BUFFER*BLOCKSIZE)));
+   pblocks = (int *)malloc((unsigned)(sizeof(int)*(n*(numBLOCKs)*BLOCKSIZE+BUFFER*BLOCKSIZE)));	
+  
+   child_flags = new bool [maxN];
+   parent_flags = new bool [maxN];
+   
+   
+   noChange = true;
+	
+   nPOP = (int)(popProfile[0].popsize.size());
+   
+   N = popProfile[0].popStart.back()+popProfile[0].popsize.back();
+   
+   currentProfile = 0;
+	
+   nextFree = 0;
+   
+   for (int i = 0; i < N; i++)
+      child_flags[i] = 0;
+     
+   int count=0;
+   for(int k=0;k<nPOP;k++){    
+	   for (int i = 0; i < nSample[k]; i++)
+	      {
+		  child_First[popProfile[0].popStart[k]+i]=0;
+		  child_Last[popProfile[0].popStart[k]+i]=numBLOCKs-1;
+
+	      child_flags[popProfile[0].popStart[k]+i] = 1;
+
+	      for (int j = 0; j < numBLOCKs; j++)
+         		children[popProfile[0].popStart[k]+i][j] = cblocks + (count * L + j * BLOCKSIZE);
+	      	
+	      for (int j = 0; j < L; j++)
+	         {
+	         genepool[nextFree] = 1;
+	         genetimes[nextFree] = 0;
+	         children[popProfile[0].popStart[k]+i][j / BLOCKSIZE][j % BLOCKSIZE] = nextFree++;
+	         }
+	         
+	         count++;
+	      }
+   }
+
+   for (int i = 0; i < L; i++)
+      fragments[i] = n;
+
+   for (int i = 0; i < N; i++)
+      parent_flags[i] = 0;
+
+   activeSegments = L;
+   
+   //for (int i = 0; i < N; i++)
+   //  for (int j = 0; j < numBLOCKs; j++)
+   //     parents[i][j] = NULL;
+      
+   pblockcounter = 0;
+ 
+}
+
+void InitializeMutation(vector< vector<bool> > &return_chromosome){
+	
+// 	   return_chromosome.clear();
+//    return_chromosome.resize(n);
+//    for(int i=0;i<n;i++){
+//  	   return_chromosome[i].reserve(max(1,(int)(mutation*(double)numChr*(double)L*(double)length*40*(double)N/(double)nPOP)));
+// 	}
+//    
+//    chromosome.resize(n);
+//    for(int i=0;i<n;i++){
+//  		chromosome[i].resize(L);
+// 		for(int j=0;j<L;j++){
+//  		   chromosome[i][j].reserve(max(1,(int)(mutation*(double)length*40*(double)N/(double)nPOP)));
+// 		}
+//    }
+	
+//    cout<<"WGCS-- chromosomes capacity="<<chromosome[0][0].capacity()<<" size="<<chromosome[0][0].size()<<endl;
+   
+   for(int i=0;i<n;i++){
+	   for(int j=0;j<L;j++){
+		   chromosome[i][j].clear();
+	   }
+   }
+   
+//    cout<<"WGCS-- after clear() chromosomes capacity="<<chromosome[0][0].capacity()<<" size="<<chromosome[0][0].size()<<endl;
+   
+   branchMutation = new int [nextFree];
+   for(int i=0;i<nextFree;i++)branchMutation[i]=0;
+
+}
+
+void PrepareMutation(){
+	
+   for(int i=0;i<L;i++)fragmentLength[i]=0;
+   
+   for(int i=0;i<nextFree;i++){
+	   geneIndex[i] = -1; 
+	   branchFragment[i] = -1;
+   }	
+	
+}
+   
+int SampleParent(int individual, int previous){
+	
+	if(fixedPOP){	   
+	   int population = individual/popProfile[0].popsize[0]; 
+	   
+	   if (previous == -1 && nPOP>1){
+	      if (globalRandom.Next() < migration){
+		      int random = globalRandom.NextInt()%(nPOP-1); 
+	          if(population<=random)population=random+1;
+	          else population=random;
+	      }
+	   }
+	      
+	   return globalRandom.NextInt() % (popProfile[0].popsize[0]) + popProfile[0].popStart[population];
+   }
+   else{
+	   
+	   if(noChange){
+		   int population=0; 
+	   
+		   for(int i=0;i<nPOP;i++){
+			   if(individual>=popProfile[currentProfile].popStart[i])population=i;
+			   else break;
+		   }
+		   
+//		   cout<<"No change: ind="<<individual<<" pop="<<population<<endl;
+		   
+		   if (previous == -1 && nPOP>1){
+		      if (globalRandom.Next() < migration){
+			      int random = globalRandom.NextInt()%(nPOP-1); 
+		          if(population<=random)population=random+1;
+		          else population=random;
+		      }
+		   }
+		   
+//		   cout<<"new pop="<<population<<endl;
+		      
+		   return globalRandom.NextInt() % (popProfile[currentProfile].popsize[population]) + popProfile[currentProfile].popStart[population];
+	   }
+	   else{
+		   int population=0; 
+	   
+		   for(int i=0;i<nPOP;i++){
+			   if(individual>=popProfile[currentProfile].popStart[i])population=i; 
+			   else break;
+		   }
+//		   cout<<"Change: ind="<<individual<<" pop="<<population<<endl;
+		   
+		   if (previous == -1 && nPOP>1){
+		      if (globalRandom.Next() < migration){
+			      int random = globalRandom.NextInt()%(nPOP-1); 
+		          if(population<=random)population=random+1;
+		          else population=random;
+		      }
+		   }
+		   
+//		   cout<<"new pop="<<population<<endl;
+		   
+		   if(popProfile[currentProfile].popChange[population].size()==1){
+//			   cout<<"To next pop="<<popProfile[currentProfile].popChange[population][0]<<endl;
+			   
+				return globalRandom.NextInt() % (popProfile[currentProfile+1].popsize[popProfile[currentProfile].popChange[population][0]]) 
+						+ popProfile[currentProfile+1].popStart[popProfile[currentProfile].popChange[population][0]];	   
+		   }
+		   else{
+		   		double selectRandom = globalRandom.Next();
+		   		for(unsigned int i=0;i<popProfile[currentProfile].popChange[population].size();i++){
+			   		if(selectRandom < popProfile[currentProfile].selectProb[population][i]){
+				   		
+//				   		cout<<"To next pop="<<popProfile[currentProfile].popChange[population][i]<<endl;
+				   		
+				   		return globalRandom.NextInt() % (popProfile[currentProfile+1].popsize[popProfile[currentProfile].popChange[population][i]]) 
+						+ popProfile[currentProfile+1].popStart[popProfile[currentProfile].popChange[population][i]];
+			   		}	
+				}
+				cout<<"Random number generator error!"<<endl;
+				exit(1);
+	   	   }
+		   
+		   
+	   }
+	   
+   }
+}
+
+void readPopProfile(string filename, unsigned int numPopulation){
+	   
+	   int popProfile_index=0;
+	   
+	   popProfile.clear();
+	   struct popProfilesStruct temp_popProfile;
+	   popProfile.push_back(temp_popProfile);	
+	   			   
+	   ifstream popFile;
+	   string fileline;
+	   char *term;
+	   char *term1;
+	   
+	   popFile.open(filename.c_str());
+	   if(!popFile)error("Population profile: %s cannot be opened!",filename.c_str());
+	   
+	   getline(popFile,fileline);					// the first line is the population sizes at generation 0
+	   
+	   term = strtok ((char*)fileline.c_str()," \t");  // the first term is the generation
+	   popProfile[popProfile_index].generation=atoi(term);
+	   
+	   if(popProfile[popProfile_index].generation!=0)error("The first line in population profile should be generation 0!");
+	   
+	   term = strtok (NULL, " \t");					   // the second term and after are the sizes of that population
+	   while(term!=NULL){
+			   popProfile[popProfile_index].popsize.push_back(atoi(term));
+			   term = strtok (NULL, " \t");
+	   }
+
+	   popProfile[popProfile_index].popStart.push_back(0); // the first population starts from 0
+	   for(unsigned int i=0;i<popProfile[popProfile_index].popsize.size()-1;i++){
+	   	   popProfile[popProfile_index].popStart.push_back(popProfile[popProfile_index].popsize[i]+popProfile[popProfile_index].popStart.back());
+	   }
+	   		   
+	   maxN = popProfile[popProfile_index].popStart.back()+popProfile[popProfile_index].popsize.back();
+	   
+	   if(popProfile[0].popsize.size()!=numPopulation)error("The number of populations specified by -pop is not consistent with generation 0 in population profile:%s",filename.c_str());
+	   
+	   // read in the population correspondence to the next population profile
+	   
+		   
+	   while(popFile.peek()!=EOF){
+		   getline(popFile,fileline);						// the even number lines are the correspondence of populations of different generation
+		   popProfile[popProfile_index].popChange.resize(popProfile[popProfile_index].popsize.size());
+		   
+		   term = strtok ((char*)fileline.c_str(),"- \t");  // the FROM population
+		   while(term!=NULL){
+			   
+			   term1 = strtok (NULL, "- \t");				// the TO population
+
+			   popProfile[popProfile_index].popChange[atoi(term)-1].push_back(atoi(term1)-1);
+			   
+			   term = strtok (NULL, "- \t");
+		   }
+		   
+		   popProfile_index++;
+		   popProfile.push_back(temp_popProfile);
+		   
+		   if(popFile.peek()==EOF)error("Error in population profile: the last line should specify the population sizes.");
+		   
+		   // read the next population profile
+		   
+	   	   getline(popFile,fileline);					// the population sizes
+   
+		   term = strtok ((char*)fileline.c_str()," \t");  // the first term is the generation
+		   popProfile[popProfile_index].generation=atoi(term);
+		   
+		   term = strtok (NULL, " \t");					   // the second term and after are the sizes of that population
+		   while(term!=NULL){
+				   popProfile[popProfile_index].popsize.push_back(atoi(term));
+				   term = strtok (NULL, " \t");
+		   }
+		   
+		   popProfile[popProfile_index].popStart.push_back(0); // the first population starts from 0
+		   for(unsigned int i=0;i<popProfile[popProfile_index].popsize.size()-1;i++){
+		   	   popProfile[popProfile_index].popStart.push_back(popProfile[popProfile_index].popsize[i]+popProfile[popProfile_index].popStart.back());
+		   }
+		   
+		   if(maxN < popProfile[popProfile_index].popStart.back()+popProfile[popProfile_index].popsize.back())
+		   			maxN = popProfile[popProfile_index].popStart.back()+popProfile[popProfile_index].popsize.back();
+	   }	
+	   
+	   for(unsigned int i=0;i<popProfile.size();i++){
+		   popProfile[i].selectProb.resize(popProfile[i].popChange.size());
+		   for(unsigned int j=0;j<popProfile[i].popChange.size();j++){
+			   float sum=0.0;
+			   for(unsigned int k=0;k<popProfile[i].popChange[j].size();k++)sum += popProfile[i+1].popsize[popProfile[i].popChange[j][k]];
+			   
+			   popProfile[i].selectProb[j].push_back(((float)popProfile[i+1].popsize[popProfile[i].popChange[j][0]])/sum);
+			   for(unsigned int k=1;k<popProfile[i].popChange[j].size();k++)
+			   		popProfile[i].selectProb[j].push_back(popProfile[i].selectProb[j].back()+popProfile[i+1].popsize[popProfile[i].popChange[j][k]]/sum);
+		   }
+	   }
+	   // check for consistence
+	   
+	   int min, max;
+	   
+	   for(unsigned int i=0;i<popProfile.size()-1;i++){
+			min=99999;max=-1;
+			for(unsigned int j=0;j<popProfile[i].popChange.size();j++){
+				if(popProfile[i].popChange[j].size()==0)error("One population does not have its parent population: profile file line %d, population %d\n",i+1,j+1);	
+				for(unsigned int k=0;k<popProfile[i].popChange[j].size();k++){
+					if(min>popProfile[i].popChange[j][k])min=popProfile[i].popChange[j][k];
+					if(max<popProfile[i].popChange[j][k])max=popProfile[i].popChange[j][k];
+				}
+			}   
+			if(min!=0 || max!=(int)(popProfile[i+1].popsize.size()-1))
+				error("The population changes does not consistent with the next population profile. line=%d, min=%d, max=%d, num POP in next line=%d\n",i+1,min+1,max+1,popProfile[i+1].popsize.size());
+	   }
+	
+}
+
+
+void readRecombination(string filename, int pieces){
+
+	   ifstream recombFile;
+	   string fileline;
+	   char *term;
+		
+	   vector<float> recDistribution;
+	   vector<float> recRate;
+	   
+	   float sum=0.0;
+	   float random;
+	   
+	   recombFile.open(filename.c_str());
+	   if(!recombFile)error("Recombination profile: %s cannot be opened!",filename.c_str());
+	   
+	   getline(recombFile,fileline);		// the first line is head of distribution table
+	   term = strtok ((char*)fileline.c_str()," \t");
+	   term = strtok (NULL, " \t");	
+	   
+	   if(strcmp(term,"Freq")==0 || strcmp(term,"freq")==0){
+		   while(recombFile.peek()!=EOF){
+			   getline(recombFile,fileline);
+			   
+			   term = strtok ((char*)fileline.c_str()," \t"); // the first number is the recombination rate
+			   recRate.push_back((float)(atof(term)));
+			   
+			   if(atof(term)<0 || atof(term)>1)error("Recombination rate should be between 0 and 1. input=%f",atof(term));
+			   
+			   term = strtok (NULL, " \t");					  // the second number is the frequency of recombination rate
+			   recDistribution.push_back((float)(atof(term)));
+			   
+			   if(atof(term)<=0)error("Frequency should be positive. input=%f",atof(term));
+			   
+			   sum += (float)(atof(term));
+		   }
+		   
+		   if(recRate.size()!=recDistribution.size() || recRate.size()==0)error("Errors in the recombination profile: %s",filename.c_str());
+		   
+		   // normalize the frequency 
+		   for(unsigned int i=0;i<recDistribution.size();i++) recDistribution[i] /= sum;
+		   
+		   // compute the cumulative distribution
+		   for(unsigned int i=1;i<recDistribution.size();i++) recDistribution[i] += recDistribution[i-1];
+		   
+		   recombVector = new float [pieces-1];
+		   
+		   for(int i=0;i<pieces-1;i++){
+			   random = (float)(globalRandom.Next());
+			   for(unsigned int j=0;j<recDistribution.size();j++){
+					if(random < recDistribution[j]){
+						recombVector[i]=recRate[j];
+						break;   
+					}
+			   }
+		   }
+		   
+			//check the input
+			cout<<"*** RECOMBINATION RATES PROFILE ***"<<endl;
+			cout<<"Recombination_Rate\tCumulative_frequency"<<endl;
+			for(unsigned int i=0;i<recRate.size();i++){
+				cout<<recRate[i]<<"\t\t\t"<<recDistribution[i]<<endl;  
+			}		
+			   
+   	   }
+   	   else if(strcmp(term,"Pos")==0 || strcmp(term,"pos")==0){
+	   	   float currentRec=0.0,nextRec=0.0;
+	   	   int nextPos=-1;
+	   	   
+	   	   if(recombFile.peek()!=EOF){
+		   	   getline(recombFile,fileline);
+		   	   term = strtok ((char*)fileline.c_str()," \t");
+	   		   currentRec=(float)(atof(term));
+	   		   
+	   		   if(recombFile.peek()!=EOF){
+		   		   getline(recombFile,fileline);
+		   	   	   term = strtok ((char*)fileline.c_str()," \t");
+	   		       nextRec=(float)(atof(term));
+	   		       term = strtok (NULL, " \t");	
+	   		       nextPos=atoi(term);
+	   		   }
+	   		   else{
+		   		   nextPos=-1;
+	   		   }
+	   		   
+   		   }
+   		   else{
+	   		   error("Recombination file does not contain recombination rate information.\n");
+   		   }
+   		   
+		   recombVector = new float [pieces-1];
+   		   
+   		   for(int i=0;i<pieces-1;i++){
+	   		   if(i+1==nextPos){
+		   		    currentRec=nextRec;
+		   		    if(recombFile.peek()!=EOF){
+		   		   		getline(recombFile,fileline);
+		   	   	   		term = strtok ((char*)fileline.c_str()," \t");
+	   		       		nextRec=(float)(atof(term));
+	   		       		term = strtok (NULL, " \t");	
+	   		       		nextPos=atoi(term);
+	   		   		}
+	   		   		else{
+		   		   		nextPos=-1;
+	   		   		}
+		   		}
+		   		
+		   		recombVector[i]=currentRec;
+		   }
+   		   
+   		   
+   		   
+   	   }
+   	   else{
+	   	   error("The second column is not \"freq\" or \"pos\".\n");
+   	   }
+   	   
+   	   
+
+	   
+	   cout<<"The recombination rates between fragments are:"<<endl<<"\t";
+	   for(int i=0;i<pieces-1;i++)cout<<"("<<i+1<<") "<<recombVector[i]<<" ";
+	   cout<<"("<<pieces<<")"<<endl;
+	   cout<<"*** END OF RECOMBINATION RATES PROFILE ***"<<endl;
+	   	
+}
+
+void drawtrees(){
+	
+	map<int,string > tree;
+	ostringstream  temp;
+	
+			   
+	for(int i=0;i<L;i++){
+		tree.clear();
+		
+		for(int j=0;j<n;j++){
+			if(tree[genepool[j*L+i]]!=""){
+				temp.str("");
+				temp <<tree[genepool[j*L+i]]<<","<<j+1<<":"<<genetimes[genepool[j*L+i]]-genetimes[j*L+i];
+				tree[genepool[j*L+i]]=	temp.str();
+			}
+			else{
+				temp.str("");
+				temp <<j+1<<":"<<genetimes[genepool[j*L+i]]-genetimes[j*L+i];
+				tree[genepool[j*L+i]]=	temp.str();
+			}
+		}
+		
+	
+		for(int j=n*L;j<nextFree;j++){
+			if(tree[j]!=""){
+				if(genepool[j]==0){
+					cout<<"The tree for fragment "<<i+1<<"= ("<<tree[j]<<");"<<endl;
+				}
+				else{							
+					if(tree[genepool[j]]!=""){
+						temp.str("");
+						temp <<tree[genepool[j]]<<",("<<tree[j]<<"):"<<genetimes[genepool[j]]-genetimes[j];
+						tree[genepool[j]]=	temp.str();
+					}
+					else{
+						temp.str("");
+						temp <<"("<<tree[j]<<"):"<<genetimes[genepool[j]]-genetimes[j];
+						tree[genepool[j]]=temp.str();
+					}
+				}
+			}
+		}
+	}
+	
+	tree.clear();
+}
+
+
+
+
+void genome(string popSize, 						// effective population size of each population
+		  int nSubPOP, 							// number of populations
+		  vector<int> & nSubSample, 			// number of samples (haplotypes) draw from each populations
+		  int numPieces, 						// number of fragments for each sample (chromosome)
+		  int pieceLen, 						// length in base pair of each fragment
+		  int numIndepRegion, 					// number of independent regions (independent chromosome)
+		  int s,								// fixed number of SNPs want to simulate, randomly place s SNPs on the genealogy, -1 means the number of mutation follows a Poisson distribution 
+		  string rec,							// recombination rate between consecutive fragments per generation
+		  double mut, 							// mutation rate per generation per base pair
+		  double mig, 							// migration rate per generation
+		  vector< vector<bool> > &return_chromosome, // the return chromosome by connecting all independent chromosome into one long chromosome
+		  long t,								// random seed 
+		  int drawtree,						// drawtree=1 output the genealogy trees for each fragment and chromosome; drawtree=0 do not output the trees	
+		  bool printparameters		   	 		// =1 print out all input parameters, =0 do not print input parameters
+		  )
+{   
+   if(atoi(popSize.c_str())>0){            	// if the population size is fixed during evoluation
+	   N = atoi(popSize.c_str())*nSubPOP;
+	   maxN = N;
+	   
+	   fixedPOP = true;
+	   
+	   popProfile.clear();
+	   struct popProfilesStruct temp_popProfile;
+	   popProfile.push_back(temp_popProfile);
+	   
+	   popProfile[0].generation = 0;
+	   
+	   for(int i=0;i<nSubPOP;i++)popProfile[0].popsize.push_back(atoi(popSize.c_str()));
+	   for(int i=0;i<nSubPOP;i++)popProfile[0].popStart.push_back(i*atoi(popSize.c_str()));
+	   
+	   cout<<"*** POPULATION PROFILE ***"<<endl;
+	   cout<<"The size of each population is fixed to "<<atoi(popSize.c_str())<<endl;
+	   cout<<"Sizes of populations = ";
+	   for(unsigned int j=0;j<popProfile[0].popsize.size();j++)cout<<popProfile[0].popsize[j]<<" ";
+	   cout<<endl;
+// 	   cout<<"POP Start index = ";
+// 	   for(int j=0;j<popProfile[0].popStart.size();j++)cout<<popProfile[0].popStart[j]<<" ";
+// 	   cout<<endl;
+	   cout<<"*** END OF POPULATION PROFILE ***"<<endl<<endl;
+			
+   }
+   else{									// popSize stores the filename of the population profile during evoluation
+ 	   fixedPOP = false;	  
+   	
+	   currentProfile = 0;
+	   
+	   readPopProfile(popSize,nSubPOP);
+	   
+	   N = popProfile[0].popStart.back()+popProfile[0].popsize.back();
+	   
+	   // print out the population profile information.
+	   
+	   cout<<"*** POPULATION PROFILE ***"<<endl;
+	   cout<<"Maximum total size of populations at all generations: N = "<<maxN<<endl;
+	   cout<<"Total size of populations at generation 0: N = "<<N<<endl<<endl;
+	   
+	   for(unsigned int i=0;i<popProfile.size();i++){
+			cout<<"Generation = "<<popProfile[i].generation<<endl;
+			cout<<"Sizes of populations= ";
+			for(unsigned int j=0;j<popProfile[i].popsize.size();j++)cout<<popProfile[i].popsize[j]<<" ";
+			cout<<endl;
+// 			cout<<"POP Start index = ";
+// 			for(int j=0;j<popProfile[i].popStart.size();j++)cout<<popProfile[i].popStart[j]<<" ";
+// 			cout<<endl;
+			if(popProfile[i].popChange.size()>0){
+				cout<<"Populations change backward in time:"<<endl;
+				for(unsigned int j=0;j<popProfile[i].popChange.size();j++){
+					cout<<" From "<<j+1<<" to ";
+					for(unsigned int k=0;k<popProfile[i].popChange[j].size();k++)cout<<popProfile[i].popChange[j][k]+1<<" ";
+					cout<<endl;
+				}
+				
+// 				cout<<"Selection cumulative probability:"<<endl;   
+// 				for(int j=0;j<popProfile[i].selectProb.size();j++){
+// 					cout<<" From "<<j+1<<" with prob = ";
+// 					for(int k=0;k<popProfile[i].selectProb[j].size();k++)cout<<popProfile[i].selectProb[j][k]<<" ";
+// 					cout<<endl;
+// 				}
+			}
+			if(i==popProfile.size()-1)cout<<"*** END OF POPULATION PROFILE ***"<<endl;
+			cout<<endl;
+	   }
+   }
+   
+   n=0;
+   nPOP=nSubPOP;
+   nSample.resize(nPOP);
+   for(int i=0;i<nPOP;i++){
+	   if(popProfile[0].popsize[i]<nSubSample[i])error("Population size should be larger than sample size!");
+	   n += nSubSample[i];
+	   nSample[i]=nSubSample[i];
+   }
+   
+   L = numPieces;
+   length = pieceLen;
+   numChr = numIndepRegion;
+   SNP=s;
+   
+   if(atof(rec.c_str())>0){				// if the recombination rate is fixed
+	   recombination = atof(rec.c_str());
+   }
+   else{						// the distribution of recombination rate is described in the file (filename stored in rec)
+   	   recombination = -1;
+	   readRecombination(rec,L);
+   }
+   
+   
+   
+   mutation = mut;
+   migration = mig;
+   
+   if(t<=0)t=(long)(time(0)); // if t >0 we set the seed to t
+   cout<<endl<<"random seed="<<t<<endl;
+   globalRandom.Reset(t);	
+   
+   if(printparameters){
+	   printf("popSize=%s, n=%d, nPOP=%d, L=%d, length=%d, numChr=%d, SNP=%d, rec=%s, mut=%.4e, mig=%.4e, drawtree=%d, printparameters=%d\n",popSize.c_str(),n,nPOP,L,length,numChr,SNP,rec.c_str(),mut,mig,drawtree,printparameters);
+	   cout<<"subSample = ";
+	   for(int i=0;i<nPOP;i++){
+		   cout<<nSample[i]<<" ";
+	   }
+	   cout<<endl;
+   }
+   
+   if(N<=0 || n<=0 || L<=0 || length <=0 || numChr <=0 || atof(rec.c_str()) <0 || mut <0 || mig<0){
+	   cout<<"Input parameters have to be positive."<<endl;
+	   printf("N=%d, n=%d, L=%d, length=%d, numChr=%d, recombination=%s, mutation=%lf, migration=%lf\n",N,n,L,length,numChr,rec.c_str(),mut,mig);
+	   exit(0);
+   }
+
+   poolsize = 2 * n * L - L;
+
+   AllocateMemory(return_chromosome);
+
+//    cout<<"memory OK!"<<endl<<endl;
+   
+	for(int chr=0;chr<numChr;chr++){
+		
+	   Initialize();
+	   
+// 	   cout<<"Initialization OK!"<<endl;
+	
+	   bool done = false;
+	   int  generation = 0, units = n * L;
+	   
+	   while (!done)
+	      {
+	//       printf("\n\nGeneration %d -- Tracking %d units of sequence [pool used %.2f%%]   %c",
+	//              generation++, units, nextFree * 100. / poolsize,
+	//              !debug ? '\r' : '\n');
+	
+		  generation++;
+		  
+		  if(!fixedPOP && currentProfile<popProfile.size()-1)
+				if(generation == popProfile[currentProfile+1].generation)noChange=false;
+		  
+	      //if (debug)
+	      //   {
+	      //   printf("F  : ");
+	      //   for (int i = 0; i < L; i++)
+	      //      printf("%3d ", fragments[i]);
+	      //   printf("\n");
+	      //   }
+	
+	      // Position of the first segment allocated to the current generation
+	      int generationStart = nextFree;
+	
+	      // For each individual in the population
+	      for (int i = 0; i < N && !done; i++)
+	         if (child_flags[i])
+	            {
+		            
+	            /*if (debug) printf("%2d : ", i);*/
+	
+	            // The position we are currently tracking and its distance
+	            // from the previous position
+	            int distance = 0, parent = -1;
+	            int startPiece = 0;
+	            float recombProb;
+           
+			    for (int block = child_First[i]; block <= child_Last[i]; block++)
+                if (children[i][block] == NULL)
+                   distance += BLOCKSIZE;
+                else
+                   {
+                   int pos = BLOCKSIZE * block;
+                   int pos_in_block;
+				   int pos_bound = (L < BLOCKSIZE * (block + 1) ? L : BLOCKSIZE * (block + 1) );
+
+                   while (pos <  pos_bound )
+                       {
+	                       
+	                   pos_in_block = pos % BLOCKSIZE;
+                       // Skip through uninformative positions along each chromosome
+                       if (children[i][block][pos_in_block] == -1)
+                          {
+                          //if (debug) printf("  . \n");
+                          distance++;
+                          pos++;
+                          continue;
+                          }
+                         
+                       // Pick an ancestor at random for the current piece of chromosome
+                       
+                       if(recombination>0){			// if the recombination rate is fixed over the genome
+	                       recombProb = (float)(pow(1.0 - recombination, distance));
+                       }
+                       else{						// if the recombination rate is specified by the profile
+	                       recombProb = 1;
+	                       for(int recV_iter=startPiece;recV_iter<pos;recV_iter++)recombProb *= (1-recombVector[recV_iter]);
+                       }
+                       
+                    
+                       if (parent == -1 ||
+                           distance > 0 && globalRandom.Next() > recombProb){
+	                            parent = SampleParent(i, parent);
+                        }
+                            
+                       
+                       TouchParent(parent, pos);
+                       
+                       
+                       // Reset distance between segments
+                       distance = 1;
+                       startPiece=pos;
+
+                       //if (debug) printf("parent=%3d pos=%d address=%d address2=%d\n", parent,pos,parents[parent][block],parents[parent][pos/BLOCKSIZE]);
+                       
+                       // Check if we sampled a new ancestral segment ...
+                       if (parents[parent][block][pos_in_block] == -1)
+                          {
+                          // We delay sampling a new gene from the pool until we know
+                          // a coalescent event occured ...
+                          parents[parent][block][pos_in_block] = children[i][block][pos_in_block];
+                          
+                          }
+                       // Or if a coalescent event occured
+                       else
+                          {
+                          // If we previously delayed sampling this parent
+                          if (parents[parent][block][pos_in_block] < generationStart)
+                            {
+	                            genetimes[nextFree] = generation;
+	                            genepool[parents[parent][block][pos_in_block]] = nextFree;
+	                            parents[parent][block][pos_in_block] = nextFree++;
+
+	                            if (nextFree == poolsize && units != 2)
+	                               error("Genepool exhausted\n");
+                            }
+
+                          genepool[children[i][block][pos_in_block]] = parents[parent][block][pos_in_block];
+                          
+                          fragments[pos]--;
+                          units--;
+
+                         if (fragments[pos] == 1)
+                            {
+                            // Reached the MRCA for this fragment
+                            MRCAtimes[pos] = generation;
+                            MRCA[pos] = parents[parent][block][pos_in_block];
+                            
+                            genepool[parents[parent][block][pos_in_block]]=0;
+                            
+                            parents[parent][block][pos_in_block] = -1;
+                            units--;
+         
+                            // One less segment to track
+                            if (--activeSegments == 0)
+                               {
+                               done = true;
+                               break;
+                               }
+                            }
+                         }
+                      pos++;
+                      }
+                   }
+	            }
+	            
+	      NewGeneration();
+	    
+	      }
+	
+	   if(drawtree){
+		   cout<<endl<<"Genealogy trees for chromosome:"<<chr+1<<endl;
+		   drawtrees();
+	   }  
+	      
+	   // assign mutations
+
+	   cout<<"Simulate mutation for Chromosome "<<chr+1<<endl;
+	   
+	   FreeCoalescentMemory();
+	   
+	   cout<<"Coalescent Memory free"<<endl;
+	   
+	   InitializeMutation(return_chromosome);
+
+	   cout<<"Mutation process initialized"<<endl;
+	    	
+// 	    cout<<"genepool:\t";
+// 		for(int i=0;i<2*n*L-L;i++)cout<<genepool[i]<<" ";
+// 		cout<<endl;
+// 		
+// 		cout<<"genetimes:\t";
+// 		for(int i=0;i<2*n*L-L;i++)cout<<genetimes[i]<<" ";
+// 		cout<<endl;
+// 	   
+// 		cout<<"nextFree="<<nextFree<<endl;
+		
+// 		cout<<"mutation:\t";
+// 		for(int i=0;i<nextFree;i++)cout<<branchMutation[i]<<" ";
+// 		cout<<endl;
+		
+	   if(SNP<=0){
+		   for(int i=0;i<nextFree;i++){
+				if(genepool[i]!=0)branchMutation[i]=poissonInt((double)length*(double)(genetimes[genepool[i]]-genetimes[i])*mutation);   
+		   }
+	   }
+	   else{
+		   unsigned long * cumulativeCount = new unsigned long [nextFree];
+		   
+		   if(genepool[0]!=0)cumulativeCount[0] = genetimes[genepool[0]]-genetimes[0];
+		   else cumulativeCount[0] = 0;
+		   
+		   for(int i=1;i<nextFree;i++){
+				if(genepool[i]!=0)cumulativeCount[i] = cumulativeCount[i-1] + genetimes[genepool[i]]-genetimes[i]; 
+				else cumulativeCount[i] = cumulativeCount[i-1];  
+		   }
+		   
+		   unsigned long totalGeneration = cumulativeCount[nextFree - 1];
+		   
+// 		   cout<<"Total Generation = "<<totalGeneration<<endl;
+// 		   
+// 		   if((unsigned long)SNP>totalGeneration){
+// 				cout<<"Required number of SNPs("<<SNP<<") is more than total number of generations("<<totalGeneration<<")!"<<endl;   
+// 		   }
+		   
+		   for(int i=0;i<SNP;i++){
+			   
+			   unsigned long randomLong = globalRandom.NextInt()%totalGeneration;
+			   
+			   int begin = 0;
+			   int end = nextFree - 1;
+			   int middle = (begin+end)/2;
+			   
+			   while(1){
+				   if(cumulativeCount[middle]<randomLong){
+						begin = middle;
+				   }
+				   else{
+						end = middle;
+				   }
+				   
+				   if(end == begin+1)break;
+				   
+				   middle = (begin+end)/2;
+			   }
+			   
+			   if(cumulativeCount[end] < randomLong || (cumulativeCount[begin] >= randomLong && begin != 0)){
+					cout<<"searching error: begin="<<begin<<" end="<<end<<" random="<<randomLong<<" c1="<<cumulativeCount[begin]<<" c2="<<cumulativeCount[end]<<endl;   
+					error("see above message!\n");
+			   }
+			   
+			   int sampleIndex;
+			   
+			   if(cumulativeCount[begin]>=randomLong)sampleIndex = begin;
+			   else sampleIndex = end;
+			   
+			   if(genepool[sampleIndex]==0){
+				   cout<<"genepool error: begin="<<begin<<" end="<<end<<" random="<<randomLong<<" c1="<<cumulativeCount[begin]<<" c2="<<cumulativeCount[end]<<endl;   
+				   error("see above message!\n");
+			   }
+			   
+			   branchMutation[sampleIndex]++;
+		   }
+		   
+		   delete [] cumulativeCount;
+		   
+	   }
+   
+	   delete [] genetimes;
+	   
+	   cout<<"Mutation assigned"<<endl;
+	    
+//  	cout<<"mutation:\t";
+// 		for(int i=0;i<nextFree;i++)cout<<branchMutation[i]<<" ";
+// 		cout<<endl;
+
+	   AllocateMutationMemory();
+	   
+	   cout<<"Mutation Memory allocated"<<endl;
+	   	
+	   PrepareMutation();
+	
+	   for(int i=0;i<n;i++){
+	   		for(int j=0;j<L;j++){
+		   		if(branchMutation[i*L+j]>0){
+			   		geneIndex[i*L+j] = fragmentLength[j];
+			   		fragmentLength[j] += branchMutation[i*L+j];
+		   		}
+		   		branchFragment[i*L+j] = j;
+	   		}
+	   }   
+	      
+	   
+	   for(int i=0;i<nextFree;i++){
+			if(genepool[genepool[i]]!=0 && genepool[i]!=0){
+			   if(branchMutation[genepool[i]]>0 && branchFragment[genepool[i]]== -1 ){
+				   geneIndex[genepool[i]] = fragmentLength[branchFragment[i]];
+				   fragmentLength[branchFragment[i]] += branchMutation[genepool[i]];
+			   }
+			   branchFragment[genepool[i]]=branchFragment[i];
+			}
+	   }
+	   
+	   //cout<<"branchfragment done"<<endl;
+	   
+//    	cout<<"branchFragment:\t";
+// 		for(int i=0;i<nextFree;i++)cout<<branchFragment[i]<<" ";
+// 		cout<<endl;
+
+//    	cout<<"geneIndex:\t";
+// 		for(int i=0;i<nextFree;i++)cout<<geneIndex[i]<<" ";
+// 		cout<<endl;
+// 	
+// 		cout<<"fragmentLength:\t";
+// 		for(int i=0;i<L;i++)cout<<fragmentLength[i]<<" ";
+// 		cout<<endl;
+		
+
+// 	    cout<<"genepool genetime mutation frag index\n";
+// 		for(int i=0;i<nextFree;i++)cout<<genepool[i]<<" "<<genetimes[i]<<" "<<branchMutation[i]<<" "<<branchFragment[i]<<" "<<geneIndex[i]<<endl;
+// 		cout<<endl;
+	   
+
+   
+       delete [] branchFragment;
+       
+	   cout<<"Internal memory free"<<endl;
+	   
+	   // print the fragment index for each SNP
+	   cout<<"SNP genetic position scaled to [0,1]:"<<endl;
+	   if(L>1){
+		   for(int i=0;i<L;i++){
+				for(int j=0;j<fragmentLength[i];j++)cout<<(float)i/(float)(L-1)<<" ";   
+		   }
+	   }
+	   else{
+		   cout<<"Only one fragment to be simulated. No recombination between SNPs.";
+	   }
+	   cout<<endl;
+	   	   
+	   int totalLength=0;
+	   if(SNP <=0){	
+		   for(int i=0;i<L;i++)totalLength += fragmentLength[i];
+	   }
+	   
+	   
+	   for(int i=0;i<n;i++){
+	   		for(int j=0;j<L;j++){
+	   			chromosome[i][j].resize(fragmentLength[j],0);
+
+	   			for(int k=0;k<branchMutation[i*L+j];k++)chromosome[i][j][geneIndex[i*L+j]+k] = 1;
+		   		
+		   		int current=i*L+j;
+  		
+		   		while(genepool[genepool[current]]!=0){
+			   		for(int k=0;k<branchMutation[genepool[current]];k++)chromosome[i][j][geneIndex[genepool[current]]+k]=1;
+			   		current=genepool[current];
+		   		}
+	   		}
+	   }
+	   
+	   cout<<"Chromosome done"<<endl;
+	   
+	   FreeMutationMemory();
+	   
+	   if(SNP <=0){	
+		   for(int i=0;i<n;i++){
+				return_chromosome[i].reserve(return_chromosome[i].size()+totalLength);
+		   }
+	   }
+	   
+	   for(int i=0;i<n;i++){
+			for(int j=0;j<L;j++){
+				return_chromosome[i].insert(return_chromosome[i].end(),chromosome[i][j].begin(),chromosome[i][j].end());
+			}
+	   }	 
+	   
+	   cout<<"Return chromosome done"<<endl;
+	   
+	   // print the chromosomes
+// 	   for(int i=0;i<n;i++){
+// 			printf("Chr %d: ",i);
+// 			for(int j=0;j<L;j++){
+// 				for(int k=0;k<chromosome[i][j].size();k++)cout<<chromosome[i][j][k];
+// 				cout<<" ";	
+// 			}
+// 			cout<<endl;
+// 	   }	   
+// 	   
+// 	   
+// 	  for(int i=0;i<n;i++){
+// 			printf("return Chr %d: ",i);
+// 			for(int j=0;j<return_chromosome[i].size();j++){
+// 				cout<<return_chromosome[i][j];
+// 			}
+// 			cout<<endl;
+// 	   }
+
+	   
+   
+   }
+   
+   FreeMemory();
+   
+   printf("\nDone simulating ARG ...\n");
+   
+   
+   
+}
+
+
diff --git a/genome.h b/genome.h
new file mode 100644
index 0000000..ec8f113
--- /dev/null
+++ b/genome.h
@@ -0,0 +1,62 @@
+////////////////////////////////////////////////////////////////////// 
+// genome.h 
+//////////////////////////////////////////////////////////////////////////////
+//              COPYRIGHT NOTICE FOR GENOME CODE
+//
+// Copyright (C) 2006 - 2009, Liming Liang and Goncalo Abecasis,
+// 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.
+//
+///////////////////////////////////////////////////////////////////////////////
+
+
+#include <vector>
+#include <string>
+
+using namespace std;
+
+void genome(string popSize, 						// effective population size of each population
+		  int nSubPOP, 							// number of populations
+		  vector<int> & nSubSample, 			// number of samples (haplotypes) draw from each populations
+		  int numPieces, 						// number of fragments for each sample (chromosome)
+		  int pieceLen, 						// length in base pair of each fragment
+		  int numIndepRegion, 					// number of independent regions (independent chromosome)
+		  int s,								// fixed number of SNPs want to simulate, randomly place s SNPs on the genealogy 
+		  string rec,							// recombination rate between consecutive fragments per generation
+		  double mut, 							// mutation rate per generation per base pair
+		  double mig, 							// migration rate per generation
+		  vector< vector<bool> > &return_chromosome, // the return chromosome by connecting all independent chromosome into one long chromosome
+		  long t,								// random seed 
+		  int drawtree,						// drawtree=1 output the genealogy trees for each fragment and chromosome; drawtree=0 do not output the trees	
+		  bool printparameters =false 	 		// =1 print out all input parameters, =0 do not print input parameters
+		  );
+
+
+
+ 
diff --git a/makefile b/makefile
new file mode 100644
index 0000000..bfe2dd4
--- /dev/null
+++ b/makefile
@@ -0,0 +1,27 @@
+# module name
+NAME = genome
+
+# switches
+SW = -O2 -Wall
+
+# libreries
+LIB = -lm
+
+# compiler
+CC = g++
+
+# main module 
+
+OBJ = Main.o genome.o  Random.o Error.o stochastic.o
+
+$(NAME): $(OBJ)
+	$(CC) -o $(NAME) $(OBJ) $(LIB) $(SW)
+
+.cpp.o : 
+	$(CC) $(SW) -o $@ -c $*.cpp 
+
+.SUFFIXES : .cpp .c .o $(SUFFIXES)
+
+clean : 
+	rm -rf *.o core*
+
diff --git a/makefile-32bitsys b/makefile-32bitsys
new file mode 100644
index 0000000..7626bae
--- /dev/null
+++ b/makefile-32bitsys
@@ -0,0 +1,27 @@
+# module name
+NAME = genome-32bit
+
+# switches
+SW = -O2 -Wall -m32
+
+# libreries
+LIB = -lm
+
+# compiler
+CC = g++
+
+# main module 
+
+OBJ = Main.o genome.o  Random.o Error.o stochastic.o
+
+$(NAME): $(OBJ)
+	$(CC) -o $(NAME) $(OBJ) $(LIB) $(SW)
+
+.cpp.o : 
+	$(CC) $(SW) -o $@ -c $*.cpp 
+
+.SUFFIXES : .cpp .c .o $(SUFFIXES)
+
+clean : 
+	rm -rf *.o core*
+
diff --git a/population.txt b/population.txt
new file mode 100644
index 0000000..e162178
--- /dev/null
+++ b/population.txt
@@ -0,0 +1,7 @@
+0	500	500	500
+1-1 2-1 3-2
+10	700	300
+1-1 1-2 2-3
+20	500	200	300
+1-1 2-1 3-1
+30	500
\ No newline at end of file
diff --git a/readme.txt b/readme.txt
new file mode 100644
index 0000000..1b62610
--- /dev/null
+++ b/readme.txt
@@ -0,0 +1,84 @@
+
+GENOME: coalescent-based whole genome simulator
+(c) 2006-2009 Liming Liang, Goncalo Abecasis
+
+
+
+DISCLAIMER
+==========
+
+This is version of GENOME is provided to you free of charge.  It comes with
+no guarantees. If you report bugs and provide helpful comments, the next
+version will be better.
+
+
+DISTRIBUTION
+============
+
+The latest version of GENOME is available at:
+
+   http://www.sph.umich.edu/csg/liang/genome
+
+If you use GENOME please register by e-mailing lianglim@umich.edu. 
+
+Code for GENOME is subject to specific
+license conditions -- please see file LICENSE.GENOME for additional
+information. 
+
+Code for the Mersenne Twister random number generator is subject to specific
+license conditions -- please see file LICENSE.twister for additional
+information.
+
+
+
+
+BRIEF INSTRUCTIONS
+==================
+
+A complete manual is available at:
+
+     http://www.sph.umich.edu/csg/liang/genome
+
+
+Parameters and default values are listed below:
+
+      -pop       number of subpopulations and size of subsamples [ 2 10 10 ] 
+      -N         effective size of each subpopulation or the filename for population profile [10000] 
+      -c         number of independent regions (chromosomes) to simulate [1] 
+      -pieces    number of fragments per independent region [100] 
+      -len       length in base of each fragment [10000] 
+      -s         fixed number of SNPs per independent region (chromosome), -1 = number of SNPs follows Poisson distribution [-1] 
+      -rec       recombination rate bewteen consecutive fragments or the filename for recombination rate distribution [0.0001] 
+      -mut       mutation rate per generation per base pair [1e-08] 
+      -mig       migration rate per generation per individual [0.00025] 
+      -seed      random seed, -1 = use time as the seed [-1] 
+      -tree      1=draw the genealogy trees, 0=do not output [0] 
+      -maf       Output SNPs with minor allele frequence greater than [0] 
+      -prop      To keep this proportion of SNPs with MAF < the value of -maf parameter [1] 
+
+
+
+An example command line is:
+
+     genome -pop 2 3 5 -N 100 -tree 1
+
+This command will simulate 3 and 5 sequences from two populations respectively.
+Each population is of size 100. Other parameters will take the default values listed in "[]".
+The genealogy trees in Newick format will be output.
+
+recombination.txt contains an example of the recombination rates profile.
+population.txt contains an example of the population history profile. 
+For detail of these two files please refer to the -N and -rec parameters in the manual.
+
+
+REFERENCE
+=========
+Liang L., Zollner S., Abecasis G.R. (2007) GENOME: a rapid coalescent-based whole genome simulator.
+Bioinformatics 23(12):1565-7
+
+
+
+Comments and suggestions are welcome!
+
+Liming Liang
+lianglim@umich.edu
diff --git a/recombination-pos.txt b/recombination-pos.txt
new file mode 100644
index 0000000..921a858
--- /dev/null
+++ b/recombination-pos.txt
@@ -0,0 +1,6 @@
+Rec		Pos
+0.0001	1
+0.001		3
+0.5	5
+0.7	6
+0.9	7
diff --git a/recombination.txt b/recombination.txt
new file mode 100644
index 0000000..5fdc153
--- /dev/null
+++ b/recombination.txt
@@ -0,0 +1,3 @@
+Rec		Freq
+0.0001	9
+0.001		1
diff --git a/stochastic.cpp b/stochastic.cpp
new file mode 100644
index 0000000..61eea36
--- /dev/null
+++ b/stochastic.cpp
@@ -0,0 +1,115 @@
+////////////////////////////////////////////////////////////////////// 
+// stochastic.cpp 
+//////////////////////////////////////////////////////////////////////////////
+//              COPYRIGHT NOTICE FOR GENOME CODE
+//
+// Copyright (C) 2006 - 2009, Liming Liang and Goncalo Abecasis,
+// 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.
+//
+///////////////////////////////////////////////////////////////////////////////
+
+
+
+#include "Random.h"
+#include <math.h>
+#include <iostream>
+
+using namespace std;
+
+#define Pi 3.1415926535897932
+
+double lngamma(double z)
+{
+	double result,sum;
+	static double c[7]={1.000000000190015, 76.18009172947146, -86.50532032941677, 24.01409824083091, -1.231739572450155, 0.1208650973866179e-2, -0.5395239384953e-5};
+	
+	sum=c[0];
+	for(int i=1;i<=6;i++)sum += c[i]/(z+i);
+	
+	result = (z+0.5)*log(z+5.5)-(z+5.5)+log(2.5066282746310005*sum/z);
+	
+	return result;
+}
+
+
+
+int poissonInt(double lambda) 
+{
+	static double mean=-1, waitTime,s,logmean,c;
+	double count, eventTime,ratio,y;
+	
+	static double maxInt=pow(2.0,(double)(8*sizeof(int)-1))-1;
+	
+	if(lambda <0){
+		cout<<"The mean of Poisson should be >=0, input="<<lambda<<endl;
+		exit(1);	
+	}
+	else if(lambda < 12){
+		if(lambda != mean){
+			mean = lambda;
+			waitTime = exp(-lambda);	
+		}
+		
+		count = 0;
+		eventTime = globalRandom.Next();
+		
+		while(eventTime > waitTime){
+			count++;
+			eventTime *= globalRandom.Next();	
+		}
+	}
+	else{
+		if(lambda != mean){
+			mean = lambda;
+			s=sqrt(2*lambda);
+			logmean=log(lambda);
+			c=lambda*log(lambda)-lngamma(lambda+1);
+		}
+		
+		do{
+			do{
+				y=tan(Pi*globalRandom.Next());
+				count = floor(s*y+lambda);
+			}while(count<0);
+			
+			ratio = 0.9*(1+y*y)*exp(count*logmean-lngamma(count+1)-c);
+			
+		}while(globalRandom.Next()>ratio);
+		
+	}
+	
+	if(count>maxInt)count=maxInt;
+	
+	return (int)count;		
+}
+
+
+
+
diff --git a/stochastic.h b/stochastic.h
new file mode 100644
index 0000000..63be72c
--- /dev/null
+++ b/stochastic.h
@@ -0,0 +1,42 @@
+////////////////////////////////////////////////////////////////////// 
+// stochastic.h 
+//////////////////////////////////////////////////////////////////////////////
+//              COPYRIGHT NOTICE FOR GENOME CODE
+//
+// Copyright (C) 2006 - 2009, Liming Liang and Goncalo Abecasis,
+// 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.
+//
+///////////////////////////////////////////////////////////////////////////////
+
+
+int poissonInt(double lambda);
+
+
+