library(MASS)
source('smgm.R')

## Rothman, Bickel, Levina and Zhu, EJS
## p dimension, ex: example number
## return: the Sigma matrix
RBLZ = function(p, ex=1)
{
  M = matrix(0,p,p)
  if(ex==0) {
    M = diag(p)
  }
  else if(ex==1) {
    for(i in 1:p)
      for(j in 1:p)
        M[i,j] = 0.7^(abs(i-j))
  }
  else if(ex==2) {
    diag(M) = 1
    for(i in 1:p)
      for(j in 1:p) {
        tmp = abs(i-j)
        if(tmp==1)
          M[i,j] = 0.4
        else if(tmp==2 | tmp==3)
          M[i,j] = 0.2
        else if(tmp==4)
          M[i,j] = 0.1
      }
    M = solve(M)
      } else {
    if(ex==3) {
      a = 0.1
    } else {
      a = 0.5 }
    M[upper.tri(M)] = 0.5*rbinom(p*(p-1)/2,1,a)
    M = M+t(M)
    diag(M) = 0
    eg = eigen(M)$values
    d = (eg[1]-p*eg[p])/(p-1)
    M=M+d*diag(p)
    M = solve(M)
  }
  M
}

n = 20
p = 20
q = 20
sig = RBLZ(p, ex=1)
psi = RBLZ(q, ex=1)

Sig.Psi = sig%x%psi

x = mvrnorm(n, rep(0,p*q), Psi.Sig)
xx = array(t(x), c(p,q,n))

fit = smgm(xx, scale=F,start='warm',method='scad')