[med-svn] [Git][med-team/plink1-9][upstream] New upstream version 1.90~b6.5-180913
Dylan Aïssi
gitlab at salsa.debian.org
Sun Sep 16 14:53:09 BST 2018
Dylan Aïssi pushed to branch upstream at Debian Med / plink1.9
Commits:
0cb99996 by Dylan Aïssi at 2018-09-16T13:38:21Z
New upstream version 1.90~b6.5-180913
- - - - -
6 changed files:
- plink.c
- plink_calc.c
- plink_data.c
- plink_glm.c
- plink_ld.c
- plink_rserve.c
Changes:
=====================================
plink.c
=====================================
@@ -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
@@ -105,10 +105,10 @@ 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
- " "
+ ""
#ifdef STABLE_BUILD
" " // (don't want this when version number has two trailing digits)
#else
=====================================
plink_calc.c
=====================================
@@ -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 {
=====================================
plink_data.c
=====================================
@@ -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)) {
=====================================
plink_glm.c
=====================================
@@ -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;
}
}
=====================================
plink_ld.c
=====================================
@@ -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++) {
@@ -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;
@@ -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
@@ -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);
=====================================
plink_rserve.c
=====================================
@@ -63,6 +63,7 @@ int32_t rserve_call(char* rplugin_fname, char* rplugin_host_or_socket, int32_t r
Rinteger* r_l;
Rinteger* r_g;
Rdouble* r_data;
+ Rexp* r_eval_dummy; // bugfix (11 Sep 2018): must delete rc->eval return-val
double* pheno_d2;
double* covar_d2 = nullptr;
double* dptr;
@@ -81,6 +82,12 @@ int32_t rserve_call(char* rplugin_fname, char* rplugin_host_or_socket, int32_t r
uint32_t cur_data_len;
uint32_t uii;
int32_t ii;
+ if (!pheno_nm_ct) {
+ // bugfix (11 Sep 2018): this was segfaulting instead of printing an error
+ // message
+ logerrprint("Error: --R only processes samples with nonmissing phenotype values, and all\nphenotype values are missing.\n");
+ goto rserve_call_ret_1;
+ }
if (fopen_checked(rplugin_fname, "r", &infile)) {
goto rserve_call_ret_OPEN_FAIL;
}
@@ -198,13 +205,19 @@ int32_t rserve_call(char* rplugin_fname, char* rplugin_host_or_socket, int32_t r
rc->assign("n", r_n);
rc->assign("PHENO", r_p);
rc->assign("CLUSTER", r_s);
- rc->eval("CLUSTER[CLUSTER==-1] <- NA");
+ r_eval_dummy = rc->eval("CLUSTER[CLUSTER==-1] <- NA");
+ if (r_eval_dummy) {
+ delete r_eval_dummy;
+ }
if (covar_ct) {
r_cov = new Rdouble(covar_d2, pheno_nm_ct * covar_ct);
rc->assign("c", r_cov);
- rc->eval("COVAR<-matrix(c,nrow=n,byrow=T)");
+ r_eval_dummy = rc->eval("COVAR<-matrix(c,nrow=n,byrow=T)");
} else {
- rc->eval("COVAR<-NA");
+ r_eval_dummy = rc->eval("COVAR<-NA");
+ }
+ if (r_eval_dummy) {
+ delete r_eval_dummy;
}
memcpy(outname_end, ".auto.R", 8);
} else {
@@ -314,11 +327,20 @@ int32_t rserve_call(char* rplugin_fname, char* rplugin_host_or_socket, int32_t r
r_g = new Rinteger(geno_int_buf, block_size * pheno_nm_ct);
rc->assign("l", r_l);
rc->assign("g", r_g);
- rc->eval("GENO<-matrix(g,nrow=n,byrow=T)");
- rc->eval("GENO[GENO==-1] <- NA");
+ r_eval_dummy = rc->eval("GENO<-matrix(g,nrow=n,byrow=T)");
+ if (r_eval_dummy) {
+ delete r_eval_dummy;
+ }
+ r_eval_dummy = rc->eval("GENO[GENO==-1] <- NA");
+ if (r_eval_dummy) {
+ delete r_eval_dummy;
+ }
delete r_l;
delete r_g;
- rc->eval(inbuf_start);
+ r_eval_dummy = rc->eval(inbuf_start);
+ if (r_eval_dummy) {
+ delete r_eval_dummy;
+ }
r_data = (Rdouble*)rc->eval("Rplink(PHENO,GENO,CLUSTER,COVAR)");
if (r_data) {
dptr = r_data->doubleArray();
@@ -463,6 +485,7 @@ int32_t rserve_call(char* rplugin_fname, char* rplugin_host_or_socket, int32_t r
retval = RET_NETWORK;
break;
}
+ rserve_call_ret_1:
fclose_cond(infile);
fclose_cond(outfile);
bigstack_reset(bigstack_mark);
View it on GitLab: https://salsa.debian.org/med-team/plink1-9/commit/0cb9999606e20c7e42de1812e1ecb3a48438dca9
--
View it on GitLab: https://salsa.debian.org/med-team/plink1-9/commit/0cb9999606e20c7e42de1812e1ecb3a48438dca9
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/20180916/f32d00f9/attachment-0001.html>
More information about the debian-med-commit
mailing list