Projthis.Rmd

---
title: "Projthis"
output: github_document
---

# Create Projthis Project
```{r}
projthis::proj_create("/Users/bropc/Documents/LMU/Master Statistics and Data Science/Masterarbeit/R Master/StARS_Simulations")
```

Load Projthis
```{r}
library("projthis")
```

Create Workflow
```{r}
proj_use_workflow("workflow", git_ignore_data = FALSE)
```

Create a workflow project
```{r}
# creates a new rmd file from a template
proj_workflow_use_rmd("00-import")
```

```{r}
install.packages("C:/R Projekte/QUIC_1.1.1.tar.gz", repos = NULL, type="source")
install.packages("/Users/bropc/Documents/LMU/Master Statistics and Data Science/Masterarbeit/R Master/QUIC_1.1.1.tar", repos = NULL, type="source")
```

```{r}
# Load necessary libraries
library(pulsar)
library(huge)
library(QUIC)
library(orca)
library(ggplot2)


# Specify configurations
configs <- list(
  list(n=800, p=40),
  list(n=400, p=100),
  list(n=200, p=200)
)


# Function to run experiment
run_hub <- function(n, p, s, J, rho, u, b, N) {
  
  # 1. Set Parameters
  # Size of each hub group
  s <- 20
  # Number of hub groups
  J <- floor(p/s) 
  # Off-diagonal elements
  rho <- (1/(s+1))
  # Diagonal elements
  u <- 0.001
  # Subsampled dataset size 
  b = (floor(10*sqrt(n)))/n
  # Number of Subsamples N
  N = 20
  
  # 2. Data Generation
  set.seed(123)
  hub <- huge.generator(n = n, d = p, graph = "hub", v = rho, u = u, g = J, vis = TRUE)
  hub_data <- hub$data
  
  true_graph <- hub$theta != 0
  true_cov <- hub$sigma
  lambda_path  <- getLamPath(max = getMaxCov(hub_data), min = 0.01, len = 40)
  lambda <- list(lambda=lambda_path)
  
  # 3. Run QUIC method
  quicr <- function(hub_data, lambda, ...) {
  S <- cov(hub_data)
  est <- QUIC(S, rho = 1, path = lambda, msg = 0, tol = 1e-2, ...)
  est$path <- lapply(seq(length(lambda)), function(i) {
  ## convert precision array to adj list
  tmp <- est$X[,,i]; diag(tmp) <- 0
  as(tmp!=0, "lMatrix")
  })
  est
  }

  # Run pulsar package
  out.p <- pulsar(
      hub_data,
      fun = quicr,
      fargs = lambda,
      criterion = c('stars', 'gcd'),
      thresh = 0.1,
      subsample.ratio = b,
      rep.num = N,
      seed = NULL,
      lb.stars = TRUE,
      ub.stars = TRUE,
      ncores = 1,
      refit = FALSE
  )
  # Get optimal index for gcd
  lam_gam <- get.opt.index(out.p, criterion = "gcd")
  # Set optimal index for gcd
  opt.index(out.p, criterion = "gcd") <- lam_gam

  fit  <- refit(out.p, criterion = c("stars", "gcd"))
  stars_graph <- fit[["refit"]][["stars"]]
  gstars_graph <- fit[["refit"]][["gcd"]]

  # 4. Run Oracle QUIC
  oracle_quic <- function(true_cov, lambda, ...) {
    S <- true_cov
    est <- QUIC(S, rho = 1, path = lambda, msg = 0, tol = 1e-2, ...)
    est$path <- lapply(seq(length(lambda)), function(i) {
        tmp <- est$X[,,i]
        diag(tmp) <- 0
        as(tmp != 0, "lMatrix")
    })
    est
  }
  
  oracle_results <- oracle_quic(true_cov, lambda_path)

  # Minimize total number of different edges between the estimated and true graph
  best_lambda_index <- which.min(sapply(1:length(lambda_path), function(i) {
    estimated_graph <- oracle_results$path[[i]]
    sum(estimated_graph != true_graph)
  }))

  best_lambda <- lambda_path[best_lambda_index]

  # Extract the oracle precision matrix for the best lambda:
  oracle_graph <- oracle_results$X[,,best_lambda_index]

  
  # 5. Calculate metrics
 
  # List of estimated graphs
  estimated_graphs <- list(stars = stars_graph, gstars = gstars_graph, oracle = oracle_graph)

  # Create empty lists to store results
  F1_scores <- list()
  jaccard_indices <- list()

  # Loop through each estimated graph and calculate metrics
  for (name in names(estimated_graphs)) {
  
  # Convert the matrices to binary (1 for edge, 0 for no edge)
  est_bin <- as.numeric(estimated_graphs[[name]] != 0)
  true_bin <- as.numeric(true_graph != 0)

  # Compute true positives, false positives, true negatives, and false negatives
  TP <- sum(est_bin == 1 & true_bin == 1)
  FP <- sum(est_bin == 1 & true_bin == 0)
  TN <- sum(est_bin == 0 & true_bin == 0)
  FN <- sum(est_bin == 0 & true_bin == 1)

  # Calculate precision, recall, F1-score, and Jaccard index
  precision <- TP / (TP + FP)
  recall <- TP / (TP + FN)
  F1 <- 2 * (precision * recall) / (precision + recall)
  jaccard_index <- TP / (TP + FP + FN)

  # Save results to lists
  F1_scores[[name]] <- F1
  jaccard_indices[[name]] <- jaccard_index
  }
  
  return(list(F1_stars = F1_scores[["stars"]], F1_gcd = F1_scores[["gstars"]], F1_oracle = F1_scores[["oracle"]]))

}


# Run the loop for each configuration
results <- list()
for (config in configs) {
  cat("Running simulation for configuration:", paste0("n=", config$n, " p=", config$p), "\n")
  res <- run_hub(config$n, config$p, s, J, rho, u, b, N)
  results[[paste0("n=", config$n, " p=", config$p)]] <- res
}

# Create a data frame to store results for plotting
df <- data.frame(np_setting = names(results), 
                 F1_stars = sapply(results, function(x) x$F1_stars),
                 F1_gcd = sapply(results, function(x) x$F1_gcd),
                 F1_oracle = sapply(results, function(x) x$F1_oracle))




```


