Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

numerical issue resulting in negative hivpop -> inconsistent ART initiations #46

Open
jeffeaton opened this issue Nov 4, 2024 · 0 comments

Comments

@jeffeaton
Copy link
Collaborator

jeffeaton commented Nov 4, 2024

This is not expected to happen in actual applications. Can't quite figure out what to do to address.

Small negative in Xartelig_15plus resulting in large ratio greater than 0 ratio for artelig_hm[hm] / Xartelig_15plus in line here:

eppasm/src/eppasm.cpp

Lines 913 to 914 in 30d4902

( (1.0 - art_alloc_mxweight) * artelig_hm[hm] / Xartelig_15plus +
art_alloc_mxweight * expect_mort_artelig_hm[hm] / expect_mort_artelig15plus);

Error occurs at t = 48, hts = 7, g = 0, ha = 0, hm = 6. Printed values of variables:

grad: -0.000000
artnum_hts: 27292.800000
artinit_hts: 21488.420514
artinit_hm: -213.749612
artelig: 2.081747e-16
Xartelig_15plus: -4.821249e-15
expect_mort_artelig15plus: 2.146363e-314

Issue discovered in Rob's ago-mort-disagg.pjnz example file when setting incidence rate to 0 (debugging discrepancy in ART coverage among entrants).


Example code:

pjnz <- "~/Downloads/ago-mort-disagg.pjnz"

fp <- prepare_directincid(pjnz)
fp$incidinput[] <- 0.0
fp$cd4_nonaids_excess_mort <- array(0.0, dim(fp$cd4_mort))
fp$art_nonaids_excess_mort <- array(0.0, dim(fp$art_mort))

modC <- simmod(fp)
modR <- simmod(fp, VERSION = "R")

round(colSums(attr(modC, "hivpop")[,,,49],,1) + colSums(attr(modC, "artpop")[,,,,49],,2),2)
round(apply(modC[,,2,49], 2, tapply, fp$ss$ag.idx, sum), 2)

round(colSums(attr(modR, "hivpop")[,,,49],,1) + colSums(attr(modR, "artpop")[,,,,49],,2),2)
round(apply(modR[,,2,49], 2, tapply, fp$ss$ag.idx, sum), 2)

Output:

> devtools::load_all()
ℹ Loading eppasm
> pjnz <- "~/Downloads/ago-mort-disagg.pjnz"
> fp <- prepare_directincid(pjnz)
> fp$incidinput[] <- 0.0
> fp$cd4_nonaids_excess_mort <- array(0.0, dim(fp$cd4_mort))
> fp$art_nonaids_excess_mort <- array(0.0, dim(fp$art_mort))
> modC <- simmod(fp)
modR <- simmod(fp, VERSION = "R")

> > round(colSums(attr(modC, "hivpop")[,,,49],,1) + colSums(attr(modC, "artpop")[,,,,49],,2),2)
         [,1]    [,2]
 [1,] 1278.34 1298.64
 [2,] 1492.60 1542.64
 [3,] 1710.75 1810.65
 [4,]  932.02 1048.01
 [5,]  305.96  349.94
 [6,] -101.56   55.19
 [7,]   -1.80    3.44
 [8,]    0.00    0.00
 [9,]    0.00    0.00
> round(apply(modC[,,2,49], 2, tapply, fp$ss$ag.idx, sum), 2)
     [,1]    [,2]
1 1278.34 1298.64
2 1492.60 1542.64
3 1710.75 1810.65
4  932.02 1048.01
5  305.96  349.94
6  -56.40   55.19
7   -0.36    3.44
8    0.00    0.00
9    0.00    0.00
> round(colSums(attr(modR, "hivpop")[,,,49],,1) + colSums(attr(modR, "artpop")[,,,,49],,2),2)
         [,1]    [,2]
 [1,] 1286.68 1298.64
 [2,] 1492.65 1542.64
 [3,] 1710.75 1810.65
 [4,]  945.62 1048.01
 [5,]  305.96  349.94
 [6,]   47.99   55.19
 [7,]    2.93    3.44
 [8,]    0.00    0.00
 [9,]    0.00    0.00
> round(apply(modR[,,2,49], 2, tapply, fp$ss$ag.idx, sum), 2)
     [,1]    [,2]
1 1286.68 1298.64
2 1492.65 1542.64
3 1710.75 1810.65
4  945.62 1048.01
5  305.96  349.94
6   47.99   55.19
7    2.93    3.44
8    0.00    0.00
9    0.00    0.00
> 
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant