From 9c4598d97f2a899a9b64bc8221bb16613a27cd06 Mon Sep 17 00:00:00 2001 From: Jim Thorson Date: Wed, 1 Jan 2025 12:14:13 -0800 Subject: [PATCH] finish initial dsemRTMB --- R/compute_nll.R | 87 +++++++++++++++++++++++++--------- R/dsemRTMB.R | 55 ---------------------- scratch/example.RDS | Bin 2774 -> 1118 bytes scratch/example2.R | 81 +++++++++++++++++++++++++++++++ scratch/example2_tmp.R | 102 ++++++++++++++++++++++++++++++++++++++++ scratch/test_dsemRTMB.R | 20 +++++++- 6 files changed, 267 insertions(+), 78 deletions(-) create mode 100644 scratch/example2.R create mode 100644 scratch/example2_tmp.R diff --git a/R/compute_nll.R b/R/compute_nll.R index c91bcf6..1ab4131 100644 --- a/R/compute_nll.R +++ b/R/compute_nll.R @@ -12,6 +12,24 @@ function( parlist, "c" <- ADoverload("c") "[<-" <- ADoverload("[<-") + # Temporary fix for solve(adsparse) returning dense-matrix + sparse_solve = function(x){ + invx = solve(x) + if( RTMB:::ad_context() ){ + out = sparseMatrix( + i = row(invx), + j = col(invx), + x = 1, + ) + out = AD(out) + out@x = invx + #out = drop0(out) # drop0 doesn't work + return(out) + }else{ + return(invx) + } + } + # Unpack parameters explicitly beta_z = parlist$beta_z delta0_j = parlist$delta0_j @@ -46,6 +64,47 @@ function( parlist, Gamma_kk = AD(ram$G_kk) V_kk = t(Gamma_kk) %*% Gamma_kk + # Rescale I-Rho and Gamma if using constant marginal variance options + if( (options[2]==1) || (options[2]==2) ){ + invIminusRho_kk = sparse_solve(IminusRho_kk) # solve(adsparse) returns dense-matrix + #print(class(invIminusRho_kk)) + + # Hadamard squared LU-decomposition + # See: https://eigen.tuxfamily.org/dox/group__QuickRefPage.html + squared_invIminusRho_kk = invIminusRho_kk + #print(class(squared_invIminusRho_kk)) + #print(invIminusRho_kk) + #REPORT( invIminusRho_kk ) + squared_invIminusRho_kk@x = squared_invIminusRho_kk@x^2 + + if( options[2] == 1 ){ + # 1-matrix + ones_k1 = matrix(1, nrow=n_k, ncol=1) + + # Calculate diag( t(Gamma) * Gamma ) + squared_Gamma_kk = Gamma_kk + squared_Gamma_kk@x = squared_Gamma_kk@x^2 + sigma2_k1 = t(squared_Gamma_kk) %*% ones_k1; + + # Rowsums + margvar_k = solve(squared_invIminusRho_kk, sigma2_k1) + + # Rescale IminusRho_kk and Gamma + invmargsd_kk = invsigma_kk = AD(Diagonal(n_k)) + invmargsd_kk@x = 1 / sqrt(margvar_k) + invsigma_kk@x = 1 / sqrt(sigma2_k1) + IminusRho_kk = invmargsd_kk %*% IminusRho_kk; + Gamma_kk = invsigma_kk %*% Gamma_kk; + }else{ + # calculate diag(Gamma)^2 + targetvar_k = diag(Gamma_kk)^2 + + # Rescale Gamma + margvar_k = solve(squared_invIminusRho_kk, targetvar_k) + diag(Gamma_kk) = sqrt(margvar_k) + } + } + # Calculate effect of initial condition -- SPARSE version delta_tj = matrix( 0, nrow=n_t, ncol=n_j ) if( length(delta0_j)>0 ){ @@ -62,34 +121,20 @@ function( parlist, # Full rank GMRF z_tj = x_tj - # Doesn't work ... mat2triplet not implemented - #V_kk = matrix(0, nrow=n_k, ncol=n_k) - #REPORT( V_kk ) - #V_kk = as.matrix(V_kk) - #invV_kk = solve(V_kk) - #Qtmp_kk = invV_kk - #Qtmp_kk = AD(Qtmp_kk) - #REPORT( Qtmp_kk ) - #Q_kk = Qtmp_kk - - # Works - invV_kk = AD(ram$G_kk) - invV_kk@x = 1 / Gamma_kk@x^2 - Q_kk = t(IminusRho_kk) %*% invV_kk %*% IminusRho_kk + # Works for diagonal + #invV_kk = AD(ram$G_kk) + #invV_kk@x = 1 / Gamma_kk@x^2 + #Q_kk = t(IminusRho_kk) %*% invV_kk %*% IminusRho_kk - # Experiment - #IminusRho_dense = as.matrix( IminusRho_kk ) - #V_dense = as.matrix( V_kk ) - #Q_RHS = solve(V_kk, IminusRho_kk) - #Q_kk = t(IminusRho_kk) %*% Q_RHS + # Works in general + invV_kk = sparse_solve( V_kk ) + Q_kk = t(IminusRho_kk) %*% invV_kk %*% IminusRho_kk # Fine from here jnll_gmrf = -1 * dgmrf( as.vector(z_tj), mu=as.vector(xhat_tj + delta_tj), Q=Q_kk, log=TRUE ) - #jnll_gmrf = 0 REPORT( Q_kk ) }else{ # Reduced rank projection .. dgmrf is lower precision than GMRF in CPP - #jnll_gmrf = -1 * dgmrf( as.vector(x_tj), mu=rep(1,n_k), Q=Diagonal(n_k), log=TRUE ) jnll_gmrf = -1 * sum( dnorm(x_tj, mean=0, sd=1, log=TRUE) ) # diff --git a/R/dsemRTMB.R b/R/dsemRTMB.R index f99b3a2..c2931d8 100644 --- a/R/dsemRTMB.R +++ b/R/dsemRTMB.R @@ -69,11 +69,6 @@ function( sem, delta0_j = rep(0, n_j), x_tj = ifelse( is.na(y_tj), 0, y_tj ) ) - #if( control$gmrf_parameterization=="separable" ){ - # Params$x_tj = ifelse( is.na(tsdata), 0, tsdata ) - #}else{ - # Params$eps_tj = ifelse( is.na(tsdata), 0, tsdata ) - #} # Turn off initial conditions if( estimate_delta0==FALSE ){ @@ -113,8 +108,6 @@ function( sem, } # Build object - #source( file.path(R'(C:\Users\James.Thorson\Desktop\Git\dsem\R)', "make_matrices.R") ) - #source( file.path(R'(C:\Users\James.Thorson\Desktop\Git\dsem\R)', "get_jnll.R") ) cmb <- function(f, ...) function(p) f(p, ...) ## Helper to make closure #f(parlist, model, tsdata, family) obj = RTMB::MakeADFun( @@ -130,36 +123,6 @@ function( sem, profile = control$profile, silent = TRUE ) - # - if( FALSE ){ - repRTMB = obj$report() - range(Rep$Gamma_kk - repRTMB$Gamma_kk) - range(Rep$Rho_kk - repRTMB$Rho_kk) - range(Rep$Q_kk - repRTMB$Q_kk) - range(Rep$mu_tj - repRTMB$mu_tj) - range(Rep$delta_tj - repRTMB$delta_tj) - range(Rep$xhat_tj - repRTMB$xhat_tj) - range(Rep$jnll_gmrf - repRTMB$jnll_gmrf) - - # - dgmrf( repRTMB$z_tj, mu=repRTMB$xhat_tj + repRTMB$delta_tj, Q=repRTMB$Q_kk, log=TRUE ) - dgmrf( Rep$z_tj, mu=Rep$xhat_tj + Rep$delta_tj, Q=Rep$Q_kk, log=TRUE ) - dgmrf( as.vector(Rep$z_tj), mu=as.vector(Rep$xhat_tj + Rep$delta_tj), Q=Rep$Q_kk, log=TRUE ) - mvtnorm::dmvnorm( as.vector(Rep$z_tj - (Rep$xhat_tj + Rep$delta_tj)), sigma=as.matrix(solve(Rep$Q_kk)), log=TRUE ) - dGMRF( as.vector(Rep$z_tj), mu=as.vector(Rep$xhat_tj + Rep$delta_tj), Q=Rep$Q_kk ) - - # - Rep = fit$obj$report( fit$obj$env$last.par ) - - # - example = list( - y = repRTMB$z_tj, - mu = repRTMB$xhat_tj + repRTMB$delta_tj, - Q = Rep$Q_kk - ) - saveRDS( example, file=file.path(R'(C:\Users\James.Thorson\Desktop\Git\dsem\scratch)','example.RDS') ) - } - if(control$quiet==FALSE) list_parameters(obj) # bundle internal = list( @@ -171,24 +134,6 @@ function( sem, covs = covs ) - # Parse priors -# if( !is.null(prior_negloglike) ){ -# # prior_negloglike = \(obj) -dnorm(obj$par[1],0,1,log=TRUE) -# prior_value = tryCatch( expr = prior_negloglike(obj) ) -# if( is.na(prior_value) ) stop("Check `prior_negloglike(obj$par)`") -# obj$fn_orig = obj$fn -# obj$gr_orig = obj$gr -# -# # BUild prior evaluator -# requireNamespace(RTMB) -# priors_obj = RTMB::MakeADFun( func = prior_negloglike, -# parameters = list(par=obj$par), -# silent = TRUE ) -# obj$fn = \(pars) obj$fn_orig(pars) + priors_obj$fn(pars) -# obj$gr = \(pars) obj$gr_orig(pars) + priors_obj$gr(pars) -# internal$priors_obj = priors_obj -# } - # Further bundle out = list( "obj" = obj, "ram" = ram, diff --git a/scratch/example.RDS b/scratch/example.RDS index 56452ce7c18e8700dff74828f9e2af3db8ec4524..2dbd1c99a9756dbde44ce082498be426422393bb 100644 GIT binary patch literal 1118 zcmV-k1flyMiwFP!000002JP5+Y!p=-!11?7FFdFqSOt-*9CBDxEEI7Usl}r}!6In6 z#x0bx-4@%5pj?VYE@=^O1cjmqT1AR-(-uYXK&b_l#Nc1X7&YR5{v-JPb?0Z6)x?;% z(Mh(+r*CKH&Ac~n-kbUTc2XN-jK_GrF5^kyqe+;MpEF>HG4A$kxQ)wmxiMZhpIydl zQlu>A-$mGA!hf#TigFfe${CP|BqT!)l2$Ijh4!%*qb*W#3EDv#yc8YK5tktibiiDJ zPUws)(FIqbE4raOdY~uL(F?uN2Yqoh`k_AtU?8qR1_t3;T!-s17&jmj41^hqEDXbq z7>=9ZgC8SsGe#mCw;%_ja4T-Z?YINE7>zL)i#&|Oc-)D6Ou${3h)K8`_n-iiF$GgG z4bx%aUIb8x8JLM8%!1sn*$AQpA(WyFb8sKZQGrU_kGXgN^DrL|;vp=+LM*~!EWyKg z1WWNKs;~@?VL4V{C01cI9>*F~qXuj71lHk6JcXz644%bvcpmHV0ybbHUc^h-gw3eM z7QBpCuobUj8@A&$ypB4&fj6TkwhV?$&bLhXFFvv36WgH^TOQ8mO>R$iu?&w)QFnOv zpQ~k~?xs#@PAUC83briFZ&~b= zGU~ycG}8)HmX+)C+eGb72$WbA>~yC_e3?>MG;*|6SzbIZ`ebvgK+r0ZcWz|;M6GPU z8ve)+T<7f9%PAAPXhhGe=hpLU8`_SxrR`~(+OD>(?Q0*jAKDk~kM>FXrG3->X&<$p zB32QOo>wChZ9}6IZBHW=ZCj%j?T1D%+9!=>w0|1eXkRtT(SB>hqt~F(k6w>PLV9f) z73p=#=d}^ntkIHQzeY|v78*tAxM)PBW2DiQj+aJSI(8a$={RZx7P8A&-LDu$j{^vfgBk$-0uYBfn}2+zuW9f0y5Ykn`cAiGm*jra&UZ3nQTp=a z`uP6%{{O{({r6e-PfB%`?o!$Mux7H)ex)DM=h^MC;@$jNu;J72|lXUu6u&rTRgI~jV4fAQ3U&H=1>|ew2XgEGC zt}23kw79$M^PZN_LvNEC+2^>HykFl=-KG1}F)n?O9#0R#wtmD~D#j(gj(LL%*hh=Q zJn{GScv6HOPg{<8G>pgbWMu2xnP2PUmSBI;IL;^>&zxj^I|j$CeW#C)!T#s*7W#Y~ zU+j7PxVRpJ{m;SqMPvTC@M~!273kw)5ucCaip6ooK~Kc|}EOb(+HGN1V?;gP#wMFBa><&nF(^5-?99 z+9X`pWXzL{^G?P*Nf@6X-DeB2-toBp349#xj}++1{Jz3CexGEl*SFnpTwZ8>us*!? z#W+92yJ4MrK<^2iKR^D^12DcH+Wr{FpR0k`9)|IQF>W~85g5mxyOC&doVIuL*75}N z1Y#V2o&%xt{raopJRh`m*Z0>pQ1`VV&hv@wUl>n#3myBn@Hk(e0GwAroqiq#xGn{M z*0&4w1=9+Dq5F&aqrI*Bthhc_+(*TW&_33ECGZztg8vIWZpHO2`V=jn|KIRM>21+} zFhBB(^$^a;)1Cgujd&72+Ir+!OKH=S!sPxGU&-Q2EoC{S`Bk(Z0N)def4f)rcvf3W z^ZAs8WfkapRcy>Dy9XsFHYGHvMhi&jre&FC_1;=|Tk zEzKWvnKs~CzqFv%zjT$>@}qU%fjhr^yQ{QX(YI6TTRt4#Ra(7Jxc%2yKxmG3^B5!* z5;v=E^qhrC#LSc04SIh7WilpWefCvi3Tc0NXz$Yr$y8Qt&>u_*Wn}~XA(SvyInZZc z{lZ!8Kz|q|0td+Y^aMyG#Si@vlxdW1=m$}zvr2_N`-(Dy(gS_=sxT8L&-$Y&QIuZj zvsVaquD#J8OPNjd%E|idRVbR$2mNsr_Vv0i`XQ9LM6bH6Kb~Tt1foBI= z`wIGS=#kL<1bqbbY0$d~dJy#K(7OxzNa!=5_Ym|^&}TyLDd?l2M?vo;=wqPIg5F!u z$3mYC-Cxjyp+`gSBk1Fx&w<`o&_kflg&rX2jQNMf(kS_qt;Cp*SS+1VK-or&`H96cD20^m#F(#GER#}1 z*+Go?i^UdESiav$jQNbkvM9xrUBsB*SZpz+gtD6$^Bs#Vp|F>MJ;ZpQve;4zd$Xw{ z_WZ~1OLpI~dz0O}?Ec(Kd5&_Na)Q`%BKOC_zXJXW_@9UWB(dj8?q2}^O87Ux{{s9? z#GXUBKMwv?@K?hBBK*z7o?E$pA^fZ1uY&(2_)igg&gK4i_}9Q+4gWs)PZN7C=Kciu z*TP=|e?9zXh&@MheAO03%&)wXg1b+_v8{vN${{~-LW#9jlqKNbEw__x4+2>$cLUK_YS4gP%ix59rI{tLujGq^t; z{sQ>7!T$>U7m2-=aDN8;h462O{|Nk-h`q*eeRc!>);3i=Q`In4-BQ&mRee&`Bvlw;VfE~^a7|dPIwM>UR)<=IYr$&FS>Zac`f*OU2CO!;3hRG}=sBOS zJ z*0%E>K-hyI{sDx4+$7v3JS03Nyd=COd?b7&{3NMLbwwmoDk`RX-=5#gc(kl>4YdJ%yPnPCqz48juYlO z!QuqX2{BHX_Y;}#0L$kRlzm|+ZF6kuv5V<1-ljO zQBbE~ufTH>&r7@@@uI{_68j|TCH70aEO9{Mpu{1G!xFDZ9Fb^{cva%4#A_1ABpM}- z3!IQRDbXa+EOAQWw8R;S7KyVG=OkJs&P!a7xF~T+;hrGqNHP1WEF@?R?Snf0>B;GxmZ~NrOjSu$b z9kM*|Ze3XY;X56E%rB3#9U0wo_~%bzW!r-fsIi;G^F*OnKLKYgEJgS6&}s zyLzO4f3fADVSdbO?#JAZxgT>s=6=lmnENsJ(=q=C__3`><`R;aY zq%HWNvO7T(A{4U{qhR4QgW&*t()Y@Dspm5*#BZm|92-`_U2~gPffh|?@`zg c5B~d;&FHw7)^LgaT=dER0BR(^ga}ap0OuyZvH$=8 diff --git a/scratch/example2.R b/scratch/example2.R new file mode 100644 index 0000000..77e34e4 --- /dev/null +++ b/scratch/example2.R @@ -0,0 +1,81 @@ + library(Matrix) + library(RTMB) + + # Load minimal example + example = readRDS( file=file.path(R'(C:\Users\James.Thorson\Desktop\Git\dsem\scratch)',"example.RDS")) + attach(example) + set.seed(101) + z = rnorm(nrow(Rho_kk)) + + # Define function with both versions + f = function(p){ + # + sparse_solve = function(x){ + invx = solve(x) + if( RTMB:::ad_context() ){ + out = sparseMatrix( + i = row(invx), + j = col(invx), + x = 1, + ) + out = AD(out) + out@x = invx + #out = drop0(out) # drop0 doesn't work + return(out) + }else{ + return(invx) + } + } + + Gamma_kk = AD( p[1] * Gamma_kk ) + Rho_kk = AD( p[2] * Rho_kk ) + V_kk = t(Gamma_kk) %*% Gamma_kk + IminusRho_kk = Diagonal(nrow(Rho_kk)) - Rho_kk + + # Option-1 ... works + invV_kk = Gamma_kk + invV_kk@x = 1 / Gamma_kk@x^2 + Q1_kk = t(IminusRho_kk) %*% invV_kk %*% IminusRho_kk + REPORT( Q1_kk ) + + # Option-2 + invV_kk = solve(V_kk) + Q2_kk = t(IminusRho_kk) %*% invV_kk %*% IminusRho_kk + REPORT( Q2_kk ) + + # Option-3 + invV_kk = solve(V_kk) + #if( RTMB:::ad_context() ){ + # temp_kk = sparseMatrix( + # i = row(invV_kk), + # j = col(invV_kk), + # x = 1, + # ) + # temp_kk = AD(temp_kk) + # temp_kk@x = invV_kk + # invV_kk = temp_kk + #} + invV_kk = sparse_solve(V_kk) + Q3_kk = t(IminusRho_kk) %*% invV_kk %*% IminusRho_kk + REPORT( Q3_kk ) + + # Option-1 WORKS for diagonal matrix only + #jnll_gmrf = -1 * dgmrf( z, mu=rep(0,nrow(Rho_kk)), Q=Q1_kk, log=TRUE ) + + # Option-2 DOES NOT WORK + #jnll_gmrf = -1 * dgmrf( z, mu=rep(0,nrow(Rho_kk)), Q=Q2_kk, log=TRUE ) + + # Option-3 WORKS generally, but requires a manual coersion from dense to sparse matrix + jnll_gmrf = -1 * dgmrf( z, mu=rep(0,nrow(Rho_kk)), Q=Q3_kk, log=TRUE ) + return(jnll_gmrf) + } + + # Works for Option-1 + f(c(1,1)) + obj = MakeADFun( f, c(1,1) ) + obj$fn(obj$par) + + # Show they're identical if compiled using Option-1 + Rep = obj$report() + identical( Rep$Q1_kk, Rep$Q2_kk ) + identical( Rep$Q1_kk, Rep$Q3_kk ) diff --git a/scratch/example2_tmp.R b/scratch/example2_tmp.R new file mode 100644 index 0000000..1e9fd66 --- /dev/null +++ b/scratch/example2_tmp.R @@ -0,0 +1,102 @@ + library(Matrix) + library(RTMB) + + # Load minimal example + example = readRDS( file=file.path(R'(C:\Users\James.Thorson\Desktop\Git\dsem\scratch)',"example.RDS")) + attach(example) + set.seed(101) + z = rnorm(nrow(Rho_kk)) + + # Define function with both versions + f = function(p){ + Gamma_kk = AD( p[1] * Gamma_kk ) + Rho_kk = AD( p[2] * Rho_kk ) + V_kk = t(Gamma_kk) %*% Gamma_kk + IminusRho_kk = Diagonal(nrow(Rho_kk)) - Rho_kk + + # Option-1 ... works + invV_kk = Gamma_kk + invV_kk@x = 1 / Gamma_kk@x^2 + Q1_kk = t(IminusRho_kk) %*% invV_kk %*% IminusRho_kk + REPORT( Q1_kk ) + + # Option-2 + invV_kk = solve(V_kk) + #trip = mat2triplet(invV_kk) + #print(trip) + #print(invV_kk@i) + #print(invV_kk@p) + #print(invV_kk@x) + #p = invV_kk@p + #n.p = length(p) + #dp <- p[-1L] - p[-n.p] + #n.dp = length(dp) + #j = rep.int(seq.int(from = 0L, length.out = n.dp), dp) + 1 + #print(j) + #invV_kk = sparseMatrix( + # i = invV_kk@i + 1, + # p = invV_kk@p, + # #j = j, + # x = invV_kk@x + # ) + if( RTMB:::ad_context() ){ + #print(invV_kk[1:10,1:10]) + #class(invV_kk) + #print(row(invV_kk)) + + #vec = as.matrix(invV_kk)[1,1] + #print(vec) + + #Dim = rep(nrow(V_kk),2) + #invV2_kk = drop0(sparseMatrix( i=1, j=1, x=0, dims=Dim )) # Make with a zero + #for(i in seq_len(nrow(invV_kk)) ){ + #for(j in seq_len(ncol(invV_kk)) ){ + # Ind_kk = sparseMatrix( + # i = i, + # j = j, + # x = 1, + # dims = Dim + # ) + # invV2_kk = invV2_kk + invV_kk[i,j] * Ind_kk + #}} + invV2_kk = sparseMatrix( + i = row(invV_kk), + j = col(invV_kk), + x = 1, + ) + invV2_kk = AD(invV2_kk) + invV2_kk@x = invV_kk + REPORT(invV2_kk) + }else{ + invV2_kk = invV_kk + } + #invV2_kk@x = invV_kk@x + #invV_kk = AD(invV_kk) + #print(length(invV2_kk@x)) + #print(length(invV_kk@x)) + #print(invV2_kk) + Q2_kk = t(IminusRho_kk) %*% invV2_kk %*% IminusRho_kk + REPORT( Q2_kk ) + #print(Q2_kk) + + # Option-3 ... try rebuilding + #Q3_kk = AD(sparseMatrix(i=Q2_kk@i, p=Q2_kk@p, x=Q2_kk@x)) + #REPORT( Q3_kk ) + + # Option-1 WORKS + #jnll_gmrf = -1 * dgmrf( z, mu=rep(0,nrow(Rho_kk)), Q=Q1_kk, log=TRUE ) + + # Option-2 DOES NOT WORK + jnll_gmrf = -1 * dgmrf( z, mu=rep(0,nrow(Rho_kk)), Q=Q2_kk, log=TRUE ) + #return(jnll_gmrf) + } + + # Works for Option-1 + f(c(1,1)) + obj = MakeADFun( f, c(1,1) ) +# obj$fn(obj$par) + + # Show they're identical if compiled using Option-1 +# Rep = obj$report() +# identical( Rep$Q1_kk, Rep$Q2_kk ) +# attr(Rep$Q1_kk,"Dim") diff --git a/scratch/test_dsemRTMB.R b/scratch/test_dsemRTMB.R index b2d1954..403286f 100644 --- a/scratch/test_dsemRTMB.R +++ b/scratch/test_dsemRTMB.R @@ -57,6 +57,7 @@ fit = dsem( sem=sem, run_model = TRUE, use_REML = TRUE, gmrf_parameterization = "projection", + constant_variance = c("conditional", "marginal", "diagonal")[2], # marginal not working map = Map, parameters = Params ) ) @@ -78,19 +79,34 @@ source( file.path(R'(C:\Users\James.Thorson\Desktop\Git\dsem\R)', "read_model.R" source( file.path(R'(C:\Users\James.Thorson\Desktop\Git\dsem\R)', "dsemRTMB.R") ) # Define prior -log_prior = function(p) dnorm( p$beta_z[1], mean=0, sd=0.01, log=TRUE) +log_prior = function(p) dnorm( p$beta_z[1], mean=0, sd=0.1, log=TRUE) fitRTMB = dsemRTMB( sem = sem, tsdata = tsdata, estimate_delta0 = TRUE, family = rep("normal",ncol(tsdata)), log_prior = log_prior, - control = dsem_control( quiet = FALSE, + control = dsem_control( quiet = TRUE, run_model = TRUE, use_REML = TRUE, gmrf_parameterization = "projection", + constant_variance = c("conditional", "marginal", "diagonal")[2], # marginal not working map = Map, parameters = Params ) ) +#Rep = fitRTMB$obj$report() +range(fit$opt$par - fitRTMB$opt$par) + +if( FALSE ){ + Rep = fitRTMB$obj$report() + dgmrf( as.vector(Rep$z_tj), mu=as.vector(Rep$xhat_tj + Rep$delta_tj), Q=Rep$Q0_kk, log=TRUE ) + + example = list( + Gamma_kk = Rep$Gamma_kk, + Rho_kk = Rep$Rho_kk + ) + saveRDS( example, file=file.path(R'(C:\Users\James.Thorson\Desktop\Git\dsem\scratch)',"example.RDS")) + # Use in example2.R +} # obj = fitRTMB$obj