text:
# functions to run StochasticGene for Bayesian parameter estimation of stochastic Markov models of gene transcription

using Dates
using DelimitedFiles

"""
    fit_rna(nchains::Int,gene::String,cell::String,datacond,G::Int,maxtime::Float64,infolder::String,resultfolder::String,datafolder,inlabel,label,nsets,runcycle::Bool=false,transient::Bool=false,samplesteps::Int=40000,warmupsteps=0,annealsteps=0,temp=100.,tempanneal=100.,root = "/home/carsonc/scrna/")
    fit_rna(nchains::Int,gene::String,cell::String,fittedparam::String,datacond,G::Int,maxtime::Float64,infolder::String,resultfolder::String,datafolder,inlabel,label,nsets,runcycle::Bool=false,transient::Bool=false,samplesteps::Int=100000,warmupsteps=20000,annealsteps=0,temp=1.,tempanneal=100.,root = "/home/carsonc/scrna/")
    fit_rna(nchains::Int,gene::String,cell::String,fittedparam::String,fixedeffects::String,datacond,G::Int,maxtime::Float64,infolder::String,resultfolder::String,datafolder,inlabel,label,nsets,runcycle::Bool=false,transient::Bool=false,samplesteps::Int=100000,warmupsteps=20000,annealsteps=0,temp=1.,tempanneal=100.,root = "/home/carsonc/scrna/")
    fit_rna(nchains::Int,gene::String,cell::String,fittedparam::Vector,fixedeffects::Tuple,datacond,G::Int,maxtime::Float64,infolder::String,resultfolder::String,datafolder,inlabel,label,nsets,runcycle::Bool=false,transient::Bool=false,samplesteps::Int=100000,warmupsteps=20000,annealsteps=100000,temp=1.,tempanneal=100.,root = "/home/carsonc/scrna/")
    fit_rna(nchains::Int,data,gene::String,decayrate::Float64,nalleles::Int,fittedparam::Vector,fixedeffects::Tuple,datacond,G::Int,maxtime::Float64,infolder::String,resultfolder::String,datafolder,inlabel,label,nsets,runcycle::Bool=false,transient::Bool=false,samplesteps::Int=100000,warmupsteps=20000,annealsteps=100000,temp=1.,tempanneal=100.,root = "/home/carsonc/scrna/")

Fit steady state or transient GM model to RNA data for a single gene, write the result (through function finalize), and return nothing.

# Arguments
- `nchains`: number of MCMC chains
- `gene`: gene name
- `cell`: cell type
- `datacond`: condition, if more than one condition is used enter as a single string separated by hyphen, e.g. "DMSO-AUXIN"
- `maxtime`: float maximum time for entire run
- `infolder`: folder pointing to results used as initial conditions
- `resultfolder`: folder where results go
- `datafolder`: folder for data
- `inlabel`: name of input files (not including gene name but including condition)
- `label`: = name of output files
- `nsets`: int number of rate sets
- `runcycle`: if true, cycle through all parameters sequentially in MCMC
- `samplesteps`: int number of samples
- `warmupsteps`: int number of warmup steps
- `annealsteps`: in number of annealing steps
- `temp`: MCMC temperature
- `tempanneal`: starting temperature for annealing
- `root`: root folder of data and Results folders
- `fittedparam`: indices of parameters to be fit (input as string of ints separated by "-")
- `fixedeffects`: string indicating which rate is fixed, e.g. "eject"

"""
function fit_rna(nchains::Int,gene::String,cell::String,datacond,G::Int,maxtime::Float64,infolder::String,resultfolder::String,datafolder,inlabel,label,nsets,runcycle::Bool=false,transient::Bool=false,samplesteps::Int=100000,warmupsteps=20000,annealsteps=0,temp=1.,tempanneal=100.,root = "/home/carsonc/scrna/")
    fittedparam = make_fittedparam(G,nsets)
    fit_rna(nchains,gene,cell,fittedparam,(),datacond,G,maxtime,infolder,resultfolder,datafolder,inlabel,label,nsets,runcycle,transient,samplesteps,warmupsteps,annealsteps,temp,tempanneal,root)
    nothing
end

function fit_rna(nchains::Int,gene::String,cell::String,fittedparam::String,datacond,G::Int,maxtime::Float64,infolder::String,resultfolder::String,datafolder,inlabel,label,nsets,runcycle::Bool=false,transient::Bool=false,samplesteps::Int=100000,warmupsteps=20000,annealsteps=0,temp=1.,tempanneal=100.,root = "/home/carsonc/scrna/")
    fittedparam = parse.(Int,split(fittedparam,"-"))
    fit_rna(nchains,gene,cell,fittedparam,(),datacond,G,maxtime,infolder,resultfolder,datafolder,inlabel,label,nsets,runcycle,transient,samplesteps,warmupsteps,annealsteps,temp,tempanneal,root)
    nothing
end

function fit_rna(nchains::Int,gene::String,cell::String,fittedparam::String,fixedeffects::String,datacond,G::Int,maxtime::Float64,infolder::String,resultfolder::String,datafolder,inlabel,label,nsets,runcycle::Bool=false,transient::Bool=false,samplesteps::Int=100000,warmupsteps=20000,annealsteps=0,temp=1.,tempanneal=100.,root = "/home/carsonc/scrna/")
    fittedparam = parse.(Int,split(fittedparam,"-"))
    fixedeffects = make_fixedeffects(fixedeffects,G,nsets)
    fit_rna(nchains,gene,cell,fittedparam,fixedeffects,datacond,G,maxtime,infolder,resultfolder,datafolder,inlabel,label,nsets,runcycle,transient,samplesteps,warmupsteps,annealsteps,temp,tempanneal,root)
    nothing
end

function fit_rna(nchains::Int,gene::String,cell::String,fittedparam::Vector,fixedeffects::Tuple,datacond,G::Int,maxtime::Float64,infolder::String,resultfolder::String,datafolder,inlabel,label,nsets,runcycle::Bool=false,transient::Bool=false,samplesteps::Int=100000,warmupsteps=20000,annealsteps=100000,temp=1.,tempanneal=100.,root = "/home/carsonc/scrna/")
    println(now())
    gene = check_genename(gene)
    println(gene)
    decayrate = decay(gene,cell,root)
    nalleles = alleles(gene,cell,root)
    if decayrate < 0
        throw("decayrate < 0")
    else
        println("decay rate: ",decayrate)
    end

    if occursin("-",datafolder)
        datafolder = string.(split(datafolder,"-"))
    end
    if occursin("-",datacond)
        datacond = string.(split(datacond,"-"))
    end
    if transient
        data = make_data(gene,datacond,datafolder,label,root,["T0","T30","T120"],[0.,30.,120.])
    else
        data = make_data(gene,datacond,datafolder,label,root)
    end
    fit_rna(nchains,data,gene,decayrate,nalleles,fittedparam,fixedeffects,datacond,G,maxtime,infolder,resultfolder,datafolder,inlabel,label,nsets,runcycle,transient,samplesteps,warmupsteps,annealsteps,temp,tempanneal,root)
end

function check_genename(gene)
    if occursin("[",gene)
        gene = replace(gene,"[" => "(")
        gene = replace(gene,"]" => ")")
    end
    return gene
end

function fit_rna(nchains::Int,data,gene::String,decayrate::Float64,nalleles::Int,fittedparam::Vector,fixedeffects::Tuple,datacond,G::Int,maxtime::Float64,infolder::String,resultfolder::String,datafolder,inlabel,label,nsets,runcycle::Bool=false,transient::Bool=false,samplesteps::Int=100000,warmupsteps=20000,annealsteps=100000,temp=1.,tempanneal=100.,root = "/home/carsonc/scrna/")

    printinfo(gene,G,datacond,datafolder,infolder,resultfolder,maxtime)
    println(data.nRNA)

    yieldprior = .05
    r = getr(gene,G,nalleles,decayrate,inlabel,infolder,nsets,root,data)
    cv = getcv(gene,G,nalleles,fittedparam,inlabel,infolder,root)
    if cv != 0.02
        warmupsteps = 0
    end
    if runcycle
        maxtime /= 2
        cv = .02
        r = cycle(nchains,data,r,G,nalleles,nsets,cv,fittedparam,fixedeffects,decayrate,yieldprior,maxtime,temp,tempanneal)
        annealsteps = 0
        warmupsteps = div(samplesteps,5)
    end

    options,model,param = make_structures(r,G,nalleles,nsets,cv,fittedparam,fixedeffects,decayrate,yieldprior,samplesteps,annealsteps,warmupsteps,maxtime,temp,tempanneal)
    initial_ll(param,data,model)
    fit,stats,waic = StochasticGene.run_mh(data,model,options,nchains);
    finalize(data,model,fit,stats,waic,temp,resultfolder,root)
    println(now())
    nothing
end

function fit_fish(nchains::Int,gene::String,datacond,G::Int,maxtime::Float64,infolder::String,resultfolder::String,datafolder,inlabel,label,nsets,runcycle::Bool=false,transient::Bool=false,samplesteps::Int=100000,warmupsteps=20000,annealsteps=0,temp=1.,tempanneal=100.,root = "/home/carsonc/scrna/")

end

"""
rna_fish(gene,cond,fishfolder,rnafolder,yield,root)

output RNA histogram and downsampled FISH histogram with loss
"""
function rna_fish(gene,cond,fishfolder,rnafolder,yield,root)
    datarna = StochasticGene.histograms_rna(StochasticGene.scRNApath(gene,cond,rnafolder,root),gene,false)
    f=reduce_fish(gene,cond,datarna[1],fishfolder,yield,root)
    return datarna[2],f
end

"""
reduce_fish(gene,cond,nhist,fishfolder,yield,root)

sample fish histogram with loss (1-yield)
"""
function reduce_fish(gene,cond,nhist,fishfolder,yield,root)
    datafish = StochasticGene.histograms_rna(StochasticGene.FISHpath(gene,cond,fishfolder,root),gene,true)
    StochasticGene.technical_loss(datafish[2],yield,nhist)
end

"""
make_structures(r,G,nalleles,nsets,cv,fittedparam,fixedeffects,decayrate,yieldprior,samplesteps,annealsteps,warmupsteps,maxtime,temp,tempanneal)

return structures options, model, and param
"""
function make_structures(r,G,nalleles,nsets,cv,fittedparam,fixedeffects,decayrate,yieldprior,samplesteps,annealsteps,warmupsteps,maxtime,temp,tempanneal)
    if length(fixedeffects) > 0
        model = StochasticGene.model_rna(r,G,nalleles,nsets,cv,fittedparam,fixedeffects,decayrate,yieldprior,0)
    else
        model = StochasticGene.model_rna(r,G,nalleles,nsets,cv,fittedparam,decayrate,yieldprior,0)
    end
    options = StochasticGene.MHOptions(samplesteps,annealsteps,warmupsteps,0,maxtime,temp,tempanneal)
    param,_ = StochasticGene.initial_proposal(model)
    return options,model,param
end

"""
cycle(nchains,data,r,G,nalleles,nsets,cv,fittedparam,fixedeffects,decayrate,yieldprior,maxtime,temp,tempanneal)

run_mh by cycling through individual parameters sequentially
"""
function cycle(nchains,data,r,G,nalleles,nsets,cv,fittedparam,fixedeffects,decayrate,yieldprior,maxtime,temp,tempanneal)
    options,model,param = make_structures(r,G,nalleles,nsets,cv,fittedparam,fixedeffects,decayrate,yieldprior,100,0,0,maxtime/10,temp,tempanneal)
    initial_ll(param,data,model,"pre-cycle ll:")
    t0 = time()
    while (time() - t0 < maxtime)
        for fp in fittedparam
            model = StochasticGene.model_rna(r,G,nalleles,nsets,cv,[fp],decayrate,yieldprior,0)
            fit,stats,waic = StochasticGene.run_mh(data,model,options,nchains);
            r = StochasticGene.get_rates(fit.parml,model)
        end
    end
    return r
end
"""
initial_ll(param,data,model,message="initial ll:")

compute and print initial loglikelihood
"""
function initial_ll(param,data,model,message="initial ll:")
    ll,_ = StochasticGene.loglikelihood(param,data,model)
    println(message,ll)
end

"""
printinfo(gene,G,cond,datafolder,infolder,resultfolder,maxtime)

print out run information
"""
function printinfo(gene,G,cond,datafolder,infolder,resultfolder,maxtime)
    println(datafolder)
    println(gene," ",G," ",cond)
    println("in: ", infolder," out: ",resultfolder)
    println(maxtime)
end

"""
finalize(data,model,fit,stats,waic,temp,resultfolder,root)

write out run results and print out final loglikelihood and deviance
"""
function finalize(data,model,fit,stats,waic,temp,resultfolder,root)
    writefile = joinpath(root,resultfolder)
    StochasticGene.writeall(writefile,fit,stats,waic,data,temp,model)
    println("final ll: ",fit.llml)
    println(fit.accept," ",fit.total)
    println("Deviance: ",StochasticGene.deviance(fit,data,model))
    println(stats.meanparam)
