text:
#Trimmed mean; amount of trimming defaults to 20% or 0.2
function tmean{S <: Real}(x::Vector{S}; tr::Real=0.2)
    if tr < 0 || tr > 0.5
        error("tr cannot be smaller than 0 or larger than 0.5")
    elseif tr == 0
        return mean(x)
    elseif tr == .5
        return median(x, checknan=false)
    else
        n   = length(x)
        lo  = floor(n*tr)+1
        hi  = n+1-lo
        return mean(sort!(copy(x))[lo:hi])
    end
end
#tmean{S <: Real}(x::DataVector{S}; tr::Real=0.2)=tmean(removeNA(x), tr=tr)

function tmean!{S <: Real}(x::Vector{S}; tr::Real=0.2)
    if tr < 0 || tr > 0.5
        error("tr cannot be smaller than 0 or larger than 0.5")
    elseif tr == 0
        return mean(x)
    elseif tr == .5
        return median!(x, checknan=false)
    else
        n   = length(x)
        lo  = floor(n*tr)+1
        hi  = n+1-lo
        return mean(sort!(x)[lo:hi])
    end
end

#Winsorize data
function winval{S <: Real}(x::Vector{S}; tr::Real=0.2)
    xcopy   = copy(x)
    n       = length(x)
    xcopy   = sort!(xcopy)
    ibot    = floor(tr*n)+1
    itop    = n-ibot+1
    xbot, xtop = xcopy[ibot], xcopy[itop]
    return  [x[i]<=xbot?xbot:(x[i]>=xtop?xtop:x[i]) for i=1:n]
end
winval{S <: Real}(x::DataVector{S}; tr::Real=0.2)=winval(removeNA(x), tr=tr)

#Winsorized mean
winmean{S <: Real}(x::Vector{S}; tr=0.2)=mean(winval(x, tr=tr))
#winmean{S <: Real}(x::DataVector{S}; tr::Real=0.2)=mean(winval(removeNA(x), tr=tr))

#Winsorized variance
winvar{S <: Real}(x::Vector{S}; tr=0.2)=var(winval(x, tr=tr))
#winvar{S <: Real}(x::DataVector{S}; tr::Real=0.2)=var(winval(remove(x), tr=tr))

#winsorized standard deviation
winstd{S <: Real}(x::Vector{S}; tr=0.2)=std(winval(x, tr=tr))
#winstd{S <: Real}(x::DataVector{S}; tr::Real=0.2)=std(winval(remove(x), tr=tr))

#Estimate the standard error of the gamma trimmed mean
trimse{S <: Real}(x::Vector{S}; tr::Real=0.2)=
    sqrt(winvar(x,tr=tr))/((1-2tr)*sqrt(length(x)))
#trimse{S <: Real}(x::DataVector{S}; tr::Real=0.2)=trimse(removeNA(x), tr=tr)

#Compute a 1-alpha confidence interval for the trimmed mean
function trimci{S <: Real}(x::Vector{S}; tr::Real=0.2, alpha::Real=0.05, nullvalue::Real=0, method=true)
    se  = trimse(x, tr=tr)
    n   = length(x)
    df::Int64   = n-2*floor(tr*n)-1
    estimate    = tmean(x, tr=tr)
    confint     = [estimate-Rmath.qt(1-alpha/2, df)*se,
                   estimate+Rmath.qt(1-alpha/2, df)*se]
    statistic   = (estimate-nullvalue)/se
    pval        = 2*(1-Rmath.pt(abs(statistic),df))
    METHOD      = method ? "1-alpha confidence interval for the trimmed mean\n": nothing
    output           = testOutput()
    output.method    = METHOD
    output.df        = df
    output.estimate  = estimate
    output.ci        = confint
    output.statistic = statistic
    output.p         = pval
    return output
end
#trimci{S <: Real}(x::DataVector{S}; tr=0.2, alpha=0.05, nullvalue=0, method=true)=
 #   trimci(removeNA(x), tr=tr, alpha=alpha, nullvalue=nullvalue, method=method)

#Stein's method
function stein1{S <: Real}(x::Vector{S}, delta::Real; alpha::Real=0.05, pow::Real=0.8, oneside::Bool=false, n=nothing, variance=nothing)
    delta    = abs(delta)
    n        = n==nothing?length(x):n
    variance = variance==nothing?var(x):variance
    df       = n-1
    alpha    = !oneside?alpha/2:alpha
    d        = (delta/(Rmath.qt(pow, df)-Rmath.qt(alpha, df)))*(delta/(Rmath.qt(pow, df)-Rmath.qt(alpha, df)))
    N        = max([n, floor(variance/d)+1])
    return integer(N)
end
#stein1{S <: Real}(x::DataVector{S}, delta::Real; alpha::Real=0.05, pow::Real=0.8, oneside::Bool=false, n=nothing, variance=nothing)=
 #   stein1(removeNA(x), delta, alpha=alpha, pow=pow, oneside=oneside, n=n, variance=variance)

function stein2{S <: Real, T <: Real}(x1::Vector{S}, x2::Vector{T}; mu0::Real=0.0, alpha::Real=0.05, method::Bool=true)
    n       = length(x1)
    df      = n-1
    N       = n+length(x2)
    xbar    = mean([x1, x2])
    test    = sqrt(N)*(xbar-mu0)/std(x1)
    crit    = Rmath.qt(1-alpha/2, df)
    lo, hi  = xbar-crit*std(x1), xbar+crit*std(x1)
    sig     = 2*(1-Rmath.pt(test, df))
    METHOD  = method?"The 2nd stage of Stein's method":nothing
    output  =testOutput()
    output.method    = METHOD
    output.df        = df
    output.estimate  = xbar
    output.ci        = [lo, hi]
    output.statistic = test
    output.crit      = crit
    output.p         = sig
    return output
end
#stein2{S <: Real, T <: Real}(x1::DataVector{S}, x2::DataVector{T}; mu0::Real=0.0, alpha::Real=0.05, method::Bool=true)=
 #   stein2(removeNA(x1), removeNA(x2), mu0=mu0, alpha=alpha, method=method)

#Ideal fourths
function idealf{S <: Real}(x::Vector{S}; method::Bool=true)
    n       = length(x)
    j       = integer(floor(n/4+5/12))
    y       = sort(x)
    g       = n/4-j+5/12
    k       = n-j+1
    METHOD  = method?"The ideal fourths":nothing
    idealfOutput((1-g).*y[j]+g.*y[j+1], (1-g).*y[k]+g.*y[k-1], METHOD)
end
#idealf{S <: Real}(x::DataVector{S}; method::Bool=true)=idealf(removeNA(x), method=method)

#Percentage bend midvariance
function pbvar{S <: Real}(x::Vector{S}; beta::Real=0.2)
    output  = 0.0
    const n = length(x)
    w       = zeros( n)
    med     = median(x, checknan=false)
    w       = [abs(x[i]-med) for i=1:n]
    w       = sort!(w)
    m       = floor((1-beta).*n+0.5)
    z       = 0.0
    omega   = w[m]
    counter = 0
    if omega > 0
        for i = 1:n
            y = (x[i]-med)/omega
            z += y>1.0?1.0:(y< -1.0?1.0:y*y)
            counter+=abs(y) < 1?1:0
        end
        output = n*(omega*omega)*z/(counter*counter)
    end
    return output
end
#pbvar{S <: Real}(x::DataVector{S}; beta::Real=0.2)=pbvar(removeNA(x), beta=beta)

#Biweight midvariance
function bivar{S <: Real}(x::Vector{S})
    m   = median(x, checknan=false)
    n   = length(x)
    u   = 0.0
    MAD = mad(x)
    av  = zeros( n)
    top = bot=0.0
    q   = Rmath.qnorm(0.75)
    for i = 1:n
        u     = abs(x[i]-m)./(9.*q.*MAD)
        av[i] = u<1.0?1:0
        top   += (n*av[i]*(x[i]-m)*(x[i]-m)*(1-u*u).^4)
        bot   += av[i]*(1-u*u)*(1-5*u*u)
    end
    top/(bot*bot)
end
#bivar{S <: Real}(x::DataVector{S})=bivar(removeNA(x))

#Compute the tau measure of location as described in Yohai and Zamar (JASA, 1988, 83, 406-413)
function tauloc{S <: Real}(x::Vector{S}; cval::Real=4.5)
    s     = Rmath.qnorm(0.75)*mad(x)
    med   = median(x, checknan=false)
    n     = length(x)
    Wnom  = Wden = 0.0
    for i in 1:n
        y = (x[i]-med)/s
        temp = (1.0-y*y/cval/cval)*(1.0-y*y/cval/cval)
        if abs(temp) <= cval
            Wnom += temp*x[i]
            Wden += temp
        end
    end
    return Wnom/Wden
end
#tauloc{S <: Real}(x::DataVector{S}; cval::Real=4.5)=tauloc(removeNA(x), cval=cval)

#Compute the tau measure of scale as described in Yohai and Zamar (JASA, 1988, 83, 406-413)
function tauvar{S <: Real}(x::Vector{S}; cval::Real=3.0)
    x     = x[:]
    s     = Rmath.qnorm(0.75)*mad(x)
    tloc  = tauloc(x)
    n     = length(x)
    W     = 0.0
    cval2 = cval*cval
    [W    += min(((x[i]-tloc)/s)*((x[i]-tloc)/s), cval2) for i=1:n]
    s*s*W/n
end
#tauvar{S <: Real}(x::DataVector{S}; cval::Real=3.0)=tauvar(removeNA(x), cval=cval)

#Outlier detection method using the boxplot rule. The ideal fourths are
#used to estimate the quartiles
function outbox{S <: Real}(x::Vector{S}; mbox::Bool=false, gval::Real=NaN, method::Bool=true)
    n    = length(x)
    temp = idealf(x, method=false)
    cl   = 0.0
    cu   = 0.0
    if mbox
        if isnan(gval)
            gval=(17.63*n-23.64)/(7.74*n-3.71)
        end
        cl = median(x, checknan=false) - gval*(temp.upper_quartile-temp.lower_quartile)
        cu = median(x, checknan=false) + gval*(temp.upper_quartile-temp.lower_quartile)
    elseif !mbox
        if isnan(gval)
            gval=1.5
        end
        cl = temp.lower_quartile - gval*(temp.upper_quartile-temp.lower_quartile)
        cu = temp.upper_quartile + gval*(temp.upper_quartile-temp.lower_quartile)
    end
    flag=falses(n)
    for i=1:n
        if (x[i].<cl) || (x[i].>cu)
            flag[i]=true
        end
    end
    vec=1:n
    outid  = sum(flag)>0?vec[flag]:[]
    keepid = vec[!flag]
    outval = x[flag]
    nout   = length(outid)
    if method && !mbox
        METHOD = "Outlier detection method using \nthe ideal-fourths based boxplot rule\n"
    elseif method && mbox
        METHOD = "Outlier detection method using \nthe ideal-fourths based boxplot rule\n(using the modification suggested by Carling (2000))\n"
    else
        METHOD = nothing
    end
    outOutput(outid, keepid, outval, nout, METHOD)
end
#outbox{S <: Real}(x::DataVector{S}; mbox=false, gval::Real=NaN, method=true)=outbox(removeNA(x), mbox=mbox, gval=gval, method=method)

#Compute adaptive kernel density estimate for univariate data
function akerd{S <: Real}(x::Vector{S}; hval::Real=NaN, aval::Real=0.5, op::Integer=1,
               fr::Real=0.8, pyhat::Bool=false, pts=NaN, plotit=true, xlab="",
               ylab="", title="", color="black", plottype="solid")
    xsort =  sort(x)
    if op == 1
        m =  mad(x)
        if m == 0.0
            temp = idealf(x, method=false)
            m = (temp.upper_quartile-temp.lower_quartile)/(Rmath.qnorm(0.75)-Rmath.qnorm(0.25))
        end
        if m  == 0.0
            m =  sqrt(winvar(x)./0.4129)
        end
        if m  == 0.0
            error("All measures of dispersion are equal to 0")
        end
        fhat  = rdplot(x, pyhat=true, plotit=false, fr=fr, pyhat=true)
        if m  >  0.0
            [fhat[i] = fhat[i]/(2.0*fr*m) for i=1:length(fhat)]
        end
    end
    if op == 2
        init = kde(x)
        fhat = init.density
        x    = init.x
    end
    n = length(x)
    if isnan(hval)
        sig  = std(x)
        temp = idealf(x, method=false)
        iqr  = (temp.upper_quartile-temp.lower_quartile)/1.34
        A    = min(sig, iqr)
        A    = A==0.0? winstd(x)/0.64:A
        hval = 1.06*A/n^0.2
    end
    gm = 0.0
    gm_int = 0
    nfhat = length(fhat)
    for i = 1:nfhat
        if fhat[i] > 0.0
            gm += log(fhat[i])
            gm_int += 1
        end
    end
    gm = exp(gm/gm_int)
    alam = zeros(n)
    [alam[i] = (fhat[i]/gm)^(0-aval) for i=1:nfhat]
    if isnan(pts)
        pts = x[:]
    end
    pts  = sort!(pts)
    dhat = akerd_loop(x, pts, hval, alam)
    if plotit
        p = FramedPlot()
        add(p, Curve(pts, dhat, "color", color, "type", plottype))
        if title != nothing
            setattr(p, "title", title)
        end
        if xlab != nothing
            setattr(p, "xlabel", xlab)
        end
        if ylab != nothing
            setattr(p, "ylabel", ylab)
        end
        Winston.tk(p)
    end
    if pyhat
        return dhat
    end
end
#akerd{S <: Real}(x::DataVector{S}; hval::Real=NaN, aval::Real=0.5,
 #                op::Integer=1, fr::Real=0.8, pyhat::Bool=false,
  #               pts=NaN, plotit=true, xlab="", ylab="", title="",
   #              color="black", plottype="solid")=akerd(removeNA(x),
    #             hval=hval, aval=aval, op=op, fr=fr, pyhat=pyhat, pts=pts,
     #            plotit=plotit, xlab=xlab, ylab=ylab, title=title, color=color,
      #           plottype=plottype)

#Expected frequency curve. fr controls amount of smoothing, theta is the azimuthal direction and
#phi the colatitude
function rdplot{S <: Real}(x::Vector{S}; fr::Real=NaN, plotit=true,
                           theta=50, phi=25, expand=0.5, pyhat=false, pts=NaN,
                           title="", xlab="", ylab="", color="black", plottype="solid")
    x = x[:]
    if isnan(fr)
        fr = 0.8
    end
    if isnan(pts)
        pts = x[:]
    end
    nx   = length(x)
    npts = length(pts)
    rmd  = [sum(near(x, pts[i], fr))*1.0 for i=1:npts]
    MAD  = mad!(x)
    if MAD != 0.0
        [rmd[i] = rmd[i]/(2*fr*MAD) for i=1:npts]
    end
    [rmd[i] = rmd[i]/nx i=1:npts]
    if plotit
        index = sortperm(pts);
        p     = FramedPlot()
        add(p, Curve(pts[index], rmd[index], "color", color, "type", plottype))
        if title != nothing
            setattr(p, "title", title)
        end
        if xlab != nothing
            setattr(p, "xlabel", xlab)
        end
        if ylab != nothing
            setattr(p, "ylabel", ylab)
        end
        Winston.tk(p)
    end
    if pyhat
        return rmd
    end
end
#rdplot{S <: Real}(x::DataVector{S}; fr::Real=NaN, plotit=true, theta=50, phi=25, expand=0.5,
 #                 pyhat=false, pts=NaN, title="", xlab="", ylab="", color="black", plottype="solid")=
#rdplot(removeNA(x), fr=fr, plotit=plotit, theta=theta, phi=phi, expand=expand,
 #   pyhat=pyhat, pts=pts, title=title, xlab=xlab, ylab=ylab, color=color, plottype=plottype)

#Determine which values in x are near pt
function near{S <: Real}(x::Vector{S}, pt::Real, fr::Real=1.0)
    m = mad!(x)
    if m == 0.0
        temp = idealf(x, method=false)
        m = (temp.upper_quartile-temp.lower_quartile)/(Rmath.qnorm(0.75)-Rmath.qnorm(.25))
    end
    if m == 0.0
        m = sqrt(winvar(x)/0.4129)
    end
    if m == 0.0
        error("All measures of dispersion are equal to 0")
    end
    n = length(x)
    dflag = falses(n)
    fr_m = fr*m
    #[dflag[i]=abs(x[i]-pt) <= fr_m ?true:false for i=1:n]
    for i = 1:n
        if abs(x[i]-pt) <= fr_m
            dflag[i] = true
        end
    end

    return dflag
end
#near{S <: Real}(x::DataVector{S}, pt::Real, fr::Real=1.0)=near(removeNA(x), pt, fr)

#Computing the confidence interval for the median using the
#Hettmansperger and Sheather method.
function sint{S <: Real}(x::Vector{S}; alpha::Real=0.05, method::Bool=true)
    n  = length(x)
    k  = Rmath.qbinom(alpha/2.0, n, 0.5)
    gk = Rmath.pbinom(length(x)-k, n, .5) - Rmath.pbinom(k-1, n, .5)
    if gk >= (1 - alpha)
        gkp1 = Rmath.pbinom(n-k-1, n, .5) - Rmath.pbinom(k, n, .5)
        kp   = k + 1
    else
        k    = k - 1
        gk   = Rmath.pbinom(n-k, n, .5)   - Rmath.pbinom(k-1, n, .5)
        gkp1 = Rmath.pbinom(n-k-1, n, .5) - Rmath.pbinom(k, n, .5)
        kp   = k+1
    end
    xsort=sort(x)
    nmk  = n-k
    nmkp = nmk+1
    ival = (gk-1+alpha)/(gk-gkp1)
    lam  = ((n-k)*ival)/(k+(n-2*k)*ival)
    low  = lam*xsort[kp]+(1-lam)*xsort[k]
    hi   = lam*xsort[nmk]+(1-lam)*xsort[nmkp]
    if method
        if duplicated(x)
            METHOD="Confidence interval for the median\nDuplicate values detected; hdpb() might have more power\n"
        else
            METHOD="Confidence interval for the median\n"
        end
    else
        METHOD=nothing
    end
    output=testOutput()
    output.method=METHOD
    output.ci=[low, hi]
    output
end
#sint{S <: Real}(x::DataVector{S}; alpha::Float64=0.05, method=true)=sint(removeNA(x), alpha=alpha, method=method)

#Computing standard error of the median using a method recommended by McKean and Shrader (1984)
function msmedse{S <: Real}(x::Vector{S})
    if duplicated(x)
        warn("Tied values detected. Estimate of standard error might be highly inaccurate, even with n large")
    end
    y = sort(x)
    n = length(y)
    av::Int = round((n+1)./2-Rmath.qnorm(.995).*sqrt(n/4))
    if av == 0
        av = 1
    end
    top::Int = n-av+1
    abs((y[top]-y[av])/(2*Rmath.qnorm(.995)))
end
#msmedse{S <: Real}(x::DataVector{S})=msmedse(removeNA(x))

#Compute a 1-alpha confidence interval for p, the probability of success for a binomial dist. using Pratt's method
#y is a vector of 1's and 0's, x is the number of successes observed among n trials.

function binomci(x::Int, n::Int; alpha::Real=0.05)
    if x > n
        error("x must be smaller than or equal to n")
    elseif x < 0
        error("x cannot be negative")
    elseif n == 1
        error("Something is wrong: number of observations is only 1")
    end
    if x != n && x != 0
        z     = Rmath.qnorm(1-alpha/2)
        A     = ((x+1)/(n-x)).*((x+1)/(n-x))
        B     = 81.*(x+1).*(n-x)-9.*n-8
        C     = (0-3).*z.*sqrt(9.*(x+1).*(n-x).*(9*n+5-z.*z)+n+1)
        D     = 81.*(x+1).*(x+1)-9.*(x+1).*(2+z*z)+1
        E     = 1+A.*((B+C)./D).*((B+C)./D).*((B+C)./D)
        upper = 1/E
        A     = (x./(n-x-1)).*(x./(n-x-1))
        B     = 81.*x.*(n-x-1)-9.*n-8
        C     = 3.*z.*sqrt(9.*x.*(n-x-1).*(9.*n+5-z.*z)+n+1)
        D     = 81.*x.*x-9.*x.*(2+z.*z)+1
        E     = 1+A.*((B+C)./D).^3
        lower = 1/E
    end
    if x == 0
        lower = 0.0
        upper = 1.0-alpha.^(1/n)
    end
    if x == 1
        lower = 1-(1-alpha/2).^(1/n)
        upper = 1-(alpha/2).^(1/n)
    end
    if x == (n-1)
        lower = (alpha/2).^(1/n)
        upper = (1-alpha/2).^(1/n)
    end
    if x == n
        lower = alpha.^(1/n)
        upper = 1
    end
    phat=x/n
    binomciOutput(phat, [lower, upper], n)
end
binomci(x::Vector{Int}; alpha::Real=0.05)=
sum([x[i]!=0 && x[i]!=1 for i=1:length(x)])>0?error("Data must be 0's and 1's"):binomci(length(x.==1), length(x), alpha=alpha)
#binomci(x::DataVector{Int}; alpha::Real=0.05)=
#sum([x[i]!=0 && x[i]!=1 for i=1:length(x)])>0?error("Data must be 0's or 1's"):binomci(removeNA(x), alpha=alpha)

#Computing a 1-alpha confidence interval for p, the probability of
#success for a binomial distribution, using a generalization of the
#Agresti-Coull  method that was studied by Brown, Cai DasGupta
#(Annals of Statistics, 2002, 30, 160-201.)
function acbinomci(x::Int, n::Int; alpha::Real=0.05)
    if x > n
        error("x must be smaller than or equal to n")
    elseif x < 0
        error("x cannot be negative")
    elseif n == 1
        error("Something is wrong: number of observations is only 1")
    end
    if x != n && x != 0
        cr    = Rmath.qnorm(1-alpha/2)
        ntil  = n+cr.*cr
        ptil  = (x+cr.*cr./2)./ntil
        lower = ptil-cr.*sqrt(ptil.*(1-ptil)./ntil)
        upper = ptil+cr.*sqrt(ptil.*(1-ptil)./ntil)
    end
    if x == 0
        lower = 0.0
        upper = 1.0-alpha.^(1/n)
    end
    if x == 1
        lower = 1-(1-alpha/2).^(1/n)
        upper = 1-(alpha/2).^(1/n)
    end
    if x == (n-1)
        lower = (alpha/2).^(1/n)
        upper = (1-alpha/2).^(1/n)
    end
    if x == n
        lower = alpha.^(1/n)
        upper = 1
    end
    phat = x/n
    binomciOutput(phat, [lower, upper], n)
end
acbinomci(x::Vector{Int}; alpha::Real=0.05)=
sum([x[i]!=0 && x[i]!=1 for i=1:length(x)])>0?error("Data must be 0's and 1's"):acbinomci(length(x.==1), length(x), alpha=alpha)
#acbinomci(x::DataVector{Int}; alpha::Real=0.05)=
#sum([x[i]!=0 && x[i]!=1 for i=1:length(x)])>0?error("Data must be 0's or 1's"):acbinomci(removeNA(x), alpha=alpha)

#Extension of Stein's method based on the trimmed mean.
function stein1_tr{S <: Real}(x::Vector{S}, del::Real; alpha::Real=0.05, pow::Real=0.8, tr::Real=0.2)
    if tr <0 || tr >= .5
        error("Argument tr must be between 0 and .5")
    end
    n =length(x)
    g::Int    = floor(tr*n)
    df::Int   = n-2*g-1
    t1        = Rmath.qt(pow, df)
    t2        = Rmath.qt(alpha./2, df)
    dv        = (del/(t1-t2))*(del/(t1-t2))
    nvec::Int = floor(trimse(x, tr=tr)./dv)+1
    n > nvec ? n : nvec
end

function stein1_tr{S <: Real}(x::Array{S, 2}, del::Real; alpha::Real=0.05, pow::Real=0.8, tr::Real=0.2)
    if tr<0||tr>=.5
        error("Argument tr must be between 0 and .5")
    end
    if size(x, 2)==1 || size(x, 1)==1
        temp=zeros(length(x))
        temp[:]=x
        return stein1_tr(temp, del, alpha=alpha, pow=pow, tr=tr)
    else
        J     = size(x,2)
        ic    = 0
        ntest::Integer = (J*J-J)/2
        n     = size(x, 1)
        g::Integer = floor(tr*n)
        df    = n-2g-1
        tdiff = Rmath.qt(pow, df)-Rmath.qt(alpha/2.0/ntest, df)
        dv    = (del/tdiff)*(del/tdiff)
        N     = zeros(Int, ntest)
        if ntest > 1
            for i = 1:J
                for j = 1:J
                    if i < j
                        N[ic+=1] = floor(trimse(x[:,i]-x[:,j], tr=tr)/dv)+1
                    end
                end
            end
            return max(vcat(N, n))
        elseif ntest == 1
            return int(floor(trimse(x[:,1], tr=tr)/dv)+1)
        end
    end
end
#stein1_tr{S <: Real}(x::DataArray{S, 2}, del::Real; alpha::Real=0.05, pow::Real=0.8, tr::Real=0.2)=
 #   stein1_tr(deleterows!(x), del, alpha=alpha, pow=pow, tr=tr)

#Extension of the second stage of Stein's method when performing all pairwise comparisons
#among J dependent groups.
function stein2_tr{S <: Real, T <: Real}(x::Vector{S}, y::Vector{T}; alpha::Real=0.05, tr::Real=0.2, method::Bool=true)
    if tr < 0 || tr >= .5
        error("Argument tr must be between 0 and .5")
    end
    m       = vcat(x, y)
    METHOD  = method? "Extension of the 2nd stage of Stein's method based on the trimmed mean\n":nothing
    output  = testOutput()
    output.method    = METHOD
    output.statistic = sqrt(length(m))*(1-2tr)*tmean(m, tr=tr)/winstd(m)
    output.crit         = Rmath.qt(1-alpha/2, floor(length(x)-2floor(tr*length(x))-1))
    return output
end

function stein2_tr{S <: Real, T <: Real}(x::Array{S}, y::Array{T}; alpha::Real=0.05, tr::Real=.2, method::Bool=true)
    if tr <0 || tr >= .5
        error("Argument tr must be between 0 and .5")
    end
    if size(x, 2) != size(y, 2)
        error("Number of variables in x not equal to number of variables in y")
    end
    if size(x, 2) == 1 || size(x, 1) ==1
        return stein2_tr(x[:], del, alpha=alpha, pow=pow, tr=tr, method=method)
    else
        g::Int  = floor(tr*size(x,1))
        df::Int = size(x,1)-2.*g-1
        ic::Int = 0
        J::Int  = size(x,2)
        ntest::Int = (J.*J-J)./2
        m       = vcat(x, y)
        nm      = size(m, 1)
        test=DataFrame({Integer, Integer, Float64}, ["grp1", "grp2", "test_stat"], ntest)
        if ntest > 1
            for i=1:J
                for j=1:J
                    if i<j
                        temp  = m[:,i]-m[:,j]
                        test[ic+=1, 1], test[ic, 2] = i, j
                        test[ic, 3] = sqrt(nm)*(1-2tr)*tmean(temp, tr=tr)/winstd(temp)
                    end
                end
            end
        else ntest == 1
            test[1,1], test[1, 2] = 1, 2
            test[1, 3] = test[ic, 3] = sqrt(nm)*(1-2tr)*tmean(m[:,1], tr=tr)/winstd(m[:,1])
        end
        if method
            METHOD="Extension of the 2nd stage of Stein's method based on the trimmed mean\n"
        else
            METHOD=nothing
        end
        output = testOutput()
        output.method    = METHOD
        output.statistic = test
        output.crit      = Rmath.qt(1-alpha/2./ntest, df)
        return output
    end
end

#stein2_tr{S <: Real, T <: Real}(x::DataVector{S}, y::DataVector{S}; alpha::Real=0.05, tr::Real=0.2, method::Bool=true)=
#    stein2_tr(removeNA(x), removeNA(y), alpha=alpha, tr=tr, method=method)
#stein2_tr{S <: Real, T <: Real}(x::Array{Int}, y::Array{Int}; alpha=0.05, tr=0.2, method=true)=
#    stein2_tr(convert(Array{Float64}, x), convert(Array{Float64}, y), alpha=alpha, tr=tr, method=method)

#Compute 1-alpha confidence interval for the median using the Hettmansperger-Sheather interpolation method
function sintv2{S <: Real}(x::Vector{S}; alpha::Real=0.05, nullval::Real=0.0, method::Bool=true)
    ci      = sint(x, alpha=alpha, method=false).ci
    alph    = [1:99]./100
    irem    = zeros(Int, 1)
    for i = 1:99
        irem  = i
        chkit = sint(x, alpha=alph[i], method=false).ci
        if chkit[1].>nullval || chkit[2].<nullval
            break
        end
    end
    pval=irem/100
    if pval.<=0.1
        iup  = (irem+1)/100
        alph = seq(0.001, iup, 0.001)
        for i = 1:length(alph)
            pval  = alph[i]
            chkit = sint(x, alpha=alph[i], method=false).ci
            if chkit[1] > nullval || chkit[2] < nullval
                break
            end
        end
    end
    if pval <= 0.001
        alph = seq(.0001,.001,.0001)
        for i = 1:length(alph)
            pval  = alph[i]
            chkit = sint(x, alpha=alph[i], method=false).ci
            if chkit[1] > nullval || chkit[2] < nullval
                break
            end
        end
    end
    if method
        if duplicated(x)
            METHOD="Confidence interval for the median\nusing the Hettmansperger-Sheather interpolation method\nDuplicate values detected; hdpb() might have more power\n"
        else
            METHOD="Confidence interval for the median\nusing the Hettmansperger-Sheather interpolation method\n"
        end
    else
        METHOD=nothing
    end
    output = testOutput()
    output.method = METHOD
    output.ci     = ci
    output.p      = pval
    output
end

#sintv2{S <: Real}(x::DataVector{S}; alpha::Real=0.05, nullval::Real=0.0, method::Bool=true)=
 #   sintv2(removeNA(x), alpha=alpha, nullval=nullval, method=method)

#Evaluate Huber's Psi function for each value in the vector x.
#The bending constant defaults to 1.28
hpsi{S <: Real}(x::Vector{S}, bend::Real=1.28)=[abs(x[i])<=bend?x[i]:bend*sign(x[i]) for i=1:length(x)]
#hpsi{S <: Real}(x::DataVector{S}, bend::Real=1.28)= hpsi(removeNA(x), bend=bend)

#Compute one-step M-estimator of location using Huber's Psi.
#The default bending constant is 1.28
function onestep{S <: Real}(x::Vector{S}, bend::Real=1.28)
    MED = median(x, checknan=false)
    MAD = mad(x)
    n   = length(x)
    y   = [(x[i]-MED)/MAD for i = 1:n]
    A   = sum(hpsi(y, bend))
    B   = 0
    [B  += abs(y[i]) <= bend ? 1:0 for i=1:n]
    return MED + MAD*A/B
end
#onestep{S <: Real}(x::DataVector{S}, bend::Real=1.28)=onestep(removeNA(x), bend=bend)

#Compute a bootstrap, .95 confidence interval for the
#measure of location corresponding to the argument est.
#By default, a one-step
#M-estimator of location based on Huber's Psi is used.
#The default number of bootstrap samples is nboot=500
#
#nv=null value when  computing a p-value

function onesampb{S <: Real}(x::Vector{S}; est::Function=onestep, alpha::Real=0.05, nboot::Integer=2000, seed::Union(Int, Bool)=2, nv::Real=0)
    if iseltype(seed, Int)
        srand(seed)
    elseif seed
        srand(2)
    end
    n        = length(x)
    bvec     = zeros(nboot)
    temp     = zeros(n)
    randid   = rand(1:n, n*nboot)
    for i = 1:nboot
        [temp[j] = x[randid[(i-1)*n+j]] for j=1:n]
        bvec[i]  = est(temp)
    end
    low      = round((alpha./2).*nboot) + 1
    up       = nboot-low + 1
    pv       = mean(bvec.>nv)+0.5.*mean(bvec.==nv)
    pv       = 2.*min(pv, 1-pv)
    bvec     = sort!(bvec)
    estimate = est(x)
    output   = testOutput()
    output.estimate = estimate
    output.ci       = [bvec[low], bvec[up]]
    output.p        = pv
    output
end
#onesampb{S <: Real}(x::DataVector{S}; est::Function=onestep, alpha::Real=0.05, nboot::Integer=2000, seed::Union(Int, Bool)=2, nv::Real=0)=
 #   onesampb(removeNA(x), est=est, alpha=alpha, nboot=nboot, seed=seed, nv=nv)

#MOM-estimator of location. The default bending constant is 2.24
function mom{S <: Real}(x::Vector{S}; bend::Real=2.24)
    n    = length(x)
    flag = trues(n)
    MED  = median(x, checknan=false)
    MAD  = mad(x)
    for i = 1:n
        if x[i] > MED+bend*MAD || x[i] < MED-bend*MAD
            flag[i]=false
        end
    end
    mean(x[flag])
end
#mom{S <: Real}(x::DataVector{S}; bend=2.24)=mom(removeNA(x), bend=bend)

function mom!{S <: Real}(x::Vector{S}; bend::Real=2.24)
    n    = length(x)
    flag = trues(n)
    MED  = median!(x, checknan=false)
    MAD  = mad!(x)
    for i = 1:n
        if x[i] > MED+bend*MAD || x[i] < MED-bend*MAD
            flag[i]=false
        end
    end
    mean(x[flag])
end

#Compute a bootstrap, .95 confidence interval for the
#MOM-estimator of location based on Huber's Psi.
#The default number of bootstrap samples is nboot=500
function momci{S <: Real}(x::Vector{S}; bend::Real=2.24, alpha::Real=0.05, nboot::Integer=2000, seed::Union(Int, Bool)=2, method=true)
    if iseltype(seed, Int)
        srand(seed)
    elseif seed
        srand(2)
    end
    n    = length(x)
    bvec = zeros(nboot)
    temp = zeros(n)
    randid=rand(1:n, n*nboot)
    for i = 1:nboot
        [temp[j] = x[randid[(i-1)*n+j]] for j=1:n]
        bvec[i]=mom!(temp, bend=bend)
    end
    low  = round((alpha./2).*nboot) + 1
    up   = nboot-low + 1
    bvec = sort!(bvec)
    METHOD = method? "Bootstrap .95 confidence interval for the MOM-\nestimator of location based on Huber's Psi\n":nothing
    output = testOutput()
    output.method = METHOD
    output.ci     = [bvec[low], bvec[up]]
    output
end
#momci{S <: Real}(x::DataVector{S}; bend::Real=2.24, alpha::Real=0.05, nboot::Integer=2000, seed::Union(Int, Bool)=2, method=true)=
 #   momci(removeNA(x), bend=bend, alpha=alpha, nboot=nboot, seed=seed, method=method)

#Contaminated normal distribution
function cnorm(n::Integer; epsilon::Real=0.1, k::Real=10)
    if epsilon > 1
        error("epsilon must be less than or equal to 1")
    elseif epsilon < 0
        error("epsilon must be greater than or equal to 0")
    end
    if k <= 0
        error("k must be greater than 0")
    end
    output  = zeros(n)
    epsilondiff = 1.0 - epsilon
    [output[i] = rand() > epsilondiff ? k*randn():randn() for i=1:n]
    return output
end

#   Compute a 1-alpha confidence interval for
#   a trimmed mean.
#
#   The default number of bootstrap samples is nboot=2000
#
#   win is the amount of Winsorizing before bootstrapping
#   when WIN=T.
#
#   Missing values are automatically removed.
#
#  nv is null value. That test hypothesis trimmed mean equals nv
#
#  plotit=TRUE gives a plot of the bootstrap values
#  pop=1 results in the expected frequency curve.
#  pop=2 kernel density estimate    NOT IMPLEMENTED
#  pop=3 boxplot                    NOT IMPLEMENTED
#  pop=4 stem-and-leaf              NOT IMPLEMENTED
#  pop=5 histogram
#  pop=6 adaptive kernel density estimate.
#
#  fr controls the amount of smoothing when plotting the bootstrap values
#  via the function rdplot. fr=NA means the function will use fr=.8
#  (When plotting bivariate data, rdplot uses fr=.6 by default.)

function trimpb{S <: Real}(x::Vector{S}; tr::Real=0.2, alpha::Real=0.05, nboot::Integer=2000,
                win::Union(Bool, Real)=false, plotit::Bool=false, pop::Int=1, nullval::Real=0.0,
                xlab="X", ylab="Density", fr::Union(Nothing, Real)=nothing, seed::Union(Bool, Integer)=2,
                method::Bool=true)
    if iseltype(win, Bool)
        if win
            x=winval(x, tr=0.1)
        end
    elseif win>tr
        error("The amount of Winsorizing must be <= to the amount of trimming")
    else
        x=winval(x, tr=win)
    end
    crit=alpha/2.0
    icl=round(crit*nboot)+1
    icu=nboot-icl
    if iseltype(seed, Bool)
        if seed
            srand(2)
        end
    else
        srand(seed)
    end
    n=length(x)
    bvec=zeros(nboot)
    randid=rand(1:n, n*nboot)
    for i=1:nboot
        temp=zeros(n)
        for j=1:n
            temp[j]=x[randid[(i-1)*n+j]]
        end
        bvec[i]=tmean!(temp, tr=tr)
    end
    bvec=sort!(bvec)
    pval1=pval2=0.0
    for i=1:nboot
        if bvec[i]<nullval
            pval1+=1./nboot
        end
        if bvec[i]==nullval
            pval2+=0.5/nboot
        end
    end
    pval=2*min(pval1+pval2, 1-pval1-pval2)
    ci=[bvec[icl], bvec[icu]]
    if method
        METHOD::String="Compute a 1-alpha confidence interval for a trimmed mean using the bootstrap percentile method.\n"
    else
        METHOD=nothing
    end
    if plotit
        if pop==1
            if fr==nothing
                rdplot(bvec, fr=0.6, xlab=xlab, ylab=ylab)
            else
                rdplot(bvec, fr=fr, xlab=xlab, ylab=ylab)
            end
        elseif pop==5
            p=FramedPlot()
            add(p, Histogram(hist(bvec)[2], 2))
            if xlab!=nothing
                setattr(p, "xlabel", xlab)
            end
            if ylab!=nothing
                setattr(p, "ylabel", ylab)
            end
            Winston.tk(p)
        elseif pop==6
            akerd(bvec, xlab=xlab, ylab=ylab)
        end
    end
    output = testOutput()
    output.method = METHOD
    output.ci     = ci
    output.p      = pval
    output
end

#trimpb{S <: Real}(x::DataVector{S}; tr::Real=0.2, alpha::Real=0.05, nboot::Integer=2000,
 #               win::Union(Bool, Real)=false, plotit::Bool=false, pop::Int=1, nullval::Real=0.0,
  #              xlab="X", ylab="Density", fr::Union(Nothing, Real)=nothing, seed::Union(Bool, Integer)=2,
   #             method::Bool=true)=trimpb(removeNA(x), tr=tr, alpha=alpha, nboot=nboot, win=win, plotit=plotit,
    #            pop=pop, nullval=nullval, xlab=xlab, ylab=ylab, fr=fr, seed=seed, method=method)

#  Compute a 1-alpha confidence interval for the trimmed mean
#  using a bootstrap percentile t method.
#
#  The default amount of trimming is tr=.2
#  side=T, for true,  indicates the symmetric two-sided method
#
#
#  Side=F yields an equal-tailed confidence interval
#
#
#  NOTE: p.value is reported when side=T only.
#

function trimcibt{S <: Real}(x::Vector{S}; tr::Real=0.2, alpha::Real=0.05, nboot::Integer=2000, side::Bool=true,
                             plotit::Bool=false, op::Integer=1, nullval::Real=0, seed::Union(Bool, Int)=2, method::Bool=true)
    if iseltype(seed, Bool)
        if seed
            srand(2)
        end
    else
        srand(seed)
    end
    n=length(x)
    test=(tmean(x, tr=tr)-nullval)./trimse(x, tr=tr)
    randid=rand(1:n, n*nboot)
    tempout=trimcibt_loop(x, n, nboot, tr, side, randid, test)
    if side
        tval=tempout[1]
        pval=tempout[2]
    end
    icrit=floor((1-alpha)*nboot+0.5)
    ibot=round(alpha*nboot/2)+1
    itop=nboot-ibot-1
    ci=zeros(2)
    if method && side
        METHOD="Bootstrap .95 confidence interval for the trimmed mean\nusing a bootstrap percentile t method\n"
    elseif method && !side
        METHOD="Bootstrap .95 confidence interval for the trimmed mean\nusing a bootstrap percentile t method\n[NOTE: p value is computed only when side=true]\n"
    else
        METHOD=nothing
    end
    if !side
        if plotit
            if op==1
                akerd(tempout)
            elseif op==2
                rdplot(tempout)
            end
        end
        ci[1]=tmean(x, tr=tr)-tempout[itop]*trimse(x, tr=tr)
        ci[2]=tmean(x, tr=tr)-tempout[ibot]*trimse(x, tr=tr)
        output = testOutput()
        output.method = METHOD
        output.estimate = tmean(x, tr=tr)
        output.ci = ci
        output.statistic = test
        return output
    else
        if plotit
            if op==1
                akerd(tval)
            elseif op==2
                rdplot(tval)
            end
        end
        ci[1]=tmean(x, tr=tr)-tval[icrit]*trimse(x, tr=tr)
        ci[2]=tmean(x, tr=tr)+tval[icrit]*trimse(x, tr=tr)
        output = testOutput()
        output.method = METHOD
        output.estimate = tmean(x, tr=tr)
        output.ci = ci
        output.statistic = test
        output.p = pval
        return output
    end
end
#trimcibt{S <: Real}(x::DataVector{S}; tr::Real=0.2, alpha::Real=0.05, nboot::Integer=2000, side::Bool=true,
 #                   plotit::Bool=false, op::Integer=1, nullval::Real=0, seed::Union(Bool, Int)=2, method::Bool=true) =
  #                  trimcibt(removeNA(x), tr=tr, alpha=alpha, nboot=nboot, side=side, plotit=plotit, op=op, nullval=nullval,
   #                 seed=seed, method=method)

#   Compute bootstrap estimate of the standard error of the
#   estimator est
#   The default number of bootstrap samples is nboot=1000

function bootse{S <: Real}(x::Vector{S}; nboot::Integer=1000, est::Function=median, seed::Union(Bool, Integer)=2)
   if iseltype(seed, Bool)
        if seed
            srand(2)
        end
    else
        srand(seed)
    end
    n=length(x)
    temp=zeros(n)
    bvec=zeros(nboot)
    randid=rand(1:n, n*nboot)
    for i=1:(nboot*n)
        if (i%n)!=0
            temp[i%n]=x[randid[i]]
        else
            temp[n]=x[randid[i]]
            bvec[div(i, n)]=est(temp)

        end
    end
    return std(bvec)
end

#bootse{S <: Real}(x::DataVector{S}; nboot::Integer=1000,
 #                 est::Function=median, seed::Union(Bool, Integer)=2)=
  #                bootse(removeNA(x), nboot=nboot, est=est, seed=seed)