### 1) Hub1: for n = 800 and p = 40
```{r}
library(huge)
library(pulsar)

#Set Parameters
n <- 800 
p <- 40 
#Size of each hub group
s <- 20
#Number of hub groups
J <- floor(p/s) 
#Off-diagonal elements
rho <- (1/(s+1))
#Diagonal elements
u <- 0.003
#Subsampled dataset size 
b = (floor(10*sqrt(n)))/n
#Number of Subsamples N
N = 20


set.seed(123)
hub <- huge.generator(n = n, d = p, graph = "hub", v = rho, u = u, g = J, vis = TRUE)
hub_data <- hub$data

#Lambda path
lambda_path  <- getLamPath(max = getMaxCov(hub_data), min = 0.01, len = 40)
lambda <- list(lambda=lambda_path)

true_graph <- as(hub$theta, "lMatrix")
true_cov <- hub$sigma

#Plot Adjacency Matrix
huge.plot(hub$theta)

true_graph <- as.matrix(hub$theta)
test <- as(hub$theta, "TsparseMatrix")
```

```{r}
#Define QUIC as method of choice for stars and gstars
library(QUIC)

quicr <- function(hub_data, lambda, ...) {
S <- cov(hub_data)
est <- QUIC(S, rho = 1, path = lambda, msg = 0, tol = 1e-2, ...)
est$path <- lapply(seq(length(lambda)), function(i) {
## convert precision array to adj list
tmp <- est$X[,,i]; diag(tmp) <- 0
as(tmp!=0, "lMatrix") 
})
est
}
```


## StARS, G-StARS and Oracle
```{r}

#Define QUIC as method of choice for stars and gstars
library(QUIC)

quicr <- function(hub_data, lambda, ...) {
  S <- cov(hub_data)
  est <- QUIC(S, rho = 1, path = lambda, msg = 0, tol = 1e-2, ...)
  est$path <- lapply(seq(length(lambda)), function(i) {
  ## convert precision array to adj list
  tmp <- est$X[,,i]; diag(tmp) <- 0
  tmp <- ifelse(tmp != 0, 1, 0)
  return(tmp)
  })
  est
}

#Run pulsar package
library(orca)
out.p <- pulsar(
      hub_data,
      fun = quicr,
      fargs = (lambda = lambda),
      criterion = c('stars', 'gcd'),
      thresh = 0.1,
      subsample.ratio = b,
      rep.num = N,
      seed = NULL,
      lb.stars = TRUE,
      ub.stars = TRUE,
      ncores = 1,
      refit = FALSE
)
#Get optimal index for gcd
lam_gam <- get.opt.index(out.p, criterion = "gcd")
#Set optimal index for gcd
opt.index(out.p, criterion = "gcd") <- lam_gam


out.p
plot(out.p, scale = T, invlam = FALSE, loglam = FALSE, legends = FALSE)

fit  <- refit(out.p, criterion = c("stars", "gcd"))
print(fit)

stars_graph <- fit[["refit"]][["stars"]]
gstars_graph <- fit[["refit"]][["gcd"]]


#Define Hamming Distance as criterium for Oracle
hamming_distance <- function(matrix1, matrix2) {
    return(sum(matrix1 != matrix2))
}

# Oracle with QUIC
oracle_quic <- function(data, lambda, ...) {
    S <- true_cov
    est <- QUIC(S, rho = 1, path = lambda, msg = 0, tol = 1e-2, ...)
    est$path <- lapply(seq(length(lambda)), function(i) {
    ## convert precision array to adj list
    tmp <- est$X[,,i]; diag(tmp) <- 0
    tmp <- ifelse(tmp != 0, 1, 0)
    return(tmp)
    })
    est
}

oracle_results <- oracle_quic(hub_data, lambda_path)

# Minimize Hamming distance between the estimated and true graph HAMMING!
# best_lambda_index is position of optimal lambda on path
best_lambda_index <- which.min(sapply(1:length(lambda_path), function(i) {
    estimated_graph <- oracle_results$path[[i]]
    hamming_distance(estimated_graph, true_graph)
}))

#Determine best lambda oracle
best_lambda <- lambda_path[best_lambda_index]

#Extract the oracle precision matrix for the best lambda:
oracle_graph <- oracle_results$X[,,best_lambda_index]
# Convert precision matrix to adjacency matrix
oracle_graph <- ifelse(oracle_graph != 0, 1, 0)
diag(oracle_graph) <- 0


#Plot
huge.plot(stars_graph)
huge.plot(gstars_graph)
huge.plot(oracle_graph)
```