end

"""
transient_rna(nchains,gene::String,fittedparam,cond::Vector,G::Int,maxtime::Float64,infolder::String,resultfolder,datafolder,inlabel,label,nsets,runcycle::Bool=false,samplesteps::Int=40000,warmupsteps=20000,annealsteps=100000,temp=1.,tempanneal=100.,root = "/home/carsonc/scrna/")

make structures for transient model fit
"""
function transient_rna(nchains,gene::String,fittedparam,cond::Vector,G::Int,maxtime::Float64,infolder::String,resultfolder,datafolder,inlabel,label,nsets,runcycle::Bool=false,samplesteps::Int=40000,warmupsteps=20000,annealsteps=100000,temp=1.,tempanneal=100.,root = "/home/carsonc/scrna/")
    data = make_data(gene,cond,G,datafolder,label,nsets,root)
    model = make_model(gene,G,fittedparam,inlabel,infolder,nsets,root,data)
    param,_ = StochasticGene.initial_proposal(model)
    return param, data, model
end

"""
steadystate_rna(gene::String,fittedparam,cond,G::Int,folder::String,datafolder,label,nsets,root)
steadystate_rna(r::Vector,gene::String,fittedparam,cond,G::Int,datafolder::String,label,nsets,root)

make structures for steady state rna fit

"""
function steadystate_rna(gene::String,fittedparam,cond,G::Int,folder::String,datafolder,label,nsets,root)
    datacond = string.(split(cond,"-"))
    data = make_data(gene,datacond,datafolder,label,root)
    model = make_model(gene,G,fittedparam,label * "_" * cond,folder,nsets,root,data,false)
    param,_ = StochasticGene.initial_proposal(model)
    return param, data, model
end

function steadystate_rna(r::Vector,gene::String,fittedparam,cond,G::Int,datafolder::String,label,nsets,root)
    datacond = string.(split(cond,"-"))
    data = make_data(gene,datacond,datafolder,label,root)
    model = make_model(gene,r,G,fittedparam,nsets,root)
    param,_ = StochasticGene.initial_proposal(model)
    return param, data, model
end

"""
make_model(gene,G,fittedparam,inlabel,infolder,nsets,root,data,verbose=true)
make_model(gene,r::Vector,G,fittedparam,nsets,root) = StochasticGene.model_rna(r,G,alleles(root,gene),nsets,.02,fittedparam,r[2*G],.1,0)

make model structure

"""
function make_model(gene,G,fittedparam,inlabel,infolder,nsets,root,data,verbose=true)
    decayrate = decay(root,gene)
    if verbose
        println(decayrate)
        if decayrate < 0
            throw("error")
        end
    end
    nalleles = alleles(root,gene)
    yieldprior = .1
    r = getr(gene,G,nalleles,decayrate,inlabel,infolder,nsets,root,data,verbose)
    cv = getcv(gene,G,nalleles,fittedparam,inlabel,infolder,root,verbose)
    model = StochasticGene.model_rna(r,G,nalleles,nsets,cv,fittedparam,decayrate,yieldprior,0)
end

make_model(gene,r::Vector,G,fittedparam,nsets,root) = StochasticGene.model_rna(r,G,alleles(root,gene),nsets,.02,fittedparam,r[2*G],.1,0)

function make_fixedeffects(fixedeffects,G,nsets)
    if fixedeffects == "eject"
        fixedeffects = ([2*G-1,2*G*(nsets-1) + 2*G-1],)
    elseif fixedeffects == "off"
        fixedeffects = ([2*G-2,2*G*(nsets-1) + 2*G-2],)
    elseif fixedeffects == "on"
        fixedeffects = ([2*G-3,2*G*(nsets-1) + 2*G-3],)
    elseif fixedeffects == "offeject" || fixedeffects == "ejectoff"
        fixedeffects = ([2*G-1,2*G*(nsets-1) + 2*G-1], [2*G-2,2*G*(nsets-1) + 2*G-2])
    elseif fixedeffects == "oneject" || fixedeffects == "ejecton"
        fixedeffects = ([2*G-1,2*G*(nsets-1) + 2*G-1],[2*G-3,2*G*(nsets-1) + 2*G-3])
    elseif fixedeffects == "onoff" || fixedeffects == "offon"
        fixedeffects = ([2*G-2,2*G*(nsets-1) + 2*G-2],[2*G-3,2*G*(nsets-1) + 2*G-3])
    end
    return fixedeffects
end
function make_fittedparam(G::Int,nsets)

    if nsets == 1
        if G == 3
            fittedparam = [1,2,3,4,5,7]
        elseif G == 2
            fittedparam = [1,2,3,5]
        elseif G == 1
            fittedparam = [1,3]
        end
    else
        if G == 3
            # fittedparam = [1,2,3,4,5,7,8,9,10,11,13]
            # fittedparam = [7,8,9,10,11,13]
            fittedparam = [7,8,9,10,11]
        elseif G == 2
            # fittedparam = [1,2,3,5,6,7,9]
            # fittedparam = [5,6,7,9]
            fittedparam = [5,6,7]
        elseif G == 1
            fittedparam = [3,5]
        end
    end
    return fittedparam
end

function make_data(gene::String,cond::String,datafolder,label,root)
    if cond == "null"
        cond = ""
    end
    datafile = StochasticGene.scRNApath(gene,cond,datafolder,root)
    StochasticGene.data_rna(datafile,label,gene,false)
end

function make_data(gene::String,cond::Array,datafolder::String,label,root)
    datafile = Array{String,1}(undef,length(cond))
    for i in eachindex(cond)
        datafile[i] = StochasticGene.scRNApath(gene,cond[i],datafolder,root)
    end
    StochasticGene.data_rna(datafile,label,gene,false)
end

function make_data(gene::String,cond::Array,datafolder::Array,label,root)
    datafile = Array{String,1}(undef,length(cond))
    for i in eachindex(cond)
        datafile[i] = StochasticGene.scRNApath(gene,cond[i],datafolder[i],root)
    end
    println(datafile)
    StochasticGene.data_rna(datafile,label,gene,false)
end

function make_data(gene::String,cond::String,datafolder,label,root,sets::Vector,time::Vector)
    if cond == "null"
        cond = ""
    end
    datafile =[]
    for set in sets
        folder = joinpath(datafolder,set)
        datafile = vcat(datafile,StochasticGene.scRNApath(gene,cond,folder,root))
    end
    StochasticGene.data_rna(datafile,label,times,gene,false)
end

function getr(gene,G,nalleles,decayrate,inlabel,infolder,nsets::Int,root,data,verbose=true)
    ratefile = StochasticGene.path_Gmodel("rates",gene,G,nalleles,inlabel,infolder,root)
    if verbose
        println(ratefile)
    end
    if isfile(ratefile)
        r = StochasticGene.readrates(ratefile,2)
        r[end] = clamp(r[end],eps(Float64),1-eps(Float64))
        if length(r) == 2*G*nsets + 1
            if verbose
                println(r)
            end
            return r
        end
    end
    println("No r")
    setr(G,decayrate,nsets,data)
end

function setr(G,decayrate,nsets,data)
    if G == 2
        r = [0.015,0.015,0.5,.01,1.]*decayrate/.01
    elseif G == 3
        r = [0.015,.2,.2,0.015,1.5,.01,1.]*decayrate/.01
    elseif G == 1
        if typeof(data.nRNA) <: Vector
            r = Array{Float64,1}(undef,0)
            for hist in data.histRNA
                mu = StochasticGene.mean_histogram(hist)
                r = vcat([10*mu,1.],r)
            end
            r *= decayrate
            r = [r;1.]
            nsets = 1
        else
            mu=StochasticGene.mean_histogram(data.histRNA)
            r = [10*mu,1.,1.]*decayrate
        end
    end
    if nsets > 1
        r = [r[1:end-1];r]
    end
    r[end] = .05
    println(r)
    return r
end

function getcv(gene,G,nalleles,fittedparam,inlabel,infolder,root,verbose = true)
    paramfile = StochasticGene.path_Gmodel("param_stats",gene,G,nalleles,inlabel,infolder,root)
    if isfile(paramfile)
        cv = StochasticGene.read_covlogparam(paramfile)
        cv = float.(cv)
        if ~StochasticGene.isposdef(cv) || size(cv)[1] != length(fittedparam)
            cv = .02
        end
    else
        cv = .02
    end
    if verbose
        println(cv)
    end
    return cv
end

function decay(gene,cell,root,col=2)
    folder = joinpath(root,"data/halflives")
    files = readdir(folder)
    a = nothing
    for file in files
        if occursin(cell,file)
            path = joinpath(folder,file)
            in = readdlm(path,',')
            ind = findfirst(in[:,1] .== gene)
            if ind != nothing
                a = in[ind,col]
            end
        end
    end
    decay(a,gene)
end

decay(a::Float64) = log(2)/a/60.

function decay(a,gene)
    if typeof(a) <: Number
        return decay(a)
    else
        println(gene," has no decay time")
        return -1.
    end
end

function alleles(gene,cell,root,col=3)
    folder = joinpath(root,"data/alleles")
    files = readdir(folder)
    a = nothing
    for file in files
        if occursin(cell,file)
            path = joinpath(folder,file)
            in = readdlm(path)
            ind = findfirst(in[:,1] .== gene)
            if ind != nothing
                a = in[ind,col]
            end
        end
    end
    if a != nothing
        return Int(a)
    else
        return 2
    end
end

function repair_measures(resultfolder,measurefile,ratefile::String,cond,n::Int,datafolder,root)
    front,back = split(measurefile,".")
    fix = front * "fix" * "." * back
    measurepath = joinpath(resultfolder,measurefile)
    ratepath = joinpath(resultfolder,ratefile)
    rates,_ = readdlm(ratepath,',',header=true)
    println(length(rates[:,1]))
    measures,head = readdlm(measurepath,',',header=true)
    f = open(joinpath(resultfolder,fix),"w")
    writedlm(f,[head "Deviance fixed" "LogML fixed" "AIC fixed"],',')
    for i in eachindex(rates[:,1])
        h,hd = histograms(rates[i,:],cond,n,datafolder,root)
        d = StochasticGene.deviance(h,hd)
        ll = StochasticGene.crossentropy(h,hd)
        a = 4*(n+1) + 2*ll
        writedlm(f,[measures[i:i,:] d ll a],',')
    end
    close(f)
end

function repair_measures(outfile,file,fittedparam,cond,G::Int,folder::String,datafolder,label,nsets,root)
    contents,header = readdlm(file,',',header=true)
    f = open(outfile,"w")
    writedlm(f,[header],',')
    for row in eachrow(contents)
        row[3],row[6],row[2] = measures(string(row[1]),fittedparam,cond,G,folder,datafolder,label,nsets,root)
        writedlm(f, [row],',')
    end
    close(f)
end

function measures(gene::String,fittedparam,cond,G::Int,folder::String,datafolder,label,nsets,root)
    param,data,model = steadystate_rna(gene,fittedparam,cond,G,folder,datafolder,label,nsets,root)
    ll = StochasticGene.loglikelihood(param,data,model)[1]
    return ll, StochasticGene.aic(length(fittedparam),ll), StochasticGene.deviance(data,model)
end

function deviance(r,cond,n,datafolder,root)
    h,hd = histograms(r,cond,n,datafolder,root)
    StochasticGene.deviance(h,hd)
end

function compute_deviance(outfile,ratefile::String,cond,n,datafolder,root)
    f = open(outfile,"w")
    rates,head = readdlm(ratefile,',',header=true)
    for r in eachrow(rates)
        d=deviance(r,cond,n,datafolder,root)
        writedlm(f,[gene d],',')
    end
    close(f)
end

function write_histograms(outfolder,ratefile::String,cond,n,datafolder,root)
    rates,head = readdlm(ratefile,',',header=true)
    for r in eachrow(rates)
        h,hd = histograms(r,cond,n,datafolder,root)
        f = open(joinpath(outfolder,r[1] * ".txt"),"w")
        writedlm(f,h)
        close(f)
    end
end

function write_histograms(outfolder,ratefile,fittedparam,datacond,G::Int,datafolder::String,label,nsets,root)
    rates,head = readdlm(ratefile,',',header=true)
    if ~isdir(outfolder)
        mkpath(outfolder)
    end
    cond = string.(split(datacond,"-"))
    for r in eachrow(rates)
        h = histograms(r,fittedparam,datacond,G,datafolder,label,nsets,root)
        for i in eachindex(cond)
            f = open(joinpath(outfolder,string(r[1]) * cond[i] * ".txt"),"w")
            writedlm(f,h[i])
            close(f)
        end
    end
end