#   Compute a .95 confidence interval for Pearson's correlation coefficient.
#
#   This function uses an adjusted percentile bootstrap method that
#   gives good results when the error term is heteroscedastic.
function pcorb{S <: Real, T <: Real}(x::Vector{S}, y::Vector{T}; seed::Union(Bool, Integer)=2, plotit::Bool=false)
   if iseltype(seed, Bool)
        if seed
            srand(2)
        end
    else
        srand(seed)
    end
    n=length(x)
    randid=rand(1:n, n*599)
    tempx=zeros(n)
    tempy=zeros(n)
    bvec=zeros(599)
    for i=1:(599*n)
        if (i%n)!=0
            tempx[i%n], tempy[i%n]=x[randid[i]], y[randid[i]]
        else
            tempx[n], tempy[n]=x[randid[i]], y[randid[i]]
            bvec[div(i, n)]=cor!(tempx, tempy)
        end
    end
    if n >= 250
        ilow, ihi = 15, 584
    elseif n >= 180
        ilow, ihi = 14, 585
    elseif n >= 80
        ilow, ihi = 11, 588
    elseif n >= 40
        ilow, ihi = 8, 592
    else
        ilow, ihi = 7, 593
    end
    bvec=sort!(bvec)
    if plotit
        akerd(bvec, title="Distribution of Bootstrap Pearson Correlation Coefficients")
    end
    r=cor!(x, y)
    output = testOutput()
    output.estimate = r
    output.ci = [bvec[ilow], bvec[ihi]]
    output
end
#pcorb{S <: Real, T <: Real}(x::DataVector{S}, y::DataVector{T}; nboot::Integer=2000, seed::Union(Bool, Integer)=2)
#= pcorb(removeNAVector{Int}, y::Vector{Float64}; nboot=2000, seed::Union(Bool, Int)=2)=

# Test the hypothesis of independence between x and y by
# testing the hypothesis that the regression surface is a horizontal plane.
# Stute et al. (1998, JASA, 93, 141-149).
#
#  flag=1 gives Kolmogorov-Smirnov test statistic
#  flag=2 gives the Cramer-von Mises test statistic
#  flag=3 causes both test statistics to be reported.
#
#  tr=0 results in the Cramer-von Mises test statistic when flag=2
#      With tr>0, a trimmed version of the test statistic is used.
#
function indt(x::AbstractArray, y::AbstractArray; flag::Int=1, nboot::Int=599,
              tr::Float64=0.0, seed::Union(Bool, Int)=2, method=true)
    if ndims(x)==1
        x=reshape(x,length(x),1)
    end
    if length(findin(flag, 1:3))==0
        error("flag must be set to 1, 2, or 3")
    end
    n=size(x,1)
    np=size(x,2)
    y=reshape(y, length(y), 1)
    if length(y)!=n
        error("Incondistent dimensions of x and y: number of x must match number of y")
    end
    mflag=indt_mflag(x)
    yhat=mean(y)
    res=zeros(n)
    [res[i]=y[i]-yhat for i=1:n]
    if iseltype(seed, Bool)
        if seed
            srand(2)
        end
    else
        srand(seed)
    end
    data=(rand(nboot, n)-0.5).*sqrt(12)
    rvalb=zeros(n, nboot)
    const sqrtn=sqrt(n)
    [rvalb[:,i]=regts1(data[i,:], yhat, res, mflag, x, 0) for i=1:nboot]
    [[rvalb[i,j]=abs(rvalb[i,j])/sqrtn for i=1:n] for j=1:nboot]

    dstatb=zeros(nboot)
    [dstatb[i]=max(rvalb[:,i]) for i=1:nboot]

    [[rvalb[i,j]=rvalb[i,j].*rvalb[i,j] for i=1:n] for j=1:nboot]
    wstatb=mean(rvalb, 1)
    rval=regts1(fill(1.0, n), yhat, res, mflag, x, tr)./sqrtn

    dstat=pval_d=wstat=pval_w=nothing
    if flag==1 || flag==3
        [rval[i]=abs(rval[i]) for i=1:n]
        dstat=max(rval)
        pval_d=0.0
        pval_d=sum(dstat.>=dstatb)./nboot
        pval_d=1.0-pval_d
    end
    if flag==2 || flag==3
        [rval[i]=rval[i].*rval[i] for i=1:n]
        wstat=tmean(rval, tr=tr)
        pval_w=0.0
        pval_w=sum(wstat.>=wstatb)./nboot
        pval_w=1.0-pval_w
    end
    if method
        METHOD::String="Test whether x and y are independent by testing the hypothesis\nthat the regression surface is a horizontal plane.\n"
    else
        METHOD=nothing
    end
    indtOutput(
        METHOD,
        dstat,
        pval_d,
        wstat,
        pval_w,
        flag
        )
end

function indirectTest{S <: Real, T <: Real, W <: Real}(dv::Vector{S}, iv::Vector{T}, m::Vector{W};
            nboot::Integer=5000, alpha::Real=0.05, scale::Bool=false, seed::Union(Bool, Int)=2, plotit::Bool=false)
    if iseltype(seed, Bool)
        if seed
            srand(2)
        end
    else
        srand(seed)
    end
    n = length(iv)
    randid  = rand(1:n, n*nboot)
    bvec    = sort!(bootindirect(iv, dv, m, nboot))
    bbar    = mean(bvec)
    bootci  = [bvec[round(alpha*nboot/2)], bvec[nboot - round(alpha*nboot/2) + 1]]
    bootest = mean(bvec)
    bootse  = std(bvec)
    p       = mean(bvec .<0) + 0.5*mean(bvec .==0)
    p       = 2*min(p, 1-p)
    data    = DataFrame(iv, m, dv)
    regfit1 = coeftable(lm( :(x3 ~ x1     ), data))
    regfit2 = coeftable(lm( :(x2 ~ x1     ), data))
    regfit3 = coeftable(lm( :(x3 ~ x1 + x2), data))
    regfit  = rbind(regfit1[2,:], regfit2[2,:], regfit3[2:3,:])

    estimate  = regfit2[2,1]*regfit3[3,1]
    sobel_se  = sqrt(regfit[4, 1]*regfit[4, 1]*regfit[2, 2]*regfit[2, 2]+
                     regfit[2, 1]*regfit[2, 1]*regfit[4, 2]*regfit[4, 2]+
                     regfit[2, 2]*regfit[2, 2]*regfit[4, 2]*regfit[4, 2])
    sobel     = DataFrame(estimate,
                          sobel_se,
                          estimate/sobel_se,
                          estimate - Rmath.qnorm(.975)*sobel_se,
                          estimate + Rmath.qnorm(.975)*sobel_se,
                          2*(1-Rmath.pnorm(abs(estimate/sobel_se))))
    #return nboot, n, regfit, sobel, bootest, bootci, p
    if plotit
        dens = kde(bvec)
        plot(dens.x, dens.density)
    end
    output = indirectTestOutput()
    output.nboot   = nboot
    output.n       = n
    output.sobel   = sobel
    output.regfit  = regfit
    output.bootest = bootest
    output.bootci  = bootci
    output.p       = p
    output.bootse  = bootse
    output
end

#   Compute the Winsorized correlation between x and y.
#
#   tr is the amount of Winsorization
#   This function also returns the Winsorized covariance
function wincor{S <: Real, T <: Real}(x::Vector{S}, y::Vector{T}; tr::Real=0.2)
    n = length(x)
    if n != length(y)
        error("`x` and `y` must agree in length")
    end
    g::Integer = floor(tr*n)
    xvec = winval(x, tr=tr)
    yvec = winval(y, tr=tr)
    wcor = cor(xvec, yvec)
    wcov = cov(xvec, yvec)

    if sum(x.==y) != n
        test = wcor*sqrt((n - 2)/(1 - wcor*wcor))
        sig  = 2*(1 - Rmath.pt(abs(test), n-2*g-2))
        return wcor, wcov, sig, n
    else
        return wcor, wcov, n
    end
end

#
#  Compare the trimmed means of two dependent random variables
#  using the data in x and y.
#  The default amount of trimming is 20%
#
#  Missing values (values stored as NA) are not allowed.
#
#  A confidence interval for the trimmed mean of x minus the
#  the trimmed mean of y is computed and returned in yuend$ci.
#  The significance level is returned in yuend$siglevel
function yuend{S <: Real, T <: Real}(x::Vector{S}, y::Vector{T}; tr::Real=0.2, alpha::Real=0.05, method::Bool=true)
    n = length(x)
    if n != length(y)
        error("`x` and `y` must agree in length")
    end
    h1::Integer = n - 2*floor(tr*n)
    q1 = (n - 1)*winvar(x, tr=tr)
    q2 = (n - 1)*winvar(y, tr=tr)
    q3 = (n - 1)*wincor(x, y, tr=tr)[2]
    df = h1 - 1
    se = sqrt((q1 + q2 - 2*q3)/(h1*(h1-1)))
    crit = Rmath.qt(1 - alpha/2, df)
    dif = tmean(x, tr=tr) - tmean(y, tr=tr)
    confint = [dif - crit*se, dif + crit*se]
    test = dif/se
    p = 2*(1 - Rmath.pt(abs(test), df))
    if method
        METHOD="Comparing the trimmed means of two dependent variables.\n"
    else
        METHOD=nothing
    end
    output = testOutput()
    output.method = METHOD
    output.ci = confint
    output.p = p
    output.estimate = dif
    output.se = se
    output.statistic = test
    output.n = n
    output.df = df
    output
end

#  A heteroscedastic one-way ANOVA for trimmed means using a generalization of Welch's method.

function t1way{S <: Real}(x::Array{S, 2}; tr::Real=0.2, method::Bool=true)
    n = size(x, 1)
    g = [1:size(x, 2)]
    grp = rep(g, rep(n, size(x, 2)))
    x = x[:]
    t1waycore(x, grp, tr, method)
end

function t1way{S <: Real}(x::Vector{S}, grp::Vector; tr::Real=0.2, method::Bool=true)
    g = unique(grp)
    grpcopy = [find(g.==grp[i])[1] for i=1:length(grp)]
    t1waycore(x, grpcopy, tr, method)
end

function pbos{S <: Real}(x::Vector{S}; beta::Real=0.2)
    temp    = sort( abs( x - median(x) ))
    nval    = length( x )
    omhatid::Integer = floor( (1 - beta)*nval )
    omhatx  = temp[ omhatid ]
    psi     = ( x - median(x) )./ omhatx
    i1      = length(psi[ psi .< -1 ])
    i2      = length(psi[ psi .> 1 ])
    sx      = 0.0
    [ sx += psi[i] < -1 ? 0 : [ psi[i] > 1 ? 0 : x[i] ]  for i=1:nval ]
    return ( sx  + omhatx * (i2 - i1))/(nval - i1 - i2)
end

#Compute the percentage bend correlation between x and y
#beta is the bending constant for omega sub N.
function pbcor{S <: Real, T <: Real}(x::Vector{S}, y::Vector{T}; beta::Real=0.2)
    nval = length(x)
    if length(y) != nval
        error("x and y do not agree in length.")
    end
    temp    = sort( abs( x - median(x) ))
    omhatid::Integer = floor( (1 - beta)*nval )
    omhatx  = temp[ omhatid ]
    temp    = sort( abs( y - median(y) ))
    omhaty  = temp[ omhatid ]
    a       = (x .- pbos(x, beta=beta) )./omhatx
    b       = (y .- pbos(y, beta=beta) )./omhaty
    for i = 1:nval
        if a[i] < -1
            a[i] = -1
        elseif a[i] > 1
            a[i] = 1
        end
        if b[i] < -1
            b[i] = -1
        elseif b[i] > 1
            b[i] = 1
        end
    end
    Pbcor   = sum( a.*b )/sqrt(sum( a.*a ) * sum( b.*b ))
    test    = Pbcor*sqrt( ( nval - 2 )/( 1 - Pbcor*Pbcor ) )
    sig     = 2*( 1 - Rmath.pt(abs(test), nval-2))

    METHOD=nothing

    output = testOutput()
    output.method = METHOD
    output.p = sig
    output.estimate = Pbcor
    output.statistic = test
    output.n = nval
    output.df = nval - 2
    output
end

function outer{S <: Real, T <: Real}(x::Vector{S}, y::Vector{T}, f::Function)
    nx      = length(x)
    ny      = length(y)
    output  = zeros(nx, ny)
    for i = 1:ny
        for j = 1:nx
            output[j, i] = f(x[j], y[i])
        end
    end
    return output
end

function hd{S <: Real}(x::Vector{S}; q::Real=0.5)
    #Compute the Theil-Sen regression estimator.
    # Only a single predictor is allowed in this version
    n    = length(x)
    m1   = ( n + 1 )*q
    m2   = ( n + 1 )*(1 - q)
    vec1 = [1:n]./n
    vec2 = ([1:n] - 1)./n
    w    = Rmath.pbeta( vec1, m1, m2 ) - Rmath.pbeta( vec2,  m1, m2 )
    return sum( w.*sort(x) )
end

function tsp1reg{S <: Real, T <: Real}(x::Vector{S}, y::Vector{T}, HD::Bool)
    order = sortperm( x )
    xsort = x[ order ]
    ysort = y[ order ]
    vec1  = outer( ysort, ysort, - )
    vec2  = outer( xsort, xsort, - )
    v1    = vec1[ vec2 .> 0 ]
    v2    = vec2[ vec2 .> 0 ]
    b1    = median!( v1./v2 )
    b0    = 0.0
    if !HD
        b0 = median( y ) - b1 * median( x )
    else
        b0 = hd( y ) - b1*hd( x )
    end
    return [b0, b1]
end

function tsreg_coef(mf::ModelFrame, HD::Bool, iter::Integer)
    y     = vector( model_response( mf ) )
    x     = ModelMatrix(mf).m[:,2:end]
    np, n = size(x, 2), size(x, 1)

    #ONE PREDICTOR
    if np == 1
        coef = tsp1reg( x[:], y, false )
    #    coef[1], coef[2] =  output.intercept, output.slope
    #    res = output.res
    else
    #MULTIPLE PREDICTORS
        coef_temp = zeros( np )
        for p = 1:np
            coef_temp[p] = tsp1reg( x[:,p], y, false )[2]
        end
        res = y - x*coef_temp
        if !HD
            b0 = median!( res )
        else
            b0 = hd( res )
        end
        #r = zeros( n )
        #coef_old = coef_temp[:]
        for i = 1:iter
            for p = 1:np
                r = y - x*coef_temp - b0 + coef_temp[p].*x[:,p]
                coef_temp[p] = tsp1reg(x[:,p], r, false)[2]
            end
            if !HD
                b0 = median!( y - x*coef_temp )
            else
                b0 = hd( y - x*coef_temp )
            end
            coef_old = coef_temp[:]
        end
        coef = [b0, coef_temp]
    end
    return coef
end