## Oracle
```{r}

#Subsample original data set (b x p)

# Function to generate N subsamples from the original dataset
subsample_data <- function(data, N, b) {
  n <- nrow(data)
  subsample_size <- floor(b * n)
  
  subsamples <- vector("list", N)
  
  for (i in 1:N) {
    subsample_indices <- sample(1:n, subsample_size)
    subsamples[[i]] <- data[subsample_indices, ]
  }
  
  return(subsamples)
}

# Use the function to generate the subsamples
oracle_data <- subsample_data(hub_data, N, b)


#Define Hamming Distance as criterium for Oracle
hamming_distance <- function(matrix1, matrix2) {
    return(sum(matrix1 != matrix2))
}

# Assuming lambda_path is defined somewhere above in your original code

# Oracle with QUIC function adjusted to take a single data sample
oracle_quic <- function(data, lambda, ...) {
    S <- cov(data)
    est <- QUIC(S, rho = 1, path = lambda, msg = 0, tol = 1e-2, ...)
    est$path <- lapply(seq(length(lambda)), function(i) {
    ## convert precision array to adj list
    tmp <- est$X[,,i]; diag(tmp) <- 0
    tmp <- ifelse(tmp != 0, 1, 0)
    return(tmp)
    })
    est
}

# List to store best lambda for each subsample
best_lambdas <- numeric(N)

# Iterate over each subsample
for(i in 1:N) {
    # Apply the QUIC method to the current subsample
    oracle_results <- oracle_quic(oracle_data[[i]], lambda_path)
    
    # Determine the best lambda for the current subsample
    best_lambda_index <- which.min(sapply(1:length(lambda_path), function(j) {
        estimated_graph <- oracle_results$path[[j]]
        hamming_distance(estimated_graph, true_graph)
    }))
    
    # Store the best lambda for this subsample
    best_lambdas[i] <- lambda_path[best_lambda_index]
}

# Determine the most frequently optimal lambda across all subsamples
best_lambda_overall <- as.numeric(names(sort(table(best_lambdas), decreasing = TRUE)[1]))

# Using the overall best lambda to run QUIC on the full dataset
oracle_results_full <- oracle_quic(hub_data, lambda_path)

# Find the index of the lambda in lambda_path that is closest to best_lambda_overall
closest_lambda_index <- which.min(abs(lambda_path - best_lambda_overall))

# Extract the precision matrix for the closest lambda to best_lambda_overall
oracle_graph <- oracle_results_full$X[,,closest_lambda_index]

# Convert precision matrix to adjacency matrix
oracle_graph <- ifelse(oracle_graph != 0, 1, 0)
diag(oracle_graph) <- 0

#Plot
huge.plot(oracle_graph)

```

```{r}
#out.p$stars$merge

#res10 = quicr(hub_data, c(out.p$stars$summary[19]))
```

