Nothing
.do_damped_WLS_spprec <- function(
sXaug, zInfo, # fixed descriptors of system to be solved
old_Vscaled_beta,
oldAPHLs,
APHLs_args,
damping,
dampingfactor=2, ## no need to change the input value
ypos,off,GLMMbool,etaFix,
lambda_est,
wranefblob, seq_n_u_h,
ZAL_scaling, # locally fixed, "resident"; only changed in return value
processed,
Trace,
phi_est, #H_global_scale,
n_u_h,
ZAL, # unscaled one
which_LevMar_step,
update_sXaug_constant_arglist,
# promise rather than argument:
low_pot=NULL,
v_infer_args=NULL, # not null for beta optimization with v_in_b optimization i.e. in *some* .do_damped_WLS_outer() call
stylefn, # in-loop stylefn for damped_WLS
stylefn_v_out= .spaMM.data$options$stylefns$v_out_last,
stylefn_v_in= .spaMM.data$options$stylefns$v_in_last, ##
outer,
IRLS_fn =get(".solve_v_h_IRLS_spprec", asNamespace("spaMM"), inherits=FALSE)
) {
ZAL_scaling <- 1 ## TAG: scaling for spprec
if (outer) {
trace <- max(0L,Trace-1L)
stylefn_v <- stylefn_v_out
} else {
trace <- max(0L,Trace-2L)
stylefn_v <- stylefn_v_in
}
if (processed$p_v_obj=="p_v" && which_LevMar_step!="v") {
objname <- "p_v"
} else {
objname <- APHLs_args$which <- "hlik"
}
oldlik <- oldAPHLs[[objname]]
initdamping <- damping
gainratio_grad <- zInfo$gainratio_grad
restarted_at_e_7 <- FALSE
first_it <- TRUE
prev_gainratio <- -Inf
if (Trace && ! is.null(v_infer_args)) {
cat(stylefn("[")) # cat(which_LevMar_step) #=> a substep of V_IN_B: "strict_v|b" or "b_&_v_in_b"
}
GLGLLM_const_w <- attr(processed$models,"GLGLLM_const_w")
while ( TRUE ) { ## loop on damping; each iteration produce blue + ul-greens + yellow
# if (trace && ! is.null(v_infer_args)) {
# cat(c(damping)) # => what is shown after [.v_h iter...]V_h IRLS...; .dampingfactor=* damping=
# }
if (processed$HL[1L]==1L) { ## ML fit
Vscaled_beta <- old_Vscaled_beta
LM_z <- zInfo["scaled_grad"] # still a list, but for clarity, emphasizes that only this element is needed
## maximize p_v wrt beta only
if ( which_LevMar_step=="v_b") { ## note tests on summand too !!!!!!!
LevMarblob <- get_from_MME(sXaug=sXaug, which="LevMar_step", LM_z=LM_z, damping=damping)
Vscaled_beta <- list(
v_h=Vscaled_beta$v_h + LevMarblob$dVscaled, ##
beta_eta=Vscaled_beta$beta_eta + LevMarblob$dbeta_eta
)
} else if ( which_LevMar_step=="v") { ## v_h estimation given beta
LevMarblob <- get_from_MME(sXaug=sXaug, which="LevMar_step_v_h", LM_z=LM_z, damping=damping)
Vscaled_beta$v_h <- Vscaled_beta$v_h + LevMarblob$dVscaled
} else if ( which_LevMar_step=="strict_v|b") {
# =: the case where we only fit v_h for the input beta_eta:
# only v_h in Vscaled_beta will be changed, by v_infer_args step below,
LevMarblob <- v_infer_args$LevMarblob ## LevMarblob$d... will be overwritten below
# is.null(LevMarblob) occurs in initial damping=Inf call.
} else if ( which_LevMar_step=="b_&_v_in_b") {
LevMarblob <- get_from_MME(sXaug=sXaug, which="LevMar_step", LM_z=LM_z, damping=damping) # template
dbeta_LevMarblob <- get_from_MME(sXaug=sXaug, which="LevMar_step_beta", LM_z=LM_z, damping=damping)
Vscaled_beta$beta_eta <- Vscaled_beta$beta_eta + dbeta_LevMarblob$dbeta_eta
# v_h in Vscaled_beta will be changed by v_infer_args step below.
v_infer_args$LevMarblob <- get_from_MME(sXaug=sXaug, which="LevMar_step_v_h", LM_z=LM_z, damping=damping)
} else if ( which_LevMar_step=="b_from_v_b") {
LevMarblob <- get_from_MME(sXaug=sXaug, which="LevMar_step", LM_z=LM_z, damping=damping)
Vscaled_beta$beta_eta <- Vscaled_beta$beta_eta + dbeta_LevMarblob$dbeta_eta
} else if ( which_LevMar_step=="b") { ## probably not used; get_from_MME() includes protection for pforpv=0
LevMarblob <- get_from_MME(sXaug=sXaug, which="LevMar_step_beta", LM_z=LM_z, damping=damping)
Vscaled_beta$beta_eta <- Vscaled_beta$beta_eta + dbeta_LevMarblob$dbeta_eta
}
if ( which_LevMar_step!="v" && ! is.null(v_infer_args)) {
if (Trace) cat(stylefn_v("[")) # blue '[' (within red or yellow '[')
v_h_blob <- .wrap_wrap_v_h_IRLS(IRLS_fn=IRLS_fn, v_h=Vscaled_beta$v_h, beta_eta=Vscaled_beta$beta_eta,
seq_n_u_h, GLMMbool, wranefblob, processed, lambda_est,
v_infer_args, Trace) # each underline green is a damping _loop_ not a damping step
if (Trace) cat(stylefn_v("]"))
if (trace) .cat_break_info( v_h_blob, stylefn_v, stylefn) ## prints, at the level of the outer damped_WLS, the results of the v_h IRLS
Vscaled_beta$v_h <- v_h_blob$v_h
LevMarblob$dVscaled_beta <- list(
v_h=Vscaled_beta$v_h - old_Vscaled_beta$v_h, ##
beta_eta=Vscaled_beta$beta_eta -old_Vscaled_beta$beta_eta
)
}
} else { ## joint hlik maximization
LevMarblob <- get_from_MME(sXaug=sXaug, which="LevMar_step", LM_z=zInfo["scaled_grad"], damping=damping)
Vscaled_beta <- list(
v_h=old_Vscaled_beta$v_h + LevMarblob$dVscaled,
beta_eta=old_Vscaled_beta$beta_eta + LevMarblob$dbeta_eta
)
}
v_h <- Vscaled_beta$v_h #### * ZAL_scaling =1 here
eta <- off + drop(processed$AUGI0_ZX$X.pv %*% Vscaled_beta$beta_eta) + drop(ZAL %id*% v_h) ## length nobs
newmuetablob <- .muetafn(eta=eta,BinomialDen=processed$BinomialDen,processed=processed, phi_est=phi_est)
neww.resid <- .calc_w_resid(newmuetablob$GLMweights,phi_est, obsInfo=processed$how$obsInfo)
newH_w.resid <- .calc_H_w.resid(neww.resid, muetablob=newmuetablob, processed=processed) # for LLF w.resid is not generally defined.
if (is.null(etaFix$v_h)) {
if (GLMMbool) {
u_h <- v_h
newwranefblob <- wranefblob ## keep input wranefblob since GLMM and lambda_est not changed
} else {
u_h <- processed$u_h_v_h_from_v_h(v_h)
if ( ! is.null(maybe <- attr(u_h,"v_h"))) v_h <- maybe
## update functions u_h,v_h
newwranefblob <- processed$updateW_ranefS(u_h=u_h,v_h=v_h, lambda=lambda_est)
}
} else newwranefblob <- wranefblob
sXaug_arglist <- c(update_sXaug_constant_arglist,
list(w.ranef=newwranefblob$w.ranef,
#weight_X=newweight_X,
H_w.resid=newH_w.resid))
if ( ! GLMMbool) {newZAL_scaling <- 1} ## TAG: scaling for spprec
####
APHLs_args$dvdu <- newwranefblob$dvdu
APHLs_args$u_h <- u_h
APHLs_args$muetablob <- newmuetablob
#
if (processed$p_v_obj=="p_v" && which_LevMar_step!="v") { ## new damping -> new weights -> new expensive computation to evaluate p_v
if (GLGLLM_const_w) {
newsXaug <- NULL
APHLs_args$sXaug <- sXaug
} else {
newsXaug <- do.call(processed$spprec_method, # ie, def_AUGI0_ZX_spprec
sXaug_arglist)
if (Trace) {
tracechar <- ifelse(.BLOB(newsXaug)$nonSPD,"!",".")
cat(stylefn(tracechar)) # yellow in V_IN_B case
}
APHLs_args$sXaug <- newsXaug
}
} else { ## sXaug_arglist will still be used after the loop !!!!!!!!!!!!!!!!!!!!
APHLs_args$sXaug <- newsXaug <- NULL # to compute hlik, no expensive matrix computation.
}
newAPHLs <- do.call(".calc_APHLs_from_ZX", APHLs_args)
#print(c(unlist(newAPHLs)))
newlik <- unlist(newAPHLs[objname]) # keep name
#
if (damping==0L) {
breakcond <- "damping=0"
break
}
if ( which_LevMar_step=="strict_v|b") {
breakcond <- "v|b_no_loop"
attr(breakcond,"v_pot4improv") <- .pot4improv("v", sXaug, gainratio_grad=zInfo$gainratio_grad, seq_n_u_h)
break
} # =: single call to .calc_APHLs_from_ZX to only fit v_h for the input beta_eta.
if (first_it) {
pot4improv <- .pot4improv(which_LevMar_step, sXaug, gainratio_grad=zInfo$m_grad_obj, seq_n_u_h)
loc_pot_tol <- .loc_pot_tol(which_LevMar_step, processed$spaMM_tol)
if (is.null(low_pot)) low_pot <- (pot4improv < loc_pot_tol)
}
if (low_pot) { # keeping the low_pot condition may be important for the "605" tests. We may always suppress it by the spaMM.options.
very_low_pot <- (pot4improv < loc_pot_tol/10)
breakcond <- "low_pot"
attr(breakcond, "pot4improv") <- pot4improv
attr(breakcond, "very_low_pot") <- very_low_pot
if (processed$HL[1L]==1L && which_LevMar_step=="v_b") {
v_pot4improv <- .pot4improv("v", sXaug, gainratio_grad=zInfo$m_grad_obj, seq_n_u_h)
attr(breakcond, "no_overfit") <- ((v_pot4improv < processed$spaMM_tol$v_pot_tol))
}
break
}
## ELSE
gainratio <- (newlik!=-Inf) ## -Inf occurred in binary probit with extreme eta...
if (gainratio) {
summand <- .calc_summand_gainratio_spprec(processed, which_LevMar_step, LevMarblob, seq_n_u_h, ZAL_scaling=1, gainratio_grad)
denomGainratio <- sum(summand)
if (denomGainratio<=0) {
## In the summand, all elements of dbeta*( LevenbergMstep_result$dampDpD * dbeta) should be positive.
# Inner product dbeta . LevenbergMstep_result$rhs should be positive too
## However, numerical error may lead to <0 or even -Inf
## Further, if there are both -Inf and +Inf elements the sum is NaN.
denomGainratio <- denomGainratio-sum(summand[summand<0]) # quite rough patch
}
dlogL <- newlik-oldlik
conv_logL <- abs(dlogL)/(1+abs(newlik))
gainratio <- 2*dlogL/denomGainratio ## cf computation in MadsenNT04 below 3.14, but generalized for D' != I ; numerator is reversed for maximization
# but if gradient is practically zero and damping ~0 we may not wish to compare ~0 to ~0...
# which is why we break, rather than stop, on (damping>1e100). MadsenNT04 have other stopping crits
} else { ## 2017/10/16 a patch to prevent a stop in this case, but covers up dubious computations (FIXME)
newlik <- -.Machine$double.xmax
dlogL <- newlik-oldlik
conv_logL <- abs(dlogL)/(1+abs(newlik))
denomGainratio <- Inf # bc denomGainratio can be tested below
}
if (trace) {
cat(stylefn(paste(" dampingfactor=",dampingfactor,#"innerj=",innerj,
"damping=",damping,"gainratio=",gainratio,"oldlik=",oldlik,"newlik=",newlik)))
}
if (is.nan(gainratio)) {
# if the initial logL is the solution logL, then damping diverges
# it is then possible that some element of dVscaled_beta =0 and some of dampDpD =Inf
# then the summand has some NaN
# At the same time not all elements of dVscaled_beta need be 0 (eg small changes in eta for mu=0 or 1 in binomial models)
# so testing dVscaled_beta is not sufficient to stop the algo
# (LevenbergM is quite slow in such cases)
breakcond <- "NaN_gain"
break
}
if (objname == "p_v") { # then the starting objective may be 'too high' and we need to handle that
div_gainratio <- (gainratio-prev_gainratio)*damping/dampingfactor
if (div_gainratio < -max(LevMarblob$dampDpD)/(1e05*damping)) {
breakcond <- "div_gain" # used to switch from "V_IN_B" to "strict_v|b"
break
}
}
if (gainratio > 0) { # gainratio may be always negative if initial ranefs better optimize logL than the correct solution does.
## cf Madsen-Nielsen-Tingleff again, and as in levmar library by Lourakis
damping <- damping * max(1/3,1-(2*gainratio-1)^3) # gainratio->0 factor-> 2; gainratio->1 factor->1/3
breakcond <- "OK_gain"
break
}
if (
conditionNS_for_restart <- (
initdamping>1e-7 &&
( ! restarted_at_e_7 ) &&
denomGainratio<loc_pot_tol/10
)
) { # then restart
damping <- 1e-7
dampingfactor <- 2
restarted_at_e_7 <- TRUE
prev_gainratio <- -Inf
if (trace) cat("-")
# and continue
} else if ( restarted_at_e_7 && ## has gone through 1e-7, # handles the cases where we started with too high damping
(( odd_condition <- (dampingfactor>4 && ## ie at least 2 iteration of the while() => prev_conv_logL is available
conv_logL <1e-8 && abs(prev_conv_logL) <1e-8)) || # possible reason for odd_condition is truely stuck because of "COMPoisson" reasons?
damping>initdamping)
) { # restarted && damping>init must be a sufficient condition for break
breakcond <- "stuck_obj"
break ## cases were we do not expect any significant improvement
} else { ## other UNsuccessful step
prev_gainratio <- gainratio
prev_conv_logL <- conv_logL
damping <- damping*dampingfactor
dampingfactor <- dampingfactor*2
if (damping>1e10) {
breakcond <- "div_damp"
break
}
}
first_it <- FALSE # : skipped if break in first iteration
} ################# end while(TRUE)
if (Trace && ! is.null(v_infer_args)) {
cat(stylefn("]"))
if (trace) {cat(stylefn(damping))}
}
if (trace) cat(breakcond)
if (is.null(newsXaug)) { ## which means that hlik is the local objective.
# For HL11, p_v will be used as oldAPHLs in the next call to .do_damped_WLS_outer() in an alternating algo;
# and sXaug may be needed to compute sscaled in .solve_v_h_IRLS()
# For PQL fits newsXaug has not been needed in the damping loop but will be needed after exiting this fn
# (e.g., for its next call -> LevMarblob <- get_from_MME(sXaug=sXaug, which="LevMar_step", LMrhs=zInfo$scaled_grad, damping=damping))
if (GLGLLM_const_w) {
APHLs_args$sXaug <- newsXaug <- sXaug
} else {
newsXaug <- do.call(processed$spprec_method, sXaug_arglist)
if (Trace) {
tracechar <- ifelse(.BLOB(newsXaug)$nonSPD,"!",".")
if (processed$p_v_obj=="p_v") { # v estimation within HL11
cat(stylefn_v(tracechar))
} else cat(stylefn(tracechar)) # PQL/L, vb extimation
}
APHLs_args$sXaug <- newsXaug
}
APHLs_args$which <- processed$p_v_obj # "p_v" #
newAPHLs <- do.call(".calc_APHLs_from_ZX", APHLs_args)
}
RESU <- list(lik=newlik, APHLs=newAPHLs, damping=damping, sXaug=newsXaug,
# fitted=fitted, ## FIXME: removed so that no shortcut a la Bates in calc_APHLs_from_ZX; reimplement the shorcut in that fn?
eta=newmuetablob$sane_eta, muetablob=newmuetablob, wranefblob=newwranefblob,
breakcond=breakcond,
v_h=v_h, u_h=u_h, w.resid=neww.resid) # newweight_X does not exists and not needed in spprec
# if ( ! first_it) { # if not break in first iteration
# RESU$conv_logL_not_first_it <- conv_logL
# }
if ( ! GLMMbool ) {
RESU$ZAL_scaling <- newZAL_scaling
# RESU$Xscal <- newXscal ## does not exist and presumably not needed.
Vscaled_beta$v_h <- v_h/newZAL_scaling ## represent solution in new scaling...
}
RESU$Vscaled_beta <- Vscaled_beta
return(RESU)
}
#copies to allow independent debug()ing
.do_damped_WLS_v_in_b_spprec <- .do_damped_WLS_spprec
.do_damped_WLS_outer_spprec <- .do_damped_WLS_spprec
# processed, Trace
.updates_from_w.resid <- function(RESU, phi_est, update_sXaug_constant_arglist,
wranefblob, processed, Trace, stylefn) {
RESU$w.resid <- .calc_w_resid(RESU$muetablob$GLMweights,phi_est, obsInfo=processed$how$obsInfo)
sXaug_arglist <- c(update_sXaug_constant_arglist, # contained H_global_scale but not longer so
list(w.ranef=wranefblob$w.ranef))
sXaug_arglist$H_w.resid <- .calc_H_w.resid(RESU$w.resid, muetablob=RESU$muetablob, processed=processed)
RESU$sXaug <- do.call(processed$spprec_method, # ie, def_AUGI0_ZX_spprec
sXaug_arglist)
if (Trace) {
tracechar <- ifelse(.BLOB(RESU$sXaug)$nonSPD,"!",".")
cat(stylefn(tracechar)) # yellow in V_IN_B case
}
RESU
}
.WLS_substitute_spprec <- function(update_sXaug_constant_arglist, Vscaled_beta, off, etaFix, mod_attr,
lambda_est, ZAL,
processed, phi_est,
wranefblob, Trace,stylefn) {
# Vscaled_beta must have been provided by somethin else than damped_WLS_blob
# drop, not as.vector(): names are then those of (final) eta and mu -> used by predict() when no new data
eta <- off + drop(processed$AUGI0_ZX$X.pv %*% Vscaled_beta$beta_eta) + drop(ZAL %id*% Vscaled_beta$v_h)
RESU <- list()
if (is.null(etaFix$v_h)) {
v_h <- Vscaled_beta$v_h ## * ZAL_scaling (=1)
if ( mod_attr$GLMMbool ) {
RESU$u_h <- RESU$v_h <- v_h ## keep input wranefblob since lambda_est not changed
} else {
RESU$u_h <- u_h <- processed$u_h_v_h_from_v_h(v_h)
if ( ! is.null(maybe <- attr(u_h,"v_h"))) v_h <- maybe
RESU$v_h <- v_h
## update functions u_h,v_h
RESU$wranefblob <- wranefblob <- processed$updateW_ranefS(u_h=u_h,v_h=v_h, lambda=lambda_est)
#if ( ! mod_attr$GLMMbool) {
RESU$ZAL_scaling <- 1 ## TAG: scaling for spprec
# }
}
}
RESU$muetablob <- .muetafn(eta=eta,BinomialDen=processed$BinomialDen,processed=processed, phi_est=phi_est)
if ( (! mod_attr$LLM_const_w) && (! mod_attr$GLGLLM_const_w) ) {
RESU <- .updates_from_w.resid(RESU, phi_est, update_sXaug_constant_arglist,
wranefblob, # do NOT try to use RESU$wranefblob bc it does not necessarily exists (if no updated)
processed, Trace, stylefn)
} ## else no need to update RESU$'s w.resid, sXaug, H_w.resid
return(RESU) ## contains only updated quantities
}
.solve_IRLS_as_spprec <-
function(
ZAL, y=processed$y,
n_u_h=length(u_h),
#H_global_scale,
lambda_est, muetablob=NULL, off=processed$off, maxit.mean, etaFix,
wranefblob, processed,
## for ! LMM
phi_est,
## supplement for LevenbergM
beta_eta,
## supplement for ! GLMM
u_h, v_h, w.resid=NULL,
H_w.resid,
for_intervals,
##
corrPars, # corrPars needed together with adjMatrix to define Qmat
verbose=processed$verbose,
LevM_HL11_method=.spaMM.data$options$LevM_HL11_method,
## ignored but for consistency of arguments with .solve_IRLS_as_ZX:
H_global_scale
) {
trace <- verbose["TRACE"]
if (trace) {
cat(">")
if (verbose["trace"]) cat(.pretty_summ_lambda(lambda_est,processed))
}
pforpv <- ncol(processed$AUGI0_ZX$X.pv)
nobs <- length(y)
seq_n_u_h <- seq_len(n_u_h)
ypos <- n_u_h+seq_len(nobs)
lcrandfamfam <- attr(processed$rand.families,"lcrandfamfam")
mod_attr <- attributes(processed[["models"]])
LMMbool <- mod_attr$LMMbool
GLMMbool <- mod_attr$GLMMbool
LevenbergM <- (processed$LevenbergM["LM_start"] && is.null(for_intervals))
is_HL1_1 <- (processed$HL[1L]==1L)
fpot_cond <- processed$spaMM_tol$fpot_cond
if (fpot_cond) fpot_tol <- processed$spaMM_tol$fpot_tol
if ( is.null(for_intervals) && is_HL1_1) {
if (pforpv==0L) { # outer beta
which_LevMar_step <- "v" # appeared with outer beta estim
v_iter <- 0L
} else which_LevMar_step <- default_b_step <- LevM_HL11_method[["b_step"]]
rescue_thr <- processed$spaMM_tol$rescue_thr
rescue_nbr <- 0L
prev1_rescued <- FALSE
} else which_LevMar_step <- "v_b"
old_relV_beta <- NULL
not_moving <- FALSE
damped_WLS_blob <- NULL
d_relV_b_tol <- processed$spaMM_tol$d_relV_b_tol
d_relV_b_tol_LM <- processed$spaMM_tol$d_relV_b_tol_LM
if ( LevenbergM) {
dampings_env <- list2env(.spaMM.data$options$spaMM_tol$dampings_env_v)
}
if ( ! LMMbool) {
checkpot_min_it <- as.integer(maxit.mean/4L) # reconsider all usage of this and possibly extend to spprec (see ref to pot4improv in test-mv-extra for a test)
constant_zAug_args <- list(n_u_h=n_u_h, nobs=nobs, pforpv=pforpv, y=y, off=off, ZAL=ZAL, processed=processed)
}
##### initial sXaug
ZAL_scaling <- 1 ## TAG: scaling for spprec
if (is.null(muetablob)) { ## NULL input eta allows NULL input muetablob
eta <- off + drop(processed$AUGI0_ZX$X.pv %*% beta_eta) + drop(ZAL %id*% v_h)
muetablob <- .muetafn(eta=eta,BinomialDen=processed$BinomialDen,processed=processed, phi_est=phi_est)
}
## varies within loop if ! LMM since at least the GLMweights in w.resid change
if ( is.null(w.resid) ) w.resid <- .calc_w_resid(muetablob$GLMweights,phi_est, obsInfo=processed$how$obsInfo)
## needs adjMatrix and corrPars to define Qmat
update_sXaug_constant_arglist <- list(AUGI0_ZX=processed$AUGI0_ZX, corrPars=corrPars,
cum_n_u_h=processed$cum_n_u_h #,H_global_scale=H_global_scale
)
sXaug_arglist <- c(update_sXaug_constant_arglist,
list(w.ranef=wranefblob$w.ranef))
sXaug_arglist$H_w.resid <- .calc_H_w.resid(w.resid, muetablob=muetablob, processed=processed)
sXaug <- do.call(processed$spprec_method, # ie, def_AUGI0_ZX_spprec
sXaug_arglist)
if (trace) {
stylefn <- switch(which_LevMar_step,
v=.spaMM.data$options$stylefns$vloop,
V_IN_B=.spaMM.data$options$stylefns$v_in_loop,
.spaMM.data$options$stylefns$betaloop )
if (LevenbergM) cat("LM")
tracechar <- ifelse(.BLOB(sXaug)$nonSPD,"!",".")
cat(stylefn(tracechar)) # yellow in V_IN_B case
}
if ( ! is.null(for_intervals)) {
Vscaled_beta <- list(v_h=v_h/ZAL_scaling, beta_eta=for_intervals$beta_eta)
fixefobjfn <- names(for_intervals$fixeflik)
} else {
Vscaled_beta <- list(v_h=v_h/ZAL_scaling,beta_eta=beta_eta)
}
# to be evaluated once when it becomes needed:
delayedAssign("constant_v_infer_args", list( # ultimately for the .solve_v_h_IRLS_spprec() call
X.pv=processed$AUGI0_ZX$X.pv,
ZAL=ZAL, y=y, n_u_h=n_u_h, #H_global_scale=H_global_scale,
lambda_est=lambda_est, off=off,maxit.mean=maxit.mean,etaFix=etaFix,
processed=processed, phi_est=phi_est,
trace=trace, corrPars=corrPars, dampings_env=dampings_env))
## Loop controls:
allow_LM_restart <- ( ! LMMbool && ! LevenbergM && is.null(for_intervals) && is.na(processed$LevenbergM["user_LM"]) )
if (allow_LM_restart) {
keep_init <- new.env()
names_keep <- c("sXaug","wranefblob","muetablob","u_h","H_w.resid","w.resid","v_h","beta_eta","Vscaled_beta")
for (st in names_keep) keep_init[[st]] <- environment()[[st]]
}
LMcond <- - 10. # also for hlik LM algo
best_HL1_lik <- -Inf
pot4improv <- NULL
################ L O O P ##############
for (innerj in 1:maxit.mean) {
# if (sanitize <- FALSE) {
# w_sane_etamo <- (sXaug$BLOB$WLS_mat_weights)*(muetablob$sane_eta-off)
# rhs_sanitized_v_b <- c(v_h*attr(sXaug,"w.ranef") + .crossprod(ZAL, w_sane_etamo),
# .crossprod(sXaug$AUGI0_ZX$X.pv, w_sane_etamo))
#
# insane_Vscaled_beta <- Vscaled_beta
# Vscaled_beta <- get_from_MME(sXaug, szAug =list(m_grad_obj=rhs_sanitized_v_b))
# Vscaled_beta <- list(v_h=Vscaled_beta$dv_h,beta_eta=Vscaled_beta$dbeta_eta)
# }
if( ! LevenbergM && allow_LM_restart) { ## FIXME the next step improvement would be
# ./. to keep track of lowest lambda that created problem and use LM by default then
if (innerj>3L) {
crit <- abs_d_relV_beta/(old_abs_d_relV_beta+1e-8)
LMcond <- LMcond + mean(sqrt(crit))^2
## cat(mean(abs_d_relV_beta/old_abs_d_relV_beta)," ")
# cat(LMcond/innerj," ")
if (LMcond/innerj>0.5 ||
(innerj> checkpot_min_it && pot4improv > max(10,old_pot4improv))
) {
if (trace) cat("!LM") # ie, LevenbergM!
for (st in names_keep) assign(st,keep_init[[st]])
LevenbergM <- TRUE
# Vscaled_beta included in keep_init hence assigned by assign(st,keep_init[[st]])
dampings_env <- list2env(.spaMM.data$options$spaMM_tol$dampings_env_v)
damped_WLS_blob <- NULL
allow_LM_restart <- FALSE
if ( which_LevMar_step=="v_b") {
## The LevM.negbin test finds "strict_v|b" poorer than "V_IN_B"" (note some divergent p_v's) which led to:
# which_LevMar_step <- "V_IN_B" ## not modified by if (... ! is.null(damped_WLS_blob) ...) before being used.
# but the optim_LevM's (update(br$fullfit,fixed=... test shows one should keep using "v_b" here.
# otherwise "!LM" differs (and is poorer) from LevM=TRUE (which indeed starts from "v_b")
# However, it's not clear why "V_IN_B" is poorer (and it's not the step on which the loop terminates) => (there is a fixme on the fit_as_ZX version...)
strictv_parent_info <- c(from=which_LevMar_step, breakcond="")
rescue_thr <- processed$spaMM_tol$rescue_thr
rescue_nbr <- 0L
prev1_rescued <- FALSE
}
}
}
if (innerj>2L) {
old_abs_d_relV_beta <- abs_d_relV_beta
old_pot4improv <- pot4improv
}
}
##### get the lik of the current state
if ( ! is.null(for_intervals)) {
loc_logLik_args <- list(sXaug=sXaug, processed=processed, phi_est=phi_est,
lambda_est=lambda_est, dvdu=wranefblob$dvdu, u_h=u_h, muetablob=muetablob)
oldlik <- unlist(do.call(".calc_APHLs_from_ZX",loc_logLik_args)[fixefobjfn]) # unlist keeps name
} else if (LevenbergM) { ## then logL is necessary to check for increase
if (is.null(damped_WLS_blob)) {
oldAPHLs <- .calc_APHLs_from_ZX(sXaug=sXaug, processed=processed, phi_est=phi_est, which=processed$p_v_obj,
lambda_est=lambda_est, dvdu=wranefblob$dvdu, u_h=u_h,
muetablob=muetablob)
} else { ## Levenberg and innerj>1
oldAPHLs <- damped_WLS_blob$APHLs
}
}
#####
##### RHS
if (LMMbool) {
zInfo <- list(z2=NULL,z1=y-off,sscaled=0)
y_eta_ <- y-muetablob$sane_eta
## the gradient for -p_v (or -h), independent of the scaling
zInfo$m_grad_obj <- .calc_m_grad_obj(zInfo,
z1_eta=y_eta_, # same as zInfo$z1 - etamo = y-off-(eta-off)
z1_sscaled_eta=y_eta_ , # same as zInfo$z1_sscaled - etamo =(sscaled=0)= zInfo$z1 - etamo=y-off-(eta-off)
GLMMbool, v_h, wranefblob,
H_w.resid=.BLOB(sXaug)$H_w.resid,
ZAL, X.pv=processed$AUGI0_ZX$X.pv)
} else {
if ( ! GLMMbool) {
# # arguments for init_resp_z_corrections_new called in calc_zAug_not_LMM
# init_z_args <- c(constant_init_z_args,
# list(w.ranef=wranefblob$w.ranef, u_h=u_h, v_h=v_h, dvdu=wranefblob$dvdu,
# sXaug=sXaug)) # H_w.resid provided by sXaug!
# z2 <- do.call(".init_resp_z_corrections_new",init_z_args)$z20
z2 <- .calc_z2(lcrandfamfam=lcrandfamfam, psi_M=processed$psi_M, cum_n_u_h=processed$cum_n_u_h, rand.families=processed$rand.families,
u_h=u_h, lambda_est=lambda_est, v_h=v_h, dvdu=wranefblob$dvdu)
} else z2 <- rep(0,n_u_h)
calc_zAug_args <- c(constant_zAug_args,
list(muetablob=muetablob, dlogWran_dv_h=wranefblob$dlogWran_dv_h,
sXaug=sXaug,
w.ranef=wranefblob$w.ranef,
w.resid=w.resid,
z2=z2) )
zInfo <- do.call(".calc_zAug_not_LMM",calc_zAug_args)
if (GLMMbool) zInfo$z2 <- NULL
etamo <- muetablob$sane_eta - off
## the gradient for -p_v (or -h), independent of the scaling
zInfo$m_grad_obj <- .calc_m_grad_obj(zInfo, z1_eta=zInfo$z1-etamo, z1_sscaled_eta=zInfo$z1_sscaled - etamo , GLMMbool, v_h, wranefblob,
H_w.resid=.BLOB(sXaug)$H_w.resid,
ZAL, X.pv=processed$AUGI0_ZX$X.pv)
}
## keep name 'w'zAug to emphasize the distinct weightings of zaug and Xaug (should have been so everywhere)
#####
##### improved Vscaled_beta
## he solver uses a 'beta first" approach to solving for d_beta and d_v... even for LMM
if (LevenbergM) {
zInfo$gainratio_grad <- zInfo$m_grad_obj
zInfo$scaled_grad <- zInfo$m_grad_obj
}
#
if ( ! is.null(for_intervals)) {
currentDy <- (for_intervals$fixeflik-oldlik)
if (LMMbool) etamo <-muetablob$sane_eta - off
intervalBlob <- .intervalStep_spprec(old_v_h_beta=Vscaled_beta,
sXaug=sXaug,zInfo=zInfo,
for_intervals=for_intervals,
currentlik=oldlik,currentDy=currentDy,
ZAL=ZAL, etamo=etamo)
damped_WLS_blob <- NULL
Vscaled_beta <- intervalBlob$v_h_beta
} else if (LevenbergM) {
if (trace>1L) {
stylefn <- switch(which_LevMar_step,
v=.spaMM.data$options$stylefns$vloop,
V_IN_B=.spaMM.data$options$stylefns$v_in_loop,
.spaMM.data$options$stylefns$betaloop )
if (pforpv) {
maxs_grad <- c(max(abs(zInfo$m_grad_obj[seq_n_u_h])),max(abs(zInfo$m_grad_obj[-seq_n_u_h])))
} else maxs_grad <- c(max(abs(zInfo$m_grad_obj[seq_n_u_h])),0) # outer beta
cat(stylefn("iter=",innerj,", max(|grad|): v=",maxs_grad[1L],"beta=",maxs_grad[2L],";\n"))
}
constant_APHLs_args <- list(processed=processed, which=processed$p_v_obj, sXaug=sXaug, phi_est=phi_est, lambda_est=lambda_est)
# the following block needs m_grad_v the new m_grad_v hence its position
if (is_HL1_1 && ! is.null(damped_WLS_blob)) {
#### Get next LM step && conditionally update old_relV_beta ####
# we assess convergence at the end of the loop by comparing old_relV_beta to relV_beta. We update old_relV_beta
# (1) here, in all cases where v has been updated;
# or (2) in one reversal case, it needs to be updated after the test at the end of the loop.
if (just_rescued <- identical(attr(damped_WLS_blob, "step"), "rescue")) {
rescue_nbr <- rescue_nbr + 1L
old_relV_beta <- relV_beta
if (pforpv) {
if (prev1_rescued || rescue_nbr > rescue_thr["V_IN_B"]) {
which_LevMar_step <- "V_IN_B" # # [yellow [cyan .cyan .....ul green...
} else which_LevMar_step <- "v_b" # outer beta
}
} else if (which_LevMar_step=="v_b" || which_LevMar_step=="b_from_v_b" ) {
if (rescue_nbr > rescue_thr["strictv"] && #rescue has been previously needed in the outer loop
damped_WLS_blob$breakcond != "low_pot" ## in particular, if OK_gain, we play safer if we know a problem occurred previously
## "stuck_obj" occurs here too...
) {
strictv_parent_info <- c(from=which_LevMar_step, breakcond=damped_WLS_blob$breakcond)
which_LevMar_step <- "strict_v|b"
# } else if (max(abs(m_grad_v)) > max(abs(old_m_grad_v))) which_LevMar_step <- "v" # test is fausse bonne idee...
} else {
v_parent_info <- c(from=which_LevMar_step, breakcond=damped_WLS_blob$breakcond)
# Tried this from [v3.12.34 -> 54] with poor effect on test-nloptr #362... (spcorr)
# if (damped_WLS_blob$breakcond=="stuck_obj") {
# which_LevMar_step <- "V_IN_B"
# } else
which_LevMar_step <- "v" # standard switch from yellow to underlined cyan. It is generally not a good idea to switch immediately to "V_IN_B"
v_iter <- 0L
}
} else if (which_LevMar_step=="v") {
if (damped_WLS_blob$breakcond == "low_pot") { ## LevMar apparently maximized h wrt v after several iterations
#cat(damped_WLS_blob$breakcond)
old_relV_beta <- relV_beta ## serves to assess convergence !!! which is thus dependent on condition ( hlik_stuck || ! need_v_step)
if (pforpv) which_LevMar_step <- default_b_step # We should not reach this line when RHS is "v_in_b" # ! outer beta
} else {
## v_h estimation not yet converged, continue with it
# But this means that "stuck_obj" on a "v" step must have dealt with elsewhere, otherwise we are looping on stuck_obj
# => special rescue case.
v_iter <- v_iter+1L
if (v_iter>10L && pforpv) which_LevMar_step <- "V_IN_B"
}
} else if (which_LevMar_step=="strict_v|b") {
old_relV_beta <- relV_beta
if (strictv_parent_info[["from"]] %in% c("v_b","b_from_v_b" )
&& strictv_parent_info[["breakcond"]] !="low_pot") { # (which implies that rescue was not just called)
which_LevMar_step <- "v_b"
} else { # strictv was called after V_IN_B (!)
if (rescue_nbr > rescue_thr["re_V_IN_B"]) {
which_LevMar_step <- "V_IN_B"
} else which_LevMar_step <- "v_b"
}
} else if (which_LevMar_step=="V_IN_B") {
## .wrap_do_damped_WLS_outer(which_LevMar_step=="V_IN_B") has run .do_damped_WLS_spprec(which_LevMar_step=="b_&_v_in_b") ie.
## a damping loop of {a beta updating followed by a v_h_IRLS}'s. We assess the resulting breakcond.
breakcond <- damped_WLS_blob$breakcond
if (breakcond=="stuck_obj" || breakcond=="div_gain") {
strictv_parent_info <- c(from=which_LevMar_step,breakcond=breakcond)
which_LevMar_step <- "strict_v|b" # the call to "strict_v|b" may seem odd but results in clean optim
#If we did that in the wrap... then we would next compare two identical "strict_v|b"
} else { # a good breakcond MAY result in switching back to v_b
old_relV_beta <- relV_beta
if (rescue_nbr > rescue_thr["re_V_IN_B"]) {
which_LevMar_step <- "V_IN_B"
} else which_LevMar_step <- "v_b"
}
} else if (default_b_step=="v_in_b") { # presumably not used
old_relV_beta <- relV_beta
} else { ## "b" or any unanticipated case # presumably not used
warning("Unexpected case in .solve_IRLS_as_ZX(): please contact the package maintainer.")
# as v_b cas:
if (rescue_nbr > rescue_thr["strictv"] && #rescue has been previously needed in the outer loop
damped_WLS_blob$breakcond != "low_pot") { #rescue has been previously needed in the outer loop
strictv_parent_info <- c(from=which_LevMar_step,breakcond=damped_WLS_blob$breakcond)
which_LevMar_step <- "strict_v|b"
} else {
v_parent_info <- c(from=which_LevMar_step, breakcond=damped_WLS_blob$breakcond)
which_LevMar_step <- "v"
v_iter <- 0L
}
}
prev1_rescued <- just_rescued
}
new_damping <- .get_new_damping(dampings_env$v[[which_LevMar_step]], which_LevMar_step)
damped_WLS_blob <- .wrap_do_damped_WLS_outer(
damped_WLS_fn = .do_damped_WLS_outer_spprec,
LevM_HL11_method=LevM_HL11_method, # contains the rescue_thr options => any possibility to simplify arguments ?
rescue= (is_HL1_1 && rescue_thr["rescue"]),
which_LevMar_step=which_LevMar_step,
old_relV_beta=old_relV_beta,
sXaug=sXaug, zInfo=zInfo,
old_Vscaled_beta=Vscaled_beta,
oldAPHLs=oldAPHLs,
APHLs_args = constant_APHLs_args,
damping=new_damping,
Trace= trace,
ypos=ypos,off=off,
GLMMbool=GLMMbool,etaFix=etaFix,
lambda_est=lambda_est,
wranefblob=wranefblob,seq_n_u_h=seq_n_u_h,
update_sXaug_constant_arglist=update_sXaug_constant_arglist,
#H_global_scale=H_global_scale,
ZAL_scaling= ZAL_scaling,
processed=processed,
phi_est=phi_est, n_u_h=n_u_h, ZAL=ZAL,
constant_v_infer_args=constant_v_infer_args,
looseness= if ( is.null(damped_WLS_blob) || ## start strict
new_damping>1e-7) {## use strict when there are trace of difficulties (in particular, failure to improve)
1 } else {processed$spaMM_tol$loose_fac},
low_pot=NULL ## explicit for clarity, but its the default
)
#old_m_grad_v <- m_grad_v
dampings_env$v[[attr(damped_WLS_blob,"step")]] <- damped_WLS_blob$damping
## LevM PQL
if (! is_HL1_1) {
if (damped_WLS_blob$lik < oldAPHLs$hlik) { ## if LevM step failed to find a damping that increases the hlik
# Tis should occur only bc of (1) numerically challenging conditions e.g mu close to bounds; or (2) optimum has been
# found and floating point innacurracies matter. We try to exclude the second case by the following test:
if ( ! ((breakcond <- damped_WLS_blob$breakcond)=="low_pot" && attr(breakcond,"very_low_pot"))) {
damped_WLS_blob <- NULL
.diagnose_conv_problem_LevM( beta_cov_info=.calc_beta_cov_info_spprec(X.pv = sXaug$AUGI0_ZX$X.pv,envir = sXaug$BLOB),
processed) # writes into 'processed'
dVscaled_beta <- get_from_MME(sXaug,szAug=zInfo) ################### FIT
Vscaled_beta <- list(v_h=Vscaled_beta$v_h+dVscaled_beta$dv_h,
beta_eta=Vscaled_beta$beta_eta+dVscaled_beta$dbeta_eta)
if (TRUE) { # not clear what is best here
break
} else {
LevenbergM <- FALSE ## desperate move...
# for (st in names_keep) assign(st,keep_init[[st]])
# Vscaled_beta <- c(v_h/ZAL_scaling ,beta_eta) ## RHS from assign(st,keep_init[[st]])
}
}
}
}
} else { ## IRLS: always accept new v_h_beta
dVscaled_beta <- get_from_MME(sXaug,szAug=zInfo) ################### FIT
devel <- FALSE
sanitize <- (FALSE || devel) # needed for devel test, speculative otherwise
# sanitize <- (TRUE || devel) # needed for devel test, speculative otherwise
if (sanitize) { # sanitize the PREVIOUS values
w_sane_etamo <- (sXaug$BLOB$WLS_mat_weights)*(muetablob$sane_eta-off)
rhs_sanitized_v_b <- c(v_h*attr(sXaug,"w.ranef") + .crossprod(ZAL, w_sane_etamo),
.crossprod(sXaug$AUGI0_ZX$X.pv, w_sane_etamo))
insane_Vscaled_beta <- Vscaled_beta
Vscaled_beta <- get_from_MME(sXaug, szAug =list(m_grad_obj=rhs_sanitized_v_b))
Vscaled_beta <- list(v_h=Vscaled_beta$dv_h,beta_eta=Vscaled_beta$dbeta_eta)
} # (but this may have its own problems; and avoid samitizing twice)
# UPDATE from sanitize or not previous value: new = old + dVscaled_beta
Vscaled_beta <- list(v_h=Vscaled_beta$v_h+dVscaled_beta$dv_h, ## DON'T REMOVE ME !!
beta_eta=Vscaled_beta$beta_eta+dVscaled_beta$dbeta_eta)
if (devel) { # devel code to check solutions by spprec vs spcorr.
# Difficult to extend to LevM case and may fail if Vscaled_beta is not sanaitized
testblob <- .spprec2spcorr(sXaug, zInfo=zInfo)
if ( ! is.null(v_h_beta_c <- testblob$v_h_beta)) {
# exclude case where spprec G is SPD but spcorr H is nonSPD => WLS_mat differ
if (diff(range(Vscaled_beta$v_h-v_h_beta_c$v_h))>1e-10) browser()
if (diff(range(Vscaled_beta$beta_eta-v_h_beta_c$beta))>1e-10) browser()
}#
}
##################
damped_WLS_blob <- NULL
if ( ! LMMbool && innerj>= checkpot_min_it) {
pot4improv <- get_from_MME(sXaug=sXaug, which="Mg_solve_g", B=zInfo$m_grad_obj)
}
}
if (trace>5L) .prompt()
##### Everything that is needed for
# (1) assessment of convergence: c(v_h*sqrt(wranefblob$w.ranef),beta_eta)
# (2) all return elements are updated as function of the latest Vscaled_beta
# (itself possible updated to new scaling by the following assign()'s).
# In particular We need muetablob and (if ! LMM) sXaug, hence a lot of stuff.
# Hence, the following code is useful whether a break occurs or not.
if ( is.null(damped_WLS_blob) ) { ## fits nothing, but updates variables in case of standard IRLS, or of intervals
WLS_blob <- .WLS_substitute_spprec(update_sXaug_constant_arglist, Vscaled_beta, off, etaFix, mod_attr=mod_attr,
lambda_est=lambda_est, ZAL,
processed, phi_est,
wranefblob, Trace=trace, stylefn)
for (st in names(WLS_blob)) assign(st,WLS_blob[[st]])
if (fpot_cond && ! LMMbool && is.null(for_intervals)) { # fpot_cond is FALSE except in possible private usage,
Mg_solve_g <- sum(.unlist(dVscaled_beta)*zInfo$m_grad_obj) # using old m_grad_obj
if (250*Mg_solve_g < fpot_tol) break
}
} else {
if (is_HL1_1 && (
(which_LevMar_step =="V_IN_B" && damped_WLS_blob$breakcond=="OK_gain") ||
(which_LevMar_step =="v" && damped_WLS_blob$breakcond=="low_pot")
) && damped_WLS_blob$APHLs$p_v>best_HL1_lik ) {
#cat(crayon::red("ICI!\n"))
best_HL1_damped_WLS_blob <- damped_WLS_blob
best_HL1_lik <- best_HL1_damped_WLS_blob$APHLs$p_v
#print(best_HL1_lik)
}
varnames <- intersect(names(damped_WLS_blob), # sXaug (and the weights) need not be present if(GLGLLM_const_w)
c("w.resid", ## !important! cf test-adjacency-corrMatrix.R
"sXaug","Vscaled_beta","wranefblob","v_h","u_h","muetablob"))
list2env(damped_WLS_blob[varnames], envir = environment())
if ( ! GLMMbool ) ZAL_scaling <- damped_WLS_blob$ZAL_scaling ## cf 'TAG:', this line in case I decide to use a ZAL_scaling !=1 later
}
# At this point all return elements are updated as function of the latest Vscaled_beta.
# In particular We need muetablob and (if ! LMM) sXaug, hence a lot of stuff.
#####
beta_eta <- Vscaled_beta$beta_eta
##### assessment of convergence
if (innerj==maxit.mean) {
if (maxit.mean>1L) {
if (LevenbergM) processed$LevenbergM["LM_start"] <- TRUE
if (trace) {
cat(crayon::red("!"))
} else if ( ! identical(processed$warned_maxit_mean, TRUE)) {
processed$warned_maxit_mean <- TRUE
if (!is.null(for_intervals)) {
message("Iterative algorithm converges slowly.")
} else message("Iterative algorithm converges slowly. See help('convergence') for suggestions.")
}
}
break
} else {
relV_beta <- c(v_h*sqrt(wranefblob$w.ranef),beta_eta) ## convergence on v_h relative to sqrt(lambda), more exactly for Gaussian
abs_d_relV_beta <- abs(relV_beta - old_relV_beta) ## for ML, comparison between estimates when ( hlik_stuck || ! need_v_step )
if (is_HL1_1 && LevenbergM) {
not_moving <- (
# exclude cases of possible p_v overfit by v_h, such as cases "b" "v_b" "b_from_v_b" "b_&_v_in_b"
(which_LevMar_step %in% c("v", "V_IN_B","strict_v|b") || default_b_step=="v_in_b") &&
( ! is.null(old_relV_beta)) &&
{
#print(c(mean(abs_d_relV_beta[seq_n_u_h]),relV_beta[-seq_n_u_h],old_relV_beta[-seq_n_u_h]))
meanmean <- mean(c(mean(abs_d_relV_beta[seq_n_u_h]), mean(abs_d_relV_beta[-seq_n_u_h])), na.rm=TRUE) # second mean may be NaN
meanmean < d_relV_b_tol_LM
}
)
} else not_moving <- (
( ! is.null(old_relV_beta)) &&
{
meanmean <- mean(c(mean(abs_d_relV_beta[seq_n_u_h]), mean(abs_d_relV_beta[-seq_n_u_h])), na.rm=TRUE) # second mean may be NaN
meanmean < d_relV_b_tol
}
) #In ! LevM, v&b are fitted simultaneously without damping
if (not_moving) {
# not_moving_Wattr <- .diagnose_coeff_not_moving(coeff_not_moving = not_moving,relV_beta, damped_WLS_blob, innerj,
# damping, is_HL1_1, oldAPHLs, Ftol=processed$spaMM_tol$Ftol_LM, trace, LevenbergM,stylefn=identity)
break
}
# More ad hoc breaks for cases where the coefficients keep moving although the total potential is low:
if ( ! is.null(damped_WLS_blob) ) {
if (is_HL1_1) {
if ( which_LevMar_step=="v_b" && damped_WLS_blob$breakcond=="low_pot" && attr(damped_WLS_blob$breakcond,"no_overfit")) {
break
} else if ( which_LevMar_step=="v" &&
damped_WLS_blob$breakcond=="low_pot" && ## essential condition
( ( ! pforpv ) || # outer beta
v_parent_info[["breakcond"]]=="low_pot" ## essential condition; otherwise poor fits (eg test-nloptr # 422 p_v_out_def )
)
) {
break
} else if ( which_LevMar_step=="v_b" &&
damped_WLS_blob$breakcond=="stuck_obj" ) {
# I removed this break from [v3.12.34 -> 54] with poor effect on test-nloptr #362... (spcorr)
break # this case occurs in the test-nloptr tests
# not an obvious termination condition, but seems OK following "v"&&"low_pot" (and stops the alternation between these two states)
# but if previous steps were "V_IN_B" & OK_gain, it may be worth returning to such a step
} else if ( which_LevMar_step=="strict_v|b" &&
# => so that the $breakcond is { "v|b_no_loop" (i.e. no damping loop in this case), with "v_pot4improv" attribute tested next:}
strictv_parent_info[["breakcond"]]=="stuck_obj" && # we don't test "low_pot" here bc this case does not occur according to previous code
attr(damped_WLS_blob$breakcond,"v_pot4improv")<1e-10) {
} else if ( which_LevMar_step =="V_IN_B" && damped_WLS_blob$breakcond=="low_pot") {
# I could further test attr(damped_WLS_blob$breakcond,"very_low_pot") here
break # motivated by poisson 'smaller' fit in test_COMPoisson_difficult with control.HLfit=list(max.iter.mean=1000),
} #else {print(which_LevMar_step); str(damped_WLS_blob$breakcond)}
} else {# tendency of PQL/L meanmean to converge very slowly despite low_pot (and indeed tiny changes in hlik) => ad hoc adjustment
if (damped_WLS_blob$breakcond=="low_pot" && attr(damped_WLS_blob$breakcond,"very_low_pot")) break
}
}
#
if ( ! (is_HL1_1 && LevenbergM)) { ## This is for the special case of reversal of LevenbergM condition from F to T in LevM PQL !!!!
old_relV_beta <- relV_beta
} ## ELSE old_relV_beta was updated when required in block for which_LevMar_step.
}
} ################ E N D LOOP ##############
#if (trace>4L) browser()
if ( ! is.null(for_intervals) && for_intervals$phi_pred_OK) {
warnobjfn <- names(for_intervals$warnlik)
warnlik <- unlist(do.call(".calc_APHLs_from_ZX",loc_logLik_args)[warnobjfn])
if ((for_intervals$warnlik-warnlik) < -1e-4 &&
(is.null(bestlik <- processed$envir$confint_best$lik) || warnlik > bestlik)) {
if (is.null(bestlik)) {
locmess <- paste("A higher",warnobjfn,"was found than for the original fit.",
"\nThis suggests the original fit did not fully maximize",warnobjfn,
"\n (REML fits, or numerical accuracy issues). Expect more information at end of computation.")
message(locmess)
}
processed$envir$confint_best$lik <- warnlik
processed$envir$confint_best$beta_eta <- .unscale(sXaug$AUGI0_ZX$X.pv, Vscaled_beta[n_u_h+seq_len(pforpv)])
processed$envir$confint_best$ranPars <- for_intervals$ranFix
processed$envir$confint_best$ranPars$lambda <- unique(lambda_est) # ugly but get_inits_from_fit() cannot be used in this context...
}
}
if (trace>1L && (LevenbergM)) {
stylefn <- .spaMM.data$options$stylefns$betalast
if (pforpv) { # outer beta
maxs_grad <- c(max(abs(zInfo$m_grad_obj[seq_n_u_h])),max(abs(zInfo$m_grad_obj[-seq_n_u_h])))
} else maxs_grad <- c(max(abs(zInfo$m_grad_obj[seq_n_u_h])),0)
cat(stylefn("iter=",innerj,", max(|grad|): v=",maxs_grad[1L],"beta=",maxs_grad[2L],";\n"))
}
if (best_HL1_lik > -Inf) {
if (best_HL1_lik> damped_WLS_blob$APHLs$p_v) {
# cat("restoring better fit")
damped_WLS_blob <- best_HL1_damped_WLS_blob
varnames <- intersect(names(damped_WLS_blob), # sXaug (and the weights) need not be present if(GLGLLM_const_w)
c("w.resid", ## !important! cf test-adjacency-corrMatrix.R
"sXaug","Vscaled_beta","wranefblob","v_h","u_h","muetablob"))
list2env(damped_WLS_blob[varnames], envir = environment())
}
}
names(beta_eta) <- colnames(processed$AUGI0_ZX$X.pv)
RESU <- list(sXaug=sXaug,
## used by calc_APHLs_from_ZX
#fitted=fitted, ## FIXME: removed so that no shortcut a la Bates in calc_APHLs_from_ZX; reimplement the shorcut in that fn?
#weight_X=NA,
nobs=nobs, pforpv=pforpv, seq_n_u_h=seq_n_u_h, u_h=u_h,
muetablob=muetablob,
lambda_est=lambda_est,
phi_est=phi_est,
## used by other code
beta_eta=beta_eta, w.resid=w.resid, wranefblob=wranefblob,
v_h=v_h, eta=muetablob$sane_eta, innerj=innerj)
## for diagnostic purposes
if ( ! LMMbool ) RESU$m_grad_obj <- zInfo$gainratio_grad ## frm ZInfo bc other copies of m_grad_obj may be missing if not LevenbergM
return(RESU)
}
.intervalStep_spprec <- function(old_v_h_beta,sXaug,zInfo,currentlik,for_intervals,currentDy,ZAL, etamo) {
## voir code avant 18/10/2014 pour une implem rustique de VenzonM pour debugage
## somewhat more robust algo (FR->FR: still improvable ?), updates according to a quadratic form of lik near max
## then target.dX = (current.dX)*sqrt(target.dY/current.dY) where dX,dY are relative to the ML x and y
## A nice thing of this conception is that if the target lik cannot be approached,
## the inferred x converges to the ML x => this x won't be recognized as a CI bound (important for locoptim)
parmcol_ZX <- for_intervals$parmcol_ZX
v_h_beta_vec <- old_v_h_beta_vec <- unlist(old_v_h_beta)
v_h_beta_vec[] <- NA
if (currentDy <0) {
v_h_beta_vec[parmcol_ZX] <- old_v_h_beta_vec[parmcol_ZX]
} else {
currentDx <- (old_v_h_beta_vec[parmcol_ZX]-for_intervals$MLparm)
targetDy <- (for_intervals$fixeflik-for_intervals$targetlik)
Dx <- currentDx*sqrt(targetDy/currentDy)
#cat(currentDx," ",targetDy," ",currentDy," ",Dx,"\n")
## pb is if Dx=0 , Dx'=0... and Dx=0 can occur while p_v is still far from the target, because other params have not converged.
## (fixme) patch:
if (currentDy<targetDy) { ## we are close to the ML: we extrapolate a bit more confidently
min_abs_Dx <- for_intervals$asympto_abs_Dparm/1000
} else min_abs_Dx <- 1e-6 ## far from ML: more cautious move our of Dx=0
Dx <- sign(currentDx)*max(abs(Dx),min_abs_Dx)
v_h_beta_vec[parmcol_ZX] <- for_intervals$MLparm + Dx
}
## gradient changes: (the get_from_MME() call only uses m_grad_obj to get rhs)
# scratch code that led o this is in devel_confint_spprec.R
parmcol_X <- for_intervals$parmcol_X
off_newparm <- sXaug$AUGI0_ZX$X.pv[,parmcol_X]*v_h_beta_vec[parmcol_ZX] # function of NEW tentative value of MLparm
# so it does not give the difference between etamo and i_etamo, which depends on OLD value of MLparm
oovb <- old_v_h_beta_vec[-(parmcol_ZX)]
seq_n_u_h <- seq_along(old_v_h_beta$v_h)
i_etamo <- drop(ZAL %*% oovb[seq_n_u_h] +
sXaug$AUGI0_ZX$X.pv[,-(parmcol_X),drop=FALSE] %*% oovb[-seq_n_u_h]) #
#
rhs <- sXaug$BLOB$WLS_mat_weights*(etamo-off_newparm- i_etamo) #_F I X M E__ this block code looks inefficient and unclear...
zInfo$m_grad_obj <- zInfo$m_grad_obj[-parmcol_ZX]+
c(drop(crossprod(ZAL,rhs)),
drop(crossprod(sXaug$AUGI0_ZX$X.pv[,-parmcol_X,drop=FALSE],rhs))) # tjrs zInfoo$m_grad_obj
#
int_sXaug <- sXaug # list with two elements: environments AUGI0_ZX and BLOB
# Modify $AUGI0_ZX:
int_AUGI0_ZX <- as.list(sXaug$AUGI0_ZX) ## as.list() local copy avoids global modifs of original sXaug$AUGI0_ZX which is an envir
int_AUGI0_ZX$X.pv <- int_AUGI0_ZX$X.pv[,-(parmcol_X),drop=FALSE]
int_AUGI0_ZX$ZeroBlock <- int_AUGI0_ZX$ZeroBlock[,-(parmcol_X),drop=FALSE]
int_sXaug$AUGI0_ZX <- int_AUGI0_ZX ## Replaces an envir by a list in the local copy;
# Modify $BLOB, leaving untouched the non-promises ("G_CHMfactor" "Gmat" "chol_Q" "pMat_G") .
# ***which will presumably be written over since they are not treated as promises***
# But any preexisting evaluated promise should in principle be written over.
# A conflict in dimension between functions of X would likely be a conflict between previously and newly evaluated promise
.init_promises_spprec(sXaug=int_sXaug) # reinit X.pv promises (ZtWX XtWX r12 qrXa DpD LZtWX...)
# and non-X.Pv promises (Md2hdv2 tcrossfac_Md2hdv2 inv_L_G_ZtsqrW invL_G.P sortPerm...), which may not be needed.
# This won't use the original fit promises anyway (as .confint_LRT_single_par() -> get_HLCorcall() -> .preprocess() -> fresh .init_promises_spprec)
# This attr(.,"pforpv") is used to determined whether there are beta coeffs to compute in .AUGI0_ZX_spprec()
attr(int_sXaug,"pforpv") <- attr(int_sXaug,"pforpv") - length(parmcol_X) # a priori -1
v_h_beta_vec[-(parmcol_ZX)] <- old_v_h_beta_vec[-(parmcol_ZX)]+ unlist(get_from_MME(int_sXaug,szAug=zInfo))
return(list(v_h_beta=relist(v_h_beta_vec,old_v_h_beta)))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.