Skip to content

Commit

Permalink
Debugged internal code (bias_coeff_admix_fit) shared by `admix_prop…
Browse files Browse the repository at this point in the history
…_1d_linear` and `admix_prop_1d_circular` for edge cases. M1mac error only (all other systems were fine before)
  • Loading branch information
alexviiia committed Aug 9, 2021
1 parent bff026a commit 99bc358
Show file tree
Hide file tree
Showing 4 changed files with 35 additions and 19 deletions.
2 changes: 1 addition & 1 deletion DESCRIPTION
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
Package: bnpsd
Title: Simulate Genotypes from the BN-PSD Admixture Model
Version: 1.3.12
Version: 1.3.13
Authors@R: c(
person(given = "Alejandro",
family = "Ochoa",
Expand Down
8 changes: 8 additions & 0 deletions NEWS.md
Original file line number Diff line number Diff line change
Expand Up @@ -252,3 +252,11 @@ New functions and bug fixes dealing with reordering tree edges and tips.
- 6th CRAN submission.
- Removed "LazyData: true" from DESCRIPTION (to avoid a new "NOTE" on CRAN).
- Fixed spelling in documentation.

# bnpsd 1.3.13 (2021-08-09)

- 6th CRAN submission, second attempt.
- Debugged internal code (`bias_coeff_admix_fit`) shared by `admix_prop_1d_linear` and `admix_prop_1d_circular` for edge cases.
- Error was only observed on M1mac architecture (previous code worked on all other systems!).
- If a bias coefficient of 1 was desired, expected sigma to be `Inf`, but instead an error was encountered.
- Previous error message: "f() values at end points not of opposite sign" (in `stats::uniroot`)
34 changes: 23 additions & 11 deletions R/bias_coeff_admix_fit.R
Original file line number Diff line number Diff line change
Expand Up @@ -60,6 +60,12 @@ bias_coeff_admix_fit <- function(
if (bias_coeff < 0)
stop('Desired `bias_coeff` must be >= 0, passed ', bias_coeff)

# tolerance should be precision-dependent, but when configured with --disable-long-double the value of .Machine$double.eps doesn't change!
# so here we test for that config difference (by testing if .Machine$sizeof.longdouble is zero) and reduce the tolerance accordingly
tol <- .Machine$double.eps
if ( .Machine$sizeof.longdouble == 0 )
tol <- sqrt( tol )

# handle an edge case that is problematic in some systems (looking at you Apple M1)
if ( bias_coeff == 1 ) {
# in the existing cases the answer is:
Expand All @@ -75,7 +81,8 @@ bias_coeff_admix_fit <- function(
)
# get actual bias coefficient of this case
bias_coeff2 <- bias_coeff_admix(admix_proportions, coanc_subpops)
if ( bias_coeff2 == bias_coeff )
# equality is not always feasible due to machine precision, but if it's close enough let's do this
if ( abs( bias_coeff2 - bias_coeff ) < tol )
return( sigma )
}

Expand All @@ -92,12 +99,6 @@ bias_coeff_admix_fit <- function(
if (bias_coeff < bias_coeff_min)
stop('Desired `bias_coeff` must be greater than ', bias_coeff_min, ' (the minimum achievable with `sigma = 0`), passed ', bias_coeff, '. Tip: This minimum depends most strongly on the input `coanc_subpops`, so the main alternative to increasing `bias_coeff` is to change the shape of `coanc_subpops` (its overall scale does not change the minimum value)')

# tolerance should be precision-dependent, but when configured with --disable-long-double the value of .Machine$double.eps doesn't change!
# so here we test for that config difference (by testing if .Machine$sizeof.longdouble is zero) and reduce the tolerance accordingly
tol <- .Machine$double.eps
if ( .Machine$sizeof.longdouble == 0 )
tol <- sqrt( tol )

# function whose zero we want!
# "sigma" is the only parameter to optimize
# need this function inside here as it depends on extra parameters
Expand All @@ -117,14 +118,25 @@ bias_coeff_admix_fit <- function(
# get bias coefficient, return difference from desired value!
delta <- bias_coeff_admix(admix_proportions, coanc_subpops) - bias_coeff

# hack to prevent issues when attempting to fit extreme values (at boundaries)
# in particular, without this (and in particular in lower precision settings, such as --disable-long-double), uniroot below complains because delta has the same sign on both extrema (due to machine error, not a real sign issue). If the sign wasn't checked then the tolerance would find that the solution is at the boundary, so here we force the check earlier.
if ( abs(delta) < tol )
delta <- 0
## # hack to prevent issues when attempting to fit extreme values (at boundaries)
## # in particular, without this (and in particular in lower precision settings, such as --disable-long-double), uniroot below complains because delta has the same sign on both extrema (due to machine error, not a real sign issue). If the sign wasn't checked then the tolerance would find that the solution is at the boundary, so here we force the check earlier.
## if ( abs(delta) < tol )
## delta <- 0

# return delta now
return( delta )
}

# to avoid uniroot dying, check opposite ends and make sure sings are opposite (which is what that function checks and dies of if it fails)
# I believe the issue is handled earlier, but just in case, directly assess this problematic case and return hack answers if that happens
v0 <- bias_coeff_admix_objective( 0 ) # sigma = 0, s = admix_prop_bias_coeff_min
v1 <- bias_coeff_admix_objective( 1 ) # sigma = Inf, s = 1
# this function is monotonically increasing, so v0 should be negative and v1 positive
# handle cases when this is not so
if ( v0 >= 0 )
return( 0 )
if ( v1 <= 0 )
return( Inf )

# default tolerance of .Machine$double.eps^0.25 (~ 1e-4) was not good enough, reduced to ~ 2e-16
obj <- stats::uniroot(
Expand Down
10 changes: 3 additions & 7 deletions cran-comments.md
Original file line number Diff line number Diff line change
@@ -1,16 +1,12 @@
## Test environments
* local R: x86_64-redhat-linux-gnu R 4.0.5 (2021-03-31)
* local R-devel: x86_64-pc-linux-gnu R Under development (unstable) (2021-07-27 r80667)
* local R-devel-noLD: x86_64-pc-linux-gnu R Under development (unstable) (2021-07-27 r80667)
* local R-devel: x86_64-pc-linux-gnu R Under development (unstable) (2021-08-09 r80724)
* local R-devel-noLD: x86_64-pc-linux-gnu R Under development (unstable) (2021-08-09 r80724)
* win-builder release: x86_64-w64-mingw32 R 4.1.0 (2021-05-18)
* win-builder devel: x86_64-w64-mingw32 R Under development (unstable) (2021-07-25 r80663)
* win-builder devel: x86_64-w64-mingw32 R Under development (unstable) (2021-08-05 r80717)
* rhub (macos): x86_64-apple-darwin17.0 R 4.1.0 (2021-05-18)
* rhub (windows): x86_64-w64-mingw32 R Under development (unstable) (2021-07-03 r80596)
* rhub (solaris): i386-pc-solaris2.10 R 4.1.0 (2021-05-18)
* rhub (debian): x86_64-pc-linux-gnu R Under development (unstable) (2021-06-27 r80567)
* rhub (debian-noLD): x86_64-pc-linux-gnu R Under development (unstable) (2021-06-27 r80567)
* rhub (ubuntu): x86_64-pc-linux-gnu R Under development (unstable) (2021-06-27 r80567)
* rhub (fedora): x86_64-pc-linux-gnu R Under development (unstable) (2021-06-27 r80567)

## R CMD check results
There were no ERRORs, WARNINGs, or NOTEs.
Expand Down

0 comments on commit 99bc358

Please sign in to comment.