## ACHTUNG: Sicherheitssafe für hub_generator
```{r}
library(MASS)
library(pulsar)
library(Matrix)
library(igraph)
library(huge)

generator <- function (n, p, graph, rho, g, vis, verbose) 
{
    gcinfo(FALSE)
    if (verbose) 
        cat("Generating data from the multivariate normal distribution with the", 
            graph, "graph structure....")
  
  g.large = p%%g #Rest Funktion
    g.small = g - g.large
    n.small = floor(p/g)
    n.large = n.small + 1
    g.list = c(rep(n.small, g.small), rep(n.large, g.large))
    g.ind = rep(c(1:g), g.list)
    rm(g.large, g.small, n.small, n.large, g.list)
    gc()
    
    theta = matrix(0, p, p) #Here Theta defined (pxp matrix with entries "0")
  
     if (graph == "hub") { #Define Hub connections
        for (i in 1:g) {
            tmp = which(g.ind == i)
            theta[tmp[1], tmp] = 1
            theta[tmp, tmp[1]] = 1
            rm(tmp)
            gc()
        }
    }
  
    diag(theta) = 0
    omega = theta * rho
    diag(omega) = abs(min(eigen(omega)$values)) + 0.1 #Ensure that the omega is positive definite by adding a small positive value
    sigma = cov2cor(solve(omega)) #Give Cov and transform it into Corr
    #omega = solve(sigma) 
    diag(omega) = 1 #Set diagonal of precision matrix to 1 (Liu et al.)
    x = mvrnorm(n, rep(0, p), sigma)
    sigmahat = cor(x)
    if (vis == TRUE) {
        fullfig = par(mfrow = c(2, 2), pty = "s", omi = c(0.3, 
            0.3, 0.3, 0.3), mai = c(0.3, 0.3, 0.3, 0.3))
        fullfig[1] = image(theta, col = gray.colors(256), main = "Adjacency Matrix")
        fullfig[2] = image(sigma, col = gray.colors(256), main = "Covariance Matrix")
        g = graph.adjacency(theta, mode = "undirected", diag = FALSE)
        layout.grid = layout.fruchterman.reingold(g)
        fullfig[3] = plot(g, layout = layout.grid, edge.color = "gray50", 
            vertex.color = "red", vertex.size = 3, vertex.label = NA, 
            main = "Graph Pattern")
        fullfig[4] = image(sigmahat, col = gray.colors(256), 
            main = "Empirical Matrix")
        rm(fullfig, g, layout.grid)
        gc()
    }
    if (verbose) 
        cat("done.\n")
    rm(vis, verbose)
    gc()
    sim = list(data = x, sigma = sigma, sigmahat = sigmahat, 
        omega = omega, theta = Matrix(theta, sparse = TRUE), 
        sparsity = sum(theta)/(p * (p - 1)), graph.type = graph)
    class(sim) = "sim"
    return(sim)
}
```


## Random Schrott
```{r}

generator_rand <- function (n, p, prob, vis = FALSE, verbose = TRUE) 
{
    gcinfo(FALSE)
    if (verbose) 
        cat("Generating data from the multivariate normal 
            distribution with the random graph structure....")
  
    prob = sqrt(prob/2) * (prob < 0.5) + 
      (1 - sqrt(0.5 - 0.5 * prob)) * (prob >= 0.5)
    
    theta = matrix(0, p, p)
        tmp = matrix(runif(p^2, 0, 0.5), p, p) #set threshold
        tmp = tmp + t(tmp)
        theta[tmp < prob] = 1 #thresholding: If prob greater than 0.5 than edge is present (1)
        rm(tmp)
        gc()
        
    diag(theta) = 0 #To ensure that there are no self-loops in the graph
    omega = theta * 0.3
    diag(omega) = abs(min(eigen(omega)$values)) + 0.1 #!!!
    sigma = cov2cor(solve(omega)) #!!!
    omega = solve(sigma)
    x = mvrnorm(n, rep(0, p), sigma)
    sigmahat = cor(x) #!!!
    
    if (vis == TRUE) {
        fullfig = par(mfrow = c(2, 2), pty = "s", omi = c(0.3, 
            0.3, 0.3, 0.3), mai = c(0.3, 0.3, 0.3, 0.3))
        fullfig[1] = image(theta, col = gray.colors(256), main = "Adjacency Matrix")
        fullfig[2] = image(sigma, col = gray.colors(256), main = "Covariance Matrix")
        g = graph.adjacency(theta, mode = "undirected", diag = FALSE)
        layout.grid = layout.fruchterman.reingold(g)
        fullfig[3] = plot(g, layout = layout.grid, edge.color = "gray50", 
            vertex.color = "red", vertex.size = 3, vertex.label = NA, 
            main = "Graph Pattern")
        fullfig[4] = image(sigmahat, col = gray.colors(256), 
            main = "Empirical Matrix")
        rm(fullfig, g, layout.grid)
        gc()
    }
    
    if (verbose) 
        cat("done.\n")
    rm(vis, verbose)
    gc()
    
    sim = list(data = x, sigma = sigma, sigmahat = sigmahat, 
    omega = omega, theta = Matrix(theta, sparse = TRUE), 
    sparsity = sum(theta)/(p * (p - 1))) #Sparsity level = number of present edges divided by number of potential edges p x p.
    class(sim) = "sim"
    return(sim)
}
```

