[med-svn] [Git][med-team/plink1-9][upstream] New upstream version 1.90~b6.21-201019
Dylan Aïssi
gitlab at salsa.debian.org
Wed Nov 11 09:32:04 GMT 2020
Dylan Aïssi pushed to branch upstream at Debian Med / plink1.9
Commits:
7b42b568 by Dylan Aïssi at 2020-11-11T10:26:23+01:00
New upstream version 1.90~b6.21-201019
- - - - -
11 changed files:
- Makefile
- − Makefile.std
- SFMT.h
- plink.c
- plink_calc.c
- plink_common.h
- plink_family.c
- plink_first_compile
- plink_glm.c
- plink_help.c
- plink_ld.c
Changes:
=====================================
Makefile
=====================================
@@ -29,8 +29,9 @@ endif
BIN ?= plink
CC ?= gcc
CXX ?= g++
-CFLAGS ?= -Wall -O2
-CXXFLAGS ?= -Wall -O2
+CFLAGS ?= -Wall -O2 -g -I../2.0/simde
+CXXFLAGS ?= -Wall -O2 -g -I../2.0/simde
+MACFLAGS ?= -mmacosx-version-min=10.9
PREFIX ?= /usr/local
DESTDIR ?= .
@@ -42,9 +43,11 @@ ifeq ($(SYS), MAC)
ifeq "$(GCC_GTEQ_43)" "1"
CFLAGS ?= -Wall -O2 -flax-vector-conversions
endif
+ CFLAGS += $(MACFLAGS)
+ CXXFLAGS += $(MACFLAGS)
BLASFLAGS ?= -framework Accelerate
LDFLAGS ?= -ldl
- ZLIB ?= ../zlib-1.2.11/libz.1.2.11.dylib
+ ZLIB ?= -L. ../zlib-1.2.11/libz-64.a
endif
ifeq ($(SYS), WIN)
@@ -55,13 +58,13 @@ ifeq ($(SYS), WIN)
# NO_LAPACK.
BLASFLAGS ?= -L. lapack/liblapack.a -L. lapack/librefblas.a
LDFLAGS ?= -lm -static-libgcc
- ZLIB ?= ../zlib-1.2.11/libz.a
+ ZLIB ?= -L. ../zlib-1.2.11/libz.a
endif
# These must appear after the MAC/WIN-specific ?= statements.
BLASFLAGS ?= -L/usr/lib64/atlas -llapack -lblas -lcblas -latlas
LDFLAGS ?= -lm -lpthread -ldl
-ZLIB ?= ../zlib-1.2.11/libz.so.1.2.11
+ZLIB ?= -L. ../zlib-1.2.11/libz.so.1.2.11
ifdef NO_LAPACK
BLASFLAGS=
@@ -75,10 +78,10 @@ OBJS = plink.o plink_assoc.o plink_calc.o plink_cluster.o plink_cnv.o plink_comm
all: $(BIN)
plink: $(OBJS)
- $(CXX) $(OBJS) $(LDFLAGS_EXTRA) $(BLASFLAGS) $(LDFLAGS) -L. $(ZLIB) -o $@
+ $(CXX) $(OBJS) $(LDFLAGS_EXTRA) $(BLASFLAGS) $(LDFLAGS) $(ZLIB) -o $@
plinkw: $(OBJS)
- gfortran -O2 $(OBJS) $(LDFLAGS_EXTRA) -Wl,-Bstatic $(BLASFLAGS) $(LDFLAGS) -L. $(ZLIB) -o $@
+ gfortran -O2 $(OBJS) $(LDFLAGS_EXTRA) -Wl,-Bstatic $(BLASFLAGS) $(LDFLAGS) $(ZLIB) -o $@
install:
$(INSTALL) -d $(DESTDIR)$(PREFIX)/bin
@@ -98,7 +101,7 @@ clean:
# includes a C++ header and exposed functions will need to be declared with
# extern "C".
%.o: %.c
- $(CXX) -c $(CPPFLAGS) $(CFLAGS) $< -o $@
+ $(CXX) -c $(CFLAGS) $< -o $@
%.o: %.cc
- $(CXX) -x c++ -c $(CPPFLAGS) $(CXXFLAGS) $< -o $@
+ $(CXX) -x c++ -c $(CXXFLAGS) $< -o $@
=====================================
Makefile.std deleted
=====================================
@@ -1,104 +0,0 @@
-# General-purpose Makefile for PLINK 1.90
-#
-# Compilation options:
-# Do not link to LAPACK NO_LAPACK
-
-# Leave blank after "=" to disable; put "= 1" to enable
-# (when enabled, "#define NOLAPACK" must be uncommented in plink_common.h)
-NO_LAPACK =
-
-# Variable that allows additional linker flags (e.g., "-L/path/to/libs") to be
-# passed into the linker; to use this, invoke 'make LDFLAGS_EXTRA="<opts>"'.
-LDFLAGS_EXTRA =
-
-.PHONY: clean install
-
-# should autodetect system
-SYS = UNIX
-ifdef SystemRoot
- SYS = WIN
-else
- UNAME := $(shell uname)
- ifeq ($(UNAME), Darwin)
- SYS = MAC
- endif
-endif
-
-# Allow these to be overridden by make arguments or env variables, so people
-# don't have to edit the Makefile to build in a different environment.
-BIN ?= plink
-CC ?= gcc
-CXX ?= g++
-CFLAGS ?= -Wall -O2
-CXXFLAGS ?= -Wall -O2
-
-PREFIX ?= /usr/local
-DESTDIR ?= .
-INSTALL ?= install
-STRIP ?= strip
-
-ifeq ($(SYS), MAC)
- GCC_GTEQ_43 := $(shell expr `g++ -dumpversion | sed -e 's/\.\([0-9][0-9]\)/\1/g' -e 's/\.\([0-9]\)/0\1/g' -e 's/^[0-9]\{3,4\}$$/&00/'` \>= 40300)
- ifeq "$(GCC_GTEQ_43)" "1"
- CFLAGS ?= -Wall -O2 -flax-vector-conversions
- endif
- BLASFLAGS ?= -framework Accelerate
- LDFLAGS ?= -ldl
- ZLIB ?= ../zlib-1.2.11/libz.1.2.11.dylib
-endif
-
-ifeq ($(SYS), WIN)
-# Note that, unlike the Linux and Mac build processes, this STATICALLY links
-# LAPACK, since we have not gotten around to trying dynamically-linked LAPACK
-# on Windows.
-# If you don't already have LAPACK built, you'll probably want to turn on
-# NO_LAPACK.
- BLASFLAGS ?= -L. lapack/liblapack.a -L. lapack/librefblas.a
- LDFLAGS ?= -lm -static-libgcc
- ZLIB ?= ../zlib-1.2.11/libz.a
-endif
-
-# These must appear after the MAC/WIN-specific ?= statements.
-BLASFLAGS ?= -L/usr/lib64/atlas -llapack -lblas -lcblas -latlas
-LDFLAGS ?= -lm -lpthread -ldl
-ZLIB ?= ../zlib-1.2.11/libz.so.1.2.11
-
-ifdef NO_LAPACK
- BLASFLAGS=
-endif
-
-OBJS = plink.o plink_assoc.o plink_calc.o plink_cluster.o plink_cnv.o plink_common.o plink_data.o plink_dosage.o plink_family.o plink_filter.o plink_glm.o plink_help.o plink_homozyg.o plink_lasso.o plink_ld.o plink_matrix.o plink_misc.o plink_perm.o plink_rserve.o plink_set.o plink_stats.o SFMT.o dcdflib.o pigz.o yarn.o Rconnection.o hfile.o bgzf.o
-
-# In the event that you are still concurrently using PLINK 1.07, we suggest
-# renaming that binary to "plink107", and this one to "plink" or "plink1".
-
-all: $(BIN)
-
-plink: $(OBJS)
- $(CXX) $(OBJS) $(LDFLAGS_EXTRA) $(BLASFLAGS) $(LDFLAGS) -L. $(ZLIB) -o $@
-
-plinkw: $(OBJS)
- gfortran -O2 $(OBJS) $(LDFLAGS_EXTRA) -Wl,-Bstatic $(BLASFLAGS) $(LDFLAGS) -L. $(ZLIB) -o $@
-
-install:
- $(INSTALL) -d $(DESTDIR)$(PREFIX)/bin
- $(INSTALL) -c $(BIN) $(DESTDIR)$(PREFIX)/bin
-
-install-strip: install
- $(STRIP) $(DESTDIR)$(PREFIX)/bin/$(BIN)
-
-clean:
- rm -f $(OBJS) plink plinkw
-
-# Pattern-based rules for compiling object (.o) files; basically identical to
-# GNU make's built-in rules, except we explicitly use "g++" instead of $(CC).
-
-# Compiling C files with a C++ compiler is deprecated, but the code needs
-# to be cleaned up before we can switch to cc. E.g. plink_rserve.c
-# includes a C++ header and exposed functions will need to be declared with
-# extern "C".
-%.o: %.c
- $(CXX) -c $(CPPFLAGS) $(CFLAGS) $< -o $@
-
-%.o: %.cc
- $(CXX) -x c++ -c $(CPPFLAGS) $(CXXFLAGS) $< -o $@
=====================================
SFMT.h
=====================================
@@ -129,7 +129,12 @@ extern "C" {
128-bit SIMD like data type for standard C
------------------------------------------*/
#ifdef __LP64__
- #include <emmintrin.h>
+ #ifdef __SSE2__
+ #include <emmintrin.h>
+ #else
+ #define SIMDE_ENABLE_NATIVE_ALIASES
+ #include "x86/sse2.h"
+ #endif
/** 128-bit data structure */
union W128_T {
=====================================
plink.c
=====================================
@@ -93,7 +93,7 @@
static const char ver_str[] =
#ifdef STABLE_BUILD
- "PLINK v1.90b6.18"
+ "PLINK v1.90b6.22"
#else
"PLINK v1.90p"
#endif
@@ -105,10 +105,10 @@ static const char ver_str[] =
#else
" 32-bit"
#endif
- " (16 Jun 2020)";
+ " (3 Nov 2020)";
static const char ver_str2[] =
// include leading space if day < 10, so character length stays the same
- ""
+ " "
#ifdef STABLE_BUILD
"" // (don't want this when version number has two trailing digits)
#else
@@ -3327,6 +3327,9 @@ int32_t main(int32_t argc, char** argv) {
uint32_t lasso_lambda_iters = 0;
uint32_t testmiss_modifier = 0;
uint32_t testmiss_mperm_val = 0;
+#ifdef USE_MKL
+ uint32_t mkl_native = 0;
+#endif
// this default limit plays well with e.g. fbstring small-string optimization
uint32_t new_id_max_allele_len = 23;
@@ -6158,14 +6161,12 @@ int32_t main(int32_t argc, char** argv) {
}
}
} else if (!memcmp(argptr2, "istance-matrix", 15)) {
- if (distance_exp != 0.0) {
- logerrprint("Error: --distance-matrix cannot be used with --distance-exp.\n");
+ if (calculation_type & CALC_DISTANCE) {
+ // not worth the maintenance burden of making a few of these subcases
+ // work
+ logerrprint("Error: --distance-matrix cannot be used with --distance.\n");
goto main_ret_INVALID_CMDLINE;
- }
- if (dist_calc_type & DISTANCE_1_MINUS_IBS) {
- logerrprint("Error: --distance-matrix flag cannot be used with \"--distance 1-ibs\".\n");
- goto main_ret_INVALID_CMDLINE_A;
- }
+ }
calculation_type |= CALC_PLINK1_DISTANCE_MATRIX;
goto main_param_zero;
} else if (!memcmp(argptr2, "ummy", 5)) {
@@ -7681,6 +7682,10 @@ int32_t main(int32_t argc, char** argv) {
glm_modifier |= GLM_INTERACTION;
goto main_param_zero;
} else if (!memcmp(argptr2, "bs-matrix", 10)) {
+ if (calculation_type & CALC_DISTANCE) {
+ logerrprint("Error: --ibs-matrix cannot be used with --distance.\n");
+ goto main_ret_INVALID_CMDLINE;
+ }
calculation_type |= CALC_PLINK1_IBS_MATRIX;
goto main_param_zero;
} else if (!memcmp(argptr2, "d-delim", 8)) {
@@ -9716,6 +9721,11 @@ int32_t main(int32_t argc, char** argv) {
} else if (!memcmp(argptr2, "oweb", 5)) {
logprint("Note: --noweb has no effect since no web check is implemented yet.\n");
goto main_param_zero;
+ } else if (!memcmp(argptr2, "ative", 6)) {
+#ifdef USE_MKL
+ mkl_native = 1;
+#endif
+ goto main_param_zero;
} else {
goto main_ret_INVALID_CMDLINE_UNRECOGNIZED;
}
@@ -13304,6 +13314,12 @@ int32_t main(int32_t argc, char** argv) {
goto main_ret_INVALID_CMDLINE_A;
}
+#ifdef USE_MKL
+ if (!mkl_native) {
+ mkl_cbwr_set(MKL_CBWR_COMPATIBLE);
+ }
+#endif
+
free_cond(subst_argv);
free_cond(script_buf);
free_cond(rerun_buf);
=====================================
plink_calc.c
=====================================
@@ -2689,6 +2689,8 @@ int32_t unrelated_herit_batch(uint32_t load_grm_bin, char* grmname, char* phenon
#endif
int32_t ibs_test_calc(pthread_t* threads, char* read_dists_fname, uintptr_t unfiltered_sample_ct, uintptr_t* sample_exclude, uintptr_t sample_ct, uintptr_t perm_ct, uintptr_t pheno_nm_ct, uintptr_t pheno_ctrl_ct, uintptr_t* pheno_nm, uintptr_t* pheno_c) {
+ // g_dists and g_half_marker_ct_recip assumed to be populated by
+ // calc_distance().
unsigned char* bigstack_mark = g_bigstack_base;
uintptr_t unfiltered_sample_ctl = BITCT_TO_WORDCT(unfiltered_sample_ct);
uintptr_t pheno_nm_ctl = BITCT_TO_WORDCT(pheno_nm_ct);
@@ -3751,6 +3753,7 @@ uint32_t distance_d_write_1mibs_sq_emitn(uint32_t overflow_ct, unsigned char* re
int32_t distance_d_write(FILE** outfile_ptr, FILE** outfile2_ptr, FILE** outfile3_ptr, int32_t dist_calc_type, char* outname, char* outname_end, double* dists, double half_marker_ct_recip, uint32_t sample_ct, int32_t first_sample_idx, int32_t end_sample_idx, int32_t parallel_idx, int32_t parallel_tot, unsigned char* membuf) {
// membuf assumed to be of at least size sample_ct * 8.
+ printf("distance_d_write\n");
uint32_t bin4 = dist_calc_type & DISTANCE_BIN4;
int32_t shape = dist_calc_type & DISTANCE_SHAPEMASK;
int32_t write_alcts = dist_calc_type & DISTANCE_ALCT;
@@ -7573,7 +7576,8 @@ int32_t calc_distance(pthread_t* threads, uint32_t parallel_idx, uint32_t parall
FILE* outfile3 = nullptr;
uintptr_t unfiltered_sample_ct4 = (unfiltered_sample_ct + 3) / 4;
uint64_t dists_alloc = 0;
- uint32_t missing_wt_needed = ((calculation_type & CALC_DISTANCE) || ((!read_dists_fname) && (calculation_type & (CALC_IBS_TEST | CALC_GROUPDIST | CALC_REGRESS_DISTANCE)))) && (!(dist_calc_type & DISTANCE_FLAT_MISSING));
+ uint32_t unweighted_and_flat = ((calculation_type & (CALC_PLINK1_DISTANCE_MATRIX | CALC_PLINK1_IBS_MATRIX)) || (dist_calc_type & DISTANCE_FLAT_MISSING)) && (!distance_wts_fname) && (distance_exp == 0.0);
+ uint32_t missing_wt_needed = ((calculation_type & CALC_DISTANCE) || ((!read_dists_fname) && (calculation_type & (CALC_IBS_TEST | CALC_GROUPDIST | CALC_REGRESS_DISTANCE)))) && (!unweighted_and_flat);
uint32_t unwt_needed = 0;
uint32_t marker_weight_sum = 0;
int32_t retval = 0;
@@ -7986,7 +7990,8 @@ int32_t calc_distance(pthread_t* threads, uint32_t parallel_idx, uint32_t parall
giptr3++;
}
}
- fill_subset_weights(subset_weights, &(main_weights[ukk]));
+ fill_subset_weights(subset_weights, &(main_weights[marker_idx - ujj + ukk]));
+ g_subset_weights_i = &(wtbuf[ukk]);
uii = is_last_block && (ukk + (MULTIPLEX_DIST_EXP / 3) >= ujj);
if (spawn_threads2(threads, &calc_wdist_thread, dist_thread_ct, uii)) {
goto calc_distance_ret_THREAD_CREATE_FAIL;
@@ -8001,6 +8006,20 @@ int32_t calc_distance(pthread_t* threads, uint32_t parallel_idx, uint32_t parall
} while (!is_last_block);
putc_unlocked('\r', stdout);
logprint("Distance matrix calculation complete.\n");
+ if (calculation_type & (CALC_DISTANCE | CALC_IBS_TEST)) {
+ // Calculate g_half_marker_ct_recip here since main_weights may be
+ // overwritten below.
+ if (main_weights) {
+ dyy = 0.0;
+ marker_uidx = 0;
+ for (marker_idx = 0; marker_idx < marker_ct; marker_idx++) {
+ dyy += main_weights[marker_idx];
+ }
+ g_half_marker_ct_recip = 0.5 / dyy;
+ } else {
+ g_half_marker_ct_recip = 0.5 / (double)marker_ct;
+ }
+ }
bigstack_reset(masks);
if (calculation_type & (CALC_PLINK1_DISTANCE_MATRIX | CALC_PLINK1_IBS_MATRIX)) {
if (bigstack_alloc_c(16 * sample_ct, &writebuf)) {
@@ -8159,35 +8178,17 @@ int32_t calc_distance(pthread_t* threads, uint32_t parallel_idx, uint32_t parall
}
}
- if (calculation_type & (CALC_DISTANCE | CALC_IBS_TEST)) {
- if ((distance_exp == 0.0) || (!(dist_calc_type & (DISTANCE_IBS | DISTANCE_1_MINUS_IBS)))) {
- g_half_marker_ct_recip = 0.5 / (double)marker_ct;
- } else {
- dyy = 0.0;
- marker_uidx = 0;
- for (marker_idx = 0; marker_idx < marker_ct; marker_uidx++, marker_idx++) {
- next_unset_ul_unsafe_ck(marker_exclude, &marker_uidx);
- dxx = set_allele_freqs[marker_uidx];
- if ((dxx > 0.0) && (dxx < 1.0)) {
- dyy += pow(2 * dxx * (1.0 - dxx), -distance_exp);
- } else {
- dyy += 1.0;
- }
- }
- g_half_marker_ct_recip = 0.5 / dyy;
- }
- if (calculation_type & CALC_DISTANCE) {
- if (!parallel_idx) {
- retval = distance_d_write_ids(outname, outname_end, dist_calc_type, unfiltered_sample_ct, sample_exclude, sample_ids, max_sample_id_len);
- if (retval) {
- goto calc_distance_ret_1;
- }
- LOGPRINTFWW("IDs written to %s .\n", outname);
- }
- retval = distance_d_write(&outfile, &outfile2, &outfile3, dist_calc_type, outname, outname_end, g_dists, g_half_marker_ct_recip, sample_ct, g_thread_start[0], g_thread_start[dist_thread_ct], parallel_idx, parallel_tot, (unsigned char*)geno);
+ if (calculation_type & CALC_DISTANCE) {
+ if (!parallel_idx) {
+ retval = distance_d_write_ids(outname, outname_end, dist_calc_type, unfiltered_sample_ct, sample_exclude, sample_ids, max_sample_id_len);
if (retval) {
- goto calc_distance_ret_1;
+ goto calc_distance_ret_1;
}
+ LOGPRINTFWW("IDs written to %s .\n", outname);
+ }
+ retval = distance_d_write(&outfile, &outfile2, &outfile3, dist_calc_type, outname, outname_end, g_dists, g_half_marker_ct_recip, sample_ct, g_thread_start[0], g_thread_start[dist_thread_ct], parallel_idx, parallel_tot, (unsigned char*)geno);
+ if (retval) {
+ goto calc_distance_ret_1;
}
}
bigstack_reset(bigstack_mark);
=====================================
plink_common.h
=====================================
@@ -187,14 +187,12 @@
// http://esr.ibiblio.org/?p=5095 ).
#ifdef __LP64__
- #ifndef __SSE2__
- // It's obviously possible to support this by writing 64-bit non-SSE2 code
- // shadowing each SSE2 intrinsic, but this almost certainly isn't worth the
- // development/testing effort until regular PLINK 2.0 development is
- // complete. No researcher has ever asked me for this feature.
- #error "64-bit builds currently require SSE2. Try producing a 32-bit build instead."
+ #ifdef __SSE2__
+ #include <emmintrin.h>
+ #else
+ #define SIMDE_ENABLE_NATIVE_ALIASES
+ #include "x86/sse2.h"
#endif
- #include <emmintrin.h>
#define VECFTYPE __m128
#define VECITYPE __m128i
=====================================
plink_family.c
=====================================
@@ -3894,7 +3894,7 @@ int32_t dfam(pthread_t* threads, FILE* bedfile, uintptr_t bed_offset, char* outn
uint32_t dfam_sample_ct;
uint32_t dfam_sample_ctl;
uint32_t dfam_sample_ctl2;
- uint32_t dfam_sample_ctv;
+ uint32_t dfam_sample_ctaw;
uint32_t chrom_fo_idx;
uint32_t chrom_end;
uint32_t chrom_idx;
@@ -4165,7 +4165,7 @@ int32_t dfam(pthread_t* threads, FILE* bedfile, uintptr_t bed_offset, char* outn
dfam_sample_ct = unfiltered_sample_ct - popcount_longs(dfam_sample_exclude, unfiltered_sample_ctl);
dfam_sample_ctl = BITCT_TO_WORDCT(dfam_sample_ct);
dfam_sample_ctl2 = QUATERCT_TO_WORDCT(dfam_sample_ct);
- dfam_sample_ctv = BITCT_TO_ALIGNED_WORDCT(dfam_sample_ct);
+ dfam_sample_ctaw = BITCT_TO_ALIGNED_WORDCT(dfam_sample_ct);
if (bigstack_alloc_ui(unfiltered_sample_ct, &sample_uidx_to_idx)) {
goto dfam_ret_NOMEM;
}
@@ -4384,7 +4384,7 @@ int32_t dfam(pthread_t* threads, FILE* bedfile, uintptr_t bed_offset, char* outn
perm_vec_ctcl8m = round_up_pow2(g_perm_vec_ct, CACHELINE_DBL);
if (bigstack_alloc_ul(dfam_sample_ct * perm_vec_wcta, &g_dfam_perm_vecst) ||
- bigstack_alloc_ul(g_perm_vec_ct * dfam_sample_ctv, &g_perm_vecs)) {
+ bigstack_alloc_ul(g_perm_vec_ct * dfam_sample_ctaw, &g_perm_vecs)) {
goto dfam_ret_NOMEM;
}
// initialize phenotype permutations.
@@ -5728,7 +5728,8 @@ int32_t qfam(pthread_t* threads, FILE* bedfile, uintptr_t bed_offset, char* outn
}
for (block_idx = 0, marker_idx = marker_idx_base; block_idx < block_size; marker_uidx++, marker_idx++) {
if (IS_SET(marker_exclude, marker_uidx)) {
- marker_uidx = next_set_ul_unsafe(marker_exclude, marker_uidx);
+ // bugfix (16 Sep 2020): set -> unset
+ marker_uidx = next_unset_ul_unsafe(marker_exclude, marker_uidx);
seek_flag = 1;
}
if (marker_uidx >= chrom_end) {
=====================================
plink_first_compile
=====================================
@@ -25,4 +25,4 @@ make
# Compile
cd ../1.9
-make -f Makefile.std plink
+make
=====================================
plink_glm.c
=====================================
@@ -1702,7 +1702,7 @@ uint32_t glm_logistic(uintptr_t cur_batch_size, uintptr_t param_ct, uintptr_t sa
uintptr_t joint_test_requested = (constraints_con_major? 1 : 0);
uintptr_t param_ctx = param_ct + joint_test_requested;
uintptr_t param_ctx_msi = param_ctx - skip_intercept;
- uintptr_t sample_validx_ctv = BITCT_TO_ALIGNED_WORDCT(sample_valid_ct + missing_ct);
+ uintptr_t sample_validx_ctaw = BITCT_TO_ALIGNED_WORDCT(sample_valid_ct + missing_ct);
uintptr_t perm_fail_ct = 0;
uintptr_t cur_word = 0;
uintptr_t perm_idx;
@@ -1811,7 +1811,7 @@ uint32_t glm_logistic(uintptr_t cur_batch_size, uintptr_t param_ct, uintptr_t sa
}
}
coef = &(coef[param_cta4]);
- perm_vecs = &(perm_vecs[sample_validx_ctv]);
+ perm_vecs = &(perm_vecs[sample_validx_ctaw]);
}
return perm_fail_ct;
}
@@ -3343,6 +3343,7 @@ THREAD_RET_TYPE glm_logistic_maxt_thread(void* arg) {
dxx = (double)coef[pidx * cur_param_cta4 + 1];
dxx *= dxx;
dxx /= (double)regression_results[pidx * param_ctx_m1];
+ printf("coef: %g denom: %g stat: %g\n", coef[pidx * cur_param_cta4 + 1], regression_results[pidx * param_ctx_m1], dxx);
if (dxx > stat_high) {
success_2incr += 2;
} else if (dxx > stat_low) {
@@ -5972,7 +5973,7 @@ int32_t glm_logistic_assoc(pthread_t* threads, FILE* bedfile, uintptr_t bed_offs
uintptr_t param_idx_end;
uintptr_t sample_valid_ct;
uintptr_t sample_valid_cta4;
- uintptr_t sample_valid_ctv;
+ uintptr_t sample_valid_ctaw;
uintptr_t sample_valid_ctv2;
uintptr_t sample_idx;
uintptr_t param_ctx_max;
@@ -6062,7 +6063,7 @@ int32_t glm_logistic_assoc(pthread_t* threads, FILE* bedfile, uintptr_t bed_offs
goto glm_logistic_assoc_ret_1;
}
sample_valid_cta4 = round_up_pow2(sample_valid_ct, 4);
- sample_valid_ctv = BITCT_TO_ALIGNED_WORDCT(sample_valid_ct);
+ sample_valid_ctaw = BITCT_TO_ALIGNED_WORDCT(sample_valid_ct);
sample_valid_ctv2 = QUATERCT_TO_ALIGNED_WORDCT(sample_valid_ct);
final_mask = get_final_mask(sample_valid_ct);
param_ct_maxa4 = round_up_pow2(param_ct_max, 4);
@@ -6403,16 +6404,18 @@ int32_t glm_logistic_assoc(pthread_t* threads, FILE* bedfile, uintptr_t bed_offs
}
}
- if (bigstack_alloc_ul(sample_valid_ctv, &pheno_c_collapsed)) {
+ if (bigstack_alloc_ul(sample_valid_ctaw, &pheno_c_collapsed)) {
goto glm_logistic_assoc_ret_NOMEM;
}
+ // bugfix (21 Sep 2020): trailing word not guaranteed to be zero
+ pheno_c_collapsed[sample_valid_ctaw - 1] = 0;
copy_bitarr_subset(pheno_c, load_mask, unfiltered_sample_ct, sample_valid_ct, pheno_c_collapsed);
- g_perm_case_ct = popcount_longs(pheno_c_collapsed, sample_valid_ctv);
+ g_perm_case_ct = popcount_longs(pheno_c_collapsed, sample_valid_ctaw);
if ((!g_perm_case_ct) || (g_perm_case_ct == sample_valid_ct)) {
goto glm_logistic_assoc_ret_PHENO_CONSTANT;
}
if (do_perms) {
- if (bigstack_alloc_ul(orig_perm_batch_size * sample_valid_ctv, &g_perm_vecs)) {
+ if (bigstack_alloc_ul(orig_perm_batch_size * sample_valid_ctaw, &g_perm_vecs)) {
goto glm_logistic_assoc_ret_NOMEM;
}
g_perm_is_1bit = 1;
@@ -7980,7 +7983,7 @@ int32_t glm_logistic_nosnp(pthread_t* threads, FILE* bedfile, uintptr_t bed_offs
double* dptr;
uintptr_t sample_valid_ct;
uintptr_t sample_valid_cta4;
- uintptr_t sample_valid_ctv;
+ uintptr_t sample_valid_ctaw;
uintptr_t sample_valid_ctv2;
uintptr_t sample_uidx_stop;
uintptr_t sample_idx;
@@ -8074,7 +8077,7 @@ int32_t glm_logistic_nosnp(pthread_t* threads, FILE* bedfile, uintptr_t bed_offs
}
}
sample_valid_cta4 = round_up_pow2(sample_valid_ct, 4);
- sample_valid_ctv = BITCT_TO_ALIGNED_WORDCT(sample_valid_ct);
+ sample_valid_ctaw = BITCT_TO_ALIGNED_WORDCT(sample_valid_ct);
sample_valid_ctv2 = QUATERCT_TO_ALIGNED_WORDCT(sample_valid_ct);
if (condition_mname || condition_fname) {
@@ -8318,11 +8321,12 @@ int32_t glm_logistic_nosnp(pthread_t* threads, FILE* bedfile, uintptr_t bed_offs
bigstack_alloc_f(param_ct * param_cta4, ¶m_2d_buf2) ||
bigstack_alloc_f(perm_batch_size * (param_ctx - uii), ®ression_results) ||
bigstack_alloc_ul(BITCT_TO_WORDCT(perm_batch_size), &perm_fails) ||
- bigstack_alloc_ul(perm_batch_size * sample_valid_ctv, &g_perm_vecs)) {
+ bigstack_alloc_ul(perm_batch_size * sample_valid_ctaw, &g_perm_vecs)) {
goto glm_logistic_nosnp_ret_NOMEM;
}
+ g_perm_vecs[sample_valid_ctaw - 1] = 0;
copy_bitarr_subset(pheno_c, load_mask, unfiltered_sample_ct, sample_valid_ct, g_perm_vecs);
- g_perm_case_ct = popcount_longs(g_perm_vecs, sample_valid_ctv);
+ g_perm_case_ct = popcount_longs(g_perm_vecs, sample_valid_ctaw);
if ((!g_perm_case_ct) || (g_perm_case_ct == sample_valid_ct)) {
goto glm_logistic_nosnp_ret_PHENO_CONSTANT;
}
@@ -8813,7 +8817,7 @@ uint32_t glm_logistic_dosage(uintptr_t sample_ct, uintptr_t* cur_samples, uintpt
return 0;
}
uintptr_t sample_valid_cta4 = round_up_pow2(sample_valid_ct, 4);
- uintptr_t sample_valid_ctv = BITCT_TO_ALIGNED_WORDCT(sample_valid_ct);
+ uintptr_t sample_valid_ctaw = BITCT_TO_ALIGNED_WORDCT(sample_valid_ct);
uintptr_t param_cta4 = round_up_pow2(param_ct, 4);
float* fptr = covars_cov_major;
uintptr_t case_ct;
@@ -8822,8 +8826,9 @@ uint32_t glm_logistic_dosage(uintptr_t sample_ct, uintptr_t* cur_samples, uintpt
uintptr_t covar_idx;
double dxx;
double dyy;
+ perm_vec[sample_valid_ctaw - 1] = 0;
copy_bitarr_subset(pheno_c, cur_samples, sample_ct, sample_valid_ct, perm_vec);
- case_ct = popcount_longs(perm_vec, sample_valid_ctv);
+ case_ct = popcount_longs(perm_vec, sample_valid_ctaw);
if ((!case_ct) || (case_ct == sample_valid_ct)) {
return 0;
}
=====================================
plink_help.c
=====================================
@@ -2153,6 +2153,9 @@ int32_t disp_help(uint32_t param_ct, char** argv) {
" --perm-batch-size <val> : Set number of permutations per batch for some\n"
" permutation tests.\n"
);
+ help_print("native", &help_ctrl, 0,
+" --native : Allow Intel MKL to use processor-dependent code paths.\n"
+ );
help_print("output-min-p", &help_ctrl, 0,
" --output-min-p <p> : Specify minimum p-value to write to reports.\n"
);
=====================================
plink_ld.c
=====================================
@@ -762,9 +762,9 @@ int32_t ld_prune(Ld_info* ldip, FILE* bedfile, uintptr_t bed_offset, uintptr_t m
uintptr_t founder_ct = popcount_longs(founder_info, unfiltered_sample_ctv2 / 2);
uintptr_t founder_ctl = BITCT_TO_WORDCT(founder_ct);
#ifdef __LP64__
- uintptr_t founder_ctv = BITCT_TO_ALIGNED_WORDCT(founder_ct);
+ uintptr_t founder_ctaw = BITCT_TO_ALIGNED_WORDCT(founder_ct);
#else
- uintptr_t founder_ctv = founder_ctl;
+ uintptr_t founder_ctaw = founder_ctl;
#endif
uintptr_t founder_ct_mld = (founder_ct + MULTIPLEX_LD - 1) / MULTIPLEX_LD;
uint32_t founder_ct_mld_m1 = ((uint32_t)founder_ct_mld) - 1;
@@ -921,7 +921,7 @@ int32_t ld_prune(Ld_info* ldip, FILE* bedfile, uintptr_t bed_offset, uintptr_t m
bigstack_alloc_ul(unfiltered_sample_ctv2, &loadbuf) ||
bigstack_alloc_ul(window_max * founder_ct_192_long, &geno) ||
bigstack_alloc_ul(window_max * founder_ct_192_long, &geno_masks) ||
- bigstack_alloc_ul(window_max * founder_ctv, &geno_mmasks) ||
+ bigstack_alloc_ul(window_max * founder_ctaw, &geno_mmasks) ||
bigstack_alloc_ui(window_max, &missing_cts) ||
bigstack_alloc_d(window_max, &sums) ||
bigstack_alloc_d(window_max, &variance_recips)) {
@@ -986,7 +986,7 @@ int32_t ld_prune(Ld_info* ldip, FILE* bedfile, uintptr_t bed_offset, uintptr_t m
if (is_haploid && hh_exists) {
haploid_fix(hh_exists, founder_include2, founder_male_include2, founder_ct, is_x, is_y, (unsigned char*)(&(geno[ulii * founder_ct_192_long])));
}
- if (!ld_process_load(&(geno[ulii * founder_ct_192_long]), &(geno_masks[ulii * founder_ct_192_long]), &(geno_mmasks[ulii * founder_ctv]), &(missing_cts[ulii]), &(sums[ulii]), &(variance_recips[ulii]), founder_ct, is_x && (!ignore_x), weighted_x, nonmale_founder_ct, founder_male_include2, nonmale_geno, nonmale_masks, ulii * founder_ct_192_long)) {
+ if (!ld_process_load(&(geno[ulii * founder_ct_192_long]), &(geno_masks[ulii * founder_ct_192_long]), &(geno_mmasks[ulii * founder_ctaw]), &(missing_cts[ulii]), &(sums[ulii]), &(variance_recips[ulii]), founder_ct, is_x && (!ignore_x), weighted_x, nonmale_founder_ct, founder_male_include2, nonmale_geno, nonmale_masks, ulii * founder_ct_192_long)) {
SET_BIT(uii, pruned_arr);
cur_exclude_ct++;
}
@@ -1037,7 +1037,7 @@ int32_t ld_prune(Ld_info* ldip, FILE* bedfile, uintptr_t bed_offset, uintptr_t m
} else {
non_missing_ct = fixed_non_missing_ct - missing_cts[ujj];
if (fixed_missing_ct && missing_cts[ujj]) {
- non_missing_ct += popcount_longs_intersect(&(geno_mmasks[uii * founder_ctv]), &(geno_mmasks[ujj * founder_ctv]), founder_ctl);
+ non_missing_ct += popcount_longs_intersect(&(geno_mmasks[uii * founder_ctaw]), &(geno_mmasks[ujj * founder_ctaw]), founder_ctl);
}
}
non_missing_ctd = (double)((int32_t)non_missing_ct);
@@ -1255,7 +1255,7 @@ int32_t ld_prune(Ld_info* ldip, FILE* bedfile, uintptr_t bed_offset, uintptr_t m
memcpy(&(nonmale_geno[uii * founder_ct_192_long]), &(nonmale_geno[ujj * founder_ct_192_long]), founder_ct_192_long * sizeof(intptr_t));
memcpy(&(nonmale_masks[uii * founder_ct_192_long]), &(nonmale_masks[ujj * founder_ct_192_long]), founder_ct_192_long * sizeof(intptr_t));
}
- memcpy(&(geno_mmasks[uii * founder_ctv]), &(geno_mmasks[ujj * founder_ctv]), founder_ctl * sizeof(intptr_t));
+ memcpy(&(geno_mmasks[uii * founder_ctaw]), &(geno_mmasks[ujj * founder_ctaw]), founder_ctl * sizeof(intptr_t));
live_indices[uii] = live_indices[ujj];
start_arr[uii] = start_arr[ujj];
missing_cts[uii] = missing_cts[ujj];
@@ -1301,7 +1301,7 @@ int32_t ld_prune(Ld_info* ldip, FILE* bedfile, uintptr_t bed_offset, uintptr_t m
if (is_haploid && hh_exists) {
haploid_fix(hh_exists, founder_include2, founder_male_include2, founder_ct, is_x, is_y, (unsigned char*)(&(geno[cur_window_size * founder_ct_192_long])));
}
- if (!ld_process_load(&(geno[cur_window_size * founder_ct_192_long]), &(geno_masks[cur_window_size * founder_ct_192_long]), &(geno_mmasks[cur_window_size * founder_ctv]), &(missing_cts[cur_window_size]), &(sums[cur_window_size]), &(variance_recips[cur_window_size]), founder_ct, is_x && (!ignore_x), weighted_x, nonmale_founder_ct, founder_male_include2, nonmale_geno, nonmale_masks, cur_window_size * founder_ct_192_long)) {
+ if (!ld_process_load(&(geno[cur_window_size * founder_ct_192_long]), &(geno_masks[cur_window_size * founder_ct_192_long]), &(geno_mmasks[cur_window_size * founder_ctaw]), &(missing_cts[cur_window_size]), &(sums[cur_window_size]), &(variance_recips[cur_window_size]), founder_ct, is_x && (!ignore_x), weighted_x, nonmale_founder_ct, founder_male_include2, nonmale_geno, nonmale_masks, cur_window_size * founder_ct_192_long)) {
SET_BIT(window_unfiltered_end, pruned_arr);
cur_exclude_ct++;
}
@@ -3742,10 +3742,10 @@ THREAD_RET_TYPE fast_epi_thread(void* arg) {
uintptr_t marker_ct = g_epi_marker_ct;
uint32_t case_ct = g_epi_case_ct;
uint32_t ctrl_ct = g_epi_ctrl_ct;
- uint32_t case_ctv3 = BITCT_TO_ALIGNED_WORDCT(case_ct);
- uint32_t ctrl_ctv3 = BITCT_TO_ALIGNED_WORDCT(ctrl_ct);
- uint32_t case_ctsplit = 3 * case_ctv3;
- uint32_t ctrl_ctsplit = 3 * ctrl_ctv3;
+ uint32_t case_ctaw3 = BITCT_TO_ALIGNED_WORDCT(case_ct);
+ uint32_t ctrl_ctaw3 = BITCT_TO_ALIGNED_WORDCT(ctrl_ct);
+ uint32_t case_ctsplit = 3 * case_ctaw3;
+ uint32_t ctrl_ctsplit = 3 * ctrl_ctaw3;
uint32_t tot_ctsplit = case_ctsplit + ctrl_ctsplit;
uint32_t is_case_only = (g_epi_flag / EPI_FAST_CASE_ONLY) & 1;
uint32_t group_ct = 2 - is_case_only;
@@ -3864,14 +3864,14 @@ THREAD_RET_TYPE fast_epi_thread(void* arg) {
nm_case_fixed = is_set_ul(zmiss1, block_idx1 * 2);
nm_ctrl_fixed = is_set_ul(zmiss1, block_idx1 * 2 + 1);
nm_fixed = nm_case_fixed & nm_ctrl_fixed;
- tot1[0] = popcount_longs(cur_geno1, case_ctv3);
- tot1[1] = popcount_longs(&(cur_geno1[case_ctv3]), case_ctv3);
- tot1[2] = popcount_longs(&(cur_geno1[2 * case_ctv3]), case_ctv3);
+ tot1[0] = popcount_longs(cur_geno1, case_ctaw3);
+ tot1[1] = popcount_longs(&(cur_geno1[case_ctaw3]), case_ctaw3);
+ tot1[2] = popcount_longs(&(cur_geno1[2 * case_ctaw3]), case_ctaw3);
if (!is_case_only) {
cur_geno1_ctrls = &(cur_geno1[case_ctsplit]);
- tot1[3] = popcount_longs(cur_geno1_ctrls, ctrl_ctv3);
- tot1[4] = popcount_longs(&(cur_geno1_ctrls[ctrl_ctv3]), ctrl_ctv3);
- tot1[5] = popcount_longs(&(cur_geno1_ctrls[2 * ctrl_ctv3]), ctrl_ctv3);
+ tot1[3] = popcount_longs(cur_geno1_ctrls, ctrl_ctaw3);
+ tot1[4] = popcount_longs(&(cur_geno1_ctrls[ctrl_ctaw3]), ctrl_ctaw3);
+ tot1[5] = popcount_longs(&(cur_geno1_ctrls[2 * ctrl_ctaw3]), ctrl_ctaw3);
if (is_boost) {
if (nm_fixed) {
cur_boost_precalc2 = &(boost_precalc2[block_idx2 * 6]);
@@ -3894,7 +3894,7 @@ THREAD_RET_TYPE fast_epi_thread(void* arg) {
cur_zmiss2 = (zmiss2[block_idx2 / BITCT2] >> (2 * (block_idx2 % BITCT2))) & 3;
cur_zmiss2_tmp = cur_zmiss2 & 1;
if (nm_case_fixed) {
- two_locus_count_table_zmiss1(cur_geno1, cur_geno2, counts, case_ctv3, cur_zmiss2_tmp);
+ two_locus_count_table_zmiss1(cur_geno1, cur_geno2, counts, case_ctaw3, cur_zmiss2_tmp);
if (cur_zmiss2_tmp) {
counts[2] = tot1[0] - counts[0] - counts[1];
counts[5] = tot1[1] - counts[3] - counts[4];
@@ -3903,7 +3903,7 @@ THREAD_RET_TYPE fast_epi_thread(void* arg) {
counts[7] = cur_tot2[1] - counts[1] - counts[4];
counts[8] = cur_tot2[2] - counts[2] - counts[5];
} else {
- two_locus_count_table(cur_geno1, cur_geno2, counts, case_ctv3, cur_zmiss2_tmp);
+ two_locus_count_table(cur_geno1, cur_geno2, counts, case_ctaw3, cur_zmiss2_tmp);
if (cur_zmiss2_tmp) {
counts[2] = tot1[0] - counts[0] - counts[1];
counts[5] = tot1[1] - counts[3] - counts[4];
@@ -3913,7 +3913,7 @@ THREAD_RET_TYPE fast_epi_thread(void* arg) {
if (!is_case_only) {
cur_zmiss2_tmp = cur_zmiss2 >> 1;
if (nm_ctrl_fixed) {
- two_locus_count_table_zmiss1(cur_geno1_ctrls, &(cur_geno2[case_ctsplit]), &(counts[9]), ctrl_ctv3, cur_zmiss2_tmp);
+ two_locus_count_table_zmiss1(cur_geno1_ctrls, &(cur_geno2[case_ctsplit]), &(counts[9]), ctrl_ctaw3, cur_zmiss2_tmp);
if (cur_zmiss2_tmp) {
counts[11] = tot1[3] - counts[9] - counts[10];
counts[14] = tot1[4] - counts[12] - counts[13];
@@ -3922,7 +3922,7 @@ THREAD_RET_TYPE fast_epi_thread(void* arg) {
counts[16] = cur_tot2[4] - counts[10] - counts[13];
counts[17] = cur_tot2[5] - counts[11] - counts[14];
} else {
- two_locus_count_table(cur_geno1_ctrls, &(cur_geno2[case_ctsplit]), &(counts[9]), ctrl_ctv3, cur_zmiss2_tmp);
+ two_locus_count_table(cur_geno1_ctrls, &(cur_geno2[case_ctsplit]), &(counts[9]), ctrl_ctaw3, cur_zmiss2_tmp);
if (cur_zmiss2_tmp) {
counts[11] = tot1[3] - counts[9] - counts[10];
counts[14] = tot1[4] - counts[12] - counts[13];
@@ -5073,8 +5073,8 @@ THREAD_RET_TYPE ld_dprime_thread(void* arg) {
uintptr_t block_idx1_end = ((tidx + 1) * g_ld_idx1_block_size) / g_ld_thread_ct;
uintptr_t marker_idx2_maxw = g_ld_marker_ctm8;
uintptr_t founder_ct = g_ld_founder_ct;
- uint32_t founder_ctv3 = BITCT_TO_ALIGNED_WORDCT(founder_ct);
- uint32_t founder_ctsplit = 3 * founder_ctv3;
+ uint32_t founder_ctaw3 = BITCT_TO_ALIGNED_WORDCT(founder_ct);
+ uint32_t founder_ctsplit = 3 * founder_ctaw3;
uintptr_t* geno1 = g_ld_geno1;
uintptr_t* zmiss1 = g_epi_zmiss1;
uintptr_t* sex_male = g_ld_sex_male;
@@ -5150,17 +5150,17 @@ THREAD_RET_TYPE ld_dprime_thread(void* arg) {
is_x1 = (block_idx1 >= xstart1) && (block_idx1 < xend1);
nm_fixed = is_set_ul(zmiss1, block_idx1);
cur_geno1 = &(geno1[block_idx1 * founder_ctsplit]);
- tot1[0] = popcount_longs(cur_geno1, founder_ctv3);
- tot1[1] = popcount_longs(&(cur_geno1[founder_ctv3]), founder_ctv3);
- tot1[2] = popcount_longs(&(cur_geno1[2 * founder_ctv3]), founder_ctv3);
+ tot1[0] = popcount_longs(cur_geno1, founder_ctaw3);
+ tot1[1] = popcount_longs(&(cur_geno1[founder_ctaw3]), founder_ctaw3);
+ tot1[2] = popcount_longs(&(cur_geno1[2 * founder_ctaw3]), founder_ctaw3);
if (is_x1 || x2_present) {
memcpy(cur_geno1_male, cur_geno1, founder_ctsplit * sizeof(intptr_t));
- bitvec_and(sex_male, founder_ctv3, cur_geno1_male);
- tot1[3] = popcount_longs(cur_geno1_male, founder_ctv3);
- bitvec_and(sex_male, founder_ctv3, &(cur_geno1_male[founder_ctv3]));
- tot1[4] = popcount_longs(&(cur_geno1_male[founder_ctv3]), founder_ctv3);
- bitvec_and(sex_male, founder_ctv3, &(cur_geno1_male[2 * founder_ctv3]));
- tot1[5] = popcount_longs(&(cur_geno1_male[2 * founder_ctv3]), founder_ctv3);
+ bitvec_and(sex_male, founder_ctaw3, cur_geno1_male);
+ tot1[3] = popcount_longs(cur_geno1_male, founder_ctaw3);
+ bitvec_and(sex_male, founder_ctaw3, &(cur_geno1_male[founder_ctaw3]));
+ tot1[4] = popcount_longs(&(cur_geno1_male[founder_ctaw3]), founder_ctaw3);
+ bitvec_and(sex_male, founder_ctaw3, &(cur_geno1_male[2 * founder_ctaw3]));
+ tot1[5] = popcount_longs(&(cur_geno1_male[2 * founder_ctaw3]), founder_ctaw3);
}
cur_geno2 = &(geno2[block_idx2 * founder_ctsplit]);
rptr = &(results[2 * block_idx1 * marker_idx2_maxw]);
@@ -5168,7 +5168,7 @@ THREAD_RET_TYPE ld_dprime_thread(void* arg) {
cur_tot2 = &(tot2[block_idx2 * 3]);
cur_zmiss2 = is_set_ul(zmiss2, block_idx2);
if (nm_fixed) {
- two_locus_count_table_zmiss1(cur_geno1, cur_geno2, counts, founder_ctv3, cur_zmiss2);
+ two_locus_count_table_zmiss1(cur_geno1, cur_geno2, counts, founder_ctaw3, cur_zmiss2);
if (cur_zmiss2) {
counts[2] = tot1[0] - counts[0] - counts[1];
counts[5] = tot1[1] - counts[3] - counts[4];
@@ -5177,7 +5177,7 @@ THREAD_RET_TYPE ld_dprime_thread(void* arg) {
counts[7] = cur_tot2[1] - counts[1] - counts[4];
counts[8] = cur_tot2[2] - counts[2] - counts[5];
} else {
- two_locus_count_table(cur_geno1, cur_geno2, counts, founder_ctv3, cur_zmiss2);
+ two_locus_count_table(cur_geno1, cur_geno2, counts, founder_ctaw3, cur_zmiss2);
if (cur_zmiss2) {
counts[2] = tot1[0] - counts[0] - counts[1];
counts[5] = tot1[1] - counts[3] - counts[4];
@@ -5186,7 +5186,7 @@ THREAD_RET_TYPE ld_dprime_thread(void* arg) {
}
is_x2 = ((block_idx2 < xend2) && (block_idx2 >= xstart2));
if (is_x1 || is_x2) {
- two_locus_count_table(cur_geno1_male, cur_geno2, &(counts[9]), founder_ctv3, cur_zmiss2);
+ two_locus_count_table(cur_geno1_male, cur_geno2, &(counts[9]), founder_ctaw3, cur_zmiss2);
if (cur_zmiss2) {
counts[11] = tot1[3] - counts[9] - counts[10];
counts[14] = tot1[4] - counts[12] - counts[13];
@@ -5244,8 +5244,8 @@ int32_t ld_report_dprime(pthread_t* threads, Ld_info* ldip, FILE* bedfile, uintp
uintptr_t unfiltered_sample_ct4 = (unfiltered_sample_ct + 3) / 4;
uintptr_t founder_ct = g_ld_founder_ct;
uintptr_t founder_ctl = BITCT_TO_WORDCT(founder_ct);
- uintptr_t founder_ctv3 = BITCT_TO_ALIGNED_WORDCT(founder_ct);
- uintptr_t founder_ctsplit = 3 * founder_ctv3;
+ uintptr_t founder_ctaw3 = BITCT_TO_ALIGNED_WORDCT(founder_ct);
+ uintptr_t founder_ctsplit = 3 * founder_ctaw3;
uintptr_t final_mask = get_final_mask(founder_ct);
uintptr_t orig_marker_ctm8 = g_ld_marker_ctm8;
uintptr_t marker_idx2_maxw = orig_marker_ctm8;
@@ -5347,8 +5347,8 @@ int32_t ld_report_dprime(pthread_t* threads, Ld_info* ldip, FILE* bedfile, uintp
goto ld_report_dprime_ret_NOMEM;
}
for (block_idx1 = 0; block_idx1 < idx1_block_size; block_idx1++) {
- g_ld_geno1[block_idx1 * founder_ctsplit + founder_ctv3 - 1] = 0;
- g_ld_geno1[block_idx1 * founder_ctsplit + 2 * founder_ctv3 - 1] = 0;
+ g_ld_geno1[block_idx1 * founder_ctsplit + founder_ctaw3 - 1] = 0;
+ g_ld_geno1[block_idx1 * founder_ctsplit + 2 * founder_ctaw3 - 1] = 0;
g_ld_geno1[block_idx1 * founder_ctsplit + founder_ctsplit - 1] = 0;
}
@@ -5375,8 +5375,8 @@ int32_t ld_report_dprime(pthread_t* threads, Ld_info* ldip, FILE* bedfile, uintp
idx2_block_size -= 4;
}
for (block_idx2 = 0; block_idx2 < idx2_block_size; block_idx2++) {
- g_ld_geno2[block_idx2 * founder_ctsplit + founder_ctv3 - 1] = 0;
- g_ld_geno2[block_idx2 * founder_ctsplit + 2 * founder_ctv3 - 1] = 0;
+ g_ld_geno2[block_idx2 * founder_ctsplit + founder_ctaw3 - 1] = 0;
+ g_ld_geno2[block_idx2 * founder_ctsplit + 2 * founder_ctaw3 - 1] = 0;
g_ld_geno2[block_idx2 * founder_ctsplit + founder_ctsplit - 1] = 0;
}
marker_uidx1 = next_unset_unsafe(marker_exclude_idx1, 0);
@@ -5546,7 +5546,7 @@ int32_t ld_report_dprime(pthread_t* threads, Ld_info* ldip, FILE* bedfile, uintp
if (is_haploid && hh_exists) {
haploid_fix(hh_exists, founder_include2, founder_male_include2, founder_ct, is_x, is_y, (unsigned char*)loadbuf);
}
- load_and_split3(nullptr, loadbuf, founder_ct, &(g_ld_geno1[block_idx1 * founder_ctsplit]), dummy_nm, dummy_nm, founder_ctv3, 0, 0, 1, &ulii);
+ load_and_split3(nullptr, loadbuf, founder_ct, &(g_ld_geno1[block_idx1 * founder_ctsplit]), dummy_nm, dummy_nm, founder_ctaw3, 0, 0, 1, &ulii);
if (ulii == 3) {
SET_BIT(block_idx1, g_epi_zmiss1);
}
@@ -5596,11 +5596,11 @@ int32_t ld_report_dprime(pthread_t* threads, Ld_info* ldip, FILE* bedfile, uintp
haploid_fix(hh_exists, founder_include2, founder_male_include2, founder_ct, is_x, is_y, (unsigned char*)loadbuf);
}
ulptr = &(g_ld_geno2[block_idx2 * founder_ctsplit]);
- load_and_split3(nullptr, loadbuf, founder_ct, ulptr, dummy_nm, dummy_nm, founder_ctv3, 0, 0, 1, &ulii);
+ load_and_split3(nullptr, loadbuf, founder_ct, ulptr, dummy_nm, dummy_nm, founder_ctaw3, 0, 0, 1, &ulii);
uiptr = &(g_epi_tot2[block_idx2 * 3]);
- uiptr[0] = popcount_longs(ulptr, founder_ctv3);
- uiptr[1] = popcount_longs(&(ulptr[founder_ctv3]), founder_ctv3);
- uiptr[2] = popcount_longs(&(ulptr[2 * founder_ctv3]), founder_ctv3);
+ uiptr[0] = popcount_longs(ulptr, founder_ctaw3);
+ uiptr[1] = popcount_longs(&(ulptr[founder_ctaw3]), founder_ctaw3);
+ uiptr[2] = popcount_longs(&(ulptr[2 * founder_ctaw3]), founder_ctaw3);
if (ulii == 3) {
SET_BIT(block_idx2, g_epi_zmiss2);
}
@@ -9366,10 +9366,10 @@ int32_t epistasis_report(pthread_t* threads, Epi_info* epi_ip, FILE* bedfile, ui
uintptr_t case_ctl2 = QUATERCT_TO_WORDCT(case_ct);
uintptr_t case_ctv2 = QUATERCT_TO_ALIGNED_WORDCT(case_ct);
uintptr_t ctrl_ctl2 = QUATERCT_TO_WORDCT(ctrl_ct);
- uintptr_t case_ctv3 = BITCT_TO_ALIGNED_WORDCT(case_ct);
- uintptr_t ctrl_ctv3 = BITCT_TO_ALIGNED_WORDCT(ctrl_ct);
- uintptr_t case_ctsplit = 3 * case_ctv3;
- uintptr_t ctrl_ctsplit = 3 * ctrl_ctv3;
+ uintptr_t case_ctaw3 = BITCT_TO_ALIGNED_WORDCT(case_ct);
+ uintptr_t ctrl_ctaw3 = BITCT_TO_ALIGNED_WORDCT(ctrl_ct);
+ uintptr_t case_ctsplit = 3 * case_ctaw3;
+ uintptr_t ctrl_ctsplit = 3 * ctrl_ctaw3;
uintptr_t pct = 1;
uintptr_t marker_uidx2 = 0;
uintptr_t marker_uidx2_trail = 0;
@@ -9728,7 +9728,7 @@ int32_t epistasis_report(pthread_t* threads, Epi_info* epi_ip, FILE* bedfile, ui
pct_thresh = tests_expected / 100;
if (is_case_only) {
g_epi_ctrl_ct = 0;
- ctrl_ctv3 = 0;
+ ctrl_ctaw3 = 0;
ctrl_ctsplit = 0;
memcpy(outname_end, ".epi.co", 8);
} else {
@@ -9798,11 +9798,11 @@ int32_t epistasis_report(pthread_t* threads, Epi_info* epi_ip, FILE* bedfile, ui
bigstack_alloc_ui(ulii, &g_epi_n_sig_ct1);
bigstack_alloc_ui(ulii, &g_epi_fail_ct1);
for (block_idx1 = 0; block_idx1 < idx1_block_size; block_idx1++) {
- g_epi_geno1[block_idx1 * tot_ctsplit + case_ctv3 - 1] = 0;
- g_epi_geno1[block_idx1 * tot_ctsplit + 2 * case_ctv3 - 1] = 0;
+ g_epi_geno1[block_idx1 * tot_ctsplit + case_ctaw3 - 1] = 0;
+ g_epi_geno1[block_idx1 * tot_ctsplit + 2 * case_ctaw3 - 1] = 0;
g_epi_geno1[block_idx1 * tot_ctsplit + case_ctsplit - 1] = 0;
- g_epi_geno1[block_idx1 * tot_ctsplit + case_ctsplit + ctrl_ctv3 - 1] = 0;
- g_epi_geno1[block_idx1 * tot_ctsplit + case_ctsplit + 2 * ctrl_ctv3 - 1] = 0;
+ g_epi_geno1[block_idx1 * tot_ctsplit + case_ctsplit + ctrl_ctaw3 - 1] = 0;
+ g_epi_geno1[block_idx1 * tot_ctsplit + case_ctsplit + 2 * ctrl_ctaw3 - 1] = 0;
g_epi_geno1[block_idx1 * tot_ctsplit + tot_ctsplit - 1] = 0;
}
if (is_triangular) {
@@ -9837,11 +9837,11 @@ int32_t epistasis_report(pthread_t* threads, Epi_info* epi_ip, FILE* bedfile, ui
idx2_block_size -= 16;
}
for (block_idx2 = 0; block_idx2 < idx2_block_size; block_idx2++) {
- g_epi_geno2[block_idx2 * tot_ctsplit + case_ctv3 - 1] = 0;
- g_epi_geno2[block_idx2 * tot_ctsplit + 2 * case_ctv3 - 1] = 0;
+ g_epi_geno2[block_idx2 * tot_ctsplit + case_ctaw3 - 1] = 0;
+ g_epi_geno2[block_idx2 * tot_ctsplit + 2 * case_ctaw3 - 1] = 0;
g_epi_geno2[block_idx2 * tot_ctsplit + case_ctsplit - 1] = 0;
- g_epi_geno2[block_idx2 * tot_ctsplit + case_ctsplit + ctrl_ctv3 - 1] = 0;
- g_epi_geno2[block_idx2 * tot_ctsplit + case_ctsplit + 2 * ctrl_ctv3 - 1] = 0;
+ g_epi_geno2[block_idx2 * tot_ctsplit + case_ctsplit + ctrl_ctaw3 - 1] = 0;
+ g_epi_geno2[block_idx2 * tot_ctsplit + case_ctsplit + 2 * ctrl_ctaw3 - 1] = 0;
g_epi_geno2[block_idx2 * tot_ctsplit + tot_ctsplit - 1] = 0;
}
marker_uidx = next_unset_ul_unsafe(marker_exclude1, marker_uidx_base);
@@ -9934,7 +9934,7 @@ int32_t epistasis_report(pthread_t* threads, Epi_info* epi_ip, FILE* bedfile, ui
goto epistasis_report_ret_READ_FAIL;
}
}
- if (load_and_split3(bedfile, loadbuf, unfiltered_sample_ct, &(g_epi_geno1[block_idx1 * tot_ctsplit]), pheno_nm, pheno_c, case_ctv3, ctrl_ctv3, IS_SET(marker_reverse, marker_uidx_tmp), is_case_only, &ulii)) {
+ if (load_and_split3(bedfile, loadbuf, unfiltered_sample_ct, &(g_epi_geno1[block_idx1 * tot_ctsplit]), pheno_nm, pheno_c, case_ctaw3, ctrl_ctaw3, IS_SET(marker_reverse, marker_uidx_tmp), is_case_only, &ulii)) {
goto epistasis_report_ret_READ_FAIL;
}
if (ulii) {
@@ -10047,18 +10047,18 @@ int32_t epistasis_report(pthread_t* threads, Epi_info* epi_ip, FILE* bedfile, ui
}
}
ulptr = &(g_epi_geno2[block_idx2 * tot_ctsplit]);
- if (load_and_split3(bedfile, loadbuf, unfiltered_sample_ct, ulptr, pheno_nm, pheno_c, case_ctv3, ctrl_ctv3, IS_SET(marker_reverse, marker_uidx2), is_case_only, &ulii)) {
+ if (load_and_split3(bedfile, loadbuf, unfiltered_sample_ct, ulptr, pheno_nm, pheno_c, case_ctaw3, ctrl_ctaw3, IS_SET(marker_reverse, marker_uidx2), is_case_only, &ulii)) {
goto epistasis_report_ret_READ_FAIL;
}
uiptr = &(g_epi_tot2[block_idx2 * tot_stride]);
- uiptr[0] = popcount_longs(ulptr, case_ctv3);
- uiptr[1] = popcount_longs(&(ulptr[case_ctv3]), case_ctv3);
- uiptr[2] = popcount_longs(&(ulptr[2 * case_ctv3]), case_ctv3);
+ uiptr[0] = popcount_longs(ulptr, case_ctaw3);
+ uiptr[1] = popcount_longs(&(ulptr[case_ctaw3]), case_ctaw3);
+ uiptr[2] = popcount_longs(&(ulptr[2 * case_ctaw3]), case_ctaw3);
if (!is_case_only) {
- ulptr = &(ulptr[case_ctv3 * 3]);
- uiptr[3] = popcount_longs(ulptr, ctrl_ctv3);
- uiptr[4] = popcount_longs(&(ulptr[ctrl_ctv3]), ctrl_ctv3);
- uiptr[5] = popcount_longs(&(ulptr[2 * ctrl_ctv3]), ctrl_ctv3);
+ ulptr = &(ulptr[case_ctaw3 * 3]);
+ uiptr[3] = popcount_longs(ulptr, ctrl_ctaw3);
+ uiptr[4] = popcount_longs(&(ulptr[ctrl_ctaw3]), ctrl_ctaw3);
+ uiptr[5] = popcount_longs(&(ulptr[2 * ctrl_ctaw3]), ctrl_ctaw3);
if (is_boost) {
boost_calc_p_bc(uiptr[0], uiptr[1], uiptr[2], uiptr[3], uiptr[4], uiptr[5], &(g_epi_boost_precalc2[block_idx2 * 6]));
}
@@ -10413,10 +10413,10 @@ int32_t indep_pairphase(Ld_info* ldip, FILE* bedfile, uintptr_t bed_offset, uint
uintptr_t unfiltered_sample_ctv2 = QUATERCT_TO_ALIGNED_WORDCT(unfiltered_sample_ct);
uintptr_t founder_ct = popcount_longs(founder_info, unfiltered_sample_ctv2 / 2);
uintptr_t founder_ctl = BITCT_TO_WORDCT(founder_ct);
- uintptr_t founder_ctv3 = BITCT_TO_ALIGNED_WORDCT(founder_ct);
+ uintptr_t founder_ctaw3 = BITCT_TO_ALIGNED_WORDCT(founder_ct);
// no actual case/control split here, but keep the variables the same to
// minimize divergence from ld_report_dprime()
- uintptr_t founder_ctsplit = 3 * founder_ctv3;
+ uintptr_t founder_ctsplit = 3 * founder_ctaw3;
uintptr_t final_mask = get_final_mask(founder_ct);
uintptr_t window_max = 1;
uintptr_t* founder_include2 = nullptr;
@@ -10533,8 +10533,8 @@ int32_t indep_pairphase(Ld_info* ldip, FILE* bedfile, uintptr_t bed_offset, uint
fill_all_bits(founder_ct, dummy_nm);
// bugfix: this loop must start at 0, not 1
for (ulii = 0; ulii < window_max; ulii++) {
- geno[ulii * founder_ctsplit + founder_ctv3 - 1] = 0;
- geno[ulii * founder_ctsplit + 2 * founder_ctv3 - 1] = 0;
+ geno[ulii * founder_ctsplit + founder_ctaw3 - 1] = 0;
+ geno[ulii * founder_ctsplit + 2 * founder_ctaw3 - 1] = 0;
geno[ulii * founder_ctsplit + founder_ctsplit - 1] = 0;
}
if ((chrom_info_ptr->xymt_codes[X_OFFSET] != -2) && is_set(chrom_info_ptr->chrom_mask, chrom_info_ptr->xymt_codes[X_OFFSET])) {
@@ -10562,10 +10562,10 @@ int32_t indep_pairphase(Ld_info* ldip, FILE* bedfile, uintptr_t bed_offset, uint
haploid_fix(hh_exists, founder_include2, founder_male_include2, founder_ct, is_x, is_y, (unsigned char*)loadbuf);
}
cur_geno1 = &(geno[ulii * founder_ctsplit]);
- load_and_split3(nullptr, loadbuf, founder_ct, cur_geno1, dummy_nm, dummy_nm, founder_ctv3, 0, 0, 1, &ulkk);
- cur_tots[ulii * 3] = popcount_longs(cur_geno1, founder_ctv3);
- cur_tots[ulii * 3 + 1] = popcount_longs(&(cur_geno1[founder_ctv3]), founder_ctv3);
- cur_tots[ulii * 3 + 2] = popcount_longs(&(cur_geno1[2 * founder_ctv3]), founder_ctv3);
+ load_and_split3(nullptr, loadbuf, founder_ct, cur_geno1, dummy_nm, dummy_nm, founder_ctaw3, 0, 0, 1, &ulkk);
+ cur_tots[ulii * 3] = popcount_longs(cur_geno1, founder_ctaw3);
+ cur_tots[ulii * 3 + 1] = popcount_longs(&(cur_geno1[founder_ctaw3]), founder_ctaw3);
+ cur_tots[ulii * 3 + 2] = popcount_longs(&(cur_geno1[2 * founder_ctaw3]), founder_ctaw3);
if ((!cur_tots[ulii * 3 + 1]) && ((!cur_tots[ulii * 3]) || (!cur_tots[ulii * 3 + 2]))) {
SET_BIT(uljj, pruned_arr);
cur_exclude_ct++;
@@ -10596,12 +10596,12 @@ int32_t indep_pairphase(Ld_info* ldip, FILE* bedfile, uintptr_t bed_offset, uint
nm_fixed = is_set_ul(zmiss, ulii);
if (is_x) {
memcpy(cur_geno1_male, cur_geno1, founder_ctsplit * sizeof(intptr_t));
- bitvec_and(sex_male_collapsed, founder_ctv3, cur_geno1_male);
- tot1[3] = popcount_longs(cur_geno1_male, founder_ctv3);
- bitvec_and(sex_male_collapsed, founder_ctv3, &(cur_geno1_male[founder_ctv3]));
- tot1[4] = popcount_longs(&(cur_geno1_male[founder_ctv3]), founder_ctv3);
- bitvec_and(sex_male_collapsed, founder_ctv3, &(cur_geno1_male[2 * founder_ctv3]));
- tot1[5] = popcount_longs(&(cur_geno1_male[2 * founder_ctv3]), founder_ctv3);
+ bitvec_and(sex_male_collapsed, founder_ctaw3, cur_geno1_male);
+ tot1[3] = popcount_longs(cur_geno1_male, founder_ctaw3);
+ bitvec_and(sex_male_collapsed, founder_ctaw3, &(cur_geno1_male[founder_ctaw3]));
+ tot1[4] = popcount_longs(&(cur_geno1_male[founder_ctaw3]), founder_ctaw3);
+ bitvec_and(sex_male_collapsed, founder_ctaw3, &(cur_geno1_male[2 * founder_ctaw3]));
+ tot1[5] = popcount_longs(&(cur_geno1_male[2 * founder_ctaw3]), founder_ctaw3);
}
for (; uljj < cur_window_size; uljj++) {
if (IS_SET(pruned_arr, live_indices[uljj])) {
@@ -10611,7 +10611,7 @@ int32_t indep_pairphase(Ld_info* ldip, FILE* bedfile, uintptr_t bed_offset, uint
cur_tot2 = &(cur_tots[uljj * 3]);
cur_zmiss2 = IS_SET(zmiss, uljj);
if (nm_fixed) {
- two_locus_count_table_zmiss1(cur_geno1, cur_geno2, counts, founder_ctv3, cur_zmiss2);
+ two_locus_count_table_zmiss1(cur_geno1, cur_geno2, counts, founder_ctaw3, cur_zmiss2);
if (cur_zmiss2) {
counts[2] = tot1[0] - counts[0] - counts[1];
counts[5] = tot1[1] - counts[3] - counts[4];
@@ -10620,7 +10620,7 @@ int32_t indep_pairphase(Ld_info* ldip, FILE* bedfile, uintptr_t bed_offset, uint
counts[7] = cur_tot2[1] - counts[1] - counts[4];
counts[8] = cur_tot2[2] - counts[2] - counts[5];
} else {
- two_locus_count_table(cur_geno1, cur_geno2, counts, founder_ctv3, cur_zmiss2);
+ two_locus_count_table(cur_geno1, cur_geno2, counts, founder_ctaw3, cur_zmiss2);
if (cur_zmiss2) {
counts[2] = tot1[0] - counts[0] - counts[1];
counts[5] = tot1[1] - counts[3] - counts[4];
@@ -10628,7 +10628,7 @@ int32_t indep_pairphase(Ld_info* ldip, FILE* bedfile, uintptr_t bed_offset, uint
}
}
if (is_x) {
- two_locus_count_table(cur_geno1_male, cur_geno2, &(counts[9]), founder_ctv3, cur_zmiss2);
+ two_locus_count_table(cur_geno1_male, cur_geno2, &(counts[9]), founder_ctaw3, cur_zmiss2);
if (cur_zmiss2) {
counts[11] = tot1[3] - counts[9] - counts[10];
counts[14] = tot1[4] - counts[12] - counts[13];
@@ -10746,10 +10746,10 @@ int32_t indep_pairphase(Ld_info* ldip, FILE* bedfile, uintptr_t bed_offset, uint
haploid_fix(hh_exists, founder_include2, founder_male_include2, founder_ct, is_x, is_y, (unsigned char*)loadbuf);
}
cur_geno1 = &(geno[cur_window_size * founder_ctsplit]);
- load_and_split3(nullptr, loadbuf, founder_ct, cur_geno1, dummy_nm, dummy_nm, founder_ctv3, 0, 0, 1, &ulkk);
- cur_tots[((uintptr_t)cur_window_size) * 3] = popcount_longs(cur_geno1, founder_ctv3);
- cur_tots[((uintptr_t)cur_window_size) * 3 + 1] = popcount_longs(&(cur_geno1[founder_ctv3]), founder_ctv3);
- cur_tots[((uintptr_t)cur_window_size) * 3 + 2] = popcount_longs(&(cur_geno1[2 * founder_ctv3]), founder_ctv3);
+ load_and_split3(nullptr, loadbuf, founder_ct, cur_geno1, dummy_nm, dummy_nm, founder_ctaw3, 0, 0, 1, &ulkk);
+ cur_tots[((uintptr_t)cur_window_size) * 3] = popcount_longs(cur_geno1, founder_ctaw3);
+ cur_tots[((uintptr_t)cur_window_size) * 3 + 1] = popcount_longs(&(cur_geno1[founder_ctaw3]), founder_ctaw3);
+ cur_tots[((uintptr_t)cur_window_size) * 3 + 2] = popcount_longs(&(cur_geno1[2 * founder_ctaw3]), founder_ctaw3);
if ((!cur_tots[((uintptr_t)cur_window_size) * 3 + 1]) && ((!cur_tots[((uintptr_t)cur_window_size) * 3]) || (!cur_tots[((uintptr_t)cur_window_size) * 3 + 2]))) {
SET_BIT(window_unfiltered_end, pruned_arr);
cur_exclude_ct++;
View it on GitLab: https://salsa.debian.org/med-team/plink1-9/-/commit/7b42b568b020db6213f36f1956e2fd4fff925ba2
--
View it on GitLab: https://salsa.debian.org/med-team/plink1-9/-/commit/7b42b568b020db6213f36f1956e2fd4fff925ba2
You're receiving this email because of your account on salsa.debian.org.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://alioth-lists.debian.net/pipermail/debian-med-commit/attachments/20201111/91646f1b/attachment-0001.html>
More information about the debian-med-commit
mailing list