function tsreg(formula::Expr, dataframe::AbstractDataFrame; varfun::Function=pbvar, corfun::Function=pbcor, HD::Bool=false, iter::Integer=10)
#  Compute Theil-Sen regression estimator
#  Gauss-Seidel algorithm is used when there is more than one predictor
    mf = ModelFrame( formula, dataframe )

    #WILL ADD THE ABILITY TO DO MULTIPLE REGRESSION
    output = regOut()
    coef   = zeros(2)
    res    = zeros(n)
    if np == 1
        output = tsp1reg( vector( mf[2] ), mr )
        coef[1], coef[2] =  output.intercept, output.slope
    else
        stop("Only 1 predictor is allowed.")
    end

    res  = temp1.res
    yhat = y - res

    epow   = pbvar(yhat)/pbvar(y)
    if epow >= 1
        epow = sqrt( pbcor(yhat, y).estimate )
    end
    stre   = sqrt(epow)
    output = DataFrame()
    output["b0"] = temp1.intercept
    output["b1"] = temp1.slope
    output["strength of assoc."] = stre
    output["explanatory power"]  = epow

    if np == 1
        lm_coef = coef(lm( formula, dataframe ))
        p     = FramedPlot()
        pts   = Points( vector(mf[2]), mr , "type", "dot")
        s1    = Slope( lm_coef[2], (0, lm_coef[1]), "type", "solid", "color", "blue")
        s2    = Slope( temp1.slope, (0, temp1.intercept), "type", "solid", "color", "red")
        add(p, pts, s1, s2)
        Winston.display(p)
    end
    output
end

#function DataFrame.describe(x::regOut)


text:
using GLMakie
using CSV
using DataFrames

genemetab = CSV.read("data/genemetab.csv", DataFrame)

scatter(log.(1 .+ genemetab.glutdegr), log.(1 .+ genemetab.glutsynth), genemetab.glutgut,
    axis=(; xlabel="degradataion", ylabel="synthesis", zlabel="concentration"))

ext: jl

max_stars_repo_name: Klepac-Ceraj-Lab/ResonanceAnalysis

max_issues_repo_name: Klepac-Ceraj-Lab/ResonanceAnalysis

max_forks_repo_name: Klepac-Ceraj-Lab/ResonanceAnalysis


text:
const MAX_VAL = 6f0
const TABLE_SIZE = 1000

sigmoid_table = zeros(Tsf, TABLE_SIZE)
logsigm_table = zeros(Tsf, TABLE_SIZE)
for i in 1:TABLE_SIZE
	sigmoid_table[i] = exp(((i-1) / TABLE_SIZE * 2 - 1) * MAX_VAL)
	sigmoid_table[i] = sigmoid_table[i] / (1. + sigmoid_table[i])
	logsigm_table[i] = log(sigmoid_table[i] + realmin(Tsf))
end

function index(x::Tsf)
	return round((x + MAX_VAL) * (TABLE_SIZE / MAX_VAL / 2) ) + 1
end

function cached_sigmoid(x::Tsf)
	if x >= MAX_VAL
		return 1f0
	elseif x <= -MAX_VAL
		return 0f0
	end

	idx = index(x)
	if idx > TABLE_SIZE
		return 1f0
	end
	return sigmoid_table[idx]
end

function cached_logsigm(x::Tsf)
	if x >= MAX_VAL
		return 0f0
	elseif x <= -MAX_VAL
		return -6.01f0
	end

	idx = index(x)
	if idx > TABLE_SIZE
		return 0f0
	end
	return logsigm_table[idx]
end


text:
using Catlab
using CombinatorialSpaces
using CombinatorialSpaces.CombinatorialMaps
using Catlab, Catlab.CategoricalAlgebra.CSets
using CombinatorialSpaces.CombinatorialMaps: TheoryCombinatorialMap
using Catlab.Permutations: cycles

export AbstractPlanarGraph, PlanarGraph,
    face, src, dst, boundary, isboundary,
    rem_half_edge!, rem_half_edges!, rem_edge!,
    add_corolla!,
    planar_rz

export incident, σ, α, ϕ, trace_faces, trace_vertices

@present TheoryPlanarGraph <: TheoryCombinatorialMap begin
    V::Ob
    F::Ob

    src::Hom(H,V)   # the source vertex of a half edge
    face::Hom(H,F)  # the face corresponding to a half edge
end

const AbstractPlanarGraph = AbstractACSetType(TheoryPlanarGraph)
const PlanarGraph = CSetType(TheoryPlanarGraph, index=[:ϕ, :src, :face])

face(g::AbstractPlanarGraph, he) = subpart(g, he, :face)
src(g::AbstractPlanarGraph, he) = subpart(g, he, :src)
dst(g::AbstractPlanarGraph, he) = src(g, α(g, he))

"""
    boundary(g::AbstractPlanarGraph)

Returns all half edges on the boundary of a planar graph `g`.
"""
function boundary(g::AbstractPlanarGraph)
    for f in trace_faces(g)
        if face(g, f[1]) == 0
            return f
        end
    end
    error("Boundary not found")
end

isboundary(g::AbstractPlanarGraph, he) = (face(g, he) == 0)

rem_half_edge!(g::AbstractPlanarGraph, he) = rem_part!(g, :H, he)
rem_half_edges!(g::AbstractPlanarGraph, hes) = rem_parts!(g, :H, hes)

function rem_edge!(g::AbstractPlanarGraph, he)

end

function add_corolla!(g::AbstractPlanarGraph, valence::Int)
    v = add_vertex!(g)
    n = nparts(g, :H)
    add_parts!(g, :H, valence; src=v, σ=circshift((n+1):(n+valence), -1))
    return g
end

function pair_half_edges!(g::AbstractPlanarGraph, h1, h2)
    set_subpart!(g, [h1;h2], :α, [h2; h1])
    return g
end

function update_phi!(g::AbstractPlanarGraph)
    hes = collect(1:nparts(g, :H))
    for he in hes
        prev_he = α(g, σ(g, he))
        set_subpart!(g, prev_he, :ϕ, he)
    end
    return g
end

function update_face!(g::AbstractPlanarGraph, he; isboundary = false)
    current_he = he
    next_he = ϕ(g, current_he)
    f = 0
    if !isboundary
        if face(g, he) == 0
            f = add_part!(g, :F)
        else
            f = face(g, he)
        end
    end
    while next_he != he
        set_subpart!(g, current_he, :face, f)
        current_he = next_he
        next_he = ϕ(g, current_he)
    end
    set_subpart!(g, current_he, :face, f)
    return g
end

function planar_rz()
    g = PlanarGraph()
    valences = [3,3,3,3,3,3,3,3,2,4,2,2,2]
    for vl in valences
        add_corolla!(g, vl)
    end
    he_pairs = [(1,26), (2,27), (3,4), (5, 30), (6,7), (8,31), (9,10), (11,33), (12,35),
        (13,36), (14,34), (15,16), (17,32), (18,19), (20,29), (21,22), (23,28), (24,25)]
    pair_half_edges!(g, [p[1] for p in he_pairs], [p[2] for p in he_pairs])
    update_phi!(g)
    update_face!(g, 3; isboundary = true)
    he2f = [2, 5, 8, 11, 20, 1]
    for he in he2f
        update_face!(g, he)
    end

    return g
end


text:
export AbstractAsyncInstrument,
    AbstractSyncInstrument,
    Measurement,
    Counter,
    ObservableCounter,
    Histogram,
    ObservableGauge,
    UpDownCounter,
    ObservableUpDownCounter

"""
    Measurement(value, [attributes=StaticAttrs()])

This is to follow [the specification](https://github.com/open-telemetry/opentelemetry-specification/blob/main/specification/metrics/api.md#measurement)

See also [StaticAttrs](@ref)
"""
struct Measurement{V,T<:StaticAttrs}
    value::V
    attributes::T
end

Measurement(v) = Measurement(v, StaticAttrs())

"""
Sync instrument usually pass its measurement immediately.
"""
abstract type AbstractSyncInstrument{T} <: AbstractInstrument{T} end

"""
Async instrument usually has a callback function (which is named *Observable* in
OpenTelemetry) to get its measurement.

!!! note

    If the return of the callback function is not a [`Measurement`](@ref),
    it will be converted into a `Measurement` with an empty
    [`StaticAttrs`](@ref) implicitly when being uploaded to the associated meter
    provider.
"""
abstract type AbstractAsyncInstrument{T} <: AbstractInstrument{T} end

function examine_instrument(
    ins::AbstractInstrument;
    max_unit_length = 63,
    max_description_length = 1024,
)
    !isnothing(match(r"[a-zA-Z][_0-9a-zA-Z\\.\\-]{0,62}$", ins.name)) ||
        throw(ArgumentError("invalid name: $(ins.name)"))
    length(ins.unit) <= max_unit_length ||
        throw(ArgumentError("length of unit should be no more than $max_unit_length"))
    length(ins.description) <= max_description_length || throw(
        ArgumentError(
            "length of description should be no more than $max_description_length",
        ),
    )
end

(ins::AbstractSyncInstrument)(amount; kw...) =
    ins(Measurement(amount, StaticAttrs(values(kw))))
(ins::AbstractSyncInstrument)(m::Measurement) = push!(ins.meter.provider, ins => m)

(ins::AbstractAsyncInstrument)() = push!(ins.meter.provider, ins => ins.callback())

#####

"""
    Counter{T}(name, meter;unit="", description="") where T

`Counter` is a [`AbstractSyncInstrument`](@ref) which supports non-negative increments.

See more details from [the specification](https://github.com/open-telemetry/opentelemetry-specification/blob/main/specification/metrics/api.md#counter).
"""
struct Counter{T} <: AbstractSyncInstrument{T}
    meter::Meter
    name::String
    unit::String
    description::String
    function Counter{T}(name, meter; unit = "", description = "") where {T}
        c = new{T}(meter, name, unit, description)
        examine_instrument(c)
        push!(meter, c)
        c
    end
end

(c::Counter{T})(; kw...) where {T} = c(one(T); kw...)

function (c::Counter)(m::Measurement)
    if m.value < 0
        # TODO: custom error
        throw(ArgumentError("amount must be non-negative, got $(m.value)"))
    end
    push!(c.meter, c => m)
end

"""
    ObservableCounter{T}(callback, name, meter; unit = "", description = "") where {T}

`ObservableCounter` is an [`AbstractAsyncInstrument`](@ref) which reports monotonically increasing value(s) when the instrument is being observed.

See more details from [the specification](https://github.com/open-telemetry/opentelemetry-specification/blob/main/specification/metrics/api.md#asynchronous-counter)
"""
struct ObservableCounter{T,F} <: AbstractAsyncInstrument{T}
    callback::F
    meter::Meter
    name::String
    unit::String
    description::String
    function ObservableCounter{T}(
        callback::F,
        name,
        meter;
        unit = "",
        description = "",
    ) where {T,F}
        c = new{T,F}(callback, meter, name, unit, description)
        examine_instrument(c)
        push!(meter, c)
        c
    end
end

function (ins::ObservableCounter{T})() where {T}
    v = ins.callback()
    m = v isa Measurement ? v : Measurement(v)
    if m.value >= zero(T)
        push!(ins.meter.provider, ins => m)
    else
        # TODO: custom error
        throw(ArgumentError("amount must be non-negative, got $v"))
    end
end

"""
    Histogram{T}(name, meter; unit = "", description = "") where {T}

`Histogram` is an `AbstractSyncInstrument` which can be used to report arbitrary values that are likely to be statistically meaningful. It is intended for statistics such as histograms, summaries, and percentile.

See more details from [the specification](https://github.com/open-telemetry/opentelemetry-specification/blob/main/specification/metrics/api.md#histogram)
"""
struct Histogram{T} <: AbstractSyncInstrument{T}
    meter::Meter
    name::String
    unit::String
    description::String
    function Histogram{T}(name, meter; unit = "", description = "") where {T}
        h = new{T}(meter, name, unit, description)
        examine_instrument(h)
        push!(meter, h)
        h
    end
end

"""
    ObservableGauge{T}(callback, name, meter; unit = "", description = "",) where {T}

`ObservableGauge` is an [`AbstractAsyncInstrument`](@ref) which reports non-additive value(s) (e.g. the room temperature - it makes no sense to report the temperature value from multiple rooms and sum them up) when the instrument is being observed.

See also the details from [the specification](https://github.com/open-telemetry/opentelemetry-specification/blob/main/specification/metrics/api.md#asynchronous-gauge)
"""
struct ObservableGauge{T,F} <: AbstractAsyncInstrument{T}
    callback::F
    meter::Meter
    name::String
    unit::String
    description::String
    function ObservableGauge{T}(
        callback::F,
        name,
        meter;
        unit = "",
        description = "",
    ) where {T,F}
        g = new{T,F}(callback, meter, name, unit, description)
        examine_instrument(g)
        push!(meter, g)
        g
    end
end

"""
    UpDownCounter{T}(name, meter; unit = "", description = "") where {T}

`UpDownCounter` is a [`AbstractSyncInstrument`](@ref) which supports increments and decrements.

See also the details from [the specification](https://github.com/open-telemetry/opentelemetry-specification/blob/main/specification/metrics/api.md#updowncounter).
"""
struct UpDownCounter{T} <: AbstractSyncInstrument{T}
    meter::Meter
    name::String
    unit::String
    description::String
    function UpDownCounter{T}(name, meter; unit = "", description = "") where {T}
        c = new{T}(meter, name, unit, description)
        examine_instrument(c)
        push!(meter, c)
        c
    end
end

"""
    ObservableUpDownCounter{T}(callback, name, meter; unit = "", description = "") where {T}

`ObservableUpDownCounter` is an [`AbstractAsyncInstrument`](@ref) which reports additive value(s) (e.g. the process heap size - it makes sense to report the heap size from multiple processes and sum them up, so we get the total heap usage) when the instrument is being observed.

See more details from [the specification](https://github.com/open-telemetry/opentelemetry-specification/blob/main/specification/metrics/api.md#asynchronous-updowncounter).
"""
struct ObservableUpDownCounter{T,F} <: AbstractAsyncInstrument{T}
    callback::F
    meter::Meter
    name::String
    unit::String
    description::String
    function ObservableUpDownCounter{T}(
        callback::F,
        name,
        meter;
        unit = "",
        description = "",
    ) where {T,F}
        c = new{T,F}(callback, meter, name, unit, description)
        examine_instrument(c)
        push!(meter, c)
        c
    end
end


text:
@testset "Ulam basis" begin

using ValidatedNumerics
using InvariantMeasures
using InvariantMeasures.BasisDefinition

B = Ulam(4)
@test nonzero_on(B, (@interval(0.1), @interval(0.3))) == (1,2)
@test nonzero_on(B, (@interval(0), @interval(1))) == (1,4)
@test nonzero_on(B, (@interval(0.3), @interval(0.31))) == (2,2)
@test nonzero_on(B, (@interval(1), @interval(1))) == (4,4)

@test is_refinement(Ulam(8), Ulam(4))
@test is_refinement(Ulam(8), Ulam(8))
@test !is_refinement(Ulam(4), Ulam(8))

end


text:
#
# This file is part of the DiscreteEvents.jl Julia package, MIT license
#
# Paul Bayer, 2019
#
# This is a Julia package for discrete event simulation
#

