1## Left Medcouple
2lmc <- function(x, na.rm = FALSE, ...) {
3    -mc(x[x <= median(x, na.rm = na.rm)], na.rm = na.rm, ...)
4}
5
6## Right Medcouple
7rmc <- function(x, na.rm = FALSE, ...) {
8    mc(x[x >= median(x, na.rm = na.rm)], na.rm = na.rm, ...)
9}
10
11
12## ## Generic function
13## mc <- function (x, ...)
14##       UseMethod("mc")
15
16## ## Default method (for numeric vectors):
17## mc.default <- function(x, na.rm = FALSE,
18mc <- function(x, na.rm = FALSE, doReflect = (length(x) <= 100)
19             , doScale = TRUE # <- chg default to 'FALSE' ?
20             , eps1 = 1e-14, eps2 = 1e-15 # << new in 0.93-2 (2018-07..)
21             , maxit = 100, trace.lev = 0
22             , full.result = FALSE
23               )
24{
25    x <- as.numeric(x)
26    ina <- is.na(x)
27    if (na.rm)
28        x <- x[!ina]
29    else if (any(ina))
30        return(NA_real_)
31    ## ==> x is NA-free from here on
32
33    ## if(length(l.. <- list(...)))
34    ##     stop("In mc(): invalid argument(s) : ",
35    ##          paste(sQuote(names(l..)), collapse=","), call. = FALSE)
36    rr <- mcComp(x, doReflect, doScale=doScale, eps1=eps1, eps2=eps2,
37                 maxit=maxit, trace.lev=trace.lev)
38
39    if(!(conv1 <- rr[["converged"]]) |
40       (doReflect && !(conv2 <- rr[["converged2"]]))) {
41        stop("mc(): not 'converged' ",
42             if(!conv1) paste("in", rr[["iter"]], "iterations"),
43             if(doReflect && !conv2)
44             paste(if(!conv1)" *and*",
45                   "'reflect part' in", rr[["iter2"]], "iterations"),
46             "; try enlarging eps1, eps2 !?\n")
47    }
48
49    m <- if (doReflect) (rr[["medc"]] - rr[["medc2"]]) / 2  else rr[["medc"]]
50    structure(m, mcComp = if(full.result) rr)
51}
52
53## eps1 = 1e-13, eps2 = eps1  <==>  original code which only had  'eps = 1e-13'
54##                                  hardcoded in C code.
55## These defaults do *not* make sense here, but in mc().
56## However, currently they are used in ../tests/mc-etc.R
57mcComp <- function(x, doReflect, doScale, eps1, eps2, maxit = 1000, trace.lev = 1)
58{
59    stopifnot(is.logical(doReflect), length(doReflect) == 1L, !is.na(doReflect),
60              is.logical(doScale),   length(doScale)   == 1L, !is.na(doScale),
61              is.1num(eps1), eps1 >= 0,
62              is.1num(eps2), eps2 >= 0,
63              length(maxit     <- as.integer(maxit)) == 1,
64              length(trace.lev <- as.integer(trace.lev)) == 1
65              )
66    ## Assumption [from caller, = mc()]: 'x' has no NAs (but can have +-Inf)
67    x <- as.numeric(x)
68    n <- as.integer(length(x))
69    eps <- as.double(c(eps1, eps2))
70    c.iter <- c(maxit, trace.lev)
71    ## NAOK=TRUE: to allow  +/- Inf to be passed
72    ans <- .C(mc_C, x, n,
73              eps = eps, iter = c.iter,
74              medc = double(1)
75            , doScale = doScale
76            , NAOK=TRUE)[c("medc", "eps", "iter")]
77    it <- ans[["iter"]]
78    ans[["converged"]] <- it[2] == 1
79    ans[["iter"]] <- it[1]
80
81    if (doReflect) { ## also compute on reflected data
82	a2 <- .C(mc_C, -x, n,
83                 eps2 = eps, iter2 = c.iter,
84                 medc2 = double(1)
85               , doScale = doScale
86               , NAOK=TRUE)[c("medc2", "iter2", "doScale")]
87        it <- a2[["iter2"]]
88        a2[["converged2"]] <- it[2] == 1
89        a2[["iter2"]] <- it[1]
90
91	c(ans, a2)
92    }
93    else ans
94}
95
96