[Genabel-commits] r1028 - in pkg/ProbABEL: . doc examples src

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Mon Nov 26 07:59:20 CET 2012


Author: lckarssen
Date: 2012-11-26 07:59:20 +0100 (Mon, 26 Nov 2012)
New Revision: 1028

Added:
   pkg/ProbABEL/examples/compare_with_old.sh
   pkg/ProbABEL/examples/compare_with_old_short.sh
   pkg/ProbABEL/examples/test_other/
   pkg/ProbABEL/src/cholesky.cpp
   pkg/ProbABEL/src/cholesky.h
   pkg/ProbABEL/src/command_line_settings.cpp
   pkg/ProbABEL/src/command_line_settings.h
   pkg/ProbABEL/src/coxph_data.cpp
   pkg/ProbABEL/src/coxph_data.h
   pkg/ProbABEL/src/data.cpp
   pkg/ProbABEL/src/eigen_mematrix.cpp
   pkg/ProbABEL/src/eigen_mematrix.h
   pkg/ProbABEL/src/gendata.cpp
   pkg/ProbABEL/src/gendata.h
   pkg/ProbABEL/src/maskedmatrix.cpp
   pkg/ProbABEL/src/maskedmatrix.h
   pkg/ProbABEL/src/phedata.cpp
   pkg/ProbABEL/src/phedata.h
   pkg/ProbABEL/src/reg1.cpp
   pkg/ProbABEL/src/regdata.cpp
   pkg/ProbABEL/src/regdata.h
   pkg/ProbABEL/src/usage.cpp
   pkg/ProbABEL/src/utilities.cpp
   pkg/ProbABEL/src/utilities.h
Modified:
   pkg/ProbABEL/configure.ac
   pkg/ProbABEL/doc/INSTALL
   pkg/ProbABEL/src/Makefile.am
   pkg/ProbABEL/src/chinv2.c
   pkg/ProbABEL/src/chsolve2.c
   pkg/ProbABEL/src/coxfit2.c
   pkg/ProbABEL/src/data.h
   pkg/ProbABEL/src/main.cpp
   pkg/ProbABEL/src/mematri1.h
   pkg/ProbABEL/src/mematrix.h
   pkg/ProbABEL/src/reg1.h
   pkg/ProbABEL/src/survS.h
   pkg/ProbABEL/src/survproto.h
   pkg/ProbABEL/src/testchol.cpp
   pkg/ProbABEL/src/usage.h
Log:
ProbABEL: Merged changes from the ProbABEL-refactoring branch to trunk. 

Thanks to Maarten Kooijman for the large amount of work he put into this branch! 




Modified: pkg/ProbABEL/configure.ac
===================================================================
--- pkg/ProbABEL/configure.ac	2012-11-20 21:57:23 UTC (rev 1027)
+++ pkg/ProbABEL/configure.ac	2012-11-26 06:59:20 UTC (rev 1028)
@@ -14,10 +14,10 @@
 AC_PROG_LN_S
 if test -z "$CXXFLAGS"; then
    # User did not set CXXFLAGS, so we can put in our own defaults
-    CXXFLAGS="-g -O3"
+    CXXFLAGS="-g -O2"
 fi
 if test -z "$CPPFLAGS"; then
-   # User did not set CXXFLAGS, so we can put in our own defaults
+   # User did not set CPPFLAGS, so we can put in our own defaults
     CPPFLAGS="-Wall"
 fi
 # If CXXFLAGS/CPPFLAGS are already set AC_PROG_CXX will not overwrite them
@@ -32,9 +32,45 @@
 
 # Checks for header files.
 AC_FUNC_ALLOCA
-AC_CHECK_HEADERS([float.h inttypes.h libintl.h limits.h stddef.h stdint.h \
-	 		  stdlib.h string.h sys/param.h wchar.h wctype.h])
+AC_CHECK_HEADERS([float.h inttypes.h libintl.h limits.h stddef.h \
+			  stdint.h stdlib.h string.h sys/param.h wchar.h wctype.h])
 
+
+# See if we want use of the Eigen library enabled (yes by default) and if so,
+# whether we can find the library files.
+AC_ARG_ENABLE([eigen],
+ [AS_HELP_STRING([--enable-eigen], [Use the Eigen template library for fast \
+				    linear algebra (Eigen can be downloaded \
+				    from eigen.tuxfamily.org). This is enabled \
+				    by default])],
+ [eigen=${enableval}],
+ [eigen=yes])
+
+if test "x$eigen" = "xyes"; then
+    AC_MSG_NOTICE([building using the Eigen headers enabled])
+
+    AC_ARG_WITH([eigen-include-path],
+	[AS_HELP_STRING([--with-eigen-include-path],
+	  [location of the Eigen headers, defaults to /usr/include/eigen3])],
+	[CXXFLAGS+=" -I${withval} -DEIGEN"],
+	[CXXFLAGS+=' -I/usr/include/eigen3 -DEIGEN'
+	CPPFLAGS+=' -I/usr/include/eigen3'])
+    dnl AC_SUBST([EIGEN_CFLAGS])
+
+    # Check for the EIGEN header files
+    AC_CHECK_HEADERS([Eigen/Dense Eigen/LU])
+
+    if test x$ac_cv_header_Eigen_Dense = xno; then
+      AC_MSG_ERROR([Could not find the Eigen header files. Did you specify \
+--with-eigen-include-path correctly? Or use --disable-eigen \
+to disable use of fast linear algebra.])
+    fi
+else
+   AC_MSG_NOTICE([not using Eigen for linear algebra])
+fi
+AM_CONDITIONAL([WITH_EIGEN], test "x$eigen" = "xyes")
+
+
 # Checks for typedefs, structures, and compiler characteristics.
 AC_HEADER_STDBOOL
 AC_C_INLINE