"""
    DiscreteEvents

A Julia package for generating and simulating discrete events. It runs on Julia
`VERSION ≥ v"1.0"`. Multithreading requires `VERSION ≥ v"1.3"`.

The current stable, registered version is installed with
```julia
pkg> add DiscreteEvents
```

The development version is installed with:
```julia
pkg> add("https://github.com/pbayer/DiscreteEvents.jl")
```
"""
module DiscreteEvents

"""
    version

Gives the package version:

```jldoctest
julia> using DiscreteEvents

julia> DiscreteEvents.version
v"0.3.0"
```
"""
const version = v"0.3.0"

using Unitful, Random, DataStructures, Logging, .Threads
import Unitful: FreeUnits, Time

include("types.jl")
include("components.jl")
include("fclosure.jl")
include("schedule.jl")
include("events.jl")
include("clock.jl")
include("process.jl")
include("threads.jl")
include("timer.jl")
include("printing.jl")
include("resources.jl")

export  Clock, PClock, RTClock, RTC, setUnit!, 𝐶,
        Action, Timing, at, after, every, before, until,
        tau, sample_time!, fun, event!, periodic!,
        incr!, run!, stop!, resume!, sync!, resetClock!,
        Prc, process!, interrupt!, delay!, wait!, now!,
        fork!, collapse!, pclock, diagnose, onthread,
        Resource

Random.seed!(123)
rng = MersenneTwister(2020)
𝐶.state == Undefined() ? init!(𝐶) : nothing

end # module


text:
@testset "L-shape" begin
    Ns = 8:8:16
    optim_prec_final = 64

    RR = RealField(1024)
    result = RR("[9.6397238440219410527114592624 +/- 1.43e-26]")
    goalradius = 3e-2

    domain = LShape(RealField(128))
    u = LShapeEigenfunction(domain)
    interval = ball(RR(9.6), RR(0.5))

    println("Computing eigenvalue for the $domain")
    λs, _ = iteratemps(
        domain,
        u,
        interval,
        Ns,
        optim_prec_final = optim_prec_final,
        optim_prec_linear = true,
        show_trace = true,
    )
    rad = Float64(radius(λs[end]))

    @test overlaps(result, λs[end])
    @test rad < goalradius

    if rad < goalradius
        @printf "radius = %e < %e\n" rad goalradius
    else
        @printf "radius = %e >= %e\n" rad goalradius
    end
end


text:
abstract type AbstractAction{T} end
abstract type AbstractInteraction{T} <: AbstractAction{T} end
abstract type AbstractSpaceAction{T,N} <: AbstractAction{T} end
abstract type AbstractSpaceInteraction{T,N} <: AbstractSpaceAction{T,N} end
abstract type AbstractPopulationP{NP,AP} end

struct NullAction{T} <: AbstractAction{T} end
(na::NullAction)(args...) = nothing

struct CompositeAction{T,P} <: AbstractAction{T}
    arr::Vector{P}
end
CompositeAction(arr::Vector{A}) where {T, A<:AbstractAction{T}} = CompositeAction{T,A}(arr)
function (cp::CompositeParameter)(args...)
    actions = [param(args...) for param in cp.arr]
    types = Union{typeof.(actions)...}
    actions = Array{types}(actions)
    CompositeAction(actions)
end
function (ca::CompositeAction)(args...)
    for action in ca.arr
        action(args...)
    end
end


text:
"""
    TEVars(target_future::Vector{Int}
        target_presentpast::Vector{Int}
        source_presentpast::Vector{Int}
        conditioned_presentpast::Vector{Int})

Which axes of the state space correspond to the future of the target,
the present/past of the target, the present/past of the source, and
the present/past of any conditioned variables? Indices correspond to variables
of the embedding or colums of a `Dataset`.

This information is used by the transfer entropy estimators to ensure
the marginal distributions are computed correctly.
"""
struct TEVars
    target_future::Vector{Int}
    target_presentpast::Vector{Int}
    source_presentpast::Vector{Int}
    conditioned_presentpast::Vector{Int}
end

"""
    TEVars(target_future::Vector{Int},
            target_presentpast::Vector{Int},
            source_presentpast::Vector{Int})

Which axes of the state space correspond to the future of the target,
the present/past of the target, and the present/past of the source?
Indices correspond to variables of the embedding or colums of a `Dataset`.

This information is used by the transfer entropy estimators to ensure
the marginal distributions are computed correctly.

This three-argument constructor assumes there will be no conditional variables.
"""
TEVars(Tf::Vector{Int}, Tpp::Vector{Int}, Spp::Vector{Int}) = TEVars(Tf, Tpp, Spp, Int[])

"""
    TEVars(;Tf = Int[], Tpp = Int[], Spp = Int[], Cpp = Int[])

Which axes of the state space correspond to the future of the target (`Tf`),
the present/past of the target (`Tpp`), the present/past of the source (`Spp`), and
the present/past of any conditioned variables (`Cpp`)? Indices correspond to variables
of the embedding or colums of a `Dataset`.

This information is used by the transfer entropy estimators to ensure
the marginal distributions are computed correctly.
"""
TEVars(;Tf::Vector{Int} = Int[],
	    Tpp::Vector{Int} = Int[],
	    Spp::Vector{Int} = Int[],
		Cpp::Vector{Int} = Int[]) =
	TEVars(Tf, Tpp, Spp, Cpp)

export TEVars


text:
using SavGol
using Test
using DelimitedFiles

@testset "SavGol.jl" begin
    # Write your own tests here.
    data  = readdlm("signal.txt")
    sg = SG(20,3)
    smoothed = apply_filter(sg[:,1],data)

end


text:
#Function for obtaining the gravitational acceleration of at a position w.r.t. the central reference frame
function GetGravAccel(sys::PlanetarySystem, Pos::Vector{T}, ep::S) where {T,S <: Real}
    g = zeros(3)
    for body in sys.Bodies
        #Obtain the position of the vehicle position w.r.t. the C.o.M. of the current celestial body
        PosCart = Pos - CelestialBodyPosition(sys.InertialFrame_Body, body, ep)

        #Compute gravitational contribution of the current celestial body in the V-frame w.r.t. the C.o.M. of the current selected body -> Transform to inertial frame
        if (typeof(body) <: abstractStar) #in the case of a star, always use central gravity
            glocalv = CentralGravity(body, PosCart)
        else
            glocalv = body.Environment.GravityModel(body, PosCart)
        end
        #Compute the rotation quaternion of the V-frame of the current celestial body to the inertial frame (J2000)
        #Because they are vectors, no translation is needed, as the inertial frame of each body is equal to the global inertial frame (in orientation)
        Rsph,lon,lat = PositionCart2Sph(PosCart)
        qv2r = Q_V2R(lon,lat)
        gbodyr = RotateVector(qv2r, glocalv)
        g += BodyFixed2Inertial3(body,ep)*gbodyr
    end
    return g
end

#Clearly a need to transform a position vector w.r.t. a body to the inertial frame.
function GetInertialFramePos(sys::PlanetarySystem, Body::T, Pos::Vector{S}, ep::P) where {T <: abstractCelestialBody, S,P <: Real}
    return (Pos + CelestialBodyPosition(sys.InertialFrame_Body, Body, ep))
end


text:
# Bogumił Kamiński, 2022

# Codes for chapter 10

# Code for section 10.1

aq = [10.0   8.04  10.0  9.14  10.0   7.46   8.0   6.58
       8.0   6.95   8.0  8.14   8.0   6.77   8.0   5.76
      13.0   7.58  13.0  8.74  13.0  12.74   8.0   7.71
       9.0   8.81   9.0  8.77   9.0   7.11   8.0   8.84
      11.0   8.33  11.0  9.26  11.0   7.81   8.0   8.47
      14.0   9.96  14.0  8.1   14.0   8.84   8.0   7.04
       6.0   7.24   6.0  6.13   6.0   6.08   8.0   5.25
       4.0   4.26   4.0  3.1    4.0   5.39  19.0  12.50
      12.0  10.84  12.0  9.13  12.0   8.15   8.0   5.56
       7.0   4.82   7.0  7.26   7.0   6.42   8.0   7.91
       5.0   5.68   5.0  4.74   5.0   5.73   8.0   6.89];

using DataFrames

# Code for listing 10.1

aq1 = DataFrame(aq, ["x1", "y1", "x2", "y2", "x3", "y3", "x4", "y4"])
DataFrame(aq, [:x1, :y1, :x2, :y2, :x3, :y3, :x4, :y4])

# Code for creating DataFrame with automatic column names

DataFrame(aq, :auto)

# Codes for creating DataFrame from vector of vectors

aq_vec = collect(eachcol(aq))
DataFrame(aq_vec, ["x1", "y1", "x2", "y2", "x3", "y3", "x4", "y4"])
DataFrame(aq_vec, :auto)

# Codes for section 10.1.2

data = (set1=(x=aq[:, 1], y=aq[:, 2]),
        set2=(x=aq[:, 3], y=aq[:, 4]),
        set3=(x=aq[:, 5], y=aq[:, 6]),
        set4=(x=aq[:, 7], y=aq[:, 8]));

data.set1.x

DataFrame(x1=data.set1.x, y1=data.set1.y,
          x2=data.set2.x, y2=data.set2.y,
          x3=data.set3.x, y3=data.set3.y,
          x4=data.set4.x, y4=data.set4.y)

DataFrame(:x1 => data.set1.x, :y1 => data.set1.y,
          :x2 => data.set2.x, :y2 => data.set2.y,
          :x3 => data.set3.x, :y3 => data.set3.y,
          :x4 => data.set4.x, :y4 => data.set4.y)

DataFrame([:x1 => data.set1.x, :y1 => data.set1.y,
           :x2 => data.set2.x, :y2 => data.set2.y,
           :x3 => data.set3.x, :y3 => data.set3.y,
           :x4 => data.set4.x, :y4 => data.set4.y]);

[(i, v) for i in 1:4 for v in [:x, :y]]

[string(v, i) for i in 1:4 for v in [:x, :y]]

[string(v, i) => getproperty(data[i], v)
        for i in 1:4 for v in [:x, :y]]

DataFrame([string(v, i) => getproperty(data[i], v)
           for i in 1:4 for v in [:x, :y]]);

data_dict = Dict([string(v, i) => getproperty(data[i], v)
                         for i in 1:4 for v in [:x, :y]])
collect(data_dict)

DataFrame(data_dict)

df1 = DataFrame(x1=data.set1.x)
df1.x1 === data.set1.x

df2 = DataFrame(x1=data.set1.x; copycols=false)
df2.x1 === data.set1.x

df = DataFrame(x=1:3, y=1)
df.x

DataFrame(x=[1], y=[1, 2, 3])

using RCall
r_df = R"data.frame(a=1:6, b=1:2, c=1:3)"
julia_df = rcopy(r_df)

# Codes for section 10.1.3

data.set1
DataFrame(data.set1)

DataFrame([(a=1, b=2), (a=3, b=4), (a=5, b=6)])

data

# Code for listing 10.2

aq2 = DataFrame(data)

# Codes for listing 10.3

data_dfs = map(DataFrame, data)

# Codes for vertical concatenation examples

vcat(data_dfs.set1, data_dfs.set2, data_dfs.set3, data_dfs.set4)

vcat(data_dfs.set1, data_dfs.set2, data_dfs.set3, data_dfs.set4;
     source="source_id")

vcat(data_dfs.set1, data_dfs.set2, data_dfs.set3, data_dfs.set4;
     source="source_id"=>string.("set", 1:4))

reduce(vcat, collect(data_dfs);
       source="source_id"=>string.("set", 1:4))

# Code for listing 10.4

df1 = DataFrame(a=1:3, b=11:13)
df2 = DataFrame(a=4:6, c=24:26)
vcat(df1, df2)
vcat(df1, df2; cols=:union)

# Code for listing 10.5

df_agg = DataFrame()
append!(df_agg, data_dfs.set1)
append!(df_agg, data_dfs.set2)

# Code for appending tables to a data frame

df_agg = DataFrame()
append!(df_agg, data.set1)
append!(df_agg, data.set2)

# Code for promote keyword argument

df1 = DataFrame(a=1:3, b=11:13)
df2 = DataFrame(a=4:6, b=[14, missing, 16])
append!(df1, df2)
append!(df1, df2; promote=true)

# Code for section 10.2.3

df = DataFrame()
push!(df, (a=1, b=2))
push!(df, (a=3, b=4))

df = DataFrame(a=Int[], b=Int[])
push!(df, [1, 2])
push!(df, [3, 4])

function sim_step(current)
    dx, dy = rand(((1,0), (-1,0), (0,1), (0,-1)))
    return (x=current.x + dx, y=current.y + dy)
end

using BenchmarkTools
@btime rand(((1,0), (-1,0), (0,1), (0,-1)));

dx, dy = (10, 20)
dx
dy

using FreqTables
using Random
Random.seed!(1234);
proptable([rand(((1,0), (-1,0), (0,1), (0,-1))) for _ in 1:10^7])

using Random
Random.seed!(6);
walk = DataFrame(x=0, y=0)
for _ in 1:10
    current = walk[end, :]
    push!(walk, sim_step(current))
end
walk

using Plots
plot(walk.x, walk.y;
     legend=false,
     series_annotations=1:11,
     xticks=range(extrema(walk.x)...),
     yticks=range(extrema(walk.y)...))

extrema(walk.y)

range(1, 5)

(3/4)^9

# Code for listing 10.6

function walk_unique() #A
    walk = DataFrame(x=0, y=0)
    for _ in 1:10
        current = walk[end, :]
        push!(walk, sim_step(current))
    end
    return nrow(unique(walk)) == nrow(walk) #B
end
Random.seed!(2);
proptable([walk_unique() for _ in 1:10^5])

# code for serialization

using Serialization
serialize("walk.bin", walk)
deserialize("walk.bin") == walk

max_stars_repo_name: Mo-Gul/JuliaForDataAnalysis

max_issues_repo_name: Mo-Gul/JuliaForDataAnalysis

max_forks_repo_name: Mo-Gul/JuliaForDataAnalysis


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# Run the distributed test outside of the main driver since it needs its own
# set of dedicated workers.

include("testenv.jl")
cmd = `$test_exename $test_exeflags distributed_exec.jl`

if !success(pipeline(cmd; stdout=STDOUT, stderr=STDERR)) && ccall(:jl_running_on_valgrind,Cint,()) == 0
    error("Distributed test failed, cmd : $cmd")
end

max_stars_repo_name: caodi/julia

max_issues_repo_name: caodi/julia

max_forks_repo_name: caodi/julia


text:
#!/usr/bin/julia

# Trizen
# 28 August 2020
# https://github.com/trizen

# Generate all the positive divisors of n.

using Primes

