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); + + +