```{r}

generator_rand <- function (n, p, prob, vis = FALSE, verbose = TRUE) 
{
    gcinfo(FALSE)
    if (verbose) 
        cat("Generating data from the multivariate normal 
            distribution with the random graph structure....")
  
    prob = sqrt(prob/2) * (prob < 0.5) + 
      (1 - sqrt(0.5 - 0.5 * prob)) * (prob >= 0.5)
    
    theta = matrix(0, p, p)
        tmp = matrix(runif(p^2, 0, 0.5), p, p) #set threshold
        tmp = tmp + t(tmp)
        theta[tmp < prob] = 1 #thresholding: If prob greater than 0.5 than edge is present (1)
        rm(tmp)
        gc()
        
    diag(theta) = 0 #To ensure that there are no self-loops in the graph
    omega = theta
    diag(omega) = abs(min(eigen(omega)$values)) + 0.1  #!!!
    sigma = solve(omega) #!!!
    omega = solve(sigma)
    x = mvrnorm(n, rep(0, p), sigma)
    sigmahat = cov(x) #!!!
    
    if (vis == TRUE) {
        fullfig = par(mfrow = c(2, 2), pty = "s", omi = c(0.3, 
            0.3, 0.3, 0.3), mai = c(0.3, 0.3, 0.3, 0.3))
        fullfig[1] = image(theta, col = gray.colors(256), main = "Adjacency Matrix")
        fullfig[2] = image(sigma, col = gray.colors(256), main = "Covariance Matrix")
        g = graph.adjacency(theta, mode = "undirected", diag = FALSE)
        layout.grid = layout.fruchterman.reingold(g)
        fullfig[3] = plot(g, layout = layout.grid, edge.color = "gray50", 
            vertex.color = "red", vertex.size = 3, vertex.label = NA, 
            main = "Graph Pattern")
        fullfig[4] = image(sigmahat, col = gray.colors(256), 
            main = "Empirical Matrix")
        rm(fullfig, g, layout.grid)
        gc()
    }
    
    if (verbose) 
        cat("done.\n")
    rm(vis, verbose)
    gc()
    
    sim = list(data = x, sigma = sigma, sigmahat = sigmahat, 
    omega = omega, theta = Matrix(theta, sparse = TRUE), 
    sparsity = sum(theta)/(p * (p - 1))) #Sparsity level = number of present edges divided by number of potential edges p x p.
    class(sim) = "sim"
    return(sim)
}
```

```{r}

#Set Parameters
n <- 800 
p <- 40 
prob <- 3/p

generator_rand <- function (n, p, prob, vis = FALSE, verbose = TRUE) 
{
    gcinfo(FALSE)
    if (verbose) 
        cat("Generating data from the multivariate normal 
            distribution with the random graph structure....")

    prob = sqrt(prob/2) * (prob < 0.5) + 
      (1 - sqrt(0.5 - 0.5 * prob)) * (prob >= 0.5)
    
    theta = matrix(0, p, p)
        tmp = matrix(runif(p^2, 0, 0.5), p, p) #set threshold
        tmp = tmp + t(tmp)
        theta[tmp < prob] = 1 #thresholding: If prob greater than 0.5 than edge is present (1)
        rm(tmp)
        gc()
        
    diag(theta) = 0 #To ensure that there are no self-loops in the graph
    
    random_matrix <- matrix(runif(p^2, -0.5, 0.5), p, p)
    omega = theta + (random_matrix + t(random_matrix)) / 2
    
    # Compute eigenvalues
    eigenvalues = eigen(omega, symmetric = TRUE)$values

    # If the smallest eigenvalue is negative, add its absolute value plus a small constant to the diagonal
    if (min(eigenvalues) < 0) {
    diag(omega) = diag(omega) - min(eigenvalues) + 0.1
    }
    
    #diag(omega) = abs(min(eigen(omega)$values)) + 0.1  #!!!
    sigma = solve(omega) #!!!
    omega = solve(sigma)
    x = mvrnorm(n, rep(0, p), sigma)
    sigmahat = cov(x) #!!!
    
    
eigenvalues    
prob
theta
omega
range(omega)
str(omega)
sigma
x
```

