solving the old dependency group problem

RBBGCMuso/R/assistantFunctions.R

@@ -150,6 +150,7 @@ readValuesFromFile  <- function(filename, linums){
 
 	return(values)
 }
+
 #' readMeasuredMuso
 #'
 #' MuSo data reader

RBBGCMuso/R/musoRand.R

@@ -8,7 +8,8 @@
 #' @importFrom limSolve xsample
 #' @export
 
-musoRand <- function(parameters, iterations=3000, fileType="epc", constrains = NULL, burnin = NULL){
+musoRand <- function(parameters, iterations=3000, fileType="epc", sourceFile=NULL, constrains = NULL, burnin = NULL){
+
     if(is.null(constrains)){
         constMatrix <- constrains
         constMatrix <- getOption("RMuso_constMatrix")[[fileType]][[as.character(getOption("RMuso_version"))]]
@@ -20,20 +21,48 @@ musoRand <- function(parameters, iterations=3000, fileType="epc", constrains = N
     constMatrix <- constMatrix[,-1]
 
     depTableMaker <- function(constMatrix,parameters){
-        # browser()
         parameters <- parameters[order(parameters[,1]),] ## BUG!!!
         selectedRows <- constMatrix[,"INDEX"] %in% parameters[,1]
+        # constMatrix[constMatrix[,"INDEX"] %in% parameters[,1],]
         rankList <- rank(constMatrix[selectedRows,2])
         constMatrix[selectedRows,c(5,6)] <- parameters[rankList,c(2,3)]
+
+
+
+
+
         logiConstrain <- (constMatrix[,"GROUP"] %in% constMatrix[constMatrix[,"INDEX"] %in% parameters[,1],"GROUP"] &
                           (constMatrix[,"GROUP"]!=0)) | ((constMatrix[,"INDEX"] %in% parameters[,1]) & (constMatrix[,"GROUP"] == 0))
         constMatrix <- constMatrix[logiConstrain,]
         constMatrix <- constMatrix[order(apply(constMatrix[,7:8],1,function(x){x[1]/10+abs(x[2])})),]
-	constMatrix
+
+
+
+        paragroups <- unique(constMatrix[constMatrix[,"GROUP"] != 0,"GROUP"])
+        missingMembers <- list()
+        for(group in paragroups){
+            groupMembers <- constMatrix[constMatrix[,"GROUP"] == group,"INDEX"]
+            missingMemberElems <- groupMembers[is.na(match(groupMembers,parameters[,1]))]
+            if(length(missingMemberElems) > 0){
+                missingMembers$indices  <- c(missingMembers$indices,
+                                             match(missingMemberElems,constMatrix[,"INDEX"]))
+                if(is.null(sourceFile)){
+                    stop(sprintf("All group members of the group (%s) have to be in parameters if sourceFile not(epc file)",group))
+                }
+
+                missingMembers$values <- c(missingMembers$values,
+                                           suppressWarnings(readValuesFromFile(sourceFile,missingMemberElems))) 
+            }
+        }
+
+        list(depTable=constMatrix,missingMembers=missingMembers)
     }
     # browser()
-    genMat0 <- function(dep){
+    genMat0 <- function(dep, missingIndices){
         numberOfVariable <- nrow(dep)
+        if(length(missingMembers) != 0){
+            dep <- dep[-missingIndices,]
+        }
         G <- rbind(diag(numberOfVariable), -1*diag(numberOfVariable))
         h <- c(dependences[,5], -1*dependences[,6])
         return(list(G=G,h=h))
@@ -70,7 +99,11 @@ musoRand <- function(parameters, iterations=3000, fileType="epc", constrains = N
             if(Range[2]==N){
                 G<-cbind(matrix(ncol=(Range[1]-1),nrow=numRowsInG,data=0),G)
             } else {
-		G <- cbind(matrix(ncol=(Range[1]-1),nrow=numRowsInG,data=0),G,matrix(ncol=(N-Range[2]),nrow=numRowsInG,data=0))
+                G <- cbind(matrix(ncol=(Range[1]-1),
+                                  nrow=numRowsInG,data=0),
+                           G,
+                           matrix(ncol=(N-Range[2]),
+                                  nrow=numRowsInG,data=0))
             }
         }
         return(list(G=-1*G,h=-1*rep(0,nrow(G))))
@@ -149,6 +182,8 @@ musoRand <- function(parameters, iterations=3000, fileType="epc", constrains = N
     }
 
     dependences <- depTableMaker(constMatrix, parameters)
+    missingMembers <- dependences$missingMembers
+    dependences <- dependences$depTable
     dependences <- cbind(dependences,1:nrow(dependences))
     colnames(dependences)[ncol(dependences)] <- "rowIndex"
     # browser()
@@ -168,7 +203,7 @@ musoRand <- function(parameters, iterations=3000, fileType="epc", constrains = N
             }
         }
 
-        Gh0<- genMat0(dependences)
+        Gh0<- genMat0(dependences, missingMembers$indices)
         G <- do.call(rbind,lapply(Gh,function(x){x$G}))
         G<- rbind(Gh0$G,G)
         h <- do.call(c,lapply(Gh,function(x){x$h}))
@@ -176,12 +211,22 @@ musoRand <- function(parameters, iterations=3000, fileType="epc", constrains = N
         E <- do.call(rbind,lapply(Ef,function(x){x$E}))
         f <- do.call(c,lapply(Ef,function(x){x$f}))
         # browser()
-        randVal <- suppressWarnings(limSolve::xsample(G=G,H=h,E=E,F=f,burninlength=burnin, iter = iterations))$X
+        if(length(missingMembers$indices)!=0){
-    } else{
+            Ep <- matrix(data=0,ncol=numberOfVariable,nrow=length(missingMembers$indices))
-        Gh0<-genMat0(dependences)
+            Ep[1:length(missingMembers$indices),missingMembers$indices] <- 1
-        randVal <- suppressWarnings(xsample(G=Gh0$G,H=Gh0$h, iter = iterations))$X
+            E <- rbind(E,Ep)
+            f <- c(f,missingMembers$values)
         }
 
+        randVal <- suppressWarnings(limSolve::xsample(G=G,H=h,E=E,F=f,burninlength=burnin, iter = iterations))$X
+    } else{
+        Gh0<-genMat0(dependences,NULL)
+        randVal <- suppressWarnings(xsample(G=Gh0$G,H=Gh0$h, iter = iterations))$X
+    }
+    if(length(missingMembers$indices)!=0){
+        return(list(INDEX = dependences$INDEX[-missingMembers$indices],
+                    randVal=randVal[,-missingMembers$indices]))
+    }
     results <- list(INDEX =dependences$INDEX, randVal=randVal)
     return(results)
 }