library(Rcpp)
sourceCpp(file="/storage11_7T/fuy/TADA-A/cell_WES/DNM/ll_sum.cpp")
sourceCpp(file="/storage11_7T/fuy/TADA-A/cell_WES/DNM/multi_annota.cpp")
sourceCpp(file="/storage11_7T/fuy/TADA-A/cell_WES/DNM/post.cpp")
sourceCpp(file="/storage11_7T/fuy/TADA-A/cell_WES/DNM/ll_sum2.cpp")

compact_data = readRDS("/storage11_7T/fuy/TADA-A/cell_WES/DNM/2021-01-19cp6_all_gene_PO_6788SNV_sf_uni_prior_compact.rds")
data=compact_data$base_info

gene_prior_file = "/storage11_7T/fuy/TADA-A/db/MS_data/prior/new_uniform_gene_prior.txt"
optimization_iteration = 2000

library(data.table)
gene_prior = fread(gene_prior_file)
gene_prior = gene_prior[order(gene_prior$genename),]

selected_annotations = c(2)
logP_Zg0 = sumall0(data)

fr_pi <- function(par){

    all_rr = par[1:length(selected_annotations)] 
    gene_prior$prior <- rep(par[length(par)],nrow(gene_prior))
    
    logP_Zg1 = sumall1(data,selected_annotations,all_rr)

    logP_table<-data.table(logP_Zg1 = logP_Zg1, logP_Zg0 = logP_Zg0, genename = names(data))
    logP_table <- logP_table[gene_prior, on = "genename"]
    logP_table <- logP_table[complete.cases(logP_table)]
    idx = match(unique(logP_table$genename),logP_table$genename)
    idx2 = c(idx,nrow(logP_table) + 1)
    pr = logP_table[idx,]$prior
    ll_sum1 <- ll_sum(idx2,logP_table$logP_Zg1,logP_table$logP_Zg0,pr)
    ll_sum1 
    }

ll.lst = c()
gm.lst = c()

for(i in seq(0.01,0.2,0.02)){
    
a = c(rep(.1,1),i)

tmm = proc.time()
gama = optim(a, fr_pi ,control=list("fnscale"=-1, "maxit" = optimization_iteration))$par
gm.lst = rbind(gm.lst,gama)
ll = fr_pi(gama)
ll.lst = c(ll.lst,ll)
proc.time() - tmm

gm.lst
ll.lst
    
    }

plot(ll.lst)

which.max(ll.lst)

gm.lst
ll.lst