[med-svn] [gwama] 02/04: Imported Upstream version 2.1
Dylan Aïssi
bob.dybian-guest at moszumanska.debian.org
Wed Nov 19 08:15:39 UTC 2014
This is an automated email from the git hooks/post-receive script.
bob.dybian-guest pushed a commit to branch master
in repository gwama.
commit 909413013fdf7b6a843bf26f295f02dd4f29ac22
Author: Dylan Aïssi <bob.dybian at gmail.com>
Date: Wed Nov 19 08:54:34 2014 +0100
Imported Upstream version 2.1
---
Makefile | 12 +
__MACOSX/._Makefile | Bin 0 -> 229 bytes
__MACOSX/._ap.cpp | Bin 0 -> 229 bytes
__MACOSX/._ap.h | Bin 0 -> 229 bytes
__MACOSX/._apvt.h | Bin 0 -> 229 bytes
__MACOSX/._chisquaredistr.cpp | Bin 0 -> 229 bytes
__MACOSX/._chisquaredistr.h | Bin 0 -> 229 bytes
__MACOSX/._cohort.cpp | Bin 0 -> 229 bytes
__MACOSX/._cohort.h | Bin 0 -> 229 bytes
__MACOSX/._commandLine.cpp | Bin 0 -> 171 bytes
__MACOSX/._commandLine.h | Bin 0 -> 229 bytes
__MACOSX/._gammaf.cpp | Bin 0 -> 229 bytes
__MACOSX/._gammaf.h | Bin 0 -> 229 bytes
__MACOSX/._global.cpp | Bin 0 -> 171 bytes
__MACOSX/._global.h | Bin 0 -> 171 bytes
__MACOSX/._igammaf.cpp | Bin 0 -> 229 bytes
__MACOSX/._igammaf.h | Bin 0 -> 229 bytes
__MACOSX/._main.cpp | Bin 0 -> 171 bytes
__MACOSX/._marker.cpp | Bin 0 -> 171 bytes
__MACOSX/._marker.h | Bin 0 -> 229 bytes
__MACOSX/._normaldistr.cpp | Bin 0 -> 229 bytes
__MACOSX/._normaldistr.h | Bin 0 -> 229 bytes
__MACOSX/._problem.cpp | Bin 0 -> 229 bytes
__MACOSX/._problem.h | Bin 0 -> 229 bytes
__MACOSX/._readFile.cpp | Bin 0 -> 171 bytes
__MACOSX/._readFile.h | Bin 0 -> 229 bytes
__MACOSX/._resource.h | Bin 0 -> 229 bytes
__MACOSX/._statistics.cpp | Bin 0 -> 229 bytes
__MACOSX/._statistics.h | Bin 0 -> 229 bytes
__MACOSX/._stdafx.h | Bin 0 -> 229 bytes
__MACOSX/._study.cpp | Bin 0 -> 229 bytes
__MACOSX/._study.h | Bin 0 -> 229 bytes
ap.cpp | 459 ++++++++++++++++++
ap.h | 582 ++++++++++++++++++++++
apvt.h | 754 ++++++++++++++++++++++++++++
chisquaredistr.cpp | 157 ++++++
chisquaredistr.h | 143 ++++++
cohort.cpp | 47 ++
cohort.h | 52 ++
commandLine.cpp | 378 +++++++++++++++
commandLine.h | 40 ++
gammaf.cpp | 338 +++++++++++++
gammaf.h | 103 ++++
global.cpp | 75 +++
global.h | 75 +++
igammaf.cpp | 430 ++++++++++++++++
igammaf.h | 156 ++++++
main.cpp | 181 +++++++
marker.cpp | 1078 +++++++++++++++++++++++++++++++++++++++++
marker.h | 151 ++++++
normaldistr.cpp | 377 ++++++++++++++
normaldistr.h | 174 +++++++
problem.cpp | 15 +
problem.h | 15 +
readFile.cpp | 585 ++++++++++++++++++++++
readFile.h | 50 ++
resource.h | 14 +
statistics.cpp | 559 +++++++++++++++++++++
statistics.h | 52 ++
stdafx.h | 0
study.cpp | 124 +++++
study.h | 56 +++
tools.cpp | 110 +++++
tools.h | 28 ++
64 files changed, 7370 insertions(+)
diff --git a/Makefile b/Makefile
new file mode 100755
index 0000000..9be1847
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,12 @@
+#GWAMA program
+
+VERSION = 2.7
+
+CC = g++
+
+DEBUGFLAGS = -Wno-deprecated -O3
+
+GWAMA: main.cpp
+
+ g++ main.cpp marker.cpp statistics.cpp study.cpp chisquaredistr.cpp normaldistr.cpp gammaf.cpp igammaf.cpp ap.cpp global.cpp problem.cpp tools.cpp cohort.cpp commandLine.cpp readFile.cpp $(DEBUGFLAGS) -o GWAMA
+
diff --git a/__MACOSX/._Makefile b/__MACOSX/._Makefile
new file mode 100644
index 0000000..f55cdb1
Binary files /dev/null and b/__MACOSX/._Makefile differ
diff --git a/__MACOSX/._ap.cpp b/__MACOSX/._ap.cpp
new file mode 100644
index 0000000..f55cdb1
Binary files /dev/null and b/__MACOSX/._ap.cpp differ
diff --git a/__MACOSX/._ap.h b/__MACOSX/._ap.h
new file mode 100644
index 0000000..f55cdb1
Binary files /dev/null and b/__MACOSX/._ap.h differ
diff --git a/__MACOSX/._apvt.h b/__MACOSX/._apvt.h
new file mode 100644
index 0000000..f55cdb1
Binary files /dev/null and b/__MACOSX/._apvt.h differ
diff --git a/__MACOSX/._chisquaredistr.cpp b/__MACOSX/._chisquaredistr.cpp
new file mode 100644
index 0000000..f55cdb1
Binary files /dev/null and b/__MACOSX/._chisquaredistr.cpp differ
diff --git a/__MACOSX/._chisquaredistr.h b/__MACOSX/._chisquaredistr.h
new file mode 100644
index 0000000..f55cdb1
Binary files /dev/null and b/__MACOSX/._chisquaredistr.h differ
diff --git a/__MACOSX/._cohort.cpp b/__MACOSX/._cohort.cpp
new file mode 100644
index 0000000..f55cdb1
Binary files /dev/null and b/__MACOSX/._cohort.cpp differ
diff --git a/__MACOSX/._cohort.h b/__MACOSX/._cohort.h
new file mode 100644
index 0000000..f55cdb1
Binary files /dev/null and b/__MACOSX/._cohort.h differ
diff --git a/__MACOSX/._commandLine.cpp b/__MACOSX/._commandLine.cpp
new file mode 100644
index 0000000..14f037e
Binary files /dev/null and b/__MACOSX/._commandLine.cpp differ
diff --git a/__MACOSX/._commandLine.h b/__MACOSX/._commandLine.h
new file mode 100644
index 0000000..f55cdb1
Binary files /dev/null and b/__MACOSX/._commandLine.h differ
diff --git a/__MACOSX/._gammaf.cpp b/__MACOSX/._gammaf.cpp
new file mode 100644
index 0000000..f55cdb1
Binary files /dev/null and b/__MACOSX/._gammaf.cpp differ
diff --git a/__MACOSX/._gammaf.h b/__MACOSX/._gammaf.h
new file mode 100644
index 0000000..f55cdb1
Binary files /dev/null and b/__MACOSX/._gammaf.h differ
diff --git a/__MACOSX/._global.cpp b/__MACOSX/._global.cpp
new file mode 100644
index 0000000..14f037e
Binary files /dev/null and b/__MACOSX/._global.cpp differ
diff --git a/__MACOSX/._global.h b/__MACOSX/._global.h
new file mode 100644
index 0000000..14f037e
Binary files /dev/null and b/__MACOSX/._global.h differ
diff --git a/__MACOSX/._igammaf.cpp b/__MACOSX/._igammaf.cpp
new file mode 100644
index 0000000..f55cdb1
Binary files /dev/null and b/__MACOSX/._igammaf.cpp differ
diff --git a/__MACOSX/._igammaf.h b/__MACOSX/._igammaf.h
new file mode 100644
index 0000000..f55cdb1
Binary files /dev/null and b/__MACOSX/._igammaf.h differ
diff --git a/__MACOSX/._main.cpp b/__MACOSX/._main.cpp
new file mode 100644
index 0000000..14f037e
Binary files /dev/null and b/__MACOSX/._main.cpp differ
diff --git a/__MACOSX/._marker.cpp b/__MACOSX/._marker.cpp
new file mode 100644
index 0000000..14f037e
Binary files /dev/null and b/__MACOSX/._marker.cpp differ
diff --git a/__MACOSX/._marker.h b/__MACOSX/._marker.h
new file mode 100644
index 0000000..f55cdb1
Binary files /dev/null and b/__MACOSX/._marker.h differ
diff --git a/__MACOSX/._normaldistr.cpp b/__MACOSX/._normaldistr.cpp
new file mode 100644
index 0000000..f55cdb1
Binary files /dev/null and b/__MACOSX/._normaldistr.cpp differ
diff --git a/__MACOSX/._normaldistr.h b/__MACOSX/._normaldistr.h
new file mode 100644
index 0000000..f55cdb1
Binary files /dev/null and b/__MACOSX/._normaldistr.h differ
diff --git a/__MACOSX/._problem.cpp b/__MACOSX/._problem.cpp
new file mode 100644
index 0000000..f55cdb1
Binary files /dev/null and b/__MACOSX/._problem.cpp differ
diff --git a/__MACOSX/._problem.h b/__MACOSX/._problem.h
new file mode 100644
index 0000000..f55cdb1
Binary files /dev/null and b/__MACOSX/._problem.h differ
diff --git a/__MACOSX/._readFile.cpp b/__MACOSX/._readFile.cpp
new file mode 100644
index 0000000..14f037e
Binary files /dev/null and b/__MACOSX/._readFile.cpp differ
diff --git a/__MACOSX/._readFile.h b/__MACOSX/._readFile.h
new file mode 100644
index 0000000..f55cdb1
Binary files /dev/null and b/__MACOSX/._readFile.h differ
diff --git a/__MACOSX/._resource.h b/__MACOSX/._resource.h
new file mode 100644
index 0000000..f55cdb1
Binary files /dev/null and b/__MACOSX/._resource.h differ
diff --git a/__MACOSX/._statistics.cpp b/__MACOSX/._statistics.cpp
new file mode 100644
index 0000000..f55cdb1
Binary files /dev/null and b/__MACOSX/._statistics.cpp differ
diff --git a/__MACOSX/._statistics.h b/__MACOSX/._statistics.h
new file mode 100644
index 0000000..f55cdb1
Binary files /dev/null and b/__MACOSX/._statistics.h differ
diff --git a/__MACOSX/._stdafx.h b/__MACOSX/._stdafx.h
new file mode 100644
index 0000000..f55cdb1
Binary files /dev/null and b/__MACOSX/._stdafx.h differ
diff --git a/__MACOSX/._study.cpp b/__MACOSX/._study.cpp
new file mode 100644
index 0000000..f55cdb1
Binary files /dev/null and b/__MACOSX/._study.cpp differ
diff --git a/__MACOSX/._study.h b/__MACOSX/._study.h
new file mode 100644
index 0000000..f55cdb1
Binary files /dev/null and b/__MACOSX/._study.h differ
diff --git a/ap.cpp b/ap.cpp
new file mode 100755
index 0000000..9154b17
--- /dev/null
+++ b/ap.cpp
@@ -0,0 +1,459 @@
+/********************************************************************
+AP Library version 1.2
+Copyright (c) 2003-2007, Sergey Bochkanov (ALGLIB project).
+See www.alglib.net or alglib.sources.ru for details.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+- Redistributions of source code must retain the above copyright
+ notice, this list of conditions and the following disclaimer.
+
+- Redistributions in binary form must reproduce the above copyright
+ notice, this list of conditions and the following disclaimer listed
+ in this license in the documentation and/or other materials
+ provided with the distribution.
+
+- Neither the name of the copyright holders nor the names of its
+ contributors may 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 "stdafx.h"
+#include "ap.h"
+
+
+/********************************************************************
+Optimized ABLAS interface
+********************************************************************/
+#ifdef AP_WIN32
+#include <windows.h>
+extern "C"
+{
+typedef double (*_ddot1)(const double*, const double*, long);
+typedef void (*_dmove1)(const double*, const double*, long);
+typedef void (*_dmoves1)(const double*, const double*, long, double);
+typedef void (*_dmoveneg1)(const double*, const double*, long);
+typedef void (*_dadd1)(const double*, const double*, long);
+typedef void (*_dadds1)(const double*, const double*, long, double);
+typedef void (*_dsub1)(const double*, const double*, long);
+typedef void (*_dmuls1)(const double*, long, double);
+}
+HINSTANCE ABLAS = LoadLibrary("ablas.dll");
+
+static _ddot1 ddot1 = ABLAS==NULL ? NULL : (_ddot1) GetProcAddress(ABLAS, "ASMDotProduct1");
+static _dmove1 dmove1 = ABLAS==NULL ? NULL : (_dmove1) GetProcAddress(ABLAS, "ASMMove1");
+static _dmoves1 dmoves1 = ABLAS==NULL ? NULL : (_dmoves1) GetProcAddress(ABLAS, "ASMMoveS1");
+static _dmoveneg1 dmoveneg1 = ABLAS==NULL ? NULL : (_dmoveneg1) GetProcAddress(ABLAS, "ASMMoveNeg1");
+static _dadd1 dadd1 = ABLAS==NULL ? NULL : (_dadd1) GetProcAddress(ABLAS, "ASMAdd1");
+static _dadds1 dadds1 = ABLAS==NULL ? NULL : (_dadds1) GetProcAddress(ABLAS, "ASMAddS1");
+static _dsub1 dsub1 = ABLAS==NULL ? NULL : (_dsub1) GetProcAddress(ABLAS, "ASMSub1");
+static _dmuls1 dmuls1 = ABLAS==NULL ? NULL : (_dmuls1) GetProcAddress(ABLAS, "ASMMulS1");
+#endif
+
+const double ap::machineepsilon = 5E-16;
+const double ap::maxrealnumber = 1E300;
+const double ap::minrealnumber = 1E-300;
+
+/********************************************************************
+ap::complex operations
+********************************************************************/
+const bool ap::operator==(const ap::complex& lhs, const ap::complex& rhs)
+{ return lhs.x==rhs.x && lhs.y==rhs.y; }
+
+const bool ap::operator!=(const ap::complex& lhs, const ap::complex& rhs)
+{ return lhs.x!=rhs.x || lhs.y!=rhs.y; }
+
+const ap::complex ap::operator+(const ap::complex& lhs)
+{ return lhs; }
+
+const ap::complex ap::operator-(const ap::complex& lhs)
+{ return ap::complex(-lhs.x, -lhs.y); }
+
+const ap::complex ap::operator+(const ap::complex& lhs, const ap::complex& rhs)
+{ ap::complex r = lhs; r += rhs; return r; }
+
+const ap::complex ap::operator+(const ap::complex& lhs, const double& rhs)
+{ ap::complex r = lhs; r += rhs; return r; }
+
+const ap::complex ap::operator+(const double& lhs, const ap::complex& rhs)
+{ ap::complex r = rhs; r += lhs; return r; }
+
+const ap::complex ap::operator-(const ap::complex& lhs, const ap::complex& rhs)
+{ ap::complex r = lhs; r -= rhs; return r; }
+
+const ap::complex ap::operator-(const ap::complex& lhs, const double& rhs)
+{ ap::complex r = lhs; r -= rhs; return r; }
+
+const ap::complex ap::operator-(const double& lhs, const ap::complex& rhs)
+{ ap::complex r = lhs; r -= rhs; return r; }
+
+const ap::complex ap::operator*(const ap::complex& lhs, const ap::complex& rhs)
+{ return ap::complex(lhs.x*rhs.x - lhs.y*rhs.y, lhs.x*rhs.y + lhs.y*rhs.x); }
+
+const ap::complex ap::operator*(const ap::complex& lhs, const double& rhs)
+{ return ap::complex(lhs.x*rhs, lhs.y*rhs); }
+
+const ap::complex ap::operator*(const double& lhs, const ap::complex& rhs)
+{ return ap::complex(lhs*rhs.x, lhs*rhs.y); }
+
+const ap::complex ap::operator/(const ap::complex& lhs, const ap::complex& rhs)
+{
+ ap::complex result;
+ double e;
+ double f;
+ if( fabs(rhs.y)<fabs(rhs.x) )
+ {
+ e = rhs.y/rhs.x;
+ f = rhs.x+rhs.y*e;
+ result.x = (lhs.x+lhs.y*e)/f;
+ result.y = (lhs.y-lhs.x*e)/f;
+ }
+ else
+ {
+ e = rhs.x/rhs.y;
+ f = rhs.y+rhs.x*e;
+ result.x = (lhs.y+lhs.x*e)/f;
+ result.y = (-lhs.x+lhs.y*e)/f;
+ }
+ return result;
+}
+
+const ap::complex ap::operator/(const double& lhs, const ap::complex& rhs)
+{
+ ap::complex result;
+ double e;
+ double f;
+ if( fabs(rhs.y)<fabs(rhs.x) )
+ {
+ e = rhs.y/rhs.x;
+ f = rhs.x+rhs.y*e;
+ result.x = lhs/f;
+ result.y = -lhs*e/f;
+ }
+ else
+ {
+ e = rhs.x/rhs.y;
+ f = rhs.y+rhs.x*e;
+ result.x = lhs*e/f;
+ result.y = -lhs/f;
+ }
+ return result;
+}
+
+const ap::complex ap::operator/(const ap::complex& lhs, const double& rhs)
+{ return ap::complex(lhs.x/rhs, lhs.y/rhs); }
+
+const double ap::abscomplex(const ap::complex &z)
+{
+ double w;
+ double xabs;
+ double yabs;
+ double v;
+
+ xabs = fabs(z.x);
+ yabs = fabs(z.y);
+ w = xabs>yabs ? xabs : yabs;
+ v = xabs<yabs ? xabs : yabs;
+ if( v==0 )
+ return w;
+ else
+ {
+ double t = v/w;
+ return w*sqrt(1+t*t);
+ }
+}
+
+const ap::complex ap::conj(const ap::complex &z)
+{ return ap::complex(z.x, -z.y); }
+
+const ap::complex ap::csqr(const ap::complex &z)
+{ return ap::complex(z.x*z.x-z.y*z.y, 2*z.x*z.y); }
+
+/********************************************************************
+BLAS functions
+********************************************************************/
+double ap::vdotproduct(const double *v1, const double *v2, int N)
+{
+#ifdef AP_WIN32
+ if( ddot1!=NULL )
+ return ddot1(v1, v2, N);
+#endif
+ return ap::_vdotproduct<double>(v1, v2, N);
+}
+
+ap::complex ap::vdotproduct(const ap::complex *v1, const ap::complex *v2, int N)
+{
+ return ap::_vdotproduct<ap::complex>(v1, v2, N);
+}
+
+void ap::vmove(double *vdst, const double* vsrc, int N)
+{
+#ifdef AP_WIN32
+ if( dmove1!=NULL )
+ {
+ dmove1(vdst, vsrc, N);
+ return;
+ }
+#endif
+ ap::_vmove<double>(vdst, vsrc, N);
+}
+
+void ap::vmove(ap::complex *vdst, const ap::complex* vsrc, int N)
+{
+ ap::_vmove<ap::complex>(vdst, vsrc, N);
+}
+
+void ap::vmoveneg(double *vdst, const double *vsrc, int N)
+{
+#ifdef AP_WIN32
+ if( dmoveneg1!=NULL )
+ {
+ dmoveneg1(vdst, vsrc, N);
+ return;
+ }
+#endif
+ ap::_vmoveneg<double>(vdst, vsrc, N);
+}
+
+void ap::vmoveneg(ap::complex *vdst, const ap::complex *vsrc, int N)
+{
+ ap::_vmoveneg<ap::complex>(vdst, vsrc, N);
+}
+
+void ap::vmove(double *vdst, const double *vsrc, int N, double alpha)
+{
+#ifdef AP_WIN32
+ if( dmoves1!=NULL )
+ {
+ dmoves1(vdst, vsrc, N, alpha);
+ return;
+ }
+#endif
+ ap::_vmove<double,double>(vdst, vsrc, N, alpha);
+}
+
+void ap::vmove(ap::complex *vdst, const ap::complex *vsrc, int N, double alpha)
+{
+ ap::_vmove<ap::complex,double>(vdst, vsrc, N, alpha);
+}
+
+void ap::vmove(ap::complex *vdst, const ap::complex *vsrc, int N, ap::complex alpha)
+{
+ ap::_vmove<ap::complex,ap::complex>(vdst, vsrc, N, alpha);
+}
+
+void ap::vadd(double *vdst, const double *vsrc, int N)
+{
+#ifdef AP_WIN32
+ if( dadd1!=NULL )
+ {
+ dadd1(vdst, vsrc, N);
+ return;
+ }
+#endif
+ ap::_vadd<double>(vdst, vsrc, N);
+}
+
+void ap::vadd(ap::complex *vdst, const ap::complex *vsrc, int N)
+{
+ ap::_vadd<ap::complex>(vdst, vsrc, N);
+}
+
+void ap::vadd(double *vdst, const double *vsrc, int N, double alpha)
+{
+#ifdef AP_WIN32
+ if( dadds1!=NULL )
+ {
+ dadds1(vdst, vsrc, N, alpha);
+ return;
+ }
+#endif
+ ap::_vadd<double,double>(vdst, vsrc, N, alpha);
+}
+
+void ap::vadd(ap::complex *vdst, const ap::complex *vsrc, int N, double alpha)
+{
+ ap::_vadd<ap::complex,double>(vdst, vsrc, N, alpha);
+}
+
+void ap::vadd(ap::complex *vdst, const ap::complex *vsrc, int N, ap::complex alpha)
+{
+ ap::_vadd<ap::complex,ap::complex>(vdst, vsrc, N, alpha);
+}
+
+void ap::vsub(double *vdst, const double *vsrc, int N)
+{
+#ifdef AP_WIN32
+ if( dsub1!=NULL )
+ {
+ dsub1(vdst, vsrc, N);
+ return;
+ }
+#endif
+ ap::_vsub<double>(vdst, vsrc, N);
+}
+
+void ap::vsub(ap::complex *vdst, const ap::complex *vsrc, int N)
+{
+ ap::_vsub<ap::complex>(vdst, vsrc, N);
+}
+
+void ap::vsub(double *vdst, const double *vsrc, int N, double alpha)
+{
+#ifdef AP_WIN32
+ if( dadds1!=NULL )
+ {
+ dadds1(vdst, vsrc, N, -alpha);
+ return;
+ }
+#endif
+ ap::_vsub<double,double>(vdst, vsrc, N, alpha);
+}
+
+void ap::vsub(ap::complex *vdst, const ap::complex *vsrc, int N, double alpha)
+{
+ ap::_vsub<ap::complex,double>(vdst, vsrc, N, alpha);
+}
+
+void ap::vsub(ap::complex *vdst, const ap::complex *vsrc, int N, ap::complex alpha)
+{
+ ap::_vsub<ap::complex,ap::complex>(vdst, vsrc, N, alpha);
+}
+
+void ap::vmul(double *vdst, int N, double alpha)
+{
+#ifdef AP_WIN32
+ if( dmuls1!=NULL )
+ {
+ dmuls1(vdst, N, alpha);
+ return;
+ }
+#endif
+ ap::_vmul<double,double>(vdst, N, alpha);
+}
+
+void ap::vmul(ap::complex *vdst, int N, double alpha)
+{
+ ap::_vmul<ap::complex,double>(vdst, N, alpha);
+}
+
+void ap::vmul(ap::complex *vdst, int N, ap::complex alpha)
+{
+ ap::_vmul<ap::complex,ap::complex>(vdst, N, alpha);
+}
+
+/********************************************************************
+standard functions
+********************************************************************/
+int ap::sign(double x)
+{
+ if( x>0 ) return 1;
+ if( x<0 ) return -1;
+ return 0;
+}
+
+double ap::randomreal()
+{
+ int i = rand();
+ while(i==RAND_MAX)
+ i =rand();
+ return double(i)/double(RAND_MAX);
+}
+
+int ap::randominteger(int maxv)
+{ return rand()%maxv; }
+
+int ap::round(double x)
+{ return int(floor(x+0.5)); }
+
+int ap::trunc(double x)
+{ return int(x>0 ? floor(x) : ceil(x)); }
+
+int ap::ifloor(double x)
+{ return int(floor(x)); }
+
+int ap::iceil(double x)
+{ return int(ceil(x)); }
+
+double ap::pi()
+{ return 3.14159265358979323846; }
+
+double ap::sqr(double x)
+{ return x*x; }
+
+int ap::maxint(int m1, int m2)
+{
+ return m1>m2 ? m1 : m2;
+}
+
+int ap::minint(int m1, int m2)
+{
+ return m1>m2 ? m2 : m1;
+}
+
+double ap::maxreal(double m1, double m2)
+{
+ return m1>m2 ? m1 : m2;
+}
+
+double ap::minreal(double m1, double m2)
+{
+ return m1>m2 ? m2 : m1;
+}
+
+/********************************************************************
+Service routines:
+********************************************************************/
+void* ap::amalloc(size_t size, size_t alignment)
+{
+ if( alignment<=1 )
+ {
+ //
+ // no alignment, just call malloc
+ //
+ void *block = malloc(sizeof(void*)+size);
+ void **p = (void**)block;
+ *p = block;
+ return (void*)((char*)block+sizeof(void*));
+ }
+ else
+ {
+ //
+ // align.
+ //
+ void *block = malloc(alignment-1+sizeof(void*)+size);
+ char *result = (char*)block+sizeof(void*);
+ //if( ((unsigned int)(result))%alignment!=0 )
+ // result += alignment - ((unsigned int)(result))%alignment;
+ if( (result-(char*)0)%alignment!=0 )
+ result += alignment - (result-(char*)0)%alignment;
+ *((void**)(result-sizeof(void*))) = block;
+ return result;
+ }
+}
+
+void ap::afree(void *block)
+{
+ void *p = *((void**)((char*)block-sizeof(void*)));
+ free(p);
+}
+
+int ap::vlen(int n1, int n2)
+{
+ return n2-n1+1;
+}
+
diff --git a/ap.h b/ap.h
new file mode 100755
index 0000000..98d5893
--- /dev/null
+++ b/ap.h
@@ -0,0 +1,582 @@
+/********************************************************************
+AP Library version 1.2.1
+
+Copyright (c) 2003-2007, Sergey Bochkanov (ALGLIB project).
+See www.alglib.net or alglib.sources.ru for details.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+- Redistributions of source code must retain the above copyright
+ notice, this list of conditions and the following disclaimer.
+
+- Redistributions in binary form must reproduce the above copyright
+ notice, this list of conditions and the following disclaimer listed
+ in this license in the documentation and/or other materials
+ provided with the distribution.
+
+- Neither the name of the copyright holders nor the names of its
+ contributors may 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 AP_H
+#define AP_H
+
+#include <stdlib.h>
+#include <string>
+#include <math.h>
+
+/********************************************************************
+Array bounds check
+********************************************************************/
+#define AP_ASSERT
+
+#ifndef AP_ASSERT //
+#define NO_AP_ASSERT // This code avoids definition of the
+#endif // both AP_ASSERT and NO_AP_ASSERT symbols
+#ifdef NO_AP_ASSERT //
+#ifdef AP_ASSERT //
+#undef NO_AP_ASSERT //
+#endif //
+#endif //
+
+
+/********************************************************************
+Current environment.
+********************************************************************/
+#ifndef AP_WIN32
+#ifndef AP_UNKNOWN
+#define AP_UNKNOWN
+#endif
+#endif
+#ifdef AP_WIN32
+#ifdef AP_UNKNOWN
+#error Multiple environments are declared!
+#endif
+#endif
+
+/********************************************************************
+This symbol is used for debugging. Do not define it and do not remove
+comments.
+********************************************************************/
+//#define UNSAFE_MEM_COPY
+
+
+/********************************************************************
+Namespace of a standard library AlgoPascal.
+********************************************************************/
+namespace ap
+{
+
+/********************************************************************
+Service routines:
+ amalloc - allocates an aligned block of size bytes
+ afree - frees block allocated by amalloc
+ vlen - just alias for n2-n1+1
+********************************************************************/
+void* amalloc(size_t size, size_t alignment);
+void afree(void *block);
+int vlen(int n1, int n2);
+
+/********************************************************************
+Exception class.
+********************************************************************/
+class ap_error
+{
+public:
+ ap_error(){};
+ ap_error(const char *s){ msg = s; };
+
+ std::string msg;
+
+ static void make_assertion(bool bClause)
+ { if(!bClause) throw ap_error(); };
+ static void make_assertion(bool bClause, const char *msg)
+ { if(!bClause) throw ap_error(msg); };
+private:
+};
+
+/********************************************************************
+Class defining a complex number with double precision.
+********************************************************************/
+class complex;
+
+class complex
+{
+public:
+ complex():x(0.0),y(0.0){};
+ complex(const double &_x):x(_x),y(0.0){};
+ complex(const double &_x, const double &_y):x(_x),y(_y){};
+ complex(const complex &z):x(z.x),y(z.y){};
+
+ complex& operator= (const double& v){ x = v; y = 0.0; return *this; };
+ complex& operator+=(const double& v){ x += v; return *this; };
+ complex& operator-=(const double& v){ x -= v; return *this; };
+ complex& operator*=(const double& v){ x *= v; y *= v; return *this; };
+ complex& operator/=(const double& v){ x /= v; y /= v; return *this; };
+
+ complex& operator= (const complex& z){ x = z.x; y = z.y; return *this; };
+ complex& operator+=(const complex& z){ x += z.x; y += z.y; return *this; };
+ complex& operator-=(const complex& z){ x -= z.x; y -= z.y; return *this; };
+ complex& operator*=(const complex& z){ double t = x*z.x-y*z.y; y = x*z.y+y*z.x; x = t; return *this; };
+ complex& operator/=(const complex& z)
+ {
+ ap::complex result;
+ double e;
+ double f;
+ if( fabs(z.y)<fabs(z.x) )
+ {
+ e = z.y/z.x;
+ f = z.x+z.y*e;
+ result.x = (z.x+z.y*e)/f;
+ result.y = (z.y-z.x*e)/f;
+ }
+ else
+ {
+ e = z.x/z.y;
+ f = z.y+z.x*e;
+ result.x = (z.y+z.x*e)/f;
+ result.y = (-z.x+z.y*e)/f;
+ }
+ *this = result;
+ return *this;
+ };
+
+ double x, y;
+};
+
+const complex operator/(const complex& lhs, const complex& rhs);
+const bool operator==(const complex& lhs, const complex& rhs);
+const bool operator!=(const complex& lhs, const complex& rhs);
+const complex operator+(const complex& lhs);
+const complex operator-(const complex& lhs);
+const complex operator+(const complex& lhs, const complex& rhs);
+const complex operator+(const complex& lhs, const double& rhs);
+const complex operator+(const double& lhs, const complex& rhs);
+const complex operator-(const complex& lhs, const complex& rhs);
+const complex operator-(const complex& lhs, const double& rhs);
+const complex operator-(const double& lhs, const complex& rhs);
+const complex operator*(const complex& lhs, const complex& rhs);
+const complex operator*(const complex& lhs, const double& rhs);
+const complex operator*(const double& lhs, const complex& rhs);
+const complex operator/(const complex& lhs, const complex& rhs);
+const complex operator/(const double& lhs, const complex& rhs);
+const complex operator/(const complex& lhs, const double& rhs);
+const double abscomplex(const complex &z);
+const complex conj(const complex &z);
+const complex csqr(const complex &z);
+
+
+/********************************************************************
+Templates for vector operations
+********************************************************************/
+#include "apvt.h"
+
+/********************************************************************
+BLAS functions
+********************************************************************/
+double vdotproduct(const double *v1, const double *v2, int N);
+complex vdotproduct(const complex *v1, const complex *v2, int N);
+
+void vmove(double *vdst, const double* vsrc, int N);
+void vmove(complex *vdst, const complex* vsrc, int N);
+
+void vmoveneg(double *vdst, const double *vsrc, int N);
+void vmoveneg(complex *vdst, const complex *vsrc, int N);
+
+void vmove(double *vdst, const double *vsrc, int N, double alpha);
+void vmove(complex *vdst, const complex *vsrc, int N, double alpha);
+void vmove(complex *vdst, const complex *vsrc, int N, complex alpha);
+
+void vadd(double *vdst, const double *vsrc, int N);
+void vadd(complex *vdst, const complex *vsrc, int N);
+
+void vadd(double *vdst, const double *vsrc, int N, double alpha);
+void vadd(complex *vdst, const complex *vsrc, int N, double alpha);
+void vadd(complex *vdst, const complex *vsrc, int N, complex alpha);
+
+void vsub(double *vdst, const double *vsrc, int N);
+void vsub(complex *vdst, const complex *vsrc, int N);
+
+void vsub(double *vdst, const double *vsrc, int N, double alpha);
+void vsub(complex *vdst, const complex *vsrc, int N, double alpha);
+void vsub(complex *vdst, const complex *vsrc, int N, complex alpha);
+
+void vmul(double *vdst, int N, double alpha);
+void vmul(complex *vdst, int N, double alpha);
+void vmul(complex *vdst, int N, complex alpha);
+
+
+/********************************************************************
+Template of a dynamical one-dimensional array
+********************************************************************/
+template<class T, bool Aligned = false>
+class template_1d_array
+{
+public:
+ template_1d_array()
+ {
+ m_Vec=0;
+ m_iVecSize = 0;
+ m_iLow = 0;
+ m_iHigh = -1;
+ };
+
+ ~template_1d_array()
+ {
+ if(m_Vec)
+ {
+ if( Aligned )
+ ap::afree(m_Vec);
+ else
+ delete[] m_Vec;
+ }
+ };
+
+ template_1d_array(const template_1d_array &rhs)
+ {
+ m_Vec=0;
+ m_iVecSize = 0;
+ m_iLow = 0;
+ m_iHigh = -1;
+ if( rhs.m_iVecSize!=0 )
+ setcontent(rhs.m_iLow, rhs.m_iHigh, rhs.getcontent());
+ };
+
+
+ const template_1d_array& operator=(const template_1d_array &rhs)
+ {
+ if( this==&rhs )
+ return *this;
+
+ if( rhs.m_iVecSize!=0 )
+ setcontent(rhs.m_iLow, rhs.m_iHigh, rhs.getcontent());
+ else
+ {
+ m_Vec=0;
+ m_iVecSize = 0;
+ m_iLow = 0;
+ m_iHigh = -1;
+ }
+ return *this;
+ };
+
+
+ const T& operator()(int i) const
+ {
+ #ifndef NO_AP_ASSERT
+ ap_error::make_assertion(i>=m_iLow && i<=m_iHigh);
+ #endif
+ return m_Vec[ i-m_iLow ];
+ };
+
+
+ T& operator()(int i)
+ {
+ #ifndef NO_AP_ASSERT
+ ap_error::make_assertion(i>=m_iLow && i<=m_iHigh);
+ #endif
+ return m_Vec[ i-m_iLow ];
+ };
+
+
+ void setbounds( int iLow, int iHigh )
+ {
+ if(m_Vec)
+ {
+ if( Aligned )
+ ap::afree(m_Vec);
+ else
+ delete[] m_Vec;
+ }
+ m_iLow = iLow;
+ m_iHigh = iHigh;
+ m_iVecSize = iHigh-iLow+1;
+ if( Aligned )
+ m_Vec = (T*)ap::amalloc(m_iVecSize*sizeof(T), 16);
+ else
+ m_Vec = new T[m_iVecSize];
+ };
+
+
+ void setcontent( int iLow, int iHigh, const T *pContent )
+ {
+ setbounds(iLow, iHigh);
+ for(int i=0; i<m_iVecSize; i++)
+ m_Vec[i] = pContent[i];
+ };
+
+
+ T* getcontent()
+ {
+ return m_Vec;
+ };
+
+ const T* getcontent() const
+ {
+ return m_Vec;
+ };
+
+
+ int getlowbound(int iBoundNum = 0) const
+ {
+ return m_iLow;
+ };
+
+
+ int gethighbound(int iBoundNum = 0) const
+ {
+ return m_iHigh;
+ };
+
+ raw_vector<T> getvector(int iStart, int iEnd)
+ {
+ if( iStart>iEnd || wrongIdx(iStart) || wrongIdx(iEnd) )
+ return raw_vector<T>(0, 0, 1);
+ else
+ return raw_vector<T>(m_Vec+iStart-m_iLow, iEnd-iStart+1, 1);
+ };
+
+
+ const_raw_vector<T> getvector(int iStart, int iEnd) const
+ {
+ if( iStart>iEnd || wrongIdx(iStart) || wrongIdx(iEnd) )
+ return const_raw_vector<T>(0, 0, 1);
+ else
+ return const_raw_vector<T>(m_Vec+iStart-m_iLow, iEnd-iStart+1, 1);
+ };
+private:
+ bool wrongIdx(int i) const { return i<m_iLow || i>m_iHigh; };
+
+ T *m_Vec;
+ long m_iVecSize;
+ long m_iLow, m_iHigh;
+};
+
+
+
+/********************************************************************
+Template of a dynamical two-dimensional array
+********************************************************************/
+template<class T, bool Aligned = false>
+class template_2d_array
+{
+public:
+ template_2d_array()
+ {
+ m_Vec=0;
+ m_iVecSize=0;
+ m_iLow1 = 0;
+ m_iHigh1 = -1;
+ m_iLow2 = 0;
+ m_iHigh2 = -1;
+ };
+
+ ~template_2d_array()
+ {
+ if(m_Vec)
+ {
+ if( Aligned )
+ ap::afree(m_Vec);
+ else
+ delete[] m_Vec;
+ }
+ };
+
+ template_2d_array(const template_2d_array &rhs)
+ {
+ m_Vec=0;
+ m_iVecSize=0;
+ m_iLow1 = 0;
+ m_iHigh1 = -1;
+ m_iLow2 = 0;
+ m_iHigh2 = -1;
+ if( rhs.m_iVecSize!=0 )
+ {
+ setbounds(rhs.m_iLow1, rhs.m_iHigh1, rhs.m_iLow2, rhs.m_iHigh2);
+ for(int i=m_iLow1; i<=m_iHigh1; i++)
+ vmove(&(operator()(i,m_iLow2)), &(rhs(i,m_iLow2)), m_iHigh2-m_iLow2+1);
+ }
+ };
+ const template_2d_array& operator=(const template_2d_array &rhs)
+ {
+ if( this==&rhs )
+ return *this;
+
+ if( rhs.m_iVecSize!=0 )
+ {
+ setbounds(rhs.m_iLow1, rhs.m_iHigh1, rhs.m_iLow2, rhs.m_iHigh2);
+ for(int i=m_iLow1; i<=m_iHigh1; i++)
+ vmove(&(operator()(i,m_iLow2)), &(rhs(i,m_iLow2)), m_iHigh2-m_iLow2+1);
+ }
+ else
+ {
+ m_Vec=0;
+ m_iVecSize=0;
+ m_iLow1 = 0;
+ m_iHigh1 = -1;
+ m_iLow2 = 0;
+ m_iHigh2 = -1;
+ }
+ return *this;
+ };
+
+ const T& operator()(int i1, int i2) const
+ {
+ #ifndef NO_AP_ASSERT
+ ap_error::make_assertion(i1>=m_iLow1 && i1<=m_iHigh1);
+ ap_error::make_assertion(i2>=m_iLow2 && i2<=m_iHigh2);
+ #endif
+ return m_Vec[ m_iConstOffset + i2 +i1*m_iLinearMember];
+ };
+
+ T& operator()(int i1, int i2)
+ {
+ #ifndef NO_AP_ASSERT
+ ap_error::make_assertion(i1>=m_iLow1 && i1<=m_iHigh1);
+ ap_error::make_assertion(i2>=m_iLow2 && i2<=m_iHigh2);
+ #endif
+ return m_Vec[ m_iConstOffset + i2 +i1*m_iLinearMember];
+ };
+
+ void setbounds( int iLow1, int iHigh1, int iLow2, int iHigh2 )
+ {
+ if(m_Vec)
+ {
+ if( Aligned )
+ ap::afree(m_Vec);
+ else
+ delete[] m_Vec;
+ }
+ int n1 = iHigh1-iLow1+1;
+ int n2 = iHigh2-iLow2+1;
+ m_iVecSize = n1*n2;
+ if( Aligned )
+ {
+ //if( n2%2!=0 )
+ while( (n2*sizeof(T))%16!=0 )
+ {
+ n2++;
+ m_iVecSize += n1;
+ }
+ m_Vec = (T*)ap::amalloc(m_iVecSize*sizeof(T), 16);
+ }
+ else
+ m_Vec = new T[m_iVecSize];
+ m_iLow1 = iLow1;
+ m_iHigh1 = iHigh1;
+ m_iLow2 = iLow2;
+ m_iHigh2 = iHigh2;
+ m_iConstOffset = -m_iLow2-m_iLow1*n2;
+ m_iLinearMember = n2;
+ };
+
+ void setcontent( int iLow1, int iHigh1, int iLow2, int iHigh2, const T *pContent )
+ {
+ setbounds(iLow1, iHigh1, iLow2, iHigh2);
+ for(int i=m_iLow1; i<=m_iHigh1; i++, pContent += m_iHigh2-m_iLow2+1)
+ vmove(&(operator()(i,m_iLow2)), pContent, m_iHigh2-m_iLow2+1);
+ };
+
+ int getlowbound(int iBoundNum) const
+ {
+ return iBoundNum==1 ? m_iLow1 : m_iLow2;
+ };
+
+ int gethighbound(int iBoundNum) const
+ {
+ return iBoundNum==1 ? m_iHigh1 : m_iHigh2;
+ };
+
+ raw_vector<T> getcolumn(int iColumn, int iRowStart, int iRowEnd)
+ {
+ if( (iRowStart>iRowEnd) || wrongColumn(iColumn) || wrongRow(iRowStart) ||wrongRow(iRowEnd) )
+ return raw_vector<T>(0, 0, 1);
+ else
+ return raw_vector<T>(&((*this)(iRowStart, iColumn)), iRowEnd-iRowStart+1, m_iLinearMember);
+ };
+
+ raw_vector<T> getrow(int iRow, int iColumnStart, int iColumnEnd)
+ {
+ if( (iColumnStart>iColumnEnd) || wrongRow(iRow) || wrongColumn(iColumnStart) || wrongColumn(iColumnEnd))
+ return raw_vector<T>(0, 0, 1);
+ else
+ return raw_vector<T>(&((*this)(iRow, iColumnStart)), iColumnEnd-iColumnStart+1, 1);
+ };
+
+ const_raw_vector<T> getcolumn(int iColumn, int iRowStart, int iRowEnd) const
+ {
+ if( (iRowStart>iRowEnd) || wrongColumn(iColumn) || wrongRow(iRowStart) ||wrongRow(iRowEnd) )
+ return const_raw_vector<T>(0, 0, 1);
+ else
+ return const_raw_vector<T>(&((*this)(iRowStart, iColumn)), iRowEnd-iRowStart+1, m_iLinearMember);
+ };
+
+ const_raw_vector<T> getrow(int iRow, int iColumnStart, int iColumnEnd) const
+ {
+ if( (iColumnStart>iColumnEnd) || wrongRow(iRow) || wrongColumn(iColumnStart) || wrongColumn(iColumnEnd))
+ return const_raw_vector<T>(0, 0, 1);
+ else
+ return const_raw_vector<T>(&((*this)(iRow, iColumnStart)), iColumnEnd-iColumnStart+1, 1);
+ };
+private:
+ bool wrongRow(int i) const { return i<m_iLow1 || i>m_iHigh1; };
+ bool wrongColumn(int j) const { return j<m_iLow2 || j>m_iHigh2; };
+
+ T *m_Vec;
+ long m_iVecSize;
+ long m_iLow1, m_iLow2, m_iHigh1, m_iHigh2;
+ long m_iConstOffset, m_iLinearMember;
+};
+
+
+typedef template_1d_array<int> integer_1d_array;
+typedef template_1d_array<double,true> real_1d_array;
+typedef template_1d_array<complex> complex_1d_array;
+typedef template_1d_array<bool> boolean_1d_array;
+
+typedef template_2d_array<int> integer_2d_array;
+typedef template_2d_array<double,true> real_2d_array;
+typedef template_2d_array<complex> complex_2d_array;
+typedef template_2d_array<bool> boolean_2d_array;
+
+
+/********************************************************************
+Constants and functions introduced for compatibility with AlgoPascal
+********************************************************************/
+extern const double machineepsilon;
+extern const double maxrealnumber;
+extern const double minrealnumber;
+
+int sign(double x);
+double randomreal();
+int randominteger(int maxv);
+int round(double x);
+int trunc(double x);
+int ifloor(double x);
+int iceil(double x);
+double pi();
+double sqr(double x);
+int maxint(int m1, int m2);
+int minint(int m1, int m2);
+double maxreal(double m1, double m2);
+double minreal(double m1, double m2);
+
+};//namespace ap
+
+
+#endif
diff --git a/apvt.h b/apvt.h
new file mode 100755
index 0000000..e98ec60
--- /dev/null
+++ b/apvt.h
@@ -0,0 +1,754 @@
+/********************************************************************
+Templates for AP Library
+Copyright (c) 2003-2008, Sergey Bochkanov (ALGLIB project).
+See www.alglib.net or alglib.sources.ru for details.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+- Redistributions of source code must retain the above copyright
+ notice, this list of conditions and the following disclaimer.
+
+- Redistributions in binary form must reproduce the above copyright
+ notice, this list of conditions and the following disclaimer listed
+ in this license in the documentation and/or other materials
+ provided with the distribution.
+
+- Neither the name of the copyright holders nor the names of its
+ contributors may 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 APVT_H
+#define APVT_H
+
+/********************************************************************
+Template defining vector in memory. It is used by the basic
+subroutines of linear algebra.
+
+Vector consists of Length elements of type T, starting from an element,
+which Data is pointed to. Interval between adjacent elements equals
+the value of Step.
+
+The class provides an access for reading only.
+********************************************************************/
+template<class T>
+class const_raw_vector
+{
+public:
+ const_raw_vector(const T *Data, int Length, int Step):
+ pData(const_cast<T*>(Data)),iLength(Length),iStep(Step){};
+
+ const T* GetData() const
+ { return pData; };
+
+ int GetLength() const
+ { return iLength; };
+
+ int GetStep() const
+ { return iStep; };
+protected:
+ T *pData;
+ int iLength, iStep;
+};
+
+
+/********************************************************************
+Template defining vector in memory, derived from const_raw_vector.
+It is used by the basic subroutines of linear algebra.
+
+Vector consists of Length elements of type T, starting from an element,
+which Data is pointed to. Interval between adjacent elements equals
+the value of Step.
+
+The class provides an access both for reading and writing.
+********************************************************************/
+template<class T>
+class raw_vector : public const_raw_vector<T>
+{
+public:
+ raw_vector(T *Data, int Length, int Step):const_raw_vector<T>(Data, Length, Step){};
+
+ T* GetData()
+ { return const_raw_vector<T>::pData; };
+};
+
+
+/********************************************************************
+Dot product
+********************************************************************/
+template<class T>
+T vdotproduct(const_raw_vector<T> v1, const_raw_vector<T> v2)
+{
+ ap_error::make_assertion(v1.GetLength()==v2.GetLength());
+ if( v1.GetStep()==1 && v2.GetStep()==1 )
+ {
+ //
+ // fast
+ //
+ T r = 0;
+ const T *p1 = v1.GetData();
+ const T *p2 = v2.GetData();
+ int imax = v1.GetLength()/4;
+ int i;
+ for(i=imax; i!=0; i--)
+ {
+ r += (*p1)*(*p2) + p1[1]*p2[1] + p1[2]*p2[2] + p1[3]*p2[3];
+ p1+=4;
+ p2+=4;
+ }
+ for(i=0; i<v1.GetLength()%4; i++)
+ r += (*(p1++))*(*(p2++));
+ return r;
+ }
+ else
+ {
+ //
+ // general
+ //
+ int offset11 = v1.GetStep(), offset12 = 2*offset11, offset13 = 3*offset11, offset14 = 4*offset11;
+ int offset21 = v2.GetStep(), offset22 = 2*offset21, offset23 = 3*offset21, offset24 = 4*offset21;
+ T r = 0;
+ const T *p1 = v1.GetData();
+ const T *p2 = v2.GetData();
+ int imax = v1.GetLength()/4;
+ int i;
+ for(i=0; i<imax; i++)
+ {
+ r += (*p1)*(*p2) + p1[offset11]*p2[offset21] + p1[offset12]*p2[offset22] + p1[offset13]*p2[offset23];
+ p1+=offset14;
+ p2+=offset24;
+ }
+ for(i=0; i<v1.GetLength()%4; i++)
+ {
+ r += (*p1)*(*p2);
+ p1+=offset11;
+ p2+=offset21;
+ }
+ return r;
+ }
+}
+
+
+/********************************************************************
+Dot product, another form
+********************************************************************/
+template<class T>
+T _vdotproduct(const T *v1, const T *v2, int N)
+{
+ T r = 0;
+ int imax = N/4;
+ int i;
+ for(i=imax; i!=0; i--)
+ {
+ r += (*v1)*(*v2) + v1[1]*v2[1] + v1[2]*v2[2] + v1[3]*v2[3];
+ v1+=4;
+ v2+=4;
+ }
+ for(i=0; i<N%4; i++)
+ r+=(*(v1++))*(*(v2++));
+ return r;
+}
+
+
+/********************************************************************
+Copy one vector into another
+********************************************************************/
+template<class T>
+void vmove(raw_vector<T> vdst, const_raw_vector<T> vsrc)
+{
+ ap_error::make_assertion(vdst.GetLength()==vsrc.GetLength());
+ if( vdst.GetStep()==1 && vsrc.GetStep()==1 )
+ {
+ //
+ // fast
+ //
+ T *p1 = vdst.GetData();
+ const T *p2 = vsrc.GetData();
+ int imax = vdst.GetLength()/2;
+ int i;
+ for(i=imax; i!=0; i--)
+ {
+ *p1 = *p2;
+ p1[1] = p2[1];
+ p1 += 2;
+ p2 += 2;
+ }
+ if(vdst.GetLength()%2 != 0)
+ *p1 = *p2;
+ return;
+ }
+ else
+ {
+ //
+ // general
+ //
+ int offset11 = vdst.GetStep(), offset12 = 2*offset11, offset13 = 3*offset11, offset14 = 4*offset11;
+ int offset21 = vsrc.GetStep(), offset22 = 2*offset21, offset23 = 3*offset21, offset24 = 4*offset21;
+ T *p1 = vdst.GetData();
+ const T *p2 = vsrc.GetData();
+ int imax = vdst.GetLength()/4;
+ int i;
+ for(i=0; i<imax; i++)
+ {
+ *p1 = *p2;
+ p1[offset11] = p2[offset21];
+ p1[offset12] = p2[offset22];
+ p1[offset13] = p2[offset23];
+ p1 += offset14;
+ p2 += offset24;
+ }
+ for(i=0; i<vdst.GetLength()%4; i++)
+ {
+ *p1 = *p2;
+ p1 += vdst.GetStep();
+ p2 += vsrc.GetStep();
+ }
+ return;
+ }
+}
+
+
+/********************************************************************
+vmove, abother form
+********************************************************************/
+template<class T>
+void _vmove(T *vdst, const T* vsrc, int N)
+{
+ int imax = N/2;
+ int i;
+ for(i=imax; i!=0; i--)
+ {
+ *vdst = *vsrc;
+ vdst[1] = vsrc[1];
+ vdst += 2;
+ vsrc += 2;
+ }
+ if(N%2!=0)
+ *vdst = *vsrc;
+}
+
+
+/********************************************************************
+Copy one vector multiplied by -1 into another.
+********************************************************************/
+template<class T>
+void vmoveneg(raw_vector<T> vdst, const_raw_vector<T> vsrc)
+{
+ ap_error::make_assertion(vdst.GetLength()==vsrc.GetLength());
+ if( vdst.GetStep()==1 && vsrc.GetStep()==1 )
+ {
+ //
+ // fast
+ //
+ T *p1 = vdst.GetData();
+ const T *p2 = vsrc.GetData();
+ int imax = vdst.GetLength()/2;
+ int i;
+ for(i=0; i<imax; i++)
+ {
+ *p1 = -*p2;
+ p1[1] = -p2[1];
+ p1 += 2;
+ p2 += 2;
+ }
+ if(vdst.GetLength()%2 != 0)
+ *p1 = -*p2;
+ return;
+ }
+ else
+ {
+ //
+ // general
+ //
+ int offset11 = vdst.GetStep(), offset12 = 2*offset11, offset13 = 3*offset11, offset14 = 4*offset11;
+ int offset21 = vsrc.GetStep(), offset22 = 2*offset21, offset23 = 3*offset21, offset24 = 4*offset21;
+ T *p1 = vdst.GetData();
+ const T *p2 = vsrc.GetData();
+ int imax = vdst.GetLength()/4;
+ int i;
+ for(i=imax; i!=0; i--)
+ {
+ *p1 = -*p2;
+ p1[offset11] = -p2[offset21];
+ p1[offset12] = -p2[offset22];
+ p1[offset13] = -p2[offset23];
+ p1 += offset14;
+ p2 += offset24;
+ }
+ for(i=0; i<vdst.GetLength()%4; i++)
+ {
+ *p1 = -*p2;
+ p1 += vdst.GetStep();
+ p2 += vsrc.GetStep();
+ }
+ return;
+ }
+}
+
+
+/********************************************************************
+vmoveneg, another form
+********************************************************************/
+template<class T>
+void _vmoveneg(T *vdst, const T *vsrc, int N)
+{
+ T *p1 = vdst;
+ const T *p2 = vsrc;
+ int imax = N/2;
+ int i;
+ for(i=0; i<imax; i++)
+ {
+ *p1 = -*p2;
+ p1[1] = -p2[1];
+ p1 += 2;
+ p2 += 2;
+ }
+ if(N%2 != 0)
+ *p1 = -*p2;
+}
+
+
+/********************************************************************
+Copy one vector multiplied by a number into another vector.
+********************************************************************/
+template<class T, class T2>
+void vmove(raw_vector<T> vdst, const_raw_vector<T> vsrc, T2 alpha)
+{
+ ap_error::make_assertion(vdst.GetLength()==vsrc.GetLength());
+ if( vdst.GetStep()==1 && vsrc.GetStep()==1 )
+ {
+ //
+ // fast
+ //
+ T *p1 = vdst.GetData();
+ const T *p2 = vsrc.GetData();
+ int imax = vdst.GetLength()/4;
+ int i;
+ for(i=imax; i!=0; i--)
+ {
+ *p1 = alpha*(*p2);
+ p1[1] = alpha*p2[1];
+ p1[2] = alpha*p2[2];
+ p1[3] = alpha*p2[3];
+ p1 += 4;
+ p2 += 4;
+ }
+ for(i=0; i<vdst.GetLength()%4; i++)
+ *(p1++) = alpha*(*(p2++));
+ return;
+ }
+ else
+ {
+ //
+ // general
+ //
+ int offset11 = vdst.GetStep(), offset12 = 2*offset11, offset13 = 3*offset11, offset14 = 4*offset11;
+ int offset21 = vsrc.GetStep(), offset22 = 2*offset21, offset23 = 3*offset21, offset24 = 4*offset21;
+ T *p1 = vdst.GetData();
+ const T *p2 = vsrc.GetData();
+ int imax = vdst.GetLength()/4;
+ int i;
+ for(i=0; i<imax; i++)
+ {
+ *p1 = alpha*(*p2);
+ p1[offset11] = alpha*p2[offset21];
+ p1[offset12] = alpha*p2[offset22];
+ p1[offset13] = alpha*p2[offset23];
+ p1 += offset14;
+ p2 += offset24;
+ }
+ for(i=0; i<vdst.GetLength()%4; i++)
+ {
+ *p1 = alpha*(*p2);
+ p1 += vdst.GetStep();
+ p2 += vsrc.GetStep();
+ }
+ return;
+ }
+}
+
+
+/********************************************************************
+vmove, another form
+********************************************************************/
+template<class T, class T2>
+void _vmove(T *vdst, const T *vsrc, int N, T2 alpha)
+{
+ T *p1 = vdst;
+ const T *p2 = vsrc;
+ int imax = N/4;
+ int i;
+ for(i=imax; i!=0; i--)
+ {
+ *p1 = alpha*(*p2);
+ p1[1] = alpha*p2[1];
+ p1[2] = alpha*p2[2];
+ p1[3] = alpha*p2[3];
+ p1 += 4;
+ p2 += 4;
+ }
+ for(i=0; i<N%4; i++)
+ *(p1++) = alpha*(*(p2++));
+}
+
+
+/********************************************************************
+Vector addition
+********************************************************************/
+template<class T>
+void vadd(raw_vector<T> vdst, const_raw_vector<T> vsrc)
+{
+ ap_error::make_assertion(vdst.GetLength()==vsrc.GetLength());
+ if( vdst.GetStep()==1 && vsrc.GetStep()==1 )
+ {
+ //
+ // fast
+ //
+ T *p1 = vdst.GetData();
+ const T *p2 = vsrc.GetData();
+ int imax = vdst.GetLength()/4;
+ int i;
+ for(i=imax; i!=0; i--)
+ {
+ *p1 += *p2;
+ p1[1] += p2[1];
+ p1[2] += p2[2];
+ p1[3] += p2[3];
+ p1 += 4;
+ p2 += 4;
+ }
+ for(i=0; i<vdst.GetLength()%4; i++)
+ *(p1++) += *(p2++);
+ return;
+ }
+ else
+ {
+ //
+ // general
+ //
+ int offset11 = vdst.GetStep(), offset12 = 2*offset11, offset13 = 3*offset11, offset14 = 4*offset11;
+ int offset21 = vsrc.GetStep(), offset22 = 2*offset21, offset23 = 3*offset21, offset24 = 4*offset21;
+ T *p1 = vdst.GetData();
+ const T *p2 = vsrc.GetData();
+ int imax = vdst.GetLength()/4;
+ int i;
+ for(i=0; i<imax; i++)
+ {
+ *p1 += *p2;
+ p1[offset11] += p2[offset21];
+ p1[offset12] += p2[offset22];
+ p1[offset13] += p2[offset23];
+ p1 += offset14;
+ p2 += offset24;
+ }
+ for(i=0; i<vdst.GetLength()%4; i++)
+ {
+ *p1 += *p2;
+ p1 += vdst.GetStep();
+ p2 += vsrc.GetStep();
+ }
+ return;
+ }
+}
+
+
+/********************************************************************
+vadd, another form
+********************************************************************/
+template<class T>
+void _vadd(T *vdst, const T *vsrc, int N)
+{
+ T *p1 = vdst;
+ const T *p2 = vsrc;
+ int imax = N/4;
+ int i;
+ for(i=imax; i!=0; i--)
+ {
+ *p1 += *p2;
+ p1[1] += p2[1];
+ p1[2] += p2[2];
+ p1[3] += p2[3];
+ p1 += 4;
+ p2 += 4;
+ }
+ for(i=0; i<N%4; i++)
+ *(p1++) += *(p2++);
+}
+
+
+/********************************************************************
+Add one vector multiplied by a number to another vector.
+********************************************************************/
+template<class T, class T2>
+void vadd(raw_vector<T> vdst, const_raw_vector<T> vsrc, T2 alpha)
+{
+ ap_error::make_assertion(vdst.GetLength()==vsrc.GetLength());
+ if( vdst.GetStep()==1 && vsrc.GetStep()==1 )
+ {
+ //
+ // fast
+ //
+ T *p1 = vdst.GetData();
+ const T *p2 = vsrc.GetData();
+ int imax = vdst.GetLength()/4;
+ int i;
+ for(i=imax; i!=0; i--)
+ {
+ *p1 += alpha*(*p2);
+ p1[1] += alpha*p2[1];
+ p1[2] += alpha*p2[2];
+ p1[3] += alpha*p2[3];
+ p1 += 4;
+ p2 += 4;
+ }
+ for(i=0; i<vdst.GetLength()%4; i++)
+ *(p1++) += alpha*(*(p2++));
+ return;
+ }
+ else
+ {
+ //
+ // general
+ //
+ int offset11 = vdst.GetStep(), offset12 = 2*offset11, offset13 = 3*offset11, offset14 = 4*offset11;
+ int offset21 = vsrc.GetStep(), offset22 = 2*offset21, offset23 = 3*offset21, offset24 = 4*offset21;
+ T *p1 = vdst.GetData();
+ const T *p2 = vsrc.GetData();
+ int imax = vdst.GetLength()/4;
+ int i;
+ for(i=0; i<imax; i++)
+ {
+ *p1 += alpha*(*p2);
+ p1[offset11] += alpha*p2[offset21];
+ p1[offset12] += alpha*p2[offset22];
+ p1[offset13] += alpha*p2[offset23];
+ p1 += offset14;
+ p2 += offset24;
+ }
+ for(i=0; i<vdst.GetLength()%4; i++)
+ {
+ *p1 += alpha*(*p2);
+ p1 += vdst.GetStep();
+ p2 += vsrc.GetStep();
+ }
+ return;
+ }
+}
+
+
+/********************************************************************
+vadd, another form
+********************************************************************/
+template<class T, class T2>
+void _vadd(T *vdst, const T *vsrc, int N, T2 alpha)
+{
+ T *p1 = vdst;
+ const T *p2 = vsrc;
+ int imax = N/4;
+ int i;
+ for(i=imax; i!=0; i--)
+ {
+ *p1 += alpha*(*p2);
+ p1[1] += alpha*p2[1];
+ p1[2] += alpha*p2[2];
+ p1[3] += alpha*p2[3];
+ p1 += 4;
+ p2 += 4;
+ }
+ for(i=0; i<N%4; i++)
+ *(p1++) += alpha*(*(p2++));
+}
+
+
+/********************************************************************
+Vector subtraction
+********************************************************************/
+template<class T>
+void vsub(raw_vector<T> vdst, const_raw_vector<T> vsrc)
+{
+ ap_error::make_assertion(vdst.GetLength()==vsrc.GetLength());
+ if( vdst.GetStep()==1 && vsrc.GetStep()==1 )
+ {
+ //
+ // fast
+ //
+ T *p1 = vdst.GetData();
+ const T *p2 = vsrc.GetData();
+ int imax = vdst.GetLength()/4;
+ int i;
+ for(i=imax; i!=0; i--)
+ {
+ *p1 -= *p2;
+ p1[1] -= p2[1];
+ p1[2] -= p2[2];
+ p1[3] -= p2[3];
+ p1 += 4;
+ p2 += 4;
+ }
+ for(i=0; i<vdst.GetLength()%4; i++)
+ *(p1++) -= *(p2++);
+ return;
+ }
+ else
+ {
+ //
+ // general
+ //
+ int offset11 = vdst.GetStep(), offset12 = 2*offset11, offset13 = 3*offset11, offset14 = 4*offset11;
+ int offset21 = vsrc.GetStep(), offset22 = 2*offset21, offset23 = 3*offset21, offset24 = 4*offset21;
+ T *p1 = vdst.GetData();
+ const T *p2 = vsrc.GetData();
+ int imax = vdst.GetLength()/4;
+ int i;
+ for(i=0; i<imax; i++)
+ {
+ *p1 -= *p2;
+ p1[offset11] -= p2[offset21];
+ p1[offset12] -= p2[offset22];
+ p1[offset13] -= p2[offset23];
+ p1 += offset14;
+ p2 += offset24;
+ }
+ for(i=0; i<vdst.GetLength()%4; i++)
+ {
+ *p1 -= *p2;
+ p1 += vdst.GetStep();
+ p2 += vsrc.GetStep();
+ }
+ return;
+ }
+}
+
+
+/********************************************************************
+vsub, another form
+********************************************************************/
+template<class T>
+void _vsub(T *vdst, const T *vsrc, int N)
+{
+ T *p1 = vdst;
+ const T *p2 = vsrc;
+ int imax = N/4;
+ int i;
+ for(i=imax; i!=0; i--)
+ {
+ *p1 -= *p2;
+ p1[1] -= p2[1];
+ p1[2] -= p2[2];
+ p1[3] -= p2[3];
+ p1 += 4;
+ p2 += 4;
+ }
+ for(i=0; i<N%4; i++)
+ *(p1++) -= *(p2++);
+}
+
+
+/********************************************************************
+Subtract one vector multiplied by a number from another vector.
+********************************************************************/
+template<class T, class T2>
+void vsub(raw_vector<T> vdst, const_raw_vector<T> vsrc, T2 alpha)
+{
+ vadd(vdst, vsrc, -alpha);
+}
+
+
+/********************************************************************
+vsub, another form
+********************************************************************/
+template<class T, class T2>
+void _vsub(T *vdst, const T *vsrc, int N, T2 alpha)
+{
+ _vadd(vdst, vsrc, N, -alpha);
+}
+
+
+/********************************************************************
+In-place vector multiplication
+********************************************************************/
+template<class T, class T2>
+void vmul(raw_vector<T> vdst, T2 alpha)
+{
+ if( vdst.GetStep()==1 )
+ {
+ //
+ // fast
+ //
+ T *p1 = vdst.GetData();
+ int imax = vdst.GetLength()/4;
+ int i;
+ for(i=imax; i!=0; i--)
+ {
+ *p1 *= alpha;
+ p1[1] *= alpha;
+ p1[2] *= alpha;
+ p1[3] *= alpha;
+ p1 += 4;
+ }
+ for(i=0; i<vdst.GetLength()%4; i++)
+ *(p1++) *= alpha;
+ return;
+ }
+ else
+ {
+ //
+ // general
+ //
+ int offset11 = vdst.GetStep(), offset12 = 2*offset11, offset13 = 3*offset11, offset14 = 4*offset11;
+ T *p1 = vdst.GetData();
+ int imax = vdst.GetLength()/4;
+ int i;
+ for(i=0; i<imax; i++)
+ {
+ *p1 *= alpha;
+ p1[offset11] *= alpha;
+ p1[offset12] *= alpha;
+ p1[offset13] *= alpha;
+ p1 += offset14;
+ }
+ for(i=0; i<vdst.GetLength()%4; i++)
+ {
+ *p1 *= alpha;
+ p1 += vdst.GetStep();
+ }
+ return;
+ }
+}
+
+
+/********************************************************************
+vmul, another form
+********************************************************************/
+template<class T, class T2>
+void _vmul(T *vdst, int N, T2 alpha)
+{
+ T *p1 = vdst;
+ int imax = N/4;
+ int i;
+ for(i=imax; i!=0; i--)
+ {
+ *p1 *= alpha;
+ p1[1] *= alpha;
+ p1[2] *= alpha;
+ p1[3] *= alpha;
+ p1 += 4;
+ }
+ for(i=0; i<N%4; i++)
+ *(p1++) *= alpha;
+}
+
+#endif
diff --git a/chisquaredistr.cpp b/chisquaredistr.cpp
new file mode 100755
index 0000000..b6d9bb6
--- /dev/null
+++ b/chisquaredistr.cpp
@@ -0,0 +1,157 @@
+/*************************************************************************
+Cephes Math Library Release 2.8: June, 2000
+Copyright 1984, 1987, 2000 by Stephen L. Moshier
+
+Contributors:
+ * Sergey Bochkanov (ALGLIB project). Translation from C to
+ pseudocode.
+
+See subroutines comments for additional copyrights.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+- Redistributions of source code must retain the above copyright
+ notice, this list of conditions and the following disclaimer.
+
+- Redistributions in binary form must reproduce the above copyright
+ notice, this list of conditions and the following disclaimer listed
+ in this license in the documentation and/or other materials
+ provided with the distribution.
+
+- Neither the name of the copyright holders nor the names of its
+ contributors may 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 <stdafx.h>
+#include "chisquaredistr.h"
+
+/*************************************************************************
+Chi-square distribution
+
+Returns the area under the left hand tail (from 0 to x)
+of the Chi square probability density function with
+v degrees of freedom.
+
+
+ x
+ -
+ 1 | | v/2-1 -t/2
+ P( x | v ) = ----------- | t e dt
+ v/2 - | |
+ 2 | (v/2) -
+ 0
+
+where x is the Chi-square variable.
+
+The incomplete gamma integral is used, according to the
+formula
+
+y = chdtr( v, x ) = igam( v/2.0, x/2.0 ).
+
+The arguments must both be positive.
+
+ACCURACY:
+
+See incomplete gamma function
+
+
+Cephes Math Library Release 2.8: June, 2000
+Copyright 1984, 1987, 2000 by Stephen L. Moshier
+*************************************************************************/
+double chisquaredistribution(double v, double x)
+{
+ double result;
+
+ ap::ap_error::make_assertion(x>=0&&v>=1, "Domain error in ChiSquareDistribution");
+ result = incompletegamma(v/2.0, x/2.0);
+ return result;
+}
+
+
+/*************************************************************************
+Complemented Chi-square distribution
+
+Returns the area under the right hand tail (from x to
+infinity) of the Chi square probability density function
+with v degrees of freedom:
+
+ inf.
+ -
+ 1 | | v/2-1 -t/2
+ P( x | v ) = ----------- | t e dt
+ v/2 - | |
+ 2 | (v/2) -
+ x
+
+where x is the Chi-square variable.
+
+The incomplete gamma integral is used, according to the
+formula
+
+y = chdtr( v, x ) = igamc( v/2.0, x/2.0 ).
+
+The arguments must both be positive.
+
+ACCURACY:
+
+See incomplete gamma function
+
+Cephes Math Library Release 2.8: June, 2000
+Copyright 1984, 1987, 2000 by Stephen L. Moshier
+*************************************************************************/
+double chisquarecdistribution(double v, double x)
+{
+ double result;
+
+ ap::ap_error::make_assertion(x>=0&&v>=1, "Domain error in ChiSquareDistributionC");
+ result = incompletegammac(v/2.0, x/2.0);
+ return result;
+}
+
+
+/*************************************************************************
+Inverse of complemented Chi-square distribution
+
+Finds the Chi-square argument x such that the integral
+from x to infinity of the Chi-square density is equal
+to the given cumulative probability y.
+
+This is accomplished using the inverse gamma integral
+function and the relation
+
+ x/2 = igami( df/2, y );
+
+ACCURACY:
+
+See inverse incomplete gamma function
+
+
+Cephes Math Library Release 2.8: June, 2000
+Copyright 1984, 1987, 2000 by Stephen L. Moshier
+*************************************************************************/
+double invchisquaredistribution(double v, double y)
+{
+ double result;
+
+ ap::ap_error::make_assertion(y>=0&&y<=1&&v>=1, "Domain error in InvChiSquareDistribution");
+ result = 2*invincompletegammac(0.5*v, y);
+ return result;
+}
+
+
+
diff --git a/chisquaredistr.h b/chisquaredistr.h
new file mode 100755
index 0000000..c8196d7
--- /dev/null
+++ b/chisquaredistr.h
@@ -0,0 +1,143 @@
+/*************************************************************************
+Cephes Math Library Release 2.8: June, 2000
+Copyright 1984, 1987, 2000 by Stephen L. Moshier
+
+Contributors:
+ * Sergey Bochkanov (ALGLIB project). Translation from C to
+ pseudocode.
+
+See subroutines comments for additional copyrights.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+- Redistributions of source code must retain the above copyright
+ notice, this list of conditions and the following disclaimer.
+
+- Redistributions in binary form must reproduce the above copyright
+ notice, this list of conditions and the following disclaimer listed
+ in this license in the documentation and/or other materials
+ provided with the distribution.
+
+- Neither the name of the copyright holders nor the names of its
+ contributors may 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 _chisquaredistr_h
+#define _chisquaredistr_h
+
+#include "ap.h"
+
+#include "gammaf.h"
+#include "normaldistr.h"
+#include "igammaf.h"
+
+
+/*************************************************************************
+Chi-square distribution
+
+Returns the area under the left hand tail (from 0 to x)
+of the Chi square probability density function with
+v degrees of freedom.
+
+
+ x
+ -
+ 1 | | v/2-1 -t/2
+ P( x | v ) = ----------- | t e dt
+ v/2 - | |
+ 2 | (v/2) -
+ 0
+
+where x is the Chi-square variable.
+
+The incomplete gamma integral is used, according to the
+formula
+
+y = chdtr( v, x ) = igam( v/2.0, x/2.0 ).
+
+The arguments must both be positive.
+
+ACCURACY:
+
+See incomplete gamma function
+
+
+Cephes Math Library Release 2.8: June, 2000
+Copyright 1984, 1987, 2000 by Stephen L. Moshier
+*************************************************************************/
+double chisquaredistribution(double v, double x);
+
+
+/*************************************************************************
+Complemented Chi-square distribution
+
+Returns the area under the right hand tail (from x to
+infinity) of the Chi square probability density function
+with v degrees of freedom:
+
+ inf.
+ -
+ 1 | | v/2-1 -t/2
+ P( x | v ) = ----------- | t e dt
+ v/2 - | |
+ 2 | (v/2) -
+ x
+
+where x is the Chi-square variable.
+
+The incomplete gamma integral is used, according to the
+formula
+
+y = chdtr( v, x ) = igamc( v/2.0, x/2.0 ).
+
+The arguments must both be positive.
+
+ACCURACY:
+
+See incomplete gamma function
+
+Cephes Math Library Release 2.8: June, 2000
+Copyright 1984, 1987, 2000 by Stephen L. Moshier
+*************************************************************************/
+double chisquarecdistribution(double v, double x);
+
+
+/*************************************************************************
+Inverse of complemented Chi-square distribution
+
+Finds the Chi-square argument x such that the integral
+from x to infinity of the Chi-square density is equal
+to the given cumulative probability y.
+
+This is accomplished using the inverse gamma integral
+function and the relation
+
+ x/2 = igami( df/2, y );
+
+ACCURACY:
+
+See inverse incomplete gamma function
+
+
+Cephes Math Library Release 2.8: June, 2000
+Copyright 1984, 1987, 2000 by Stephen L. Moshier
+*************************************************************************/
+double invchisquaredistribution(double v, double y);
+
+
+#endif
diff --git a/cohort.cpp b/cohort.cpp
new file mode 100755
index 0000000..c8a3c85
--- /dev/null
+++ b/cohort.cpp
@@ -0,0 +1,47 @@
+/*************************************************************************
+GWAMA software: May, 2010
+
+Contributors:
+ * Andrew P Morris amorris at well.ox.ac.uk
+ * Reedik Magi reedik at well.ox.ac.uk
+
+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 "cohort.h"
+#include <iostream>
+#include <stdlib.h>
+using namespace std;
+cohort::cohort(void)
+{
+ _name = "N";
+ _columnMarkerName = -9;
+ _columnEffectAllele = -9;
+ _columnOtherAllele = -9;
+ _columnEffectAlleleFreq = -9;
+ _columnStrand = -9;
+ _columnBeta = -9;
+ _columnSE = -9;
+ _columnOR = -9;
+ _columnOR_95L = -9;
+ _columnOR_95U = -9;
+ _columnN = -9;
+ _columnImputed = -9;
+ _directLambda = 1;
+ _imputedLambda = 1;
+ _directCount = 0;
+ _imputedCount = 0;
+ _columnCount = -9;
+}
+
+
diff --git a/cohort.h b/cohort.h
new file mode 100755
index 0000000..29bc10e
--- /dev/null
+++ b/cohort.h
@@ -0,0 +1,52 @@
+/*************************************************************************
+GWAMA software: May, 2010
+
+Contributors:
+ * Andrew P Morris amorris at well.ox.ac.uk
+ * Reedik Magi reedik at well.ox.ac.uk
+
+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 _COHORT_H_
+#define _COHORT_H_
+#include <iostream>
+#include <stdlib.h>
+using namespace std;
+class cohort
+{
+public:
+ cohort(void);
+ string _name;
+ int _columnMarkerName;
+ int _columnEffectAllele;
+ int _columnOtherAllele;
+ int _columnEffectAlleleFreq;
+ int _columnStrand;
+ int _columnBeta;
+ int _columnSE;
+ int _columnOR;
+ int _columnOR_95L;
+ int _columnOR_95U;
+ int _columnN;
+ int _columnImputed;
+ int _columnCount;
+ double _directLambda;
+ double _imputedLambda;
+ int _directCount;
+ int _imputedCount;
+
+
+};
+#endif
+
diff --git a/commandLine.cpp b/commandLine.cpp
new file mode 100755
index 0000000..7e0acd4
--- /dev/null
+++ b/commandLine.cpp
@@ -0,0 +1,378 @@
+/*************************************************************************
+GWAMA software: May, 2010
+
+Contributors:
+ * Andrew P Morris amorris at well.ox.ac.uk
+ * Reedik Magi reedik at well.ox.ac.uk
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*************************************************************************/
+
+
+#include <iostream>
+#include <map>
+#include <vector>
+#include <stdlib.h>
+#include <fstream>
+#include "marker.h"
+#include "study.h"
+#include <math.h>
+#include <cctype> // std::toupper
+
+#include "global.h"
+#include "tools.h"
+
+using namespace std;
+
+void printVersion(string version)
+{
+ cout << "GWAMA version: " << version << endl;
+ exit(0);
+}
+void printHelp(string version)
+{
+ string v;
+ double len =17 - version.size()-1;
+ int f = (int) ((len/2.0)+0.6);
+ int r = (int) ((len/2.0));
+ for (int i=0; i < f;i++)v+=" ";
+ v+='v';
+ v+=version;
+ for (int i=0; i < r;i++)v+=" ";
+ cout << endl;
+ cout << "@----------------------------------------------------------@" << endl;
+ cout << "| GWAMA |" << v << "| May, 2010 |" << endl;
+ cout << "|----------------------------------------------------------|" << endl;
+ cout << "| (C) 2008 Reedik Magi & Andrew P Morris |" << endl;
+ cout << "| GNU General Public License, v2 |" << endl;
+ cout << "|----------------------------------------------------------|" << endl;
+ cout << "| For documentation, citation & bug-report instructions: |" << endl;
+ cout << "| http://www.well.ox.ac.uk/GWAMA/ |" << endl;
+ cout << "@----------------------------------------------------------@" << endl;
+ cout << endl;
+ cout << endl;
+ cout << "Command line options:"<<endl;
+ cout << endl;
+ cout << "GWAMA --filelist {filename} or -i {filename} Specify "<<endl;
+ cout << " studies' result files. Default = "<<endl;
+ cout << " gwama.in "<<endl<<endl;
+ cout << " --output {fileroot} or -o {fileroot} Specify file "<<endl;
+ cout << " root for output of analysis. Default "<<endl;
+ cout << " <gwama> (gwama.out) "<<endl<<endl;
+ cout << " --random or -r Use random effect correction "<<endl<<endl;;
+ cout << " --genomic_control or -gc Use genomic control for "<<endl;
+ cout << " adjusting studies' result files "<<endl<<endl;
+ cout << " --genomic_control_output or -gco Use genomic "<<endl;
+ cout << " control on meta-analysis summary. "<<endl;
+ cout << " (i.e. results of meta- "<<endl;
+ cout << " analysis are corrected for gc) "<<endl<<endl;
+ cout << " --quantitative or -qt Use this option, if trait is"<<endl;
+ cout << " quantitative (columns BETA & SE). "<<endl;
+ cout << " Default is binary trait (columns OR, "<<endl;
+ cout << " OR95_U, OR_95_L) "<<endl<<endl;
+ cout << " --threshold {0-1} or -t {0-1} The p-value threshold "<<endl;
+ cout << " for showing direction. Default = 1 "<<endl<<endl;
+ cout << " --map {filename} or -m {filename} Marker position "<<endl;
+ cout << " file for chromosome and position info"<<endl<<endl;
+ cout << " --no_alleles No allele information has been given."<<endl;
+ cout << " Expecting always the same EA. "<<endl<<endl;
+ cout << " --indel_alleles Allele labes might contain more "<<endl;
+ cout << " than single letter. No strand checks."<<endl<<endl;
+ cout << " --sex Run gender-differentiated and gender-"<<endl;
+ cout << " heterogeneity analysis. Gender info "<<endl;
+ cout << " must be provided in filelist file. "<<endl;
+ cout << " (second column after file names is "<<endl;
+ cout << " either M or F) More info in tutorial."<<endl<<endl;
+ cout << " --name_marker alternative header to marker name col"<<endl<<endl;
+ cout << " --name_strand alternative header to strand column "<<endl<<endl;
+ cout << " --name_n alternative header to sample size col"<<endl<<endl;
+ cout << " --name_eaf alternative header to EAF column "<<endl<<endl;
+ cout << " --name_beta alternative header to beta column "<<endl<<endl;
+ cout << " --name_se alternative header to std. err. col "<<endl<<endl;
+ cout << " --name_or alternative header to OR column "<<endl<<endl;
+ cout << " --name_or_95l alternative header to OR 95L column "<<endl<<endl;
+ cout << " --name_or_95u alternative header to OR 95U column "<<endl<<endl;
+ cout << " --help or -h Print this help "<<endl<<endl;
+ cout << " --version or -v Print GWAMA version number "<<endl<<endl;
+ exit (0);
+}
+
+
+bool
+readCommandline (int argc, char *argv[], global & _g)
+{
+ for (int i = 1; i < argc; i++)
+ {
+ if (string(argv[i]).compare("--version")==0 || string(argv[i]).compare("-v")==0)
+ {
+ printVersion(_g._version);
+ return 1;
+ }
+ else if (string(argv[i]).compare("--help")==0 || string(argv[i]).compare("-h")==0)
+ {
+ printHelp(_g._version);
+ return 1;
+ }
+ else if (string(argv[i]).compare("--genomic_control")==0 || string(argv[i]).compare("-gc")==0)
+ {
+ _g._genomicControl=true;
+ cout << "Genomic control enabled" << endl;
+ }
+ else if (string(argv[i]).compare("--no_alleles")==0)
+ {
+ _g._noAlleles=true;
+ cout << "All effects are expected to be according to the same allele." << endl;
+ }
+
+ else if (string(argv[i]).compare("--random")==0 || string(argv[i]).compare("-r")==0)
+ {
+ _g._randomEffect=true;
+ cout << "Using random effect correction" << endl;
+ }
+ else if (string(argv[i]).compare("--quantitative")==0 || string(argv[i]).compare("-qt")==0)
+ {
+ _g._binaryTrait = false;
+ cout << "Quantitative trait (BETA+SE)" << endl;
+
+ }
+ else if (string(argv[i]).compare("--genomic_control_output")==0 || string(argv[i]).compare("-gco")==0)
+ {
+ _g._genomicControlOutput=true;
+ cout << "Genomic control for meta-analysis output enabled" << endl;
+ }
+ else if (string(argv[i]).compare("--sex")==0)
+ {
+ _g._genderHet=true;
+ cout << "Gender-specific heterogeneity analysis enabled" << endl;
+ }
+ else if (string(argv[i]).compare("--indel_alleles")==0)
+ {
+ _g._indel=true;
+ cout << "Indel allele names enabled. No strand checks" << endl;
+ }
+ else if (string(argv[i]).compare("--name_marker")==0)
+ {
+ if (i != argc-1)
+ {
+ cout<<"Alternative marker column name: " << uc(string(argv[i+1])) << endl;;
+ _g._alternativeMarkerNames.push_back(uc(string(argv[i+1])));
+ }
+ if (i == argc-1 )
+ {
+ cout<<"ERROR: Missing a value for --name_marker\n";
+ return 0;
+ }
+ i++;
+ }
+ else if (string(argv[i]).compare("--name_ea")==0)
+ {
+ if (i != argc-1)
+ {
+ cout<<"Alternative effect allele column name: " << uc(string(argv[i+1])) << endl;;
+ _g._alternativeEffectAlleles.push_back(uc(string(argv[i+1])));
+ }
+ if (i == argc-1 )
+ {
+ cout<<"ERROR: Missing a value for --name_ea\n";
+ return 0;
+ }
+ i++;
+ }
+ else if (string(argv[i]).compare("--name_nea")==0)
+ {
+ if (i != argc-1)
+ {
+ cout<<"Alternative other allele column name: " << uc(string(argv[i+1])) << endl;;
+ _g._alternativeOtherAlleles.push_back(uc(string(argv[i+1])));
+ }
+ if (i == argc-1 )
+ {
+ cout<<"ERROR: Missing a value for --name_nea\n";
+ return 0;
+ }
+ i++;
+ }
+ else if (string(argv[i]).compare("--name_or")==0)
+ {
+ if (i != argc-1)
+ {
+ cout<<"Alternative OR column name: " << uc(string(argv[i+1])) << endl;;
+ _g._alternativeORs.push_back(uc(string(argv[i+1])));
+ }
+ if (i == argc-1 )
+ {
+ cout<<"ERROR: Missing a value for --name_or\n";
+ return 0;
+ }
+ i++;
+ }
+ else if (string(argv[i]).compare("--name_or_95l")==0)
+ {
+ if (i != argc-1)
+ {
+ cout<<"Alternative lower CI of OR column name: " << uc(string(argv[i+1])) << endl;;
+ _g._alternativeOR_95Ls.push_back(uc(string(argv[i+1])));
+ }
+ if (i == argc-1 )
+ {
+ cout<<"ERROR: Missing a value for --name_or_95l\n";
+ return 0;
+ }
+ i++;
+ }
+ else if (string(argv[i]).compare("--name_or_95u")==0)
+ {
+ if (i != argc-1)
+ {
+ cout<<"Alternative upper CI of OR column name: " << uc(string(argv[i+1])) << endl;;
+ _g._alternativeOR_95Us.push_back(uc(string(argv[i+1])));
+ }
+ if (i == argc-1 )
+ {
+ cout<<"ERROR: Missing a value for --name_or_95u\n";
+ return 0;
+ }
+ i++;
+ }
+ else if (string(argv[i]).compare("--name_beta")==0)
+ {
+ if (i != argc-1)
+ {
+ cout<<"Alternative beta column name: " << uc(string(argv[i+1])) << endl;;
+ _g._alternativeBetas.push_back(uc(string(argv[i+1])));
+ }
+ if (i == argc-1 )
+ {
+ cout<<"ERROR: Missing a value for --name_beta\n";
+ return 0;
+ }
+ i++;
+ }
+ else if (string(argv[i]).compare("--name_se")==0)
+ {
+ if (i != argc-1)
+ {
+ cout<<"Alternative std. error column name: " << uc(string(argv[i+1])) << endl;;
+ _g._alternativeSEs.push_back(uc(string(argv[i+1])));
+ }
+ if (i == argc-1 )
+ {
+ cout<<"ERROR: Missing a value for --name_se\n";
+ return 0;
+ }
+ i++;
+ }
+ else if (string(argv[i]).compare("--name_n")==0)
+ {
+ if (i != argc-1)
+ {
+ cout<<"Alternative sample size column name: " << uc(string(argv[i+1])) << endl;;
+ _g._alternativeNs.push_back(uc(string(argv[i+1])));
+ }
+ if (i == argc-1 )
+ {
+ cout<<"ERROR: Missing a value for --name_n\n";
+ return 0;
+ }
+ i++;
+ }
+ else if (string(argv[i]).compare("--name_eaf")==0)
+ {
+ if (i != argc-1)
+ {
+ cout<<"Alternative effect allele frequency column name: " << uc(string(argv[i+1])) << endl;;
+ _g._alternativeEffectAlleleFreqs.push_back(uc(string(argv[i+1])));
+ }
+ if (i == argc-1 )
+ {
+ cout<<"ERROR: Missing a value for --name_eaf\n";
+ return 0;
+ }
+ i++;
+ }
+ else if (string(argv[i]).compare("--name_strand")==0)
+ {
+ if (i != argc-1)
+ {
+ cout<<"Alternative strand column name: " << uc(string(argv[i+1])) << endl;;
+ _g._alternativeStrands.push_back(uc(string(argv[i+1])));
+ }
+ if (i == argc-1 )
+ {
+ cout<<"ERROR: Missing a value for --name_strand\n";
+ return 0;
+ }
+ i++;
+ }
+
+ else if (string(argv[i]).compare("--filelist")==0 || string(argv[i]).compare("-i")==0)
+ {
+ if (i != argc-1)
+ {
+ _g._fileList = string(argv[i+1]);
+ }
+ if (i == argc-1 )
+ {
+ cout<<"ERROR: Missing a value for --filelist\n";
+ return 0;
+ }
+ i++;
+ }
+ else if (string(argv[i]).compare("--output")==0 || string(argv[i]).compare("-o")==0)
+ {
+ if (i != argc-1)
+ {
+ _g._outputRoot = string(argv[i+1]);
+ }
+ if (i == argc-1 )
+ {
+ cout<<"ERROR: Missing a value for --output\n";
+ return 0;
+ }
+ i++;
+ }
+ else if (string(argv[i]).compare("--map")==0 || string(argv[i]).compare("-m")==0)
+ {
+ if (i != argc-1)
+ {
+ _g._markerMap = string(argv[i+1]);
+ }
+ if (i == argc-1 )
+ {
+ cout<<"ERROR: Missing a value for --map\n";
+ return 0;
+ }
+ i++;
+ }
+ else if (string(argv[i]).compare("--threshold")==0 || string(argv[i]).compare("-t")==0)
+ {
+ if (i != argc-1)
+ {
+ if (atof(argv[i+1]) >0 && atof(argv[i+1]) <= 1)
+ _g._thresholdPValDir = atof(argv[i+1]);
+ }
+ if (i == argc-1)
+ {
+ cout<<"ERROR: Missing or errogeneous value for --threshold\n";
+ return 0;
+ }
+ i++;
+ }
+ else
+ {
+ cout<<"Unknown command line option: " << string(argv[i]) << ". Exit program.\n";
+ return false;
+ }
+ }
+ return true;
+}
diff --git a/commandLine.h b/commandLine.h
new file mode 100755
index 0000000..d8bdf8f
--- /dev/null
+++ b/commandLine.h
@@ -0,0 +1,40 @@
+/*************************************************************************
+GWAMA software: May, 2010
+
+Contributors:
+ * Andrew P Morris amorris at well.ox.ac.uk
+ * Reedik Magi reedik at well.ox.ac.uk
+
+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 _COMMANDLINE_H_
+#define _COMMANDLINE_H_
+
+
+#include <iostream>
+#include <map>
+#include <vector>
+#include <stdlib.h>
+#include <fstream>
+#include "marker.h"
+#include "study.h"
+#include <math.h>
+#include <cctype> // std::toupper
+#include "global.h"
+using namespace std;
+bool readCommandline(int argc, char *argv[], global & _g);
+void printVersion(string version);
+void printHelp(string version);
+
+#endif
+
diff --git a/gammaf.cpp b/gammaf.cpp
new file mode 100755
index 0000000..060767e
--- /dev/null
+++ b/gammaf.cpp
@@ -0,0 +1,338 @@
+/*************************************************************************
+Cephes Math Library Release 2.8: June, 2000
+Copyright by Stephen L. Moshier
+
+Contributors:
+ * Sergey Bochkanov (ALGLIB project). Translation from C to
+ pseudocode.
+
+See subroutines comments for additional copyrights.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+- Redistributions of source code must retain the above copyright
+ notice, this list of conditions and the following disclaimer.
+
+- Redistributions in binary form must reproduce the above copyright
+ notice, this list of conditions and the following disclaimer listed
+ in this license in the documentation and/or other materials
+ provided with the distribution.
+
+- Neither the name of the copyright holders nor the names of its
+ contributors may 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 <stdafx.h>
+#include "gammaf.h"
+
+double gammastirf(double x);
+
+/*************************************************************************
+Gamma function
+
+Input parameters:
+ X - argument
+
+Domain:
+ 0 < X < 171.6
+ -170 < X < 0, X is not an integer.
+
+Relative error:
+ arithmetic domain # trials peak rms
+ IEEE -170,-33 20000 2.3e-15 3.3e-16
+ IEEE -33, 33 20000 9.4e-16 2.2e-16
+ IEEE 33, 171.6 20000 2.3e-15 3.2e-16
+
+Cephes Math Library Release 2.8: June, 2000
+Original copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier
+Translated to AlgoPascal by Bochkanov Sergey (2005, 2006, 2007).
+*************************************************************************/
+double gamma(double x)
+{
+ double result;
+ double p;
+ double pp;
+ double q;
+ double qq;
+ double z;
+ int i;
+ double sgngam;
+
+ sgngam = 1;
+ q = fabs(x);
+ if( q>33.0 )
+ {
+ if( x<0.0 )
+ {
+ p = ap::ifloor(q);
+ i = ap::round(p);
+ if( i%2==0 )
+ {
+ sgngam = -1;
+ }
+ z = q-p;
+ if( z>0.5 )
+ {
+ p = p+1;
+ z = q-p;
+ }
+ z = q*sin(ap::pi()*z);
+ z = fabs(z);
+ z = ap::pi()/(z*gammastirf(q));
+ }
+ else
+ {
+ z = gammastirf(x);
+ }
+ result = sgngam*z;
+ return result;
+ }
+ z = 1;
+ while(x>=3)
+ {
+ x = x-1;
+ z = z*x;
+ }
+ while(x<0)
+ {
+ if( x>-0.000000001 )
+ {
+ result = z/((1+0.5772156649015329*x)*x);
+ return result;
+ }
+ z = z/x;
+ x = x+1;
+ }
+ while(x<2)
+ {
+ if( x<0.000000001 )
+ {
+ result = z/((1+0.5772156649015329*x)*x);
+ return result;
+ }
+ z = z/x;
+ x = x+1.0;
+ }
+ if( x==2 )
+ {
+ result = z;
+ return result;
+ }
+ x = x-2.0;
+ pp = 1.60119522476751861407E-4;
+ pp = 1.19135147006586384913E-3+x*pp;
+ pp = 1.04213797561761569935E-2+x*pp;
+ pp = 4.76367800457137231464E-2+x*pp;
+ pp = 2.07448227648435975150E-1+x*pp;
+ pp = 4.94214826801497100753E-1+x*pp;
+ pp = 9.99999999999999996796E-1+x*pp;
+ qq = -2.31581873324120129819E-5;
+ qq = 5.39605580493303397842E-4+x*qq;
+ qq = -4.45641913851797240494E-3+x*qq;
+ qq = 1.18139785222060435552E-2+x*qq;
+ qq = 3.58236398605498653373E-2+x*qq;
+ qq = -2.34591795718243348568E-1+x*qq;
+ qq = 7.14304917030273074085E-2+x*qq;
+ qq = 1.00000000000000000320+x*qq;
+ result = z*pp/qq;
+ return result;
+ return result;
+}
+
+
+/*************************************************************************
+Natural logarithm of gamma function
+
+Input parameters:
+ X - argument
+
+Result:
+ logarithm of the absolute value of the Gamma(X).
+
+Output parameters:
+ SgnGam - sign(Gamma(X))
+
+Domain:
+ 0 < X < 2.55e305
+ -2.55e305 < X < 0, X is not an integer.
+
+ACCURACY:
+arithmetic domain # trials peak rms
+ IEEE 0, 3 28000 5.4e-16 1.1e-16
+ IEEE 2.718, 2.556e305 40000 3.5e-16 8.3e-17
+The error criterion was relative when the function magnitude
+was greater than one but absolute when it was less than one.
+
+The following test used the relative error criterion, though
+at certain points the relative error could be much higher than
+indicated.
+ IEEE -200, -4 10000 4.8e-16 1.3e-16
+
+Cephes Math Library Release 2.8: June, 2000
+Copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier
+Translated to AlgoPascal by Bochkanov Sergey (2005, 2006, 2007).
+*************************************************************************/
+double lngamma(double x, double& sgngam)
+{
+ double result;
+ double a;
+ double b;
+ double c;
+ double p;
+ double q;
+ double u;
+ double w;
+ double z;
+ int i;
+ double logpi;
+ double ls2pi;
+ double tmp;
+
+ sgngam = 1;
+ logpi = 1.14472988584940017414;
+ ls2pi = 0.91893853320467274178;
+ if( x<-34.0 )
+ {
+ q = -x;
+ w = lngamma(q, tmp);
+ p = ap::ifloor(q);
+ i = ap::round(p);
+ if( i%2==0 )
+ {
+ sgngam = -1;
+ }
+ else
+ {
+ sgngam = 1;
+ }
+ z = q-p;
+ if( z>0.5 )
+ {
+ p = p+1;
+ z = p-q;
+ }
+ z = q*sin(ap::pi()*z);
+ result = logpi-log(z)-w;
+ return result;
+ }
+ if( x<13 )
+ {
+ z = 1;
+ p = 0;
+ u = x;
+ while(u>=3)
+ {
+ p = p-1;
+ u = x+p;
+ z = z*u;
+ }
+ while(u<2)
+ {
+ z = z/u;
+ p = p+1;
+ u = x+p;
+ }
+ if( z<0 )
+ {
+ sgngam = -1;
+ z = -z;
+ }
+ else
+ {
+ sgngam = 1;
+ }
+ if( u==2 )
+ {
+ result = log(z);
+ return result;
+ }
+ p = p-2;
+ x = x+p;
+ b = -1378.25152569120859100;
+ b = -38801.6315134637840924+x*b;
+ b = -331612.992738871184744+x*b;
+ b = -1162370.97492762307383+x*b;
+ b = -1721737.00820839662146+x*b;
+ b = -853555.664245765465627+x*b;
+ c = 1;
+ c = -351.815701436523470549+x*c;
+ c = -17064.2106651881159223+x*c;
+ c = -220528.590553854454839+x*c;
+ c = -1139334.44367982507207+x*c;
+ c = -2532523.07177582951285+x*c;
+ c = -2018891.41433532773231+x*c;
+ p = x*b/c;
+ result = log(z)+p;
+ return result;
+ }
+ q = (x-0.5)*log(x)-x+ls2pi;
+ if( x>100000000 )
+ {
+ result = q;
+ return result;
+ }
+ p = 1/(x*x);
+ if( x>=1000.0 )
+ {
+ q = q+((7.9365079365079365079365*0.0001*p-2.7777777777777777777778*0.001)*p+0.0833333333333333333333)/x;
+ }
+ else
+ {
+ a = 8.11614167470508450300*0.0001;
+ a = -5.95061904284301438324*0.0001+p*a;
+ a = 7.93650340457716943945*0.0001+p*a;
+ a = -2.77777777730099687205*0.001+p*a;
+ a = 8.33333333333331927722*0.01+p*a;
+ q = q+a/x;
+ }
+ result = q;
+ return result;
+}
+
+
+double gammastirf(double x)
+{
+ double result;
+ double y;
+ double w;
+ double v;
+ double stir;
+
+ w = 1/x;
+ stir = 7.87311395793093628397E-4;
+ stir = -2.29549961613378126380E-4+w*stir;
+ stir = -2.68132617805781232825E-3+w*stir;
+ stir = 3.47222221605458667310E-3+w*stir;
+ stir = 8.33333333333482257126E-2+w*stir;
+ w = 1+w*stir;
+ y = exp(x);
+ if( x>143.01608 )
+ {
+ v = pow(x, 0.5*x-0.25);
+ y = v*(v/y);
+ }
+ else
+ {
+ y = pow(x, x-0.5)/y;
+ }
+ result = 2.50662827463100050242*y*w;
+ return result;
+}
+
+
+
diff --git a/gammaf.h b/gammaf.h
new file mode 100755
index 0000000..5abb786
--- /dev/null
+++ b/gammaf.h
@@ -0,0 +1,103 @@
+/*************************************************************************
+Cephes Math Library Release 2.8: June, 2000
+Copyright by Stephen L. Moshier
+
+Contributors:
+ * Sergey Bochkanov (ALGLIB project). Translation from C to
+ pseudocode.
+
+See subroutines comments for additional copyrights.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+- Redistributions of source code must retain the above copyright
+ notice, this list of conditions and the following disclaimer.
+
+- Redistributions in binary form must reproduce the above copyright
+ notice, this list of conditions and the following disclaimer listed
+ in this license in the documentation and/or other materials
+ provided with the distribution.
+
+- Neither the name of the copyright holders nor the names of its
+ contributors may 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 _gammaf_h
+#define _gammaf_h
+
+#include "ap.h"
+
+/*************************************************************************
+Gamma function
+
+Input parameters:
+ X - argument
+
+Domain:
+ 0 < X < 171.6
+ -170 < X < 0, X is not an integer.
+
+Relative error:
+ arithmetic domain # trials peak rms
+ IEEE -170,-33 20000 2.3e-15 3.3e-16
+ IEEE -33, 33 20000 9.4e-16 2.2e-16
+ IEEE 33, 171.6 20000 2.3e-15 3.2e-16
+
+Cephes Math Library Release 2.8: June, 2000
+Original copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier
+Translated to AlgoPascal by Bochkanov Sergey (2005, 2006, 2007).
+*************************************************************************/
+double gamma(double x);
+
+
+/*************************************************************************
+Natural logarithm of gamma function
+
+Input parameters:
+ X - argument
+
+Result:
+ logarithm of the absolute value of the Gamma(X).
+
+Output parameters:
+ SgnGam - sign(Gamma(X))
+
+Domain:
+ 0 < X < 2.55e305
+ -2.55e305 < X < 0, X is not an integer.
+
+ACCURACY:
+arithmetic domain # trials peak rms
+ IEEE 0, 3 28000 5.4e-16 1.1e-16
+ IEEE 2.718, 2.556e305 40000 3.5e-16 8.3e-17
+The error criterion was relative when the function magnitude
+was greater than one but absolute when it was less than one.
+
+The following test used the relative error criterion, though
+at certain points the relative error could be much higher than
+indicated.
+ IEEE -200, -4 10000 4.8e-16 1.3e-16
+
+Cephes Math Library Release 2.8: June, 2000
+Copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier
+Translated to AlgoPascal by Bochkanov Sergey (2005, 2006, 2007).
+*************************************************************************/
+double lngamma(double x, double& sgngam);
+
+
+#endif
diff --git a/global.cpp b/global.cpp
new file mode 100755
index 0000000..d725a9d
--- /dev/null
+++ b/global.cpp
@@ -0,0 +1,75 @@
+/*************************************************************************
+GWAMA software: May, 2010
+
+Contributors:
+ * Andrew P Morris amorris at well.ox.ac.uk
+ * Reedik Magi reedik at well.ox.ac.uk
+
+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 "global.h"
+#include <iostream>
+#include <map>
+#include <vector>
+#include <stdlib.h>
+#include <fstream>
+using namespace std;
+global::global()
+{
+ _errNR = 1; //count of errors and warnings found
+ _thresholdPValDir = 1; //p-value threshold for directions
+ _markerCount = 1;
+ _studyCount = 0;
+ _genomicControl = false; //genomic control for input files
+ _genomicControlOutput = false; //genomic control for output
+ _binaryTrait = true; //by default binary trait
+ _randomEffect = false; //use random effect
+ _genderHet = false; //use gender-specific analysis
+ _noAlleles = false; //dont use EA and NEA columns
+ _indel = false; //dont use long allele names
+ _fileList = "gwama.in"; //gwama.in file
+ _markerMap = "N"; //map file name
+ _outputRoot = "gwama"; //output file root (w/o file extension)
+ _version = "2.1";
+
+ outputLambda=1;
+ outputLambda_male=1;
+ outputLambda_female=1;
+
+
+
+ _alternativeMarkerNames.push_back("SNP");
+ _alternativeMarkerNames.push_back("RS");
+ _alternativeMarkerNames.push_back("SNPID");
+ _alternativeMarkerNames.push_back("ID");
+ _alternativeMarkerNames.push_back("MARKER");
+ _alternativeMarkerNames.push_back("MARKERNAME");
+ _alternativeStrands.push_back("STRAND");
+ _alternativeEffectAlleles.push_back("EA");
+ _alternativeEffectAlleles.push_back("EFFECT_ALLELE");
+ _alternativeOtherAlleles.push_back("NON_EFFECT_ALLELE");
+ _alternativeOtherAlleles.push_back("OTHER_ALLELE");
+ _alternativeOtherAlleles.push_back("NEA");
+ _alternativeEffectAlleleFreqs.push_back("EFFECT_ALLELE_FREQUENCY");
+ _alternativeEffectAlleleFreqs.push_back("EAF");
+ _alternativeEffectAlleleFreqs.push_back("EFFECT_ALLELE_FREQ");
+ _alternativeBetas.push_back("BETA");
+ _alternativeBetas.push_back("EFFECT");
+ _alternativeSEs.push_back("SE");
+ _alternativeSEs.push_back("STDERR");
+ _alternativeORs.push_back("OR");
+ _alternativeOR_95Ls.push_back("OR_95L");
+ _alternativeOR_95Us.push_back("OR_95U");
+ _alternativeNs.push_back("N");
+ _alternativeImputeds.push_back("IMPUTED");
+}
diff --git a/global.h b/global.h
new file mode 100755
index 0000000..049ef6b
--- /dev/null
+++ b/global.h
@@ -0,0 +1,75 @@
+/*************************************************************************
+GWAMA software: May, 2010
+
+Contributors:
+ * Andrew P Morris amorris at well.ox.ac.uk
+ * Reedik Magi reedik at well.ox.ac.uk
+
+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 _GLOBAL_H_
+#define _GLOBAL_H_
+
+#include <iostream>
+#include <map>
+#include <vector>
+#include <stdlib.h>
+#include <fstream>
+
+using namespace std;
+class global
+{
+public:
+ int _errNR; //count of errors and warnings found
+ double _thresholdPValDir; //p-value threshold for directions
+ bool _genomicControl; //genomic control for input files
+ bool _genomicControlOutput; //genomic control for output
+ bool _binaryTrait; //by default binary trait
+ bool _randomEffect; //use random effect
+ bool _genderHet; //use gender-specific analysis
+ bool _noAlleles; //use gender-specific analysis
+ bool _indel; //use indel alleles
+ vector <string> _alternativeMarkerNames;
+ vector <string> _alternativeEffectAlleles;
+ vector <string> _alternativeOtherAlleles;
+ vector <string> _alternativeEffectAlleleFreqs;
+ vector <string> _alternativeStrands;
+ vector <string> _alternativeBetas;
+ vector <string> _alternativeSEs;
+ vector <string> _alternativeORs;
+ vector <string> _alternativeOR_95Ls;
+ vector <string> _alternativeOR_95Us;
+ vector <string> _alternativeNs;
+ vector <string> _alternativeImputeds;
+ string _fileList; //gwama.in file
+ string _markerMap; //map file name
+ string _outputRoot; //output file root (w/o file extension)
+ string _version;
+
+ int _markerCount;
+ int _studyCount;
+ vector <string> studyList;
+ vector <string> studyGenders;
+
+ double outputLambda;
+ double outputLambda_male;;
+ double outputLambda_female;
+
+
+ global(); //constructor
+};
+
+
+
+#endif
+
diff --git a/igammaf.cpp b/igammaf.cpp
new file mode 100755
index 0000000..8baa331
--- /dev/null
+++ b/igammaf.cpp
@@ -0,0 +1,430 @@
+/*************************************************************************
+Cephes Math Library Release 2.8: June, 2000
+Copyright by Stephen L. Moshier
+
+Contributors:
+ * Sergey Bochkanov (ALGLIB project). Translation from C to
+ pseudocode.
+
+See subroutines comments for additional copyrights.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+- Redistributions of source code must retain the above copyright
+ notice, this list of conditions and the following disclaimer.
+
+- Redistributions in binary form must reproduce the above copyright
+ notice, this list of conditions and the following disclaimer listed
+ in this license in the documentation and/or other materials
+ provided with the distribution.
+
+- Neither the name of the copyright holders nor the names of its
+ contributors may 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 <stdafx.h>
+#include "igammaf.h"
+
+/*************************************************************************
+Incomplete gamma integral
+
+The function is defined by
+
+ x
+ -
+ 1 | | -t a-1
+ igam(a,x) = ----- | e t dt.
+ - | |
+ | (a) -
+ 0
+
+
+In this implementation both arguments must be positive.
+The integral is evaluated by either a power series or
+continued fraction expansion, depending on the relative
+values of a and x.
+
+ACCURACY:
+
+ Relative error:
+arithmetic domain # trials peak rms
+ IEEE 0,30 200000 3.6e-14 2.9e-15
+ IEEE 0,100 300000 9.9e-14 1.5e-14
+
+Cephes Math Library Release 2.8: June, 2000
+Copyright 1985, 1987, 2000 by Stephen L. Moshier
+*************************************************************************/
+double incompletegamma(double a, double x)
+{
+ double result;
+ double igammaepsilon;
+ double ans;
+ double ax;
+ double c;
+ double r;
+ double tmp;
+
+ igammaepsilon = 0.000000000000001;
+ if( x<=0||a<=0 )
+ {
+ result = 0;
+ return result;
+ }
+ if( x>1&&x>a )
+ {
+ result = 1-incompletegammac(a, x);
+ return result;
+ }
+ ax = a*log(x)-x-lngamma(a, tmp);
+ if( ax<-709.78271289338399 )
+ {
+ result = 0;
+ return result;
+ }
+ ax = exp(ax);
+ r = a;
+ c = 1;
+ ans = 1;
+ do
+ {
+ r = r+1;
+ c = c*x/r;
+ ans = ans+c;
+ }
+ while(c/ans>igammaepsilon);
+ result = ans*ax/a;
+ return result;
+}
+
+
+/*************************************************************************
+Complemented incomplete gamma integral
+
+The function is defined by
+
+
+ igamc(a,x) = 1 - igam(a,x)
+
+ inf.
+ -
+ 1 | | -t a-1
+ = ----- | e t dt.
+ - | |
+ | (a) -
+ x
+
+
+In this implementation both arguments must be positive.
+The integral is evaluated by either a power series or
+continued fraction expansion, depending on the relative
+values of a and x.
+
+ACCURACY:
+
+Tested at random a, x.
+ a x Relative error:
+arithmetic domain domain # trials peak rms
+ IEEE 0.5,100 0,100 200000 1.9e-14 1.7e-15
+ IEEE 0.01,0.5 0,100 200000 1.4e-13 1.6e-15
+
+Cephes Math Library Release 2.8: June, 2000
+Copyright 1985, 1987, 2000 by Stephen L. Moshier
+*************************************************************************/
+double incompletegammac(double a, double x)
+{
+ double result;
+ double igammaepsilon;
+ double igammabignumber;
+ double igammabignumberinv;
+ double ans;
+ double ax;
+ double c;
+ double yc;
+ double r;
+ double t;
+ double y;
+ double z;
+ double pk;
+ double pkm1;
+ double pkm2;
+ double qk;
+ double qkm1;
+ double qkm2;
+ double tmp;
+
+ igammaepsilon = 0.000000000000001;
+ igammabignumber = 4503599627370496.0;
+ igammabignumberinv = 2.22044604925031308085*0.0000000000000001;
+ if( x<=0||a<=0 )
+ {
+ result = 1;
+ return result;
+ }
+ if( x<1||x<a )
+ {
+ result = 1-incompletegamma(a, x);
+ return result;
+ }
+ ax = a*log(x)-x-lngamma(a, tmp);
+ if( ax<-709.78271289338399 )
+ {
+ result = 0;
+ return result;
+ }
+ ax = exp(ax);
+ y = 1-a;
+ z = x+y+1;
+ c = 0;
+ pkm2 = 1;
+ qkm2 = x;
+ pkm1 = x+1;
+ qkm1 = z*x;
+ ans = pkm1/qkm1;
+ do
+ {
+ c = c+1;
+ y = y+1;
+ z = z+2;
+ yc = y*c;
+ pk = pkm1*z-pkm2*yc;
+ qk = qkm1*z-qkm2*yc;
+ if( qk!=0 )
+ {
+ r = pk/qk;
+ t = fabs((ans-r)/r);
+ ans = r;
+ }
+ else
+ {
+ t = 1;
+ }
+ pkm2 = pkm1;
+ pkm1 = pk;
+ qkm2 = qkm1;
+ qkm1 = qk;
+ if( fabs(pk)>igammabignumber )
+ {
+ pkm2 = pkm2*igammabignumberinv;
+ pkm1 = pkm1*igammabignumberinv;
+ qkm2 = qkm2*igammabignumberinv;
+ qkm1 = qkm1*igammabignumberinv;
+ }
+ }
+ while(t>igammaepsilon);
+ result = ans*ax;
+ return result;
+}
+
+
+/*************************************************************************
+Inverse of complemented imcomplete gamma integral
+
+Given p, the function finds x such that
+
+ igamc( a, x ) = p.
+
+Starting with the approximate value
+
+ 3
+ x = a t
+
+ where
+
+ t = 1 - d - ndtri(p) sqrt(d)
+
+and
+
+ d = 1/9a,
+
+the routine performs up to 10 Newton iterations to find the
+root of igamc(a,x) - p = 0.
+
+ACCURACY:
+
+Tested at random a, p in the intervals indicated.
+
+ a p Relative error:
+arithmetic domain domain # trials peak rms
+ IEEE 0.5,100 0,0.5 100000 1.0e-14 1.7e-15
+ IEEE 0.01,0.5 0,0.5 100000 9.0e-14 3.4e-15
+ IEEE 0.5,10000 0,0.5 20000 2.3e-13 3.8e-14
+
+Cephes Math Library Release 2.8: June, 2000
+Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier
+*************************************************************************/
+double invincompletegammac(double a, double y0)
+{
+ double result;
+ double igammaepsilon;
+ double iinvgammabignumber;
+ double x0;
+ double x1;
+ double x;
+ double yl;
+ double yh;
+ double y;
+ double d;
+ double lgm;
+ double dithresh;
+ int i;
+ int dir;
+ double tmp;
+
+ igammaepsilon = 0.000000000000001;
+ iinvgammabignumber = 4503599627370496.0;
+ x0 = iinvgammabignumber;
+ yl = 0;
+ x1 = 0;
+ yh = 1;
+ dithresh = 5*igammaepsilon;
+ d = 1/(9*a);
+ y = 1-d-invnormaldistribution(y0)*sqrt(d);
+ x = a*y*y*y;
+ lgm = lngamma(a, tmp);
+ i = 0;
+ while(i<10)
+ {
+ if( x>x0||x<x1 )
+ {
+ d = 0.0625;
+ break;
+ }
+ y = incompletegammac(a, x);
+ if( y<yl||y>yh )
+ {
+ d = 0.0625;
+ break;
+ }
+ if( y<y0 )
+ {
+ x0 = x;
+ yl = y;
+ }
+ else
+ {
+ x1 = x;
+ yh = y;
+ }
+ d = (a-1)*log(x)-x-lgm;
+ if( d<-709.78271289338399 )
+ {
+ d = 0.0625;
+ break;
+ }
+ d = -exp(d);
+ d = (y-y0)/d;
+ if( fabs(d/x)<igammaepsilon )
+ {
+ result = x;
+ return result;
+ }
+ x = x-d;
+ i = i+1;
+ }
+ if( x0==iinvgammabignumber )
+ {
+ if( x<=0 )
+ {
+ x = 1;
+ }
+ while(x0==iinvgammabignumber)
+ {
+ x = (1+d)*x;
+ y = incompletegammac(a, x);
+ if( y<y0 )
+ {
+ x0 = x;
+ yl = y;
+ break;
+ }
+ d = d+d;
+ }
+ }
+ d = 0.5;
+ dir = 0;
+ i = 0;
+ while(i<400)
+ {
+ x = x1+d*(x0-x1);
+ y = incompletegammac(a, x);
+ lgm = (x0-x1)/(x1+x0);
+ if( fabs(lgm)<dithresh )
+ {
+ break;
+ }
+ lgm = (y-y0)/y0;
+ if( fabs(lgm)<dithresh )
+ {
+ break;
+ }
+ if( x<=0.0 )
+ {
+ break;
+ }
+ if( y>=y0 )
+ {
+ x1 = x;
+ yh = y;
+ if( dir<0 )
+ {
+ dir = 0;
+ d = 0.5;
+ }
+ else
+ {
+ if( dir>1 )
+ {
+ d = 0.5*d+0.5;
+ }
+ else
+ {
+ d = (y0-yl)/(yh-yl);
+ }
+ }
+ dir = dir+1;
+ }
+ else
+ {
+ x0 = x;
+ yl = y;
+ if( dir>0 )
+ {
+ dir = 0;
+ d = 0.5;
+ }
+ else
+ {
+ if( dir<-1 )
+ {
+ d = 0.5*d;
+ }
+ else
+ {
+ d = (y0-yl)/(yh-yl);
+ }
+ }
+ dir = dir-1;
+ }
+ i = i+1;
+ }
+ result = x;
+ return result;
+}
+
+
+
diff --git a/igammaf.h b/igammaf.h
new file mode 100755
index 0000000..e86cf9d
--- /dev/null
+++ b/igammaf.h
@@ -0,0 +1,156 @@
+/*************************************************************************
+Cephes Math Library Release 2.8: June, 2000
+Copyright by Stephen L. Moshier
+
+Contributors:
+ * Sergey Bochkanov (ALGLIB project). Translation from C to
+ pseudocode.
+
+See subroutines comments for additional copyrights.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+- Redistributions of source code must retain the above copyright
+ notice, this list of conditions and the following disclaimer.
+
+- Redistributions in binary form must reproduce the above copyright
+ notice, this list of conditions and the following disclaimer listed
+ in this license in the documentation and/or other materials
+ provided with the distribution.
+
+- Neither the name of the copyright holders nor the names of its
+ contributors may 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 _igammaf_h
+#define _igammaf_h
+
+#include "ap.h"
+
+#include "gammaf.h"
+#include "normaldistr.h"
+
+
+/*************************************************************************
+Incomplete gamma integral
+
+The function is defined by
+
+ x
+ -
+ 1 | | -t a-1
+ igam(a,x) = ----- | e t dt.
+ - | |
+ | (a) -
+ 0
+
+
+In this implementation both arguments must be positive.
+The integral is evaluated by either a power series or
+continued fraction expansion, depending on the relative
+values of a and x.
+
+ACCURACY:
+
+ Relative error:
+arithmetic domain # trials peak rms
+ IEEE 0,30 200000 3.6e-14 2.9e-15
+ IEEE 0,100 300000 9.9e-14 1.5e-14
+
+Cephes Math Library Release 2.8: June, 2000
+Copyright 1985, 1987, 2000 by Stephen L. Moshier
+*************************************************************************/
+double incompletegamma(double a, double x);
+
+
+/*************************************************************************
+Complemented incomplete gamma integral
+
+The function is defined by
+
+
+ igamc(a,x) = 1 - igam(a,x)
+
+ inf.
+ -
+ 1 | | -t a-1
+ = ----- | e t dt.
+ - | |
+ | (a) -
+ x
+
+
+In this implementation both arguments must be positive.
+The integral is evaluated by either a power series or
+continued fraction expansion, depending on the relative
+values of a and x.
+
+ACCURACY:
+
+Tested at random a, x.
+ a x Relative error:
+arithmetic domain domain # trials peak rms
+ IEEE 0.5,100 0,100 200000 1.9e-14 1.7e-15
+ IEEE 0.01,0.5 0,100 200000 1.4e-13 1.6e-15
+
+Cephes Math Library Release 2.8: June, 2000
+Copyright 1985, 1987, 2000 by Stephen L. Moshier
+*************************************************************************/
+double incompletegammac(double a, double x);
+
+
+/*************************************************************************
+Inverse of complemented imcomplete gamma integral
+
+Given p, the function finds x such that
+
+ igamc( a, x ) = p.
+
+Starting with the approximate value
+
+ 3
+ x = a t
+
+ where
+
+ t = 1 - d - ndtri(p) sqrt(d)
+
+and
+
+ d = 1/9a,
+
+the routine performs up to 10 Newton iterations to find the
+root of igamc(a,x) - p = 0.
+
+ACCURACY:
+
+Tested at random a, p in the intervals indicated.
+
+ a p Relative error:
+arithmetic domain domain # trials peak rms
+ IEEE 0.5,100 0,0.5 100000 1.0e-14 1.7e-15
+ IEEE 0.01,0.5 0,0.5 100000 9.0e-14 3.4e-15
+ IEEE 0.5,10000 0,0.5 20000 2.3e-13 3.8e-14
+
+Cephes Math Library Release 2.8: June, 2000
+Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier
+*************************************************************************/
+double invincompletegammac(double a, double y0);
+
+
+#endif
diff --git a/main.cpp b/main.cpp
new file mode 100755
index 0000000..61c32cf
--- /dev/null
+++ b/main.cpp
@@ -0,0 +1,181 @@
+/*************************************************************************
+GWAMA software: May, 2009
+
+Contributors:
+ * Andrew P Morris amorris at well.ox.ac.uk
+ * Reedik Magi reedik at well.ox.ac.uk
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*************************************************************************/
+
+
+#include <iostream>
+#include <map>
+#include <vector>
+#include <stdlib.h>
+#include <fstream>
+#include "marker.h"
+#include "study.h"
+#include <math.h>
+#include <cctype> // std::toupper
+#include <zlib.h>
+#include "global.h"
+#include "commandLine.h"
+#include "readFile.h"
+#include "tools.h"
+
+using namespace std;
+
+
+
+
+
+int
+main(int argc, char *argv[])
+{
+ string gwamaversion = "2.1";
+ global _g;
+
+ char sb [11];
+ vector <string> fileGenders;
+ vector <string> snpNames;
+ vector <study> studies;
+ bool proc;
+
+ proc = readCommandline( argc, argv, _g);
+ if (!proc){cerr << "Error reading command line options! Exit program." << endl; exit(1);}
+ string gwamagc = _g._outputRoot + ".gc.out";
+ string gwamaout = _g._outputRoot + + ".out";
+ string gwamaerr = _g._outputRoot + + ".err.out";
+ string gwamalog = _g._outputRoot + + ".log.out";
+ ofstream LAMBDA_FILE (gwamagc.c_str());
+ ofstream ERR (gwamaerr.c_str());
+ ofstream LOG (gwamalog.c_str());
+
+ LAMBDA_FILE << "Cohort\tdirectly_genotyped_markers_lambda\t\tdirectly_genotyped_markers_count\timputed_markers_lambda\timputed_markers_count\n";
+ LOG << "Running GWAMA " << gwamaversion << endl;
+ proc = readFilelist(_g, ERR, LOG);
+ if (!proc){cerr << "Error reading file list! Exit program." << endl; exit(1);}
+
+ map <string, int> markerPosition;
+ vector <marker> markerList;
+ for (int i = 0; i < _g._studyCount; i++)
+ {
+ cout << "------------------------------------\nReading file: " << _g.studyList[i] << endl;
+ if (_g._genderHet && _g.studyGenders.size()>i) cout << "Gender status: " << _g.studyGenders[i] << endl;
+ else if (_g._genderHet){ cerr << "Error reading gender status from filelist. Exit program<< "; exit(1);}
+ proc = readCohort(i,_g,markerPosition,markerList,ERR,LOG, LAMBDA_FILE);
+ }
+ ofstream OUT (gwamaout.c_str());
+
+ if (_g._markerCount==1) {cerr << "No markers found! Exit program!"<<endl; exit(1);}
+
+ if (_g._markerMap != "N")
+ {
+ cout << "------------------------------------\nReading map file"<<endl;
+ proc = readMapFile(_g._markerMap,markerList,markerPosition, ERR, _g._errNR, LOG);
+ }
+ cout << "------------------------------------\nPreparing output..."<<endl;
+ if (_g._genomicControlOutput)
+ {
+ vector <double> _chiStat;
+ vector <double> _chiStat_male;
+ vector <double> _chiStat_female;
+ for (int i = 0; i < _g._markerCount-1; i++)
+ {
+ if (markerList[i]._studyCount>0)
+ {_chiStat.push_back((markerList[i].sumBeta(0)*markerList[i].sumBeta(0))/(markerList[i].sumSE(0)*markerList[i].sumSE(0)));}
+ if (_g._genderHet)
+ {
+ if (markerList[i].male_studyCount>0)
+ {_chiStat_male.push_back((markerList[i].male_sumBeta(0)*markerList[i].male_sumBeta(0))/(markerList[i].male_sumSE(0)*markerList[i].male_sumSE(0)));}
+ if (markerList[i].female_studyCount>0)
+ {_chiStat_female.push_back((markerList[i].female_sumBeta(0)*markerList[i].female_sumBeta(0))/(markerList[i].female_sumSE(0)*markerList[i].female_sumSE(0)));}
+
+ }
+ }
+
+
+ double _median = 0;
+ sortVec( _chiStat, _chiStat.size());
+ if ((_chiStat.size())%2!=0) _median = _chiStat[((_chiStat.size())/2)-1];
+ else _median = (_chiStat[(_chiStat.size())/2 -1] + _chiStat[(_chiStat.size())/2])/2;
+ _g.outputLambda = _median/0.4549364; //median of chi-sq from R ... qchisq(0.5, df= 1)
+ if (_g.outputLambda>1)LOG << "Output p-value is corrected with lambda=" << _g.outputLambda << endl;
+ else LOG << "Output p-value is not corrected. Lambda=" << _g.outputLambda << endl;
+ if (_g.outputLambda>1)cout <<"------------------------------------\nOutput p-value is corrected with lambda=" << _g.outputLambda << endl;
+ else cout <<"------------------------------------\nOutput p-value is not corrected. Lambda=" << _g.outputLambda << endl;
+ if (_g._genderHet)
+ {
+ _median = 0;
+ sortVec( _chiStat_male, _chiStat_male.size());
+ if ((_chiStat_male.size())%2!=0) _median = _chiStat_male[((_chiStat_male.size())/2)-1];
+ else _median = (_chiStat_male[(_chiStat_male.size()-1)/2 -1] + _chiStat_male[(_chiStat_male.size()-1)/2])/2;
+ _g.outputLambda_male = _median/0.4549364; //median of chi-sq from R ... qchisq(0.5, df= 1)
+ if (_g.outputLambda_male>1)LOG << "Output p-value is corrected with lambda=" << _g.outputLambda_male << endl;
+ else LOG << "Output p-value is not corrected. Lambda=" << _g.outputLambda_male << endl;
+ if (_g.outputLambda_male>1)cout <<"------------------------------------\nOutput p-value is corrected with lambda=" << _g.outputLambda_male << endl;
+ else cout <<"------------------------------------\nOutput p-value is not corrected. Lambda=" << _g.outputLambda_male << endl;
+
+ _median = 0;
+ sortVec( _chiStat_female, _chiStat_female.size());
+ if ((_chiStat_female.size())%2!=0) _median = _chiStat_female[((_chiStat_female.size())/2)-1];
+ else _median = (_chiStat_female[(_chiStat_female.size()-1)/2 -1] + _chiStat_female[(_chiStat_female.size()-1)/2])/2;
+ _g.outputLambda_female = _median/0.4549364; //median of chi-sq from R ... qchisq(0.5, df= 1)
+ if (_g.outputLambda_female>1)LOG << "Output p-value is corrected with lambda=" << _g.outputLambda_female << endl;
+ else LOG << "Output p-value is not corrected. Lambda=" << _g.outputLambda_female << endl;
+ if (_g.outputLambda_female>1)cout <<"------------------------------------\nOutput p-value is corrected with lambda=" << _g.outputLambda_female << endl;
+ else cout <<"------------------------------------\nOutput p-value is not corrected. Lambda=" << _g.outputLambda_female << endl;
+
+
+ }
+
+ }
+ if (_g._markerMap != "N")
+ {OUT <<"chromosome\tposition\t";}
+ if (_g._binaryTrait)
+ {OUT << "rs_number\treference_allele\tother_allele\teaf\tOR\tOR_se\tOR_95L\tOR_95U\tz\tp-value\t_-log10_p-value\tq_statistic\tq_p-value\ti2\tn_studies\tn_samples\teffects";}
+ else
+ {OUT << "rs_number\treference_allele\tother_allele\teaf\tbeta\tse\tbeta_95L\tbeta_95U\tz\tp-value\t_-log10_p-value\tq_statistic\tq_p-value\ti2\tn_studies\tn_samples\teffects";}
+
+ if (_g._genderHet)
+ {
+ if (_g._binaryTrait)
+ {OUT << "\tmale_eaf\tmale_OR\tmale_OR_se\tmale_OR_95L\tmale_OR_95U\tmale_z\tmale_p-value\tmale_n_studies\tmale_n_samples\tfemale_eaf\tfemale_OR\tfemale_OR_se\tfemale_OR_95L\tfemale_OR_95U\tfemale_z\tfemale_p-value\tfemale_n_studies\tfemale_n_samples\tgender_differentiated_p-value\tgender_heterogeneity_p-value";}
+ else
+ {OUT << "\tmale_eaf\tmale_beta\tmale_se\tmale_beta_95L\tmale_beta_95U\tmale_z\tmale_p-value\tmale_n_studies\tmale_n_samples\tfemale_eaf\tfemale_beta\tfemale_se\tfemale_beta_95L\tfemale_beta_95U\tfemale_z\tfemale_p-value\tfemale_n_studies\tfemale_n_samples\tgender_differentiated_p-value\tgender_heterogeneity_p-value";}
+ }
+ OUT << endl;
+
+
+ for (int i = 0; i < _g._markerCount-1; i++)
+ {
+ if (_g._randomEffect)markerList[i].doRandom(_g);
+ OUT << markerList[i].printMarker( _g) << endl;
+ }
+
+ OUT.close();
+ LOG << "Analysis finished." << endl;
+
+ ERR.close();
+ LOG.close();
+ cout << "------------------------------------\nGWAMA program finished current job successfully!" << endl;
+ cout << "Please check " << _g._outputRoot << ".log.out for full information." << endl;
+ cout << "Analysis finished." << endl;
+ return 0;
+}
+
+
+
+
+
diff --git a/marker.cpp b/marker.cpp
new file mode 100755
index 0000000..c753855
--- /dev/null
+++ b/marker.cpp
@@ -0,0 +1,1078 @@
+/*************************************************************************
+GWAMA software: May, 2009
+
+Contributors:
+ * Andrew P Morris amorris at well.ox.ac.uk
+ * Reedik Magi reedik at well.ox.ac.uk
+
+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 "marker.h"
+#include <iostream>
+#include <map>
+#include <vector>
+#include <stdlib.h>
+#include <fstream>
+#include <math.h>
+#include "statistics.h"
+#include "study.h"
+#include "chisquaredistr.h"
+#include "global.h"
+#include "tools.h"
+
+using namespace std;
+
+
+double
+abs_d(double x)
+{
+if (x < 0) return x * -1;
+return x;
+
+}
+
+marker::marker(std::string name, int popCount)
+{
+ _name = name;
+ for (int i = 0; i < popCount; i++)
+ {
+ _direction += "?";
+ }
+ _effectAllele="N";
+ _otherAllele="N";
+ _studyCount = 0;
+ male_studyCount = 0;
+ female_studyCount = 0;
+ _position = -9;
+ _chromosome = -9;
+ _N=0;
+ _Nok=true;
+ _EAF=0;
+ _B=0;
+ _W=0;
+ _W2=0;
+ _Q=0;
+
+ _Nmale=0;
+ _EAFmale=0;
+ _Bmale=0;
+ _Wmale=0;
+ _W2male=0;
+ _Qmale=0;
+ _Brmale=0;
+ _Wrmale=0;
+ _Qrmale=0;
+
+ _Nfemale=0;
+ _EAFfemale=0;
+ _Bfemale=0;
+ _Wfemale=0;
+ _W2female=0;
+ _Qfemale=0;
+ _Brfemale=0;
+ _Wrfemale=0;
+ _Qrfemale=0;
+
+}
+
+marker::~marker(void)
+{
+}
+
+bool
+marker::addMap(int chromosome, int position)
+{
+ _chromosome=chromosome;
+ _position=position;
+ return true;
+}
+
+
+int
+marker::addCohort(int fileNr, global & g, double directLambda, double imputedLambda,
+ string myMarker,string myEffectAllele, string myOtherAllele, double myEffectAlleleFreq,
+ bool myStrand, double myBeta, double mySE, bool myImputed, int myN, ofstream & ERR,
+ ofstream & LOG, problem & markerProblems)
+{
+ //READING DATA
+ char sb [11];
+ double _se = 0;
+ if (g._genomicControl)
+ {
+ if (imputedLambda<1)imputedLambda = 1;
+ if (directLambda<1)directLambda = 1;
+ if (myImputed) _se = mySE * sqrt(imputedLambda);
+ else _se = mySE * sqrt(directLambda);
+ }
+ else _se = mySE;
+
+ if (_effectAllele.compare("N")==0)_effectAllele = myEffectAllele;
+ if (_otherAllele.compare("N")==0)_otherAllele = myOtherAllele;
+
+ double _beta = 0;
+ double _eaf = 0;
+ if (_effectAllele == myEffectAllele && _otherAllele == myOtherAllele)
+ {
+ _beta = myBeta;
+ if (myEffectAlleleFreq>0) _eaf = myEffectAlleleFreq;
+ }
+ else if (_effectAllele == myOtherAllele && _otherAllele == myEffectAllele)//we need to align according to other allele
+ {
+ _beta = myBeta * -1;
+ if (myEffectAlleleFreq>0) _eaf = 1 - myEffectAlleleFreq;
+ }
+ else //something is wrong with alleles - i'll switch strand and check if that helps
+ {
+ if (_effectAllele == flip(myEffectAllele) && _otherAllele == flip(myOtherAllele))
+ {
+ _beta = myBeta;
+ if (myEffectAlleleFreq>0) _eaf = myEffectAlleleFreq;
+ sprintf(sb, "W%.9d" , g._errNR); g._errNR++;
+ LOG << "Warning: Marker " << myMarker << " has wrong strand. (" << sb << ")" <<endl;
+ ERR << sb << " " << g.studyList[fileNr] << " warning: Marker " << _name << " has wrong strand. Strand flipped." <<endl;
+ markerProblems.problemStrand++;
+ }
+ else if (_effectAllele == flip(myOtherAllele) && _otherAllele == flip(myEffectAllele))
+ {
+ _beta = myBeta * -1;
+ if (myEffectAlleleFreq>0) _eaf = 1 - myEffectAlleleFreq;
+ sprintf(sb, "W%.9d" , g._errNR); g._errNR++;
+ LOG << "Warning: Marker " << myMarker << " has wrong strand. (" << sb << ")" <<endl;
+ ERR << sb << " " << g.studyList[fileNr] << " warning: Marker " << _name << " has wrong strand. Strand flipped." <<endl;
+ markerProblems.problemStrand++;
+ }
+ else //cannot resolve the problem - marker dropped
+ {
+ sprintf(sb, "E%.9d" , g._errNR); g._errNR++;
+ LOG << "Error: Marker " << myMarker << " has wrong alleles. (" << sb << ")" <<endl;
+ ERR << sb << " " << g.studyList[fileNr] << " error: Marker " << _name << " has wrong alleles. Marker dropped." <<endl;
+ markerProblems.wrongAlleles++;
+ return 1;
+ }
+ }
+ markerProblems.markersOK++;
+ //ADD STATS
+ double _wstat = 1/ (_se*_se);
+ _B+=_beta*_wstat;
+ _W+=_wstat;
+ _W2+=_wstat*_wstat;
+ _Q+=(_beta*_beta)*_wstat;
+
+
+ if (g._randomEffect)
+ {
+ all_beta.push_back(_beta);
+ all_se.push_back(_se);
+ }
+
+ if (g._genderHet)
+ {
+ if (g.studyGenders[fileNr]=="M")
+ {
+ _Bmale+=_beta*_wstat;
+ _Wmale+=_wstat;
+ _W2male+=_wstat*_wstat;
+ _Qmale+=(_beta*_beta)*_wstat;
+ if (g._randomEffect)
+ {
+ male_beta.push_back(_beta);
+ male_se.push_back(_se);
+ }
+
+ }
+ else if
+ (g.studyGenders[fileNr]=="F")
+ {
+ _Bfemale+=_beta*_wstat;
+ _Wfemale+=_wstat;
+ _W2female+=_wstat*_wstat;
+ _Qfemale+=(_beta*_beta)*_wstat;
+ if (g._randomEffect)
+ {
+ female_beta.push_back(_beta);
+ female_se.push_back(_se);
+ }
+
+ }
+ }
+
+
+
+ //ADD N, POPCOUNT & WEIGHTED EAF
+ _studyCount++;
+
+
+ if (abs_d(_EAF - _eaf)>0.3 && _EAF >0 && _eaf >0) //large diversity of allele frequencies - printing warning to log file
+ {
+ sprintf(sb, "W%.9d" , g._errNR); g._errNR++;
+ LOG << "Warning: Marker " << _name << " has large effect allele frequency discrepancy. (" << sb << ")" <<endl;
+
+ ERR << sb << " Warning: Marker " << myMarker << " has large effect allele frequency discrepancy." <<endl;
+ ERR << sb << " Weighted EAF: " << _EAF << ", but file " << g.studyList[fileNr] << " has EAF: " << _eaf <<endl;
+ markerProblems.problemStrand++;
+ }
+
+
+
+
+ if (myN>0 && _eaf>0)
+ {
+ if (_N>=0) _EAF = ((_EAF*_N)+(_eaf*myN))/(_N+myN);
+ if (_N>=0)_N+=myN;
+ if (g._genderHet)
+ {
+
+ if (g.studyGenders[fileNr]=="M")
+ {
+ male_studyCount++;
+ if (_Nmale>=0) _EAFmale = ((_EAFmale*_Nmale)+(_eaf*myN))/(_Nmale+myN);
+ if (_Nmale>=0)_Nmale+=myN;
+ }
+ else if
+ (g.studyGenders[fileNr]=="F")
+ {
+ female_studyCount++;
+ if (_Nfemale>=0) _EAFfemale = ((_EAFfemale*_Nfemale)+(_eaf*myN))/(_Nfemale+myN);
+ if (_Nfemale>=0)_Nfemale+=myN;
+ }
+ }
+ }
+ else if (myN>0)
+ {
+ if (_N>=0) _N+=myN;
+ if (g._genderHet)
+ {
+
+ if (g.studyGenders[fileNr]=="M")
+ {
+ male_studyCount++;
+ if (_Nmale>=0) _Nmale+=myN;
+ }
+ else if
+ (g.studyGenders[fileNr]=="F")
+ {
+ female_studyCount++;
+ if (_Nfemale>=0) _Nfemale+=myN;
+ }
+ }
+ }
+ else //stop using N as one of populations didnt had number
+ {
+ if (_eaf>0)
+ {
+ if (_EAF>0)_EAF = (_EAF + _eaf)/2;
+ else _EAF = _eaf;
+ }
+ if (g._genderHet)
+ {
+
+ if (g.studyGenders[fileNr]=="M")
+ {
+ male_studyCount++;
+ }
+ else if
+ (g.studyGenders[fileNr]=="F")
+ {
+ female_studyCount++;
+ }
+ }
+ _Nok=false;
+ }
+
+ //direction
+ double x = abs_d(_beta * _beta * (1/(_se * _se)));
+ if (bdm_chi2(x, 1)<=g._thresholdPValDir)
+ {
+ if (_beta>0)_direction[fileNr]='+';
+ else if (_beta<0)_direction[fileNr]='-';
+ else _direction[fileNr]='0';
+ }
+
+ return 0;
+}
+
+
+std::string marker::printMarker(global & g)
+{
+ char sb [1024]; //temporary char buffer
+ std::string x=""; //starting to collect all information into this string
+ double sbd;
+ bool useRandom = g._randomEffect;
+ bool map = (g._markerMap!="N");
+ bool useMetaLambda = g._genomicControlOutput;
+ bool _bt = g._binaryTrait;
+ double outputLambda = g.outputLambda;
+ double outputLambda_male = g.outputLambda_male;
+ double outputLambda_female = g.outputLambda_female;
+
+if (map)
+{
+ sprintf(sb, "%d\t" , _chromosome);
+ x.append(sb);
+ sprintf(sb, "%d\t" , _position);
+ x.append(sb);
+}
+
+
+ x.append(_name); // printing SNP name
+
+ x.append("\t"); //printing reference allele
+ x+=_effectAllele;
+
+ x.append("\t"); //printing other allele
+ x+=_otherAllele;
+
+ x.append("\t"); // EAF
+ if (_EAF>0)sprintf(sb, "%.6f" , _EAF);
+ else sprintf(sb, "%d" , -9);
+ x.append(sb);
+
+ if (_bt)
+ {
+ x.append("\t");//printing OR
+ sbd = OR(useRandom);
+ sprintf(sb, "%.6f" , sbd);
+ x.append(sb);
+
+ x.append("\t");//printing se
+ sbd = (OR(useRandom)-Lower95OR(useRandom))/1.96;
+ sprintf(sb, "%.6f" , sbd);
+ x.append(sb);
+
+
+ x.append("\t");//printing L95 OR
+ sbd = Lower95OR(useRandom);
+ sprintf(sb, "%.6f" , sbd);
+ x.append(sb);
+
+ x.append("\t");//printing U95 OR
+ sbd = Upper95OR(useRandom);
+ sprintf(sb, "%.6f" , sbd);
+ x.append(sb);
+ }
+ else
+ {
+ x.append("\t");//printing average beta
+ sbd = sumBeta(useRandom);
+ sprintf(sb, "%.6f" , sbd);
+ x.append(sb);
+
+ x.append("\t");//printing se
+ sbd = sumSE(useRandom);
+ sprintf(sb, "%.6f" , sbd);
+ x.append(sb);
+
+
+ x.append("\t");//printing L95 beta
+ sbd = Lower95Beta(useRandom);
+ sprintf(sb, "%.6f" , sbd);
+ x.append(sb);
+
+ x.append("\t");//printing U95 beta
+ sbd = Upper95Beta(useRandom);
+ sprintf(sb, "%.6f" , sbd);
+ x.append(sb);
+ }
+
+ x.append("\t");//printing z
+ sbd = z(useRandom);
+ sprintf(sb, "%.6f" , sbd);
+ x.append(sb);
+
+ x.append("\t");//p-value
+ if (useMetaLambda) sbd = pValue(outputLambda, useRandom);
+ else sbd = pValue(1, useRandom);
+ if (sbd>=0.001) {sprintf(sb, "%.6f" , sbd);}
+ else {sprintf(sb, "%.2E" , sbd);}
+ x.append(sb);
+
+ x.append("\t");//-10 log p-value
+ sprintf(sb, "%.6f" , -log10(sbd));
+ x.append(sb);
+
+ x.append("\t");//q-statistic
+ sbd = qStatistic();
+ sprintf(sb, "%.6f" , sbd);
+ x.append(sb);
+
+ x.append("\t");//qPValue
+ sbd = qPValue();
+ if (sbd>=0.001) {sprintf(sb, "%.6f" , sbd);}
+ else if (sbd<0.001){sprintf(sb, "%.2E" , sbd);}
+ else {sprintf(sb, "1");}
+ x.append(sb);
+
+ x.append("\t");//i2
+ sbd = i2();
+ sprintf(sb, "%.6f" , sbd);
+ x.append(sb);
+
+ x.append("\t"); // printing study nr
+ sprintf(sb, "%d" , _studyCount);
+ x.append(sb);
+
+ x.append("\t"); // printing sample nr
+ if (_Nok)sprintf(sb, "%d" , _N);
+ else sprintf(sb, "%d" , -9);
+ x.append(sb);
+
+ x.append("\t"); //printing effect directions
+ x.append(_direction);
+
+/////////////
+ if (g._genderHet)
+ {
+
+/////////////MALE
+ if (male_studyCount>0)
+ {
+ x.append("\t"); // EAF
+ if (_EAFmale>0)sprintf(sb, "%.6f" , _EAFmale);
+ else sprintf(sb, "%d" , -9);
+ x.append(sb);
+
+ if (_bt)
+ {
+ x.append("\t");//printing OR
+ sbd = male_OR(useRandom);
+ sprintf(sb, "%.6f" , sbd);
+ x.append(sb);
+
+ x.append("\t");//printing se
+ sbd = (male_OR(useRandom)-male_Lower95OR(useRandom))/1.96;
+ sprintf(sb, "%.6f" , sbd);
+ x.append(sb);
+
+
+ x.append("\t");//printing L95 OR
+ sbd = male_Lower95OR(useRandom);
+ sprintf(sb, "%.6f" , sbd);
+ x.append(sb);
+
+ x.append("\t");//printing U95 OR
+ sbd = male_Upper95OR(useRandom);
+ sprintf(sb, "%.6f" , sbd);
+ x.append(sb);
+ }
+ else
+ {
+ x.append("\t");//printing average beta
+ sbd = male_sumBeta(useRandom);
+ sprintf(sb, "%.6f" , sbd);
+ x.append(sb);
+
+ x.append("\t");//printing se
+ sbd = male_sumSE(useRandom);
+ sprintf(sb, "%.6f" , sbd);
+ x.append(sb);
+
+
+ x.append("\t");//printing L95 beta
+ sbd = male_Lower95Beta(useRandom);
+ sprintf(sb, "%.6f" , sbd);
+ x.append(sb);
+
+ x.append("\t");//printing U95 beta
+ sbd = male_Upper95Beta(useRandom);
+ sprintf(sb, "%.6f" , sbd);
+ x.append(sb);
+ }
+
+ x.append("\t");//printing z
+ sbd = male_z(useRandom);
+ sprintf(sb, "%.6f" , sbd);
+ x.append(sb);
+
+ x.append("\t");//p-value
+ if (useMetaLambda) sbd = male_pValue(g, useRandom);
+ else sbd = male_pValue(g, useRandom);
+ if (sbd>=0.001) {sprintf(sb, "%.6f" , sbd);}
+ else {sprintf(sb, "%.2E" , sbd);}
+ x.append(sb);
+
+ x.append("\t"); // printing study nr
+ sprintf(sb, "%d" , male_studyCount);
+ x.append(sb);
+
+ x.append("\t"); // printing sample nr
+ if (_Nok)sprintf(sb, "%d" , _Nmale);
+ else sprintf(sb, "%d" , -9);
+ x.append(sb);
+ }
+ else
+ {
+ x.append("\t-9\t-9\t-9\t-9\t-9\t-9\t-9\t-9\t-9");
+ }
+
+///////////FEMALE
+ if (female_studyCount>0)
+ {
+ x.append("\t"); // EAF
+ if (_EAFfemale>0)sprintf(sb, "%.6f" , _EAFfemale);
+ else sprintf(sb, "%d" , -9);
+ x.append(sb);
+
+ if (_bt)
+ {
+ x.append("\t");//printing OR
+ sbd = female_OR(useRandom);
+ sprintf(sb, "%.6f" , sbd);
+ x.append(sb);
+
+ x.append("\t");//printing se
+ sbd = (female_OR(useRandom) - female_Lower95OR(useRandom))/1.96;
+ sprintf(sb, "%.6f" , sbd);
+ x.append(sb);
+
+
+ x.append("\t");//printing L95 OR
+ sbd = female_Lower95OR(useRandom);
+ sprintf(sb, "%.6f" , sbd);
+ x.append(sb);
+
+ x.append("\t");//printing U95 OR
+ sbd = female_Upper95OR(useRandom);
+ sprintf(sb, "%.6f" , sbd);
+ x.append(sb);
+ }
+ else
+ {
+ x.append("\t");//printing average beta
+ sbd = female_sumBeta(useRandom);
+ sprintf(sb, "%.6f" , sbd);
+ x.append(sb);
+
+ x.append("\t");//printing se
+ sbd = female_sumSE(useRandom);
+ sprintf(sb, "%.6f" , sbd);
+ x.append(sb);
+
+
+ x.append("\t");//printing L95 beta
+ sbd = female_Lower95Beta(useRandom);
+ sprintf(sb, "%.6f" , sbd);
+ x.append(sb);
+
+ x.append("\t");//printing U95 beta
+ sbd = female_Upper95Beta(useRandom);
+ sprintf(sb, "%.6f" , sbd);
+ x.append(sb);
+ }
+
+ x.append("\t");//printing z
+ sbd = female_z(useRandom);
+ sprintf(sb, "%.6f" , sbd);
+ x.append(sb);
+
+ x.append("\t");//p-value
+ if (useMetaLambda) sbd = female_pValue(g, useRandom);
+ else sbd = female_pValue(g, useRandom);
+ if (sbd>=0.001) {sprintf(sb, "%.6f" , sbd);}
+ else {sprintf(sb, "%.2E" , sbd);}
+ x.append(sb);
+
+ x.append("\t"); // printing study nr
+ sprintf(sb, "%d" , female_studyCount);
+ x.append(sb);
+
+ x.append("\t"); // printing sample nr
+ if (_Nok)sprintf(sb, "%d" , _Nfemale);
+ else sprintf(sb, "%d" , -9);
+ x.append(sb);
+ }
+ else
+ {
+ x.append("\t-9\t-9\t-9\t-9\t-9\t-9\t-9\t-9\t-9");
+ }
+
+ if (male_studyCount>0 && female_studyCount>0)
+ {
+ x.append("\t");//p-value
+ if (useMetaLambda) sbd = genderdiff_pValue(g, useRandom);
+ else sbd = genderdiff_pValue(g, useRandom);
+ if (sbd>=0.001) {sprintf(sb, "%.6f" , sbd);}
+ else {sprintf(sb, "%.2E" , sbd);}
+ x.append(sb);
+
+ x.append("\t");//p-value
+ if (useMetaLambda) sbd = genderhet_pValue(g, useRandom);
+ else sbd = genderhet_pValue(g, useRandom);
+ if (sbd>=0.001) {sprintf(sb, "%.6f" , sbd);}
+ else {sprintf(sb, "%.2E" , sbd);}
+ x.append(sb);
+
+ }
+ else
+ {
+ x.append("\t-9\t-9");
+ }
+ }
+/////////////////////
+
+
+
+
+
+
+ return x;
+
+}
+
+bool
+marker::doRandom(global & g)
+{
+ double statQ = _Q - _W*((_B / _W)*(_B / _W));
+ double tau2 =0;
+ if ((_W-(_W2/_W))!=0)tau2=(statQ -(_studyCount-1))/(_W-(_W2/_W));
+ if (tau2<0)tau2=0;
+
+ _Wr = 0;
+ _Br = 0;
+ _Qr = 0;
+ for (int i = 0; i < all_se.size(); i++)
+ {
+ double se2;
+ se2 = all_se[i]*all_se[i];
+ _Wr += (1 / (se2 + tau2));
+ _Br += all_beta[i] * (1/(se2 + tau2));
+ _Qr += all_beta[i] * all_beta[i] * (1/se2);
+ }
+
+ if (g._genderHet)
+ {
+/////////////////
+ statQ = _Qmale - _Wmale*((_Bmale / _Wmale)*(_Bmale / _Wmale));
+ tau2 =0;
+ if ((_Wmale-(_W2male/_Wmale))!=0)tau2=(statQ -(male_studyCount-1))/(_Wmale-(_W2male/_Wmale));
+ if (tau2<0)tau2=0;
+
+ _Wrmale = 0;
+ _Brmale = 0;
+ _Qrmale = 0;
+ for (int i = 0; i < male_se.size(); i++)
+ {
+ double se2;
+ se2 = male_se[i]*male_se[i];
+ _Wrmale += (1 / (se2 + tau2));
+ _Brmale += male_beta[i] * (1/(se2 + tau2));
+ _Qrmale += male_beta[i] * male_beta[i] * (1/se2);
+ }
+
+
+
+//////////////////
+ statQ = _Qfemale - _Wfemale*((_Bfemale / _Wfemale)*(_Bfemale / _Wfemale));
+ tau2 =0;
+ if ((_Wfemale-(_W2female/_Wfemale))!=0)tau2=(statQ -(female_studyCount-1))/(_Wfemale-(_W2female/_Wfemale));
+ if (tau2<0)tau2=0;
+
+ _Wrfemale = 0;
+ _Brfemale = 0;
+ _Qrfemale = 0;
+ for (int i = 0; i < female_se.size(); i++)
+ {
+ double se2;
+ se2 = female_se[i]*female_se[i];
+ _Wrfemale += (1 / (se2 + tau2));
+ _Brfemale += female_beta[i] * (1/(se2 + tau2));
+ _Qrfemale += female_beta[i] * female_beta[i] * (1/se2);
+ }
+
+
+/////////////////
+
+
+
+ }
+return true;
+
+}
+
+
+
+double marker::sumBeta(bool useRandom)
+{
+ if (useRandom)return _Br / _Wr;
+ return _B / _W;
+}
+double marker::sumSE(bool useRandom)
+{
+ if (useRandom)return (1/(sqrt(_Wr)));
+ return (1/(sqrt(_W)));
+}
+
+double marker::Lower95Beta(bool useRandom)
+{
+ if (useRandom)return (_Br / _Wr) - (1.96/(sqrt(_Wr)));
+ return (_B / _W) - (1.96/(sqrt(_W)));
+}
+double marker::Upper95Beta(bool useRandom)
+{
+ if (useRandom)return (_Br / _Wr) + (1.96/(sqrt(_Wr)));
+ return (_B / _W) + (1.96/(sqrt(_W)));
+}
+double marker::OR(bool useRandom)
+{
+ if (useRandom)return exp(_Br / _Wr);
+ return exp(_B / _W);
+}
+double marker::Upper95OR(bool useRandom)
+{
+ if (useRandom)return exp((_Br / _Wr) + (1.96/(sqrt(_Wr))));
+ return exp((_B / _W) + (1.96/(sqrt(_W))));
+}
+double marker::Lower95OR(bool useRandom)
+{
+ if (useRandom) return exp((_Br / _Wr) - (1.96/(sqrt(_Wr)))) ;
+ return exp((_B / _W) - (1.96/(sqrt(_W)))) ;
+}
+double marker::z(bool useRandom)
+{
+ if (useRandom) return (_Br / _Wr) * sqrt(_Wr);
+ return (_B / _W) * sqrt(_W);
+}
+double marker::pValue(double outputLambda, bool useRandom)
+{
+ double x;
+ double beta;
+ double se;
+ if (outputLambda>1)
+ {
+ if (useRandom)
+ {
+ beta = _Br / _Wr;
+ se = (1/(sqrt(_Wr)));
+ }
+ else
+ {
+ beta = _B / _W;
+ se = (1/(sqrt(_W)));
+ }
+ se = se * sqrt(outputLambda);
+ x = (beta/se)*(beta/se);
+ }
+ else
+ {
+ if (useRandom)x = ((_Br / _Wr) * sqrt(_Wr))*((_Br / _Wr) * sqrt(_Wr));
+ else x = ((_B / _W) * sqrt(_W))*((_B / _W) * sqrt(_W));
+ }
+ double z = abs_d(x);
+ return bdm_chi2(z, 1);
+ return 1;
+}
+
+double
+marker::genderdiff_pValue(global & g, bool useRandom)
+{
+ double x;
+ double beta;
+ double se;
+ if (g.outputLambda_male>1)
+ {
+ if (useRandom)
+ {
+ beta = _Brmale / _Wrmale;
+ se = (1/(sqrt(_Wrmale)));
+ }
+ else
+ {
+ beta = _Bmale / _Wmale;
+ se = (1/(sqrt(_Wmale)));
+ }
+ se = se * sqrt(g.outputLambda_male);
+ x = (beta/se)*(beta/se);
+ }
+ else
+ {
+ if (useRandom)x = ((_Brmale / _Wrmale) * sqrt(_Wrmale))*((_Brmale / _Wrmale) * sqrt(_Wrmale));
+ else x = ((_Bmale / _Wmale) * sqrt(_Wmale))*((_Bmale / _Wmale) * sqrt(_Wmale));
+ }
+ double z = abs_d(x);
+
+//FEMALE
+ if (g.outputLambda_female>1)
+ {
+ if (useRandom)
+ {
+ beta = _Brfemale / _Wrfemale;
+ se = (1/(sqrt(_Wrfemale)));
+ }
+ else
+ {
+ beta = _Bfemale / _Wfemale;
+ se = (1/(sqrt(_Wfemale)));
+ }
+ se = se * sqrt(g.outputLambda_female);
+ x = (beta/se)*(beta/se);
+ }
+ else
+ {
+ if (useRandom)x = ((_Brfemale / _Wrfemale) * sqrt(_Wrfemale))*((_Brfemale / _Wrfemale) * sqrt(_Wrfemale));
+ else x = ((_Bfemale / _Wfemale) * sqrt(_Wfemale))*((_Bfemale / _Wfemale) * sqrt(_Wfemale));
+ }
+
+ z += abs_d(x);
+ return bdm_chi2(z, 2);
+
+ return 1;
+}
+double
+marker::genderhet_pValue(global & g, bool useRandom)
+{
+ double x;
+ double beta;
+ double se;
+ if (g.outputLambda>1)
+ {
+ if (useRandom)
+ {
+ beta = _Br / _Wr;
+ se = (1/(sqrt(_Wr)));
+ }
+ else
+ {
+ beta = _B / _W;
+ se = (1/(sqrt(_W)));
+ }
+ se = se * sqrt(g.outputLambda);
+ x = (beta/se)*(beta/se);
+ }
+ else
+ {
+ if (useRandom)x = ((_Br / _Wr) * sqrt(_Wr))*((_Br / _Wr) * sqrt(_Wr));
+ else x = ((_B / _W) * sqrt(_W))*((_B / _W) * sqrt(_W));
+ }
+ double z = abs_d(x);
+
+
+ if (g.outputLambda_male>1)
+ {
+ if (useRandom)
+ {
+ beta = _Brmale / _Wrmale;
+ se = (1/(sqrt(_Wrmale)));
+ }
+ else
+ {
+ beta = _Bmale / _Wmale;
+ se = (1/(sqrt(_Wmale)));
+ }
+ se = se * sqrt(g.outputLambda_male);
+ x = (beta/se)*(beta/se);
+ }
+ else
+ {
+ if (useRandom)x = ((_Brmale / _Wrmale) * sqrt(_Wrmale))*((_Brmale / _Wrmale) * sqrt(_Wrmale));
+ else x = ((_Bmale / _Wmale) * sqrt(_Wmale))*((_Bmale / _Wmale) * sqrt(_Wmale));
+ }
+ double z2 = abs_d(x);
+
+//FEMALE
+ if (g.outputLambda_female>1)
+ {
+ if (useRandom)
+ {
+ beta = _Brfemale / _Wrfemale;
+ se = (1/(sqrt(_Wrfemale)));
+ }
+ else
+ {
+ beta = _Bfemale / _Wfemale;
+ se = (1/(sqrt(_Wfemale)));
+ }
+ se = se * sqrt(g.outputLambda_female);
+ x = (beta/se)*(beta/se);
+ }
+ else
+ {
+ if (useRandom)x = ((_Brfemale / _Wrfemale) * sqrt(_Wrfemale))*((_Brfemale / _Wrfemale) * sqrt(_Wrfemale));
+ else x = ((_Bfemale / _Wfemale) * sqrt(_Wfemale))*((_Bfemale / _Wfemale) * sqrt(_Wfemale));
+ }
+
+ z2 += abs_d(x);
+
+ return bdm_chi2(abs_d(z2-z), 1);
+
+ return 1;
+}
+
+
+
+double marker::qStatistic()
+{
+ return _Q - _W*((_B / _W)*(_B / _W));
+}
+double marker::qPValue()
+{
+ double x = _Q - _W*((_B / _W)*(_B / _W));
+ if (x>0 && _studyCount > 1) return 1-chisquaredistribution(_studyCount-1, x);
+ return 1;
+}
+double marker::i2()
+{
+ double h2 = (_Q-_W*((_B / _W)*(_B / _W))) /(_studyCount-1);
+ if (h2<1) return 0;
+ return (h2-1)/h2;
+}
+
+
+////////////////////////////
+//MALE
+double marker::male_sumBeta(bool useRandom)
+{
+ if (useRandom)return _Brmale / _Wrmale;
+ return _Bmale / _Wmale;
+}
+double marker::male_sumSE(bool useRandom)
+{
+ if (useRandom)return (1/(sqrt(_Wrmale)));
+ return (1/(sqrt(_Wmale)));
+}
+
+double marker::male_Lower95Beta(bool useRandom)
+{
+ if (useRandom)return (_Brmale / _Wrmale) - (1.96/(sqrt(_Wrmale)));
+ return (_Bmale / _Wmale) - (1.96/(sqrt(_Wmale)));
+}
+double marker::male_Upper95Beta(bool useRandom)
+{
+ if (useRandom)return (_Brmale / _Wrmale) + (1.96/(sqrt(_Wrmale)));
+ return (_Bmale / _Wmale) + (1.96/(sqrt(_Wmale)));
+}
+double marker::male_OR(bool useRandom)
+{
+ if (useRandom)return exp(_Brmale / _Wrmale);
+ return exp(_Bmale / _Wmale);
+}
+double marker::male_Upper95OR(bool useRandom)
+{
+ if (useRandom)return exp((_Brmale / _Wrmale) + (1.96/(sqrt(_Wrmale))));
+ return exp((_Bmale / _Wmale) + (1.96/(sqrt(_Wmale))));
+}
+double marker::male_Lower95OR(bool useRandom)
+{
+ if (useRandom) return exp((_Brmale / _Wrmale) - (1.96/(sqrt(_Wrmale)))) ;
+ return exp((_Bmale / _Wmale) - (1.96/(sqrt(_Wmale)))) ;
+}
+double marker::male_z(bool useRandom)
+{
+ if (useRandom) return (_Brmale / _Wrmale) * sqrt(_Wrmale);
+ return (_Bmale / _Wmale) * sqrt(_Wmale);
+}
+double marker::male_pValue(global & g, bool useRandom)
+{
+ double x;
+ double beta;
+ double se;
+ double outputLambda = g.outputLambda_male;
+ if (outputLambda>1)
+ {
+ if (useRandom)
+ {
+ beta = _Brmale / _Wrmale;
+ se = (1/(sqrt(_Wrmale)));
+ }
+ else
+ {
+ beta = _Bmale / _Wmale;
+ se = (1/(sqrt(_Wmale)));
+ }
+ se = se * sqrt(outputLambda);
+ x = (beta/se)*(beta/se);
+ }
+ else
+ {
+ if (useRandom)x = ((_Brmale / _Wrmale) * sqrt(_Wrmale))*((_Brmale / _Wrmale) * sqrt(_Wrmale));
+ else x = ((_Bmale / _Wmale) * sqrt(_Wmale))*((_Bmale / _Wmale) * sqrt(_Wmale));
+ }
+ double z = abs_d(x);
+ return bdm_chi2(z, 1);
+ return 1;
+}
+
+////////////////////////////////
+//FEMALE
+
+double marker::female_sumBeta(bool useRandom)
+{
+ if (useRandom)return _Brfemale / _Wrfemale;
+ return _Bfemale / _Wfemale;
+}
+double marker::female_sumSE(bool useRandom)
+{
+ if (useRandom)return (1/(sqrt(_Wrfemale)));
+ return (1/(sqrt(_Wfemale)));
+}
+
+double marker::female_Lower95Beta(bool useRandom)
+{
+ if (useRandom)return (_Brfemale / _Wrfemale) - (1.96/(sqrt(_Wrfemale)));
+ return (_Bfemale / _Wfemale) - (1.96/(sqrt(_Wfemale)));
+}
+double marker::female_Upper95Beta(bool useRandom)
+{
+ if (useRandom)return (_Brfemale / _Wrfemale) + (1.96/(sqrt(_Wrfemale)));
+ return (_Bfemale / _Wfemale) + (1.96/(sqrt(_Wfemale)));
+}
+double marker::female_OR(bool useRandom)
+{
+ if (useRandom)return exp(_Brfemale / _Wrfemale);
+ return exp(_Bfemale / _Wfemale);
+}
+double marker::female_Upper95OR(bool useRandom)
+{
+ if (useRandom)return exp((_Brfemale / _Wrfemale) + (1.96/(sqrt(_Wrfemale))));
+ return exp((_Bfemale / _Wfemale) + (1.96/(sqrt(_Wfemale))));
+}
+double marker::female_Lower95OR(bool useRandom)
+{
+ if (useRandom) return exp((_Brfemale / _Wrfemale) - (1.96/(sqrt(_Wrfemale)))) ;
+ return exp((_Bfemale / _Wfemale) - (1.96/(sqrt(_Wfemale)))) ;
+}
+double marker::female_z(bool useRandom)
+{
+ if (useRandom) return (_Brfemale / _Wrfemale) * sqrt(_Wrfemale);
+ return (_Bfemale / _Wfemale) * sqrt(_Wfemale);
+}
+double marker::female_pValue(global & g, bool useRandom)
+{
+ double x;
+ double beta;
+ double se;
+ double outputLambda = g.outputLambda_female;
+ if (outputLambda>1)
+ {
+ if (useRandom)
+ {
+ beta = _Brfemale / _Wrfemale;
+ se = (1/(sqrt(_Wrfemale)));
+ }
+ else
+ {
+ beta = _Bfemale / _Wfemale;
+ se = (1/(sqrt(_Wfemale)));
+ }
+ se = se * sqrt(outputLambda);
+ x = (beta/se)*(beta/se);
+ }
+ else
+ {
+ if (useRandom)x = ((_Brfemale / _Wrfemale) * sqrt(_Wrfemale))*((_Brfemale / _Wrfemale) * sqrt(_Wrfemale));
+ else x = ((_Bfemale / _Wfemale) * sqrt(_Wfemale))*((_Bfemale / _Wfemale) * sqrt(_Wfemale));
+ }
+ double z = abs_d(x);
+ return bdm_chi2(z, 1);
+ return 1;
+}
diff --git a/marker.h b/marker.h
new file mode 100755
index 0000000..a96b17e
--- /dev/null
+++ b/marker.h
@@ -0,0 +1,151 @@
+/*************************************************************************
+GWAMA software: May, 2009
+
+Contributors:
+ * Andrew P Morris amorris at well.ox.ac.uk
+ * Reedik Magi reedik at well.ox.ac.uk
+
+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.
+*************************************************************************/
+
+
+#pragma once
+#include <iostream>
+#include <map>
+#include <vector>
+#include <stdlib.h>
+#include <fstream>
+#include "study.h"
+#include "global.h"
+#include "problem.h"
+
+using namespace std;
+
+class marker
+{
+
+private:
+ std::string _name;
+
+ std::string _direction;
+ string _effectAllele;
+ string _otherAllele;
+ int _chromosome;
+ int _position;
+ int _N;
+ bool _Nok; //if any population is missing N, then it turns false
+ double _EAF;
+
+ double _B;
+ double _W;
+ double _W2;
+ double _Q;
+ double _Br;
+ double _Wr;
+ double _Qr;
+
+ //gender specific
+ int _Nmale;
+ double _EAFmale;
+
+ double _Bmale;
+ double _Wmale;
+ double _W2male;
+ double _Qmale;
+ double _Brmale;
+ double _Wrmale;
+ double _Qrmale;
+
+ int _Nfemale;
+ double _EAFfemale;
+
+ double _Bfemale;
+ double _Wfemale;
+ double _W2female;
+ double _Qfemale;
+ double _Brfemale;
+ double _Wrfemale;
+ double _Qrfemale;
+
+ vector <double> all_se;
+ vector <double> male_se;
+ vector <double> female_se;
+ vector <double> all_beta;
+ vector <double> male_beta;
+ vector <double> female_beta;
+
+/*
+private:
+ std::vector <bool> _imputed;
+ std::vector <double> _effectAlleleFreq;
+ std::vector <char> _effectAlleleVec;
+ std::vector <char> _nonEffectAlleleVec;
+ std::vector <double> _beta;
+ std::vector <double> _se;
+ std::vector <int> _studyNr;
+ std::vector <int> _n;
+*/
+
+public:
+ int _studyCount;
+ int male_studyCount;
+ int female_studyCount;
+
+ marker(std::string name, int popCount);
+// int addStudy(bool imp, char strand, char effectAllele, char nonEffectAllele, double effectAlleleFreq,
+// double beta, double se, int studyNr,vector <study> & studies,int n, ofstream & ERR, int & errNR, ofstream & LOG);
+ int addCohort(int fileNr, global & g, double directLambda, double imputedLambda,
+ string myMarker,string myEffectAllele, string myOtherAllele, double myEffectAlleleFreq,
+ bool myStrand, double myBeta, double mySE, bool myImputed, int myN, ofstream & ERR, ofstream & LOG, problem & markerProblems);
+
+ bool addMap(int chromosome, int position);
+ bool calculateBWQ(vector <study> &, bool useLambda,bool useRandom, double threshold ,ofstream & ERR, int & errNR, ofstream & LOG);
+ bool doRandom(global & g);
+
+ double sumBeta(bool useRandom);
+ double sumSE(bool useRandom);
+ double Lower95Beta(bool useRandom);
+ double Upper95Beta(bool useRandom);
+ double OR(bool useRandom);
+ double Lower95OR(bool useRandom);
+ double Upper95OR(bool useRandom);
+ double pValue(double outputLambda,bool useRandom);
+ double genderdiff_pValue(global & g,bool useRandom);
+ double genderhet_pValue(global & g,bool useRandom);
+ double male_sumBeta(bool useRandom);
+ double male_sumSE(bool useRandom);
+ double male_Lower95Beta(bool useRandom);
+ double male_Upper95Beta(bool useRandom);
+ double male_OR(bool useRandom);
+ double male_Lower95OR(bool useRandom);
+ double male_Upper95OR(bool useRandom);
+ double male_pValue(global & g,bool useRandom);
+double male_z(bool useRandom);
+ double female_sumBeta(bool useRandom);
+ double female_sumSE(bool useRandom);
+ double female_Lower95Beta(bool useRandom);
+ double female_Upper95Beta(bool useRandom);
+ double female_OR(bool useRandom);
+ double female_Lower95OR(bool useRandom);
+ double female_Upper95OR(bool useRandom);
+ double female_pValue(global & g,bool useRandom);
+double female_z(bool useRandom);
+
+ double z(bool useRandom);
+ double qStatistic();
+ double qPValue();
+ double i2();
+ std::string printMarker(global & g);
+
+ ~marker(void);
+};
diff --git a/normaldistr.cpp b/normaldistr.cpp
new file mode 100755
index 0000000..fcb29c1
--- /dev/null
+++ b/normaldistr.cpp
@@ -0,0 +1,377 @@
+/*************************************************************************
+Cephes Math Library Release 2.8: June, 2000
+Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier
+
+Contributors:
+ * Sergey Bochkanov (ALGLIB project). Translation from C to
+ pseudocode.
+
+See subroutines comments for additional copyrights.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+- Redistributions of source code must retain the above copyright
+ notice, this list of conditions and the following disclaimer.
+
+- Redistributions in binary form must reproduce the above copyright
+ notice, this list of conditions and the following disclaimer listed
+ in this license in the documentation and/or other materials
+ provided with the distribution.
+
+- Neither the name of the copyright holders nor the names of its
+ contributors may 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 <stdafx.h>
+#include "normaldistr.h"
+
+/*************************************************************************
+Error function
+
+The integral is
+
+ x
+ -
+ 2 | | 2
+ erf(x) = -------- | exp( - t ) dt.
+ sqrt(pi) | |
+ -
+ 0
+
+For 0 <= |x| < 1, erf(x) = x * P4(x**2)/Q5(x**2); otherwise
+erf(x) = 1 - erfc(x).
+
+
+ACCURACY:
+
+ Relative error:
+arithmetic domain # trials peak rms
+ IEEE 0,1 30000 3.7e-16 1.0e-16
+
+Cephes Math Library Release 2.8: June, 2000
+Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier
+*************************************************************************/
+double erf(double x)
+{
+ double result;
+ double xsq;
+ double s;
+ double p;
+ double q;
+
+ s = ap::sign(x);
+ x = fabs(x);
+ if( x<0.5 )
+ {
+ xsq = x*x;
+ p = 0.007547728033418631287834;
+ p = 0.288805137207594084924010+xsq*p;
+ p = 14.3383842191748205576712+xsq*p;
+ p = 38.0140318123903008244444+xsq*p;
+ p = 3017.82788536507577809226+xsq*p;
+ p = 7404.07142710151470082064+xsq*p;
+ p = 80437.3630960840172832162+xsq*p;
+ q = 0.0;
+ q = 1.00000000000000000000000+xsq*q;
+ q = 38.0190713951939403753468+xsq*q;
+ q = 658.070155459240506326937+xsq*q;
+ q = 6379.60017324428279487120+xsq*q;
+ q = 34216.5257924628539769006+xsq*q;
+ q = 80437.3630960840172826266+xsq*q;
+ result = s*1.1283791670955125738961589031*x*p/q;
+ return result;
+ }
+ if( x>=10 )
+ {
+ result = s;
+ return result;
+ }
+ result = s*(1-erfc(x));
+ return result;
+}
+
+
+/*************************************************************************
+Complementary error function
+
+ 1 - erf(x) =
+
+ inf.
+ -
+ 2 | | 2
+ erfc(x) = -------- | exp( - t ) dt
+ sqrt(pi) | |
+ -
+ x
+
+
+For small x, erfc(x) = 1 - erf(x); otherwise rational
+approximations are computed.
+
+
+ACCURACY:
+
+ Relative error:
+arithmetic domain # trials peak rms
+ IEEE 0,26.6417 30000 5.7e-14 1.5e-14
+
+Cephes Math Library Release 2.8: June, 2000
+Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier
+*************************************************************************/
+double erfc(double x)
+{
+ double result;
+ double p;
+ double q;
+
+ if( x<0 )
+ {
+ result = 2-erfc(-x);
+ return result;
+ }
+ if( x<0.5 )
+ {
+ result = 1.0-erf(x);
+ return result;
+ }
+ if( x>=10 )
+ {
+ result = 0;
+ return result;
+ }
+ p = 0.0;
+ p = 0.5641877825507397413087057563+x*p;
+ p = 9.675807882987265400604202961+x*p;
+ p = 77.08161730368428609781633646+x*p;
+ p = 368.5196154710010637133875746+x*p;
+ p = 1143.262070703886173606073338+x*p;
+ p = 2320.439590251635247384768711+x*p;
+ p = 2898.0293292167655611275846+x*p;
+ p = 1826.3348842295112592168999+x*p;
+ q = 1.0;
+ q = 17.14980943627607849376131193+x*q;
+ q = 137.1255960500622202878443578+x*q;
+ q = 661.7361207107653469211984771+x*q;
+ q = 2094.384367789539593790281779+x*q;
+ q = 4429.612803883682726711528526+x*q;
+ q = 6089.5424232724435504633068+x*q;
+ q = 4958.82756472114071495438422+x*q;
+ q = 1826.3348842295112595576438+x*q;
+ result = exp(-ap::sqr(x))*p/q;
+ return result;
+}
+
+
+/*************************************************************************
+Normal distribution function
+
+Returns the area under the Gaussian probability density
+function, integrated from minus infinity to x:
+
+ x
+ -
+ 1 | | 2
+ ndtr(x) = --------- | exp( - t /2 ) dt
+ sqrt(2pi) | |
+ -
+ -inf.
+
+ = ( 1 + erf(z) ) / 2
+ = erfc(z) / 2
+
+where z = x/sqrt(2). Computation is via the functions
+erf and erfc.
+
+
+ACCURACY:
+
+ Relative error:
+arithmetic domain # trials peak rms
+ IEEE -13,0 30000 3.4e-14 6.7e-15
+
+Cephes Math Library Release 2.8: June, 2000
+Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier
+*************************************************************************/
+double normaldistribution(double x)
+{
+ double result;
+
+ result = 0.5*(erf(x/1.41421356237309504880)+1);
+ return result;
+}
+
+
+/*************************************************************************
+Inverse of the error function
+
+Cephes Math Library Release 2.8: June, 2000
+Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier
+*************************************************************************/
+double inverf(double e)
+{
+ double result;
+
+ result = invnormaldistribution(0.5*(e+1))/sqrt(double(2));
+ return result;
+}
+
+
+/*************************************************************************
+Inverse of Normal distribution function
+
+Returns the argument, x, for which the area under the
+Gaussian probability density function (integrated from
+minus infinity to x) is equal to y.
+
+
+For small arguments 0 < y < exp(-2), the program computes
+z = sqrt( -2.0 * log(y) ); then the approximation is
+x = z - log(z)/z - (1/z) P(1/z) / Q(1/z).
+There are two rational functions P/Q, one for 0 < y < exp(-32)
+and the other for y up to exp(-2). For larger arguments,
+w = y - 0.5, and x/sqrt(2pi) = w + w**3 R(w**2)/S(w**2)).
+
+ACCURACY:
+
+ Relative error:
+arithmetic domain # trials peak rms
+ IEEE 0.125, 1 20000 7.2e-16 1.3e-16
+ IEEE 3e-308, 0.135 50000 4.6e-16 9.8e-17
+
+Cephes Math Library Release 2.8: June, 2000
+Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier
+*************************************************************************/
+double invnormaldistribution(double y0)
+{
+ double result;
+ double expm2;
+ double s2pi;
+ double x;
+ double y;
+ double z;
+ double y2;
+ double x0;
+ double x1;
+ int code;
+ double p0;
+ double q0;
+ double p1;
+ double q1;
+ double p2;
+ double q2;
+
+ expm2 = 0.13533528323661269189;
+ s2pi = 2.50662827463100050242;
+ if( y0<=0 )
+ {
+ result = -ap::maxrealnumber;
+ return result;
+ }
+ if( y0>=1 )
+ {
+ result = ap::maxrealnumber;
+ return result;
+ }
+ code = 1;
+ y = y0;
+ if( y>1.0-expm2 )
+ {
+ y = 1.0-y;
+ code = 0;
+ }
+ if( y>expm2 )
+ {
+ y = y-0.5;
+ y2 = y*y;
+ p0 = -59.9633501014107895267;
+ p0 = 98.0010754185999661536+y2*p0;
+ p0 = -56.6762857469070293439+y2*p0;
+ p0 = 13.9312609387279679503+y2*p0;
+ p0 = -1.23916583867381258016+y2*p0;
+ q0 = 1;
+ q0 = 1.95448858338141759834+y2*q0;
+ q0 = 4.67627912898881538453+y2*q0;
+ q0 = 86.3602421390890590575+y2*q0;
+ q0 = -225.462687854119370527+y2*q0;
+ q0 = 200.260212380060660359+y2*q0;
+ q0 = -82.0372256168333339912+y2*q0;
+ q0 = 15.9056225126211695515+y2*q0;
+ q0 = -1.18331621121330003142+y2*q0;
+ x = y+y*y2*p0/q0;
+ x = x*s2pi;
+ result = x;
+ return result;
+ }
+ x = sqrt(-2.0*log(y));
+ x0 = x-log(x)/x;
+ z = 1.0/x;
+ if( x<8.0 )
+ {
+ p1 = 4.05544892305962419923;
+ p1 = 31.5251094599893866154+z*p1;
+ p1 = 57.1628192246421288162+z*p1;
+ p1 = 44.0805073893200834700+z*p1;
+ p1 = 14.6849561928858024014+z*p1;
+ p1 = 2.18663306850790267539+z*p1;
+ p1 = -1.40256079171354495875*0.1+z*p1;
+ p1 = -3.50424626827848203418*0.01+z*p1;
+ p1 = -8.57456785154685413611*0.0001+z*p1;
+ q1 = 1;
+ q1 = 15.7799883256466749731+z*q1;
+ q1 = 45.3907635128879210584+z*q1;
+ q1 = 41.3172038254672030440+z*q1;
+ q1 = 15.0425385692907503408+z*q1;
+ q1 = 2.50464946208309415979+z*q1;
+ q1 = -1.42182922854787788574*0.1+z*q1;
+ q1 = -3.80806407691578277194*0.01+z*q1;
+ q1 = -9.33259480895457427372*0.0001+z*q1;
+ x1 = z*p1/q1;
+ }
+ else
+ {
+ p2 = 3.23774891776946035970;
+ p2 = 6.91522889068984211695+z*p2;
+ p2 = 3.93881025292474443415+z*p2;
+ p2 = 1.33303460815807542389+z*p2;
+ p2 = 2.01485389549179081538*0.1+z*p2;
+ p2 = 1.23716634817820021358*0.01+z*p2;
+ p2 = 3.01581553508235416007*0.0001+z*p2;
+ p2 = 2.65806974686737550832*0.000001+z*p2;
+ p2 = 6.23974539184983293730*0.000000001+z*p2;
+ q2 = 1;
+ q2 = 6.02427039364742014255+z*q2;
+ q2 = 3.67983563856160859403+z*q2;
+ q2 = 1.37702099489081330271+z*q2;
+ q2 = 2.16236993594496635890*0.1+z*q2;
+ q2 = 1.34204006088543189037*0.01+z*q2;
+ q2 = 3.28014464682127739104*0.0001+z*q2;
+ q2 = 2.89247864745380683936*0.000001+z*q2;
+ q2 = 6.79019408009981274425*0.000000001+z*q2;
+ x1 = z*p2/q2;
+ }
+ x = x0-x1;
+ if( code!=0 )
+ {
+ x = -x;
+ }
+ result = x;
+ return result;
+}
+
+
+
diff --git a/normaldistr.h b/normaldistr.h
new file mode 100755
index 0000000..fa55d48
--- /dev/null
+++ b/normaldistr.h
@@ -0,0 +1,174 @@
+/*************************************************************************
+Cephes Math Library Release 2.8: June, 2000
+Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier
+
+Contributors:
+ * Sergey Bochkanov (ALGLIB project). Translation from C to
+ pseudocode.
+
+See subroutines comments for additional copyrights.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+- Redistributions of source code must retain the above copyright
+ notice, this list of conditions and the following disclaimer.
+
+- Redistributions in binary form must reproduce the above copyright
+ notice, this list of conditions and the following disclaimer listed
+ in this license in the documentation and/or other materials
+ provided with the distribution.
+
+- Neither the name of the copyright holders nor the names of its
+ contributors may 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 _normaldistr_h
+#define _normaldistr_h
+
+#include "ap.h"
+
+/*************************************************************************
+Error function
+
+The integral is
+
+ x
+ -
+ 2 | | 2
+ erf(x) = -------- | exp( - t ) dt.
+ sqrt(pi) | |
+ -
+ 0
+
+For 0 <= |x| < 1, erf(x) = x * P4(x**2)/Q5(x**2); otherwise
+erf(x) = 1 - erfc(x).
+
+
+ACCURACY:
+
+ Relative error:
+arithmetic domain # trials peak rms
+ IEEE 0,1 30000 3.7e-16 1.0e-16
+
+Cephes Math Library Release 2.8: June, 2000
+Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier
+*************************************************************************/
+double erf(double x);
+
+
+/*************************************************************************
+Complementary error function
+
+ 1 - erf(x) =
+
+ inf.
+ -
+ 2 | | 2
+ erfc(x) = -------- | exp( - t ) dt
+ sqrt(pi) | |
+ -
+ x
+
+
+For small x, erfc(x) = 1 - erf(x); otherwise rational
+approximations are computed.
+
+
+ACCURACY:
+
+ Relative error:
+arithmetic domain # trials peak rms
+ IEEE 0,26.6417 30000 5.7e-14 1.5e-14
+
+Cephes Math Library Release 2.8: June, 2000
+Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier
+*************************************************************************/
+double erfc(double x);
+
+
+/*************************************************************************
+Normal distribution function
+
+Returns the area under the Gaussian probability density
+function, integrated from minus infinity to x:
+
+ x
+ -
+ 1 | | 2
+ ndtr(x) = --------- | exp( - t /2 ) dt
+ sqrt(2pi) | |
+ -
+ -inf.
+
+ = ( 1 + erf(z) ) / 2
+ = erfc(z) / 2
+
+where z = x/sqrt(2). Computation is via the functions
+erf and erfc.
+
+
+ACCURACY:
+
+ Relative error:
+arithmetic domain # trials peak rms
+ IEEE -13,0 30000 3.4e-14 6.7e-15
+
+Cephes Math Library Release 2.8: June, 2000
+Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier
+*************************************************************************/
+double normaldistribution(double x);
+
+
+/*************************************************************************
+Inverse of the error function
+
+Cephes Math Library Release 2.8: June, 2000
+Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier
+*************************************************************************/
+double inverf(double e);
+
+
+/*************************************************************************
+Inverse of Normal distribution function
+
+Returns the argument, x, for which the area under the
+Gaussian probability density function (integrated from
+minus infinity to x) is equal to y.
+
+
+For small arguments 0 < y < exp(-2), the program computes
+z = sqrt( -2.0 * log(y) ); then the approximation is
+x = z - log(z)/z - (1/z) P(1/z) / Q(1/z).
+There are two rational functions P/Q, one for 0 < y < exp(-32)
+and the other for y up to exp(-2). For larger arguments,
+w = y - 0.5, and x/sqrt(2pi) = w + w**3 R(w**2)/S(w**2)).
+
+ACCURACY:
+
+ Relative error:
+arithmetic domain # trials peak rms
+ IEEE 0.125, 1 20000 7.2e-16 1.3e-16
+ IEEE 3e-308, 0.135 50000 4.6e-16 9.8e-17
+
+Cephes Math Library Release 2.8: June, 2000
+Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier
+*************************************************************************/
+double invnormaldistribution(double y0);
+
+
+#endif
diff --git a/problem.cpp b/problem.cpp
new file mode 100755
index 0000000..bc20f90
--- /dev/null
+++ b/problem.cpp
@@ -0,0 +1,15 @@
+#include "problem.h"
+
+problem::problem(void)
+{
+ markersAll = 0; //count of all markers
+ markersOK = 0; //count of ok markers
+ problemStrand = 0; //strand problem
+ wrongAlleles = 0; //wrong alleles problem
+ problemMulti = 0; //multiple occurances of SNP
+ problemEffect = 0;
+}
+
+problem::~problem(void)
+{
+}
diff --git a/problem.h b/problem.h
new file mode 100755
index 0000000..69773ed
--- /dev/null
+++ b/problem.h
@@ -0,0 +1,15 @@
+#pragma once
+
+class problem
+{
+public:
+ int markersAll; //count of all markers
+ int markersOK; //count of ok markers
+ int problemStrand; //strand problem
+ int wrongAlleles; //wrong alleles problem
+ int problemMulti; //multiple occurances of SNP
+ int problemEffect; //effect is missing or problematic
+
+ problem(void);
+ ~problem(void);
+};
diff --git a/readFile.cpp b/readFile.cpp
new file mode 100755
index 0000000..ed15d9f
--- /dev/null
+++ b/readFile.cpp
@@ -0,0 +1,585 @@
+/*************************************************************************
+GWAMA software: May, 2010
+
+Contributors:
+ * Andrew P Morris amorris at well.ox.ac.uk
+ * Reedik Magi reedik at well.ox.ac.uk
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*************************************************************************/
+
+#include <iostream>
+#include <map>
+#include <vector>
+#include <stdlib.h>
+#include <fstream>
+#include "marker.h"
+#include "study.h"
+#include <math.h>
+#include <cctype> // std::toupper
+#include "global.h"
+#include "tools.h"
+#include "cohort.h"
+#include "problem.h"
+
+using namespace std;
+
+
+bool
+readFilelist(global & _g, ofstream & ERR, ofstream & LOG)
+{
+ ifstream IF (_g._fileList.c_str());
+ if (IF.is_open())
+ {
+ LOG << "Reading input file list:" << endl;
+ while (! IF.eof() )
+ {
+ string line;
+ vector<string> tokens;
+ getline (IF,line);
+ int n = Tokenize(string(line), tokens, " ");
+ if (n)
+ {
+ LOG << "\t" << string(tokens[0]) << endl;
+ _g.studyList.push_back(string(tokens[0]));
+ if (_g._genderHet && n>1)
+ {
+ if (uc(string(tokens[1]))=="M"
+ || uc(string(tokens[1]))=="F")
+ _g.studyGenders.push_back(uc(string(tokens[1])));
+ else
+ {
+ cerr <<"Cannot read gender status for file : " << string(tokens[0]) << ". Exit program!" << endl;
+ exit (1);
+ }
+ }
+ else if
+ (_g._genderHet)
+ {
+ cerr <<"Cannot read gender status for file : " << string(tokens[0]) << ". Exit program!" << endl;
+ exit (1);
+
+ }
+ }
+ }
+ LOG << "END-OF-FILE" << endl;
+ }
+ IF.close();
+ _g._studyCount = _g.studyList.size();
+ LOG << "Study count: " << _g._studyCount << endl;
+ cout << "File " <<_g._fileList << " contained "<< _g._studyCount << " studies." << endl;
+ return true;
+}
+
+bool
+getLambda(string name, double & directLambda, double & imputedLambda,
+ int & directCount, int & imputedCount, int columnBeta,
+ int columnSE, int columnImputed, int columnCount, global & g)
+{
+ int currentLine =0;
+ int _imputedCount = 0;
+ int _directCount = 0;
+ vector <double> _imputedChiStat;
+ vector <double> _directChiStat;
+ ifstream F (name.c_str());
+ if (F.is_open())
+ {
+ while (! F.eof() )
+ {
+ string line;
+ vector<string> tokens;
+ getline (F,line);
+ int n = Tokenize(string(line), tokens, " ");
+ if (n)
+ {
+ if (currentLine>0) // skip headers
+ {
+ double beta=0;
+ double se=0;
+ if (g._binaryTrait)
+ {
+ double oratio = atof(tokens[columnBeta].c_str());
+ double oratio_95l = atof(tokens[columnSE].c_str());
+ if (oratio> 0 && oratio > oratio_95l && oratio_95l>0)
+ {
+ beta=log(oratio);
+ se = ((log(oratio)-log(oratio_95l))/1.96);
+ }
+ }
+ else
+ {
+ beta = atof(tokens[columnBeta].c_str());
+ se = atof(tokens[columnSE].c_str());
+ }
+ int imputed = 0;
+ if (columnImputed!=-9)
+ {
+ if (string(tokens[columnImputed])=="1" || string(tokens[columnImputed])=="0")
+ {
+ if (atoi(tokens[columnImputed].c_str())==1)imputed=1;
+ }
+
+ }
+ if (se>0)
+ {
+ if (imputed)
+ {
+ _imputedCount++;
+ _imputedChiStat.push_back((beta*beta)/(se*se));
+ }
+ else
+ {
+ _directCount++;
+ _directChiStat.push_back((beta*beta)/(se*se));
+ }
+ }
+ }
+ }
+ if (int(currentLine/10000)==currentLine/10000.0) cout << "Finding lambda. Line: " << currentLine << "\r";
+
+ currentLine++;
+ }
+ }
+
+ if (_imputedCount>0)
+ {
+ double _medianI = 0;
+ sortVec( _imputedChiStat, _imputedCount);
+ if (_imputedCount%2!=0) _medianI = _imputedChiStat[((_imputedCount+1)/2)-1];
+ else _medianI = (_imputedChiStat[_imputedCount/2 -1] + _imputedChiStat[_imputedCount/2])/2;
+ imputedLambda = _medianI/0.4549364; //median of chi-sq from R ... qchisq(0.5, df= 1)
+
+ }
+
+ if (_directCount>0)
+ {
+ double _medianD = 0;
+ sortVec( _directChiStat, _directCount);
+ if (_directCount%2!=0) _medianD = _directChiStat[((_directCount+1)/2)-1];
+ else _medianD = (_directChiStat[_directCount/2 -1] + _directChiStat[_directCount/2])/2;
+ directLambda = _medianD/0.4549364; //median of chi-sq from R ... qchisq(0.5, df= 1)
+ }
+ directCount = _directCount;
+ imputedCount = _imputedCount;
+
+ return true;
+}
+
+bool
+readMapFile(string markermap, vector <marker> & markerlist,
+ map<string, int> & markerPosition,
+ ofstream & ERR, int & errNR, ofstream & LOG)
+{
+ char sb [11];
+ cout << "Reading map: " << markermap << endl;
+ int currentLine = 0;
+ ifstream F (markermap.c_str());
+ if (F.is_open())
+ {
+ while (! F.eof() )
+ {
+ string line;
+ vector<string> tokens;
+ getline (F,line);
+ int chr = 0;
+ int pos = 0;
+ string currentmarker = "";
+ int n = Tokenize(string(line), tokens, " ");
+ if (n>2)
+ {
+ if (currentLine>0)
+ {
+ chr = atoi(tokens[0].c_str());
+ currentmarker = string(tokens[1]);
+ if (n==3) pos = atoi(tokens[2].c_str());
+ else pos = atoi(tokens[3].c_str());
+ if (markerPosition[currentmarker]!=0)
+ {
+ markerlist[markerPosition[currentmarker]-1].addMap(chr,pos);
+ }
+ }
+ }
+ currentLine ++;
+ }
+ }
+ else {
+ cerr << "Error: Cannot find or access " << markermap << ". Positions not added." << endl;
+ sprintf(sb, "E%.9d" , errNR); errNR++;
+ LOG << "Error: Cannot find or access " << markermap << ". Positions not added! (" << sb << ")" <<endl;
+ ERR << sb << " Error: Cannot find or access " << markermap << ". Positions not added!" <<endl;
+ ERR << sb << " Error: Make sure that this file is in given path and is readable." <<endl;
+ ERR << sb << " Error: Default 'marker.map' file can be downloaded from GWAMA web page." <<endl;
+ return false;
+ }
+ return true;
+}
+
+
+bool
+readCohort(int fileNr, global & g,map<string, int> & markerPosition,vector <marker> & markerlist, ofstream & ERR, ofstream & LOG, ofstream & LAMBDA_FILE)
+{
+ int currentLine = 0;
+ int _countGoodMarkers=0;
+ int _countBadMarkers=0;
+ char sb [11];
+ cohort thisCohort;
+ problem _markerProblems;
+
+ map<string, int> _markernames; //collect all OK marker names for checking for duplicates
+ int _markerNamesCount=1;
+
+ thisCohort._name=g.studyList[fileNr];
+// cout << "Reading file: " << g.studyList[fileNr] << endl;
+ LOG << " [" << fileNr+1 << "]" << " Reading file: " << g.studyList[fileNr] << endl;
+ ifstream F (g.studyList[fileNr].c_str());
+ if (F.is_open())
+ {
+ while (! F.eof() )
+ {
+ string line;
+ vector<string> tokens;
+ getline (F,line);
+ int n = Tokenize(string(line), tokens, " ");
+ if (n)
+ {
+ if (currentLine==0) // get headers
+ {
+ for (int i=0; i<n; i++)
+ {
+ //markername
+ for (unsigned int j=0; j<g._alternativeMarkerNames.size(); j++)
+ {
+ if (uc(string(tokens[i])).compare(g._alternativeMarkerNames[j])==0){thisCohort._columnMarkerName=i;}
+ }
+ //effect allele
+ for (unsigned int j=0; j<g._alternativeEffectAlleles.size(); j++)
+ {
+ if (uc(string(tokens[i])).compare(g._alternativeEffectAlleles[j])==0)thisCohort._columnEffectAllele=i;
+ }
+ //other allele
+ for (unsigned int j=0; j<g._alternativeOtherAlleles.size(); j++)
+ {
+ if (uc(string(tokens[i])).compare(g._alternativeOtherAlleles[j])==0)thisCohort._columnOtherAllele=i;
+ }
+ //effect allele freq
+ for (unsigned int j=0; j<g._alternativeEffectAlleleFreqs.size(); j++)
+ {
+ if (uc(string(tokens[i])).compare(g._alternativeEffectAlleleFreqs[j])==0)thisCohort._columnEffectAlleleFreq=i;
+ }
+ //strand
+ for (unsigned int j=0; j<g._alternativeStrands.size(); j++)
+ {
+ if (uc(string(tokens[i])).compare(g._alternativeStrands[j])==0)thisCohort._columnStrand=i;
+ }
+ //n
+ for (unsigned int j=0; j<g._alternativeNs.size(); j++)
+ {
+ if (uc(string(tokens[i])).compare(g._alternativeNs[j])==0)thisCohort._columnN=i;
+ }
+ //beta
+ for (unsigned int j=0; j<g._alternativeBetas.size(); j++)
+ {
+ if (uc(string(tokens[i])).compare(g._alternativeBetas[j])==0)thisCohort._columnBeta=i;
+ }
+ //se
+ for (unsigned int j=0; j<g._alternativeSEs.size(); j++)
+ {
+ if (uc(string(tokens[i])).compare(g._alternativeSEs[j])==0)thisCohort._columnSE=i;
+ }
+ //or
+ for (unsigned int j=0; j<g._alternativeORs.size(); j++)
+ {
+ if (uc(string(tokens[i])).compare(g._alternativeORs[j])==0)thisCohort._columnOR=i;
+ }
+ //or 95L
+ for (unsigned int j=0; j<g._alternativeOR_95Ls.size(); j++)
+ {
+ if (uc(string(tokens[i])).compare(g._alternativeOR_95Ls[j])==0)thisCohort._columnOR_95L=i;
+ }
+ //or 95U
+ for (unsigned int j=0; j<g._alternativeOR_95Us.size(); j++)
+ {
+ if (uc(string(tokens[i])).compare(g._alternativeOR_95Us[j])==0)thisCohort._columnOR_95U=i;
+ }
+ //imputed
+ for (unsigned int j=0; j<g._alternativeImputeds.size(); j++)
+ {
+ if (uc(string(tokens[i])).compare(g._alternativeImputeds[j])==0)thisCohort._columnImputed=i;
+ }
+
+ }
+ thisCohort._columnCount = n;
+ if (g._noAlleles)
+ {
+ thisCohort._columnEffectAllele=-9;
+ thisCohort._columnOtherAllele=-9;
+ }
+ }//header lines
+ else //data lines
+ {
+ bool currentMarkerIsOK=true;
+ if (thisCohort._columnCount!=n) // lets check if current row has right number of tokens
+ {
+ sprintf(sb, "E%.9d" , g._errNR); g._errNR++;
+ LOG << g.studyList[fileNr] << " row " << currentLine << " has different number of columns compared to header line ( " << n << ", " << thisCohort._columnCount << ") Skipping line!(" << sb << ")" << endl;
+ ERR << sb << " Error: "<< g.studyList[fileNr] << " row " << currentLine << " has different number of columns compared to header line ( " << n << ", " << thisCohort._columnCount << ") Skipping line!" << endl;
+
+ }
+ if (currentLine==1) //lets check if all necessary columns are present
+ {
+ if (thisCohort._columnMarkerName==-9 ||
+ (thisCohort._columnBeta==-9 && g._binaryTrait==false) ||
+ (thisCohort._columnSE==-9 && g._binaryTrait==false) ||
+ (thisCohort._columnOR==-9 && g._binaryTrait==true) ||
+ (thisCohort._columnOR_95L==-9 && g._binaryTrait==true) ||
+ (thisCohort._columnEffectAllele==-9 && g._noAlleles==false) ||
+ (thisCohort._columnOtherAllele==-9 && g._noAlleles==false))
+ {
+ cout << "\tmissing mandatory column! Skipping file." << endl;
+ sprintf(sb, "E%.9d" , g._errNR); g._errNR++;
+ LOG << g.studyList[fileNr] << " misses some of the mandatory columns. Skipping file!(" << sb << ")" << endl;
+ ERR << sb << " " << g.studyList[fileNr] << " misses some of the mandatory columns. Skipping file!" << endl;
+ if (thisCohort._columnMarkerName==-9)ERR << sb << " Markername column is missing" << endl;
+ if (thisCohort._columnBeta==-9 && g._binaryTrait==false)ERR << sb << " Beta column is missing" << endl;
+ if (thisCohort._columnSE==-9 && g._binaryTrait==false)ERR << sb << " SE column is missing" << endl;
+ if (thisCohort._columnOR==-9 && g._binaryTrait==true)ERR << sb << " OR column is missing" << endl;
+ if (thisCohort._columnOR_95L==-9 && g._binaryTrait==true)ERR << sb << " OR _95L column is missing" << endl;
+ if (thisCohort._columnOR_95U==-9 && g._binaryTrait==true)ERR << sb << " OR _95U column is missing" << endl;
+ if (thisCohort._columnEffectAllele==-9 && g._noAlleles==false)ERR << sb << " Effect allele column is missing. If all effects are according to the same allele, then please use --no_alleles option" << endl;
+ if (thisCohort._columnOtherAllele==-9 && g._noAlleles==false)ERR << sb << " Other allele column is missing. If all effects are according to the same allele, then please use --no_alleles option" << endl;
+ return 0;
+ }
+ if (thisCohort._columnStrand==-9)
+ {
+ cout << "Strand column missing! Expecting always positive strand."<< endl;
+ LOG << "Strand column missing! Expecting always positive strand."<< endl;
+ }
+ if (g._genomicControl) //calculating lambdas for the file
+ {
+ bool lambdaSuccess;
+ if (g._binaryTrait==false)lambdaSuccess = getLambda(thisCohort._name, thisCohort._directLambda,
+ thisCohort._imputedLambda, thisCohort._directCount, thisCohort._imputedCount,
+ thisCohort._columnBeta, thisCohort._columnSE, thisCohort._columnImputed,
+ thisCohort._columnCount, g);
+ else lambdaSuccess = getLambda(thisCohort._name, thisCohort._directLambda, thisCohort._imputedLambda,
+ thisCohort._directCount, thisCohort._imputedCount,thisCohort._columnOR, thisCohort._columnOR_95L, thisCohort._columnImputed,
+ thisCohort._columnCount, g);
+
+ if (lambdaSuccess)
+ {
+ cout << "GC lambda genotyped: " << thisCohort._directLambda << " (" << thisCohort._directCount<< ") imputed: " << thisCohort._imputedLambda << " (" << thisCohort._imputedCount << ")"<< endl;
+
+ LAMBDA_FILE << thisCohort._name;
+ char sb [1024]; //temporary char buffer
+ std::string x=""; //starting to collect all information into this string
+
+ x.append("\t");
+ sprintf(sb, "%.4f" , thisCohort._directLambda);
+ x.append(sb);
+
+ x.append("\t");
+ sprintf(sb, "%d" , thisCohort._directCount);
+ x.append(sb);
+
+ x.append("\t");
+ sprintf(sb, "%.4f" , thisCohort._imputedLambda);
+ x.append(sb);
+
+ x.append("\t");
+ sprintf(sb, "%d" , thisCohort._imputedCount);
+ x.append(sb);
+ LAMBDA_FILE << x <<endl;
+
+ }
+ }
+
+ }
+ //everything seems to be OK..lets read marker data
+ //marker name
+ string myMarker = string(tokens[thisCohort._columnMarkerName]);
+ _markerProblems.markersAll++;
+ //alleles
+ string myEffectAllele, myOtherAllele;
+ if (thisCohort._columnEffectAllele!=-9 && thisCohort._columnOtherAllele!=-9)
+ {
+ myEffectAllele = uc(string(tokens[thisCohort._columnEffectAllele]));
+ myOtherAllele = uc(string(tokens[thisCohort._columnOtherAllele]));
+ if (!checkAlleles(myEffectAllele, myOtherAllele) && !g._indel) //problem with alleles - reporting
+ {
+ sprintf(sb, "E%.9d" , g._errNR); g._errNR++;
+ LOG << g.studyList[fileNr] << " has problem with alleles for marker " << myMarker << "!(" << sb << ")" << endl;
+ ERR << sb << " " << g.studyList[fileNr] << " has problem with alleles for marker " << myMarker << "!"<< endl;
+ ERR << sb << " " << g.studyList[fileNr] << " Given alleles: "<< string(tokens[thisCohort._columnEffectAllele]) << "/" << string(tokens[thisCohort._columnOtherAllele]) << ". Skipping marker!" << endl;
+ currentMarkerIsOK=false;
+ _markerProblems.wrongAlleles++;
+ }
+ }
+ else {myEffectAllele = "N";myOtherAllele = "N";}
+ //strand
+ bool myStrand = true;
+ if (thisCohort._columnStrand!=-9){if (string(tokens[thisCohort._columnStrand]).compare("-")==0){myStrand = false;}}
+ if (!myStrand && !g._indel)
+ {
+ myEffectAllele = flip(myEffectAllele);
+ myOtherAllele = flip(myOtherAllele);
+ }
+ //eaf
+ double myEaf=-9;
+ if (thisCohort._columnEffectAlleleFreq!=-9)
+ {
+ if (atof(tokens[thisCohort._columnEffectAlleleFreq].c_str())>0 &&
+ atof(tokens[thisCohort._columnEffectAlleleFreq].c_str())<1)
+ {
+ myEaf = atof(tokens[thisCohort._columnEffectAlleleFreq].c_str());
+ }
+ else //eaf out of range - report it
+ {
+ sprintf(sb, "E%.9d" , g._errNR); g._errNR++;
+ LOG << g.studyList[fileNr] << " has problem with effect allele frequency for marker " << myMarker << "!(" << sb << ")" << endl;
+ ERR << sb << " " << g.studyList[fileNr] << " has problem with effect allele frequency " << myMarker << "!"<< endl;
+ ERR << sb << " " << g.studyList[fileNr] << " Given value: "<< string(tokens[thisCohort._columnEffectAlleleFreq]) << ". Value not used!" << endl;
+ }
+ }
+ //effect+se (or+ci)
+ double myBeta=-9;
+ double mySE=-9;
+ if (g._binaryTrait)
+ {
+ double oratio = atof(tokens[thisCohort._columnOR].c_str());
+ double oratio_95l = atof(tokens[thisCohort._columnOR_95L].c_str());
+ if (oratio>0 && oratio_95l< oratio && oratio_95l>0)
+ {
+ myBeta=log(oratio);
+ mySE = ((log(oratio)-log(oratio_95l))/1.96);
+ }
+ else // problem with or value
+ {
+ sprintf(sb, "E%.9d" , g._errNR); g._errNR++;
+ LOG << g.studyList[fileNr] << " has problem with odds ratio and its 95_L for marker " << myMarker << "!(" << sb << ")" << endl;
+ ERR << sb << " " << g.studyList[fileNr] << " has problem with odds ratio or its CI of " << myMarker << "!"<< endl;
+ ERR << sb << " " << g.studyList[fileNr] << " Given values: OR="<< string(tokens[thisCohort._columnOR]) << " CI_95L="<<string(tokens[thisCohort._columnOR_95L]) << ". Marker not used!" << endl;
+ _markerProblems.problemEffect++;
+ currentMarkerIsOK=false;
+
+ }
+ }
+ else
+ {
+ myBeta = atof(tokens[thisCohort._columnBeta].c_str());
+ mySE = atof(tokens[thisCohort._columnSE].c_str());
+ if (mySE<=0)
+ {
+ sprintf(sb, "E%.9d" , g._errNR); g._errNR++;
+ LOG << g.studyList[fileNr] << " has problem with beta and se for marker " << myMarker << "!(" << sb << ")" << endl;
+ ERR << sb << " " << g.studyList[fileNr] << " has problem with odds ratio or its CI of " << myMarker << "!"<< endl;
+ ERR << sb << " " << g.studyList[fileNr] << " Given values: BETA="<< string(tokens[thisCohort._columnBeta]) << " SE="<<string(tokens[thisCohort._columnSE]) << ". Marker not used!" << endl;
+ _markerProblems.problemEffect++;
+ currentMarkerIsOK=false;
+ }
+ }
+ //imputed
+ bool myImputed = false;
+ if (thisCohort._columnImputed!=-9)
+ {
+ if (string(tokens[thisCohort._columnImputed])=="1" || string(tokens[thisCohort._columnImputed])=="0")
+ {
+ if (atoi(tokens[thisCohort._columnImputed].c_str())==1){myImputed = true;}
+ }
+ else
+ {
+ sprintf(sb, "E%.9d" , g._errNR); g._errNR++;
+ LOG << g.studyList[fileNr] << " has problem with imputation status for marker " << myMarker << "!(" << sb << ")" << endl;
+ ERR << sb << " " << g.studyList[fileNr] << " has problem with imputation status for " << myMarker << "!"<< endl;
+ ERR << sb << " " << g.studyList[fileNr] << " Given values: IMPUTED="<< string(tokens[thisCohort._columnBeta]) << "Value must be 1-imputed, 0-directly genotyped." << endl;
+ }
+ }
+ //n
+ int myN = -9;
+ if (thisCohort._columnN!=-9)
+ {
+ if (atoi(tokens[thisCohort._columnN].c_str())>0){myN = atoi(tokens[thisCohort._columnN].c_str());}
+ else
+ {
+ sprintf(sb, "E%.9d" , g._errNR); g._errNR++;
+ LOG << g.studyList[fileNr] << " has problem with sample size for marker " << myMarker << "!(" << sb << ")" << endl;
+ ERR << sb << " " << g.studyList[fileNr] << " has problem with sample size for " << myMarker << "!"<< endl;
+ ERR << sb << " " << g.studyList[fileNr] << " Given values: N="<< string(tokens[thisCohort._columnN]) << "Value must be larger than zero." << endl;
+ }
+ }
+ //all data read - will add marker to marker list
+ if (currentMarkerIsOK)
+ {
+ if (_markernames[myMarker]==0)
+ {
+ _markernames[myMarker]=_markerNamesCount;
+ _markerNamesCount++;
+
+ if (markerPosition[myMarker]==0)
+ {
+ markerlist.push_back(marker(myMarker, g._studyCount));
+ if (markerlist[g._markerCount-1].addCohort(fileNr, g,
+ thisCohort._directLambda, thisCohort._imputedLambda,
+ myMarker,myEffectAllele, myOtherAllele, myEaf,
+ myStrand, myBeta, mySE, myImputed, myN, ERR, LOG, _markerProblems)==0)
+ // markerlist[markerCount-1].addStudy(currentimputed,
+ // currentstrand, currenteffectallele, currentnoneffectallele, currenteffectallelefreq,
+ // currentbeta, currentse, i, studies,n, ERR, errNR, LOG) ;
+ {
+ markerPosition[myMarker]=g._markerCount;
+ g._markerCount++;
+ }
+ }
+ else //marker is already existing in database
+ {
+ int x = markerPosition[myMarker]; //this is the line number of current marker
+ markerlist[x-1].addCohort(fileNr, g,
+ thisCohort._directLambda, thisCohort._imputedLambda,
+ myMarker,myEffectAllele, myOtherAllele, myEaf,
+ myStrand, myBeta, mySE, myImputed, myN, ERR, LOG, _markerProblems);
+
+ // markerlist[x-1].addStudy(currentimputed,
+ // currentstrand, currenteffectallele, currentnoneffectallele, currenteffectallelefreq,
+ // currentbeta, currentse, i, studies,n, ERR, errNR, LOG);
+
+ }
+ }
+ else
+ {
+ _markerProblems.problemMulti++;
+ }
+ _countGoodMarkers++;
+ }
+ else
+ {
+ _countBadMarkers++;
+ }
+
+ } //data lines end here
+
+ }
+ if (int(currentLine/10000)==currentLine/10000.0) cout << "Line: " << currentLine << "\r";
+ currentLine++;
+ }//while lines
+ F.close();
+ } //if is open
+ else
+ {
+ cerr << "Cannot open file. Exit program!"<< endl;
+ exit(1);
+ }
+ cout << "Marker count: " << _markerProblems.markersAll << " Markers passing sanity check: " << _markerProblems.markersOK << endl;
+ cout << "Strand problems: " << _markerProblems.problemStrand << " Wrong alleles: " << _markerProblems.wrongAlleles << endl;
+ cout << "Effect problems: " << _markerProblems.problemEffect << " Multiple occurances: " << _markerProblems.problemMulti << endl;
+ return true;
+}
+
diff --git a/readFile.h b/readFile.h
new file mode 100755
index 0000000..fcdfa13
--- /dev/null
+++ b/readFile.h
@@ -0,0 +1,50 @@
+/*************************************************************************
+GWAMA software: May, 2010
+
+Contributors:
+ * Andrew P Morris amorris at well.ox.ac.uk
+ * Reedik Magi reedik at well.ox.ac.uk
+
+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 _READFILE_H_
+#define _READFILE_H_
+
+#include <iostream>
+#include <map>
+#include <vector>
+#include <stdlib.h>
+#include <fstream>
+#include "marker.h"
+#include "study.h"
+#include <math.h>
+#include <cctype> // std::toupper
+#include "global.h"
+#include "tools.h"
+
+
+using namespace std;
+
+
+bool readFilelist(global & _g, ofstream & ERR, ofstream & LOG);
+bool
+getLambda(string name, double & directLambda, double & imputedLambda, int & directCount, int & imputedCount, int columnBeta, int columnSE, int columnImputed, int columnCount, global & g);
+bool
+readCohort(int fileNr, global & g,map<string, int> & markerPosition,vector <marker> & markerlist, ofstream & ERR, ofstream & LOG, ofstream & LAMBDA_FILE);
+bool
+readMapFile(string markermap, vector <marker> & markerlist,
+ map<string, int> & markerPosition,
+ ofstream & ERR, int & errNR, ofstream & LOG);
+
+
+#endif
diff --git a/resource.h b/resource.h
new file mode 100755
index 0000000..dba181e
--- /dev/null
+++ b/resource.h
@@ -0,0 +1,14 @@
+//{{NO_DEPENDENCIES}}
+// Microsoft Visual C++ generated include file.
+// Used by gwama_2.rc
+
+// Next default values for new objects
+//
+#ifdef APSTUDIO_INVOKED
+#ifndef APSTUDIO_READONLY_SYMBOLS
+#define _APS_NEXT_RESOURCE_VALUE 101
+#define _APS_NEXT_COMMAND_VALUE 40001
+#define _APS_NEXT_CONTROL_VALUE 1001
+#define _APS_NEXT_SYMED_VALUE 101
+#endif
+#endif
diff --git a/statistics.cpp b/statistics.cpp
new file mode 100755
index 0000000..1800007
--- /dev/null
+++ b/statistics.cpp
@@ -0,0 +1,559 @@
+/*************************************************************************
+GWAMA software: May, 2009
+
+Contributors:
+ * Andrew P Morris amorris at well.ox.ac.uk
+ * Reedik Magi reedik at well.ox.ac.uk
+
+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 __BDM_STATISTICS_CPP__
+
+//
+// SNP Assistant
+//
+// Statistical methods
+//
+// Authors:
+// Reedik M�gi <reedik at ut.ee>
+// Lauris Kaplinski <lauris at kaplinski.com>
+//
+// Copyright (C) 2002-2003 BioData, Ltd.
+//
+
+//
+// This is pure C
+//
+
+#define _USE_MATH_DEFINES
+#include <math.h>
+
+//#include "macros.h"
+//#include "compat.h"
+#include "statistics.h"
+
+
+double
+MAX (double a, double b)
+{
+if (a>b) return a;
+return b;
+}
+
+int
+MAX (int a, int b)
+{
+if (a>b) return a;
+return b;
+}
+
+
+
+
+
+static double bdm_prob_X (int x, int n11, int n12, int n22);
+
+
+double
+bdm_t2x2_chi (double n11, double n12, double n21, double n22)
+{
+ double nt = n11 + n12 + n21 + n22;
+ if (nt <= 1e-9) return -1.0;
+ double rnt = 1.0 / nt;
+ double e11 = rnt * (n11 + n12) * (n11 + n21);
+ double e12 = rnt * (n11 + n12) * (n12 + n22);
+ double e21 = rnt * (n21 + n22) * (n11 + n21);
+ double e22 = rnt * (n21 + n22) * (n12 + n22);
+ double x11 = (e11 != 0.0) ? ((n11 - e11) * (n11 - e11)) / e11 : 0;
+ double x12 = (e12 != 0.0) ? ((n12 - e12) * (n12 - e12)) / e12 : 0;
+ double x21 = (e21 != 0.0) ? ((n21 - e21) * (n21 - e21)) / e21 : 0;
+ double x22 = (e22 != 0.0) ? ((n22 - e22) * (n22 - e22)) / e22 : 0;
+ double xt = x11 + x12 + x21 + x22;
+ if ((e11 < 5.0) || (e12 < 5.0) || (e21 < 5.0) || (e22 < 5.0)) return -1.0;
+ return xt;
+}
+
+double
+bdm_t2x3_chi (double n11, double n12, double n13, double n21, double n22, double n23)
+{
+ double N[3][4], E[2][3];
+ double chi2;
+ int r, c, lowc;
+ // 2x3 Matrix
+ N[0][0] = n11;
+ N[0][1] = n12;
+ N[0][2] = n13;
+ N[1][0] = n21;
+ N[1][1] = n22;
+ N[1][2] = n23;
+ // Rows
+ for (r = 0; r < 2; r++) N[r][3] = N[r][0] + N[r][1] + N[r][2];
+ // Columns
+ for (c = 0; c < 3; c++) N[2][c] = N[0][c] + N[1][c];
+ // Total
+ N[2][3] = N[0][3] + N[1][3];
+ // Expected
+ lowc = 0;
+ for (r = 0; r < 2; r++) {
+ for (c = 0; c < 3; c++) {
+ E[r][c] = (N[r][3] * N[2][c]) / N[2][3];
+ if (E[r][c] <= 0.0) return -1.0;
+ if (E[r][c] < 5.0) lowc += 1;
+ }
+ }
+ if (lowc > 1) return -1.0;
+ chi2 = 0;
+ for (r = 0; r < 2; r++) {
+ for (c = 0; c < 3; c++) {
+ chi2 += (E[r][c] - N[r][c]) * (E[r][c] - N[r][c]) / E[r][c];
+ }
+ }
+ return chi2;
+}
+
+double
+bdm_t2x2_add (int n11, int n12, int n13, int n21, int n22, int n23)
+{
+ double x11 = n13 + n12 / 2.0;
+ double x12 = n11 + n12 / 2.0;
+ double x21 = n23 + n22 / 2.0;
+ double x22 = n21 + n22 / 2.0;
+
+ double x1_ = x11 + x12;
+ double x2_ = x21 + x22;
+ double x_1 = x11 + x21;
+ double x_2 = x12 + x22;
+ double x__ = x_1 + x_2;
+ if (x__ != 0.0) {
+ double expect_x11 = x_1 * x1_ / x__;
+ double expect_x12 = x_2 * x1_ / x__;
+ double expect_x21 = x_1 * x2_ / x__;
+ double expect_x22 = x_2 * x2_ / x__;
+ if ((expect_x11 >= 5.0) && (expect_x12 >= 5.0) && (expect_x21 >= 5.0) && (expect_x22 >= 5.0)) {
+ double vahe = x11 * x22 - x21 * x12;
+ double korrutis = x1_ * x2_ * x_1 * x_2;
+ return (vahe * vahe * x__) / korrutis;
+ }
+ }
+ return -1.0;
+}
+
+double
+bdm_t2x2_hz (int n11, int n12, int n13, int n21, int n22, int n23)
+{
+ double x11 = n11 + n13;
+ double x12 = n12;
+ double x21 = n21 + n23;
+ double x22 = n22;
+
+ double x1_ = x11 + x12;
+ double x2_ = x21 + x22;
+ double x_1 = x11 + x21;
+ double x_2 = x12 + x22;
+ double x__ = x_1 + x_2;
+ if (x__ != 0.0) {
+ double expect_x11 = x_1 * x1_ / x__;
+ double expect_x12 = x_2 * x1_ / x__;
+ double expect_x21 = x_1 * x2_ / x__;
+ double expect_x22 = x_2 * x2_ / x__;
+ if ((expect_x11 >= 5.0) && (expect_x12 >= 5.0) && (expect_x21 >= 5.0) && (expect_x22 >= 5.0)) {
+ double vahe = x11 * x22 - x21 * x12;
+ double korrutis = x1_ * x2_ * x_1 * x_2;
+ return (vahe * vahe * x__) / korrutis;
+ }
+ }
+ return -1.0;
+}
+
+double
+bdm_t2x2_nr (int n11, int n12, int n13, int n21, int n22, int n23)
+{
+ double x11 = n11;
+ double x12 = n12 + n13;
+ double x21 = n21;
+ double x22 = n22 + n23;
+
+ double x1_ = x11 + x12;
+ double x2_ = x21 + x22;
+ double x_1 = x11 + x21;
+ double x_2 = x12 + x22;
+ double x__ = x_1 + x_2;
+ if (x__ != 0.0) {
+ double expect_x11 = x_1 * x1_ / x__;
+ double expect_x12 = x_2 * x1_ / x__;
+ double expect_x21 = x_1 * x2_ / x__;
+ double expect_x22 = x_2 * x2_ / x__;
+ if ((expect_x11 >= 5.0) && (expect_x12 >= 5.0) && (expect_x21 >= 5.0) && (expect_x22 >= 5.0)) {
+ double vahe = x11 * x22 - x21 * x12;
+ double korrutis = x1_ * x2_ * x_1 * x_2;
+ return (vahe * vahe * x__) / korrutis;
+ }
+ }
+ return -1.0;
+}
+
+double
+bdm_t2x2_dr (int n11, int n12, int n13, int n21, int n22, int n23)
+{
+ double x11 = n13;
+ double x12 = n11 + n12;
+ double x21 = n23;
+ double x22 = n21 + n22;
+
+ double x1_ = x11 + x12;
+ double x2_ = x21 + x22;
+ double x_1 = x11 + x21;
+ double x_2 = x12 + x22;
+ double x__ = x_1 + x_2;
+ if (x__ != 0.0) {
+ double expect_x11 = x_1 * x1_ / x__;
+ double expect_x12 = x_2 * x1_ / x__;
+ double expect_x21 = x_1 * x2_ / x__;
+ double expect_x22 = x_2 * x2_ / x__;
+ if ((expect_x11 >= 5.0) && (expect_x12 >= 5.0) && (expect_x21 >= 5.0) && (expect_x22 >= 5.0)) {
+ double vahe = x11 * x22 - x21 * x12;
+ double korrutis = x1_ * x2_ * x_1 * x_2;
+ return (vahe * vahe * x__) / korrutis;
+ }
+ }
+ return -1.0;
+}
+
+double
+bdm_chi2 (double X, int k)
+{
+ if ((X <= 0.0) || (k < 1)) return 1.0;
+
+ if (k == 1) {
+ return 2 * bdm_prob_norm (sqrt (X));
+ }
+
+ if (k <= 30) {
+ if (k & 1) {
+ /* Odd */
+ /* Calculate R */
+ int k_1_2 = k >> 2;
+ double DEN = 1.0;
+ double R = 0.0;
+ int r;
+ for (r = 1; r <= (k_1_2); r++) {
+ DEN *= (2 * r - 1);
+ R += pow (X, r - 0.5) / DEN;
+ }
+ /* Calculate Q */
+ return 2 * bdm_prob_norm (sqrt (X)) + R * sqrt (2 / M_PI) * exp (-X / 2.0);
+ } else {
+ /* Even */
+ int k_2_2 = (k >> 2) - 1;
+ double DEN = 1.0;
+ double R = 0.0;
+ int r;
+ for (r = 1; r < k_2_2; r++) {
+ DEN *= (2 * r);
+ R += pow (X, r) / DEN;
+ }
+ /* Calculate Q */
+ return exp (-X / 2.0) * (1.0 + R);
+ }
+ }
+
+ if (X >= 150.0) {
+ return 0.0;
+ }
+
+ double Z;
+
+ if (X == (k - 1.0)) {
+ Z = -(1 / 3.0 + 0.08 / k) / sqrt (2.0 * k - 2);
+ } else {
+ double d;
+ d = X - k + 2.0 / 3.0 - 0.08 / k;
+ Z = d * sqrt ((k - 1) * log ((k - 1) / X) + X - (k - 1)) / fabs (X - (k - 1));
+ }
+
+ if (Z < 0.0) {
+// return 1 - bdm_prob_norm (Z);
+ return bdm_prob_norm (Z);
+ }
+
+ return bdm_prob_norm (Z);
+}
+
+double
+bdm_prob_norm_our (double x)
+{
+ double p;
+ double xa = fabs (x);
+ if (xa > 12.0) {
+ p = 0;
+ } else {
+ p = 1 / (1 + 0.33267 * xa);
+ p = 0.39894228 * exp (-0.5 * (xa * xa)) * p * (0.4361836 + p * (0.937298 * p - 0.1201676));
+ }
+
+ return (x >= 0.0) ? p : (1.0 - p);
+}
+
+double
+bdm_prob_norm (double x)
+{
+ long double p;
+ long double xa = fabs (x);
+ if (xa > 120.0) {
+ p = 0;
+ } else {
+ p = 1 / (1 + 0.33267 * xa);
+ p = 0.39894228 * exp (-0.5 * (xa * xa)) * p * (0.4361836 + p * (0.937298 * p -
+0.1201676));
+ }
+
+ return (x >= 0.0) ? p : (1.0 - p);
+}
+
+
+
+
+
+
+
+
+
+
+double
+bdm_prob_chi2 (double x2, int ndf)
+{
+ if ((x2 > 0.0) && (ndf > 0)) {
+ if (ndf == 1) {
+ double x = sqrt (x2);
+ return 2 * bdm_prob_norm (x);
+ } else {
+ if (x2 > 169) {
+ double ch = 2.0 / (9.0 * ndf);
+ double x = (exp (log (x2 / ndf) / 3.0) - 1.0 + ch) / sqrt (ch);
+ return bdm_prob_norm (x);
+ } else {
+ if (ndf == 2) {
+ return exp (-0.5 * x2);
+ } else {
+ int N1 = (ndf - 1) / 2;
+ int N2 = (ndf - 2) / 2;
+ if (N1 == N2) {
+ double sum = 1.0;
+ double re = 0.5 * x2;
+ double ch = 1.0;
+ int i;
+ for (i = 1; i <= N2; i++) {
+ ch = ch * re / i;
+ sum += ch;
+ }
+ return exp (-re) * sum;
+ } else {
+ double ch = sqrt (x2);
+ double z = 0.39894228 * exp (-0.5 * x2);
+ double p = bdm_prob_norm (ch);
+ if (ndf == 3) {
+ return 2 * (p + z * ch);
+ } else {
+ double chp = ch;
+ double re = 1.0;
+ int i;
+ for (i = 2; i <= N1; i++) {
+ re += 2;
+ chp *= (x2 / re);
+ ch += chp;
+ }
+ return 2 * (p + z * ch);
+ }
+ }
+ }
+ }
+ }
+ }
+
+ return -1.0;
+}
+
+double
+bdm_signif_min (int f11, int f12, int f21, int f22)
+{
+ static double *L = (double *) 0;
+ static double *N = (double *) 0;
+ static int size = 0;
+
+ int len = 2 * (f11 + f12 + f21 + f22);
+ if (len > size) {
+ // Need to allocate more
+ if (L) delete[] L;
+ if (N) delete[] N;
+ size = 2 * size;
+ size = MAX (size, len);
+ L = new double[size];
+ N = new double[size];
+ }
+
+ int f_1 = f11 + f21;
+ int f_2 = f12 + f22;
+ int f1_ = f11 + f12;
+ int f2_ = f21 + f22;
+ int f__ = f_1 + f_2;
+
+ double prob_f11 = 0.0;
+
+ while ((f22 >= 0.0) && (f11 >= 0.0)) {
+ double prob_f11_liidetav = 1.0;
+ int g = 0;
+ int k;
+ for (k = 0; k < f1_; k++) L[g++] = k + 1;
+ for (k = 0; k < f2_; k++) L[g++] = k + 1;
+ for (k = 0; k < f_2; k++) L[g++] = k + 1;
+ for (k = 0; k < f_1; k++) L[g++] = k + 1;
+
+ g = 0;
+ for (k = 0; k < f__; k++) N[g++] = k + 1;
+ for (k = 0; k < f11; k++) N[g++] = k + 1;
+ for (k = 0; k < f12; k++) N[g++] = k + 1;
+ for (k = 0; k < f21; k++) N[g++] = k + 1;
+ for (k = 0; k < f22; k++) N[g++] = k + 1;
+
+ len = g;
+ for (g = 0; g < len; g++) {
+ double jagatis = L[g] / N[g];
+ prob_f11_liidetav *= jagatis;
+ }
+
+ prob_f11 += prob_f11_liidetav;
+
+ f11 -= 1;
+ f12 += 1;
+ f21 += 1;
+ f22 -= 1;
+ }
+
+ return prob_f11;
+}
+
+double
+bdm_signif_max (int f11, int f12, int f21, int f22)
+{
+ static double *L = (double *) 0;
+ static double *N = (double *) 0;
+ static int size = 0;
+
+ int len = 2 * (f11 + f12 + f21 + f22);
+ if (len > size) {
+ // Need to allocate more
+ if (L) delete[] L;
+ if (N) delete[] N;
+ size = 2 * size;
+ size = MAX (size, len);
+ L = new double[size];
+ N = new double[size];
+ }
+
+ int f_1 = f11 + f21;
+ int f_2 = f12 + f22;
+ int f1_ = f11 + f12;
+ int f2_ = f21 + f22;
+ int f__ = f_1 + f_2;
+
+ double prob_f11 = 0.0;
+
+ while ((f12 >= 0.0) && (f21 >= 0.0)) {
+ double prob_f11_liidetav = 1.0;
+ int g = 0;
+ int k;
+ for (k = 0; k < f1_; k++) L[g++] = k + 1;
+ for (k = 0; k < f2_; k++) L[g++] = k + 1;
+ for (k = 0; k < f_2; k++) L[g++] = k + 1;
+ for (k = 0; k < f_1; k++) L[g++] = k + 1;
+
+ g = 0;
+ for (k = 0; k < f__; k++) N[g++] = k + 1;
+ for (k = 0; k < f11; k++) N[g++] = k + 1;
+ for (k = 0; k < f12; k++) N[g++] = k + 1;
+ for (k = 0; k < f21; k++) N[g++] = k + 1;
+ for (k = 0; k < f22; k++) N[g++] = k + 1;
+
+ len = g;
+ for (g = 0; g < len; g++) {
+ double jagatis = L[g] / N[g];
+ prob_f11_liidetav *= jagatis;
+ }
+
+ prob_f11 += prob_f11_liidetav;
+
+ f11 += 1;
+ f12 -= 1;
+ f21 -= 1;
+ f22 += 1;
+ }
+
+ return prob_f11;
+}
+
+#ifdef BDM_DEBUG
+#include <stdio.h>
+#endif
+
+static double
+bdm_prob_X (int x, int n11, int n12, int n22)
+{
+ int N = n11 + n12 + n22;
+ int N_2 = 2 * N;
+ int c1 = n11 + (n12 - x) / 2;
+ int c2 = n22 + (n12 - x) / 2;
+ int c3 = 2 * n11 + n12;
+ int c4 = N_2 - c3;
+
+ double p = 1.0;
+ for (int k = 1; k <= N_2; k++) {
+ int e = -1;
+ if (k <= N) e += 1;
+ if (k <= c1) e -= 1;
+ if (k <= x) e -= 1;
+ if (k <= c2) e -= 1;
+ if (k <= c3) e += 1;
+ if (k <= c4) e += 1;
+ if (e > 0) {
+ for (int i = 0; i < e; i++) p *= k;
+ } else if (e < 0) {
+ for (int i = 0; i > e; i--) p /= k;
+ }
+ if (k <= x) p *= 2;
+ }
+
+#ifdef BDM_DEBUG
+ double q = p;
+ p = 1.0;
+ for (int k = 1; k <= N_2; k++) {
+ if (k <= N) p *= k;
+ if (k <= c1) p /= k;
+ if (k <= x) p /= k;
+ if (k <= x) p *= 2;
+ if (k <= c2) p /= k;
+ p /= k;
+ if (k <= c3) p *= k;
+ if (k <= c4) p *= k;
+ }
+ if (fabs (1.0 - p/q) > 1e-9) {
+ printf ("ERERERERE\n");
+ }
+#endif
+
+ return p;
+}
+
+
+
diff --git a/statistics.h b/statistics.h
new file mode 100755
index 0000000..0e723eb
--- /dev/null
+++ b/statistics.h
@@ -0,0 +1,52 @@
+/*************************************************************************
+GWAMA software: May, 2009
+
+Contributors:
+ * Andrew P Morris amorris at well.ox.ac.uk
+ * Reedik Magi reedik at well.ox.ac.uk
+
+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 __BDM_STATISTICS_CPP__
+
+//
+// SNP Assistant
+//
+// Statistical methods
+//
+// Authors:
+// Reedik M�gi <reedik at ut.ee>
+// Lauris Kaplinski <lauris at kaplinski.com>
+//
+// Copyright (C) 2002-2003 BioData, Ltd.
+//
+
+//
+// This is pure C
+//
+
+#define _USE_MATH_DEFINES
+#include <math.h>
+
+
+//static double bdm_prob_X (int x, int n11, int n12, int n22);
+double bdm_chi2 (double X, int k);
+double bdm_prob_norm (double x);
+double bdm_prob_chi2 (double x2, int ndf);
+double bdm_signif_min (int f11, int f12, int f21, int f22);
+double bdm_signif_max (int f11, int f12, int f21, int f22);
+static double bdm_prob_X (int x, int n11, int n12, int n22);
+double MAX(double a, double b);
+int MAX(int a, int b);
diff --git a/stdafx.h b/stdafx.h
new file mode 100755
index 0000000..e69de29
diff --git a/study.cpp b/study.cpp
new file mode 100755
index 0000000..3c93d4d
--- /dev/null
+++ b/study.cpp
@@ -0,0 +1,124 @@
+/*************************************************************************
+GWAMA software: May, 2009
+
+Contributors:
+ * Andrew P Morris amorris at well.ox.ac.uk
+ * Reedik Magi reedik at well.ox.ac.uk
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*************************************************************************/
+
+
+#include <iostream>
+#include <map>
+#include <vector>
+#include <stdlib.h>
+#include <fstream>
+#include "study.h"
+#include <algorithm>
+
+using namespace std;
+study::study(std::string name)
+{
+ _name = name;
+ _imputedLambda = 1;
+ _directLambda = 1;
+ _imputedCount = 0;
+ _directCount = 0;
+}
+
+study::~study(void)
+{
+}
+
+std::string
+study::getName(void)
+{
+ return _name;
+
+}
+
+
+void
+study::addImputed(double x)
+{
+ _imputedChiStat.push_back(x);
+ _imputedCount ++;
+}
+
+void
+study::addDirect(double x)
+{
+ _directChiStat.push_back(x);
+ _directCount ++;
+
+}
+bool
+study::calculateLambdas(void)
+{
+ if (_imputedCount>0)
+ {
+ double _medianI = 0;
+ study::sortVec( _imputedChiStat, _imputedCount);
+ if (_imputedCount%2!=0) _medianI = _imputedChiStat[((_imputedCount+1)/2)-1];
+ else _medianI = (_imputedChiStat[_imputedCount/2 -1] + _imputedChiStat[_imputedCount/2])/2;
+ _imputedLambda = _medianI/0.4549364; //median of chi-sq from R ... qchisq(0.5, df= 1)
+
+ }
+
+ if (_directCount>0)
+ {
+ double _medianD = 0;
+ study::sortVec( _directChiStat, _directCount);
+ if (_directCount%2!=0) _medianD = _directChiStat[((_directCount+1)/2)-1];
+ else _medianD = (_directChiStat[_directCount/2 -1] + _directChiStat[_directCount/2])/2;
+ _directLambda = _medianD/0.4549364; //median of chi-sq from R ... qchisq(0.5, df= 1)
+ }
+ return true;
+}
+double
+study::getImputedLambda(void)
+{
+ return _imputedLambda;
+}
+
+double
+study::getDirectLambda(void)
+{
+ return _directLambda;
+}
+
+int
+study::getImputedCount(void)
+{
+ return _imputedCount;
+}
+
+int
+study::getDirectCount(void)
+{
+ return _directCount;
+}
+
+void
+study::printStudy()
+{
+ cout << "Study: " << _name.c_str() << endl;
+
+}
+
+void
+study::sortVec(vector <double>& x, int size)
+{
+ std::sort(x.begin(), x.end());
+}
diff --git a/study.h b/study.h
new file mode 100755
index 0000000..cce62a8
--- /dev/null
+++ b/study.h
@@ -0,0 +1,56 @@
+/*************************************************************************
+GWAMA software: May, 2009
+
+Contributors:
+ * Andrew P Morris amorris at well.ox.ac.uk
+ * Reedik Magi reedik at well.ox.ac.uk
+
+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.
+*************************************************************************/
+
+
+#pragma once
+#include <iostream>
+#include <map>
+#include <vector>
+#include <stdlib.h>
+#include <fstream>
+
+using namespace std;
+class study
+{
+private:
+ int _imputedCount;
+ int _directCount;
+ std::string _name;
+ vector <double> _imputedChiStat;
+ vector <double> _directChiStat;
+ double _imputedLambda;
+ double _directLambda;
+
+public:
+ study(std::string name);
+ ~study(void);
+ void addImputed(double x);
+ void addDirect(double x);
+ bool calculateLambdas(void);
+ double getImputedLambda(void);
+ double getDirectLambda(void);
+
+ int getImputedCount(void);
+ int getDirectCount(void);
+
+ void printStudy(void);
+ void sortVec(vector <double> &vec, int size);
+ std::string getName(void);
+};
diff --git a/tools.cpp b/tools.cpp
new file mode 100755
index 0000000..c94af81
--- /dev/null
+++ b/tools.cpp
@@ -0,0 +1,110 @@
+/*************************************************************************
+GWAMA software: May, 2009
+
+Contributors:
+ * Lauris Kaplinski lauris at kaplinski.com
+ * Reedik Magi reedik at well.ox.ac.uk
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*************************************************************************/
+
+
+#include <iostream>
+#include <map>
+#include <vector>
+#include <stdlib.h>
+#include <fstream>
+#include <math.h>
+#include <cctype> // std::toupper
+#include <zlib.h>
+#include <string>
+#include <algorithm>
+#include "tools.h"
+using namespace std;
+
+
+
+int Tokenize(const string& str1,
+ vector<string>& tokens,
+ const string& delimiters = " ")
+{
+ int cnt = 0;
+ string emptyStr = "";
+ string str = str1;
+ std::replace( str.begin(), str.end(), '\r', ' ' );
+ std::replace( str.begin(), str.end(), '\t', ' ' );
+
+
+
+ // Skip delimiters at beginning.
+ string::size_type lastPos = str.find_first_not_of(delimiters, 0);
+ // Find first "non-delimiter".
+ string::size_type pos = str.find_first_of(delimiters, lastPos);
+
+ while (string::npos != pos || string::npos != lastPos )
+ {
+ if (str.substr(lastPos, pos - lastPos) != emptyStr)
+ {
+ // Found a token, add it to the vector.
+ tokens.push_back(str.substr(lastPos, pos - lastPos));
+ // Skip delimiters. Note the "not_of"
+ lastPos = str.find_first_not_of(delimiters, pos);
+ // Find next "non-delimiter"
+ pos = str.find_first_of(delimiters, lastPos);
+ //
+ cnt++;
+ }
+ }
+ return cnt;
+}
+void
+sortVec(vector <double>& x, int size)
+{
+ std::sort(x.begin(), x.end());
+}
+
+
+
+string uc(string s)
+{
+ const int length = s.length();
+ for(int i=0; i!=length ; ++i)
+ {
+ s[i] = std::toupper(s[i]);
+ }
+ return s;
+}
+
+bool checkAlleles(string & s1, string & s2)
+{
+ if (s1.compare("1")==0) s1 = "A";
+ if (s1.compare("2")==0) s1 = "C";
+ if (s1.compare("3")==0) s1 = "G";
+ if (s1.compare("4")==0) s1 = "T";
+ if (s2.compare("1")==0) s2 = "A";
+ if (s2.compare("2")==0) s2 = "C";
+ if (s2.compare("3")==0) s2 = "G";
+ if (s2.compare("4")==0) s2 = "T";
+ if (!(s1.compare("A")==0 || s1.compare("C")==0 || s1.compare("G")==0 || s1.compare("T")==0))return false;
+ if (!(s2.compare("A")==0 || s2.compare("C")==0 || s2.compare("G")==0 || s2.compare("T")==0))return false;
+ if (s1==s2)return false;
+ return true;
+}
+string flip(string s)
+{
+ if (s.compare("A")==0) return "T";
+ if (s.compare("C")==0) return "G";
+ if (s.compare("G")==0) return "C";
+ if (s.compare("T")==0) return "A";
+ return "N";
+}
diff --git a/tools.h b/tools.h
new file mode 100755
index 0000000..856d417
--- /dev/null
+++ b/tools.h
@@ -0,0 +1,28 @@
+#ifndef _TOOLS_H_
+#define _TOOLS_H_
+
+#include <iostream>
+#include <map>
+#include <vector>
+#include <stdlib.h>
+#include <fstream>
+#include <math.h>
+#include <cctype> // std::toupper
+#include <string>
+#include <algorithm>
+using namespace std;
+
+void sortVec(vector <double>& x, int size);
+
+int Tokenize(const string& str1,
+ vector<string>& tokens,
+ const string& delimiters);
+string uc(string s); //uppercase
+bool checkAlleles(string & s1, string & s2); //check if alleles are ok and change numbers to letters if necessary
+string flip(string s); //flip the alleles if
+
+
+
+
+
+#endif
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-med/gwama.git
More information about the debian-med-commit
mailing list