library(flagr); #Use the package "flagr"

#=========================Generate a data set

J = 20;
N = 1000;
K = 1;

A <- matrix(0.4, J, K);


S <- matrix(0, J, J);
for(i in 1:(J-1)){
  S[i,i+1] <- 0.5;
  S[i+1,i] <- 0.5;
}
diag(S) <- -4;

Q <- A!=0
supp <- S!=0 #supp specifies the graph

Sigma <- matrix(1, K, K);
response <- FLaG.sample(A, S, Sigma, N); #Sample N respondents from the model specified above.

#===========================Fit the model when the graph is given
#Choose some initial value for A and S. 

A0 <- matrix(0.2, J, K);
S0 <- diag(-4, J);
Sigma0 <- 1;

res.supp <- confFLaG.supp(A0, S0, Sigma0, Q, supp, response, tol = 1e-05);
res.supp$A;
res.supp$S;
res.supp$Sigma;
res.supp$pseudoLoglik;
hist(res.supp$scoreEst);

#============================Regularized estimator
gamma = 0.2/sqrt(N); #Gamma is the tuning parameter
res.reg <- confFLaG(A0, S0, Sigma0, Q, gamma,  response, tol = 1e-05)

res.reg$S
res.reg$A
res.reg$Sigma;

#=============================Use a sequence of tuning parameters 
#=============================and then choose the model by BIC. 

A0 <- matrix(0.2, J, K);
S0 <- diag(-4, J);
Sigma0 <- 1;

gamma.vec <- seq(0.05, 1, by = 0.05)/sqrt(N);
result.list <- list();
BIC <- rep(0, length(gamma.vec));

for(i in 1:length(gamma.vec)){
  print(i);
  print("======================================================")
  
  gamma <- gamma.vec[i]
  
  temp <- confFLaG(A0, S0, Sigma0, Q, gamma, response, tol = 1e-5, print_proc =F)
  
  supp.new <- abs(temp$S) > 1e-4;  
  temp2 <- confFLaG.supp(temp$A, temp$S, temp$Sigma, Q, supp.new, response, tol = 1e-5, print_proc =F)  
  result.list <- c(result.list, list(temp2));
  
  num.par <- sum(Q) + J + (K-1)*K/2 + sum(supp.new[lower.tri(supp.new, diag = F)])
  
  BIC[i] <- -2 * temp2$pseudoLoglik + log(N) * num.par;
}

ind <- which.min(BIC);
sel.model <- result.list[[ind]];

sel.model$A;
sel.model$S;
sel.model$pseudoLoglik;