[med-svn] [plink1.9] 01/02: Imported Upstream version 1.90~b3.44-161117
Dylan Aïssi
bob.dybian-guest at moszumanska.debian.org
Fri Nov 18 22:49:06 UTC 2016
This is an automated email from the git hooks/post-receive script.
bob.dybian-guest pushed a commit to branch master
in repository plink1.9.
commit b7402f1541a3f355b3a859be290cd6710e9dde25
Author: Dylan Aïssi <bob.dybian at gmail.com>
Date: Fri Nov 18 23:48:32 2016 +0100
Imported Upstream version 1.90~b3.44-161117
---
Makefile | 29 +++--
plink.c | 6 +-
plink_assoc.c | 5 +-
plink_common.c | 4 +-
plink_dosage.c | 2 +-
plink_family.c | 8 +-
plink_filter.c | 325 +++++----------------------------------------------------
plink_glm.c | 28 ++---
plink_ld.c | 5 +-
9 files changed, 75 insertions(+), 337 deletions(-)
diff --git a/Makefile b/Makefile
index 0235f0e..3529166 100644
--- a/Makefile
+++ b/Makefile
@@ -7,6 +7,11 @@
# (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 LINKFLAGS_EXTRA="<opts>"'.
+LINKFLAGS_EXTRA =
+
+.PHONY: clean
# should autodetect system
SYS = UNIX
@@ -49,16 +54,28 @@ ifdef NO_LAPACK
BLASFLAGS=
endif
-SRC = plink.c plink_assoc.c plink_calc.c plink_cluster.c plink_cnv.c plink_common.c plink_data.c plink_dosage.c plink_family.c plink_filter.c plink_glm.c plink_help.c plink_homozyg.c plink_lasso.c plink_ld.c plink_matrix.c plink_misc.c plink_perm.c plink_rserve.c plink_set.c plink_stats.c SFMT.c dcdflib.c pigz.c yarn.c Rconnection.cc hfile.c bgzf.c
+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 "plink1". (Previously,
# "plink1"/"plink2" was suggested here; that also works for now, but it may
# lead to minor problems when PLINK 2.0 is released.)
-plink: $(SRC)
- g++ $(CFLAGS) $(SRC) -o plink $(BLASFLAGS) $(LINKFLAGS) -L. $(ZLIB)
+plink: $(OBJS)
+ g++ $^ $(LINKFLAGS_EXTRA) $(BLASFLAGS) $(LINKFLAGS) -L. $(ZLIB) -o $@
+
+plinkw: $(OBJS)
+ gfortran -O2 $^ $(LINKFLAGS_EXTRA) -Wl,-Bstatic $(BLASFLAGS) $(LINKFLAGS) -L. $(ZLIB) -o $@
+
+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).
+
+%.o: %.c
+ g++ -c $(CFLAGS) $< -o $@
-plinkw: $(SRC)
- g++ $(CFLAGS) $(SRC) -c
- gfortran -O2 $(OBJ) -o plink -Wl,-Bstatic $(BLASFLAGS) $(LINKFLAGS) -L. $(ZLIB)
+%.o: %.cc
+ g++ -x c++ -c $(CFLAGS) $< -o $@
diff --git a/plink.c b/plink.c
index beca41d..0ea2355 100644
--- a/plink.c
+++ b/plink.c
@@ -93,7 +93,7 @@
static const char ver_str[] =
#ifdef STABLE_BUILD
- "PLINK v1.90b3.43"
+ "PLINK v1.90b3.44"
#else
"PLINK v1.90p"
#endif
@@ -105,7 +105,7 @@ static const char ver_str[] =
#else
" 32-bit"
#endif
- " (10 Oct 2016)";
+ " (17 Nov 2016)";
static const char ver_str2[] =
// include leading space if day < 10, so character length stays the same
""
@@ -2693,10 +2693,12 @@ int32_t rerun(uint32_t rerun_argv_pos, uint32_t rerun_parameter_present, int32_t
while (0) {
rerun_ret_NOMEM:
print_ver();
+ fputs(errstr_nomem, stderr);
retval = RET_NOMEM;
break;
rerun_ret_OPEN_FAIL:
print_ver();
+ fputs("Error: Failed to open --rerun file.\n", stderr);
retval = RET_OPEN_FAIL;
break;
rerun_ret_LONG_LINE:
diff --git a/plink_assoc.c b/plink_assoc.c
index f41f0ad..75bb8e3 100644
--- a/plink_assoc.c
+++ b/plink_assoc.c
@@ -9543,13 +9543,14 @@ int32_t gxe_assoc(FILE* bedfile, uintptr_t bed_offset, char* outname, char* outn
fflush(stdout);
}
}
-
+ if (fclose_null(&outfile)) {
+ goto gxe_assoc_ret_WRITE_FAIL;
+ }
if (pct >= 10) {
putc_unlocked('\b', stdout);
}
fputs("\b\b", stdout);
logprint("done.\n");
- fclose_cond(outfile);
while (0) {
gxe_assoc_ret_NOMEM:
diff --git a/plink_common.c b/plink_common.c
index 8399099..8fb3a7b 100644
--- a/plink_common.c
+++ b/plink_common.c
@@ -4368,7 +4368,6 @@ void init_species(uint32_t species_code, Chrom_info* chrom_info_ptr) {
const uint32_t species_autosome_ct[] = {22, 29, 38, 31, 19, 12, 26};
const uint32_t species_max_code[] = {26, 33, 42, 33, 21, 12, 28};
fill_ulong_zero(CHROM_MASK_WORDS, chrom_info_ptr->chrom_mask);
- fill_ulong_zero(CHROM_MASK_WORDS, chrom_info_ptr->haploid_mask);
chrom_info_ptr->output_encoding = 0;
chrom_info_ptr->zero_extra_chroms = 0;
chrom_info_ptr->species = species_code;
@@ -4377,6 +4376,9 @@ void init_species(uint32_t species_code, Chrom_info* chrom_info_ptr) {
g_species_plural = species_plural_constants[species_code];
if (species_code != SPECIES_UNKNOWN) {
// these are assumed to be already initialized in the SPECIES_UNKNOWN case
+
+ // bugfix: haploid_mask was being cleared in --chr-set case
+ fill_ulong_zero(CHROM_MASK_WORDS, chrom_info_ptr->haploid_mask);
memcpy(chrom_info_ptr->xymt_codes, &(species_xymt_codes[species_code * XYMT_OFFSET_CT]), XYMT_OFFSET_CT * sizeof(int32_t));
chrom_info_ptr->autosome_ct = species_autosome_ct[species_code];
chrom_info_ptr->max_code = species_max_code[species_code];
diff --git a/plink_dosage.c b/plink_dosage.c
index 03f9995..b2e976e 100644
--- a/plink_dosage.c
+++ b/plink_dosage.c
@@ -1052,7 +1052,7 @@ int32_t plink1_dosage(Dosage_info* doip, char* famname, char* mapname, char* out
for (sample_uidx = 0, sample_idx = 0; sample_idx < sample_ct; sample_uidx++, sample_idx++) {
next_unset_unsafe_ck(sample_exclude, &sample_uidx);
if (is_set(sex_nm, sample_uidx)) {
- covar_d[sample_idx] = (double)((int32_t)is_set(sex_male, sample_idx));
+ covar_d[sample_idx] = (double)((int32_t)is_set(sex_male, sample_uidx));
} else {
CLEAR_BIT(sample_idx, covar_nm);
covar_d[sample_idx] = missing_phenod;
diff --git a/plink_family.c b/plink_family.c
index df0286d..fcb23b8 100644
--- a/plink_family.c
+++ b/plink_family.c
@@ -5099,7 +5099,7 @@ void flip_precalc(uint32_t lm_ct, double* qfam_w, double* pheno_d2, uintptr_t* n
// * geno_ssq is constant, so we precompute it.
// * The inner loop can also skip W=0 samples. So we clear those nm_lm bits.
// * We can also precompute geno_sum and qt_g_prod under the assumption of no
- // flips, and then
+ // flips, and then patch them for each flip
double geno_sum = 0.0;
double geno_ssq = 0.0;
double qt_g_prod = 0.0;
@@ -5625,7 +5625,7 @@ int32_t qfam(pthread_t* threads, FILE* bedfile, uintptr_t bed_offset, char* outn
bigstack_alloc_d(round_up_pow2(fss_ct, CACHELINE_DBL) * qfam_thread_ct, &qfam_b) ||
bigstack_alloc_d(round_up_pow2(lm_ct, CACHELINE_DBL) * qfam_thread_ct, &qfam_w) ||
bigstack_alloc_ui(fss_ct, &dummy_perm) ||
- bigstack_alloc_ul(fss_ctl, &dummy_flip)) {
+ bigstack_alloc_ul(flip_ctl, &dummy_flip)) {
goto qfam_ret_NOMEM;
}
for (uii = 0, ujj = 0, ukk = 0; ujj < sample_ct; uii++, ujj++) {
@@ -5652,7 +5652,7 @@ int32_t qfam(pthread_t* threads, FILE* bedfile, uintptr_t bed_offset, char* outn
for (uii = 0; uii < fss_ct; uii++) {
dummy_perm[uii] = uii;
}
- fill_ulong_zero(fss_ctl, dummy_flip);
+ fill_ulong_zero(flip_ctl, dummy_flip);
LOGPRINTFWW("--qfam-%s: Permuting %" PRIuPTR " families/singletons, and including %u %s in linear regression.\n", flag_suffix, fss_ct, lm_ct, g_species_plural);
LOGPRINTFWW5("Writing report to %s ... ", outname);
@@ -5698,7 +5698,7 @@ int32_t qfam(pthread_t* threads, FILE* bedfile, uintptr_t bed_offset, char* outn
}
ulptr = &(ulptr[lm_ctl]);
}
- fill_ulong_zero(fss_ctl, dummy_flip);
+ fill_ulong_zero(lm_ctl, dummy_flip);
} else {
for (ulii = 0; ulii < cur_perm_ct; ulii++) {
uint32_permute(&(g_qfam_permute[ulii * fss_ct]), &(precomputed_mods[-1]), &g_sfmt, fss_ct);
diff --git a/plink_filter.c b/plink_filter.c
index d198347..a7a3692 100644
--- a/plink_filter.c
+++ b/plink_filter.c
@@ -1708,209 +1708,6 @@ int32_t mind_filter(FILE* bedfile, uintptr_t bed_offset, char* outname, char* ou
return retval;
}
-#ifdef __LP64__
-void freq_hwe_haploid_count_120v(__m128i* vptr, __m128i* vend, __m128i* maskvp, uint32_t* ct_nmp, uint32_t* ct_hmajp) {
- const __m128i m2 = {0x3333333333333333LLU, 0x3333333333333333LLU};
- const __m128i m4 = {0x0f0f0f0f0f0f0f0fLLU, 0x0f0f0f0f0f0f0f0fLLU};
- const __m128i m8 = {0x00ff00ff00ff00ffLLU, 0x00ff00ff00ff00ffLLU};
- __m128i loader;
- __m128i loader2;
- __m128i loader3;
- __m128i to_ct_nm1;
- __m128i to_ct_hmaj1;
- __m128i to_ct_nm2;
- __m128i to_ct_hmaj2;
- __univec acc_nm;
- __univec acc_hmaj;
-
- acc_nm.vi = _mm_setzero_si128();
- acc_hmaj.vi = _mm_setzero_si128();
- do {
- loader = *vptr++;
- loader3 = _mm_srli_epi64(loader, 1);
- loader2 = _mm_xor_si128(loader, loader3); // inverted
- loader = _mm_and_si128(loader, loader3);
- loader3 = *maskvp++;
- to_ct_nm1 = _mm_andnot_si128(loader2, loader3);
- to_ct_hmaj1 = _mm_and_si128(loader, loader3);
-
- loader = *vptr++;
- loader3 = _mm_srli_epi64(loader, 1);
- loader2 = _mm_xor_si128(loader, loader3); // inverted
- loader = _mm_and_si128(loader, loader3);
- loader3 = *maskvp++;
- to_ct_nm1 = _mm_add_epi64(to_ct_nm1, _mm_andnot_si128(loader2, loader3));
- to_ct_hmaj1 = _mm_add_epi64(to_ct_hmaj1, _mm_and_si128(loader, loader3));
-
- loader = *vptr++;
- loader3 = _mm_srli_epi64(loader, 1);
- loader2 = _mm_xor_si128(loader, loader3); // inverted
- loader = _mm_and_si128(loader, loader3);
- loader3 = *maskvp++;
- to_ct_nm1 = _mm_add_epi64(to_ct_nm1, _mm_andnot_si128(loader2, loader3));
- to_ct_hmaj1 = _mm_add_epi64(to_ct_hmaj1, _mm_and_si128(loader, loader3));
-
- to_ct_nm1 = _mm_add_epi64(_mm_and_si128(to_ct_nm1, m2), _mm_and_si128(_mm_srli_epi64(to_ct_nm1, 2), m2));
- to_ct_hmaj1 = _mm_add_epi64(_mm_and_si128(to_ct_hmaj1, m2), _mm_and_si128(_mm_srli_epi64(to_ct_hmaj1, 2), m2));
-
- loader = *vptr++;
- loader3 = _mm_srli_epi64(loader, 1);
- loader2 = _mm_xor_si128(loader, loader3); // inverted
- loader = _mm_and_si128(loader, loader3);
- loader3 = *maskvp++;
- to_ct_nm2 = _mm_andnot_si128(loader2, loader3);
- to_ct_hmaj2 = _mm_and_si128(loader, loader3);
-
- loader = *vptr++;
- loader3 = _mm_srli_epi64(loader, 1);
- loader2 = _mm_xor_si128(loader, loader3); // inverted
- loader = _mm_and_si128(loader, loader3);
- loader3 = *maskvp++;
- to_ct_nm2 = _mm_add_epi64(to_ct_nm2, _mm_andnot_si128(loader2, loader3));
- to_ct_hmaj2 = _mm_add_epi64(to_ct_hmaj2, _mm_and_si128(loader, loader3));
-
- loader = *vptr++;
- loader3 = _mm_srli_epi64(loader, 1);
- loader2 = _mm_xor_si128(loader, loader3); // inverted
- loader = _mm_and_si128(loader, loader3);
- loader3 = *maskvp++;
- to_ct_nm2 = _mm_add_epi64(to_ct_nm2, _mm_andnot_si128(loader2, loader3));
- to_ct_hmaj2 = _mm_add_epi64(to_ct_hmaj2, _mm_and_si128(loader, loader3));
-
- to_ct_nm1 = _mm_add_epi64(to_ct_nm1, _mm_add_epi64(_mm_and_si128(to_ct_nm2, m2), _mm_and_si128(_mm_srli_epi64(to_ct_nm2, 2), m2)));
- to_ct_hmaj1 = _mm_add_epi64(to_ct_hmaj1, _mm_add_epi64(_mm_and_si128(to_ct_hmaj2, m2), _mm_and_si128(_mm_srli_epi64(to_ct_hmaj2, 2), m2)));
-
- acc_nm.vi = _mm_add_epi64(acc_nm.vi, _mm_add_epi64(_mm_and_si128(to_ct_nm1, m4), _mm_and_si128(_mm_srli_epi64(to_ct_nm1, 4), m4)));
- acc_hmaj.vi = _mm_add_epi64(acc_hmaj.vi, _mm_add_epi64(_mm_and_si128(to_ct_hmaj1, m4), _mm_and_si128(_mm_srli_epi64(to_ct_hmaj1, 4), m4)));
- } while (vptr < vend);
- acc_nm.vi = _mm_add_epi64(_mm_and_si128(acc_nm.vi, m8), _mm_and_si128(_mm_srli_epi64(acc_nm.vi, 8), m8));
- acc_hmaj.vi = _mm_add_epi64(_mm_and_si128(acc_hmaj.vi, m8), _mm_and_si128(_mm_srli_epi64(acc_hmaj.vi, 8), m8));
- *ct_nmp += ((acc_nm.u8[0] + acc_nm.u8[1]) * 0x1000100010001LLU) >> 48;
- *ct_hmajp += ((acc_hmaj.u8[0] + acc_hmaj.u8[1]) * 0x1000100010001LLU) >> 48;
-}
-#else
-void freq_hwe_haploid_count_12(uintptr_t* lptr, uintptr_t* maskp, uint32_t* ct_nmp, uint32_t* ct_hmajp) {
- uintptr_t loader = *lptr++;
- uintptr_t loader3 = loader >> 1;
- uintptr_t loader2 = loader ^ (~loader3);
- uint32_t to_ct_nm1;
- uint32_t to_ct_hmaj1;
- uint32_t to_ct_nm2;
- uint32_t to_ct_hmaj2;
- uintptr_t partial_nm;
- uintptr_t partial_hmaj;
- loader &= loader3;
- loader3 = *maskp++;
- to_ct_nm1 = loader2 & loader3;
- to_ct_hmaj1 = loader & loader3;
-
- loader = *lptr++;
- loader3 = loader >> 1;
- loader2 = loader ^ (~loader3);
- loader &= loader3;
- loader3 = *maskp++;
- to_ct_nm1 += loader2 & loader3;
- to_ct_hmaj1 += loader & loader3;
-
- loader = *lptr++;
- loader3 = loader >> 1;
- loader2 = loader ^ (~loader3);
- loader &= loader3;
- loader3 = *maskp++;
- to_ct_nm1 += loader2 & loader3;
- to_ct_hmaj1 += loader & loader3;
-
- loader = *lptr++;
- loader3 = loader >> 1;
- loader2 = loader ^ (~loader3);
- loader &= loader3;
- loader3 = *maskp++;
- to_ct_nm2 = loader2 & loader3;
- to_ct_hmaj2 = loader & loader3;
-
- loader = *lptr++;
- loader3 = loader >> 1;
- loader2 = loader ^ (~loader3);
- loader &= loader3;
- loader3 = *maskp++;
- to_ct_nm2 += loader2 & loader3;
- to_ct_hmaj2 += loader & loader3;
-
- loader = *lptr++;
- loader3 = loader >> 1;
- loader2 = loader ^ (~loader3);
- loader &= loader3;
- loader3 = *maskp++;
- to_ct_nm2 += loader2 & loader3;
- to_ct_hmaj2 += loader & loader3;
-
- to_ct_nm1 = (to_ct_nm1 & 0x33333333) + ((to_ct_nm1 >> 2) & 0x33333333);
- to_ct_nm1 += (to_ct_nm2 & 0x33333333) + ((to_ct_nm2 >> 2) & 0x33333333);
- partial_nm = (to_ct_nm1 & 0x0f0f0f0f) + ((to_ct_nm1 >> 4) & 0x0f0f0f0f);
- to_ct_hmaj1 = (to_ct_hmaj1 & 0x33333333) + ((to_ct_hmaj1 >> 2) & 0x33333333);
- to_ct_hmaj1 += (to_ct_hmaj2 & 0x33333333) + ((to_ct_hmaj2 >> 2) & 0x33333333);
- partial_hmaj = (to_ct_hmaj1 & 0x0f0f0f0f) + ((to_ct_hmaj1 >> 4) & 0x0f0f0f0f);
-
- loader = *lptr++;
- loader3 = loader >> 1;
- loader2 = loader ^ (~loader3);
- loader &= loader3;
- loader3 = *maskp++;
- to_ct_nm1 = loader2 & loader3;
- to_ct_hmaj1 = loader & loader3;
-
- loader = *lptr++;
- loader3 = loader >> 1;
- loader2 = loader ^ (~loader3);
- loader &= loader3;
- loader3 = *maskp++;
- to_ct_nm1 += loader2 & loader3;
- to_ct_hmaj1 += loader & loader3;
-
- loader = *lptr++;
- loader3 = loader >> 1;
- loader2 = loader ^ (~loader3);
- loader &= loader3;
- loader3 = *maskp++;
- to_ct_nm1 += loader2 & loader3;
- to_ct_hmaj1 += loader & loader3;
-
- loader = *lptr++;
- loader3 = loader >> 1;
- loader2 = loader ^ (~loader3);
- loader &= loader3;
- loader3 = *maskp++;
- to_ct_nm2 = loader2 & loader3;
- to_ct_hmaj2 = loader & loader3;
-
- loader = *lptr++;
- loader3 = loader >> 1;
- loader2 = loader ^ (~loader3);
- loader &= loader3;
- loader3 = *maskp++;
- to_ct_nm2 += loader2 & loader3;
- to_ct_hmaj2 += loader & loader3;
-
- loader = *lptr;
- loader3 = loader >> 1;
- loader2 = loader ^ (~loader3);
- loader &= loader3;
- loader3 = *maskp;
- to_ct_nm2 += loader2 & loader3;
- to_ct_hmaj2 += loader & loader3;
-
- to_ct_nm1 = (to_ct_nm1 & 0x33333333) + ((to_ct_nm1 >> 2) & 0x33333333);
- to_ct_nm1 += (to_ct_nm2 & 0x33333333) + ((to_ct_nm2 >> 2) & 0x33333333);
- partial_nm += (to_ct_nm1 & 0x0f0f0f0f) + ((to_ct_nm1 >> 4) & 0x0f0f0f0f);
- to_ct_hmaj1 = (to_ct_hmaj1 & 0x33333333) + ((to_ct_hmaj1 >> 2) & 0x33333333);
- to_ct_hmaj1 += (to_ct_hmaj2 & 0x33333333) + ((to_ct_hmaj2 >> 2) & 0x33333333);
- partial_hmaj += (to_ct_hmaj1 & 0x0f0f0f0f) + ((to_ct_hmaj1 >> 4) & 0x0f0f0f0f);
-
- *ct_nmp += (partial_nm * 0x01010101) >> 24;
- *ct_hmajp += (partial_hmaj * 0x01010101) >> 24;
-}
-#endif
-
static inline void single_marker_freqs_and_hwe(uintptr_t unfiltered_sample_ctl2, uintptr_t* lptr, uintptr_t* sample_include2, uintptr_t* founder_include2, uintptr_t* founder_ctrl_include2, uintptr_t* founder_case_include2, uintptr_t sample_ct, uint32_t* ll_ctp, uint32_t* lh_ctp, uint32_t* hh_ctp, uint32_t sample_f_ct, uint32_t* ll_ctfp, uint32_t* lh_ctfp, uint32_t* hh_ctfp, uint32_t hwe_or_geno_needed, uintptr_t sample_f_ctrl_ct, uint32_t* ll_hwep, uint32_t* lh_hwep, uint32_t* hh_hwep, i [...]
// This is best understood from the bottom third up (which is the order it
// was written). It's way overkill for just determining genotype
@@ -2078,86 +1875,6 @@ static inline uint32_t nonmissing_present_diff(uintptr_t unfiltered_sample_ctl2,
return 0;
}
-static inline void haploid_single_marker_freqs(uintptr_t unfiltered_sample_ct, uintptr_t unfiltered_sample_ctl2, uintptr_t* lptr, uintptr_t* sample_include2, uintptr_t* founder_include2, uintptr_t sample_ct, uint32_t* ll_ctp, uint32_t* hh_ctp, uint32_t sample_f_ct, uint32_t* ll_ctfp, uint32_t* hh_ctfp, uint32_t* hethap_incr_ptr) {
- // Here, we interpret heterozygotes as missing.
- // Nonmissing: (genotype ^ (~(genotype >> 1))) & 0x5555...
- // Homozygote major: (genotype & (genotype >> 1)) & 0x5555...
- uint32_t tot_a = 0;
- uint32_t tot_b = 0;
- uint32_t tot_hmaj = 0;
- uint32_t tot_nm_f = 0;
- uint32_t tot_hmaj_f = 0;
- uintptr_t* lptr_end = &(lptr[unfiltered_sample_ctl2]);
- uint32_t hethap_incr;
- uint32_t tot_nm;
- uint32_t uii;
- uintptr_t loader;
- uintptr_t loader2;
- uintptr_t loader3;
- uintptr_t loader4;
-#ifdef __LP64__
- uintptr_t cur_decr = 120;
- uintptr_t* lptr_12x_end;
- unfiltered_sample_ctl2 -= unfiltered_sample_ctl2 % 12;
- while (unfiltered_sample_ctl2 >= 120) {
- haploid_single_marker_freqs_loop:
- lptr_12x_end = &(lptr[cur_decr]);
- // Given a buffer with PLINK binary genotypes for a single marker, let
- // A := genotype & 0x5555...
- // B := (genotype >> 1) & 0x5555...
- // C := A & B
- // Then,
- // popcount(C) = homozyg major ct
- // popcount(B) = het ct + homozyg major ct
- // popcount(A) = missing_ct + homozyg major ct
- // = sample_ct - homozyg minor ct - het ct
- count_3freq_1920b((__m128i*)lptr, (__m128i*)lptr_12x_end, (__m128i*)sample_include2, &tot_a, &tot_b, &tot_hmaj);
- freq_hwe_haploid_count_120v((__m128i*)lptr, (__m128i*)lptr_12x_end, (__m128i*)founder_include2, &tot_nm_f, &tot_hmaj_f);
- lptr = lptr_12x_end;
- sample_include2 = &(sample_include2[cur_decr]);
- founder_include2 = &(founder_include2[cur_decr]);
- unfiltered_sample_ctl2 -= cur_decr;
- }
- if (unfiltered_sample_ctl2) {
- cur_decr = unfiltered_sample_ctl2;
- goto haploid_single_marker_freqs_loop;
- }
-#else
- uintptr_t* lptr_twelve_end = &(lptr[unfiltered_sample_ctl2 - unfiltered_sample_ctl2 % 12]);
- while (lptr < lptr_twelve_end) {
- count_3freq_48b(lptr, sample_include2, &tot_a, &tot_b, &tot_hmaj);
- freq_hwe_haploid_count_12(lptr, founder_include2, &tot_nm_f, &tot_hmaj_f);
- lptr = &(lptr[12]);
- sample_include2 = &(sample_include2[12]);
- founder_include2 = &(founder_include2[12]);
- }
-#endif
- tot_nm = 2 * tot_hmaj + sample_ct - tot_a - tot_b;
- hethap_incr = tot_b - tot_hmaj;
- while (lptr < lptr_end) {
- loader = *lptr++;
- loader3 = loader >> 1;
- loader4 = *sample_include2++;
- // different from tot_nm_f because of +sample_ct above
- tot_nm -= popcount2_long(loader & loader4);
- hethap_incr += popcount2_long(loader3 & (~loader) & loader4);
- loader2 = loader ^ (~loader3); // nonmissing?
- loader &= loader3; // homozyg A2?
- uii = popcount2_long(loader & loader4);
- tot_nm += 2 * uii - popcount2_long(loader3 & loader4);
- // tot_nm += popcount2_long(loader2 & loader4);
- tot_hmaj += uii;
- loader3 = *founder_include2++;
- tot_nm_f += popcount2_long(loader2 & loader3);
- tot_hmaj_f += popcount2_long(loader & loader3);
- }
- *hh_ctp = tot_hmaj;
- *ll_ctp = tot_nm - tot_hmaj;
- *hh_ctfp = tot_hmaj_f;
- *ll_ctfp = tot_nm_f - tot_hmaj_f;
- *hethap_incr_ptr = hethap_incr;
-}
-
int32_t calc_freqs_and_hwe(FILE* bedfile, char* outname, char* outname_end, uintptr_t unfiltered_marker_ct, uintptr_t* marker_exclude, uintptr_t marker_ct, char* marker_ids, uintptr_t max_marker_id_len, uintptr_t unfiltered_sample_ct, uintptr_t* sample_exclude, uintptr_t sample_exclude_ct, char* sample_ids, uintptr_t max_sample_id_len, uintptr_t* founder_info, int32_t nonfounders, int32_t maf_succ, double* set_allele_freqs, uintptr_t bed_offset, uint32_t hwe_needed, uint32_t hwe_all, uin [...]
FILE* hhfile = nullptr;
uintptr_t unfiltered_sample_ct4 = (unfiltered_sample_ct + 3) / 4;
@@ -2243,7 +1960,6 @@ int32_t calc_freqs_and_hwe(FILE* bedfile, char* outname, char* outname_end, uint
uintptr_t ulii;
uint32_t sample_male_ct;
uint32_t sample_f_male_ct;
- uint32_t hethap_incr;
uint32_t nonmales_needed;
uint32_t males_needed;
uint32_t uii;
@@ -2555,39 +2271,46 @@ int32_t calc_freqs_and_hwe(FILE* bedfile, char* outname, char* outname_end, uint
} else if (!nonmissing_nonmale_y) {
nonmissing_nonmale_y = nonmissing_present_diff(unfiltered_sample_ctv2, loadbuf, sample_include2, sample_male_include2);
}
- haploid_single_marker_freqs(unfiltered_sample_ct, unfiltered_sample_ctv2, loadbuf, sample_male_include2, founder_male_include2, sample_male_ct, &ll_ct, &hh_ct, sample_f_male_ct, &ll_ctf, &hh_ctf, &hethap_incr);
+ // Quasi-bugfix, 13 Nov 2016:
+ // There was a haploid_single_marker_freqs() function called by
+ // calc_freqs_and_hwe(), which treated heterozygous haploids as
+ // missing. This was inconsistent with PLINK 1.07, and also caused
+ // --geno to be based on a different set of missing calls than
+ // --mind.
+ single_marker_freqs_and_hwe(unfiltered_sample_ctv2, loadbuf, sample_male_include2, founder_male_include2, nullptr, nullptr, sample_male_ct, &ll_ct, &lh_ct, &hh_ct, sample_f_male_ct, &ll_ctf, &lh_ctf, &hh_ctf, 0, 0, nullptr, nullptr, nullptr, 0, 0, nullptr, nullptr, nullptr);
if ((is_x || sample_male_ct) && (sample_ct - cur_oblig_missing)) {
if (is_x) {
if (!cur_oblig_missing) {
- cur_genotyping_rate = ((int32_t)(ll_ct + hh_ct + ukk)) * sample_ct_recip;
+ cur_genotyping_rate = ((int32_t)(ll_ct + lh_ct + hh_ct + ukk)) * sample_ct_recip;
} else {
- cur_genotyping_rate = ((int32_t)(ll_ct + hh_ct + ukk)) / ((double)((int32_t)(sample_ct - cur_oblig_missing)));
+ cur_genotyping_rate = ((int32_t)(ll_ct + lh_ct + hh_ct + ukk)) / ((double)((int32_t)(sample_ct - cur_oblig_missing)));
}
} else {
if (!cur_oblig_missing) {
- cur_genotyping_rate = ((int32_t)(ll_ct + hh_ct)) * male_ct_recip;
+ cur_genotyping_rate = ((int32_t)(ll_ct + lh_ct + hh_ct)) * male_ct_recip;
} else {
- cur_genotyping_rate = ((int32_t)(ll_ct + hh_ct)) / ((double)((int32_t)(sample_male_ct - cur_oblig_missing)));
+ cur_genotyping_rate = ((int32_t)(ll_ct + lh_ct + hh_ct)) / ((double)((int32_t)(sample_male_ct - cur_oblig_missing)));
}
}
} else {
+ lh_ct = 0;
cur_genotyping_rate = 0;
nonmissing_rate_tot_max -= 1;
}
} else {
- haploid_single_marker_freqs(unfiltered_sample_ct, unfiltered_sample_ctv2, loadbuf, sample_include2, founder_include2, sample_ct, &ll_ct, &hh_ct, sample_f_ct, &ll_ctf, &hh_ctf, &hethap_incr);
+ single_marker_freqs_and_hwe(unfiltered_sample_ctv2, loadbuf, sample_include2, founder_include2, nullptr, nullptr, sample_ct, &ll_ct, &lh_ct, &hh_ct, sample_f_ct, &ll_ctf, &lh_ctf, &hh_ctf, 0, 0, nullptr, nullptr, nullptr, 0, 0, nullptr, nullptr, nullptr);
if (!cur_oblig_missing) {
- cur_genotyping_rate = ((int32_t)(ll_ct + hh_ct)) * sample_ct_recip;
+ cur_genotyping_rate = ((int32_t)(ll_ct + lh_ct + hh_ct)) * sample_ct_recip;
} else {
if (sample_ct - cur_oblig_missing) {
- cur_genotyping_rate = ((int32_t)(ll_ct + hh_ct)) / ((double)((int32_t)(sample_ct - cur_oblig_missing)));
+ cur_genotyping_rate = ((int32_t)(ll_ct + lh_ct + hh_ct)) / ((double)((int32_t)(sample_ct - cur_oblig_missing)));
} else {
cur_genotyping_rate = 0;
nonmissing_rate_tot_max -= 1;
}
}
}
- if (hethap_incr) {
+ if (lh_ct) {
if (!hhfile) {
memcpy(outname_end, ".hh", 4);
if (fopen_checked(outname, "w", &hhfile)) {
@@ -2631,7 +2354,7 @@ int32_t calc_freqs_and_hwe(FILE* bedfile, char* outname, char* outname_end, uint
if (ferror(hhfile)) {
goto calc_freqs_and_hwe_ret_WRITE_FAIL;
}
- hethap_ct += hethap_incr;
+ hethap_ct += lh_ct;
}
hwe_hapl_allfs[marker_uidx] = ll_ctf;
hwe_haph_allfs[marker_uidx] = hh_ctf;
@@ -2748,9 +2471,9 @@ int32_t write_missingness_reports(FILE* bedfile, uintptr_t bed_offset, char* out
uint32_t chrom_idx;
uint32_t chrom_end;
uint32_t cur_tot;
- uint32_t is_x;
+ // uint32_t is_x;
uint32_t is_y;
- uint32_t is_haploid;
+ // uint32_t is_haploid;
uint32_t om_ycorr;
uint32_t clidx;
uint32_t uii;
@@ -2847,9 +2570,9 @@ int32_t write_missingness_reports(FILE* bedfile, uintptr_t bed_offset, char* out
chrom_end = chrom_info_ptr->chrom_fo_vidx_start[chrom_fo_idx + 1];
marker_uidx = next_unset(marker_exclude, chrom_info_ptr->chrom_fo_vidx_start[chrom_fo_idx], chrom_end);
if (marker_uidx < chrom_end) {
- is_x = (((int32_t)chrom_idx) == chrom_info_ptr->xymt_codes[X_OFFSET]);
+ // is_x = (((int32_t)chrom_idx) == chrom_info_ptr->xymt_codes[X_OFFSET]);
is_y = (((int32_t)chrom_idx) == chrom_info_ptr->xymt_codes[Y_OFFSET]);
- is_haploid = is_set(chrom_info_ptr->haploid_mask, chrom_idx);
+ // is_haploid = is_set(chrom_info_ptr->haploid_mask, chrom_idx);
if (!is_y) {
cur_nm = sample_include2;
cur_tot = sample_ct;
@@ -2870,9 +2593,9 @@ int32_t write_missingness_reports(FILE* bedfile, uintptr_t bed_offset, char* out
if (load_raw(unfiltered_sample_ct4, bedfile, loadbuf)) {
goto write_missingness_reports_ret_READ_FAIL;
}
- if (is_haploid) {
- haploid_fix(hh_exists, sample_include2, sample_male_include2, unfiltered_sample_ct, is_x, is_y, (unsigned char*)loadbuf);
- }
+ // if (is_haploid) {
+ // haploid_fix(hh_exists, sample_include2, sample_male_include2, unfiltered_sample_ct, is_x, is_y, (unsigned char*)loadbuf);
+ // }
lptr = loadbuf;
lptr2 = cur_nm;
cptr2 = fw_strcpy(plink_maxsnp, &(marker_ids[marker_uidx * max_marker_id_len]), cptr);
diff --git a/plink_glm.c b/plink_glm.c
index ebd5b06..09f9a8c 100644
--- a/plink_glm.c
+++ b/plink_glm.c
@@ -1792,7 +1792,7 @@ uint32_t glm_logistic(uintptr_t cur_batch_size, uintptr_t param_ct, uintptr_t sa
return perm_fail_ct;
}
-uint32_t glm_fill_design(uintptr_t* loadbuf_collapsed, double* fixed_covars_cov_major, uintptr_t sample_valid_ct, uint32_t cur_param_ct, uint32_t coding_flags, uint32_t xchr_model, uintptr_t condition_list_start_idx, uintptr_t interaction_start_idx, uintptr_t sex_start_idx, uintptr_t* active_params, uintptr_t* haploid_params, uint32_t include_sex, uint32_t male_x_01, uintptr_t* sex_male_collapsed, uint32_t is_nonx_haploid, double* cur_covars_cov_major, double* cur_covars_sample_major, ui [...]
+uint32_t glm_fill_design(uintptr_t* loadbuf_collapsed, double* fixed_covars_cov_major, uintptr_t sample_valid_ct, uint32_t cur_param_ct, uint32_t coding_flags, uintptr_t condition_list_start_idx, uintptr_t interaction_start_idx, uintptr_t sex_start_idx, uintptr_t* active_params, uintptr_t* haploid_params, uint32_t include_sex, uint32_t male_x_01, uintptr_t* sex_male_collapsed, uint32_t is_nonx_haploid, double* cur_covars_cov_major, double* cur_covars_sample_major, uint32_t standard_beta) {
double* dptr = cur_covars_cov_major;
uintptr_t* ulptr_end_init = &(loadbuf_collapsed[sample_valid_ct / BITCT2]);
uintptr_t fixed_covar_nonsex_ct = interaction_start_idx - condition_list_start_idx;
@@ -2210,7 +2210,7 @@ uint32_t glm_fill_design(uintptr_t* loadbuf_collapsed, double* fixed_covars_cov_
return missing_ct;
}
-uint32_t glm_fill_design_float(uintptr_t* loadbuf_collapsed, float* fixed_covars_cov_major, uintptr_t sample_valid_ct, uint32_t cur_param_ct, uint32_t coding_flags, uint32_t xchr_model, uintptr_t condition_list_start_idx, uintptr_t interaction_start_idx, uintptr_t sex_start_idx, uintptr_t* active_params, uintptr_t* haploid_params, uint32_t include_sex, uint32_t male_x_01, uintptr_t* sex_male_collapsed, uint32_t is_nonx_haploid, float* cur_covars_cov_major) {
+uint32_t glm_fill_design_float(uintptr_t* loadbuf_collapsed, float* fixed_covars_cov_major, uintptr_t sample_valid_ct, uint32_t cur_param_ct, uint32_t coding_flags, uintptr_t condition_list_start_idx, uintptr_t interaction_start_idx, uintptr_t sex_start_idx, uintptr_t* active_params, uintptr_t* haploid_params, uint32_t include_sex, uint32_t male_x_01, uintptr_t* sex_male_collapsed, uint32_t is_nonx_haploid, float* cur_covars_cov_major) {
// rows of cur_covars_cov_major must be 16-byte aligned
float* fptr = cur_covars_cov_major;
uintptr_t* ulptr_end_init = &(loadbuf_collapsed[sample_valid_ct / BITCT2]);
@@ -2690,7 +2690,6 @@ static uintptr_t g_cur_param_ct;
static uintptr_t g_cur_constraint_ct;
static uint32_t g_standard_beta;
static uint32_t g_coding_flags;
-static uint32_t g_glm_xchr_model;
static uintptr_t g_condition_list_start_idx;
static uintptr_t g_interaction_start_idx;
static uintptr_t g_sex_start_idx;
@@ -2742,7 +2741,6 @@ THREAD_RET_TYPE glm_linear_adapt_thread(void* arg) {
__CLPK_integer dgels_lwork = g_dgels_lwork;
uint32_t standard_beta = g_standard_beta;
uint32_t coding_flags = g_coding_flags;
- uint32_t glm_xchr_model = g_glm_xchr_model;
uintptr_t condition_list_start_idx = g_condition_list_start_idx;
uintptr_t interaction_start_idx = g_interaction_start_idx;
uintptr_t sex_start_idx = g_sex_start_idx;
@@ -2812,7 +2810,7 @@ THREAD_RET_TYPE glm_linear_adapt_thread(void* arg) {
stat_high = orig_stats[marker_idx] + EPSILON;
stat_low = orig_stats[marker_idx] - EPSILON;
loadbuf_ptr = &(loadbuf[marker_bidx * sample_valid_ctv2]);
- cur_missing_ct = glm_fill_design(loadbuf_ptr, fixed_covars_cov_major, sample_valid_ct, cur_param_ct, coding_flags, glm_xchr_model, condition_list_start_idx, interaction_start_idx, sex_start_idx, active_params, haploid_params, include_sex, male_x_01, sex_male_collapsed, is_nonx_haploid, cur_covars_cov_major, cur_covars_sample_major, standard_beta);
+ cur_missing_ct = glm_fill_design(loadbuf_ptr, fixed_covars_cov_major, sample_valid_ct, cur_param_ct, coding_flags, condition_list_start_idx, interaction_start_idx, sex_start_idx, active_params, haploid_params, include_sex, male_x_01, sex_male_collapsed, is_nonx_haploid, cur_covars_cov_major, cur_covars_sample_major, standard_beta);
cur_sample_valid_ct = sample_valid_ct - cur_missing_ct;
dgels_m = (int32_t)((uint32_t)cur_sample_valid_ct);
dgels_ldb = dgels_m;
@@ -2921,7 +2919,6 @@ THREAD_RET_TYPE glm_logistic_adapt_thread(void* arg) {
uintptr_t cur_param_cta4 = round_up_pow2(cur_param_ct, 4);
uintptr_t cur_constraint_ct = g_cur_constraint_ct;
uint32_t coding_flags = g_coding_flags;
- uint32_t glm_xchr_model = g_glm_xchr_model;
uintptr_t condition_list_start_idx = g_condition_list_start_idx;
uintptr_t interaction_start_idx = g_interaction_start_idx;
uintptr_t sex_start_idx = g_sex_start_idx;
@@ -2988,7 +2985,7 @@ THREAD_RET_TYPE glm_logistic_adapt_thread(void* arg) {
stat_high = orig_stats[marker_idx] + EPSILON;
stat_low = orig_stats[marker_idx] - EPSILON;
loadbuf_ptr = &(loadbuf[marker_bidx * sample_valid_ctv2]);
- cur_missing_ct = glm_fill_design_float(loadbuf_ptr, fixed_covars_cov_major, sample_valid_ct, cur_param_ct, coding_flags, glm_xchr_model, condition_list_start_idx, interaction_start_idx, sex_start_idx, active_params, haploid_params, include_sex, male_x_01, sex_male_collapsed, is_nonx_haploid, cur_covars_cov_major);
+ cur_missing_ct = glm_fill_design_float(loadbuf_ptr, fixed_covars_cov_major, sample_valid_ct, cur_param_ct, coding_flags, condition_list_start_idx, interaction_start_idx, sex_start_idx, active_params, haploid_params, include_sex, male_x_01, sex_male_collapsed, is_nonx_haploid, cur_covars_cov_major);
cur_sample_valid_ct = sample_valid_ct - cur_missing_ct;
success_2incr = 0;
cur_fail_ct = 0;
@@ -3073,7 +3070,6 @@ THREAD_RET_TYPE glm_linear_maxt_thread(void* arg) {
__CLPK_integer dgels_lwork = g_dgels_lwork;
uint32_t standard_beta = g_standard_beta;
uint32_t coding_flags = g_coding_flags;
- uint32_t glm_xchr_model = g_glm_xchr_model;
uintptr_t condition_list_start_idx = g_condition_list_start_idx;
uintptr_t interaction_start_idx = g_interaction_start_idx;
uintptr_t sex_start_idx = g_sex_start_idx;
@@ -3146,7 +3142,7 @@ THREAD_RET_TYPE glm_linear_maxt_thread(void* arg) {
stat_high = orig_stats[marker_idx] + EPSILON;
stat_low = orig_stats[marker_idx] - EPSILON;
loadbuf_ptr = &(loadbuf[marker_bidx * sample_valid_ctv2]);
- cur_missing_ct = glm_fill_design(loadbuf_ptr, fixed_covars_cov_major, sample_valid_ct, cur_param_ct, coding_flags, glm_xchr_model, condition_list_start_idx, interaction_start_idx, sex_start_idx, active_params, haploid_params, include_sex, male_x_01, sex_male_collapsed, is_nonx_haploid, cur_covars_cov_major, cur_covars_sample_major, standard_beta);
+ cur_missing_ct = glm_fill_design(loadbuf_ptr, fixed_covars_cov_major, sample_valid_ct, cur_param_ct, coding_flags, condition_list_start_idx, interaction_start_idx, sex_start_idx, active_params, haploid_params, include_sex, male_x_01, sex_male_collapsed, is_nonx_haploid, cur_covars_cov_major, cur_covars_sample_major, standard_beta);
cur_sample_valid_ct = sample_valid_ct - cur_missing_ct;
dgels_m = (int32_t)((uint32_t)cur_sample_valid_ct);
dgels_ldb = dgels_m;
@@ -3243,7 +3239,6 @@ THREAD_RET_TYPE glm_logistic_maxt_thread(void* arg) {
uintptr_t cur_param_cta4 = round_up_pow2(cur_param_ct, 4);
uintptr_t cur_constraint_ct = g_cur_constraint_ct;
uint32_t coding_flags = g_coding_flags;
- uint32_t glm_xchr_model = g_glm_xchr_model;
uintptr_t condition_list_start_idx = g_condition_list_start_idx;
uintptr_t interaction_start_idx = g_interaction_start_idx;
uintptr_t sex_start_idx = g_sex_start_idx;
@@ -3312,7 +3307,7 @@ THREAD_RET_TYPE glm_logistic_maxt_thread(void* arg) {
stat_high = orig_stats[marker_idx] + EPSILON;
stat_low = orig_stats[marker_idx] - EPSILON;
loadbuf_ptr = &(loadbuf[marker_bidx * sample_valid_ctv2]);
- cur_missing_ct = glm_fill_design_float(loadbuf_ptr, fixed_covars_cov_major, sample_valid_ct, cur_param_ct, coding_flags, glm_xchr_model, condition_list_start_idx, interaction_start_idx, sex_start_idx, active_params, haploid_params, include_sex, male_x_01, sex_male_collapsed, is_nonx_haploid, cur_covars_cov_major);
+ cur_missing_ct = glm_fill_design_float(loadbuf_ptr, fixed_covars_cov_major, sample_valid_ct, cur_param_ct, coding_flags, condition_list_start_idx, interaction_start_idx, sex_start_idx, active_params, haploid_params, include_sex, male_x_01, sex_male_collapsed, is_nonx_haploid, cur_covars_cov_major);
cur_sample_valid_ct = sample_valid_ct - cur_missing_ct;
success_2incr = 0;
// todo: try better starting position
@@ -3378,7 +3373,6 @@ THREAD_RET_TYPE glm_linear_set_thread(void* arg) {
__CLPK_integer dgels_lwork = g_dgels_lwork;
uint32_t standard_beta = g_standard_beta;
uint32_t coding_flags = g_coding_flags;
- uint32_t glm_xchr_model = g_glm_xchr_model;
uintptr_t condition_list_start_idx = g_condition_list_start_idx;
uintptr_t interaction_start_idx = g_interaction_start_idx;
uintptr_t sex_start_idx = g_sex_start_idx;
@@ -3420,7 +3414,7 @@ THREAD_RET_TYPE glm_linear_set_thread(void* arg) {
for (; marker_bidx < marker_bceil; marker_bidx++) {
msa_ptr = &(g_mperm_save_all[marker_bidx * perm_vec_ct]);
loadbuf_ptr = &(loadbuf[marker_bidx * sample_valid_ctv2]);
- cur_missing_ct = glm_fill_design(loadbuf_ptr, fixed_covars_cov_major, sample_valid_ct, cur_param_ct, coding_flags, glm_xchr_model, condition_list_start_idx, interaction_start_idx, sex_start_idx, active_params, haploid_params, include_sex, male_x_01, sex_male_collapsed, is_nonx_haploid, cur_covars_cov_major, cur_covars_sample_major, standard_beta);
+ cur_missing_ct = glm_fill_design(loadbuf_ptr, fixed_covars_cov_major, sample_valid_ct, cur_param_ct, coding_flags, condition_list_start_idx, interaction_start_idx, sex_start_idx, active_params, haploid_params, include_sex, male_x_01, sex_male_collapsed, is_nonx_haploid, cur_covars_cov_major, cur_covars_sample_major, standard_beta);
cur_sample_valid_ct = sample_valid_ct - cur_missing_ct;
dgels_m = (int32_t)((uint32_t)cur_sample_valid_ct);
dgels_ldb = dgels_m;
@@ -3481,7 +3475,6 @@ THREAD_RET_TYPE glm_logistic_set_thread(void* arg) {
uintptr_t param_ct_m1 = cur_param_ct - 1;
uintptr_t cur_param_cta4 = round_up_pow2(cur_param_ct, 4);
uint32_t coding_flags = g_coding_flags;
- uint32_t glm_xchr_model = g_glm_xchr_model;
uintptr_t condition_list_start_idx = g_condition_list_start_idx;
uintptr_t interaction_start_idx = g_interaction_start_idx;
uintptr_t sex_start_idx = g_sex_start_idx;
@@ -3515,7 +3508,7 @@ THREAD_RET_TYPE glm_logistic_set_thread(void* arg) {
for (; marker_bidx < marker_bceil; marker_bidx++) {
msa_ptr = &(g_mperm_save_all[marker_bidx * perm_vec_ct]);
loadbuf_ptr = &(loadbuf[marker_bidx * sample_valid_ctv2]);
- cur_missing_ct = glm_fill_design_float(loadbuf_ptr, fixed_covars_cov_major, sample_valid_ct, cur_param_ct, coding_flags, glm_xchr_model, condition_list_start_idx, interaction_start_idx, sex_start_idx, active_params, haploid_params, include_sex, male_x_01, sex_male_collapsed, is_nonx_haploid, cur_covars_cov_major);
+ cur_missing_ct = glm_fill_design_float(loadbuf_ptr, fixed_covars_cov_major, sample_valid_ct, cur_param_ct, coding_flags, condition_list_start_idx, interaction_start_idx, sex_start_idx, active_params, haploid_params, include_sex, male_x_01, sex_male_collapsed, is_nonx_haploid, cur_covars_cov_major);
cur_sample_valid_ct = sample_valid_ct - cur_missing_ct;
// todo: try better starting position
fill_float_zero(cur_param_cta4 * perm_vec_ct, coef);
@@ -3595,7 +3588,6 @@ int32_t glm_common_init(FILE* bedfile, uintptr_t bed_offset, uint32_t glm_modifi
}
g_standard_beta = standard_beta;
g_coding_flags = glm_modifier & (GLM_HETHOM | GLM_DOMINANT | GLM_RECESSIVE);
- g_glm_xchr_model = glm_xchr_model;
g_include_sex = 0;
g_male_x_01 = 0;
if (!glm_xchr_model) {
@@ -5131,7 +5123,7 @@ int32_t glm_linear_assoc(pthread_t* threads, FILE* bedfile, uintptr_t bed_offset
if (marker_idx_to_uidx) {
marker_idx_to_uidx[marker_idx3] = marker_uidx2;
}
- cur_missing_ct = glm_fill_design(loadbuf_ptr, g_fixed_covars_cov_major, sample_valid_ct, cur_param_ct, glm_modifier & (GLM_HETHOM | GLM_DOMINANT | GLM_RECESSIVE), glm_xchr_model, condition_list_start_idx, interaction_start_idx, sex_start_idx, active_params, haploid_params, g_include_sex, g_male_x_01, sex_male_collapsed, g_min_ploidy_1 && (!g_is_x), g_linear_mt[0].cur_covars_cov_major, g_linear_mt[0].cur_covars_sample_major, standard_beta);
+ cur_missing_ct = glm_fill_design(loadbuf_ptr, g_fixed_covars_cov_major, sample_valid_ct, cur_param_ct, glm_modifier & (GLM_HETHOM | GLM_DOMINANT | GLM_RECESSIVE), condition_list_start_idx, interaction_start_idx, sex_start_idx, active_params, haploid_params, g_include_sex, g_male_x_01, sex_male_collapsed, g_min_ploidy_1 && (!g_is_x), g_linear_mt[0].cur_covars_cov_major, g_linear_mt[0].cur_covars_sample_major, standard_beta);
cur_sample_valid_ct = sample_valid_ct - cur_missing_ct;
g_nm_cts[marker_idx3] = cur_sample_valid_ct;
if ((cur_sample_valid_ct > cur_param_ct) && (!glm_check_vif(glm_vif_thresh, cur_param_ct, cur_sample_valid_ct, g_linear_mt[0].cur_covars_cov_major, g_linear_mt[0].param_2d_buf, g_linear_mt[0].mi_buf, g_linear_mt[0].param_2d_buf2))) {
@@ -6598,7 +6590,7 @@ int32_t glm_logistic_assoc(pthread_t* threads, FILE* bedfile, uintptr_t bed_offs
if (marker_idx_to_uidx) {
marker_idx_to_uidx[marker_idx3] = marker_uidx2;
}
- cur_missing_ct = glm_fill_design_float(loadbuf_ptr, g_fixed_covars_cov_major_f, sample_valid_ct, cur_param_ct, glm_modifier & (GLM_HETHOM | GLM_DOMINANT | GLM_RECESSIVE), glm_xchr_model, condition_list_start_idx, interaction_start_idx, sex_start_idx, active_params, haploid_params, g_include_sex, g_male_x_01, sex_male_collapsed, g_min_ploidy_1 && (!g_is_x), g_logistic_mt[0].cur_covars_cov_major);
+ cur_missing_ct = glm_fill_design_float(loadbuf_ptr, g_fixed_covars_cov_major_f, sample_valid_ct, cur_param_ct, glm_modifier & (GLM_HETHOM | GLM_DOMINANT | GLM_RECESSIVE), condition_list_start_idx, interaction_start_idx, sex_start_idx, active_params, haploid_params, g_include_sex, g_male_x_01, sex_male_collapsed, g_min_ploidy_1 && (!g_is_x), g_logistic_mt[0].cur_covars_cov_major);
cur_sample_valid_ct = sample_valid_ct - cur_missing_ct;
g_nm_cts[marker_idx3] = cur_sample_valid_ct;
if (cur_sample_valid_ct > cur_param_ct) {
diff --git a/plink_ld.c b/plink_ld.c
index 3b59217..2c9f10a 100644
--- a/plink_ld.c
+++ b/plink_ld.c
@@ -3282,7 +3282,8 @@ void fepi_counts_to_joint_effects_stats(uint32_t group_ct, uint32_t* counts, dou
}
}
dptr2 = ivv;
- for (uii = 0; uii < group_ct; uii++) {
+ uii = 0;
+ do {
dptr = &(dcounts[uii * 9]);
dptr3 = &(invcounts[uii * 9]);
dxx = dptr[8];
@@ -3290,7 +3291,7 @@ void fepi_counts_to_joint_effects_stats(uint32_t group_ct, uint32_t* counts, dou
*dptr2++ = dxx * dptr[1] * dptr3[2] * dptr3[7];
*dptr2++ = dxx * dptr[3] * dptr3[5] * dptr3[6];
*dptr2++ = dxx * dptr[4] * dptr3[5] * dptr3[7];
- }
+ } while (++uii < group_ct);
use_reg_stat = (ivv[3] > 0.5) && ((group_ct == 1) || (ivv[7] > 0.5));
if (use_reg_stat) {
dptr2 = xiv;
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-med/plink1.9.git
More information about the debian-med-commit
mailing list