From d49beaff6c97583411fd8c554217da48b86d0e43 Mon Sep 17 00:00:00 2001 From: Don Armstrong 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 +#include +#include +#include +#include "genome.h" +#include "Random.h" +#include +#include + +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 nSubSample(nSubPOP); + for(int i=0;i > 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 [ " <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;i0.0 && prop < 1.0){ // filter out the SNPs with MAF < maf + + cout<<"Rare SNPs filtered according to parameters."< keep(m,true); + + float p; + for(int i=0;iprop)keep[i]=0; + } + + + int index=0; + for(int k=0;k= 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 + +//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 +#include +#include +#include +#include +#include +#include +#include +#include + +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 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 > > chromosome; +static int * fragmentLength; +static int * branchFragment; +static int * geneIndex; +static int * branchMutation; + +static float * recombVector=NULL; + +struct popProfilesStruct{ + int generation; + vector popsize; + vector popStart; + vector< vector > selectProb; + vector< vector > popChange; +}; + +static vector popProfile; + +void print(){ + cout<<"\n*****************\nParents:\t"; + for(int i=0;i > &return_chromosome){ + + +// if(SNP<=0){ +// cout<<"Reserve size for return_chromosome="< > &return_chromosome){ + +// return_chromosome.clear(); +// return_chromosome.resize(n); +// for(int i=0;ipopProfile[i].popChange[j][k])min=popProfile[i].popChange[j][k]; + if(max recDistribution; + vector 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 tree; + ostringstream temp; + + + for(int i=0;i & 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 > &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;i0){ + cout<<"Populations change backward in time:"<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<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 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<= randomLong && begin != 0)){ + cout<<"searching error: begin="<=randomLong)sampleIndex = begin; + else sampleIndex = end; + + if(genepool[sampleIndex]==0){ + cout<<"genepool error: begin="<0){ + geneIndex[i*L+j] = fragmentLength[j]; + fragmentLength[j] += branchMutation[i*L+j]; + } + branchFragment[i*L+j] = j; + } + } + + + for(int i=0;i0 && branchFragment[genepool[i]]== -1 ){ + geneIndex[genepool[i]] = fragmentLength[branchFragment[i]]; + fragmentLength[branchFragment[i]] += branchMutation[genepool[i]]; + } + branchFragment[genepool[i]]=branchFragment[i]; + } + } + + //cout<<"branchfragment done"<1){ + for(int i=0;i +#include + +using namespace std; + +void genome(string popSize, // effective population size of each population + int nSubPOP, // number of populations + vector & 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 > &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 +#include + +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="< 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.2