function histograms(r,cond,n,datafolder,root)
    gene = String(r[1])
    datafile = StochasticGene.scRNApath(gene,cond,datafolder,root)
    hd = StochasticGene.read_scrna(datafile)
    h = StochasticGene.steady_state(r[2:2*n+3],r[end],n,length(hd),alleles(root,gene))
    return h,hd
end

function histograms(rin,fittedparam,cond,G::Int,datafolder,label,nsets,root)
    gene = string(rin[1])
    r = float.(rin[2:end])
    param,data,model = steadystate_rna(r,gene,fittedparam,cond,G,datafolder,label,nsets,root)
    StochasticGene.likelihoodarray(r,data,model)
end

function write_burst_stats(outfile,infile::String,G::String,folder,cond,root)
    folder = joinpath(root,folder)
    condarray = split(cond,"-")
    g = parse(Int,G)
    lr = 2*g
    lc = 2*g-1
    freq = Array{Float64,1}(undef,2*length(condarray))
    burst = similar(freq)
    f = open(joinpath(folder,outfile),"w")
    contents,head = readdlm(joinpath(folder,infile),',',header=true)
    label = Array{String,1}(undef,0)
    for c in condarray
        label = vcat(label, "Freq " * c, "sd","Burst Size " * c, "sd")
    end
    writedlm(f,["gene" reshape(label,1,length(label))],',')
    for r in eachrow(contents)
        gene = String(r[1])
        rates = r[2:end]
        rdecay = decay(root,gene)
        cov = StochasticGene.read_covparam(joinpath(folder,getfile("param_stats",gene,G,folder,cond)[1]))
        # mu = StochasticGene.readmean(joinpath(folder,getfile("param_stats",gene,G,folder,cond)[1]))
        if size(cov,2) < 2
            println(gene)
        end
        for i in eachindex(condarray)
            j = i-1
            freq[2*i-1], freq[2*i] = frequency(rates[1+lr*(i-1)],sqrt(cov[1+lc*j,1+lc*j]),rdecay)
            burst[2*i-1], burst[2*i] = burstsize(rates[3+lr*j],rates[2+lr*j],cov[3+lc*j,3+lc*j],cov[2+lc*j,2+lc*j],cov[2+lc*j,3+lc*j])
        end
        writedlm(f,[gene freq[1] freq[2] burst[1] burst[2] freq[3] freq[4] burst[3] burst[4]],',')
        flush(f)
    end
    close(f)
end

frequency(ron,sd,rdecay) = (ron/rdecay, sd/rdecay)

function burstsize(reject::Float64,roff,covee,covoo,coveo::Float64)
        v = StochasticGene.var_ratio(reject,roff,covee,covoo,coveo)
        return reject/roff, sqrt(v)
end

function ratestats(gene,G,folder,cond)
    filestats=joinpath(folder,getfile("param_stats",gene,G,folder,cond)[1])
    filerates = joinpath(folder,getratefile(gene,G,folder,cond)[1])
    rates = StochasticGene.readrates(filerates)
    cov = StochasticGene.read_covparam(filestats)
    # mu = StochasticGene.readmean(filestats)
    # sd = StochasticGene.readsd(filestats)
    return rates, cov
end

meanofftime(gene::String,infile,n,method,root) = sum(1 .- offtime(gene,infile,n,method,root))

function meanofftime(r::Vector,n::Int,method::Int)
    if n == 1
        return 1/r[1]
    else
        return sum(1 .- offtime(r,n,method))
    end
end

function offtime(r::Vector,n::Int,method::Int)
    _,_,TI = StochasticGene.mat_G_DT(r,n)
    vals,_ = StochasticGene.eig_decompose(TI)
    minval = min(minimum(abs.(vals[vals.!=0])),.2)
    StochasticGene.offtimeCDF(collect(1.:5/minval),r,n,TI,method)
end

function offtime(gene::String,infile,n,method,root)
    contents,head = readdlm(infile,',',header=true)
    r = float.(contents[gene .== contents[:,1],2:end-1])[1,:]
    offtime(r,n,method)

end

function write_moments(outfile,genelist,cond,datafolder,root)
    f = open(outfile,"w")
    writedlm(f,["Gene" "Expression Mean" "Expression Variance"],',')
    for gene in genelist
        datafile = StochasticGene.scRNApath(gene,cond,datafolder,root)
        # data = StochasticGene.data_rna(datafile,label,gene,false)
        len,h = StochasticGene.histograms_rna(datafile,gene,false)
        # h,hd = histograms(r,cond,n,datafolder,root)
        writedlm(f,[gene StochasticGene.mean_histogram(h) StochasticGene.var_histogram(h)],',')
    end
    close(f)
end

function join_files(file1::String,file2::String,outfile::String,addlabel::Bool=true)
    contents1,head1 = readdlm(file1,',',header=true)   # model G=2
    contents2,head2 = readdlm(file2,',',header=true)   # model G=3
    f = open(outfile,"w")
    if addlabel
        header = vcat(String.(head1[2:end]) .* "_G2",String.(head2[2:end]) .* "_G3")
    else
        header = vcat(String.(head1[2:end]),String.(head2[2:end]))
    end
    header = reshape(permutedims(header),(1,length(head1)+length(head2)-2))
    header = hcat(head1[1],header)
    println(header)
    writedlm(f,header,',')
    for row in 1:size(contents1,1)
        if contents1[row,1] == contents2[row,1]
            contents = hcat(contents1[row:row,2:end],contents2[row:row,2:end])
            contents = reshape(permutedims(contents),(1,size(contents1,2)+size(contents2,2)-2))
            contents = hcat(contents1[row,1],contents)
            writedlm(f,contents,',')
        end
    end
    close(f)
end

function join_files(models::Array,files::Array,outfile::String,addlabel::Bool=true)
    m = length(files)
    contents = Array{Array,1}(undef,m)
    headers = Array{Array,1}(undef,m)
    len = 0
    for i in 1:m
        contents[i],headers[i] = readdlm(files[i],',',header=true)
        len += length(headers[i][2:end])
    end
    f = open(outfile,"w")
    header = Array{String,1}(undef,0)
    for i in 1:m
        if addlabel
            header = vcat(header,String.(headers[i][2:end]) .* "_G$(models[i])")
        else
            header = vcat(header,String.(headers[i][2:end]))
        end
    end
    header = reshape(permutedims(header),(1,len))
    header = hcat(headers[1][1],header)
    println(header)
    writedlm(f,header,',')
    for row in 1:size(contents[1],1)
        content = contents[1][row:row,2:end]
        for i in 1:m-1
            if contents[i][row,1] == contents[i+1][row,1]
                content = hcat(content,contents[i+1][row:row,2:end])
                # content = reshape(permutedims(content),(1,len))
            end
        end
        content = hcat(contents[1][row:row,1],content)
        writedlm(f,[content],',')
    end
    close(f)
end

function best_AIC(outfile,infile)
    contents,head = readdlm(infile,',',header=true)
    head = vec(head)
    ind = occursin.("AIC",string.(head)) .& .~ occursin.("WAIC",string.(head))
    f = open(outfile,"w")
    labels = "Gene"
    for i in 1:sum(ind)
        labels = vcat(labels,"Model $(i)")
    end
    labels = vcat(labels,"Winning Model")
    writedlm(f,[reshape(labels,1,length(labels))],',')
    for row in eachrow(contents)
        writedlm(f,[row[1] reshape(row[ind],1,length(row[ind])) argmin(float.(row[ind]))],',')
    end
    close(f)
end

function sample_non1_genes(infile,n)
    contents,head = readdlm(infile,',',header=true)
    list = Array{String,1}(undef,0)
    for c in eachrow(contents)
        if c[5] != 1
            push!(list,c[1])
        end
    end
    a = sample(list,n,replace=false)
end

function best_model(file::String)
    contents,head = readdlm(file,',',header=true)
    f = open(file,"w")
    head = hcat(head,"Winning Model")
    writedlm(f,head,',')
    for row in eachrow(contents)
        if abs(row[11] - row[4]) > row[5] + row[12]
            if row[11] < row[4]
                writedlm(f, hcat(permutedims(row),3),',')
            else
                writedlm(f, hcat(permutedims(row),2),',')
            end
        else
            if row[13] < row[6]
                writedlm(f, hcat(permutedims(row),3),',')
            else
                writedlm(f, hcat(permutedims(row),2),',')
            end
        end
    end
    close(f)
end

function best_waic(folder,root)
    folder = joinpath(root,folder)
    files =readdir(folder)
    lowest = Inf
    winner = ""
    for file in files
        if occursin("measures",file)
            contents,head = readdlm(joinpath(folder,file),',',header=true)
            waic = mean(float.(contents[:,4]))
            println(mean(float.(contents[:,2]))," ",waic," ",median(float.(contents[:,4]))," ",file)
            if waic < lowest
                lowest = waic
                winner = file
            end
        end
    end
    return winner,lowest
end

function bestmodel(measures2,measures3)
    m2,head = readdlm(infile,',',header=true)
    m3,head = readdlm(infile,',',header=true)
end

function add_best_burst(filein,fileout,filemodel2,filemodel3)
    contents,head = readdlm(filein,',',header=true)
    burst2,head2 = readdlm(filemodel2,',',header=true)
    burst3,head3 = readdlm(filemodel3,',',header=true)
    f = open(fileout,"w")
    head = hcat(head,["mean off period" "bust size"])
    writedlm(f,head,',')
    for row in eachrow(contents)
        if Int(row[end]) == 2
            writedlm(f, hcat(permutedims(row),permutedims(burst2[findfirst(burst2[:,1] .== row[1]),2:3])),',')
        else
            writedlm(f, hcat(permutedims(row),permutedims(burst3[findfirst(burst3[:,1] .== row[1]),2:3])),',')
        end
    end
    close(f)
end

function add_best_occupancy(filein,fileout,filemodel2,filemodel3)
    contents,head = readdlm(filein,',',header=true)
    occupancy2,head2 = readdlm(filemodel2,',',header=true)
    occupancy3,head3 = readdlm(filemodel3,',',header=true)
    f = open(fileout,"w")
    head = hcat(head,["Off -2" "Off -1" "On State" ])
    writedlm(f,head,',')
    for row in eachrow(contents)
        if Int(row[end-2]) == 2
            writedlm(f, hcat(permutedims(row),hcat("NA",permutedims(occupancy2[findfirst(occupancy2[:,1] .== row[1]),2:end]))),',')
        else
            writedlm(f, hcat(permutedims(row),permutedims(occupancy3[findfirst(occupancy3[:,1] .== row[1]),2:end])),',')
        end
    end
    close(f)
end

function plot_model(r,n,nhist,nalleles)
    h= StochasticGene.steady_state(r[1:2*n+2],r[end],n,nhist,nalleles)
    plot(h)
    return h
end

function prune_file(list,file,outfile,header=true)
    contents,head = readdlm(file,',',header=header)
    f = open(outfile,"w")
    for c in eachrow(contents)
        if c[1] in list
            writedlm(f,[c],',')
        end
    end
    close(f)
end

# precompile(fit_rna,(String,String,Int,Float64,String,String,String,String,String,Int,String,Bool))
# fit_rna(nchains,gene::String,cond::String,G::Int,maxtime::Float64,infolder::String,resultfolder,datafolder,inlabel,label,nsets,runcycle::Bool=false,sample::Int=40000,warmup=20000,anneal=100000,temp=1.,tempanneal=100.,root = "/home/carsonc/scrna/")

max_stars_repo_name: nih-niddk-mbs/Transcription

max_issues_repo_name: nih-niddk-mbs/Transcription

max_forks_repo_name: nih-niddk-mbs/Transcription


text:
"""
The carbon index bounds are static per year and defined in the
[Methodology](https://carbon-intensity.github.io/api-definitions/#carbon-intensity-api-v2-0-0) paper published by NG-ESO
"""
const CARBON_INDEX_BOUNDS = OrderedDict(
    "Very Low" => (0, 44),
    "Low" => (45, 129),
    "Moderate" => (130, 209),
    "High" => (210, 310),
    "Very High" => (310, 400)
)

"""
Dict mapping region names to region ids.
"""
const AVAILABLE_REGIONS = Dict(
    "North Scotland" => 1,
    "South Scotland" => 2,
    "North West England" => 3,
    "North East England" => 4,
    "Yorkshire" => 5,
    "North Wales & Merseyside" => 6,
    "South Wales" => 7,
    "West Midlands" => 8,
    "East Midlands" => 9,
    "East England" => 10,
    "South West England" => 11,
    "South England" => 12,
    "London" => 13,
    "South East England" => 14,
    "England" => 15,
    "Scotland" => 16,
    "Wales" => 17,
    "GB" => 18
)

ext: jl


