convert it to tensorflow,  and python

###########
########### Build network and reproduce connections across datasets
########### Each connected subgraph is a dimension along which patients could be clustered
########### Check whether the same dimensions can be extracted via text encoder

require(CCA)

#### function trainN/test split

CVfun = function(DATA){
  train = sample(c(1:nrow(DATA)), nrow(DATA)*0.7)
  test  = setdiff(c(1:nrow(DATA)), train)
  return(list(train, test))
}

### function CCA
CCAfunc = function(X, Z, Xtest, Ztest, export="cor"){
  
  m1 = cancor(X, Z)
  m1cx = m1$xcoef[match(colnames(Xtest), rownames(m1$xcoef)),1]
  m1cy = m1$ycoef[match(colnames(Ztest), rownames(m1$ycoef)),1]
  
  m1cx[is.na(m1cx)] = 0
  m1cy[is.na(m1cy)] = 0
  
  xt = as.matrix(Xtest) %*% m1cx
  zt = as.matrix(Ztest) %*% m1cy
  
  tcor = cor(xt, zt)
  if(export == "cor"){
    return(tcor)
  }
  if(export == "vars"){
    return( prcomp( cbind(xt, zt), scale=TRUE)$x[,1] )
  }
  
}

d1 = read.csv("/zi/flstorage/HITKIP/datasets/IMAGEN/FU2/psytools/IMAGEN-IMGN_ADRS_CHILD_FU2-IMAGEN_DIGEST.csv")
d2 = read.csv("/zi/flstorage/HITKIP/datasets/IMAGEN/FU2/psytools/IMAGEN-IMGN_ANXDX_CHILD_FU2-IMAGEN_DIGEST.csv")
d3 = read.csv("/zi/flstorage/HITKIP/datasets/IMAGEN/FU2/psytools/IMAGEN-IMGN_AUDIT_CHILD_FU2-IMAGEN_DIGEST.csv")
d4 = read.csv("/zi/flstorage/HITKIP/datasets/IMAGEN/FU2/psytools/IMAGEN-IMGN_BIS_CHILD_FU2-IMAGEN_DIGEST.csv")
d5 = read.csv("/zi/flstorage/HITKIP/datasets/IMAGEN/FU2/psytools/IMAGEN-IMGN_CAPE_CHILD_FU2-IMAGEN_DIGEST.csv")
d6 = read.csv("/zi/flstorage/HITKIP/datasets/IMAGEN/FU2/psytools/IMAGEN-IMGN_CSI_CHILD_FU2-IMAGEN_DIGEST.csv")
d7 = read.csv("/zi/flstorage/HITKIP/datasets/IMAGEN/FU2/psytools/IMAGEN-IMGN_CTQ_CHILD_FU2-IMAGEN_DIGEST.csv")
d8 = read.csv("/zi/flstorage/HITKIP/datasets/IMAGEN/FU2/psytools/IMAGEN-IMGN_EDEQ_CHILD_FU2-IMAGEN_DIGEST.csv")
d9 = read.csv("/zi/flstorage/HITKIP/datasets/IMAGEN/FU2/psytools/IMAGEN-IMGN_ESPAD_CHILD_FU2-IMAGEN_DIGEST.csv")
d10= read.csv("/zi/flstorage/HITKIP/datasets/IMAGEN/FU2/psytools/IMAGEN-IMGN_IRI_CHILD_FU2-IMAGEN_DIGEST.csv")
d11= read.csv("/zi/flstorage/HITKIP/datasets/IMAGEN/FU2/psytools/IMAGEN-IMGN_JVQ_CHILD_FU2-IMAGEN_DIGEST.csv")
d12= read.csv("/zi/flstorage/HITKIP/datasets/IMAGEN/FU2/psytools/IMAGEN-IMGN_PAAQ_CHILD_FU2-IMAGEN_DIGEST.csv")
d13= read.csv("/zi/flstorage/HITKIP/datasets/IMAGEN/FU2/psytools/IMAGEN-IMGN_RRS_CHILD_FU2-IMAGEN_DIGEST.csv")
d14= read.csv("/zi/flstorage/HITKIP/datasets/IMAGEN/FU2/psytools/IMAGEN-IMGN_TCI_CHILD_FU2-IMAGEN_DIGEST.csv")
d15= read.csv("/zi/flstorage/HITKIP/datasets/IMAGEN/FU2/psytools/IMAGEN-IMGN_TFEQ_CHILD_FU2-IMAGEN_DIGEST.csv")

                
                
