update make_dsem_ram to use additive-kronecker

NAMESPACE

@@ -22,7 +22,10 @@ export(make_dsem_ram)
 export(parse_path)
 export(stepwise_selection)
 importFrom(Matrix,Cholesky)
+importFrom(Matrix,drop0)
+importFrom(Matrix,mat2triplet)
 importFrom(Matrix,solve)
+importFrom(Matrix,sparseMatrix)
 importFrom(TMB,MakeADFun)
 importFrom(TMB,compile)
 importFrom(TMB,dynlib)

R/dsem.R

@@ -25,7 +25,7 @@
 #'
 #' @importFrom TMB compile dynlib MakeADFun sdreport summary.sdreport
 #' @importFrom stats AIC sd .preformat.ts na.omit nlminb optimHess pnorm rnorm simulate time tsp<-
-#' @importFrom Matrix solve Cholesky
+#' @importFrom Matrix solve Cholesky sparseMatrix mat2triplet drop0
 #' @importFrom sem sem
 #' @importFrom igraph plot.igraph graph_from_data_frame with_sugiyama layout_
 #' @importFrom ggraph ggraph geom_edge_arc create_layout rectangle geom_node_text theme_graph

R/make_dsem_ram.R

@@ -332,22 +332,46 @@ function( sem,
   }
 
   # Loop through paths
+  G_kk = P_kk = sparseMatrix( i=1, j=1, x=0, dims=rep(length(variables)*length(times),2) )   # Make with a zero
+  #P_kk = new("dgCMatrix")
+  #P_kk = Matrix()
+  #P_kk@Dim <- as.integer(rep(length(variables)*length(times),2))
+  #P_kk@p = integer(length(variables)*length(times)+1L)
+  #G_kk = P_kk
   for( i in seq_len(nrow(model)) ){
-  for( t in seq_along(times) ){
     lag = as.numeric(model[i,2])
-    par.no = par.nos[i]
-    # Get index for "from"
-    from = c( times[t], model[i,'first'] )
-    from_index = match_row( Q_names, from )
-    from_index = ifelse( length(from_index)==0, NA, from_index )
-    # Get index for "to"
-    to = c( times[t+lag], model[i,'second'] )
-    to_index = match_row( Q_names, to )
-    to_index = ifelse( length(to_index)==0, NA, to_index )
-    ram_new = data.frame( "heads"=abs(as.numeric(model[i,'direction'])), "to"=to_index, "from"=from_index, "parameter"=par.no, "start"=startvalues[i] )
-    ram = rbind( ram, ram_new )
-  }}
-  rownames(ram) = NULL
+    L_tt = sparseMatrix( i = seq(lag+1,length(times)),
+                         j = seq(1,length(times)-lag),
+                         x = 1,
+                         dims = rep(length(times),2) )
+    P_jj = sparseMatrix( i = match(model[i,'second'],variables),
+                         j = match(model[i,'first'],variables),
+                         x = 1,
+                         dims = rep(length(variables),2) )
+    tmp_kk = kronecker(P_jj, L_tt)
+    if(abs(as.numeric(model[i,'direction']))==1){
+      P_kk = P_kk + tmp_kk * par.nos[i]
+    }else{
+      G_kk = G_kk + tmp_kk * par.nos[i]
+    }
+    #for( t in seq_along(times) ){
+    #  # Get index for "from"
+    #  from = c( times[t], model[i,'first'] )
+    #  from_index = match_row( Q_names, from )
+    #  from_index = ifelse( length(from_index)==0, NA, from_index )
+    #  # Get index for "to"
+    #  to = c( times[t+lag], model[i,'second'] )
+    #  to_index = match_row( Q_names, to )
+    #  to_index = ifelse( length(to_index)==0, NA, to_index )
+    #  ram_new = data.frame( "heads"=abs(as.numeric(model[i,'direction'])), "to"=to_index, "from"=from_index, "parameter"=par.no, "start"=startvalues[i] )
+    #  ram = rbind( ram, ram_new )
+    #}
+  }
+  #rownames(ram) = NULL
+  ram = rbind( cbind(1, sapply(mat2triplet(drop0(P_kk)),cbind)),
+               cbind(2, sapply(mat2triplet(drop0(G_kk)),cbind)) )
+  ram = data.frame( ram, startvalues[ram[,4]] )
+  colnames(ram) = c("heads", "to", "from", "parameter", "start")
 
   #
   if( isTRUE(remove_na) ){

tests/testthat/test-platform.R

@@ -3,23 +3,23 @@ test_that("dsem example is working ", {
   #skip_on_ci()
   sem = "
     Profits -> Consumption, 0, a2
-    Profits -> Consumption, -1, a3
+    Profits -> Consumption, 1, a3
     Priv_wage -> Consumption, 0, a4
     Gov_wage -> Consumption, 0, a4
     Consumption <-> Consumption, 0, v1
-    Consumption -> Consumption, -1, ar1
-    Consumption -> Consumption, -2, ar2
+    Consumption -> Consumption, 1, ar1
+    Consumption -> Consumption, 2, ar2
     Profits -> Investment, 0, b2
-    Profits -> Investment, -1, b3
-    Capital_stock -> Investment, -1, b4
+    Profits -> Investment, 1, b3
+    Capital_stock -> Investment, 1, b4
     Investment <-> Investment, 0, v2
     neg_Gov_wage <-> neg_Gov_wage, 0, v3
     GNP -> Priv_wage, 0, c2
     Taxes -> Priv_wage, 0, c2
     neg_Gov_wage -> Priv_wage, 0, c2
-    GNP -> Priv_wage, -1, c3
-    Taxes -> Priv_wage, -1, c3
-    neg_Gov_wage -> Priv_wage, -1, c3
+    GNP -> Priv_wage, 1, c3
+    Taxes -> Priv_wage, 1, c3
+    neg_Gov_wage -> Priv_wage, 1, c3
     Time -> Priv_wage, 0, c4
     Priv_wage <-> Priv_wage, 0, v4
     GNP <-> GNP, 0, v5