function divisors(n)

    d = Int64[1]

    for (p,e) in factor(n)
        t = Int64[]
        r = 1

        for i in 1:e
            r *= p
            for u in d
                push!(t, u*r)
            end
        end

        append!(d, t)
    end

    return sort(d)
end

println(divisors(5040))

max_stars_repo_name: trizen/julia-scripts

max_issues_repo_name: trizen/julia-scripts

max_forks_repo_name: trizen/julia-scripts


text:
function error_estimate(LML)
    lm = value(LML)
    p_posterior = exp.(LML.p_posterior .- lm)
    p_proposal = exp.(LML.p_proposal .- lm)
    g_posterior = exp.(LML.g_posterior)
    g_proposal = exp.(LML.g_proposal)

    n1 = length(p_posterior)
    n2 = length(g_proposal)
    s1 = n1 / (n1+n2)
    s2 = n2 / (n1+n2)
    f1 = p_proposal ./ (s1 * p_proposal + s2 * g_proposal)
    f2 = g_posterior ./ (s1 * p_posterior + s2 * g_posterior)

    ρ_f2 = spectrum0(f2)

    re2 = 1/n2 * var( f1 ) / mean( f1 )^2
    re2 += ρ_f2/n1 * var( f2 ) / mean( f2 )^2

    return (rmse = re2, cv = sqrt(re2), percent = sqrt(re2)*100)
end

## Spectrum of the autoregressive model at frequency 0

function spectrum0(x)
    σ², B, order = fit_auto_regressive(x)
    return σ² / (1-sum(B))^2
end

