flowClust segfault on a given hardcoded prior

[mikejiang]

here is the reproducible data

library(flowClust)
load("flowClust_debug.rda")#see the attached file
tmix_results <- flowClust(fr, varNames = channels, K = 5, trans = 0, usePrior = "yes", prior = prior)
Using the serial version of flowClust

 *** caught segfault ***
address (nil), cause 'unknown'

Traceback:
 1: .C("flowClust", as.double(t(y)), as.integer(ly), as.integer(py),     as.integer(K[i]), w = rep(0, K[i]), mu = rep(0, K[i] * py),     precision = initprec, lambda = as.double(rep(lambda, length.out = (if (trans >         1) K[i] else 1))), nu = as.double(rep(nu, K[i])), z = rep(0,         ly * K[i]), u = rep(0, ly * K[i]), as.integer(if (M ==         0) label else maxLabel[[M]]), uncertainty = double(ly),     as.double(rep(u.cutoff, K[i])), as.double(z.cutoff), flagOutliers = integer(ly),     as.integer(B), as.double(tol), as.integer(trans), as.integer(nu.est),     logLike = as.double(0), as.integer(control$B.lambda), as.integer(control$B.brent),     as.double(control$tol.brent), as.double(control$xLow), as.double(control$xUp),     as.double(control$nuLow), as.double(control$nuUp), mu0 = as.double(t(Mu0)),     as.double(kappa0), nu0 = as.double(nu0), lambda0 = as.double(Lambda0),     omega0 = as.double(Omega0), w0 = as.double(w0), as.integer(.model),     oorder = as.integer(oorder), PACKAGE = "flowClust")
 2: doTryCatch(return(expr), name, parentenv, handler)
 3: tryCatchOne(expr, names, parentenv, handlers[[1L]])
 4: tryCatchList(expr, classes, parentenv, handlers)
 5: tryCatch(expr, error = function(e) {    call <- conditionCall(e)    if (!is.null(call)) {        if (identical(call[[1L]], quote(doTryCatch)))             call <- sys.call(-4L)        dcall <- deparse(call)[1L]        prefix <- paste("Error in", dcall, ": ")        LONG <- 75L        msg <- conditionMessage(e)        sm <- strsplit(msg, "\n")[[1L]]        w <- 14L + nchar(dcall, type = "w") + nchar(sm[1L], type = "w")        if (is.na(w))             w <- 14L + nchar(dcall, type = "b") + nchar(sm[1L],                 type = "b")        if (w > LONG)             prefix <- paste0(prefix, "\n  ")    }    else prefix <- "Error : "    msg <- paste0(prefix, conditionMessage(e), "\n")    .Internal(seterrmessage(msg[1L]))    if (!silent && identical(getOption("show.error.messages"),         TRUE)) {        cat(msg, file = outFile)        .Internal(printDeferredWarnings())    }    invisible(structure(msg, class = "try-error", condition = e))})
 6: try(.C("flowClust", as.double(t(y)), as.integer(ly), as.integer(py),     as.integer(K[i]), w = rep(0, K[i]), mu = rep(0, K[i] * py),     precision = initprec, lambda = as.double(rep(lambda, length.out = (if (trans >         1) K[i] else 1))), nu = as.double(rep(nu, K[i])), z = rep(0,         ly * K[i]), u = rep(0, ly * K[i]), as.integer(if (M ==         0) label else maxLabel[[M]]), uncertainty = double(ly),     as.double(rep(u.cutoff, K[i])), as.double(z.cutoff), flagOutliers = integer(ly),     as.integer(B), as.double(tol), as.integer(trans), as.integer(nu.est),     logLike = as.double(0), as.integer(control$B.lambda), as.integer(control$B.brent),     as.double(control$tol.brent), as.double(control$xLow), as.double(control$xUp),     as.double(control$nuLow), as.double(control$nuUp), mu0 = as.double(t(Mu0)),     as.double(kappa0), nu0 = as.double(nu0), lambda0 = as.double(Lambda0),     omega0 = as.double(Omega0), w0 = as.double(w0), as.integer(.model),     oorder = as.integer(oorder), PACKAGE = "flowClust"))
 7: FUN(X[[i]], ...)
 8: lapply(as.list(1:length(K)), .flowClustK, y, expName = expName,     varNames = varNames, K = K, criterion = criterion, nu = nu,     lambda = lambda, trans = trans, min.count = min.count, max.count = max.count,     min = min, max = max, randomStart = randomStart, include = include,     rm.max, rm.min, prior, usePrior, ...)
 9: flowClust(fr, varNames = channels, K = 5, trans = 0, usePrior = "yes",     prior = prior)

[SamGG]

Hi Mike! How did you build the prior? I didn't get any information in the help of the flowClust function.
BTW, the usePrior = "vague" leads to an error, "vague" being not recognized as a value, whereas the doc says it's ok.

[mikejiang]

@SamGG see ?openCyto::prior_flowClust

[mikejiang]

I am not able to pinpoint the exact line in flowClust C code because under gdb it runs through entire flowClust c function and sefault at the stage of dotcode.c from R. Likely memory is corrupted somewhere before R takes it over. Here is the output of valgrind, not sure if it is helpful at all