```{r}

#Set Parameters
n <- 800 
p <- 40 
prob <- 3/p

    prob = sqrt(prob/2) * (prob < 0.5) + 
      (1 - sqrt(0.5 - 0.5 * prob)) * (prob >= 0.5)
    
    theta = matrix(0, p, p)
        tmp = matrix(runif(p^2, 0, 0.5), p, p) #set threshold
        tmp = tmp + t(tmp)
        theta[tmp < prob] = 1 #thresholding: If prob greater than 0.5 than edge is present (1)
        rm(tmp)
        gc()
        
    diag(theta) = 0 #To ensure that there are no self-loops in the graph
    
    random_matrix <- matrix(runif(p^2, -0.5, 0.5), p, p)
    omega = theta + (random_matrix + t(random_matrix)) / 2
    
    # Compute eigenvalues
    eigenvalues = eigen(omega, symmetric = TRUE)$values

    # If the smallest eigenvalue is negative, add its absolute value plus a small constant to the diagonal
    if (min(eigenvalues) < 0) {
    diag(omega) = diag(omega) - min(eigenvalues) + 0.1
    }
    
    #diag(omega) = abs(min(eigen(omega)$values)) + 0.1  #!!!
    sigma = solve(omega) #!!!
    omega = solve(sigma)
    x = mvrnorm(n, rep(0, p), sigma)
    sigmahat = cov(x) #!!!


if (vis == TRUE) {
        fullfig = par(mfrow = c(2, 2), pty = "s", omi = c(0.3, 
            0.3, 0.3, 0.3), mai = c(0.3, 0.3, 0.3, 0.3))
        fullfig[1] = image(theta, col = gray.colors(256), main = "Adjacency Matrix")
        fullfig[2] = image(sigma, col = gray.colors(256), main = "Covariance Matrix")
        g = graph.adjacency(theta, mode = "undirected", diag = FALSE)
        layout.grid = layout.fruchterman.reingold(g)
        fullfig[3] = plot(g, layout = layout.grid, edge.color = "gray50", 
            vertex.color = "red", vertex.size = 3, vertex.label = NA, 
            main = "Graph Pattern")
        fullfig[4] = image(sigmahat, col = gray.colors(256), 
            main = "Empirical Matrix")
        rm(fullfig, g, layout.grid)
        gc()
    }
    
    if (verbose) 
        cat("done.\n")
    rm(vis, verbose)
    gc()
    
    sim = list(data = x, sigma = sigma, sigmahat = sigmahat, 
    omega = omega, theta = Matrix(theta, sparse = TRUE), 
    sparsity = sum(theta)/(p * (p - 1))) #Sparsity level = number of present edges divided by number of potential edges p x p.
    class(sim) = "sim"
    return(sim)
}


#Set Parameters
n <- 800 
p <- 40 

#Subsampled dataset size 
b = (floor(10*sqrt(n)))/n
#Number of Subsamples N
N = 20

prob <- 3/p #For p = 40, Prob = 0.075

ER <- generator_rand(n = n, p = p, prob = prob, vis = TRUE)
ER_data <- ER$data

#ER_data 
ER$omega
range(ER$omega)
ER$theta
ER$sparsity

huge.plot(ER$theta)
```






```{r}

#Set Parameters
n <- 800 
p <- 40 

#Subsampled dataset size 
b = (floor(10*sqrt(n)))/n
#Number of Subsamples N
N = 20

prob <- 3/p #For p = 40, Prob = 0.075

ER <- generator_rand(n = n, p = p, prob = prob, vis = TRUE)
ER_data <- ER$data

#ER_data 
ER$omega
range(ER$omega)
ER$theta
ER$sparsity

huge.plot(ER$theta)


```


```{r}
library(MASS)
library(pulsar)
library(Matrix)
library(igraph)
library(huge)

generator <- function (n, p, graph, rho, g, vis, verbose) 
{
    gcinfo(FALSE)
    if (verbose) 
        cat("Generating data from the multivariate normal distribution with the", 
            graph, "graph structure....")
  
  g.large = p%%g #Rest Funktion
    g.small = g - g.large
    n.small = floor(p/g)
    n.large = n.small + 1
    g.list = c(rep(n.small, g.small), rep(n.large, g.large))
    g.ind = rep(c(1:g), g.list)
    rm(g.large, g.small, n.small, n.large, g.list)
    gc()
    
    theta = matrix(0, p, p) #Here Theta defined (pxp matrix with entries "0")
  
     if (graph == "hub") { #Define Hub connections
        for (i in 1:g) {
            tmp = which(g.ind == i)
            theta[tmp[1], tmp] = 1
            theta[tmp, tmp[1]] = 1
            rm(tmp)
            gc()
        }
    }
  
    diag(theta) = 0
    omega = theta * rho
    #diag(omega) <-  abs(min(eigen(omega)$values)) + 0.1
    diag(omega) <- 1
    #sigma = cov2cor(solve(omega)) #
    sigma = solve(omega) 
    omega = solve(sigma) 
    x = mvrnorm(n, rep(0, p), sigma)
    sigmahat = cov(x)
    if (vis == TRUE) {
        fullfig = par(mfrow = c(2, 2), pty = "s", omi = c(0.3, 
            0.3, 0.3, 0.3), mai = c(0.3, 0.3, 0.3, 0.3))
        fullfig[1] = image(theta, col = gray.colors(256), main = "Adjacency Matrix")
        fullfig[2] = image(sigma, col = gray.colors(256), main = "Covariance Matrix")
        g = graph.adjacency(theta, mode = "undirected", diag = FALSE)
        layout.grid = layout.fruchterman.reingold(g)
        fullfig[3] = plot(g, layout = layout.grid, edge.color = "gray50", 
            vertex.color = "red", vertex.size = 3, vertex.label = NA, 
            main = "Graph Pattern")
        fullfig[4] = image(sigmahat, col = gray.colors(256), 
            main = "Empirical Matrix")
        rm(fullfig, g, layout.grid)
        gc()
    }
    if (verbose) 
        cat("done.\n")
    rm(vis, verbose)
    gc()
    sim = list(data = x, sigma = sigma, sigmahat = sigmahat, 
        omega = omega, theta = Matrix(theta, sparse = TRUE), 
        sparsity = sum(theta)/(p * (p - 1)), graph.type = graph)
    class(sim) = "sim"
    return(sim)
}
```


