Loading plink.c +3 −3 Original line number Diff line number Diff line Loading @@ -93,7 +93,7 @@ static const char ver_str[] = #ifdef STABLE_BUILD "PLINK v1.90b6.4" "PLINK v1.90b6.5" #else "PLINK v1.90p" #endif Loading @@ -105,7 +105,7 @@ static const char ver_str[] = #else " 32-bit" #endif " (7 Aug 2018)"; " (13 Sep 2018)"; static const char ver_str2[] = // include leading space if day < 10, so character length stays the same "" Loading plink_calc.c +1 −1 Original line number Diff line number Diff line Loading @@ -3028,7 +3028,7 @@ int32_t groupdist_calc(pthread_t* threads, uint32_t unfiltered_sample_ct, uintpt dxx = g_reg_tot_x / dww; dhh_sd = sqrt((dhh_ssq / dww - dxx * dxx) / (dww - 1.0)); } if (!(g_case_ct * g_ctrl_ct)) { if (!(g_case_ct && g_ctrl_ct)) { dyy = 0.0; dhl_sd = 0.0; } else { Loading plink_data.c +1 −1 Original line number Diff line number Diff line Loading @@ -6184,7 +6184,7 @@ int32_t ped_to_bed(char* pedname, char* mapname, char* outname, char* outname_en if (bigstack_left() >= marker_ct * sample_ct4) { markers_per_pass = marker_ct; sprintf(g_logbuf, "Performing single-pass .bed write (%" PRIuPTR " variant%s, %" PRIuPTR " %s).\n", marker_ct, (marker_ct == 1)? "" : "s", sample_ct, species_str(sample_ct)); pass_ct = (marker_ct * sample_ct4)? 1 : 0; pass_ct = (marker_ct && sample_ct4)? 1 : 0; } else { if (!map_is_unsorted) { if (bigstack_alloc_ll(sample_ct, &line_starts)) { Loading plink_glm.c +17 −0 Original line number Diff line number Diff line Loading @@ -1658,12 +1658,29 @@ uint32_t logistic_regression(uint32_t sample_ct, uint32_t param_ct, float* vv, f return 1; } if (iteration >= 15) { // quasi-bugfix (11 Sep 2018): if fabsf(any coefficient) > 8e3, this is // almost certainly a form of convergence failure that didn't get // caught by the (delta_coef > 20.0) check due to a precision quirk. // (8e3 threshold ~= 1e-4 * 2^23, since floats have 23 bits of // precision) for (param_idx = 0; param_idx < param_ct; param_idx++) { if (fabsf(coef[param_idx]) > 8e3) { return 1; } } return 0; } } // Pons reported that 1.1e-3 was dangerous, so I agree with the decision to // tighten this threshold from 1e-3 to 1e-4. if (delta_coef < 1e-4) { // Be more conservative in throwing out results when we don't hit the // iteration limit. for (param_idx = 0; param_idx < param_ct; param_idx++) { if (fabsf(coef[param_idx]) > 6e4) { return 1; } } return 0; } } Loading plink_ld.c +16 −0 Original line number Diff line number Diff line Loading @@ -3233,6 +3233,7 @@ void fepi_counts_to_joint_effects_stats(uint32_t group_ct, uint32_t* counts, dou uint32_t uii; uint32_t ujj; uint32_t ukk; tot_inv_v[0] = 0.0; // gcc7 maybe-uninitialized warning dptr = dcounts; if (counts[0] && counts[1] && counts[2] && counts[3] && counts[4] && counts[5] && counts[6] && counts[7] && counts[8] && ((group_ct == 1) || (counts[9] && counts[10] && counts[11] && counts[12] && counts[13] && counts[14] && counts[15] && counts[16] && counts[17]))) { for (uii = 0; uii < group_ct; uii++) { Loading Loading @@ -4712,11 +4713,15 @@ THREAD_RET_TYPE epi_logistic_thread(void* arg) { cur_word2 = cur_word2 & (cur_word2 >> 1); cur_sample_ct -= popcount2_long((cur_word1 | cur_word2) & FIVEMASK); } unsigned char geno_pair_present[12]; if (cur_sample_ct <= 4) { goto epi_logistic_thread_regression_fail; } // 2. now populate covariate-major matrix with 16-byte-aligned, // trailing-entries-zeroed rows // quasi-bugfix (13 Sep 2018): reliably detect when this matrix is not // of full rank, and skip the regression in that case. memset(geno_pair_present, 0, 12); cur_sample_cta4 = round_up_pow2(cur_sample_ct, 4); for (widx = 0; widx < pheno_nm_ctl2; widx++) { sample_idx = widx * BITCT2; Loading @@ -4731,6 +4736,7 @@ THREAD_RET_TYPE epi_logistic_thread(void* arg) { uljj = cur_word2 & (3 * ONELU); if ((ulii != 3) && (uljj != 3)) { *fptr = 1.0; geno_pair_present[ulii + uljj * 4] = 1; fxx = (float)((intptr_t)ulii); fyy = (float)((intptr_t)uljj); // maybe this is faster with continuous writes instead of Loading @@ -4745,6 +4751,16 @@ THREAD_RET_TYPE epi_logistic_thread(void* arg) { cur_word2 >>= 2; } } if (!geno_pair_present[5]) { // not full rank if any 2x2 square in the 3x3 contingency table is // empty. if (((!geno_pair_present[0]) && (!geno_pair_present[1]) && (!geno_pair_present[4])) || ((!geno_pair_present[1]) && (!geno_pair_present[2]) && (!geno_pair_present[6])) || ((!geno_pair_present[4]) && (!geno_pair_present[8]) && (!geno_pair_present[9])) || ((!geno_pair_present[6]) && (!geno_pair_present[9]) && (!geno_pair_present[10]))) { goto epi_logistic_thread_regression_fail; } } if (cur_sample_ct < cur_sample_cta4) { loop_end = cur_sample_cta4 - cur_sample_ct; fill_float_zero(loop_end, fptr); Loading Loading
plink.c +3 −3 Original line number Diff line number Diff line Loading @@ -93,7 +93,7 @@ static const char ver_str[] = #ifdef STABLE_BUILD "PLINK v1.90b6.4" "PLINK v1.90b6.5" #else "PLINK v1.90p" #endif Loading @@ -105,7 +105,7 @@ static const char ver_str[] = #else " 32-bit" #endif " (7 Aug 2018)"; " (13 Sep 2018)"; static const char ver_str2[] = // include leading space if day < 10, so character length stays the same "" Loading
plink_calc.c +1 −1 Original line number Diff line number Diff line Loading @@ -3028,7 +3028,7 @@ int32_t groupdist_calc(pthread_t* threads, uint32_t unfiltered_sample_ct, uintpt dxx = g_reg_tot_x / dww; dhh_sd = sqrt((dhh_ssq / dww - dxx * dxx) / (dww - 1.0)); } if (!(g_case_ct * g_ctrl_ct)) { if (!(g_case_ct && g_ctrl_ct)) { dyy = 0.0; dhl_sd = 0.0; } else { Loading
plink_data.c +1 −1 Original line number Diff line number Diff line Loading @@ -6184,7 +6184,7 @@ int32_t ped_to_bed(char* pedname, char* mapname, char* outname, char* outname_en if (bigstack_left() >= marker_ct * sample_ct4) { markers_per_pass = marker_ct; sprintf(g_logbuf, "Performing single-pass .bed write (%" PRIuPTR " variant%s, %" PRIuPTR " %s).\n", marker_ct, (marker_ct == 1)? "" : "s", sample_ct, species_str(sample_ct)); pass_ct = (marker_ct * sample_ct4)? 1 : 0; pass_ct = (marker_ct && sample_ct4)? 1 : 0; } else { if (!map_is_unsorted) { if (bigstack_alloc_ll(sample_ct, &line_starts)) { Loading
plink_glm.c +17 −0 Original line number Diff line number Diff line Loading @@ -1658,12 +1658,29 @@ uint32_t logistic_regression(uint32_t sample_ct, uint32_t param_ct, float* vv, f return 1; } if (iteration >= 15) { // quasi-bugfix (11 Sep 2018): if fabsf(any coefficient) > 8e3, this is // almost certainly a form of convergence failure that didn't get // caught by the (delta_coef > 20.0) check due to a precision quirk. // (8e3 threshold ~= 1e-4 * 2^23, since floats have 23 bits of // precision) for (param_idx = 0; param_idx < param_ct; param_idx++) { if (fabsf(coef[param_idx]) > 8e3) { return 1; } } return 0; } } // Pons reported that 1.1e-3 was dangerous, so I agree with the decision to // tighten this threshold from 1e-3 to 1e-4. if (delta_coef < 1e-4) { // Be more conservative in throwing out results when we don't hit the // iteration limit. for (param_idx = 0; param_idx < param_ct; param_idx++) { if (fabsf(coef[param_idx]) > 6e4) { return 1; } } return 0; } } Loading
plink_ld.c +16 −0 Original line number Diff line number Diff line Loading @@ -3233,6 +3233,7 @@ void fepi_counts_to_joint_effects_stats(uint32_t group_ct, uint32_t* counts, dou uint32_t uii; uint32_t ujj; uint32_t ukk; tot_inv_v[0] = 0.0; // gcc7 maybe-uninitialized warning dptr = dcounts; if (counts[0] && counts[1] && counts[2] && counts[3] && counts[4] && counts[5] && counts[6] && counts[7] && counts[8] && ((group_ct == 1) || (counts[9] && counts[10] && counts[11] && counts[12] && counts[13] && counts[14] && counts[15] && counts[16] && counts[17]))) { for (uii = 0; uii < group_ct; uii++) { Loading Loading @@ -4712,11 +4713,15 @@ THREAD_RET_TYPE epi_logistic_thread(void* arg) { cur_word2 = cur_word2 & (cur_word2 >> 1); cur_sample_ct -= popcount2_long((cur_word1 | cur_word2) & FIVEMASK); } unsigned char geno_pair_present[12]; if (cur_sample_ct <= 4) { goto epi_logistic_thread_regression_fail; } // 2. now populate covariate-major matrix with 16-byte-aligned, // trailing-entries-zeroed rows // quasi-bugfix (13 Sep 2018): reliably detect when this matrix is not // of full rank, and skip the regression in that case. memset(geno_pair_present, 0, 12); cur_sample_cta4 = round_up_pow2(cur_sample_ct, 4); for (widx = 0; widx < pheno_nm_ctl2; widx++) { sample_idx = widx * BITCT2; Loading @@ -4731,6 +4736,7 @@ THREAD_RET_TYPE epi_logistic_thread(void* arg) { uljj = cur_word2 & (3 * ONELU); if ((ulii != 3) && (uljj != 3)) { *fptr = 1.0; geno_pair_present[ulii + uljj * 4] = 1; fxx = (float)((intptr_t)ulii); fyy = (float)((intptr_t)uljj); // maybe this is faster with continuous writes instead of Loading @@ -4745,6 +4751,16 @@ THREAD_RET_TYPE epi_logistic_thread(void* arg) { cur_word2 >>= 2; } } if (!geno_pair_present[5]) { // not full rank if any 2x2 square in the 3x3 contingency table is // empty. if (((!geno_pair_present[0]) && (!geno_pair_present[1]) && (!geno_pair_present[4])) || ((!geno_pair_present[1]) && (!geno_pair_present[2]) && (!geno_pair_present[6])) || ((!geno_pair_present[4]) && (!geno_pair_present[8]) && (!geno_pair_present[9])) || ((!geno_pair_present[6]) && (!geno_pair_present[9]) && (!geno_pair_present[10]))) { goto epi_logistic_thread_regression_fail; } } if (cur_sample_ct < cur_sample_cta4) { loop_end = cur_sample_cta4 - cur_sample_ct; fill_float_zero(loop_end, fptr); Loading