From d49beaff6c97583411fd8c554217da48b86d0e43 Mon Sep 17 00:00:00 2001
From: Don Armstrong <don@donarmstrong.com>
Date: Fri, 7 Feb 2014 15:18:51 -0800
Subject: [PATCH 1/1] import initial version of genome
---
Error.cpp | 88 +++
Error.h | 55 ++
LICENSE.GENOME | 36 ++
LICENSE.twister | 37 ++
Main.cpp | 211 +++++++
MathConstant.h | 73 +++
Random.cpp | 387 ++++++++++++
Random.h | 115 ++++
genome.cpp | 1396 +++++++++++++++++++++++++++++++++++++++++
genome.h | 62 ++
makefile | 27 +
makefile-32bitsys | 27 +
population.txt | 7 +
readme.txt | 84 +++
recombination-pos.txt | 6 +
recombination.txt | 3 +
stochastic.cpp | 115 ++++
stochastic.h | 42 ++
18 files changed, 2771 insertions(+)
create mode 100644 Error.cpp
create mode 100644 Error.h
create mode 100644 LICENSE.GENOME
create mode 100644 LICENSE.twister
create mode 100644 Main.cpp
create mode 100644 MathConstant.h
create mode 100644 Random.cpp
create mode 100644 Random.h
create mode 100644 genome.cpp
create mode 100644 genome.h
create mode 100644 makefile
create mode 100644 makefile-32bitsys
create mode 100644 population.txt
create mode 100644 readme.txt
create mode 100644 recombination-pos.txt
create mode 100644 recombination.txt
create mode 100644 stochastic.cpp
create mode 100644 stochastic.h
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);
+
+
+
--
2.39.5