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