library(nnls)
library(beepr)
WR<-read.csv('/Users/davidcolby/Library/CloudStorage/OneDrive-Nexus365/DPhil/Corbetti/Data/CIPW/Boset_Glass.csv',sep=",",fill=TRUE,row.names = 1)
  A <- read.csv('/Users/davidcolby/Library/CloudStorage/OneDrive-Nexus365/DPhil/Corbetti/Data/CIPW/Boset_MB_Minerals.csv',sep=",",fill=TRUE,row.names=1)

FC_reverse <- function(WR, mins, Parent, Daughter, x, P) {
  pm.name <- Parent
  fm.name <- Daughter
  mins <- as.matrix(mins)
  mins[is.na(mins)] <- 0
  
  C0 <- WR[pm.name, colnames(mins)]
  CL <- WR[fm.name, colnames(mins)]
  C1 <- t(rbind(CL, mins))
  C0 <- t(C0)
  model <- nnls(C1, C0)
  r.squared <- sum(model$residuals^2)
  f <- model$x[1]
  
  cum.prop <- model$x[-1]/(1-f)
  
  if (P == 1){
      cat(round(100*(1-f), 3), "% fc", "\n")
      print(row.names(mins))
      print(round(cum.prop, 3))
      cat("Quality of the model (squared residuals):", round(r.squared, 4), "\n")
      #return(round(r.squared,3))
    }else{
      return(round(r.squared,3))
    
  }

}

library(future.apply)
plan(multisession)
library(furrr)
# Define your function to be executed in parallel
process_combinations <- function(j, A) {
  row_combinations <- combn(nrow(A), j)
  X <- c()
  min_residual <- Inf
  best_combination <- NULL
  for (i in seq_len(ncol(row_combinations))) {
    residual <- FC_reverse(WR, A[row_combinations[,i], ], Parent, Daughter, 1,2)
    X <- c(X, residual)
    if (residual < min_residual) {
      min_residual <- residual
      best_combination <- row_combinations[, i]
    }
  }
  return(list(X = X, min_residual = min_residual, best_combination = best_combination))
  
  
}

###### run from here
for (P in 30){#index of parent compositions
  for (D in 35){#index of daughter compositions
    X <- list()
    j_vals<-2:6
    cat('Parent:',P, ', Daughter:',D,', ')
    cat('Results: ')
    Parent <- P #Chose which rows is the parental melt
    Daughter<-D #Chose Daughter melt
    min_residual <- Inf
    best_combination <- NULL
    X <- future_map(j_vals, process_combinations, A = A)
    residuals <- sapply(X, function(x) x$X)
    min_residual <- min(sapply(X, function(x) x$min_residual))
    best_combination <- X[[which.min(sapply(X, function(x) x$min_residual))]]$best_combination
    
    best_A <- A[best_combination, ]
    FC_reverse(WR,best_A,Parent, Daughter,1,1)
    #beep(2)
    #beep() 
    
  }
}
beep(3)