library(stats)
library(tidyr)
library(hier.part)
library(dplyr)

K<- 100  #value transaction
gamma1 <- 0.90 ##proba basic
gamma <- c()
gamma[1] <- gamma1 
alpha <- 0.50 ## depletude factor
n <- 3 ##number of transaction


##only work for t=3
for (t in c(2:n)){
  gamma[t]= gamma[t-1]*alpha
  
}

sol1 <- -3*(1-gamma[1])*K
sol2 <- (-K*(1-gamma[1])*gamma[1]+(-2*K+K*gamma[2])*(1-gamma[1]))*gamma[1]+((-2*K+K*gamma[2])*gamma[2]+(-3*K+K*gamma[3])*(1-gamma[2]))*(1-gamma[1])




##general case

###function to order my binary matrix
to.value <- function(vect){
  L <- length(vect)
  value <- 0
  for (i in L:1){
    value <- value + 2^(L-i)*vect[i]
  }
  return(value)
}

# matrix is a square binary matrix
# returns the values of the columns in a vector

column.values <- function(matrix){
  result <- numeric(ncol(matrix))
  for (i in 1:length(result)){
    result[i] <- to.value(matrix[,i])
  }
  return(result)
}

# matrix is a square binary matrix
# returns the matrix in the prescribed order

get.ordered.matrix <- function(matrix){
  vals <- column.values(matrix)
  return(matrix[,rev(order(vals))])
}


n <-3
K <- 100
gamma1 <- 0.99 ##proba bas
alpha <- 0.20 ## depletude factor

gamma <- data.frame(proba_value=matrix( 1,n))
gamma[1,] <- gamma1
for (t in c(2:n)){
  gamma[t,]= gamma[t-1,]*alpha
  
}
gamma$gamma <- as.character(c(1:n))

######prepa last even
scenarios <-  as.matrix(expand.grid(rep(list(1:0), n)))
scenarios <- as.data.frame(get.ordered.matrix(scenarios))
names(scenarios) <- paste("event",c(1:n),sep="_")
names(scenarios)[length(names(scenarios))]<-"last_event" 
scenarios_prev <- as.data.frame(scenarios[,c(1:n-1)])
scenarios$counting_0_prev <- as.character(pmax(rowSums(scenarios_prev==0)+1,1))##counting how many fail before

scenarios <- merge(scenarios,gamma,by.x = "counting_0_prev", by.y = "gamma")###associate rv

scenarios$proba_scen <- scenarios$last_event*scenarios$proba_value + (1-scenarios$last_event)*(1-scenarios$proba_value) #gamma for success and 1-gamma for fail
scenarios$rv <- -(as.numeric(scenarios$counting_0_prev)-1)*K- (1-scenarios$last_event)*K   ###value of loss at the end
scenarios$lost_proba <- scenarios$rv*scenarios$proba_scen


#####now we loop to arrive at the end$

for( j in c(2:n-1)){
  print(j)
  #j<-2 #for test purpose
  k <- n-j

namevar <- names(scenarios)[!names(scenarios) %in% c("counting_0_prev","last_event","proba_value","proba_scen","rv")]
namevartogroup <- names(scenarios)[!names(scenarios) %in% c("counting_0_prev","last_event","proba_value","proba_scen","rv","lost_proba")]
new_scenarios <- scenarios[,namevar]%>%
  group_by_at(namevartogroup)%>%
  summarize_each(fun=sum)



# scenarios <-  as.matrix(expand.grid(rep(list(1:0), k)))
# scenarios <- as.data.frame(get.ordered.matrix(scenarios))
# names(scenarios) <- paste("event",c(1:k),sep="_")
scenarios <- new_scenarios 
names(scenarios)[length(names(scenarios))-1]<-"last_event" 
scenarios_prev <- as.data.frame(scenarios[,c(1:k-1)])
scenarios$counting_0_prev <- as.character(pmax(rowSums(scenarios_prev==0)+1,1))##counting how many fail before
scenarios <- merge(scenarios,gamma,by.x = "counting_0_prev", by.y = "gamma")
scenarios$proba_scen <- scenarios$last_event*scenarios$proba_value + (1-scenarios$last_event)*(1-scenarios$proba_value)
scenarios$rv <- scenarios$lost_proba  ###value of loss at the end
scenarios$lost_proba <- scenarios$rv*scenarios$proba_scen

}

##finalization
newsol <- scenarios[scenarios$last_event==0,"lost_proba"]+scenarios[scenarios$last_event==1,"lost_proba"]

sol1 <- -n*(1-gamma1)*K
explisol3 <- (-K*(1-gamma1)*gamma1+(-2*K+K*gamma[2,"proba_value"])*(1-gamma1))*gamma1+((-2*K+K*gamma[2,"proba_value"])*gamma[2,"proba_value"]+(-3*K+K*gamma[3,"proba_value"])*(1-gamma[2,"proba_value"]))*(1-gamma1)



###############################old




scenarios <- combos(t)$binary
x <- 3
mat <- as.matrix(expand.grid(rep(list(1:0), x)))
mat[order(mat),]


ix<-order(mat[,1],mat[,2],mat[,3])
ix
mat[ix,]