> tmix_results <- flowClust(fr, varNames = channels, K = 5, trans = 0, usePrior = "yes", prior = prior)
Using the serial version of flowClust
==7673== Invalid read of size 1
==7673==    at 0x4C2E0E2: strlen (in /usr/lib/valgrind/vgpreload_memcheck-amd64-linux.so)
==7673==    by 0x4F32FD8: Rf_mkChar (envir.c:3722)
==7673==    by 0x4F194B8: do_dotCode (dotcode.c:2504)
==7673==    by 0x4F54D00: Rf_eval (eval.c:728)
==7673==    by 0x4F55005: forcePromise (eval.c:520)
==7673==    by 0x4F53A5E: bcEval (eval.c:4709)
==7673==    by 0x4F546DF: Rf_eval (eval.c:624)
==7673==    by 0x4F55005: forcePromise (eval.c:520)
==7673==    by 0x4F53A5E: bcEval (eval.c:4709)
==7673==    by 0x4F546DF: Rf_eval (eval.c:624)
==7673==    by 0x4F55005: forcePromise (eval.c:520)
==7673==    by 0x4F53A5E: bcEval (eval.c:4709)
==7673==  Address 0x7ff8000000000000 is not stack'd, malloc'd or (recently) free'd
==7673== 

 *** caught segfault ***

[mikejiang]

Here is the given prior

> prior
$Mu0
           [,1]       [,2]
[1,] 0.10359586 0.03751837
[2,] 0.07607299 0.25967254
[3,] 0.03840841 0.28397826
[4,] 0.13033492 0.01071706
[5,] 0.01800955 0.12326509

$Lambda0
, , 1

            [,1]       [,2]
[1,] 413.4382049  0.1614980
[2,]  -0.6939923 40.5273378
[3,]   4.8931884  7.0903981
[4,]  -5.0695791  0.9772082
[5,]   0.9906558  0.4186971

, , 2

           [,1]       [,2]
[1,] -0.6939923  4.2322180
[2,] 13.1844080  7.0903981
[3,] -5.0695791 21.8155970
[4,] 70.8559750  0.4186971
[5,]  0.1614980 78.0131340


$Omega0
, , 1

             [,1]         [,2]
[1,] 4.442715e-07 0.000000e+00
[2,] 0.000000e+00 1.046215e-06
[3,] 7.673965e-08 0.000000e+00
[4,] 0.000000e+00 1.674699e-08
[5,] 3.078914e-08 0.000000e+00

, , 2

             [,1]         [,2]
[1,] 0.000000e+00 3.078914e-08
[2,] 4.442715e-07 0.000000e+00
[3,] 0.000000e+00 1.046215e-06
[4,] 7.673965e-08 0.000000e+00
[5,] 0.000000e+00 1.674699e-08


$nu0
[1]  5471.417 15281.385  8143.787  5255.161 22821.750

$w0
[1] 10942.83 30562.77 16287.57 10510.32 45643.50

> prior$nu0
[1]  5471.417 15281.385  8143.787  5255.161 22821.750

and what the data looks like

[SamGG]

The help page of prior_flowClust of my current version of openCyto (1.12.1) does not give me much information about the value returned.
I turned to the flowClust2Prior function. I ran flowClust on the data, then transform the result into prior using flowClust2Prior and finally ran again flowClust using those prior. No crash.
I copied the slots of your given prior to the result of flowClust2Prior and ran again flowClust. Crash.
I think the crash is due to a badly handled numerical configuration in the C code or a linked library. Not easy to debug. May be those prior are not an adequate start point to help flowClust in converging rapidly. I hope you will find out what is happening.

[mikejiang]

I haven't verified it, but I suspect the invalid input value of precision could be the cause. here it is expecting a double array, but the actual input contains the error string

Browse[2]> initprec
, , 1

     [,1]                 [,2]                 
[1,] "0.0491791239549737" "0.00806431364015211"
[2,] "0"                  "0.486006746491884"  

, , 2

     [,1]                                                                                                        
[1,] "Error in chol.default(solve(Lambda0[, , j])) : \n  the leading minor of order 1 is not positive definite\n"
[2,] "Error in chol.default(solve(Lambda0[, , j])) : \n  the leading minor of order 1 is not positive definite\n"
     [,2]                                                                                                        
[1,] "Error in chol.default(solve(Lambda0[, , j])) : \n  the leading minor of order 1 is not positive definite\n"
[2,] "Error in chol.default(solve(Lambda0[, , j])) : \n  the leading minor of order 1 is not positive definite\n"

which is captured here.
Apparently the Choleski Decomposition failed on one of Lambda0 value from given prior.
@gfinak may have more detailed explanation on why chol failed.

But as far as I concern, ideally there should be a validity check (type and length) on each one of 20 arguments of this .C call since these pointers are vulnerable to any kind of violations of assumption.
It may take some efforts initially, but will pay off in the long run by eliminating the pains of chasing segfault

[mikejiang]

@gfinak , looks like you were already planning on dealing with this. Maybe this is the time for it.