text:
@testset "C functions {$ST}" begin

 v_ptr = PETSc.C.VecCreate(ST)
 chk(C.VecSetType(v_ptr, PETSc.C.VECMPI))
 b = Vec{ST}(v_ptr)
 resize!(b, mlocal=sys_size)
 global_indices = localpart(b) - 1  # zero based
  for i=1:sys_size
    idxm = [global_indices[i] ]
    val = [ rhs[i] ]
    PETSc.C.SetValues(b.p, idxm, val)
  end

  PETSc.C.AssemblyBegin(b.p)
  PETSc.C.AssemblyEnd(b.p)

   for i=1:sys_size
    idxm = [global_indices[i] ]
    val = ST[ 0.0 ]
    PETSc.C.GetValues(b.p, idxm, val)
    @test val[1] == rhs[i]
  end

  # matrix
  A = Mat(ST, mlocal=sys_size, nlocal=sys_size)
  row_range, col_range = localranges(A)

  for i=1:sys_size
    for j=1:sys_size
      row_idx = PETSc.C.PetscInt[ row_range[i] - 1 ]
      col_idx = PETSc.C.PetscInt[ row_range[j] - 1 ]
      val = ST[ A_julia[i, j] ]
      PETSc.C.SetValues(A.p, row_idx, col_idx, val)
    end
  end

  PETSc.C.AssemblyBegin(A.p)
  PETSc.C.AssemblyEnd(A.p)

  for i=1:sys_size
    for j=1:sys_size
      row_idx = PETSc.C.PetscInt[ row_range[i] - 1 ]
      col_idx = PETSc.C.PetscInt[ row_range[j] - 1 ]
      val = ST[ 0.0 ]
      PETSc.C.GetValues(A.p, row_idx, col_idx, val)
      @test val[1] == A_julia[i, j]
    end
  end

  idx = collect(row_range - 1)
  vals = collect(1:(sys_size*sys_size))
  vals2 = convert(Array{ST, 1}, vals)
  PETSc.C.SetValues(A.p, idx, idx, vals2)
  PETSc.C.AssemblyBegin(A.p)
  PETSc.C.AssemblyEnd(A.p)

  for i=1:sys_size
    for j=1:sys_size
      row_idx = PETSc.C.PetscInt[ row_range[i] - 1 ]
      col_idx = PETSc.C.PetscInt[ row_range[j] - 1 ]
      val = ST[ 0.0 ]
      PETSc.C.GetValues(A.p, row_idx, col_idx, val)
      @test val[1] == (i + (j-1)*sys_size)
    end
  end

  # block matrix
  B = Mat(ST, mlocal=sys_size, nlocal=sys_size, bs=sys_size, mtype=PETSc.C.MATMPIBAIJ)
  idx = PETSc.C.PetscInt[comm_rank]
  PETSc.C.SetValues(B.p, idx, idx, vals2)
  PETSc.C.AssemblyBegin(B.p)
  PETSc.C.AssemblyEnd(B.p)
  for i=1:sys_size
    for j=1:sys_size
      row_idx = PETSc.C.PetscInt[ row_range[i] - 1 ]
      col_idx = PETSc.C.PetscInt[ row_range[j] - 1 ]
      val = ST[ 0.0 ]
      PETSc.C.GetValues(A.p, row_idx, col_idx, val)
      @test val[1] == (i + (j-1)*sys_size)
    end
  end

  # shell matrix
  if ST == Float64
    ctx = (1, 2, 3)
    ctx_ptr = pointer_from_objref(ctx)
    c_ptr = PETSc.C.MatCreateShell(sys_size, sys_size, PETSc.C.PETSC_DETERMINE, PETSc.C.PETSC_DETERMINE, ctx_ptr)
    C = Mat{ST}(c_ptr)

    f_ptr = cfunction(mymult, PETSc.C.PetscErrorCode, (PETSc.C.Mat{ST}, PETSc.C.Vec{ST}, PETSc.C.Vec{ST}))
    PETSc.C.MatShellSetOperation(C.p, PETSc.C.MATOP_MULT, f_ptr)

    ctx_ret = PETSc.C.MatShellGetContext(C.p)
    @test unsafe_pointer_to_objref(ctx_ret) == ctx

    x = Vec(ST[1:sys_size;])
    xlocal = LocalVector(x)
    b = Vec(zeros(ST, sys_size))
    *(C, x, b)
    gc()  # avoid a finalizer problem
    blocal = LocalVector(b)
    for i=1:length(blocal)
      @test blocal[i] == ST(i*i)
    end

    restore(blocal)
  end
end


