RosettaCodeData/Task/P-value-correction/Wren/p-value-correction.wren

import "./dynamic" for Enum
import "./fmt" for Fmt
import "./seq" for Lst
import "./math" for Nums
import "./sort" for Sort

var Direction = Enum.create("Direction", ["UP", "DOWN"])

// test also for 'Unknown' correction type
var types = [
    "Benjamini-Hochberg", "Benjamini-Yekutieli", "Bonferroni", "Hochberg",
    "Holm", "Hommel", "Šidák", "Unknown"
]

var pFormat = Fn.new { |p, cols|
    var i = -cols
    var fmt = "$1.10f"
    return Lst.chunks(p, cols).map { |chunk|
        i = i + cols
        return Fmt.swrite("[$2d  $s", i, chunk.map { |v| Fmt.swrite(fmt, v) }.join(" "))
    }.join("\n")
}

var check = Fn.new { |p|
    if (p.count == 0 || Nums.min(p) < 0 || Nums.max(p) > 1) {
        Fiber.abort("p-values must be in range 0 to 1")
    }
    return p
}

var ratchet = Fn.new { |p, dir|
    var pp = p.toList
    var m = pp[0]
    if (dir == Direction.UP) {
        for (i in 1...pp.count) {
            if (pp[i] > m) pp[i] = m
            m = pp[i]
        }
    } else {
        for (i in 1...pp.count) {
            if (pp[i] < m) pp[i] = m
            m = pp[i]
        }
    }
    return pp.map { |v| (v < 1) ? v : 1 }.toList
}

var schwartzian = Fn.new { |p, mult, dir|
    var size = p.count
    var pwi = List.filled(size, null)
    for (i in 0...size) pwi[i] = [i, p[i]]
    var cmp = (dir == Direction.UP) ? Fn.new { |a, b| (b[1] - a[1]).sign } :
                                      Fn.new { |a, b| (a[1] - b[1]).sign }
    var order = Sort.merge(pwi, cmp).map { |e| e[0] }.toList
    var pa = List.filled(size, 0)
    for (i in 0...size) pa[i] = mult[i] * p[order[i]]
    pa = ratchet.call(pa, dir)
    var owi = List.filled(order.count, null)
    for (i in 0...order.count) owi[i] = [i, order[i]]
    cmp = Fn.new { |a, b| (a[1] - b[1]).sign }
    var order2 = Sort.merge(owi, cmp).map { |e| e[0] }.toList
    var res = List.filled(size, 0)
    for (i in 0...size) res[i] = pa[order2[i]]
    return res
}

var adjust = Fn.new { |p, type|
    var size = p.count
    if (size == 0) Fiber.abort("List cannot be empty.")
    if (type == "Benjamini-Hochberg") {
        var mult = List.filled(size, 0)
        for (i in 0...size) mult[i] = size / (size - i)
        return schwartzian.call(p, mult, Direction.UP)

    } else if (type == "Benjamini-Yekutieli") {
        var q = (1..size).reduce { |acc, i| acc + 1/i }
        var mult = List.filled(size, 0)
        for (i in 0...size) mult[i] = q * size / (size - i)
        return schwartzian.call(p, mult, Direction.UP)

    } else if (type == "Bonferroni") {
        return p.map { |v| (v * size).min(1) }.toList

    } else if (type == "Hochberg") {
        var mult = List.filled(size, 0)
        for (i in 0...size) mult[i] = i + 1
        return schwartzian.call(p, mult, Direction.UP)

    } else if (type == "Holm") {
        var mult = List.filled(size, 0)
        for (i in 0...size) mult[i] = size - i
        return schwartzian.call(p, mult, Direction.DOWN)

    } else if (type == "Hommel") {
        var pwi = List.filled(size, null)
        for (i in 0...size) pwi[i] = [i, p[i]]
        var cmp = Fn.new { |a, b| (a[1] - b[1]).sign }
        var order = Sort.merge(pwi, cmp).map { |e| e[0] }.toList
        var s = List.filled(size, 0)
        for (i in 0...size) s[i] = p[order[i]]
        var m = List.filled(size, 0)
        for (i in 0...size) m[i] = s[i] * size / (i + 1)
        var min = Nums.min(m)
        var q = List.filled(size, min)
        var pa = List.filled(size, min)
        for (j in size-1..2) {
            var lower = List.filled(size - j + 1, 0)                // lower indices
            for (i in 0...lower.count) lower[i] = i
            var upper = List.filled(j - 1, 0)                       // upper indices
            for (i in 0...upper.count) upper[i] = size - j + 1 + i
            var qmin = j * s[upper[0]] / 2
            for (i in 1...upper.count) {
                var temp = s[upper[i]] * j / (2 + i)
                if (temp < qmin) qmin = temp
            }
            for (i in 0...lower.count) {
                q[lower[i]] = qmin.min(s[lower[i]] * j)
            }
            for (i in 0...upper.count) q[upper[i]] = q[size - j]
            for (i in 0...size) if (pa[i] < q[i]) pa[i] = q[i]
        }
        var owi = List.filled(order.count, null)
        for (i in 0...order.count) owi[i] = [i, order[i]]
        var order2 = Sort.merge(owi, cmp).map { |e| e[0] }.toList
        var res = List.filled(size, 0)
        for (i in 0...size) res[i] = pa[order2[i]]
        return res

    } else if (type == "Šidák") {
        return p.map { |v| 1 - (1 - v).pow(size) }.toList

    } else {
        System.print("\nSorry, do not know how to do '%(type)' correction.\n" +
                     "Perhaps you want one of these?:\n" +
                     types[0...-1].map { |t| "  %(t)" }.join("\n")
        )
        Fiber.suspend()
    }
}

var adjusted = Fn.new { |p, type| "\n%(type)\n%(pFormat.call(adjust.call(check.call(p), type), 5))" }

var pValues = [
    4.533744e-01, 7.296024e-01, 9.936026e-02, 9.079658e-02, 1.801962e-01,
    8.752257e-01, 2.922222e-01, 9.115421e-01, 4.355806e-01, 5.324867e-01,
    4.926798e-01, 5.802978e-01, 3.485442e-01, 7.883130e-01, 2.729308e-01,
    8.502518e-01, 4.268138e-01, 6.442008e-01, 3.030266e-01, 5.001555e-02,
    3.194810e-01, 7.892933e-01, 9.991834e-01, 1.745691e-01, 9.037516e-01,
    1.198578e-01, 3.966083e-01, 1.403837e-02, 7.328671e-01, 6.793476e-02,
    4.040730e-03, 3.033349e-04, 1.125147e-02, 2.375072e-02, 5.818542e-04,
    3.075482e-04, 8.251272e-03, 1.356534e-03, 1.360696e-02, 3.764588e-04,
    1.801145e-05, 2.504456e-07, 3.310253e-02, 9.427839e-03, 8.791153e-04,
    2.177831e-04, 9.693054e-04, 6.610250e-05, 2.900813e-02, 5.735490e-03
]
types.each { |type| System.print(adjusted.call(pValues, type)) }