@@ -44,12 +80,12 @@
 # Checks for library functions.
 AC_FUNC_ERROR_AT_LINE
 AC_FUNC_MALLOC
-AC_CHECK_FUNCS([pow putenv sqrt strdup strncasecmp])
+AC_CHECK_FUNCS([pow putenv sqrt strdup strncasecmp floor])
 
 # Check if we can use large (>2GB) files on 32 bit platforms
 AC_SYS_LARGEFILE
 if test $ac_cv_sys_file_offset_bits != 'no'; then
-  CPPFLAGS="$CPPFLAGS -D_FILE_OFFSET_BITS=$ac_cv_sys_file_offset_bits"
+    CPPFLAGS="$CPPFLAGS -D_FILE_OFFSET_BITS=$ac_cv_sys_file_offset_bits"
 fi
 AC_FUNC_FSEEKO
 

Modified: pkg/ProbABEL/doc/INSTALL
===================================================================
--- pkg/ProbABEL/doc/INSTALL	2012-11-20 21:57:23 UTC (rev 1027)
+++ pkg/ProbABEL/doc/INSTALL	2012-11-26 06:59:20 UTC (rev 1028)
@@ -1,11 +1,30 @@
 These instructions show how to build ProbABEL.
 
+* Dependencies
+  ProbABEL can be compiled without depending on other
+  libraries. However, when the Eigen library is present
+  (http://eigen.tuxfamily.org) matrix manipulation is much more
+  efficient and fast. On Debian/Ubuntu systems the Eigen library can
+  be installed in the following way:
+apt-get install libeigen3-dev
+
+  If you install the eigen library by hand (e.g. after downloading the
+  files from the aforementioned website) you have to specify the
+  location of the header files when running ./configure (see below):
+./configure --with-eigen-include-path=/your/path/to/eigen
+  The default for this option is /usr/include/eigen3. This is the
+  directory in which files like Eigen/Dense and Eigen/Cholesky can be
+  found.
+
+  If ./configure cannot find the Eigen library it will say so. To
+  disable the use of Eigen (even if the library is present) use:
+./configure --disable-eigen
+
+
 * Compiling for Linux
-
-  If you downloaded the source from SVN (this is not
-  necessary when installing from the distributed .tar.gz file), run:
+  If you downloaded the source from SVN (this is not necessary when
+  installing from the distributed .tar.gz file), run:
 autoreconf -i
-
   to generate the ./configure script and some other files.
 
   Now, you are ready to compile and install the package. Run the
@@ -56,15 +75,15 @@
 
   Running
 make uninstall
-
   will uninstall the files previously installed by 'make install'.
 
+
 * Making packages for Linux distributions
   Instructions for creating a package for your favourite Linux
   distribution can be found in doc/packaging.txt
 
+
 * Cross-compiling for Windows on Linux
-
   The following steps show how to compile 32bit windows binaries on a
   Linux machine.
   First, install the necessary MINGW32 packages. On Debian/Ubuntu this
@@ -97,7 +116,7 @@
   subdirectories like 'share/' and 'etc/'.
 
   Now you can run
-make CFLAGS+="-O3 -static-libgcc -static-libstdc++" CXXFLAGS+="-O3 -static-libgcc -static-libstdc++"
+make CFLAGS+="-O2 -static-libgcc -static-libstdc++" CXXFLAGS+="-O2 -static-libgcc -static-libstdc++"
 make install
 
   This creates the binaries pacoxph.exe, palinear.exe and

Copied: pkg/ProbABEL/examples/compare_with_old.sh (from rev 1027, branches/ProbABEL-refactoring/ProbABEL/examples/compare_with_old.sh)
===================================================================
--- pkg/ProbABEL/examples/compare_with_old.sh	                        (rev 0)
+++ pkg/ProbABEL/examples/compare_with_old.sh	2012-11-26 06:59:20 UTC (rev 1028)
@@ -0,0 +1,64 @@
+echo height_base_add.out.txt && diff height_base_add.out.txt  test_other/height_base_add.out.txt
+echo height_base_fv_add.out.txt && diff height_base_fv_add.out.txt  test_other/height_base_fv_add.out.txt
+echo height_allcov_add.out.txt && diff height_allcov_add.out.txt  test_other/height_allcov_add.out.txt
+echo height_allcov_fv_add.out.txt && diff height_allcov_fv_add.out.txt  test_other/height_allcov_fv_add.out.txt
+echo height_int1_add.out.txt && diff height_int1_add.out.txt  test_other/height_int1_add.out.txt
+echo height_int1_fv_add.out.txt && diff height_int1_fv_add.out.txt  test_other/height_int1_fv_add.out.txt
+echo height_robust_add.out.txt && diff height_robust_add.out.txt  test_other/height_robust_add.out.txt
+echo height_robust_fv_add.out.txt && diff height_robust_fv_add.out.txt  test_other/height_robust_fv_add.out.txt
+echo height_robust_int1_add.out.txt && diff height_robust_int1_add.out.txt  test_other/height_robust_int1_add.out.txt
+echo height_robust_int1_fv_add.out.txt && diff height_robust_int1_fv_add.out.txt  test_other/height_robust_int1_fv_add.out.txt
+echo height_ngp2_add.out.txt && diff height_ngp2_add.out.txt  test_other/height_ngp2_add.out.txt
+echo height_ngp2_fv_add.out.txt && diff height_ngp2_fv_add.out.txt  test_other/height_ngp2_fv_add.out.txt
+echo height_ngp2_domin.out.txt && diff height_ngp2_domin.out.txt  test_other/height_ngp2_domin.out.txt
+echo height_ngp2_fv_domin.out.txt && diff height_ngp2_fv_domin.out.txt  test_other/height_ngp2_fv_domin.out.txt
+echo height_ngp2_over_domin.out.txt && diff height_ngp2_over_domin.out.txt  test_other/height_ngp2_over_domin.out.txt
+echo height_ngp2_fv_over_domin.out.txt && diff height_ngp2_fv_over_domin.out.txt  test_other/height_ngp2_fv_over_domin.out.txt
+echo height_ngp2_recess.out.txt && diff height_ngp2_recess.out.txt  test_other/height_ngp2_recess.out.txt
+echo height_ngp2_fv_recess.out.txt && diff height_ngp2_fv_recess.out.txt  test_other/height_ngp2_fv_recess.out.txt
+echo height_ngp2_2df.out.txt && diff height_ngp2_2df.out.txt  test_other/height_ngp2_2df.out.txt
+echo height_ngp2_fv_2df.out.txt && diff height_ngp2_fv_2df.out.txt  test_other/height_ngp2_fv_2df.out.txt
+echo height_ngp2_allcov_add.out.txt && diff height_ngp2_allcov_add.out.txt  test_other/height_ngp2_allcov_add.out.txt
+echo height_ngp2_allcov_fv_add.out.txt && diff height_ngp2_allcov_fv_add.out.txt  test_other/height_ngp2_allcov_fv_add.out.txt
+echo height_ngp2_allcov_domin.out.txt && diff height_ngp2_allcov_domin.out.txt  test_other/height_ngp2_allcov_domin.out.txt
+echo height_ngp2_allcov_fv_domin.out.txt && diff height_ngp2_allcov_fv_domin.out.txt  test_other/height_ngp2_allcov_fv_domin.out.txt
+echo height_ngp2_allcov_over_domin.out.txt && diff height_ngp2_allcov_over_domin.out.txt  test_other/height_ngp2_allcov_over_domin.out.txt
+echo height_ngp2_allcov_fv_over_domin.out.txt && diff height_ngp2_allcov_fv_over_domin.out.txt  test_other/height_ngp2_allcov_fv_over_domin.out.txt
+echo height_ngp2_allcov_recess.out.txt && diff height_ngp2_allcov_recess.out.txt  test_other/height_ngp2_allcov_recess.out.txt
+echo height_ngp2_allcov_fv_recess.out.txt && diff height_ngp2_allcov_fv_recess.out.txt  test_other/height_ngp2_allcov_fv_recess.out.txt
+echo height_ngp2_allcov_2df.out.txt && diff height_ngp2_allcov_2df.out.txt  test_other/height_ngp2_allcov_2df.out.txt
+echo height_ngp2_allcov_fv_2df.out.txt && diff height_ngp2_allcov_fv_2df.out.txt  test_other/height_ngp2_allcov_fv_2df.out.txt
+echo height_ngp2_int1_add.out.txt && diff height_ngp2_int1_add.out.txt  test_other/height_ngp2_int1_add.out.txt
+echo height_ngp2_int1_fv_add.out.txt && diff height_ngp2_int1_fv_add.out.txt  test_other/height_ngp2_int1_fv_add.out.txt
+echo height_ngp2_int1_domin.out.txt && diff height_ngp2_int1_domin.out.txt  test_other/height_ngp2_int1_domin.out.txt
+echo height_ngp2_int1_fv_domin.out.txt && diff height_ngp2_int1_fv_domin.out.txt  test_other/height_ngp2_int1_fv_domin.out.txt
+echo height_ngp2_int1_over_domin.out.txt && diff height_ngp2_int1_over_domin.out.txt  test_other/height_ngp2_int1_over_domin.out.txt
+echo height_ngp2_int1_fv_over_domin.out.txt && diff height_ngp2_int1_fv_over_domin.out.txt  test_other/height_ngp2_int1_fv_over_domin.out.txt
+echo height_ngp2_int1_recess.out.txt && diff height_ngp2_int1_recess.out.txt  test_other/height_ngp2_int1_recess.out.txt
+echo height_ngp2_int1_fv_recess.out.txt && diff height_ngp2_int1_fv_recess.out.txt  test_other/height_ngp2_int1_fv_recess.out.txt
+echo height_ngp2_int1_2df.out.txt && diff height_ngp2_int1_2df.out.txt  test_other/height_ngp2_int1_2df.out.txt
+echo height_ngp2_int1_fv_2df.out.txt && diff height_ngp2_int1_fv_2df.out.txt  test_other/height_ngp2_int1_fv_2df.out.txt
+echo height_ngp2_robust_add.out.txt && diff height_ngp2_robust_add.out.txt  test_other/height_ngp2_robust_add.out.txt
+echo height_ngp2_robust_fv_add.out.txt && diff height_ngp2_robust_fv_add.out.txt  test_other/height_ngp2_robust_fv_add.out.txt
+echo height_ngp2_robust_domin.out.txt && diff height_ngp2_robust_domin.out.txt  test_other/height_ngp2_robust_domin.out.txt
+echo height_ngp2_robust_fv_domin.out.txt && diff height_ngp2_robust_fv_domin.out.txt  test_other/height_ngp2_robust_fv_domin.out.txt
+echo height_ngp2_robust_over_domin.out.txt && diff height_ngp2_robust_over_domin.out.txt  test_other/height_ngp2_robust_over_domin.out.txt
+echo height_ngp2_robust_fv_over_domin.out.txt && diff height_ngp2_robust_fv_over_domin.out.txt  test_other/height_ngp2_robust_fv_over_domin.out.txt
+echo height_ngp2_robust_recess.out.txt && diff height_ngp2_robust_recess.out.txt  test_other/height_ngp2_robust_recess.out.txt
+echo height_ngp2_robust_fv_recess.out.txt && diff height_ngp2_robust_fv_recess.out.txt  test_other/height_ngp2_robust_fv_recess.out.txt
+echo height_ngp2_robust_2df.out.txt && diff height_ngp2_robust_2df.out.txt  test_other/height_ngp2_robust_2df.out.txt
+echo height_ngp2_robust_fv_2df.out.txt && diff height_ngp2_robust_fv_2df.out.txt  test_other/height_ngp2_robust_fv_2df.out.txt
+echo height_ngp2_robust_int1_add.out.txt && diff height_ngp2_robust_int1_add.out.txt  test_other/height_ngp2_robust_int1_add.out.txt
+echo height_ngp2_robust_int1_fv_add.out.txt && diff height_ngp2_robust_int1_fv_add.out.txt  test_other/height_ngp2_robust_int1_fv_add.out.txt
+echo height_ngp2_robust_int1_domin.out.txt && diff height_ngp2_robust_int1_domin.out.txt  test_other/height_ngp2_robust_int1_domin.out.txt
+echo height_ngp2_robust_int1_fv_domin.out.txt && diff height_ngp2_robust_int1_fv_domin.out.txt  test_other/height_ngp2_robust_int1_fv_domin.out.txt
+echo height_ngp2_robust_int1_over_domin.out.txt && diff height_ngp2_robust_int1_over_domin.out.txt  test_other/height_ngp2_robust_int1_over_domin.out.txt
+echo height_ngp2_robust_int1_fv_over_domin.out.txt && diff height_ngp2_robust_int1_fv_over_domin.out.txt  test_other/height_ngp2_robust_int1_fv_over_domin.out.txt
+echo height_ngp2_robust_int1_recess.out.txt && diff height_ngp2_robust_int1_recess.out.txt  test_other/height_ngp2_robust_int1_recess.out.txt
+echo height_ngp2_robust_int1_fv_recess.out.txt && diff height_ngp2_robust_int1_fv_recess.out.txt  test_other/height_ngp2_robust_int1_fv_recess.out.txt
+echo height_ngp2_robust_int1_2df.out.txt && diff height_ngp2_robust_int1_2df.out.txt  test_other/height_ngp2_robust_int1_2df.out.txt
+echo height_ngp2_robust_int1_fv_2df.out.txt && diff height_ngp2_robust_int1_fv_2df.out.txt  test_other/height_ngp2_robust_int1_fv_2df.out.txt
+echo mmscore_add.out.txt && diff mmscore_add.out.txt  test_other/mmscore_add.out.txt
+echo mmscore_fv_add.out.txt && diff mmscore_fv_add.out.txt test_other/mmscore_fv_add.out.txt
+echo mmscore_prob_add.out.txt && diff mmscore_prob_add.out.txt test_other/mmscore_prob_add.out.txt 
+echo mmscore_prob_fv_add.out.txt && diff mmscore_prob_fv_add.out.txt test_other/mmscore_prob_fv_add.out.txt

Copied: pkg/ProbABEL/examples/compare_with_old_short.sh (from rev 1027, branches/ProbABEL-refactoring/ProbABEL/examples/compare_with_old_short.sh)
===================================================================
--- pkg/ProbABEL/examples/compare_with_old_short.sh	                        (rev 0)
+++ pkg/ProbABEL/examples/compare_with_old_short.sh	2012-11-26 06:59:20 UTC (rev 1028)
@@ -0,0 +1,60 @@
+echo height_base_add.out.txt && diff height_base_add.out.txt  test_other/height_base_add.out.txt |wc -l
+echo height_base_fv_add.out.txt && diff height_base_fv_add.out.txt  test_other/height_base_fv_add.out.txt |wc -l
+echo height_allcov_add.out.txt && diff height_allcov_add.out.txt  test_other/height_allcov_add.out.txt |wc -l
+echo height_allcov_fv_add.out.txt && diff height_allcov_fv_add.out.txt  test_other/height_allcov_fv_add.out.txt |wc -l
+echo height_int1_add.out.txt && diff height_int1_add.out.txt  test_other/height_int1_add.out.txt |wc -l
+echo height_int1_fv_add.out.txt && diff height_int1_fv_add.out.txt  test_other/height_int1_fv_add.out.txt |wc -l
+echo height_robust_add.out.txt && diff height_robust_add.out.txt  test_other/height_robust_add.out.txt |wc -l
+echo height_robust_fv_add.out.txt && diff height_robust_fv_add.out.txt  test_other/height_robust_fv_add.out.txt |wc -l
+echo height_robust_int1_add.out.txt && diff height_robust_int1_add.out.txt  test_other/height_robust_int1_add.out.txt |wc -l
+echo height_robust_int1_fv_add.out.txt && diff height_robust_int1_fv_add.out.txt  test_other/height_robust_int1_fv_add.out.txt |wc -l
+echo height_ngp2_add.out.txt && diff height_ngp2_add.out.txt  test_other/height_ngp2_add.out.txt |wc -l
+echo height_ngp2_fv_add.out.txt && diff height_ngp2_fv_add.out.txt  test_other/height_ngp2_fv_add.out.txt |wc -l
+echo height_ngp2_domin.out.txt && diff height_ngp2_domin.out.txt  test_other/height_ngp2_domin.out.txt |wc -l
+echo height_ngp2_fv_domin.out.txt && diff height_ngp2_fv_domin.out.txt  test_other/height_ngp2_fv_domin.out.txt |wc -l
+echo height_ngp2_over_domin.out.txt && diff height_ngp2_over_domin.out.txt  test_other/height_ngp2_over_domin.out.txt |wc -l
+echo height_ngp2_fv_over_domin.out.txt && diff height_ngp2_fv_over_domin.out.txt  test_other/height_ngp2_fv_over_domin.out.txt |wc -l
+echo height_ngp2_recess.out.txt && diff height_ngp2_recess.out.txt  test_other/height_ngp2_recess.out.txt |wc -l
+echo height_ngp2_fv_recess.out.txt && diff height_ngp2_fv_recess.out.txt  test_other/height_ngp2_fv_recess.out.txt |wc -l
+echo height_ngp2_2df.out.txt && diff height_ngp2_2df.out.txt  test_other/height_ngp2_2df.out.txt |wc -l
+echo height_ngp2_fv_2df.out.txt && diff height_ngp2_fv_2df.out.txt  test_other/height_ngp2_fv_2df.out.txt |wc -l
+echo height_ngp2_allcov_add.out.txt && diff height_ngp2_allcov_add.out.txt  test_other/height_ngp2_allcov_add.out.txt |wc -l
+echo height_ngp2_allcov_fv_add.out.txt && diff height_ngp2_allcov_fv_add.out.txt  test_other/height_ngp2_allcov_fv_add.out.txt |wc -l
+echo height_ngp2_allcov_domin.out.txt && diff height_ngp2_allcov_domin.out.txt  test_other/height_ngp2_allcov_domin.out.txt |wc -l
+echo height_ngp2_allcov_fv_domin.out.txt && diff height_ngp2_allcov_fv_domin.out.txt  test_other/height_ngp2_allcov_fv_domin.out.txt |wc -l
+echo height_ngp2_allcov_over_domin.out.txt && diff height_ngp2_allcov_over_domin.out.txt  test_other/height_ngp2_allcov_over_domin.out.txt |wc -l
+echo height_ngp2_allcov_fv_over_domin.out.txt && diff height_ngp2_allcov_fv_over_domin.out.txt  test_other/height_ngp2_allcov_fv_over_domin.out.txt |wc -l
+echo height_ngp2_allcov_recess.out.txt && diff height_ngp2_allcov_recess.out.txt  test_other/height_ngp2_allcov_recess.out.txt |wc -l
+echo height_ngp2_allcov_fv_recess.out.txt && diff height_ngp2_allcov_fv_recess.out.txt  test_other/height_ngp2_allcov_fv_recess.out.txt |wc -l
+echo height_ngp2_allcov_2df.out.txt && diff height_ngp2_allcov_2df.out.txt  test_other/height_ngp2_allcov_2df.out.txt |wc -l
+echo height_ngp2_allcov_fv_2df.out.txt && diff height_ngp2_allcov_fv_2df.out.txt  test_other/height_ngp2_allcov_fv_2df.out.txt |wc -l
+echo height_ngp2_int1_add.out.txt && diff height_ngp2_int1_add.out.txt  test_other/height_ngp2_int1_add.out.txt |wc -l
+echo height_ngp2_int1_fv_add.out.txt && diff height_ngp2_int1_fv_add.out.txt  test_other/height_ngp2_int1_fv_add.out.txt |wc -l
+echo height_ngp2_int1_domin.out.txt && diff height_ngp2_int1_domin.out.txt  test_other/height_ngp2_int1_domin.out.txt |wc -l
+echo height_ngp2_int1_fv_domin.out.txt && diff height_ngp2_int1_fv_domin.out.txt  test_other/height_ngp2_int1_fv_domin.out.txt |wc -l
+echo height_ngp2_int1_over_domin.out.txt && diff height_ngp2_int1_over_domin.out.txt  test_other/height_ngp2_int1_over_domin.out.txt |wc -l
+echo height_ngp2_int1_fv_over_domin.out.txt && diff height_ngp2_int1_fv_over_domin.out.txt  test_other/height_ngp2_int1_fv_over_domin.out.txt |wc -l
+echo height_ngp2_int1_recess.out.txt && diff height_ngp2_int1_recess.out.txt  test_other/height_ngp2_int1_recess.out.txt |wc -l
+echo height_ngp2_int1_fv_recess.out.txt && diff height_ngp2_int1_fv_recess.out.txt  test_other/height_ngp2_int1_fv_recess.out.txt |wc -l
+echo height_ngp2_int1_2df.out.txt && diff height_ngp2_int1_2df.out.txt  test_other/height_ngp2_int1_2df.out.txt |wc -l
+echo height_ngp2_int1_fv_2df.out.txt && diff height_ngp2_int1_fv_2df.out.txt  test_other/height_ngp2_int1_fv_2df.out.txt |wc -l
+echo height_ngp2_robust_add.out.txt && diff height_ngp2_robust_add.out.txt  test_other/height_ngp2_robust_add.out.txt |wc -l
+echo height_ngp2_robust_fv_add.out.txt && diff height_ngp2_robust_fv_add.out.txt  test_other/height_ngp2_robust_fv_add.out.txt |wc -l
+echo height_ngp2_robust_domin.out.txt && diff height_ngp2_robust_domin.out.txt  test_other/height_ngp2_robust_domin.out.txt |wc -l
+echo height_ngp2_robust_fv_domin.out.txt && diff height_ngp2_robust_fv_domin.out.txt  test_other/height_ngp2_robust_fv_domin.out.txt |wc -l
+echo height_ngp2_robust_over_domin.out.txt && diff height_ngp2_robust_over_domin.out.txt  test_other/height_ngp2_robust_over_domin.out.txt |wc -l
+echo height_ngp2_robust_fv_over_domin.out.txt && diff height_ngp2_robust_fv_over_domin.out.txt  test_other/height_ngp2_robust_fv_over_domin.out.txt |wc -l
+echo height_ngp2_robust_recess.out.txt && diff height_ngp2_robust_recess.out.txt  test_other/height_ngp2_robust_recess.out.txt |wc -l
+echo height_ngp2_robust_fv_recess.out.txt && diff height_ngp2_robust_fv_recess.out.txt  test_other/height_ngp2_robust_fv_recess.out.txt |wc -l
+echo height_ngp2_robust_2df.out.txt && diff height_ngp2_robust_2df.out.txt  test_other/height_ngp2_robust_2df.out.txt |wc -l
+echo height_ngp2_robust_fv_2df.out.txt && diff height_ngp2_robust_fv_2df.out.txt  test_other/height_ngp2_robust_fv_2df.out.txt |wc -l
+echo height_ngp2_robust_int1_add.out.txt && diff height_ngp2_robust_int1_add.out.txt  test_other/height_ngp2_robust_int1_add.out.txt |wc -l
+echo height_ngp2_robust_int1_fv_add.out.txt && diff height_ngp2_robust_int1_fv_add.out.txt  test_other/height_ngp2_robust_int1_fv_add.out.txt |wc -l
+echo height_ngp2_robust_int1_domin.out.txt && diff height_ngp2_robust_int1_domin.out.txt  test_other/height_ngp2_robust_int1_domin.out.txt |wc -l
+echo height_ngp2_robust_int1_fv_domin.out.txt && diff height_ngp2_robust_int1_fv_domin.out.txt  test_other/height_ngp2_robust_int1_fv_domin.out.txt |wc -l
+echo height_ngp2_robust_int1_over_domin.out.txt && diff height_ngp2_robust_int1_over_domin.out.txt  test_other/height_ngp2_robust_int1_over_domin.out.txt |wc -l
+echo height_ngp2_robust_int1_fv_over_domin.out.txt && diff height_ngp2_robust_int1_fv_over_domin.out.txt  test_other/height_ngp2_robust_int1_fv_over_domin.out.txt |wc -l
+echo height_ngp2_robust_int1_recess.out.txt && diff height_ngp2_robust_int1_recess.out.txt  test_other/height_ngp2_robust_int1_recess.out.txt |wc -l
+echo height_ngp2_robust_int1_fv_recess.out.txt && diff height_ngp2_robust_int1_fv_recess.out.txt  test_other/height_ngp2_robust_int1_fv_recess.out.txt |wc -l
+echo height_ngp2_robust_int1_2df.out.txt && diff height_ngp2_robust_int1_2df.out.txt  test_other/height_ngp2_robust_int1_2df.out.txt |wc -l
+echo height_ngp2_robust_int1_fv_2df.out.txt && diff height_ngp2_robust_int1_fv_2df.out.txt  test_other/height_ngp2_robust_int1_fv_2df.out.txt |wc -l
\ No newline at end of file

Modified: pkg/ProbABEL/src/Makefile.am
===================================================================
--- pkg/ProbABEL/src/Makefile.am	2012-11-20 21:57:23 UTC (rev 1027)
+++ pkg/ProbABEL/src/Makefile.am	2012-11-26 06:59:20 UTC (rev 1028)
@@ -3,7 +3,12 @@
 ## Using wildcards in these lists doesn't work. Also GNU make's ($wildcard,) doesn't
 ## work. It gives warning message about portability, but in the end doesn't work,
 ## I tried :-).
-REGFILES = data.h mematrix.h reg1.h usage.h main.cpp mematri1.h
+REGFILES = data.h data.cpp gendata.h gendata.cpp mematrix.h mematri1.h	\
+ command_line_settings.h command_line_settings.cpp reg1.h usage.h		\
+ usage.cpp main.cpp utilities.h utilities.cpp phedata.h phedata.cpp	\
+ cholesky.h cholesky.cpp regdata.h regdata.cpp 	\
+ maskedmatrix.cpp maskedmatrix.h eigen_mematrix.h 	\
+ eigen_mematrix.cpp reg1.cpp
 
 ## Filevector files:
 FVSRC = fvlib/AbstractMatrix.cpp fvlib/CastUtils.cpp			\
@@ -42,7 +47,7 @@
 palogist_CPPFLAGS = $(AM_CPPFLAGS)
 
 pacoxph_SOURCES = $(COXSRC) $(REGFILES) $(FVSRC) $(FVHEADERS)	\
- $(RHEADERS) survS.h survproto.h
+ $(RHEADERS) survS.h survproto.h coxph_data.h coxph_data.cpp
 pacoxph_CXXFLAGS = -DCOXPH -I $(top_srcdir)/src/include $(AM_CXXFLAGS)
 pacoxph_CPPFLAGS = $(AM_CPPFLAGS)
 pacoxph_CFLAGS = -DCOXPH -I $(top_srcdir)/src/include $(AM_CFLAGS)

Modified: pkg/ProbABEL/src/chinv2.c
===================================================================
--- pkg/ProbABEL/src/chinv2.c	2012-11-20 21:57:23 UTC (rev 1027)
+++ pkg/ProbABEL/src/chinv2.c	2012-11-26 06:59:20 UTC (rev 1028)
@@ -1,54 +1,64 @@
 /*  SCCS @(#)chinv2.c	5.3 07/15/99
-** matrix inversion, given the FDF' cholesky decomposition
-**
-** input  **matrix, which contains the chol decomp of an n by n
-**   matrix in its lower triangle.
-**
-** returned: the upper triangle + diagonal contain (FDF')^{-1}
-**            below the diagonal will be F inverse
-**
-**  Terry Therneau
-*/
+ ** matrix inversion, given the FDF' cholesky decomposition
+ **
+ ** input  **matrix, which contains the chol decomp of an n by n
+ **   matrix in its lower triangle.
+ **
+ ** returned: the upper triangle + diagonal contain (FDF')^{-1}
+ **            below the diagonal will be F inverse
+ **
+ **  Terry Therneau
+ */
 #include "survS.h"
 #include "survproto.h"
 
-void chinv2(double **matrix , int n)
-     {
-     register double temp;
-     register int i,j,k;
+void chinv2(double **matrix, int n)
+{
+    register double temp;
+    register int i, j, k;
 
-     /*
+    /*
      ** invert the cholesky in the lower triangle
      **   take full advantage of the cholesky's diagonal of 1's
      */
-     for (i=0; i<n; i++){
-	  if (matrix[i][i] >0) {
-	      matrix[i][i] = 1/matrix[i][i];   /*this line inverts D */
-	      for (j= (i+1); j<n; j++) {
-		   matrix[j][i] = -matrix[j][i];
-		   for (k=0; k<i; k++)     /*sweep operator */
-			matrix[j][k] += matrix[j][i]*matrix[i][k];
-		   }
-	      }
-	  }
+    for (i = 0; i < n; i++)
+    {
+        if (matrix[i][i] > 0)
+        {
+            matrix[i][i] = 1 / matrix[i][i]; /*this line inverts D */
+            for (j = (i + 1); j < n; j++)
+            {
+                matrix[j][i] = -matrix[j][i];
+                for (k = 0; k < i; k++) /*sweep operator */
+                    matrix[j][k] += matrix[j][i] * matrix[i][k];
+            }
+        }
+    }
 
-     /*
+    /*
      ** lower triangle now contains inverse of cholesky
      ** calculate F'DF (inverse of cholesky decomp process) to get inverse
      **   of original matrix
      */
-     for (i=0; i<n; i++) {
-	  if (matrix[i][i]==0) {  /* singular row */
-		for (j=0; j<i; j++) matrix[j][i]=0;
-		for (j=i; j<n; j++) matrix[i][j]=0;
-		}
-	  else {
-	      for (j=(i+1); j<n; j++) {
-		   temp = matrix[j][i]*matrix[j][j];
-		   if (j!=i) matrix[i][j] = temp;
-		   for (k=i; k<j; k++)
-			matrix[i][k] += temp*matrix[j][k];
-		   }
-	      }
-	  }
-     }
+    for (i = 0; i < n; i++)
+    {
+        if (matrix[i][i] == 0)
+        { /* singular row */
+            for (j = 0; j < i; j++)
+                matrix[j][i] = 0;
+            for (j = i; j < n; j++)
+                matrix[i][j] = 0;
+        }
+        else
+        {
+            for (j = (i + 1); j < n; j++)
+            {
+                temp = matrix[j][i] * matrix[j][j];
+                if (j != i)
+                    matrix[i][j] = temp;
+                for (k = i; k < j; k++)
+                    matrix[i][k] += temp * matrix[j][k];
+            }
+        }
+    }
+}

Copied: pkg/ProbABEL/src/cholesky.cpp (from rev 1027, branches/ProbABEL-refactoring/ProbABEL/src/cholesky.cpp)
===================================================================
--- pkg/ProbABEL/src/cholesky.cpp	                        (rev 0)
+++ pkg/ProbABEL/src/cholesky.cpp	2012-11-26 06:59:20 UTC (rev 1028)
@@ -0,0 +1,152 @@
+/*
+ * cholesky.cpp
+ *
+ *  Created on: Mar 15, 2012
+ *      Author: mkooyman
+ */
+#if EIGEN
+#include "eigen_mematrix.h"
+#include "eigen_mematrix.cpp"
+#else
+#include "mematrix.h"
+#include "mematri1.h"
+#endif
+
+#include <string>
+#include <cstdarg>
+#include <cstdio>
+#include <cstdlib>
+
+/*  SCCS @(#)cholesky2.c    5.2 10/27/98
+ ** subroutine to do Cholesky decompostion on a matrix: C = FDF'
+ **   where F is lower triangular with 1's on the diagonal, and D is diagonal
+ **
+ ** arguments are:
+ **     n         the size of the matrix to be factored
+ **     **matrix  a ragged array containing an n by n submatrix to be factored
+ **     toler     the threshold value for detecting "singularity"
+ **
+ **  The factorization is returned in the lower triangle, D occupies the
+ **    diagonal and the upper triangle is left undisturbed.
+ **    The lower triangle need not be filled in at the start.
+ **
+ **  Return value:  the rank of the matrix (non-negative definite), or -rank
+ **     it not SPD or NND
+ **
+ **  If a column is deemed to be redundant, then that diagonal is set to zero.
+ **
+ **   Terry Therneau
+ */
+
+// modified Yurii Aulchenko 2008-05-20
+int cholesky2_mm(mematrix<double> &matrix, double toler)
+{
+    if (matrix.ncol != matrix.nrow)
+    {
+        fprintf(stderr, "cholesky2_mm: matrix should be square\n");
+        exit(1);
+    }
+    int n = matrix.ncol;
+    double temp;
+    int i, j, k;
+    double eps, pivot;
+    int rank;
+    int nonneg;
+
+    nonneg = 1;
+    eps = 0;
+    for (i = 0; i < n; i++)
+    {
+        if (matrix[i * n + i] > eps)
+            eps = matrix[i * n + i];
+        for (j = (i + 1); j < n; j++)
+            matrix[j * n + i] = matrix[i * n + j];
+    }
+    eps *= toler;
+
+    rank = 0;
+    for (i = 0; i < n; i++)
+    {
+        pivot = matrix[i * n + i];
+        if (pivot < eps)
+        {
+            matrix[i * n + i] = 0;
+            if (pivot < -8 * eps)
+                nonneg = -1;
+        } else
+        {
+            rank++;
+            for (j = (i + 1); j < n; j++)
+            {
+                temp = matrix[j * n + i] / pivot;
+                matrix[j * n + i] = temp;
+                matrix[j * n + j] -= temp * temp * pivot;
+                for (k = (j + 1); k < n; k++)
+                    matrix[k * n + j] -= temp * matrix[k * n + i];
+            }
+        }
+    }
+    return (rank * nonneg);
+}
+
+void chinv2_mm(mematrix<double> &matrix)
+{
+    if (matrix.ncol != matrix.nrow)
+    {
+        fprintf(stderr, "cholesky2_mm: matrix should be square\n");
+        exit(1);
+    }
+
+    int n = matrix.ncol;
+    register double temp;
+    register int i, j, k;
+
+    /*
+     ** invert the cholesky in the lower triangle
+     **   take full advantage of the cholesky's diagonal of 1's
+     */
+    for (i = 0; i < n; i++)
+    {
+        if (matrix[i * n + i] > 0)
+        {
+            matrix[i * n + i] = 1 / matrix[i * n + i]; /*this line inverts D */
+            for (j = (i + 1); j < n; j++)
+            {
+                matrix[j * n + i] = -matrix[j * n + i];
+                for (k = 0; k < i; k++) /*sweep operator */
+                    matrix[j * n + k] += matrix[j * n + i] * matrix[i * n + k];
+            }
+        }
+    }
+
+    /*
+     ** lower triangle now contains inverse of cholesky
+     ** calculate F'DF (inverse of cholesky decomp process) to get inverse
+     **   of original matrix
+     */
+    for (i = 0; i < n; i++)
+    {
+        if (matrix[i * n + i] == 0)
+        { /* singular row */
+            for (j = 0; j < i; j++)
+                matrix[j * n + i] = 0;
+            for (j = i; j < n; j++)
+                matrix[i * n + j] = 0;
+        } else
+        {
+            for (j = (i + 1); j < n; j++)
+            {
+                temp = matrix[j * n + i] * matrix[j * n + j];
+                if (j != i)
+                    matrix[i * n + j] = temp;
+                for (k = i; k < j; k++)
+                    matrix[i * n + k] += temp * matrix[j * n + k];
+            }
+        }
+    }
+    // ugly fix to return only inverse
+    for (int col = 1; col < n; col++)
+        for (int row = 0; row < col; row++)
+            matrix[col * n + row] = matrix[row * n + col];
+}
+

Copied: pkg/ProbABEL/src/cholesky.h (from rev 1027, branches/ProbABEL-refactoring/ProbABEL/src/cholesky.h)
===================================================================
--- pkg/ProbABEL/src/cholesky.h	                        (rev 0)
+++ pkg/ProbABEL/src/cholesky.h	2012-11-26 06:59:20 UTC (rev 1028)
@@ -0,0 +1,21 @@
+/*
+ * cholesky.h
+ *
+ *  Created on: Mar 15, 2012
+ *      Author: mkooyman
+ */
+
+#ifndef CHOLESKY_H_
+#define CHOLESKY_H_
+
+#if EIGEN
+#include "eigen_mematrix.h"
+#include "eigen_mematrix.cpp"
+#else
+#include "mematrix.h"
+#endif
+
+int cholesky2_mm(mematrix<double> &matrix, double toler);
+void chinv2_mm(mematrix<double> &matrix);
+
+#endif /* CHOLESKY_H_ */

Modified: pkg/ProbABEL/src/chsolve2.c
===================================================================
--- pkg/ProbABEL/src/chsolve2.c	2012-11-20 21:57:23 UTC (rev 1027)
+++ pkg/ProbABEL/src/chsolve2.c	2012-11-26 06:59:20 UTC (rev 1028)
@@ -1,42 +1,46 @@
 /*  SCCS @(#)chsolve2.c	5.2 10/27/98
-** Solve the equation Ab = y, where the cholesky decomposition of A and y
-**   are the inputs.
-**
-** Input  **matrix, which contains the chol decomp of an n by n
-**   matrix in its lower triangle.
-**        y[n] contains the right hand side
-**
-**  y is overwriten with b
-**
-**  Terry Therneau
-*/
+ ** Solve the equation Ab = y, where the cholesky decomposition of A and y
+ **   are the inputs.
+ **
+ ** Input  **matrix, which contains the chol decomp of an n by n
+ **   matrix in its lower triangle.
+ **        y[n] contains the right hand side
+ **
+ **  y is overwriten with b
+ **
+ **  Terry Therneau
+ */
 #include "survS.h"
 #include "survproto.h"
 
 void chsolve2(double **matrix, int n, double *y)
-     {
-     register int i,j;
-     register double temp;
+{
+    register int i, j;
+    register double temp;
 
-     /*
+    /*
      ** solve Fb =y
      */
-     for (i=0; i<n; i++) {
-	  temp = y[i] ;
-	  for (j=0; j<i; j++)
-	       temp -= y[j] * matrix[i][j] ;
-	  y[i] = temp ;
-	  }
-     /*
+    for (i = 0; i < n; i++)
+    {
+        temp = y[i];
+        for (j = 0; j < i; j++)
+            temp -= y[j] * matrix[i][j];
+        y[i] = temp;
+    }
+    /*
      ** solve DF'z =b
      */
-     for (i=(n-1); i>=0; i--) {
[TRUNCATED]

To get the complete diff run:
    svnlook diff /svnroot/genabel -r 1028


More information about the Genabel-commits mailing list