From dbab78436c0e4123b735facdf73b69c561e44cc8 Mon Sep 17 00:00:00 2001 From: Nicholas Brazeau Date: Tue, 20 Jun 2023 11:11:08 -0400 Subject: [PATCH 1/4] fix test for nonverbose --- tests/testthat/test-modelruns.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/testthat/test-modelruns.R b/tests/testthat/test-modelruns.R index a637e5f..c7f3451 100644 --- a/tests/testthat/test-modelruns.R +++ b/tests/testthat/test-modelruns.R @@ -21,6 +21,6 @@ test_that("model runs", { normalize_geodist = T, report_progress = T, return_verbose = F) - testthat::expect_length(mod, 4) + testthat::expect_length(mod, 6) }) From 0ee222a6f56554a41b3b8b82d97442dca198e193 Mon Sep 17 00:00:00 2001 From: Nicholas Brazeau Date: Sun, 3 Dec 2023 19:07:42 -0500 Subject: [PATCH 2/4] no geodist in data --- data-raw/IBD_simulation_data.R | 3 +++ data/IBD_simulation_data.rda | Bin 890 -> 738 bytes 2 files changed, 3 insertions(+) diff --git a/data-raw/IBD_simulation_data.R b/data-raw/IBD_simulation_data.R index 0023d46..4b779b0 100644 --- a/data-raw/IBD_simulation_data.R +++ b/data-raw/IBD_simulation_data.R @@ -131,6 +131,9 @@ IBD_simulation_data <- sim_IBDIBD(demesize = c(3,3,4), distmat = matrix(c(0,500, 100,750,0), nrow = 3), rate = 1e-3, Ft = 0.3) +# drop within deme +IBD_simulation_data <- IBD_simulation_data[IBD_simulation_data$geodist > 0, ] + #............................................................ # out #........................................................... diff --git a/data/IBD_simulation_data.rda b/data/IBD_simulation_data.rda index b97c7993867527627f84216dbd891ca68f89c833..661eda3b790047b4bc4c1ac1fa72ecbfccb2475b 100644 GIT binary patch literal 738 zcmV<80v-KAT4*^jL0KkKS%;_Oy8r@ffB*mc%JzJfIq}|W@>~DE->2QD(FG;e{d@oa zYWuhHJB7dj{0)i&Q7Pp#^%`VD6Uu1OlT0Q@r1XKK2-71^M8Q0!20-;aL_DGD0gQu5 zr>2C{LrsxAOwqKOP3M?5#TSTbAnHMCc;EEWYKamI-YG z@&al~z@SM;>D~qA%PU9wG;Q`b_iioq+fiQVy4Q=RSZbJ2l7B(Zex)Lel9C)l-wDYMmz}h2U1h$4YQ|wTXO%RN* zDgo=uktZ)cbpUi3h$Q+#U^+nX#o_*;nT7Bjs>#*EMXLp%!~+-V<|s)91L)+)zd-_o z&@y-;WB~{Pp&DCbXCeu87;qjs7Xyd7aOvR2UxH~U>PO+$(jXhKHTON%(#pu~oS-Sb@zybxq(TC5tQd1xDW-PhW>TM0 zlC?Fzr`|}WV0nrp6i^mHRxJ5Q6)1{Qao3z`i6&KBi&`{X74NsLp%H7BFmBe;U&ic+ Umj%IjVpIMuh^xuJ@MY}^I!kJ|M%Ue^#`Zl{eS=e ze*6FOJH@~PTm`mh!YES_>S;Eoq6A?g|p9-s|9Kw<`(G6sf#0004@pfUgr0000D z003wKpaJR{009F-A)%0H&dd31R!D(1tAJ%Gr%t)0+2`&oWdD~(Fj55tH2`Xg6O3H>kV=M zAesn>92`qcfD)MxAUQT_6Rnb^zz|TNu|@$qIZmO;cmxw*KmcOCL!J==JnnhS>cv=gYH4a1;TQoA3OJ;-o?Oo9wz^ho3H)@XlJfA!ed{-} zmhN{J1(Yxd)j9QWvY>)qaF_^8R8g|68)*W+W8Y2xYZe3Pnu5i$=$P0X*)e2QAGeW+Uhb{8l*@PYVXkL%`46R zFuS}VcK?mvqleFNE%Tq-Y`sR}WS1%@RCB<45Llt0tsgl~O8<~x0-?Znz!bPGSey!7 znsm4f$jc!#H6Ox23jwIQBKUxg)9e`(2VTL6^jbLpNMy@s6AmB}1)$ovcr3s=uquoV zo8SsXfcHwk8IK^LR-rx&6$SOs6b1#SMuP)yr1fu_Qm}HAl9p^~2m#8(hW&*}E&cm^ zl`l>@fC4w*U1odzzV-kkvgHvGE6``QaGJ|3#5>K{JB()PerJvi!>rC@hFv)!5oaJQ zvqFwM>-L^Wm^UD&N1YkC*Ki+_hW!iuBMvOKz^lg^h4^ZJqbpq}hArr3A4?n*95rEb zU$sK1Xc&h4FcOjZqyCWDBqgdkNp7uDv^m#DRV#(ViPiQ~ZM0Qg%^beOTHozih_$sU Q)E)m9az!{$kp07?P40Ppwg3PC From 0f1c9cb431085e680c31fbfd520bd63cb1ebcf1a Mon Sep 17 00:00:00 2001 From: Nicholas Brazeau Date: Sun, 3 Dec 2023 19:07:47 -0500 Subject: [PATCH 3/4] add grad test --- tests/testthat/test-consistentgradient.R | 67 ++++++++++++++++++++++++ 1 file changed, 67 insertions(+) create mode 100644 tests/testthat/test-consistentgradient.R diff --git a/tests/testthat/test-consistentgradient.R b/tests/testthat/test-consistentgradient.R new file mode 100644 index 0000000..e938ed5 --- /dev/null +++ b/tests/testthat/test-consistentgradient.R @@ -0,0 +1,67 @@ +test_that("model has deterministic results from same start", { + + # sim data + data("IBD_simulation_data", package = "discent") + inputdisc <- IBD_simulation_data + # start params + our_start_params <- rep(0.2, 3) + names(our_start_params) <- 1:3 + our_start_params <- c(our_start_params, "m" = 1e3) + + # given same start parameters and data and iters + # model should always follow the same gradient/trajectory + mod1 <- deme_inbreeding_spcoef(discdat = inputdisc, + start_params = our_start_params, + f_learningrate = 1e-5, + m_learningrate = 1e-1, + b1 = 0.9, + b2 = 0.999, + e = 1e-8, + steps = 1e4, + normalize_geodist = F, + report_progress = F, + return_verbose = F) + + mod2 <- deme_inbreeding_spcoef(discdat = inputdisc, + start_params = our_start_params, + f_learningrate = 1e-5, + m_learningrate = 1e-1, + b1 = 0.9, + b2 = 0.999, + e = 1e-8, + steps = 1e4, + normalize_geodist = F, + report_progress = F, + return_verbose = F) + + mod3 <- deme_inbreeding_spcoef(discdat = inputdisc, + start_params = our_start_params, + f_learningrate = 1e-5, + m_learningrate = 1e-1, + b1 = 0.9, + b2 = 0.999, + e = 1e-8, + steps = 1e4, + normalize_geodist = F, + report_progress = F, + return_verbose = F) + + mod4 <- deme_inbreeding_spcoef(discdat = inputdisc, + start_params = our_start_params, + f_learningrate = 1e-5, + m_learningrate = 1e-1, + b1 = 0.9, + b2 = 0.999, + e = 1e-8, + steps = 5e4, # different number of steps!! + normalize_geodist = F, + report_progress = F, + return_verbose = F) + + testthat::expect_equal(mod1,mod2) + testthat::expect_equal(mod1,mod3) + testthat::expect_equal(mod2,mod2) # trans, so know true, but for complete + testthat::expect_false(identical(mod2,mod4)) # more steps, different results +}) + + From 5240b339be723b876185b1d1fcf3288013339dba Mon Sep 17 00:00:00 2001 From: Nicholas Brazeau Date: Sun, 3 Dec 2023 19:07:55 -0500 Subject: [PATCH 4/4] make adam M --- R/main.R | 2 +- src/main.cpp | 18 +++--- tests/testthat/test-AdamByHand.R | 104 +++++++++++++++---------------- 3 files changed, 62 insertions(+), 62 deletions(-) diff --git a/R/main.R b/R/main.R index f2c58ce..6ff4a2f 100644 --- a/R/main.R +++ b/R/main.R @@ -162,7 +162,7 @@ deme_inbreeding_spcoef <- function(discdat, # catch accidental bad M start if user is standardizing distances if (normalize_geodist & (start_params[names(start_params) == "m"] > 500) ) { warning("You have selected to normalize geographic distances, but your - migration rate parameter is large. Please consider placing it on a + migration rate start parameter is large. Please consider placing it on a similar scale to your normalized geographic distances for stability.") } diff --git a/src/main.cpp b/src/main.cpp index e7bad15..b4bd509 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -184,16 +184,16 @@ Rcpp::List deme_inbreeding_coef_cpp(Rcpp::List args, Rcpp::List args_functions, fi_gradtraj[step][i] = fgrad[i]; } - // // get M moments for Adam - // m1t_m[step] = b1 * m1t_m[step-1] + (1-b1) * mgrad; - // v2t_m[step] = b2 * v2t_m[step-1] + (1-b2) * pow(mgrad, 2); - // m1t_m_hat = m1t_m[step] / (1-pow(b1, step)); - // v2t_m_hat = v2t_m[step] / (1-pow(b2, step)); - // - // // calculate and apply the update for M - // m = m - m_learningrate * (m1t_m_hat / (sqrt(v2t_m_hat) + e)); + // get M moments for Adam + m1t_m[step] = b1 * m1t_m[step-1] + (1-b1) * mgrad; + v2t_m[step] = b2 * v2t_m[step-1] + (1-b2) * pow(mgrad, 2); + m1t_m_hat = m1t_m[step] / (1-pow(b1, step)); + v2t_m_hat = v2t_m[step] / (1-pow(b2, step)); + + // calculate and apply the update for M + m = m - m_learningrate * (m1t_m_hat / (sqrt(v2t_m_hat) + e)); // vanilla GD - m = m - m_learningrate * mgrad; + // m = m - m_learningrate * mgrad; // check bounds on m // will reflect with normal based on magnitude + standard normal sd it is off to proper interval; NB also want to bound m so that it can only explore distance isolation (repulsion versus attraction) if (m < m_lowerbound) { diff --git a/tests/testthat/test-AdamByHand.R b/tests/testthat/test-AdamByHand.R index 9a23b14..249aa91 100644 --- a/tests/testthat/test-AdamByHand.R +++ b/tests/testthat/test-AdamByHand.R @@ -54,55 +54,55 @@ test_that("Fi adam by hand", { }) -# test_that("M adam by hand", { -# # sim data -# data("IBD_simulation_data", package = "discent") -# dat <- IBD_simulation_data -# -# # start params -# our_start_params <- rep(0.2, 3) -# names(our_start_params) <- 1:3 -# our_start_params <- c(our_start_params, "m" = 1e3) -# -# # tidy date and calculate gradient for F1 -# input <- dat %>% -# dplyr::filter(deme1 != deme2) %>% -# dplyr::filter(deme1 == 1 | deme2 == 1) -# -# -# input <- input %>% -# dplyr::mutate(gendist = ifelse(gendist > 0.999, 0.999, -# ifelse(gendist < 0.001, 0.001, -# gendist))) %>% # reasonable bounds on logit -# dplyr::mutate(gendist = logit(gendist)) -# -# # now run model -# inputdisc <- dat %>% -# dplyr::filter(deme1 != deme2) -# ret <- deme_inbreeding_spcoef(discdat = inputdisc, -# start_params = our_start_params, -# f_learningrate = 1e-3, -# m_learningrate = 1e-5, -# b1 = 0.9, -# b2 = 0.999, -# e = 1e-8, -# steps = 1e3, -# normalize_geodist = F, -# report_progress = T, -# return_verbose = T) -# # back out for m adam -# b1 <- 0.9 -# b2 <- 0.999 -# e <- 1e-8 -# m_learningrate <- 1e-5 -# mt_m <- b1 * 0 + (1-b1) * ret$fi_gradtraj[2,1] -# vt_m <- b2 * 0 + (1-b2) * (ret$fi_gradtraj[2,1]^2) -# mt_mhat <- mt_m / (1 - b1^1) -# vt_mhat <- vt_m / (1 - b2^1) -# mnew1 <- ret$m_run[1] - m_learningrate * (mt_mhat/(sqrt(vt_mhat) + e)) -# ret$m_run[1] -# -# # test out -# testthat::expect_equal(mnew1, ret$m_run[1]) -# -# }) +test_that("M adam by hand", { + # sim data + data("IBD_simulation_data", package = "discent") + dat <- IBD_simulation_data + + # start params + our_start_params <- rep(0.2, 3) + names(our_start_params) <- 1:3 + our_start_params <- c(our_start_params, "m" = 1e3) + + # tidy date and calculate gradient for F1 + input <- dat %>% + dplyr::filter(deme1 != deme2) %>% + dplyr::filter(deme1 == 1 | deme2 == 1) + + + input <- input %>% + dplyr::mutate(gendist = ifelse(gendist > 0.999, 0.999, + ifelse(gendist < 0.001, 0.001, + gendist))) %>% # reasonable bounds on logit + dplyr::mutate(gendist = logit(gendist)) + + # now run model + inputdisc <- dat %>% + dplyr::filter(deme1 != deme2) + ret <- deme_inbreeding_spcoef(discdat = inputdisc, + start_params = our_start_params, + f_learningrate = 1e-3, + m_learningrate = 1e-5, + b1 = 0.9, + b2 = 0.999, + e = 1e-8, + steps = 1e3, + normalize_geodist = F, + report_progress = T, + return_verbose = T) + # back out for m adam + b1 <- 0.9 + b2 <- 0.999 + e <- 1e-8 + m_learningrate <- 1e-5 + mt_m <- b1 * 0 + (1-b1) * ret$fi_gradtraj[2,1] + vt_m <- b2 * 0 + (1-b2) * (ret$fi_gradtraj[2,1]^2) + mt_mhat <- mt_m / (1 - b1^1) + vt_mhat <- vt_m / (1 - b2^1) + mnew1 <- ret$m_run[1] - m_learningrate * (mt_mhat/(sqrt(vt_mhat) + e)) + ret$m_run[1] + + # test out + testthat::expect_equal(mnew1, ret$m_run[1]) + +})