### identify overlapping subjects

d7$User.code = gsub("-I", "-C", d7$User.code)

ucode.comb = table(c(d1$User.code, d2$User.code, d3$User.code, d4$User.code, d5$User.code,
                     d6$User.code, d7$User.code, d8$User.code, d9$User.code, d10$User.code,
                     d11$User.code,d12$User.code,d13$User.code,d14$User.code,d15$User.code))
ucode.comb = names(ucode.comb[which(ucode.comb == 15)])
length(ucode.comb)

d1Mat = d1[match(ucode.comb, d1$User.code),grep("adrs",colnames(d1))]
d2Mat = d2[match(ucode.comb, d2$User.code),grep("ANXDX",colnames(d2))]
d3Mat = d3[match(ucode.comb, d3$User.code),grep("audit",colnames(d3))]
d4Mat = d4[match(ucode.comb, d4$User.code),grep("BIS",colnames(d4))]
d5Mat = d5[match(ucode.comb, d5$User.code),grep("CAPE42",colnames(d5))]
d6Mat = d6[match(ucode.comb, d6$User.code),grep("item",colnames(d6))]
d7Mat = d7[match(ucode.comb, d7$User.code),grep("CTQ",colnames(d7))]
d8Mat = d8[match(ucode.comb, d8$User.code),grep("EDEQ",colnames(d8))]
d9Mat = d9[match(ucode.comb, d9$User.code),grep("SCID",colnames(d9))]
d10Mat= d10[match(ucode.comb, d10$User.code),grep("IRI",colnames(d10))]
d11Mat= d11[match(ucode.comb, d11$User.code),grep("JVQ",colnames(d11))]
d12Mat= d12[match(ucode.comb, d12$User.code),grep("PAAQ",colnames(d12))]
d13Mat= d13[match(ucode.comb, d13$User.code),grep("RRS",colnames(d13))]
d14Mat= d14[match(ucode.comb, d14$User.code),grep("tci",colnames(d14))]
d15Mat= d15[match(ucode.comb, d15$User.code),grep("TFEQ",colnames(d15))]

d12Mat = apply(d12Mat,2,function(x){as.numeric(as.factor(x))})

#d1Mat = d1Mat[,-ncol(d1Mat)]

dList = list(d1Mat, d2Mat, d3Mat, d4Mat, d5Mat, d6Mat, d7Mat, d8Mat, d9Mat, d10Mat,
             d11Mat, d12Mat, d13Mat, d14Mat, d15Mat)

dList = lapply(dList, function(x){
          x = apply(x,2,function(y){as.numeric(as.vector(y))})
          x[is.na(x)] = 0
          x[x==""] = 0 
          remzero = which( apply(x,2,function(y){ sum(as.numeric(y),na.rm=TRUE)}) == 0)
          if(length(remzero) > 0) x = x[,-remzero]
          return(x)
})



#########
######### load metadata, match and correct for covariates

meta = read.csv("/zi/flstorage/HITKIP/datasets/IMAGEN/FU3/participants/IMAGEN_demographics.csv")

ucode.combMOD = ucode.comb
ucode.combMOD = as.numeric(gsub("-C", "", ucode.combMOD))

meta = meta[match(ucode.combMOD, meta[,1]),]

dList = lapply(dList, function(x){
  x = apply(x,2,function(y){ lm(y ~ as.factor(meta[,2]) + as.factor(meta[,3]) )$res })
  return(x)
})



#################################
#################################


## unique combinations of ratings scales

rtcmb = combn(c(1:15), 2)


## run repeated train test splits
## to evaluate median correlations of canonical variates between rating scales in test data 