function fit_auto_regressive(x)
    n = length(x)
    maxorder = Int(floor(10*log(10, n)))
    coeffs, varpred = levinson_fit(x, maxorder)
    order = argmin(vec(n * (log.(varpred') .+ 1) .+ 2 * ((1:maxorder) .+ 1)))

    b = -coeffs[order,1:order]
    σ² = varpred[order]
    σ² *= n / (n - order - 1)
    return σ², b, order
end

function levinson_fit(x::AbstractVector, p::Int)
    # Copyright 2012-2021 DSP.jl contributors : https://github.com/JuliaDSP/DSP.jl/blob/master/src/lpc.jl
    R_xx = autocov(x)
    a = zeros(p,p)
    prediction_err = zeros(1,p)

    # for m = 1
    a[1,1] = -R_xx[2]/R_xx[1]
    prediction_err[1] = R_xx[1]*(1-a[1,1]^2)

    # for m = 2,3,4,..p
    for m = 2:p
        a[m,m] = (-(R_xx[m+1] + dot(a[m-1,1:m-1],R_xx[m:-1:2]))/prediction_err[m-1])
        a[m,1:m-1] = a[m-1,1:m-1] + a[m,m] * a[m-1,m-1:-1:1]
        prediction_err[m] = prediction_err[m-1]*(1-a[m,m]^2)
    end

    # Return autocorrelation coefficients and error estimate
    a, prediction_err
end


text:
# extracts vocabularies and emojis per region
# run as:
# JULIA_PROJECT=. julia -p64 -Lsrc/vocabulary-and-emojis.jl
# and then main(lang)

using CSV, Glob, DataFrames, JLD2, CodecZlib, TextSearch, InvertedFiles
using InvertedFiles: add!

include("io.jl")

function compute_vocabulary(filename, mindocs=5)
    voc = Dict{String,Int}()
    sizehint!(voc, 100_000)
    @info "computing vocabulary for $filename"

    n = 0
    doc = Set{String}()
    for t in eachline(filename)
        n += 1
        empty!(doc)

        for w in split(t)
            if !(w in doc)
                voc[w] = get(voc, w, 0) + 1
                push!(doc, w)
            end
        end
    end

    m = length(voc)

    for (k, freq) in voc
        if freq < mindocs
            delete!(voc, k)
        end
    end

    @info "computed vocabulary for $filename, n=$n, m=$m, final voc: $(length(voc))"
    voc, m, n, filename
end

function vocabulary_table(lang)
    vocfile = "data/$lang/voc.tsv.gz"
    vocstatsfile = "data/$lang/voc-stats.tsv.gz"
    if isfile(vocfile) && isfile(vocstatsfile)
        @info "vocabulary files already exists"
        return
    end
    VOC = DataFrame(token=String[], country_code=String[], ndocs=Int[], idf=Float64[])
    STATS = DataFrame(country_code=String[], n=Int[], rawvoc=Int[], voc=Int[])
    ALL = Dict{String,Int}()
    N = 0

    L = pmap(compute_vocabulary, glob("data/$lang/messages/*.txt"))

    for (voc, m, n, filename) in L
        cc = basename(filename)[1:2]
        N += n
        push!(STATS, (cc, n, m, length(voc)))
        for (token, ndocs) in voc
            prob = (ndocs + 1) / n
            push!(VOC, (token, cc, ndocs, log(1 / prob)))
            ALL[token] = get(ALL, token, 0) + ndocs
        end
    end

    for (token, ndocs) in ALL
        prob = (ndocs + 1) / N
        push!(VOC, (token, "ALL", ndocs, log(1 / prob)))
    end

    M = length(ALL)
    push!(STATS, ("ALL", N, M, M))
    sort!(VOC, :token)

    @info "saving voc files"
    open(vocfile, "w") do f
        gz = GzipCompressorStream(f)
        CSV.write(gz, VOC, delim='\t')
        close(gz)
    end

    open(vocstatsfile, "w") do f
        gz = GzipCompressorStream(f)
        CSV.write(gz, STATS, delim='\t')
        close(gz)
    end
end

function count_emojis(filename)
    voc = load(filename, "voc")

    D = Dict{Char,Int}()
    for (token, freq) in voc
        # there are some cases when the tokenizer considers emojis as part of a word, e.g.,
        # when they start the word
        lst = emojis(token)
        if length(lst) > 0
            for e in lst
                D[e] = get(D, e, 0) + row.ndocs
            end
        end
    end

    replace(basename(filename), "voc-" => "", ".jld2" => "") => D
end

function emoji_table(lang)
    outname = "data/$lang/emojis.tsv.gz"
    isfile(outname) && return

    V = CSV.read("data/$lang/voc.tsv.gz", DataFrame, delim='\t')
    E = DataFrame(emoji=Char[], country_code=[], freq=Int[])

    D = Dict()

    for row in eachrow(V)
        cc = row.country_code

        for e in emojis(row.token)
            D[(e, cc)] = get(D, (e, cc), 0) + row.ndocs
        end
    end

    for (k, freq) in D
        e, cc = k
        push!(E, (e, cc, freq))
    end

    sort!(E, :emoji)

    open(outname, "w") do f
        gz = GzipCompressorStream(f)
        CSV.write(gz, E, delim='\t')
        close(gz)
    end
end

function normalize_text_from_csv(filename)
    config = TextConfig(lc=true, del_diac=true, del_dup=true, del_punc=false, group_num=true)
    tok = Tokenizer(config)
    @info filename
    D = CSV.read(filename, DataFrame, delim='\t')
    T = Vector{String}(undef, length(D.text))

    for (i, t) in enumerate(D.text)
        tokens = tokenize(tok, t)
        T[i] = join(decode.(tok, tokens), ' ')
    end

    T
end

function extract_text(ccpath)
    outname = rstrip(ccpath, '/') * ".txt"
    isfile(outname) && return
    L = pmap(normalize_text_from_csv, glob(joinpath(ccpath, "*.tsv.gz")))

    open(outname * ".tmp", "w") do f
        for lst in L
            for t in lst
                println(f, t)
            end
        end
    end

    mv(outname * ".tmp", outname; force=true)
end

function extract_id(ccpath)
    outname = rstrip(ccpath, '/') * "-id-list.txt.gz"
    isfile(outname) && return

    open(outname * ".tmp", "w") do f
        gz = GzipCompressorStream(f)
        for filename in glob(joinpath(ccpath, "*.tsv.gz"))
            @info "extracting identifiers from $filename"
            D = CSV.read(filename, DataFrame, delim='\t')
            for id in D.id
                println(gz, id)
            end
        end

        close(gz)
    end

    mv(outname * ".tmp", outname; force=true)
end

function main(lang)

    for ccpath in glob("data/$lang/messages/*/")
        extract_text(ccpath)
    end

    @sync @distributed for ccpath in glob("data/$lang/messages/*/")
        extract_id(ccpath)
    end

    vocabulary_table(lang)
    emoji_table(lang)
end

max_stars_repo_name: sadit/RegionalSpanish

max_issues_repo_name: sadit/RegionalSpanish

max_forks_repo_name: sadit/RegionalSpanish


text:
struct ContextPayloadLogger <: AbstractLogger
    logger::AbstractLogger
    ctxvars::Any
end

function _get_task_ctxvars()
    logger = Logging.current_logger()
    if logger isa ContextPayloadLogger
        return logger.ctxvars
    end
    return nothing
end

function with_task_ctxvars(f, ctx)
    @nospecialize
    return Logging.with_logger(f, ContextPayloadLogger(current_logger(), ctx))
end

# Forward actual logging interface:
Logging.handle_message(payload::ContextPayloadLogger, args...; kwargs...) =
    Logging.handle_message(payload.logger, args...; kwargs...)
Logging.shouldlog(payload::ContextPayloadLogger, args...) =
    Logging.shouldlog(payload.logger, args...)
Logging.min_enabled_level(payload::ContextPayloadLogger, args...) =
    Logging.min_enabled_level(payload.logger, args...)
Logging.catch_exceptions(payload::ContextPayloadLogger, args...) =
    Logging.catch_exceptions(payload.logger, args...)

"""
    ContextVariablesX.with_logger(f, logger::AbstractLogger)

Like `Logging.with_logger` but properly propagate the context variables.
"""
function with_logger(f, logger::AbstractLogger)
    @nospecialize
    cpl = Logging.current_logger()
    if cpl isa ContextPayloadLogger
        ctx = cpl.ctxvars
    else
        ctx = nothing
    end
    return Logging.with_logger(f, ContextPayloadLogger(logger, ctx))
end

"""
    ContextVariablesX.current_logger() -> logger::AbstractLogger

Like `Logging.current_logger` but unwraps `ContextPayloadLogger`.
"""
function current_logger()
    logger = Logging.current_logger()
    if logger isa ContextPayloadLogger
        return logger.logger
    else
        return logger
    end
end


text:
using SparseArrays
using LinearAlgebra
using Plots

function vekit(A, y, e)
    l=[0];
    while true
        # Berechnung gemäß Algorithmus (Arbeitsblatt)
        y_h = A*y
        l_h = dot(y_h,y)
        # neue Nährung Eigenwert zu Liste Nährungen Eigenwert hinzufügen
        l=[l;l_h]
        y=y_h/norm(y_h)
        # Abbruchkriterium
        if (norm(A*y-l_h*y)<e)
            # Liste Nährungen Eigenwert und Eigenvektor zurückgeben
            return l,y;
        end
    end
end

max_stars_repo_name: SoulKathi/vekit

max_issues_repo_name: SoulKathi/vekit

max_forks_repo_name: SoulKathi/vekit


text:
function loadpuzzle(s::String)
    s = replace(s, r"\s" => "")
    s = map(c -> c in "-._" ? 0 : parse(Int, c), s |> collect)
    reshape(s, 9, 9) |> permutedims
end
loadpuzzle(mat::Matrix{<:Integer}) = copy(mat)
function possinit!(P::Matrix{UInt16}, grid::Matrix)
    P .= 0
    for r in 1:9
        for c in 1:9
            p = (r - 1) ÷ 3 * 3 + (c - 1) ÷ 3 + 1
            v = grid[r, c]
            @assert 0 <= v <= 9
            if v > 0
                bm = 1 << (v-1)
                P[1, r] |= bm
                P[2, c] |= bm
                P[3, p] |= bm
            end
        end
    end
end
function intlog2(v)
    if v == 1 << 4
        return 4
    elseif v < 1 << 4
        if v == 1 << 2
            return 2
        elseif v > 1 << 2
            return 3
        else
            return v - 1
        end
    else
        if v == 1 << 7
            return 7
        elseif v > 1 << 7
            return 8
        elseif v < 1 << 6
            return 5
        else
            return 6
        end
    end
end
# 朴素深搜（顺序搜索，逐列）
function sudokudfs!(grid::Matrix, P::Matrix{UInt16}=zeros(UInt16, 3, 9), id=1)
    if id == 82
        return true
    end
    r = (id - 1) % 9 + 1
    c = (id - 1) ÷ 9 + 1
    p = (r - 1) ÷ 3 * 3 + (c - 1) ÷ 3 + 1
    if grid[r, c] == 0
        @inbounds mask = P[1, r] | P[2, c] | P[3, p]
        while mask < 1 << 9 - 1
            mask2 = mask | (mask+1)
            bm = mask2 - mask
            @inbounds P[1, r] |= bm
            @inbounds P[2, c] |= bm
            @inbounds P[3, p] |= bm
            if sudokudfs!(grid, P, id + 1)
                grid[r, c] = intlog2(bm) + 1
                return true
            end
            @inbounds P[1, r] &= ~bm
            @inbounds P[2, c] &= ~bm
            @inbounds P[3, p] &= ~bm
            mask = mask2
        end
        return false
    else
        return sudokudfs!(grid, P, id + 1)
    end
end
function naivesolver!(grid, P::Matrix{UInt16}=zeros(UInt16, 3, 9); check=true)
    possinit!(P, grid)
    valid = sudokudfs!(grid, P)
    return check ? (valid ? grid : nothing) : grid
end

# 重排搜索顺序
function sudokudfs!(grid::Matrix{<:Integer}, P::Matrix{UInt16}, order::Vector, id=1)
    if id == 82
        return true
    end
    id2 = order[id]
    r = (id2 - 1) % 9 + 1
    c = (id2 - 1) ÷ 9 + 1
    p = (r - 1) ÷ 3 * 3 + (c - 1) ÷ 3 + 1
    if grid[r, c] == 0
        @inbounds mask = P[1, r] | P[2, c] | P[3, p]
        while mask < 1 << 9 - 1
            mask2 = mask | (mask+1)
            bm = mask2 - mask
            @inbounds P[1, r] |= bm
            @inbounds P[2, c] |= bm
            @inbounds P[3, p] |= bm
            if sudokudfs!(grid, P, id + 1)
                grid[r, c] = intlog2(bm) + 1
                return true
            end
            @inbounds P[1, r] &= ~bm
            @inbounds P[2, c] &= ~bm
            @inbounds P[3, p] &= ~bm
            mask = mask2
        end
        return false
    else
        return sudokudfs!(grid, P, id + 1)
    end

end
function presort()
    ids = [[1:9...], [9:-1:1...],[1:9...], [9:-1:1...],
    [1:9...], [9:-1:1...],[1:9...], [9:-1:1...],[1:9...]]# 一条龙序
    ids = vcat(broadcast((a, b) -> a .+ b, ids, 0:9:72)...)
end
const ORDER = presort()
function reordersolver!(grid, P::Matrix{UInt16}=zeros(UInt16, 3, 9), order=ORDER; check=true)
    possinit!(P, grid)
    valid = sudokudfs!(grid, P, order)
    return check ? (valid ? grid : nothing) : grid
end

# 优先深搜
function possinit!(P::AbstractArray{T,3}, b::AbstractArray{W,2}) where {T,W}
    P .= 0
    for r in 1:9
        for c in 1:9
            if b[r, c] > 0
                hold!(P, r, c, b[r, c])
                P[r, c, :] .+= 1
                P[r, c, b[r, c]] -= 1
            end
        end
    end
end
function hold!(P::AbstractArray{T,3}, r, c, v) where T
    for i in 1:c - 1
        P[r, i, v] += 1
    end
    for i in c + 1:9
        P[r, i, v] += 1
    end
    for i in 1:r - 1
        P[i, c, v] += 1
    end
    for i in r + 1:9
        P[i, c, v] += 1
    end
    rr = (r - 1) % 3 + 1
    rc = (c - 1) % 3 + 1
    s = ((1, 2), (-1, 1), (-2, -1))
    for i in s[rr]
        for j in s[rc]
            P[r + i, c + j, v] += 1
        end
    end
end

function release!(P::AbstractArray{T,3}, r, c, v) where T
    for i in 1:c - 1
        P[r, i, v] -= 1
    end
    for i in c + 1:9
        P[r, i, v] -= 1
    end
    for i in 1:r - 1
        P[i, c, v] -= 1
    end
    for i in r + 1:9
        P[i, c, v] -= 1
    end
    rr = (r - 1) % 3 + 1
    rc = (c - 1) % 3 + 1
    s = ((1, 2), (-1, 1), (-2, -1))
    for i in s[rr]
        for j in s[rc]
            P[r + i, c + j, v] -= 1
        end
    end
end
mutable struct KVpair{KT,VT}
    key::KT
    value::VT
end
Base.first(p::KVpair) = p.key
Base.last(p::KVpair) = p.value

function heapinit!(P::AbstractArray{T,3}) where T
    pocount = sum(P .== 0, dims=3)[:]
    heap = tracHeap([KVpair(v, id) for (id, v) in enumerate(pocount)])
    buildminheap!(heap)
    return heap
end
function increasehold(P::AbstractArray{T,3}, H::tracHeap, r, c, v) where T
    if P[r, c, v] == 0
        id = (c - 1) * 9 + r
        lc = H.locL[id]
#         @assert last(H[lc])==id
        if lc <= length(H)
            oldK = first(H[lc])
            H[lc].key -= 1
            heapchangekey!(H, lc, H[lc], oldK)
        end

    end
P[r, c, v] += 1
end

function decreasehold(P::AbstractArray{T,3}, H::tracHeap, r, c, v) where T # for release!
    if P[r, c, v] == 1
        id = (c - 1) * 9 + r
        lc = H.locL[id]
#         @assert last(H[lc])==id
        if lc <= length(H)
#             oldK = first(H[lc])
            H[lc].key += 1
#             heapchangekey!(H, lc, H[lc]) #交给undo!去完成结构调整
        end
    end
P[r, c, v] -= 1
end

function hold!(P::AbstractArray{T,3}, r, c, v, H::tracHeap) where T
    for i in 1:c - 1
        increasehold(P, H, r, i, v)
    end
    for i in c + 1:9
        increasehold(P, H, r, i, v)
    end
    for i in 1:r - 1
        increasehold(P, H, i, c, v)
    end
    for i in r + 1:9
        increasehold(P, H, i, c, v)
    end
    rr = (r - 1) % 3 + 1
    rc = (c - 1) % 3 + 1
    s = ((1, 2), (-1, 1), (-2, -1))
    for i in s[rr]
        for j in s[rc]
            increasehold(P, H, r + i, c + j, v)
        end
    end
end

function release!(P::AbstractArray{T,3}, r, c, v, H::tracHeap) where T
    for i in 1:c - 1
        decreasehold(P, H, r, i, v)
    end
    for i in c + 1:9
        decreasehold(P, H, r, i, v)
    end
    for i in 1:r - 1
        decreasehold(P, H, i, c, v)
    end
    for i in r + 1:9
        decreasehold(P, H, i, c, v)
    end
    rr = (r - 1) % 3 + 1
    rc = (c - 1) % 3 + 1
    s = ((1, 2), (-1, 1), (-2, -1))
    for i in s[rr]
        for j in s[rc]
            decreasehold(P, H, r + i, c + j, v)
        end
    end
end
function sudokudfs!(P::AbstractArray{T,3}, H::tracHeap, id=last(heapextractmin!(H))) where T
    r = (id - 1) % 9 + 1
    c = (id - 1) ÷ 9 + 1
    if length(H) == 0
        for v in 1:9
            if P[r, c, v] == 0
                return true
            end
        end
        return false
    end
    for v in 1:9
        if P[r, c, v] == 0
            snapshot1 = historylen(H)
            hold!(P, r, c, v, H)
            snapshot2 = historylen(H)
            id = last(heapextractmin!(H))
            if sudokudfs!(P, H, id)
#                 @show id
                return true
            end
            undo!(H, tostep=snapshot2)
            release!(P, r, c, v, H)
            undo!(H, tostep=snapshot1)
        end
    end
    return false
end
function prioritysolver!(b, P::AbstractArray{T,3}=zeros(Int, 9, 9, 9); check=true) where T
    possinit!(P, b)
    H = heapinit!(P)
    valid = sudokudfs!(P, H)
    ans = reshape(mapslices(argmin, P, dims=3), 9, 9)
    return check ? (valid ? ans : nothing) : ans
end
function solvesudoku!(sudoku::Matrix{<:Integer}, args...; solver=naivesolver!, kargs...)
    solver(sudoku, args...; kargs...)
end
solvesudoku(sudoku, args...; kargs...) = solvesudoku!(loadpuzzle(sudoku), args...; kargs...)


text:
"""
    tabTruncMln_serial

Table containing the truncated response matrices from size `3*1` to `3*N`.
"""
tabTruncMln_serial = Vector{Array{Float64,2}}([zeros(Float64,3*k,3*k) for k=1:N+1])

"""
    tabEigValsMln_serial

Table containing the eigenvalues of the truncated response matrices from size `3*1` to `3*N`.
"""
tabEigValsMln_serial = Vector{Vector{Complex{Float64}}}([zeros(Float64,3*k) for k=1:N+1])

"""
    tabTruncMln!([args])

Fills the table `tabTruncMln` with the truncated matrices.
"""
function tabTruncMln!(tabTruncMln=tabTruncMln_serial, tabAln=tabAln_serial,
    tabBln=tabBln_serial, tabCln=tabCln_serial, tabDln=tabDln_serial,
    tabFln=tabFln_serial, tabGln=tabGln_serial, tabHln=tabHln_serial)
    for k=1:N+1
        for l=1:k
            for n=1:k
                tabTruncMln[k][l,n] = tabAln[l,n]
                tabTruncMln[k][l,k+n] = tabBln[l,n]
                tabTruncMln[k][l,2*k+n] = tabCln[l,n]
                tabTruncMln[k][k+l,n] = tabDln[l,n]
                tabTruncMln[k][k+l,k+n] = tabAln[l,n]
                tabTruncMln[k][k+l,2*k+n] = tabFln[l,n]
                tabTruncMln[k][2*k+l,n] = tabGln[l,n]
                tabTruncMln[k][2*k+l,k+n] = tabHln[l,n]
                tabTruncMln[k][2*k+l,2*k+n] = tabAln[l,n]
            end
        end
    end
end

"""
    tabEigValsMln!([args])

Fills the table `tabTruncMln` with the truncated matrices.
"""
function tabEigValsMln!(tabTruncMln=tabTruncMln_serial,
    tabEigValsMln=tabEigValsMln_serial)
    for k=1:N+1
        tabEigValsMln[k] = eigvals(tabTruncMln[k])
    end
end

"""
    getPhysicalEigvals([args])

Computes the physical eigenvalues of the response matrix.
"""
function getPhysicalEigvals(tabEigValsMln=tabEigValsMln_serial)
    liEgv = zeros(Complex{Float64},3*(N+1))
    liEigvals = tabEigValsMln[N+1]
    liEgvCvCurrent = zeros(Complex{Float64},N)

    for iegv=1:3*(N+1) # take one of the eigenvalues of Mln(N+1)
        egv = liEigvals[iegv]
        for k=1:N # compare to those of its truncatures
            dist_min = -1
            egvSelect = 0
            for kegv=1:3*k # at a given truncature, take the closest one
                egvTrunc = tabEigValsMln[k][kegv]
                dist = abs(egv-egvTrunc)
                if ((dist_min<0) || (dist<dist_min))
                    egvSelect = egvTrunc
                    dist_min = dist
                end
            end
            if ((dist_min>=0) && (dist_min<error_dis)) # check if the closest one is within the accepted error margin
                liEgvCvCurrent[k] = egvSelect
            else # if not, put Inf to discard
                liEgvCvCurrent[k] = Inf
            end
        end

        nb_egv = 1 #(1 comes from initial separate egv at N+1)
        for k=1:N # count number of accepted truncature eigenvalues
            if (liEgvCvCurrent[k] != Inf)
                nb_egv += 1
            end
        end

        moy_egv = 0
        var_egv = 0

        for k=1:N # compute mean
            if (liEgvCvCurrent[k] != Inf)
                moy_egv += liEgvCvCurrent[k]
            end
        end
        moy_egv = (moy_egv+egv)/nb_egv

        for k=1:N # compute variance
            if (liEgvCvCurrent[k] != Inf)
                var_egv += abs2(liEgvCvCurrent[k]-moy_egv)
            end
        end
        var_egv = (var_egv+abs2(egv-moy_egv))/nb_egv

        if ((var_egv < error_var) && (nb_egv >= threshold_nb_egv)) # recover mean physical eigenvalue
            liEgv[iegv] = egv
        else
            liEgv[iegv] = Inf
        end

        for k=1:N
            liEgvCvCurrent[k] = 0.0
        end

    end

    nb_phys_egv = 0
    for iegv=1:3*(N+1)
        if (liEgv[iegv] != Inf)
            nb_phys_egv += 1
        end
    end

    liPhysicalEgv = zeros(Complex{Float64},nb_phys_egv)
    index = 1
    for iegv=1:3*(N+1)
        if (liEgv[iegv] != Inf)
            liPhysicalEgv[index] = liEgv[iegv]
            index += 1
        end
    end

    return liPhysicalEgv
end

"""
    get_max(li)

Returns the tuple `(kEgv, egv)` where `egv` is the element of the list `li` with the highest
imaginary part and `kEgv` its position on the list.
"""
function get_max(li)
    len = length(li)
    max = -Inf
    egv = -Inf
    kEgv = 0
    for k=1:len
        if (imag(li[k])>max)
            egv = li[k]
            max = imag(egv)
            kEgv = k
        end
    end
    return kEgv, egv
end

######################################
# Clear table
######################################

"""
    tabTruncMln_clear!([args])

Clears the table containing the truncated response matrices.
"""
function tabTruncMln_clear!(tabTruncMln=tabTruncMln_serial)
    for k=1:N+1
        for l=1:k
            for n=1:k
                tabTruncMln[k][l,n] = 0.0
                tabTruncMln[k][l,k+n] = 0.0
                tabTruncMln[k][l,2*k+n] = 0.0
                tabTruncMln[k][k+l,n] = 0.0
                tabTruncMln[k][k+l,k+n] = 0.0
                tabTruncMln[k][k+l,2*k+n] = 0.0
                tabTruncMln[k][2*k+l,n] = 0.0
                tabTruncMln[k][2*k+l,k+n] = 0.0
                tabTruncMln[k][2*k+l,2*k+n] = 0.0
            end
        end
    end
end

"""
    tabEigValsMln_clear!([args])

Clears the table containing the eigenvalues.
"""
function tabEigValsMln_clear!(tabEigValsMln=tabEigValsMln_serial)
    for k=1:N+1
        tabEigValsMln[k] = zeros(Float64,3*k)
    end
end

max_stars_repo_name: KerwannTEP/LiRA

max_issues_repo_name: KerwannTEP/LiRA

max_forks_repo_name: KerwannTEP/LiRA


text:
#	cylinders and toroidal surface by Julia's implementation of PLaSM operators

using ViewerGL; GL = ViewerGL
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation

#	1D complex in 2D generated by Lar and PLaSM operators
circle = Lar.circle()()
GL.VIEW([ GL.GLFrame2, GL.GLGrid( circle...,GL.COLORS[1],1 ) ]);

#	1D complex in 1D generated by Lar and PLaSM operators
intervals = Lar.qn(5)([2,-.5])
GL.VIEW([ GL.GLFrame2, GL.GLGrid( intervals...,GL.COLORS[1],1 ) ]);

#	2D complex in 3D generated by Lar and PLaSM operators
cylinders = Lar.larModelProduct([ circle, intervals ])
GL.VIEW([ GL.GLFrame2, GL.GLGrid( cylinders...,GL.COLORS[1],1 ) ]);

#	other visualizations
V,FV = cylinders
FW = [[f] for f in FV]
GL.VIEW( push!( GL.GLExplode( V,FW, 1.2,1.2,1.2, 1,0.2 ), GL.GLFrame2) );
GL.VIEW( push!( GL.GLExplode( V,FW, 1.,1.,1., 99,1 ), GL.GLFrame2) );

# see
function schlegel3D(t::Float64,d::Float64)
    function schlegel3D0(V)
        W = Array{Float64,2}(undef,size(V))
        for k=1:size(V,2)
            W[4,k] = V[4,k] + t # translation
            W[1,k] = V[1,k]*d/W[4,k] # projection
            W[2,k] = V[2,k]*d/W[4,k] # projection
            W[3,k] = V[3,k]*d/W[4,k] # projection
        end
        return W[1:3,:]
    end
    return schlegel3D0
end

# 2D complex in 4D generated as topological product of circles
toroidal = Lar.larModelProduct([ Lar.circle(8)(), Lar.circle(2)() ])
V,FV = toroidal

# removal of 4-th coodinate: ortho-projection 4D -> 3D
GL.VIEW([ GL.GLFrame2, GL.GLGrid( V[1:3,:], FV, GL.COLORS[1],0.2 ) ]);

# schlegel projection 4D -> 3D
W = schlegel3D(-6.,2.)(V)
GL.VIEW([ GL.GLFrame2, GL.GLGrid( W, FV, GL.COLORS[1],0.2 ) ]);