```{r}
library(MASS)
library(pulsar)
library(Matrix)
library(igraph)
library(huge)

generator_rand <- function (n, p, prob, vis = FALSE, verbose = TRUE) 
{
    gcinfo(FALSE)
    if (verbose) 
        cat("Generating data from the multivariate normal 
            distribution with the random graph structure....")

    #prob = sqrt(prob/2) * (prob < 0.5) + 
     # (1 - sqrt(0.5 - 0.5 * prob)) * (prob >= 0.5)
    
    theta = matrix(0, p, p)
        tmp = matrix(runif(p^2, 0, 0.5), p, p) #set threshold
        tmp = tmp + t(tmp)
        theta[tmp < prob] = 1 #thresholding: If prob greater than 0.5 than edge is present (1)
        rm(tmp)
        gc()
    
    # Ensure no self-loops
    diag(theta) = 0 
    # Create a symmetric matrix with random values in the range [-1, 1] for off-diagonals
    random_matrix <- matrix(runif(p^2, -1, 1), p, p) # increase range
    omega = theta * (random_matrix + t(random_matrix)) / 2
    eigen_omega = eigen(omega) # Calculate eigenvalues of omega
    adjustment_factor = 1 # Initialize an adjustment factor
    current_condition_number = kappa(omega) # Calculate current condition number
    desired_condition_number = 100 # max(eigenvalue)/ min(eigenvalue) = 10
    
    #diagonal to 1, then off-diagonal between -1, 1 and then we rerun it until positve definite this!!!
    # set diagonal based on condition number let the rest as it is (diagonal greater than one)
    #matrix, change diagonal such that condition = 10, 100 (newton-raphson alg)

    # Adjust the diagonal elements
    while (current_condition_number > desired_condition_number) {
      # Adjust the smallest eigenvalue to change the condition number
      eigen_omega$values = pmax(eigen_omega$values, adjustment_factor * max(eigen_omega$values) / desired_condition_number)
      omega = eigen_omega$vectors %*% diag(eigen_omega$values) %*% t(eigen_omega$vectors) # Reconstruct omega with the adj eigenvalues
      current_condition_number = kappa(omega) # Recalculate condition number
      adjustment_factor = adjustment_factor * 1.01 # slightly increase the adjustment_factor to ensure the loop will terminate
    }
    
    omega = abs(omega) * (1 - diag(nrow(omega))) + diag(nrow(omega))
    
    # Inverse of omega to get sigma (covariance matrix)
    sigma = solve(omega) 
    # Generate multivariate normal samples
    x = MASS::mvrnorm(n, mu = rep(0, p), Sigma = sigma)
    # Estimate the covariance matrix from the generated samples
    sigmahat = cov(x)

  if (vis == TRUE) {
          fullfig = par(mfrow = c(2, 2), pty = "s", omi = c(0.3, 
              0.3, 0.3, 0.3), mai = c(0.3, 0.3, 0.3, 0.3))
          fullfig[1] = image(theta, col = gray.colors(256), main = "Adjacency Matrix")
          fullfig[2] = image(sigma, col = gray.colors(256), main = "Covariance Matrix")
          g = graph.adjacency(theta, mode = "undirected", diag = FALSE)
          layout.grid = layout.fruchterman.reingold(g)
          fullfig[3] = plot(g, layout = layout.grid, edge.color = "gray50", 
              vertex.color = "red", vertex.size = 3, vertex.label = NA, 
              main = "Graph Pattern")
          fullfig[4] = image(sigmahat, col = gray.colors(256), 
              main = "Empirical Matrix")
          rm(fullfig, g, layout.grid)
          gc()
    }
    
    if (verbose) 
        cat("done.\n")
    rm(vis, verbose)
    gc()
    
    sim = list(data = x, sigma = sigma, sigmahat = sigmahat, 
    omega = omega, theta = Matrix(theta, sparse = FALSE), #Matrix configuration
    sparsity = sum(theta)/(p * (p - 1))) #Sparsity level = number of present edges divided by number of potential edges p x p.
    class(sim) = "sim"
    return(sim)
}

```