text:
let
    setprecision(BigFloat, 256)

    tol = 1e-20
    p = 3
    n = 20
    t = map(BigFloat, (1:n)//n)
    Ut,Vt = egrss.generators(t,p);
    d = BigFloat("1e-8")*ones(BigFloat,n)
    Wt,c = egrss.potrf(Ut,Vt,d)
    Yt,Zt = egrss.trtri(Ut,Wt,c)

    K = egrss.full(Ut,Vt,d)
    Lref = cholesky(K).L
    L = egrss.full_tril(Ut,Wt,c)
    L_inv = egrss.full_tril(Yt,Zt,1.0./c)

    @test ( norm(I-L*L_inv) < n*tol )
    @test ( norm(L-Lref) < norm(Lref)*tol )
end

max_stars_repo_name: martinandersen/EGRSS

max_issues_repo_name: martinandersen/EGRSS

max_forks_repo_name: martinandersen/EGRSS


text:
# Conversion to and from Base types
using StationXML, Test
using InteractiveUtils: subtypes

@testset "Conversion" begin
    @testset "NumberType" begin
        @testset "$T" for T in subtypes(StationXML.NumberType)
            val = rand()
            @test convert(T, val).value == val
            @test convert(Float64, T(val)) == val
            @test convert(Complex{Float32}, T(val)) == Float32(val) + 0im
        end
    end

    @testset "EnumeratedStruct" begin
        @testset "$T" for T in subtypes(StationXML.EnumeratedStruct)
            val = first(StationXML.permitted_values(T))
            @test convert(T, val).value == val
            @test convert(String, T(val)) == val
        end
    end

    @testset "Identifier" begin
        let value = "some random text"
            @test convert(String, StationXML.Identifier(value)) == value
            @test convert(StationXML.Identifier, value).value == value
        end
    end

    @testset "Email" begin
        let address = "test.name@example.com"
            @test convert(String, StationXML.Email(address)) == address
            @test convert(StationXML.Email, address).value == address
            @test_throws ArgumentError convert(StationXML.Email, "not an email address")
        end
    end

    @testset "Comment" begin
        let comment = "Some text comment"
            @test convert(StationXML.Comment, comment).value == comment
        end
    end

    @testset "Units" begin
        let unit = "AMPERES"
            @test convert(StationXML.Units, unit).name == unit
            @test convert(String, StationXML.Units(unit)) == unit
        end
    end
end


text:
using GridPotentialFlow
using Plots
Δx = 0.02
Lx = 4.0
xlim = (-Lx/2,Lx/2)
ylim = (-Lx/2,Lx/2)
g = PhysicalGrid(xlim,ylim,Δx);

Rc = Lx/4
Δs = 2Δx
body = PotentialFlowBody(Circle(Rc,Δs))

Rv = 3/2*Rc
Γv = 1.0
v = Vortex(Rv,0.0,Γv);

setΓ(body,-Γv)

model = VortexModel(g,vortices=[v],bodies=[body]);
sol = solve(model);

r1 = sqrt.((body.points.x.-Rv).^2 .+ (body.points.y).^2) # distance from vortex to points on circle
r2 = sqrt.((body.points.x.-Rc^2/Rv).^2 .+ (body.points.y).^2) # distance from image vortex to points on circle
θ1 = π.+atan.((body.points.y)./(body.points.x.-Rv)) # angle from vortex to points on circle
θ2 = atan.((body.points.y),(body.points.x.-Rc^2/Rv)) # angle from image vortex to points on circle
v1x = -v.Γ./(2π*r1).*cos.(θ1.-π/2) # x velocity on circle induced by vortex
v1y = -v.Γ./(2π*r1).*sin.(θ1.-π/2) # y velocity on circle induced by vortex
v2x = v.Γ./(2π*r2).*cos.(θ2.-π/2) # x velocity on circle induced by image vortex
v2y = v.Γ./(2π*r2).*sin.(θ2.-π/2) # y velocity on circle induced by image vortex
V = sqrt.((v1x.+v2x).^2+(v1y.+v2y).^2) # velocity magnitude on circle
γ = -V; # bound vortex sheet strength on circle (velocity on circle is clockwise if positive vortex is to the right of it)

plot(sol.ψ,g)
plot!(body,fillcolor=:black,fillrange=0,fillalpha=0.25,linecolor=:black,linewidth=2)
scatter!([v.x],[v.y],color=:black,markersize=2,xlabel="x",ylabel="y")

plot(sol.f./Δs,label="f/ds",xlabel="body point index")
plot!(γ,label="gamma")

Vv = Γv/(2π*(Rv-Rc^2/Rv))
T = 2π*Rv/Vv

import OrdinaryDiffEq
function rhs(X,model,t)
    setvortexpositions!(model,X)
    Ẋ = vortexvelocities!(model)
    return Ẋ
end
X = getvortexpositions(model)
prob = OrdinaryDiffEq.ODEProblem(rhs,X,(0.0,T),model);
sol = OrdinaryDiffEq.solve(prob,dt=0.1,OrdinaryDiffEq.RK4(),dense=false,adaptive=false);
plot(body,fillcolor=:black,fillrange=0,fillalpha=0.25,linecolor=:black,linewidth=2)
plot!(map(s->s.u[1],sol.u),map(s->s.v[1],sol.u))
scatter!([v.x],[v.y],color=:black,markersize=2,xlabel="x",ylabel="y")

@testset "Vortex near cylinder" begin
    import OrdinaryDiffEq
#jl

    Rv = sqrt(v.x^2+v.y^2)
    Vθ = -v.Γ/(2π*(Rv-Rc^2/Rv))
    Tv = 2π*Rv/abs(Vθ)
#jl

    tspan = (0,Tv)
    model = VortexModel(g,bodies=[body],vortices=[v])
#jl
    function rhs(X,model,t)
        setvortexpositions!(model,X)
        Ẋ = vortexvelocities!(model)
        return Ẋ
    end
#jl
    X = getvortexpositions(model)
    prob = OrdinaryDiffEq.ODEProblem(rhs,X,tspan,model);
    sol = OrdinaryDiffEq.solve(prob,dt=0.1,OrdinaryDiffEq.RK4(),dense=false,adaptive=false);
#jl
    @test isapprox(sol.u[end].u[1], 1.5; atol = 1e-2)
    @test isapprox(sol.u[end].v[1], 0.0; atol = 1e-1)
end

body = PotentialFlowBody(Circle(Rc,Δs), U=(1.0,0.0))
model = VortexModel(g,bodies=[body]);
sol = solve(model);

θ = atan.(body.points.y,body.points.x)
γ = 2*sin.(θ);

plot(sol.ψ,g)
plot!(body,fillcolor=:black,fillrange=0,fillalpha=0.25,linecolor=:black,linewidth=2,xlabel="x",ylabel="y")

plot(sol.f./Δs,label="f/ds",xlabel="body point index")
plot!(γ,label="gamma")

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
using EzXML

# <svg version="1.0" id="Livello_1" xmlns="http://www.w3.org/2000/svg" xmlns:xlink="http://www.w3.org/1999/xlink" x="0px" y="0px" viewBox="0 0 283.5 283.5" enable-background="new 0 0 283.5 283.5" xml:space="preserve">

fn = "test.xml"
# fn = "lar4.svg"
println(isfile(fn))
doc = EzXML.readxml(fn)
# Select <a/> links that contain a non-blank text node.
# links = findall("//book/@category", doc)
links = findall("//path/@d", doc)
# links = findall("//path/@d", doc)
println(links)
# width = ndigits(length(links))
for (i, link) in enumerate(links)
#     println(lpad(i, width), ": ", strip(nodecontent(link)), " -- ", link["href"])
    println(": ", strip(nodecontent(link)))
end

# [strip(nodecontent(link)) for link in links]
# links = findall("//year/@tag", doc)
# links = findall("//book[@category and normalize-space(text()) != '']", doc)


text:
## Method definitions for generic eltypes - will make copies

# Direct and transposed algorithms

function materialize!(L::Ldiv{<:BandedColumnMajor})
    A, B = L.A, L.B
    checksquare(A)
    ldiv!(factorize(A), B)
end

function copyto!(dest::AbstractVecOrMat, L::Ldiv{<:BandedRowMajor})
    A, B = L.A, L.B
    copyto!(dest, Ldiv(transpose(factorize(transpose(A))), B))
end

function copyto!(dest::AbstractVecOrMat, L::Ldiv{<:ConjLayout{<:BandedRowMajor}})
    A, B = L.A, L.B
    copyto!(dest, Ldiv(factorize(A')', B))
end

# Direct and transposed algorithms
function ldiv!(A::BandedLU{T,<:BandedMatrix}, B::StridedVecOrMat{T}) where {T<:BlasFloat}
    m = size(A.factors,1)
    l,u = bandwidths(A.factors)
    data = bandeddata(A.factors)
    LAPACK.gbtrs!('N', l, u-l, m, data, A.ipiv, B)
end

function ldiv!(A::BandedLU{T}, B::AbstractVecOrMat{Complex{T}}) where T<:Real
    # c2r = reshape(transpose(reinterpret(T, reshape(B, (1, length(B))))), size(B, 1), 2*size(B, 2))
    a = real.(B)
    b = imag.(B)
    ldiv!(A, a)
    ldiv!(A, b)
    B .= a .+ im.*b
end

function ldiv!(transA::Transpose{T,<:BandedLU{T,<:BandedMatrix}}, B::StridedVecOrMat{T}) where {T<:BlasFloat}
    A = transA.parent
    m = size(A.factors,1)
    l,u = bandwidths(A.factors)
    data = bandeddata(A.factors)
    LAPACK.gbtrs!('T', l, u-l, m, data, A.ipiv, B)
end

function ldiv!(transA::Transpose{<:Any,<:BandedLU{<:Any,<:BandedMatrix}}, B::AbstractVecOrMat)
    A = transA.parent
    ldiv!(transpose(UnitLowerTriangular(A.factors)), ldiv!(transpose(UpperTriangular(A.factors)), B))
    _apply_inverse_ipiv_rows!(A, B)
end

ldiv!(adjF::Adjoint{T,<:BandedLU{T,<:BandedMatrix}}, B::AbstractVecOrMat{T}) where {T<:Real} =
    (F = adjF.parent; ldiv!(transpose(F), B))
function ldiv!(adjA::Adjoint{T,<:BandedLU{T,<:BandedMatrix}}, B::StridedVecOrMat{T}) where {T<:BlasComplex}
    A = adjA.parent
    m = size(A.factors,1)
    l,u = bandwidths(A.factors)
    data = bandeddata(A.factors)
    LAPACK.gbtrs!('C', l, u-l, m, data, A.ipiv, B)
end

function ldiv!(adjA::Adjoint{<:Any,<:BandedLU{<:Any,<:BandedMatrix}}, B::AbstractVecOrMat)
    error("Implement")
    A = adjA.parent
    ldiv!(adjoint(UnitLowerTriangular(A.factors)), ldiv!(adjoint(UpperTriangular(A.factors)), B))
    _apply_inverse_ipiv!(A, B)
end

\(A::Adjoint{<:Any,<:BandedLU}, B::Adjoint{<:Any,<:AbstractVecOrMat}) = A \ copy(B)
\(A::Transpose{<:Any,<:BandedLU}, B::Transpose{<:Any,<:AbstractVecOrMat}) = A \ copy(B)

factorize(A::BandedMatrix) = lu(A)


text:
using DiffEqFinancial, StochasticDiffEq
using Test

# write your own tests here
u0 = [1.0; 0.5]
σ = 0.25
prob = HestonProblem(1.0, 1.0, σ, 1.0, 1.0, u0, (0.0, 1.0))
sol = solve(prob, SRIW1(), adaptive = false, dt = 1 / 10)

prob = BlackScholesProblem((t) -> t^2, (t, u) -> 1.0, σ, 0.5, (0.0, 1.0))
sol = solve(prob, SRIW1())


text:
using Test
using OpenQuantumSystems
using Random, SparseArrays, LinearAlgebra, StableRNGs
import QuantumOpticsBase

import DelayDiffEq

@testset "master" begin

    Random.seed!(StableRNG(0), 1)

    D(op1::Array, op2::Array) = abs(norm(op1 - op2))
    D(x1::StateVector, x2::StateVector) = norm(x2 - x1)
    D(op1::AbstractOperator, op2::AbstractOperator) =
        abs(tracedistance_nh(dense(op1), dense(op2)))
    D(op1::AbstractSuperOperator, op2::AbstractSuperOperator) =
        abs(tracedistance_nh(dense(op1), dense(op2)))

    mode1 = Mode(0.2, 1.0)
    mode2 = Mode(0.3, 2.0)
    Energy = [0.0, 200.0]
    mol1 = Molecule([mode1], 3, [2.0, 200.0])
    mol2 = Molecule([mode2], 3, [3.0, 300.0])
    aggCore = AggregateCore([mol1, mol2])
    aggCore.coupling[2, 3] = 50
    aggCore.coupling[3, 2] = 50
    agg = setupAggregate(aggCore)
    aggTools = agg.tools
    aggOperators = agg.operators

    Ham_B = agg.operators.Ham_B
    Ham_I = agg.operators.Ham_I
    Ham_0 = agg.operators.Ham_0
    Ham = agg.operators.Ham

    Ham_0_lambda, Ham_0_S = eigen(Ham_0.data)
    Ham_0_Sinv = inv(Ham_0_S)
    Ham_0_lambda = diagm(Ham_0_lambda)

    basis = agg.tools.basis
    indicesLen = agg.tools.bSize
    indices = agg.tools.indices
    indicesMap = agg.tools.indicesMap
    FCFact = agg.tools.FCfactors
    FCProd = agg.tools.FCproduct

    data = Matrix(Hermitian(rand(ComplexF64, indicesLen, indicesLen)))
    W0 = DenseOperator(basis, basis, data)
    rho0 = trace_bath(W0, aggCore, aggTools)
    W0_bath = get_rho_bath(W0, aggCore, aggTools)
    normalize!(rho0)
    # tests have to be quick enough
    t_max = 0.001
    t_count = 10
    t0 = 0.
    t_step = (t_max - t0) / (t_count)
    tspan = [t0:t_step:t_max;]

    _, W_t = QME_SS_exact(
        W0,
        tspan,
        agg;
        reltol = 1e-3,
        abstol = 1e-3,
        int_reltol = 1e-4,
        int_abstol = 1e-4,
        alg = DelayDiffEq.MethodOfSteps(DelayDiffEq.Tsit5()),
    )

    _, W_int_t = QME_SI_exact(
        W0,
        tspan,
        agg;
        reltol = 1e-3,
        abstol = 1e-3,
        int_reltol = 1e-4,
        int_abstol = 1e-4,
        alg = DelayDiffEq.MethodOfSteps(DelayDiffEq.Tsit5()),
    )

    _, rho_t = QME_sS_exact(
        W0,
        tspan,
        agg;
        reltol = 1.0e-3,
        abstol = 1.0e-3,
        int_reltol = 1.0e-4,
        int_abstol = 1.0e-4,
        alg = DelayDiffEq.MethodOfSteps(DelayDiffEq.Tsit5())
    )

    _, rho_int_t = QME_sI_exact(
        W0,
        tspan,
        agg;
        reltol = 1.0e-3,
        abstol = 1.0e-3,
        int_reltol = 1.0e-4,
        int_abstol = 1.0e-4,
        alg = DelayDiffEq.MethodOfSteps(DelayDiffEq.Tsit5())
    )

end # testset

max_stars_repo_name: detrin/OpenQuantumSystems

max_issues_repo_name: detrin/OpenQuantumSystems

max_forks_repo_name: detrin/OpenQuantumSystems


text:
@inline VectorizationBase.stridedpointer(A::PtrArray) = A.ptr
@inline VectorizationBase.stridedpointer(A::StrideArray) = A.ptr.ptr

ArrayInterface.device(::AbstractStrideArray) = ArrayInterface.CPUPointer()

# AbstractStrideArray{S,D,T,N,C,B,R,X,O}

# contiguous_axis(A), contiguous_batch_size(A), stride_rank(A), bytestrides(A), offsets(A))

ArrayInterface.contiguous_axis(::Type{<:AbstractStrideArray{S,D,T,N,C}}) where {S,D,T,N,C} = ArrayInterface.Contiguous{C}()
ArrayInterface.contiguous_batch_size(::Type{<:AbstractStrideArray{S,D,T,N,C,B}}) where {S,D,T,N,C,B} = ArrayInterface.ContiguousBatch{B}()
ArrayInterface.stride_rank(::Type{<:AbstractStrideArray{S,D,T,N,C,B,R}}) where {S,D,T,N,C,B,R} = ArrayInterface.StrideRank{R}()
ArrayInterface.dense_dims(::Type{<:AbstractStrideArray{S,D}}) where {S,D} = DenseDims{D}()


text:
using CompScienceMeshes, BEAST

Γ = readmesh(joinpath(dirname(@__FILE__),"sphere2.in"))
X = raviartthomas(Γ)
Y = buffachristiansen(Γ)
Z = raviartthomas(Γ, BEAST.Continuity{:none})

ϵ, μ, ω = 1.0, 1.0, 1.0; κ = ω * √(ϵ*μ)
K, N = Maxwell3D.doublelayer(wavenumber=κ), NCross()
NK, Id = BEAST.DoubleLayerRotatedMW3D(im*κ), Identity()
E = Maxwell3D.planewave(direction=ẑ, polarization=x̂, wavenumber=κ)
H = -1/(im*μ*ω)*curl(E)
h = (n × H) × n
q = n × H

@hilbertspace j
@hilbertspace m
@hilbertspace g

mfie1 = @discretise (K+0.5N)[m,j] == h[m]  j∈X m∈Y
mfie2 = @discretise (NK-0.5Id)[g,j] == q[m]  j∈Z g∈Z

u1 = gmres(mfie1)
u2 = gmres(mfie2)

fcr1, _ = facecurrents(u1, X)
fcr2, _ = facecurrents(u2, Z)

using LinearAlgebra
using Plots

plot(title="Compare current density")
plot!(norm.(fcr1), label="MxMIFE")
scatter!(norm.(fcr2), label="DGMFIE", c=:red)

include("utils/postproc.jl")
include("utils/plotresults.jl")


text:
module Events

using Base: UUID
import UUIDs
import Dates
import Random
import Redis
import CBOR
import Mongoc
import Mongoc: BSON

export Event, EventDispatcher, EventException, StreamInfo, get_dispatcher

struct EventException <: Exception
    msg::AbstractString
end

@enum ServiceState stopped=0 running

include("./events/Event.jl")
include("./events/EventDispatcher.jl")
include("./events/EventStore.jl")

function get_dispatcher(group_id, consumer_id; claim_time=Dates.Minute(15))
    # TODO: Parse these from a URI
    host = get(ENV, "REDIS_DISPATCHER_HOST", "127.0.0.1")
    port = parse(UInt16, get(ENV, "REDIS_DISPATCHER_PORT", "6379"))
    password = get(ENV, "REDIS_DISPATCHER_PASSWORD", "")
    db = parse(UInt, get(ENV, "REDIS_DISPATCHER_DB", "0"))
    RedisEventDispatcher(group_id, consumer_id; host, port, password, db, claim_time)
end

function get_event_store(;auto_snapshot=AUTO_SNAPSHOT)
    uri = get(ENV, "MONGO_EVENTSTORE_URI", DEFAULT_MONGO_URI)
    MongoEventStore(;uri, auto_snapshot)
end

end


text:
@testset "sampling_functions" begin
    @testset "sampling input $Tx" for Tx in [Matrix, ColVecs, RowVecs]
        @testset "bayesian_linear_regression" begin
            rng, N, D = MersenneTwister(123456), 11, 5
            X, f, Σy = generate_toy_problem(rng, N, D, Tx)

            g = rand(rng, f)
            @test g(X) == g(X)  # check the sample doesn't change between evaluations

            if X isa Matrix
                Xc = ColVecs(X)
                Xr = RowVecs(X')
                @test g(X) == g(Xc)
                @test g(X) == g(Xr)
            end

            # test the Random interface
            @test rand(rng, Random.Sampler(rng, f, Val(Inf))) isa
                BayesianLinearRegressors.BLRFunctionSample

            samples1, samples2 = 1000, 1000
            samples = samples1 * samples2
            gs = rand(rng, f, samples1, samples2)
            @test size(gs) == (samples1, samples2)

            # test statistical properties of the sampled functions
            let
                Y = reduce(hcat, map(h -> h(X), reshape(gs, :)))
                m_empirical = mean(Y; dims=2)
                Σ_empirical = (Y .- mean(Y; dims=2)) * (Y .- mean(Y; dims=2))' ./ samples
                @test mean(f(X, Σy)) ≈ m_empirical atol = 1e-3 rtol = 1e-3
                @test cov(f(X, Σy)) ≈ Σ_empirical + Σy atol = 1e-2 rtol = 1e-2
            end

            # test statistical properties of in-place rand
            let
                A = Array{BayesianLinearRegressors.BLRFunctionSample,2}(
                    undef, samples1, samples2
                )
                A = rand!(rng, A, f)
                Y = reduce(hcat, map(h -> h(X), reshape(gs, :)))
                m_empirical = mean(Y; dims=2)
                Σ_empirical = (Y .- mean(Y; dims=2)) * (Y .- mean(Y; dims=2))' ./ samples
                @test mean(f(X, Σy)) ≈ m_empirical atol = 1e-3 rtol = 1e-3
                @test cov(f(X, Σy)) ≈ Σ_empirical + Σy atol = 1e-2 rtol = 1e-2
            end

            @testset "Zygote (everything dense)" begin
                function test_rand_funcs_adjoints(sample_function)
                    rng, N, D = MersenneTwister(123456), 11, 5
                    X, f, _ = generate_toy_problem(rng, N, D, Tx)
                    mw, A_Λw = f.mw, 0.1 .* randn(rng, D, D)

                    # Run the model forwards and check that output agrees with non-Zygote.
                    z, back = Zygote.pullback(sample_function, X, mw, A_Λw)
                    @test z == sample_function(X, mw, A_Λw)

                    # Compute adjoints using Zygote.
                    z̄ = randn(rng, size(z))
                    dX, dmw, dA_Λw = back(z̄)

                    # Verify adjoints via finite differencing.
                    fdm = central_fdm(5, 1)
                    @test dX ≈ first(j′vp(fdm, X -> sample_function(X, mw, A_Λw), z̄, X))
                    @test dmw ≈ first(j′vp(fdm, mw -> sample_function(X, mw, A_Λw), z̄, mw))
                    @test dA_Λw ≈
                        first(j′vp(fdm, A_Λw -> sample_function(X, mw, A_Λw), z̄, A_Λw))
                end

                function rand_funcs_single(X, mw, A_Λw)
                    Λw = Symmetric(A_Λw * A_Λw' + I)
                    f = BayesianLinearRegressor(mw, Λw)
                    g = rand(MersenneTwister(123456), f)
                    return g(X)
                end

                function rand_funcs_multi(X, mw, A_Λw)
                    Λw = Symmetric(A_Λw * A_Λw' + I)
                    f = BayesianLinearRegressor(mw, Λw)
                    gs = rand(MersenneTwister(123456), f, 1, 1)
                    return reduce(hcat, map(h -> h(X), reshape(gs, :)))
                end

                test_rand_funcs_adjoints(rand_funcs_single)
                test_rand_funcs_adjoints(rand_funcs_multi)
            end
        end

        @testset "basis_function_regression" begin
            rng, N, D = MersenneTwister(123456), 11, 2
            X, f, Σy = generate_toy_problem(rng, N, D, Tx)

            f_bf = BasisFunctionRegressor(f, ϕ)

            samples = 1_000_000
            gs = rand(rng, f_bf, samples)

            # test statistical properties of the sampled functions
            let
                Y = reduce(hcat, map(h -> h(X), reshape(gs, :)))
                m_empirical = mean(Y; dims=2)
                Σ_empirical = (Y .- mean(Y; dims=2)) * (Y .- mean(Y; dims=2))' ./ samples
                @test mean(f_bf(X, Σy)) ≈ m_empirical atol = 1e-3 rtol = 1e-3
                @test cov(f_bf(X, Σy)) ≈ Σ_empirical + Σy atol = 1e-2 rtol = 1e-2
            end
        end
    end
end


text:
using RecipesBase

@recipe function f(y::AbstractArray{<:NeuralQuantum.Measurement})
    yerror := uncertainty.(y)
    value.(y)
end

@recipe function f(x::AbstractArray, y::AbstractArray{<:NeuralQuantum.Measurement})
    yerror := uncertainty.(y)
    x, value.(y)
end

@recipe function f(x::AbstractArray{<:NeuralQuantum.Measurement}, y::AbstractArray)
    xerror := uncertainty.(x)
    value(x), y
end


text:
atreplinit() do repl
    try
        @eval using Pkg
        haskey(Pkg.installed(), "Revise") || @eval Pkg.add("Revise")
    catch
    end
    try
        @eval using Revise
        @async Revise.wait_steal_repl_backend()
    catch
    end
end

max_stars_repo_name: RePsychLing/Docker

max_issues_repo_name: RePsychLing/Docker

max_forks_repo_name: RePsychLing/Docker


text:
"""
    countne(x::UVAL_COLLECTION_TYPES, y::UVAL_COLLECTION_TYPES, n::Int)

Estimate a `n`-member distribution on the number of indices at which the elements of
two collections of uncertain values are not equal.

This is done by repeating the following procedure `n` times:

1. Draw a length-`L` realisation of `x` by drawing one random
    number from  each uncertain value furnishing the dataset. The draws are
    independent, so that no element-wise dependencies (e.g. sequential
    correlations) that are not already present in the data are introduced in
    the realisation.
2. Draw a length-`L` realisation of `y` in the same manner.
3. Count the number of indices at which the elements of the two length-`L`
    draws are not equal.

This yields `n` counts of non-equal values between `n` pairs of independent
realisations of `x` and `y`. The `n`-member distribution of nonequal-value counts
is returned as a vector.
"""
StatsBase.countne(x::UVAL_COLLECTION_TYPES, y::UVAL_COLLECTION_TYPES, n::Int) =
    resample(StatsBase.countne, x, y, n)

"""
    counteq(x::UVAL_COLLECTION_TYPES, y::UVAL_COLLECTION_TYPES, n::Int)

Estimate a `n`-member distribution on the number of indices at which the elements of
two collections of uncertain values are equal.

This is done by repeating the following procedure `n` times:

1. First, draw a length-`L` realisation of `x` by drawing one random
    number from  each uncertain value furnishing the dataset. The draws are
    independent, so that no element-wise dependencies (e.g. sequential
    correlations) that are not already present in the data are introduced in
    the realisation.
2. Second, draw a length-`L` realisation of `y` in the same manner.
3. Count the number of indices at which the elements of the two length-`L`
    draws are equal.

This yields `n` counts of non-equal values between `n` pairs of independent
realisations of `x` and `y`. The `n`-member distribution of equal-value counts
is returned as a vector.
"""
StatsBase.counteq(x::UVAL_COLLECTION_TYPES, y::UVAL_COLLECTION_TYPES, n::Int) =
    resample(StatsBase.counteq, x, y, n)

"""
    corkendall(x::UVAL_COLLECTION_TYPES, y::UVAL_COLLECTION_TYPES, n::Int)

Estimate a `n`-member distribution on Kendalls's rank correlation
coefficient between two collections of uncertain values.

This is done by repeating the following procedure `n` times:

1. First, draw a length-`L` realisation of `x` by drawing one random
    number from  each uncertain value furnishing the dataset. The draws are
    independent, so that no element-wise dependencies (e.g. sequential
    correlations) that are not already present in the data are introduced in
    the realisation.
2. Second, draw a length-`L` realisation of `y` in the same manner.
3. Compute Kendall's rank correlation coefficient between the two length-`L`
    draws.

This yields `n` computations of Kendall's rank correlation coefficient
between `n` independent pairs of realisations of `x` and `y`. The
`n`-member distribution of correlation estimates is returned as a vector.
"""
StatsBase.corkendall(x::UVAL_COLLECTION_TYPES, y::UVAL_COLLECTION_TYPES, n::Int) =
    resample(StatsBase.corspearman, x, y, n)

"""
    corspearman(x::UVAL_COLLECTION_TYPES, y::UVAL_COLLECTION_TYPES, n::Int)

Estimate a `n`-member distribution on Spearman's rank correlation
coefficient between two collections of uncertain values.

This is done by repeating the following procedure `n` times:

1. First, draw a length-`L` realisation of `x` by drawing one random
    number from  each uncertain value furnishing the dataset. The draws are
    independent, so that no element-wise dependencies (e.g. sequential
    correlations) that are not already present in the data are introduced in
    the realisation.
2. Second, draw a length-`L` realisation of `y` in the same manner.
3. Compute Spearman's rank correlation coefficient between the two length-`L`
    draws.

This yields `n` estimates of Spearman's rank correlation coefficient
between `n` independent pairs of realisations of `x` and `y`. The
`n`-member distribution of correlation estimates is returned as a vector.
"""
StatsBase.corspearman(x::UVAL_COLLECTION_TYPES, y::UVAL_COLLECTION_TYPES, n::Int) =
    resample(StatsBase.corspearman, x, y, n)

"""
    cor(x::UVAL_COLLECTION_TYPES, y::UVAL_COLLECTION_TYPES, n::Int)

Estimate a distribution on Pearson's rank correlation coefficient between
two collections of uncertain values.

This is done by repeating the following procedure `n` times:

1. First, draw a length-`L` realisation of `x` by drawing one random
    number from  each uncertain value furnishing the dataset. The draws are
    independent, so that no element-wise dependencies (e.g. sequential
    correlations) that are not already present in the data are introduced in
    the realisation.
2. Second, draw a length-`L` realisation of `y` in the same manner.
3. Compute Pearson's rank correlation coefficient between the two length-`L`
    draws.

This yields `n` estimates of Pearson's rank correlation coefficient
between `n` independent pairs of realisations of `x` and `y`. The
`n`-member distribution of correlation estimates is returned as a vector.
"""
StatsBase.cor(x::UVAL_COLLECTION_TYPES, y::UVAL_COLLECTION_TYPES, n::Int) =
    resample(StatsBase.cor, x, y, n)

"""
    cov(x::UVAL_COLLECTION_TYPES, y::UVAL_COLLECTION_TYPES, n::Int; corrected::Bool = true)

Obtain a distribution on the covariance between two collections of
uncertain values.

This is done by repeating the following procedure `n` times:

1. First, draw a length-`L` realisation of `x` by drawing one random
    number from  each uncertain value furnishing the dataset. The draws are
    independent, so that no element-wise dependencies (e.g. sequential
    correlations) that are not already present in the data are introduced in
    the realisation.
2. Second, draw a length-`L` realisation of `y` in the same manner.
3. Compute the covariance between the two length-`L` draws.

This yields `n` estimates of the covariance between `n` independent pairs
of realisations of `x` and `y`. The `n`-member distribution of covariance
estimates is returned as a vector.

If `corrected` is `true` (the default) then the sum is scaled with `n - 1` for
each pair of draws, whereas the sum is scaled with `n` if `corrected` is `false`
where `n = length(x)`.
"""
StatsBase.cov(x::UVAL_COLLECTION_TYPES, y::UVAL_COLLECTION_TYPES, n::Int; corrected::Bool = true) =
    resample(StatsBase.cov, x, y, n; corrected = corrected)

"""
    crosscor(x::UVAL_COLLECTION_TYPES, y::UVAL_COLLECTION_TYPES, [lags], n::Int; demean = true)

Obtain a distribution over the cross correlation between two collections of
uncertain values.

This is done by repeating the following procedure `n` times:

1. First, draw a length-`L` realisation of `x` by drawing one random
    number from each uncertain value furnishing the dataset. The draws are
    independent, so that no element-wise dependencies (e.g. sequential
    correlations) that are not already present in the data are introduced in
    the realisation.
2. Second, draw a length-`L` realisation of `y` in the same manner.
3. Compute the cross correlation between the two length-`L` draws.

This yields `n` estimates of the cross correlation between `n` independent pairs
of realisations of `x` and `y`. The `n`-member distribution of cross correlation
estimates is returned as a vector.

`demean` specifies whether, at each iteration, the respective means of the draws
should be subtracted from them before computing their cross correlation.

When left unspecified, the `lags` used are `-min(n-1, 10*log10(n))` to `min(n, 10*log10(n))`.

The output is normalized by `sqrt(var(x_draw)*var(y_draw))`. See `crosscov` for the unnormalized form.
"""
function StatsBase.crosscor(x::UVAL_COLLECTION_TYPES, y::UVAL_COLLECTION_TYPES, n::Int;
        demean = true)

    crosscor_estimates = Vector{Vector{Float64}}(undef, n)

    for i = 1:n
        draw_x = float.(resample(x))
        draw_y = float.(resample(y))
        crosscor_estimates[i] = StatsBase.crosscor(draw_x, draw_y, demean = demean)
    end

    return crosscor_estimates
end

function StatsBase.crosscor(x::UVAL_COLLECTION_TYPES, y::UVAL_COLLECTION_TYPES, lags, n::Int;
        demean = true)

    crosscor_estimates = Vector{Vector{Float64}}(undef, n)

    for i = 1:n
        draw_x = float.(resample(x))
        draw_y = float.(resample(y))
        crosscor_estimates[i] = StatsBase.crosscor(draw_x, draw_y, lags, demean = demean)
    end

    return crosscor_estimates
end

"""
    crosscov(x::UVAL_COLLECTION_TYPES, y::UVAL_COLLECTION_TYPES, [lags], n::Int; demean = true)

Obtain a distribution over the cross covariance function (CCF) between two
collections of uncertain values.

This is done by repeating the following procedure `n` times:

1. First, draw a length-`L` realisation of `x` by drawing one random
    number from each uncertain value furnishing the dataset. The draws are
    independent, so that no element-wise dependencies (e.g. sequential
    correlations) that are not already present in the data are introduced in
    the realisation.
2. Second, draw a length-`L` realisation of `y` in the same manner.
3. Compute the CCF between the two length-`L` draws.

This yields `n` estimates of the CCF between `n` independent pairs
of realisations of `x` and `y`. The `n`-member distribution of CCF estimates
 is returned as a vector.

`demean` specifies whether, at each iteration, the respective means of the draws
should be subtracted from them before computing their CCF.

When left unspecified, the `lags` used are `-min(n-1, 10*log10(n))` to `min(n, 10*log10(n))`.

The output is not normalized. See `crosscor` for a function with normalization.
"""
function StatsBase.crosscov(x::UVAL_COLLECTION_TYPES, y::UVAL_COLLECTION_TYPES, n::Int;
        demean = true)

    crosscov_estimates = Vector{Vector{Float64}}(undef, n)

    for i = 1:n
        draw_x = float.(resample(x))
        draw_y = float.(resample(y))
        crosscov_estimates[i] = StatsBase.crosscov(draw_x, draw_y, demean = demean)
    end

    return crosscov_estimates
end

function StatsBase.crosscov(x::UVAL_COLLECTION_TYPES, y::UVAL_COLLECTION_TYPES, lags, n::Int;
        demean = true)

    crosscov_estimates = Vector{Vector{Float64}}(undef, n)

    for i = 1:n
        draw_x = float.(resample(x))
        draw_y = float.(resample(y))
        crosscov_estimates[i] = StatsBase.crosscov(draw_x, draw_y, lags, demean = demean)
    end

    return crosscov_estimates
end

"""
    gkldiv(x::UVAL_COLLECTION_TYPES, y::UVAL_COLLECTION_TYPES, n::Int)

Obtain a distribution over the generalized Kullback-Leibler divergence between two
collections of uncertain values.

This is done by repeating the following procedure `n` times:

1. First, draw a length-`L` realisation of `x` by drawing one random
    number from each uncertain value furnishing the dataset. The draws are
    independent, so that no element-wise dependencies (e.g. sequential
    correlations) that are not already present in the data are introduced in
    the realisation.
2. Second, draw a length-`L` realisation of `y` in the same manner.
3. Compute the  generalized Kullback-Leibler divergence between the two
    length-`L` draws.

This yields `n` estimates of the  generalized Kullback-Leibler divergence
between `n` independent pairs of realisations of `x` and `y`. The `n`-member
distribution of generalized Kullback-Leibler divergence estimates is
returned as a vector.
"""
StatsBase.gkldiv(x::UVAL_COLLECTION_TYPES, y::UVAL_COLLECTION_TYPES, n::Int) =
    resample(StatsBase.gkldiv, x, y, n)

"""
    kldivergence(x::UVAL_COLLECTION_TYPES, y::UVAL_COLLECTION_TYPES, [b], n::Int)

Obtain a distribution over the Kullback-Leibler divergence between two
collections of uncertain values.

This is done by repeating the following procedure `n` times:

1. First, draw a length-`L` realisation of `x` by drawing one random
    number from each uncertain value furnishing the dataset. The draws are
    independent, so that no element-wise dependencies (e.g. sequential
    correlations) that are not already present in the data are introduced in
    the realisation.
2. Second, draw a length-`L` realisation of `y` in the same manner.
3. Compute the Kullback-Leibler divergence between the two
    length-`L` draws.

This yields `n` estimates of the Kullback-Leibler divergence
between `n` independent pairs of realisations of `x` and `y`. The `n`-member
distribution of Kullback-Leibler divergence estimates is
returned as a vector.

Optionally a real number `b` can be specified such that the divergence is
scaled by `1/log(b)`.
"""
StatsBase.kldivergence(x::UVAL_COLLECTION_TYPES, y::UVAL_COLLECTION_TYPES, n::Int) =
    resample(StatsBase.kldivergence, x, y, n)

StatsBase.kldivergence(x::UVAL_COLLECTION_TYPES, y::UVAL_COLLECTION_TYPES, b, n::Int) =
    resample(StatsBase.kldivergence, x, y, n, b)

"""
    maxad(x::UVAL_COLLECTION_TYPES, y::UVAL_COLLECTION_TYPES, n::Int)

Obtain a distribution over the maximum absolute deviation between two
collections of uncertain values.

This is done by repeating the following procedure `n` times:

1. First, draw a length-`L` realisation of `x` by drawing one random
    number from each uncertain value furnishing the dataset. The draws are
    independent, so that no element-wise dependencies (e.g. sequential
    correlations) that are not already present in the data are introduced in
    the realisation.
2. Second, draw a length-`L` realisation of `y` in the same manner.
3. Compute the maximum absolute deviation between the two
    length-`L` draws.

This yields `n` estimates of the maximum absolute deviation
between `n` independent pairs of realisations of `x` and `y`. The `n`-member
distribution of maximum absolute deviation estimates is
returned as a vector.
"""
StatsBase.maxad(x::UVAL_COLLECTION_TYPES, y::UVAL_COLLECTION_TYPES, n::Int) =
    resample(StatsBase.maxad, x, y, n)

"""
    meanad(x::UVAL_COLLECTION_TYPES, y::UVAL_COLLECTION_TYPES, n::Int)

Obtain a distribution over the mean absolute deviation between two
collections of uncertain values.

This is done by repeating the following procedure `n` times:

1. First, draw a length-`L` realisation of `x` by drawing one random
    number from each uncertain value furnishing the dataset. The draws are
    independent, so that no element-wise dependencies (e.g. sequential
    correlations) that are not already present in the data are introduced in
    the realisation.
2. Second, draw a length-`L` realisation of `y` in the same manner.
3. Compute the mean absolute deviation between the two
    length-`L` draws.

This yields `n` estimates of the mean absolute deviation
between `n` independent pairs of realisations of `x` and `y`. The `n`-member
distribution of mean absolute deviation estimates is
returned as a vector.
"""
StatsBase.meanad(x::UVAL_COLLECTION_TYPES, y::UVAL_COLLECTION_TYPES, n::Int) =
    resample(StatsBase.meanad, x, y, n)

"""
    msd(x::UVAL_COLLECTION_TYPES, y::UVAL_COLLECTION_TYPES, n::Int)

Obtain a distribution over the mean squared deviation between two
collections of uncertain values.

This is done by repeating the following procedure `n` times:

1. First, draw a length-`L` realisation of `x` by drawing one random
    number from each uncertain value furnishing the dataset. The draws are
    independent, so that no element-wise dependencies (e.g. sequential
    correlations) that are not already present in the data are introduced in
    the realisation.
2. Second, draw a length-`L` realisation of `y` in the same manner.
3. Compute the mean squared deviation between the two
    length-`L` draws.

This yields `n` estimates of the mean squared deviation
between `n` independent pairs of realisations of `x` and `y`. The `n`-member
distribution of mean squared deviation estimates is
returned as a vector.
"""
StatsBase.msd(x::UVAL_COLLECTION_TYPES, y::UVAL_COLLECTION_TYPES, n::Int) =
    resample(StatsBase.msd, x, y, n)

"""
    psnr(x::UVAL_COLLECTION_TYPES, y::UVAL_COLLECTION_TYPES, maxv, n::Int)

Obtain a distribution over the peak signal-to-noise ratio (PSNR) between two
collections of uncertain values.

This is done by repeating the following procedure `n` times:

1. First, draw a length-`L` realisation of `x` by drawing one random
    number from each uncertain value furnishing the dataset. The draws are
    independent, so that no element-wise dependencies (e.g. sequential
    correlations) that are not already present in the data are introduced in
    the realisation.
2. Second, draw a length-`L` realisation of `y` in the same manner.
3. Compute the PSNR between the two length-`L` draws.

This yields `n` estimates of the peak signal-to-noise ratio
between `n` independent pairs of realisations of `x` and `y`. The `n`-member
distribution of PSNR estimates is returned as a vector.

The PSNR is computed as `10 * log10(maxv^2 / msd(x_draw, y_draw))`, where `maxv` is
the maximum possible value `x` or `y` can take
"""
StatsBase.psnr(x::UVAL_COLLECTION_TYPES, y::UVAL_COLLECTION_TYPES, maxv, n::Int) =
    resample(StatsBase.psnr, x, y, n, maxv)

"""
    rmsd(x::UVAL_COLLECTION_TYPES, y::UVAL_COLLECTION_TYPES, n::Int, normalize = false)

Obtain a distribution over the root mean squared deviation between two
collections of uncertain values.

This is done by repeating the following procedure `n` times:

1. First, draw a length-`L` realisation of `x` by drawing one random
    number from each uncertain value furnishing the dataset. The draws are
    independent, so that no element-wise dependencies (e.g. sequential
    correlations) that are not already present in the data are introduced in
    the realisation.
2. Second, draw a length-`L` realisation of `y` in the same manner.
3. Compute the root mean squared deviation between the two
    length-`L` draws.

This yields `n` estimates of the root mean squared deviation
between `n` independent pairs of realisations of `x` and `y`. The `n`-member
distribution of root mean squared deviation estimates is
returned as a vector.

The root mean squared deviation is computed as `sqrt(msd(x_draw, y_draw))`
at each iteration. Optionally, `x_draw` and `y_draw` may be normalised.
"""
StatsBase.rmsd(x::UVAL_COLLECTION_TYPES, y::UVAL_COLLECTION_TYPES, n::Int; normalize = false) =
    resample(StatsBase.rmsd, x, y, n; normalize = normalize)

"""
    sqL2dist(x::UVAL_COLLECTION_TYPES, y::UVAL_COLLECTION_TYPES, n::Int)

Obtain a distribution over the squared L2 distance between two
collections of uncertain values.

This is done by repeating the following procedure `n` times:

1. First, draw a length-`L` realisation of `x` by drawing one random
    number from each uncertain value furnishing the dataset. The draws are
    independent, so that no element-wise dependencies (e.g. sequential
    correlations) that are not already present in the data are introduced in
    the realisation.
2. Second, draw a length-`L` realisation of `y` in the same manner.
3. Compute the squared L2 distance between the two
    length-`L` draws.

This yields `n` estimates of the squared L2 distance
between `n` independent pairs of realisations of `x` and `y`. The `n`-member
distribution of squared L2 distance estimates is returned as a vector.

The squared L2 distance is computed as ``\\sum_{i=1}^n |x_i - y_i|^2``.
"""
StatsBase.sqL2dist(x::UVAL_COLLECTION_TYPES, y::UVAL_COLLECTION_TYPES, n::Int) =
    resample(StatsBase.sqL2dist, x, y, n)


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

# Filter
b = [1., 2., 3., 4.]
x = [1., 1., 0., 1., 1., 0., 0., 0.]
@test filt(b, 1., x)  == [1., 3., 5., 8., 7., 5., 7., 4.]
@test filt(b, [1., -0.5], x)  == [1., 3.5, 6.75, 11.375, 12.6875, 11.34375, 12.671875, 10.3359375]
# With ranges
@test filt(b, 1., 1.0:10.0) == [1., 4., 10., 20., 30., 40., 50., 60., 70., 80.]
@test filt(1.:4., 1., 1.0:10.0) == [1., 4., 10., 20., 30., 40., 50., 60., 70., 80.]
# Across an array is the same as channel-by-channel
@test filt(b, 1., [x 1.0:8.0]) == [filt(b, 1., x) filt(b, 1., 1.0:8.0)]
@test filt(b, [1., -0.5], [x 1.0:8.0]) == [filt(b, [1., -0.5], x) filt(b, [1., -0.5], 1.0:8.0)]
si = zeros(3)
@test filt(b, 1., [x 1.0:8.0], si) == [filt(b, 1., x, si) filt(b, 1., 1.0:8.0, si)]
@test si == zeros(3) # Will likely fail if/when arrayviews are implemented
si = [zeros(3) ones(3)]
@test filt(b, 1., [x 1.0:8.0], si) == [filt(b, 1., x, zeros(3)) filt(b, 1., 1.0:8.0, ones(3))]
# With initial conditions: a lowpass 5-pole butterworth filter with W_n = 0.25,
# and a stable initial filter condition matched to the initial value
b = [0.003279216306360201,0.016396081531801006,0.03279216306360201,0.03279216306360201,0.016396081531801006,0.003279216306360201]
a = [1.0,-2.4744161749781606,2.8110063119115782,-1.703772240915465,0.5444326948885326,-0.07231566910295834]
si = [0.9967207836936347,-1.4940914728163142,1.2841226760316475,-0.4524417279474106,0.07559488540931815]
@test filt(b, a, ones(10), si) ≈ ones(10) # Shouldn't affect DC offset

@test_throws ArgumentError filt!([1, 2], [1], [1], [1])
@test xcorr([1, 2], [3, 4]) == [4, 11, 6]

@test fftshift([1 2 3]) == [3 1 2]
@test fftshift([1, 2, 3]) == [3, 1, 2]
@test fftshift([1 2 3; 4 5 6]) == [6 4 5; 3 1 2]
@test ifftshift([1 2 3]) == [2 3 1]
@test ifftshift([1, 2, 3]) == [2, 3, 1]
@test ifftshift([1 2 3; 4 5 6]) == [5 6 4; 2 3 1]

# Convolution
a = [1., 2., 1., 2.]
b = [1., 2., 3.]
@test conv(a, b) ≈ [1., 4., 8., 10., 7., 6.]
@test conv(complex(a, ones(4)), complex(b)) ≈ complex([1., 4., 8., 10., 7., 6.], [1., 3., 6., 6., 5., 3.])

# Discrete cosine transform (DCT) tests

if Base.fftw_vendor() != :mkl
    a = rand(8,11) + im*rand(8,11)
    @test norm(idct(dct(a)) - a) < 1e-8

    X = reshape([1,2,7,2,1,5,9,-1,3,4,6,9],3,4)
    Y = rand(17,14)
    Y[3:5,9:12] = X
    sX = view(Y,3:5,9:12)

    true_Xdct = [  13.856406460551018  -3.863239728836245   2.886751345948129  -0.274551994240164; -2.828427124746190  -2.184015211898548  -4.949747468305834   3.966116180118245; 4.898979485566356  -0.194137576915510  -2.857738033247041   2.731723009609389 ]

    true_Xdct_1 = [    5.773502691896258   4.618802153517007   6.350852961085884  10.969655114602890; -4.242640687119286  -2.121320343559643   4.242640687119286  -3.535533905932738; 1.632993161855452   2.041241452319315   5.715476066494083   0.408248290463863 ]

    true_Xdct_2 = [    8.  -3.854030797826254  -3.0  3.761176226848022;
        4.0  -2.071929829606556   4.0  -2.388955165168770; 12.  -0.765366864730179   4.0  -1.847759065022573 ]

    Xdct = dct(X)
    Xdct! = float(X); dct!(Xdct!)
    Xdct_1 = dct(X,1)
    Xdct!_1 = float(X); dct!(Xdct!_1,1)
    Xdct_2 = dct(X,2)
    Xdct!_2 = float(X); dct!(Xdct!_2,2)

    Xidct = idct(true_Xdct)
    Xidct! = copy(true_Xdct); idct!(Xidct!)
    Xidct_1 = idct(true_Xdct_1,1)
    Xidct!_1 = copy(true_Xdct_1); idct!(Xidct!_1,1)
    Xidct_2 = idct(true_Xdct_2,2)
    Xidct!_2 = copy(true_Xdct_2); idct!(Xidct!_2,2)

    pXdct = plan_dct(X)*(X)
    pXdct! = float(X); plan_dct!(pXdct!)*(pXdct!)
    pXdct_1 = plan_dct(X,1)*(X)
    pXdct!_1 = float(X); plan_dct!(pXdct!_1,1)*(pXdct!_1)
    pXdct_2 = plan_dct(X,2)*(X)
    pXdct!_2 = float(X); plan_dct!(pXdct!_2,2)*(pXdct!_2)

    pXidct = plan_idct(true_Xdct)*(true_Xdct)
    pXidct! = copy(true_Xdct); plan_idct!(pXidct!)*(pXidct!)
    pXidct_1 = plan_idct(true_Xdct_1,1)*(true_Xdct_1)
    pXidct!_1 = copy(true_Xdct_1); plan_idct!(pXidct!_1,1)*(pXidct!_1)
    pXidct_2 = plan_idct(true_Xdct_2,2)*(true_Xdct_2)
    pXidct!_2 = copy(true_Xdct_2); plan_idct!(pXidct!_2,2)*(pXidct!_2)

    sXdct = dct(sX)
    psXdct = plan_dct(sX)*(sX)
    sYdct! = copy(Y); sXdct! = view(sYdct!,3:5,9:12); dct!(sXdct!)
    psYdct! = copy(Y); psXdct! = view(psYdct!,3:5,9:12); plan_dct!(psXdct!)*(psXdct!)

    for i = 1:length(X)
        @test Xdct[i] ≈ true_Xdct[i]
        @test Xdct![i] ≈ true_Xdct[i]
        @test Xdct_1[i] ≈ true_Xdct_1[i]
        @test Xdct!_1[i] ≈ true_Xdct_1[i]
        @test Xdct_2[i] ≈ true_Xdct_2[i]
        @test Xdct!_2[i] ≈ true_Xdct_2[i]

        @test pXdct[i] ≈ true_Xdct[i]
        @test pXdct![i] ≈ true_Xdct[i]
        @test pXdct_1[i] ≈ true_Xdct_1[i]
        @test pXdct!_1[i] ≈ true_Xdct_1[i]
        @test pXdct_2[i] ≈ true_Xdct_2[i]
        @test pXdct!_2[i] ≈ true_Xdct_2[i]

        @test Xidct[i] ≈ X[i]
        @test Xidct![i] ≈ X[i]
        @test Xidct_1[i] ≈ X[i]
        @test Xidct!_1[i] ≈ X[i]
        @test Xidct_2[i] ≈ X[i]
        @test Xidct!_2[i] ≈ X[i]

        @test pXidct[i] ≈ X[i]
        @test pXidct![i] ≈ X[i]
        @test pXidct_1[i] ≈ X[i]
        @test pXidct!_1[i] ≈ X[i]
        @test pXidct_2[i] ≈ X[i]
        @test pXidct!_2[i] ≈ X[i]

        @test sXdct[i] ≈ true_Xdct[i]
        @test psXdct[i] ≈ true_Xdct[i]
        @test sXdct![i] ≈ true_Xdct[i]
        @test psXdct![i] ≈ true_Xdct[i]
    end
end # fftw_vendor() != :mkl

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia


text:
@testset "Coordinate system transformations" begin
    # Global coordinates
    lla = LLA(lat=-24.007712762068806, lon=96.44002819812094, alt=1374.7804632852078)
    lla2 = LLA(lat=-24.1, lon=96.5, alt=1659.18948659927)

    ecef_lla = ECEFfromLLA(wgs84)
    lla_ecef = LLAfromECEF(wgs84)

    ecef = ecef_lla(lla)
    ecef2 = ecef_lla(lla2)

    @test ecef ≈ ECEF(-654007.2768949205,5.794061748224952e6,-2.5796231395720555e6)
    @test lla ≈ lla_ecef(ecef)

    lla_utm = LLAfromUTM(utm_zone(lla)..., wgs84)
    utm_lla = UTMfromLLA(utm_zone(lla)..., wgs84)

    utm = utm_lla(lla)
    utm2 = utm_lla(lla2)

    @test utm ≈ UTM(239579.67583179142,7.342551466042985e6,1374.7804632852078)
    @test lla ≈ lla_utm(utm)
    # Test also the poles
    @test UTMfromLLA(0,true,wgs84)(LLA(89.0, 89.0, 89.0)) ≈ UTM(2.111009610242531e6, 1.9980623200455606e6, 89.0)
    @test LLAfromUTM(0,true,wgs84)(UTM(2.111009610242531e6, 1.9980623200455606e6, 89.0)) ≈ LLA(89.0, 89.0, 89.0)

    lla_utmz = LLAfromUTMZ(wgs84)
    utmz_lla = UTMZfromLLA(wgs84)

    utmz = utmz_lla(lla)
    utmz2 = utmz_lla(lla2)
    # Test also the poles go to UPS
    @test utmz_lla(LLA(89.0, 89.0, 89.0)) ≈ UTMZ(2.111009610242531e6, 1.9980623200455606e6, 89.0, 0, true)
    @test lla_utmz(UTMZ(2.111009610242531e6, 1.9980623200455606e6, 89.0, 0, true)) ≈ LLA(89.0, 89.0, 89.0)

    @test utmz ≈ UTMZ(239579.67583179142,7.342551466042985e6,1374.7804632852078, UInt8(47), false)
    @test lla ≈ lla_utmz(utmz)

    utmz_utm = UTMZfromUTM(47, false, wgs84)
    utm_utmz = UTMfromUTMZ(47, false, wgs84)

    @test utmz == utmz_utm(utm)
    @test utm == utm_utmz(utmz)
    # Sometimes this has to do something non-trivial - make sure it does
    utm_utmz2 = UTMfromUTMZ(46, false, wgs84)
    @test utm_utmz2(utmz) ≈ UTM(850009.7418418773, 7.340641097689279e6, 1374.7804632852078)

    # ENU coordinates
    enu_ecef = ENUfromECEF(lla2, wgs84)
    ecef_enu = ECEFfromENU(lla2, wgs84)

    enu = enu_ecef(ecef)

    @test enu_ecef ≈ ENUfromECEF(ecef2, wgs84) # (Can supply origin an any format)
    @test enu_ecef ≈ ENUfromECEF(utm2, 47, false, wgs84)
    @test enu_ecef ≈ ENUfromECEF(utmz2, wgs84)

    @test enu ≈ ENU(-6103.186938282723, 10222.54767562731, -295.55857190545794)
    @test ecef_enu(ENU(-6103.186938282723, 10222.54767562731, -295.55857190545794)) ≈ ecef

    # https://github.com/JuliaGeo/Geodesy.jl/issues/49
    @test LLA(ECEF(1, 0, 0), wgs84) ≈ LLA(89.99866260445, 0.00000000000, -6356752.314234)

    # Inverses
    @test inv(ecef_lla) == lla_ecef
    @test inv(lla_ecef) == ecef_lla
    @test inv(utm_lla) == lla_utm
    @test inv(lla_utm) == utm_lla
    @test inv(utmz_lla) == lla_utmz
    @test inv(lla_utmz) == utmz_lla
    @test inv(utmz_utm) == utm_utmz
    @test inv(utm_utmz) == utmz_utm
    @test inv(enu_ecef) == ecef_enu
    @test inv(ecef_enu) == enu_ecef

    # Web Mercator
    points = [LLA(0, 0, 0),
              LLA(49, 2, 0),
			  LLA(-19, -128, 4),
			  LLA(36, -77, 6)]
    # Reference generated using Proj4
    # [Proj4.transform(Proj4.Projection("+proj=longlat +datum=WGS84"),
    #                  Proj4.Projection("+proj=webmerc +datum=WGS84"),
    #                  [point.lon, point.lat, point.alt]) for point in points]
    points_wm = [[0.0, 0.0, 0.0],
                 [222638.98158654713, 6.274861394006577e6, 0.0],
                 [-1.4248894821539015e7, -2.1549359150858945e6, 4.0],
                 [-8.571600791082064e6, 4.300621372044271e6, 6.0]]
    wm_lla = WebMercatorfromLLA()
    lla_wm = LLAfromWebMercator()
    @testset "WebMercator" for (lla,wm) in zip(points, points_wm)
        @test wm_lla(lla) ≈ wm
        @test lla_wm(wm) ≈ lla
        ll = LatLon(lla.lat, lla.lon)
        wm_2d = wm[1:2]
        @test wm_lla(ll) ≈ wm_2d
        @test lla_wm(wm_2d) ≈ ll
	end

    # Composed transformation functions
    @test LLAfromENU(ecef, wgs84) == LLAfromECEF(wgs84) ∘ ECEFfromENU(ecef, wgs84)
    @test ENUfromLLA(ecef, wgs84) == ENUfromECEF(ecef, wgs84) ∘ ECEFfromLLA(wgs84)

    @test UTMfromECEF(47, false, wgs84) == UTMfromLLA(47, false, wgs84) ∘ LLAfromECEF(wgs84)
    @test ECEFfromUTM(47, false, wgs84) == ECEFfromLLA(wgs84) ∘ LLAfromUTM(47, false, wgs84)

    @test UTMZfromECEF(wgs84) == UTMZfromLLA(wgs84) ∘ LLAfromECEF(wgs84)
    @test ECEFfromUTMZ(wgs84) == ECEFfromLLA(wgs84) ∘ LLAfromUTMZ(wgs84)

    @test ENUfromUTMZ(ecef, wgs84) == ENUfromLLA(ecef, wgs84) ∘ LLAfromUTMZ(wgs84)
    @test UTMZfromENU(ecef, wgs84) == UTMZfromLLA(wgs84) ∘ LLAfromENU(ecef, wgs84)

    @test ENUfromUTM(ecef, 47, false, wgs84) == ENUfromLLA(ecef, wgs84) ∘ LLAfromUTM(47, false, wgs84)
    @test UTMfromENU(ecef, 47, false, wgs84) == UTMfromLLA(47, false, wgs84) ∘ LLAfromENU(ecef, wgs84)
end