repCOR = numeric()
for(repi in 1:20){

    ## 70/30 train test split
    traintest = CVfun(dList[[1]])
    train = traintest[[1]]
    test  = traintest[[2]]
    
    ## run for each pair
    
    corvec = numeric()
    for(i in 1:ncol(rtcmb)){
      curpair = rtcmb[,i]
      ip = curpair[1]
      is = curpair[2]
      curcca = CCAfunc(dList[[ip]][train,], dList[[is]][train,], 
                       dList[[ip]][test,], dList[[is]][test,], export = "cor")  
      corvec = c(corvec, curcca)
    }

    repCOR = rbind(repCOR, corvec)

} ## end of repeated train/ test split


finalCor = apply(repCOR, 2, median)
finalCor

##############
############## convert to matrix and cluster

corMat = matrix(0, 15, 15)
corMat[lower.tri(corMat)] = finalCor
corMat = corMat+ t(corMat)

corMat[corMat < 0.5] = 0  ## link thresholding
corMat [corMat != 0] = 1

library(igraph)

corgr = graph_from_adjacency_matrix(corMat, mode="undirected")

clsol = clusters(corgr)

## compress clusters into first PC

PCsol = numeric()
for(i in 1:clsol$no){
  ### identify non-zero elements in adjacency matrix
  curclIND = unique(clsol$membership)[i]
  curCor = corMat
  curCor[-which(clsol$membership == curclIND), -which(clsol$membership == curclIND)] = 0
  
  if( max(curCor) == 0) next ### no canonical correlations abova threshold
  
  curcmb = t(which(curCor != 0, arr.ind = TRUE))
  
  ### calculate canonical variates via bootstrap on for the entire IMAGEN sample
  
  curcmbRES = numeric() ### store variates for each cluster
  for(ii in 1:ncol(curcmb)){
    curpair = curcmb[,ii]
    ip = curpair[1]
    is = curpair[2]
    
    ### repeated traintest split
    
    ccaPCcur = numeric()
    for(repi in 1:100){
        ## 70/30 train test split
        reptraintest = CVfun(dList[[1]])
        reptrain = reptraintest[[1]]
        reptest  = reptraintest[[2]]
    
        curcca = CCAfunc(dList[[ip]][reptrain,], dList[[is]][reptrain,], 
                     dList[[ip]][reptest,], dList[[is]][reptest,], export = "vars")   ### extract first PC of both canonical variates
        curcca = curcca[match(seq(1:nrow(dList[[1]])), reptest)] ## match PC back to original identifiers
        
        if(repi > 1){   ## align PCs to first PC directionaly
          if(cor(curcca, ccaPCcur[,1], use="pairwise.complete.obs") < 0) curcca = curcca * -1
        }
        
        ccaPCcur = cbind(ccaPCcur, curcca)
    }    
    #cor(ccaPCcur, use="pairwise.complete.obs")       ## check that correlation structure OK
        
    ### determine median variate
    
    ccaPCcurMed = apply(ccaPCcur,1,function(x){median(x,na.rm=TRUE)})
    
    if(ii > 1){   ## align median variate to first median variate directionaly
      if(cor(ccaPCcurMed, curcmbRES[,1], use="pairwise.complete.obs") < 0) ccaPCcurMed = ccaPCcurMed * -1
    }
    
    curcmbRES = cbind(curcmbRES, ccaPCcurMed)
  }
  #cor(curcmbRES, use="pairwise.complete.obs")        ## check that correlation structure OK
  
  ### determine first PC for the canonical variates
  
  varPCA = prcomp(curcmbRES, scale=TRUE)$x[,1]
  PCsol = cbind(PCsol, varPCA)
}

dim(PCsol)
head(PCsol)


####################
####################

#################### Explore whether autoencoder captures the canonical variate dimensions

### combine data into dataframe


dListMat = do.call(cbind, dList)

library(h2o) 

h2o.no_progress()  # turn off progress bars
h2o.init(max_mem_size = "5g")  # initialize H2O instance

features <- as.h2o(dListMat)


ae1 <- h2o.deeplearning(
  x = seq_along(features),
  training_frame = features,
  autoencoder = TRUE,
  hidden = c(50, 10, 50),
  activation = 'Tanh',
  sparse = FALSE
)


# Extract the deep features
ae1_codings <- h2o.deepfeatures(ae1, features, layer = 2)

head(ae1_codings)

round(cor(PCsol, as.matrix(ae1_codings)) ,2)