# To Do List
## Implement sparsity levels in plot
## Check confidence intervals
## Relocate settings file

# 
```{r}

'
Oracles Objective: In your Oracle method, the goal is to choose the lambda value that minimizes the discrepancy (via Hamming distance) between the estimated and true graphs. This method inherently treats edges and non-edges equally.

High-Dimensional Settings: In high-dimensional settings (like n=200,p=200), where the number of potential edges (p(p−1)/2) is very large and the true graph is sparse, the majority of the entries in the adjacency matrix are zeros (non-edges). Thus, a method that equally weighs edges and non-edges might be biased towards predicting non-edges.
'

#Define Hamming Distance as criterium for Oracle
hamming_distance <- function(matrix1, matrix2) {
    return(sum(matrix1[lower.tri(matrix1)] != matrix2[lower.tri(matrix2)]))
}

# Minimize total number of different edges between the estimated and true graph
best_lambda_index <- which.min(sapply(1:length(lambda_path), function(i) {
  estimated_graph <- oracle_results$path[[i]]
  hamming_distance(estimated_graph, true_graph)
}))

```

```{r}

library(MASS)
library(Matrix)
library(igraph)
library(huge)

generator_geom <- function (n, p, rho, vis, verbose) 
{
    gcinfo(FALSE)
    if (verbose) 
        cat("Generating data from the multivariate normal distribution with the neighboorhod graph structure....")
  
    prob_function <- function(distance) exp(-4 * distance^2) / sqrt(2 * pi)

    # Step 1: Sample points from a unit square
    y <- matrix(runif(n * 2), ncol = 2)

    # Step 2: Initialize omega with diagonal elements set to 1
    omega <- diag(1, p)

    # Step 3: Populate off-diagonal elements of omega
    for (i in 1:(p-1)) {
      for (j in (i+1):p) {
        
        # Checks current row have fewer non-zero off-diagonal elements than the allowed
        if (sum(omega[i, -i] != 0) < floor(1/rho) && sum(omega[-j, j] != 0) < floor(1/rho)) { 
          
          # Euclidean distance between points in the unit square.
          distance <- sqrt(sum((y[i, ] - y[j, ])^2)) 
          
          # Randomly determines place an edge (non-zero value) between nodes 
          if (runif(1) < prob_function(distance)) {
            omega[i, j] <- omega[j, i] <- rho
          }
        }
      }
    }
    
    # Step 4: Generate multivariate normal data
    sigma <- solve(omega)  
    theta <- matrix(0, p, p)
    
    # Populate theta by checking non-zero off-diagonal elements of omega
    for (i in 1:p) {
      for (j in 1:p) {
        if (i != j && omega[i, j] != 0) {
          theta[i, j] <- 1
        }
      }
    }

    x = mvrnorm(n, mu = rep(0, p), Sigma = sigma)
    sigmahat = cov(x)
    
    if (vis == TRUE) {
        fullfig = par(mfrow = c(2, 2), pty = "s", omi = c(0.3, 
            0.3, 0.3, 0.3), mai = c(0.3, 0.3, 0.3, 0.3))
        fullfig[1] = image(theta, col = gray.colors(256), main = "Adjacency Matrix")
        fullfig[2] = image(sigma, col = gray.colors(256), main = "Covariance Matrix")
        g = graph.adjacency(theta, mode = "undirected", diag = FALSE)
        layout.grid = layout.fruchterman.reingold(g)
        fullfig[3] = plot(g, layout = layout.grid, edge.color = "gray50", 
            vertex.color = "red", vertex.size = 3, vertex.label = NA, 
            main = "Graph Pattern")
        fullfig[4] = image(sigmahat, col = gray.colors(256), 
            main = "Empirical Matrix")
        rm(fullfig, g, layout.grid)
        gc()
    }
    if (verbose) 
        cat("done.\n")
    rm(vis, verbose)
    gc()
    sim = list(data = x, sigma = sigma, sigmahat = sigmahat, 
        omega = omega, theta = Matrix(theta, sparse = FALSE), 
        act_sparsity = sum(theta)/(p * (p - 1)), graph.type = graph)
    class(sim) = "sim"
    return(sim)
}
```