text:
module CommonToAll
using PyPlot, StatsBase, Statistics, Distances, LinearAlgebra,
      DelimitedFiles, ..AbstractOperator, NearestNeighbors, CSV, DataFrames

import ..Options, ..OptionsStat, ..OptionsNonstat, ..OptionsNuisance,
       ..history, ..GP.κ, ..calcfstar!, ..AbstractOperator.Sounding

export trimxft, assembleTat1, gettargtemps, checkns, getchi2forall, nicenup,
        plot_posterior, make1Dhist, make1Dhists, setupz, zcontinue, makezρ, plotdepthtransforms,
        unwrap, getn, geomprogdepth, assemblemodelsatT, getstats, gethimage,
        assemblenuisancesatT, makenuisancehists, stretchexists,
        makegrid, whichislast, makesummarygrid, makearray, plotNEWSlabels,
        plotprofile, gridpoints, splitsoundingsbyline, dfn2hdr

function trimxft(opt::Options, burninfrac::Float64, temperaturenum::Int)
    x_ft = assembleTat1(opt, :x_ftrain, burninfrac=burninfrac, temperaturenum=temperaturenum)
    n = assembleTat1(opt, :nodes, burninfrac=burninfrac, temperaturenum=temperaturenum)
    x, ft = zeros(size(opt.xall, 1), sum(n)), zeros(size(opt.fbounds, 1), sum(n))
    nlast = 0
    for (i, xft) in enumerate(x_ft)
        x[:,nlast+1:nlast+n[i]]  = xft[1:size(opt.xall, 1), 1:n[i]]
        ft[:,nlast+1:nlast+n[i]] = xft[size(opt.xall, 1)+1:end, 1:n[i]]
        nlast += n[i]
    end
    x, ft, n
end

function assembleTat1(optin::Options, stat::Symbol; burninfrac=0.5, temperaturenum=-1)
    @assert temperaturenum!=-1 "Please explicitly specify which temperature number"
    if stat == :fstar && temperaturenum!= 1
        return assemblemodelsatT(optin, burninfrac=burninfrac, temperaturenum=temperaturenum)
    end
    isns = checkns(optin)
    @assert 0.0<=burninfrac<=1.0
    Tacrosschains = gettargtemps(optin)
    temps = sort(unique(Tacrosschains))
    niters = size(Tacrosschains, 1)
    start = round(Int, niters*burninfrac)
    start == 0 && (start = 1)
    ttarget = temps[temperaturenum]
    nmodels = sum((Tacrosschains[start:end,:] .== ttarget))
    if stat == :fstar || stat == :x_ftrain
        mat1 = Array{Array{Float64}, 1}(undef, nmodels)
    else
        mat1 = Array{Real, 1}(undef, nmodels)
    end
    opt = deepcopy(optin)
    imodel = 0
    for ichain in 1:size(Tacrosschains, 2)
        opt.fstar_filename = "models_"*opt.fdataname*isns*"_$ichain.bin"
        opt.x_ftrain_filename = "points_"*opt.fdataname*isns*"_$ichain.bin"
        opt.costs_filename = "misfits_"*opt.fdataname*isns*"_$ichain.bin"
        if stat == :fstar
            at1idx = findall(Tacrosschains[:,ichain].==ttarget) .>= start
        else
            at1idx = Tacrosschains[:,ichain].==ttarget
            at1idx[1:start-1] .= false
        end
        ninchain = sum(at1idx)
        @info "chain $ichain has $ninchain models"
        ninchain == 0 && continue
        mat1[imodel+1:imodel+ninchain] .= history(opt, stat=stat)[at1idx]
        imodel += ninchain
    end
    iters = history(opt, stat=:iter)
    @info "obtained models $(iters[start]) to $(iters[end]) at T=$ttarget"
    mat1
end

function assemblemodelsatT(opt::OptionsStat; burninfrac=0.9, temperaturenum=-1)
    @assert temperaturenum!=-1 "Please explicitly specify which temperature number"
    x_ft = assembleTat1(opt, :x_ftrain, burninfrac=burninfrac, temperaturenum=temperaturenum)
    n = assembleTat1(opt, :nodes, burninfrac=burninfrac, temperaturenum=temperaturenum)
    nmodels = length(n)
    matT = Array{Array{Float64}, 1}(undef, nmodels)
    K_y = zeros(opt.nmax, opt.nmax)
    Kstar = zeros(Float64, size(opt.xall,2), opt.nmax)
    xtest = opt.xall
    for imodel = 1:nmodels
        xtrain = @view x_ft[imodel][1:size(opt.xall, 1), 1:n[imodel]]
        ftrain = x_ft[imodel][size(opt.xall, 1)+1:end, :]
        ky = view(K_y, 1:n[imodel], 1:n[imodel])
        map!(x->κ(opt.K, x),ky,pairwise(WeightedEuclidean(1 ./opt.λ² ), xtrain, dims=2))
        K_y[diagind(K_y)] .+= opt.δ^2
        ks = view(Kstar, :, 1:n[imodel])
        map!(x->κ(opt.K, x),Kstar,pairwise(WeightedEuclidean(1 ./opt.λ² ), xtest, xtrain, dims=2))
        matT[imodel] = zeros(size(opt.fbounds, 1), size(opt.xall, 2))
        fstar = matT[imodel]
        calcfstar!(fstar, ftrain, opt, K_y, Kstar, n[imodel])
    end
    matT
end

function assemblenuisancesatT(optn::OptionsNuisance;
    burninfrac = 0.5, temperaturenum = -1)
    @assert temperaturenum != -1 "Please specify temperature idx explicitly"
    @assert 0.0 <= burninfrac < 1.0
    Tacrosschains = gettargtemps(optn)
    #this is probably insanely inefficient
    #Θ(niters*nchains) to run the unique() call.
    #however, we need to iterate over the entire array
    #at some point to build the ensemble so it's
    #fine. Possible improvement: store temperatures
    #as the "temperature number" to start with and keep
    #a separate mapping from temperature index to the float
    #value. This means our sortedTs is just 1:nchains and
    #also avoids any issues with equality testing of floats.
    sortedTs = sort(unique(Tacrosschains))
    niters = size(Tacrosschains,1)
    #this will never give a bounds error
    #because of the assert above
    firsti = round(Int, niters*burninfrac)
    firsti == 0 && (start = 1)
    # firsti = 1 + floor(Int, niters*burninfrac)
    ttarg = sortedTs[temperaturenum]
    nmodels = sum(Tacrosschains[firsti:end,:] .== ttarg)
    fdataname = optn.fdataname
    #drop iteration number
    mvals = zeros(nmodels,optn.nnu)
    imodel = 0
    for ichain = 1:size(Tacrosschains,2)
        at1idx = Tacrosschains[:,ichain].==ttarg
        at1idx[1:firsti-1] .= false
        ninchain = sum(at1idx)
        @info "chain $ichain has $ninchain models"
        ninchain == 0 && continue
        vals_filename = "values_nuisance_"*fdataname*"$ichain.bin"
        ndat = readdlm(vals_filename, ' ', Float64)[:,2:end]
        mvals[imodel+1:imodel+ninchain,:] .= ndat[at1idx,:]
        imodel += ninchain
    end
    mvals
end

function getnchains(costs_filename)
c = 0
    for fname in readdir(pwd())
           r = Regex(costs_filename)
           c += !isa(match(r, fname), Nothing)
    end
    c
end

function gettargtemps(opt_in::Options)
    isns = checkns(opt_in)
    opt = deepcopy(opt_in)
    costs_filename = "misfits_"*opt.fdataname*isns
    nchains = getnchains(costs_filename)
    @info "Number of chains is $nchains"
    # now look at any chain to get how many iterations
    opt.costs_filename    = costs_filename*"_1.bin"
    iters          = history(opt, stat=:iter)
    niters         = length(iters)
    @info "McMC has run for $(iters[end]) iterations"
    # then create arrays of unsorted by temperature T
    Tacrosschains  = zeros(Float64, niters, nchains)
    # get the values into the arrays
    for ichain in 1:nchains
        opt.costs_filename = costs_filename*"_$ichain.bin"
        Tacrosschains[:,ichain] = history(opt, stat=:T)
    end
    Tacrosschains
end

function gettargtemps(optn_in::OptionsNuisance)
    optn = deepcopy(optn_in)
    costs_filename = "misfits_"*optn.fdataname*"nuisance"
    nchains = getnchains(costs_filename)
    @info "Number of chains is $nchains"
    optn.costs_filename = costs_filename*"_1.bin"
    iters = history(optn, stat=:iter)
    niters = length(iters)
    @info "MCMC has run for $(iters[end]) iterations"
    Tacrosschains = zeros(Float64, niters, nchains)
    for ichain in 1:nchains
        optn.costs_filename = costs_filename*"_$ichain.bin"
        Tacrosschains[:,ichain] = history(optn, stat=:T)
    end
    Tacrosschains
end

function getstats(optin::Options;
                  figsize=(5,6), fontsize=12,
                  nxticks=5, nyticks=5, alpha=0.6)
    isns = checkns(optin)
    opt = deepcopy(optin)
    costs_filename = "misfits_"*opt.fdataname*isns
    nchains = getnchains(costs_filename)
    chains = 1:nchains
    @info "Number of chains is $nchains"
    opt.costs_filename    = costs_filename*"_1.bin"
    iters          = history(opt, stat=:iter)
    statnames = [:acceptanceRateBirth, :acceptanceRateDeath,
                 :acceptanceRatePosition, :acceptanceRateProperty, :acceptanceRateDC]
    f,ax = plt.subplots(length(statnames), 1,
                        sharex=true, sharey=true, figsize=figsize)
    maxar = 0
    for (ichain, chain) in enumerate(chains)
        opt.costs_filename = costs_filename*"_$chain.bin"
        for (istat, statname) in enumerate(statnames)
            ar = history(opt, stat=statname)
            mx = maximum(ar[.!isnan.(ar)])
            mx > maxar && (maxar = mx)
            ax[istat].plot(iters, ar, alpha=alpha)
            if ichain == nchains
                ax[istat].set_title("$statname "*isns)
                ax[istat].set_ylabel("acc. %")
                if istat == length(statnames)
                    ax[istat].set_xlabel("mcmc step no.")
                    ax[istat].set_xticks(LinRange(iters[1], iters[end], nxticks))
                    ax[istat].set_xlim(extrema(iters))
                    ax[istat].set_yticks(LinRange(0, maxar, nyticks))
                    ax[istat].set_ylim(0, maxar)
                end
            end
        end
    end
    nicenup(f, fsize=fontsize)
end

function getstats(optin::OptionsNuisance;
                  figsize=(5,6), fontsize=12,
                  nxticks=5, nyticks=5, alpha=0.6)
    opt = deepcopy(optin)
    costs_filename = "misfits_"*opt.fdataname*"nuisance"
    nchains = getnchains(costs_filename)
    chains = 1:nchains
    @info "Number of chains is $nchains"
    opt.costs_filename = costs_filename*"_1.bin"
    iters          = history(opt, stat=:iter)
    statname = :acceptanceRate
    f,ax = plt.subplots(length(optin.idxnotzero), 1, sharex=true, sharey=true, figsize=figsize)
    maxar = 0
    for (ichain, chain) in enumerate(chains)
        opt.costs_filename = costs_filename*"_$chain.bin"
        ar = history(opt, stat=statname)[:,optin.idxnotzero]
        for (i, idx) in enumerate(optin.idxnotzero)
            mx = maximum(ar[.!isnan.(ar)])
            mx > maxar && (maxar = mx)
            ax[i].plot(iters, ar[:,i], alpha=alpha)
            if ichain == nchains
                ax[i].set_title("Nuisance #$idx")
                ax[i].set_ylabel("acc. %")
                if i == length(optin.idxnotzero)
                    ax[i].set_xlabel("mcmc step no.")
                    ax[i].set_xticks(LinRange(iters[1], iters[end], nxticks))
                    ax[i].set_xlim(extrema(iters))
                end
            end
        end
    end
    nicenup(f, fsize=fontsize)
end

function checkns(optin::Options)
    isns = typeof(optin) == OptionsNonstat
    @info "ns is $isns"
    ns = "ns"
    isns || (ns="s")
    return ns
end

function getchi2forall(optn_in::OptionsNuisance;
                       nchains = 1,
                       figsize = (6,4),
                       fsize = 8,
                       alpha = 0.25,
                       nxticks = 0,
                       gridon = false)
    optn = deepcopy(optn_in)
    fdataname = optn.fdataname
    costs_filename = "misfits_"*fdataname*"nuisance"
    if nchains == 1
        nchains = getnchains(costs_filename)
    end
    optn.costs_filename = costs_filename*"_1.bin"
    iters = history(optn, stat=:iter)
    niters = length(iters)
    Tacrosschains = zeros(Float64, niters, nchains)
    χ2acrosschains = zeros(Float64, niters, nchains)
    for chain in 1:nchains
        optn.costs_filename = costs_filename*"_$(chain).bin"
        Tacrosschains[:,chain] = history(optn, stat=:T)
        χ2acrosschains[:,chain] = history(optn, stat=:misfit)
    end
    Torder = sort([(i,j) for i=1:niters, j=1:nchains],
                    by = ix->Tacrosschains[ix...],
                    dims = 2)
    χ2sorted = [χ2acrosschains[Torder[i,j]...] for i=1:niters, j=1:nchains]
    fig, ax = subplots(3,1, sharex=true, figsize=figsize)
    ax[1].plot(iters, χ2acrosschains, alpha=alpha)
    ax[1].set_xlim(extrema(iters)...)
    ax[1].set_ylim(0,100)
    ax[1].set_title("unsorted χ^2 misfit")
    ax[2].plot(iters, Tacrosschains, alpha=alpha)
    ax[2].set_title("temperature")
    ax[3].plot(iters, χ2sorted, alpha=alpha)
    ax[3].set_ylim(0,100)
    ax[3].set_title("sorted χ^2 misfit")
    ax[3].set_xlabel("iterations")

    nxticks == 0 || ax[3].set_xticks(iters[1]:div(iters[end],nxticks):iters[end])
    nicenup(fig, fsize = fsize)
end

function getchi2forall(opt_in::Options;
                        nchains         = 1,
                        figsize         = (6,4),
                        fsize           = 8,
                        alpha           = 0.25,
                        nxticks         = 0,
                        gridon          = false,
                        omittemp        = false,
                      )
    # now look at any chain to get how many iterations
    isns = checkns(opt_in)
    opt = deepcopy(opt_in)
    costs_filename = "misfits_"*opt.fdataname*isns
    if nchains == 1 # then actually find out how many chains there are saved
        nchains = getnchains(costs_filename)
    end
    opt.costs_filename    = costs_filename*"_1.bin"
    iters          = history(opt, stat=:iter)
    niters         = length(iters)
    # then create arrays of unsorted by temperature T, k, and chi2
    Tacrosschains  = zeros(Float64, niters, nchains)
    kacrosschains  = zeros(Int, niters, nchains)
    X2by2inchains  = zeros(Float64, niters, nchains)
    # get the values into the arrays
    for ichain in 1:nchains
        opt.costs_filename = costs_filename*"_$ichain.bin"
        Tacrosschains[:,ichain] = history(opt, stat=:T)
        kacrosschains[:,ichain] = history(opt, stat=:nodes)
        X2by2inchains[:,ichain] = history(opt, stat=:U)
    end

    f, ax = plt.subplots(3,1, sharex=true, figsize=figsize)
    ax[1].plot(iters, kacrosschains, alpha=alpha)
    ax[1].set_xlim(extrema(iters)...)
    ax[1].set_title(isns*" unsorted by temperature")
    gridon && ax[1].grid()
    ax[1].set_ylabel("# nodes")
    ax[2].plot(iters, X2by2inchains, alpha=alpha)
    gridon && ax[2].grid()
    ax[2].set_ylabel("-Log L")
    gridon && ax[3].grid()
    ax[3].plot(iters, Tacrosschains, alpha=alpha)
    ax[3].set_ylabel("Temperature")
    ax[3].set_xlabel("iterations")
    nxticks == 0 || ax[3].set_xticks(iters[1]:div(iters[end],nxticks):iters[end])
    nicenup(f, fsize=fsize)

    Tunsorted = copy(Tacrosschains)
    for jstep = 1:niters
        sortidx = sortperm(vec(Tacrosschains[jstep,:]))
        X2by2inchains[jstep,:] = X2by2inchains[jstep,sortidx]
        kacrosschains[jstep,:] = kacrosschains[jstep,sortidx]
        Tacrosschains[jstep,:] = Tacrosschains[jstep,sortidx]
    end
    nrows = omittemp ? 2 : 3
    f, ax = plt.subplots(nrows, 1, sharex=true, figsize=figsize)
    nchainsatone = sum(Tacrosschains[1,:] .== 1)
    ax[1].plot(iters, kacrosschains, alpha=alpha)
    ax[1].set_xlim(extrema(iters)...)
    ax[1].set_ylabel("# nuclei")
    ax[1].set_title(isns*" sorted by temperature")
    ax[1].plot(iters, kacrosschains[:,1:nchainsatone], "k", alpha=alpha)
    gridon && ax[1].grid()
    ax[2].plot(iters, X2by2inchains, alpha=alpha)
    ax[2].plot(iters, X2by2inchains[:,1:nchainsatone], "k", alpha=alpha)
    ax[2].set_ylabel("-Log L")
    gridon && ax[2].grid()
    if !omittemp
        ax[3].plot(iters, Tunsorted, alpha=alpha, color="gray")
        ax[3].set_ylabel("Temperature")
        gridon && ax[3].grid()
    end
    nxticks == 0 || ax[3].set_xticks(iters[1]:div(iters[end],nxticks):iters[end])
    ax[nrows].set_xlabel("iterations")
    nicenup(f, fsize=fsize)

end

function geomprogdepth(n, dy, c)
    dy*(1.0-c^n)/(1-c)
end

function getn(z, dy, c)
    log(1 - z/dy*(1-c))/log(c)
end

function unwrap(v, inplace=false)
  # currently assuming an array
  unwrapped = inplace ? v : copy(v)
  for i in 2:length(v)
    while unwrapped[i] - unwrapped[i-1] >= pi
      unwrapped[i] -= 2pi
    end
    while unwrapped[i] - unwrapped[i-1] <= -pi
      unwrapped[i] += 2pi
    end
  end
  return unwrapped
end

function setupz(zstart, extendfrac;n=100, dz=10, showplot=false, forextension=false)
    @assert extendfrac>1
    znrange        = 1.0:n
    zboundaries    = [zstart, zstart .+ geomprogdepth.(znrange, dz, extendfrac)...]
    thickness      = diff(zboundaries)
    zall           = zboundaries[1:end-1] + thickness/2
    znall          = getn.(zall .- zstart, dz, extendfrac)
    showplot && plotdepthtransforms(zall, znall, zboundaries)
    zbreturn = forextension ? zboundaries : zboundaries[1:end-1]
    return zall, znall, zbreturn
end

function zcontinue(;zall=nothing, znall=nothing, zboundaries=nothing, n=nothing,
                extendfrac=nothing, dz=nothing, showplot=false)
    isa(dz, Nothing) && (dz=diff(zboundaries)[end])
    zall_, znall_, zboundaries_ = setupz(zboundaries[end], extendfrac, n=n, dz=dz, showplot=false)
    zall, znall, zboundaries = [zall;zall_], [znall; znall[end] .+ znall_], [zboundaries[1:end-1]; zboundaries_]
    zlast = zboundaries[end] + diff(zboundaries)[end]
    showplot && plotdepthtransforms(zall, znall, [zboundaries;zlast])
    @assert all(diff(zall).>0)
    @assert all(diff(znall).>0)
    @assert all(diff(zboundaries).>0)
    zall, znall, zboundaries
end

function setupz(zstart;n=100, dz=10)
    zboundaries    = range(zstart, step=dz, length=n)
    zall           = zboundaries[1:end-1] .+ dz/2
    znall          = 1.5:1:n-0.5
    plotdepthtransforms(zall, znall, zboundaries)
    return zall, znall, zboundaries[1:end-1]
end

function plotdepthtransforms(zall, znall, zboundaries)
    figure()
    plot(znall, zall)
    xlabel("depth index")
    ylabel("depth associated m")
    grid()
    nicenup(gcf())

    f, ax = plt.subplots(1, 1, figsize=(10,5))
    ax.stem(zboundaries[1:end-1], zboundaries[1:end-1], markerfmt="")
    ax.stem(zall, zall, "k--", markerfmt=" ")
    ax.set_xlabel("depth m")
    ax.set_ylabel("depth m")
    nicenup(gcf())

    thickness = diff(zboundaries)
    f, ax = plt.subplots(1, 2, sharey=true, figsize=(11,5))
    ax[1].stem(zall,thickness, "k--", markerfmt=" ")
    ax[1].set_ylabel("thickness m")
    ax[1].set_xlabel("depth m")
    ax[1].yaxis.grid(which="major")
    ax[2].stem(znall,thickness, "k--", markerfmt=" ")
    ax[2].set_xlabel("depth index")
    ax[2].yaxis.grid(which="major")
    nicenup(gcf())
end

function makezρ(zboundaries::Array{Float64, 1};
          zfixed      = [-1e6,    0.],
          ρfixed      = [1e12     0.3])
    @assert length(zfixed) == length(ρfixed)
    z = [zfixed..., zboundaries...]
    @assert all(diff(z).>0)
    ρ = zeros(length(z))
    nfixed = length(ρfixed)
    ρ[1:nfixed] .= ρfixed
    z, ρ, nfixed
end

function nicenup(g::PyPlot.Figure;fsize=16)
    for ax in gcf().axes
        ax.tick_params("both",labelsize=fsize)
        ax.xaxis.label.set_fontsize(fsize)
        ax.yaxis.label.set_fontsize(fsize)
        any(keys(ax) .== :zaxis) && ax.zaxis.label.set_fontsize(fsize)
        ax.title.set_fontsize(fsize)

        if typeof(ax.get_legend_handles_labels()[1]) != Array{Any,1}
            ax.legend(loc="best", fontsize=fsize)
        end
    end
    g.tight_layout()
end

function plot_posterior(F::Operator1D,
                        optns::OptionsNonstat,
                        opts::OptionsStat;
                        temperaturenum = 1,
                        nbins = 50,
                        burninfrac=0.5,
                        qp1=0.05,
                        qp2=0.95,
                        vmaxpc=1.0,
                        cmappdf = "inferno",
                        figsize=(12,5),
                        pdfnormalize=false,
                        fsize=14,
                        showlscale1vd=false,
                        istothepow=false,
                        doplot = true,
                        a = 0.25,
                        lw = 1)
    @assert 0<vmaxpc<=1
    M = assembleTat1(optns, :fstar, burninfrac=burninfrac, temperaturenum=temperaturenum)
    himage_ns, edges_ns, CI_ns, meanimage_ns,
    meandiffimage_ns, sdslope_ns, = gethimage(F, M, optns, temperaturenum=temperaturenum,
                            nbins=nbins, qp1=qp1, qp2=qp2, istothepow=istothepow,
                            islscale=false, pdfnormalize=pdfnormalize)
    M = assembleTat1(opts, :fstar, burninfrac=burninfrac, temperaturenum=temperaturenum)
    himage, edges, CI, meanimage,
    meandiffimage, sdslope, = gethimage(F, M, opts, temperaturenum=temperaturenum,
                            nbins=nbins, qp1=qp1, qp2=qp2, istothepow=false,
                            islscale=true, pdfnormalize=pdfnormalize)
    f,ax = plt.subplots(1, 3+convert(Int, showlscale1vd), sharey=true, figsize=figsize)
    xall = opts.xall
    diffs = diff(xall[:])
    xmesh = vcat(xall[1:end-1] - diffs/2, xall[end]-diffs[end]/2, xall[end])
    vmin, vmax = extrema(himage_ns)
    vmax = vmin+vmaxpc*(vmax-vmin)
    im1 = ax[1].pcolormesh(edges_ns[:], xmesh, himage_ns, cmap=cmappdf, vmax=vmax)
    cb1 = colorbar(im1, ax=ax[1])
    cb1.ax.set_xlabel("pdf \nns")
    @show propmin, propmax = getbounds(CI_ns, optns.fbounds)
    ax[1].set_xlim(propmin, propmax)
    ax[2].fill_betweenx(xall[:],meandiffimage_ns[:]-sdslope_ns[:],meandiffimage_ns[:]+sdslope_ns[:],alpha=a, color="gray")
    ax[2].plot(meandiffimage_ns, xall[:], linewidth=lw, color="k")
    ax[2].set_xlabel("slope")

    vmin, vmax = extrema(himage)
    vmax = vmin+vmaxpc*(vmax-vmin)
    im3 = ax[3].pcolormesh(edges[:], xmesh, himage, cmap=cmappdf, vmax=vmax)
    ax[3].set_ylim(extrema(xall)...)
    propmin, propmax = getbounds(CI, opts.fbounds)
    ax[3].set_xlim(propmin, propmax)
    ax[3].invert_yaxis()
    cb3 = colorbar(im3, ax=ax[3])
    cb3.ax.set_xlabel("pdf \nstationary")

    if showlscale1vd
        ax[4].plot(meandiffimage[:], xall[:], linewidth=2, color="g", alpha=0.5)
        ax[4].fill_betweenx(xall[:], meandiffimage[:]-sdslope[:],meandiffimage[:]+sdslope[:], alpha=.5)
        ax[4].set_xlabel("slope")
    end

    ax[1].plot(CI_ns, xall[:], linewidth=lw, color="w", alpha=a)
    # ax[1].plot(CI_ns, xall[:], linewidth=1, color="c", linestyle="--", alpha=0.5)
    # ax[1].plot(meanimage_ns[:], xall[:], linewidth=2, color="w", alpha=0.5)
    # ax[1].plot(meanimage_ns[:], xall[:], linewidth=2, color="k", linestyle="--", alpha=0.5)
    ax[3].plot(CI, xall[:], linewidth=lw, color="w", alpha=a)
    # ax[3].plot(CI, xall[:], linewidth=1, color="c", linestyle="--", alpha=0.5)
    # ax[3].plot(meanimage[:], xall[:], linewidth=2, color="w", alpha=0.5)
    # ax[3].plot(meanimage[:], xall[:], linewidth=2, color="k", linestyle="--", alpha=0.5)
    ax[1].set_xlabel(L"\log_{10} \rho")
    ax[1].set_ylabel("depth (m)")
    ax[3].set_xlabel(L"\log_{10} λ")
    nicenup(f, fsize=fsize)
end

function plot_posterior(F::Operator1D,
                        opt::OptionsStat;
                    temperaturenum = 1,
                    nbins = 50,
                    burninfrac=0.5,
                    qp1=0.05,
                    qp2=0.95,
                    vmaxpc=1.0,
                    cmappdf = "inferno",
                    figsize=(5,5),
                    pdfnormalize=false,
                    istothepow=false,
                    fsize=14,
                    plotCI = true,
                    plotmean = true,
                    alpha = 1.0,
                    CIcolor = ["w", "k"],
                    meancolor = ["m", "r"],
                    lwidth = 2,
                    doplot = true)
    @assert 0<vmaxpc<=1

    M = assembleTat1(opt, :fstar, burninfrac=burninfrac, temperaturenum=temperaturenum)
    himage, edges, CI, meanimage, meandiffimage, sdslope, = gethimage(F, M, opt, burninfrac=burninfrac, temperaturenum=temperaturenum,
                nbins=nbins, qp1=qp1, qp2=qp2, istothepow=istothepow,
                islscale=false, pdfnormalize=pdfnormalize)

    if doplot
        f,ax = plt.subplots(1,2, sharey=true, figsize=figsize)
        xall = opt.xall
        diffs = diff(xall[:])
        xmesh = vcat(xall[1:end-1] - diffs/2, xall[end]-diffs[end]/2, xall[end])
        vmin, vmax = extrema(himage)
        vmax = vmin+vmaxpc*(vmax-vmin)
        im1 = ax[1].pcolormesh(edges[:], xmesh, himage, cmap=cmappdf, vmax=vmax)
        ax[1].set_ylim(extrema(xall)...)
        ax[1].invert_yaxis()
        cb1 = colorbar(im1, ax=ax[1])
        cb1.ax.set_xlabel("pdf \nstationary")
        plotCI && ax[1].plot(CI, xall[:], linewidth=lwidth, color=CIcolor[1], alpha=alpha)
        plotCI && ax[1].plot(CI, xall[:], linewidth=lwidth, color=CIcolor[2], linestyle="--", alpha=alpha)
        plotmean && ax[1].plot(meanimage[:], xall[:], linewidth=lwidth, color=meancolor[1], alpha=alpha)
        plotmean && ax[1].plot(meanimage[:], xall[:], linewidth=lwidth, color=meancolor[2], linestyle="--", alpha=alpha)
        ax[1].set_xlabel(L"\log_{10} \rho")
        ax[1].set_ylabel("depth (m)")
        bounds = copy(opt.fbounds)
        if stretchexists(F)
            bounds = [minimum(F.low) maximum(F.low + F.Δ)]
        end
        istothepow && (bounds .= 10 .^ bounds)
        propmin, propmax = getbounds(CI, bounds)
        ax[1].set_xlim(propmin, propmax)
        ax[2].plot(meandiffimage[:], xall[:], linewidth=2, color="k", linestyle="-")
        zeroside = meandiffimage[:]-sdslope[:]
        zeroside[zeroside .< 0] .= 0
        ax[2].fill_betweenx(xall[:],zeroside,meandiffimage[:]+sdslope[:],alpha=.25)
        ax[2].set_xlabel("mean slope")
        nicenup(f, fsize=fsize)
    end
    CI[:,1], CI[:,2], CI[:,3], meanimage, meandiffimage, sdslope
end

function getbounds(CI, bounds)
    propmin = min(minimum(CI), minimum(bounds))
    propmax = max(maximum(CI), maximum(bounds))
    propmin, propmax
end

function plot_posterior(operator::Operator1D,
                        optn::OptionsNuisance;
                        temperaturenum = 1,
                        nbins = 50,
                        burninfrac=0.5,
                        figsize=(8,16),
                        pdfnormalize=false,
                        fsize=14,
                        qp1 = 0.05,
                        qp2 = 0.95,
                        doplot=true)
    hists, CI = makenuisancehists(optn, qp1, qp2, burninfrac = burninfrac,
                                     nbins = nbins, temperaturenum = temperaturenum,
                                     pdfnormalize = pdfnormalize)
    if doplot
        fig,ax = subplots(length(hists), 1, figsize=figsize)
        length(hists) == 1 && (ax = [ax])
        for (i, h) = enumerate(hists)
            bwidth = diff(h.edges[1])[1]
            bx = h.edges[1][1:end-1] .+ bwidth/2
            ax[i].bar(bx, h.weights, width=bwidth, edgecolor="black")
        end
        nicenup(fig, fsize=fsize)
    end
    hists, CI
end

function makenuisancehists(optn::OptionsNuisance, qp1, qp2; burninfrac = 0.5, nbins = 50,
    temperaturenum = -1, pdfnormalize = false)
    @assert temperaturenum != -1 "Please explicitly specify the temperature index"
    nuisanceatT = assemblenuisancesatT(optn, burninfrac = burninfrac,
                                    temperaturenum = temperaturenum)
    nnu = length(optn.idxnotzero)
    CI = zeros(Float64, nnu, 3)
    h = Vector{Histogram}(undef, nnu)
    for (i, islice) in enumerate(optn.idxnotzero)
        numin, numax = extrema(optn.bounds[islice,:])
        edges = LinRange(numin, numax, nbins+1)
        h[i] = fit(Histogram, nuisanceatT[:,islice], edges)
        CI[i,:] = quantile(nuisanceatT[:,islice], [qp1, 0.5, qp2])
    end
    h, CI
end

function firstderiv(x)
    fd = similar(x)
    fd[1]        =  x[2]     - x[1]
    fd[end]      =  x[end]   - x[end-1]
    fd[2:end-1] .= (x[3:end] - x[1:end-2])/2
    abs.(fd)
    # fd
end

function secondderiv(x)
    sd = similar(x)
    sd[1]        = 2x[1] - 5x[2] + 4x[3] - x[4]
    sd[2:end-1] .= (x[3:end] - x[2:end-1] + x[1:end-2])
    sd[end]      = 2x[end] - 5x[end-1] + 4x[end-2] - x[end-3]
    abs.(sd)
end

function gethimage(F::Operator, M::AbstractArray, opt::Options;
                burninfrac = 0.5,
                nbins = 50,
                rhomin=Inf,
                rhomax=-Inf,
                qp1=0.05,
                qp2=0.95,
                islscale = false,
                pdfnormalize=false,
                istothepow = false,
                temperaturenum=1)
    T = x->x
    if (rhomin == Inf) && (rhomax == -Inf)
        if !stretchexists(F) # if no stretch
            for (i,mm) in enumerate(M)
                rhomin_mm = minimum(mm)
                rhomax_mm = maximum(mm)
                rhomin_mm < rhomin && (rhomin = rhomin_mm)
                rhomax_mm > rhomax && (rhomax = rhomax_mm)
            end
            if (typeof(opt) == OptionsStat && opt.needλ²fromlog) && islscale
                T = x->0.5log10(x)
            end
        else # there is a stretch
            rhomin = minimum(F.low)
            rhomax = maximum(F.low + F.Δ)
        end
    else
        @assert rhomin < rhomax
    end
    if istothepow && !islscale
        T = x->10. ^x
    end
    rhomin, rhomax = map(x->T(x), (rhomin, rhomax))
    edges = LinRange(rhomin, rhomax, nbins+1)
    himage = zeros(Float64, length(M[1]), nbins)
    CI = zeros(Float64, length(M[1]), 3)
    meanimage = zeros(length(M[1]))
    for ilayer=1:length(M[1])
        if !stretchexists(F) # if no stretch
            mthislayer = [T(m[ilayer]) for m in M]
            himage[ilayer,:] = fit(Histogram, mthislayer, edges).weights
            CI[ilayer,:] = [quantile(mthislayer,(qp1, 0.5, qp2))...]
            meanimage[ilayer] = mean(vec(mthislayer))
        else # there is an affine stretch with depth
            mthislayer = [m[ilayer] for m in M]
            expandedthislayer = T.(F.low[ilayer] .+ mthislayer*F.Δ[ilayer])
            himage[ilayer,:] = fit(Histogram, expandedthislayer, edges).weights
            CI[ilayer,:] = [quantile(expandedthislayer,(qp1, 0.5, qp2))...]
            meanimage[ilayer] = mean(vec(expandedthislayer))
        end
        himage[ilayer,:] = himage[ilayer,:]/sum(himage[ilayer,:])/(diff(edges)[1])
        pdfnormalize && (himage[ilayer,:] = himage[ilayer,:]/maximum(himage[ilayer,:]))
    end
    if !stretchexists(F)
        Tm = [T.(m) for m in M]
    else
        Tm = [T.(F.low + vec(m).*F.Δ) for m in M]
    end
    Mslope = mean(firstderiv.(Tm))
    sdevslope = std(firstderiv.(Tm))
    himage, edges, CI, meanimage, Mslope, sdevslope
end

function stretchexists(F::Operator)
    in(:stretch, fieldnames(typeof(F))) && F.stretch
end

# plotting codes for 2D sections in AEM

function splitsoundingsbyline(soundings::Array{S, 1}) where S<:Sounding
    alllines = [s.linenum for s in soundings]
    linenum  = unique(alllines)
    nlines   = length(linenum)
    linestartidx = zeros(Int, nlines)
    for i = 1:nlines
        linestartidx[i] = findfirst(alllines .== linenum[i])
    end
    linestartidx
end

function makegrid(vals::AbstractArray, soundings::Array{S, 1};
    dr=10, zall=[NaN], dz=-1) where S<:Sounding
    @assert all(.!isnan.(zall))
    @assert dz>0
    X = [s.X for s in soundings]
    Y = [s.Y for s in soundings]
    x0, y0 = X[1], Y[1]
    R  = sqrt.((X .- x0).^2 + (Y .- y0).^2)
    rr, zz = [r for z in zall, r in R], [z for z in zall, r in R]
    topo = [s.Z for s in soundings]
    zz = topo' .- zz # mAHD
    kdtree = KDTree([rr[:]'; zz[:]'])
    gridr = range(R[1], R[end], step=dr)
    gridz = reverse(range(extrema(zz)..., step=dz))
    rr, zz = [r for z in gridz, r in gridr], [z for z in gridz, r in gridr]
    idxs, = nn(kdtree, [rr[:]'; zz[:]'])
    img = zeros(size(rr))
    for i = 1:length(img)
        img[i] = vals[idxs[i]]
    end
    topofine = gridpoints(R, gridr, topo)
    img[zz .>topofine'] .= NaN
    img, gridr, gridz, topofine, R
end

function makegrid(vals::AbstractArray, soundings::Array{S, 1}, XYZrho;
    dr=10, zall=[NaN], dz=-1) where S<:Sounding
    @assert all(.!isnan.(zall))
    @assert dz>0
    X = [s.X for s in soundings]
    Y = [s.Y for s in soundings]
    x0, y0 = X[1], Y[1]
    R  = sqrt.((X .- x0).^2 + (Y .- y0).^2)
    rr, zz = [r for z in zall, r in R], [z for z in zall, r in R]
    topo = [s.Z for s in soundings]
    zz = topo' .- zz # mAHD
    kdtree = KDTree([rr[:]'; zz[:]'])
    gridr = range(R[1], R[end], step=dr)
    gridz = reverse(range(extrema(zz)..., step=dz))
    rr, zz = [r for z in gridz, r in gridr], [z for z in gridz, r in gridr]
    idxs, = nn(kdtree, [rr[:]'; zz[:]'])
    # comparison with ZYZrho
    x0, y0 = XYZrho[1,1], XYZrho[1,2]
    rrcompare = sqrt.((XYZrho[:,1] .- x0).^2 + (XYZrho[:,2] .- y0).^2)
    kdtreecompare = KDTree([rrcompare[:]'; XYZrho[:,3]'])
    idxscompare, = nn(kdtreecompare, [rr[:]'; zz[:]'])
    img, imgcompare = zeros(size(rr)), zeros(size(rr))
    for i = 1:length(img)
        img[i] = vals[idxs[i]]
        imgcompare[i] = XYZrho[idxscompare[i], 4]
    end
    topofine = gridpoints(R, gridr, topo)
    img[zz .>topofine'] .= NaN
    imgcompare[zz .> topofine'] .= NaN
    img, gridr, gridz, topofine, R, imgcompare
end

function gridpoints(R, gridr, points::Array{Float64, 1})
    @assert length(R) == length(points)
    kdtree = KDTree(R[:]')
    idxs, = nn(kdtree, gridr[:]')
    pointsgrid = [points[idxs[i]] for i = 1:length(gridr)]
end

function gridpoints(R, gridr, points::Array{Float64, 2})
    nsound = length(R)
    nprops = size(points, 1)
    @assert nsound == size(points, 2)
    pointsgrid = zeros(length(gridr), nprops)
    for i = 1:nprops
        pointsgrid[:,i] = gridpoints(R, gridr, vec(points[i,:]))
    end
    pointsgrid
end

function whichislast(soundings::AbstractArray)
    X = [s.X for s in soundings]
    Y = [s.Y for s in soundings]
    Eislast, Nislast = true, true
    X[1]>X[end] && (Eislast = false)
    Y[1]>Y[end] && (Nislast = false)
    Eislast, Nislast
end

function makesummarygrid(soundings, pl, pm, ph, ρmean, vdmean, vddev, zall, dz; dr=10)
    # first flip ρ to σ and so pl and ph are interchanged
    phgrid, gridx, gridz, topofine, R = makegrid(-pl, soundings, zall=zall, dz=dz, dr=dr)
    plgrid,                           = makegrid(-ph, soundings, zall=zall, dz=dz, dr=dr)
    pmgrid,                           = makegrid(-pm, soundings, zall=zall, dz=dz, dr=dr)
    σmeangrid,                        = makegrid(-ρmean, soundings, zall=zall, dz=dz, dr=dr)
    ∇zmeangrid,                       = makegrid(vdmean, soundings, zall=zall, dz=dz, dr=dr)
    ∇zsdgrid,                         = makegrid(vddev, soundings, zall=zall, dz=dz, dr=dr)
    Z = [s.Z for s in soundings]
    phgrid, plgrid, pmgrid, σmeangrid, ∇zmeangrid, ∇zsdgrid, gridx, gridz, topofine, R, Z
end

function makexyzrho(soundings, zall)
    linename = "_line_$(soundings[1].linenum)_summary.txt"
    fnames = ["rho_low", "rho_mid", "rho_hi", "rho_avg"].*linename
    dlow, dmid, dhigh, davg = map(x->readdlm(x), fnames)
    ndepths = length(zall)
    for (i, d) in enumerate([dlow, dmid, dhigh, davg])
        @assert ndims(d) == 2
        @assert size(d, 2) == length(soundings)
        @assert size(d, 1) == ndepths
        xyzrho = makearray(soundings, d, zall)
        writedlm(fnames[i][1:end-4]*"_xyzrho.txt", xyzrho)
    end
end

function makearray(soundings, d, zall)
    outarray = zeros(length(d), 4)
    ndepths = length(zall)
        for (i, s) in enumerate(soundings)
            z = s.Z .- zall # convert depths to mAHD
            x, y = s.X, s.Y
            idx = (i-1)*ndepths+1:i*ndepths
            outarray[idx,:] = [[x y].*ones(ndepths) z d[:,i]]
        end
    outarray
end

function plotNEWSlabels(Eislast, Nislast, gridx, gridz, axarray)
    for s in axarray
        Eislast ? s.text(gridx[1], minimum(s.get_ylim()), "W", backgroundcolor="w") : s.text(gridx[1], minimum(s.get_ylim()), "E", backgroundcolor="w")
        Nislast ? s.text(gridx[end], minimum(s.get_ylim()), "N", backgroundcolor="w") : s.text(gridx[end], minimum(s.get_ylim()), "S", backgroundcolor="w")
    end
end

function plotprofile(ax, idxs, Z, R)
    for idx in idxs
        ax.plot(R[idx]*[1,1], [ax.get_ylim()[1], Z[idx]], "-w")
        ax.plot(R[idx]*[1,1], [ax.get_ylim()[1], Z[idx]], "--k")
    end
end

function selectwithin1Dinterval(M::AbstractVector, z,
                                zbounds, vbounds,
                                cond = :mean)
    @assert(any(cond .== [:all, :any, :median, :mean]))
    @assert size(vbounds, 1) == size(zbounds, 1)
    @assert length(z) == length(M[1])
    @assert all(vbounds[:,1] .< vbounds[:,2])
    @assert all(zbounds[:,1] .< zbounds[:,2])
    nconditions = size(vbounds, 1)
    idx = zeros(Bool, length(M), nconditions)
    for (i, m) in enumerate(M)
       for j in 1:nconditions
            idxdepth = zbounds[j,1] .< z .< zbounds[j,2]
            if any(cond .== [:median, :mean])
                if vbounds[j,1] < eval(cond)(m[idxdepth]) < vbounds[j,2]
                    idx[i,j] = true
                end
            else
                if eval(cond)(vbounds[j,1] .< (m[idxdepth]) .< vbounds[j,2])
                    idx[i,j] = true
                end
            end
       end
    end
    vec(reduce(&, idx, dims=2))
end

function block1Dvalues(M::AbstractVector, z, zbounds, cond = :median)
    @assert(any(cond .== [:median, :mean]))
    @assert length(z) == length(M[1])
    @assert all(zbounds[:,1] .< zbounds[:,2])
    nconditions = size(zbounds, 1)
    Mblock = zeros(length(M), nconditions)
    for (i, m) in enumerate(M)
        for j in 1:nconditions
            idxdepth = zbounds[j,1] .< z .<= zbounds[j,2]
            Mblock[i,j] = eval(cond)(m[idxdepth])
        end
     end
    Mblock
end

function correlationplot(M::Array{Float64, 2}; figsize=(5,5), nbins=25)
    f, ax = plt.subplots(size(M, 2), size(M, 2), figsize=figsize, sharex=true, sharey=true)
    nrows = size(M, 2)
    for j = 1:nrows
        for i = j:nrows
            h = normalize(fit(Histogram, (M[:,i], M[:,j]), nbins=nbins))
            ax[nrows*(j-1)+i].pcolormesh(h.edges[2], h.edges[1], h.weights)
        end
    end
    ax[1].set_aspect(1)
end

function plot1Dpatches(ax, zlist, xlist; alpha=0.25, fc="red", ec="blue", lw=2)
    @assert length(zlist) == length(xlist)
    for i in 1:size(zlist, 1)
        x0, z0 = xlist[i,1], zlist[i,1]
        delx = xlist[i,2] - x0
        delz = zlist[i,2] - z0
        ax.add_patch(matplotlib.patches.Rectangle((x0,z0),delx,delz, alpha=alpha, fc=fc, ec=fc))
        ax.add_patch(matplotlib.patches.Rectangle((x0,z0),delx,delz, fill=false, ec=ec, lw=lw, linestyle="--"))
    end
end

#filtering the projection added part from oasis
function filter_proj!(df)
    filter!(row -> !contains(row.combined, "RT=PROJ"), df)
end

#filtering the transformation added part from oasis
function filter_trans!(df)
    filter!(row -> !contains(row.combined,"RT=TRNS"), df)
end

#function to read the *dfn file and extract the column number and column names as a *.txt file
function dfn2hdr(dfnfile::String)
    df = CSV.read(dfnfile, DataFrame; header=false);
    df = coalesce.(df,"")

    colname = []
    for col in 1:length(names(df))
        push!(colname,names(df)[col])
    end
    t = df[!, colname[1]]
    for i in colname[2:end]
		t = t .* df[!, i]
	end
    df[!,"combined"] = t
    filter_proj!(df)
    filter_trans!(df)
    row_number = nrow(df)
    insertcols!(df, 1, :inc => 0)
    insertcols!(df, 1, :first => 0)
    insertcols!(df,1, :second => 0)

    cumulative_columns = 0  #this will set up a cumulative variable

    for i in 2:size(df,1)   #this will do operations in row level
        regex_literal1 = r"NAME=" #few parameters have their names after this keyword
        regex_literal2 = r"RT="  #four parameters have their names after this keyword
        if !contains(string(df[i,:].combined),regex_literal1)
            continue
        end

        if contains(df[i,:].combined,"RT:A4F")
            df[i,:].combined = replace(df[i,:].combined, ":A4F" => "AA")
        end

        inc_regex = r":([0-9]+)F"
        inc_match = match(inc_regex, df[i,:].combined)  #matching the keyword with the string
        if isnothing(inc_match)
            inc = 0
            idx2 =  idx2 = split(string(df[i,:].combined), regex_literal1)
            idx2_r = first.(split.(idx2[2], ":"))

        else
            inc = parse(Int64, inc_match.captures[1]) #the outcome is int64
            idx2 =  idx2 = split(string(df[i,:].combined), regex_literal2)
            idx2_r = first.(split.(idx2[2], ":"))

        end

        df[i,:].inc = inc
        df[i,:].first = cumulative_columns + 1
        df[i,:].second = df[i,:].first + inc
        cumulative_columns += inc + 1

        if occursin(";END DEFN", idx2_r)  #type is string
            idx2_r = first.(split.(idx2_r,";"))
        else
            idx2_r = first.(split.(idx2_r,"END"))
        end

        #start writing to a file here
        io = open("header.txt","a")
        if (df[i,:].first != df[i,:].second)
            writedlm(io,[df[i,:].first df[i,:].second idx2_r])
        else
            writedlm(io, [df[i,:].first idx2_r])
        end
        close(io)
    end
end

end # module CommonToAll


text:
struct MediaConfig
    series_path::String
    movies_path::String
end

struct Episode
    number::Union{Int,UnitRange}
    location::String
end

struct Season
    number::Int
    episodes::Vector{Episode}
end

struct TVShow
    name::String
    seasons::Vector{Season}
end

padded(x::Number) = lpad(x, 2, '0')

plex_denomination(ep::Episode) = join(["e"] .* padded.(collect(ep.number)), "-")
plex_denomination(s::Season, is_dir::Bool = false) = string(is_dir ? "Season " : "s", padded(s.number))
plex_denomination(s::Season, ep::Episode) = string(plex_denomination(s), plex_denomination(ep))
plex_denomination(show::TVShow, s::Season, ep::Episode) = string(show.name, " - ", plex_denomination(s, ep))

Base.:(==)(x::Season, y::Season) = x.number == y.number && length(x.episodes) == length(y.episodes) && all(x.episodes .== y.episodes)

Season(number::Int, location::AbstractString) = Season(number, location, 1)
Season(number::Int, locations::Vector{<:AbstractString}, args...) = Season(number, locations[number], args...)

function transfer(config::MediaConfig, show::TVShow; dest=nothing, force::Bool=false, confirm::Bool=true)
    dest = isnothing(dest) ? joinpath(config.series_path, show.name) : dest
    transfer(show, dest; force, confirm)
end

transfer(show::TVShow, dest::AbstractString; force::Bool=false, confirm::Bool=true) = map(x -> transfer(show, x, dest; force, confirm), show.seasons)

function transfer(src::AbstractString, dest::AbstractString; force::Bool=false)
    if force || src ≠ dest
        mv(src, dest; force)
    else
        @warn "Ignored existing destination $dest. Set 'force' to true to replace."
    end
end

function discover_seasons(name::AbstractString, src::AbstractString)
    name_parts = lowercase.(filter(x -> !isnothing(match(r"^\w+$", x)), split(basename(name), ['.', ' ', '_', '-'])))
    season_dirs = filter(x -> all(occursin(part, lowercase(x)) for part ∈ name_parts) && isdir(x), readdir(src, join=true))
    matches = filter(x -> !isnothing(x.first), map(x -> match(r"(?:Season\W*|Saison\W*|S|s)(\d{1,2})", x) => x, season_dirs))
    map(x -> Season(parse(Int, x.first.captures[1]), discover_episodes(x.second)), matches)
end

const ignored_extensions = [".jl", ".srt", ".jpeg", ".png", ".jpg", ".bmp"]

function discover_episodes(src::AbstractString)::Vector{Episode}
    episode_candidates = filter(!isnothing, find_media.(readdir(src, join=true)))
    matches = filter(x -> !isnothing(x.first), map(x -> something(match(r"\/(\d+)\s\-\s[\dA-Z]", x), match(r"(?:s|S)\d{2}(?:\-|x|e|E)(\d{2})", x), Some(nothing)) => x, episode_candidates))
    map(x -> Episode(parse(Int, x.first.captures[1]), x.second), matches)
end

function find_media(src::AbstractString)
    if isfile(src)
        last(splitext(src)) ∉ ignored_extensions && stat(src).size > 50e6 && return src # filesize more than 50 MB
        nothing
    elseif isdir(src)
        medias = filter(!isnothing, find_media.(readdir(src, join=true)))
        if length(medias) > 1
            @warn "Found more than one candidate while searching for an episode. Taking the bigger candidate.\nCandidates found: $medias"
        end
        if !isempty(medias)
            medias[argmax(stat.(medias))]
        else
            nothing
        end
    else
        nothing
    end
end

function transfer(show::TVShow, season::Season, dest::AbstractString; force::Bool=false, confirm::Bool=true)
    source_files = getproperty.(season.episodes, :location)
    pardirs = unique(dirname.(source_files))
    dest = joinpath(dest, plex_denomination(season, true))
    dest ∈ pardirs && error("Source and destination folders are the same for a non-copying transfer operation.")

    dests = map(x -> joinpath(dest, string(plex_denomination(show, season, x), last(splitext(x.location)))), season.episodes)

    printstyled(show.name, " ", bold=true, color=:yellow)
    printstyled(plex_denomination(season, true), bold=false, color=:yellow)
    println()
    show_transfer_recap.(source_files, dests)

    n = length(source_files)

    if n == 0
        @info "No files found for transfer."
        return
    end

    if length(pardirs) == 1
        @info "$n files ready for transfer from \"$(first(pardirs))\" to \"$dest\""
    else
        @info "$n files ready for transfer from multiple sources to \"$dest\""
    end

    if confirm
        @warn "The files will be moved to the new destination. Continue? (Y/n)"
        await_confirmation(true) || return
    end

    mkpath(dest)
    transfer.(source_files, dests; force)
    @info "Transfer successful."
end

show_transfer_recap(src::AbstractString, dst::AbstractString; padding=70) = begin
    printstyled(src, color=:red)
    printstyled(" ", rpad("=>", padding - length(src) - 1, '='), " ")
    printstyled(dst, color=:green)
    println()
end

function await_confirmation(default)
    try
        while true
            answer = strip(readline(), [' ', '\n'])
            if isempty(answer) && !isnothing(default)
                return default
            elseif lowercase(answer) ∉ ["y", "n"]
                println("Answer not understood. Transfer files? (y/n)")
            elseif answer == "y"
                return true
            else
                @info "Operation aborted."
                return false
            end
            default = nothing
        end
    catch e
        if e isa InterruptException
            @info "Operation aborted."
            exit(0)
        else
            rethrow(e)
        end
    end
end

ext: jl


text:
module SecondOrderSeawaterPolynomials

export
    SecondOrderSeawaterPolynomial,
    RoquetSeawaterPolynomial,
    RoquetEquationOfState

using SeawaterPolynomials: AbstractSeawaterPolynomial, BoussinesqEquationOfState

import SeawaterPolynomials: ρ′, thermal_expansion, haline_contraction

"""
    struct SecondOrderSeawaterPolynomial{FT} <: AbstractSeawaterPolynomial

Container of coefficients for a second-order polynomial function of
absolute salinity `Sᴬ`, conservative temperature `Θ`, and geopotential
depth `Z` for seawater density.

The coefficients have the form

```math
R_ᵦᵪᵩ ,
```

where ``β, χ, φ`` denote the order of the term to which the coefficent corresponds:
``β`` is the polynomial order of absolute salinity, ``Sᴬ``, ``χ`` is the polynomial order
of conservative temperature, ``Θ``, and ``φ`` is the order of geopotential height, ``Z``.

For a `SecondOrderSeawaterPolynomial`, ``β + χ + φ < 3``.

The coefficient ``R₁₁₀`` arises in the seawater polynomial as

```
seawater_polynomial(Θ, Sᴬ, Z) = ⋯ + R₁₁₀ * Sᴬ * Θ + ⋯
```
"""
@Base.kwdef struct SecondOrderSeawaterPolynomial{FT} <: AbstractSeawaterPolynomial
    R₁₀₀ :: FT = 0
    R₀₁₀ :: FT = 0
    R₁₀₁ :: FT = 0
    R₀₁₁ :: FT = 0
    R₁₁₀ :: FT = 0
    R₀₂₀ :: FT = 0
    R₂₀₀ :: FT = 0
end

const EOS₂ = BoussinesqEquationOfState{<:SecondOrderSeawaterPolynomial}

@inline ρ′(Θ, Sᴬ, Z, eos::EOS₂) = (  eos.seawater_polynomial.R₁₀₀ * Sᴬ
                                   + eos.seawater_polynomial.R₀₁₀ * Θ
                                   + eos.seawater_polynomial.R₀₂₀ * Θ^2
                                   - eos.seawater_polynomial.R₀₁₁ * Θ * Z
                                   + eos.seawater_polynomial.R₂₀₀ * Sᴬ^2
                                   - eos.seawater_polynomial.R₁₀₁ * Sᴬ * Z
                                   + eos.seawater_polynomial.R₁₁₀ * Sᴬ * Θ )

@inline thermal_expansion(Θ, Sᴬ, Z, eos::EOS₂) = (      eos.seawater_polynomial.R₀₁₀
                                                  + 2 * eos.seawater_polynomial.R₀₂₀ * Θ
                                                  -     eos.seawater_polynomial.R₀₁₁ * Z
                                                  +     eos.seawater_polynomial.R₁₁₀ * Sᴬ )

@inline haline_contraction(Θ, Sᴬ, Z, eos::EOS₂) = (      eos.seawater_polynomial.R₁₀₀
                                                   + 2 * eos.seawater_polynomial.R₂₀₀ * Sᴬ
                                                   -     eos.seawater_polynomial.R₁₀₁ * Z
                                                   +     eos.seawater_polynomial.R₁₁₀ * Θ )

"""
    RoquetSeawaterPolynomial([FT=Float64,] coefficient_set=:SecondOrder)

Returns a `SecondOrderSeawaterPolynomial` with coefficients optimized by

> Roquet et al., "Defining a Simplified yet 'Realistic' Equation of State for Seawater", Journal of Physical Oceanography (2015).

The `coefficient_set` is a symbol or string that selects one of the "sets" of
optimized second order coefficients.

Coefficient sets
================

    - :Linear: a linear equation of state, ρ = ρᵣ + R₁₀₀ * Θ + R₀₁₀ * Sᴬ

    - :Cabbeling: includes quadratic temperature term,
                  ρ = ρᵣ + R₁₀₀ * Θ + R₀₁₀ * Sᴬ + R₀₂₀ * Θ^2

    - :CabbelingThermobaricity: includes 'thermobaricity' term,
                                ρ = ρᵣ + R₁₀₀ * Θ + R₀₁₀ * Sᴬ + R₀₂₀ * Θ^2 - R₀₁₁ * Θ * Z

    - :Freezing: same as :cabbeling_thermobaricity with modified constants to increase
                   accuracy near freezing

    - :SecondOrder: includes quadratic salinity, halibaricity, and thermohaline term,
                    ρ = ρᵣ + R₁₀₀ * Θ + R₀₁₀ * Sᴬ + R₀₂₀ * Θ^2 - R₀₁₁ * T * Z
                           + R₂₀₀ * Sᴬ^2 - R₁₀₁ * Sᴬ * Z + R₁₁₀ * Sᴬ * Θ

The optimized coefficients are reported in Table 3 of Roquet et al., "Defining a Simplified
yet 'Realistic' Equation of State for Seawater", Journal of Physical Oceanography (2015), and
further discussed around equations (12)--(15). The optimization minimizes errors in estimated
horizontal density gradient estiamted from climatological temperature and salinity distributions
between the 5 simplified forms chosen by Roquet et. al and the full-fledged
[TEOS-10](http://www.teos-10.org) equation of state.
"""
RoquetSeawaterPolynomial(FT::DataType, coefficient_set=:SecondOrder) =
    eval(Symbol(coefficient_set, :RoquetSeawaterPolynomial))(FT)

# For Float64 default
RoquetSeawaterPolynomial(coefficient_set=:SecondOrder) =
    eval(Symbol(coefficient_set, :RoquetSeawaterPolynomial))(Float64)

"""
    RoquetEquationOfState([FT=Float64,] coefficient_set=:SecondOrder; reference_density=1024.6)

Returns an `BoussinesqEquationOfState` with a `RoquetSeawaterPolynomial` corresponding to
`coefficient_set` and with `reference density = 1024.6 kg m⁻³`, the average density of seawater
at the surface of the world ocean.

Type

```julia
help?> RoquetSeawaterPolynomial
```

for options for the `coefficient_set`. The optimzed coefficient sets for the
`RoquetSeawaterPolynomial` are listed in Table 3 in

> Roquet et al., "Defining a Simplified yet 'Realistic' Equation of State for Seawater", Journal of Physical Oceanography (2015).
"""
RoquetEquationOfState(FT::DataType, coefficient_set=:SecondOrder; reference_density=1024.6) =
    BoussinesqEquationOfState(RoquetSeawaterPolynomial(FT, coefficient_set), reference_density)

RoquetEquationOfState(coefficient_set=:SecondOrder; reference_density=1024.6) =
    RoquetEquationOfState(Float64, coefficient_set, reference_density=1024.6)

"""
    LinearRoquetSeawaterPolynomial(FT=Float64)

Parameters for a linear equation of state optimized for the 'current' oceanic
temperature and salinity distribution.

For more information, type `help?> RoquetSeawaterPolynomial`.
"""
LinearRoquetSeawaterPolynomial(FT=Float64) =
    SecondOrderSeawaterPolynomial{FT}(
                                      R₀₁₀ = - 1.775e-1,
                                      R₁₀₀ =   7.718e-1,
                                     )

"""
    CabbelingRoquetSeawaterPolynomial(FT=Float64)

Parameters for a minimal equation of state that describes cabbeling,
optimized for the 'current' oceanic temperature and salinity distribution.

For more information, type `help?> RoquetSeawaterPolynomial`.
"""
CabbelingRoquetSeawaterPolynomial(FT=Float64) =
    SecondOrderSeawaterPolynomial{FT}(
                                      R₀₁₀ = - 0.844e-1,
                                      R₁₀₀ =   7.718e-1,
                                      R₀₂₀ = - 4.561e-3,
                                     )

"""
    CabbelingThermobaricityRoquetSeawaterPolynomial(FT=Float64)

Parameters for a minimal equation of state that describes cabbeling and thermobaric
effects on sewater density, optimized for the 'current' oceanic temperature and salinity
distribution.

For more information, type `help?> RoquetSeawaterPolynomial`.
"""
CabbelingThermobaricityRoquetSeawaterPolynomial(FT=Float64) =
    SecondOrderSeawaterPolynomial{FT}(
                                      R₀₁₀ = - 0.651e-1,
                                      R₁₀₀ =   7.718e-1,
                                      R₀₂₀ = - 5.027e-3,
                                      R₀₁₁ = - 2.5681e-5,
                                     )

"""
    FreezingRoquetSeawaterPolynomial(FT=Float64)

Parameters for a minimal equation of state that describes seawater density near its
freezing point, optimized for the 'current' oceanic temperature and salinity
distribution.

For more information, type `help?> RoquetSeawaterPolynomial`.
"""
FreezingRoquetSeawaterPolynomial(FT=Float64) =
    SecondOrderSeawaterPolynomial{FT}(
                                      R₀₁₀ = - 0.491e-1,
                                      R₁₀₀ =   7.718e-1,
                                      R₀₂₀ = - 5.027e-3,
                                      R₀₁₁ = - 2.5681e-5,
                                     )

"""
    SecondOrderRoquetSeawaterPolynomial(FT=Float64)

Parameters for a fully second-order equation of state for seawater,
optimized for the 'current' oceanic temperature and salinity  distribution.

For more information, type `help?> RoquetSeawaterPolynomial`.
"""
SecondOrderRoquetSeawaterPolynomial(FT=Float64) =
    SecondOrderSeawaterPolynomial{FT}(
                                      R₀₁₀ =   0.182e-1,
                                      R₁₀₀ =   8.078e-1,
                                      R₀₂₀ = - 4.937e-3,
                                      R₀₁₁ = - 2.4677e-5,
                                      R₂₀₀ = - 1.115e-4,
                                      R₁₀₁ = - 8.241e-6,
                                      R₁₁₀ = - 2.446e-3
                                     )

end # module


text:
#=
   MOOMIN is a method for interpreting differential
   expression data, see publication for further details
   https://doi.org/10.1093/bioinformatics/btz584
   Author: T. Pusa
   Date:   06/2020
=#
module moomin

using MAT
using DelimitedFiles
using JuMP

include("data.jl")
include("weight.jl")
include("genes.jl")
include("reactions.jl")
include("model.jl")
include("MILP.jl")
include("output.jl")

"""
    runMoomin(pathToData, pathToModel, optimizer;
               delimiter='\\t', geneIDhead="GeneID", PPDEhead="PPDE", logFChead="logFC",
               modelStr="model",
               pThresh=0.9, alpha=3, precision=7,
               stoichiometry=true, enumerate=1,
               printLevel=1, timeLimit=1000)

   Main function to run the full algorithm. Input is a DE data file and metabolic
   network.

#  INPUTS

- `pathToData`:      A String that gives the location of file containing DE results
- `pathToModel`:     A String that gives the location of a .mat file containing
                      the metabolic network to be used
- `optimizer`:       An Optimizer object acting as the solver handle

# OPTIONAL INPUTS

- `delimiter`:      A Char that is used as the delimiter in `pathToData` (default: '\\t')
- `geneIDhead`:     A String that gives the header of the column in `pathToData`
                     that contains gene IDs (default: "GeneID")
- `PPDEhead`:       A String that gives the header of the column in `pathToData`
                     that contains PPDEs (default: "PPDE")
- `logFChead`:      A String that gives the header of the column in `pathToData`
                     that contains log fold changes (default: "logFC")
- `modelStr`:       A String that is the name of the model structure in `pathToModel`
                     Default usually works (default: "model")
- `pThresh`:        Double to be used as the threshold for a gene to be DE
                     (default: 0.9)
- `alpha`:          Double, parameter of the weight function, a higher value means less
                    evidence is needed for a change to be inferred (default: 3)
- `precision`:      Int to determine with how many significant digits reaction weights
                     are given
- `stoichiometry`:  Boolean to determine whether to use stoichiometric (true) or
                     simpler, topological constraints (false) (default: true)
- `enumerate`:      Int, how many alternative solutions are sought at most
- `printLevel`:     Int, determine how much progress info to print (0-2, default: 1)
- `timeLimit`:      Int, time limit for the MILP solver, in seconds (default: 1000)

"""
function runMoomin(pathToData, pathToModel, optimizer;
               delimiter='\t', geneIDhead="GeneID", PPDEhead="PPDE", logFChead="logFC",
               modelStr="model",
            	pThresh=0.9, alpha=3, precision=7,
               stoichiometry=true, enumerate=1,
               printLevel=1, timeLimit=1000)

   model = readModel(pathToData, pathToModel;
   	delimiter=delimiter, geneIDhead=geneIDhead, PPDEhead=PPDEhead, logFChead=logFChead,
   	modelStr=modelStr,
   	pThresh=pThresh, alpha=alpha, precision=precision)

   solve!(model, optimizer, stoichiometry=stoichiometry, enumerate=enumerate,
            printLevel=printLevel, timeLimit=timeLimit)

   return model
end

export Data, Genes, Reactions, MoominModel, nReactions, readData,
   mapDataToGenes!, weight, colour!,
   getAssGenes, readModel, createMILP, solve!, interpretSolution, runMoomin,
   fillOutputs!, writeOutput

end

max_stars_repo_name: htpusa/moominJL

max_issues_repo_name: htpusa/moominJL

max_forks_repo_name: htpusa/moominJL


text:
# CODATA 2018 NIST SP 961
const C0 = 299792458.0  # c, speed of light in vacuum, m/s
const U0 = 1.25663706212e-6  # μ₀, vacuum permeability, H/m
const E0 = 8.8541878128e-12  # ϵ₀, vacuum permittivity, F/m
const Z0 = 376.730313668  # Z₀, vacuum impedance, Ω

const QE = -1.602176634e-19  # qₑ, charge of an electron, C
const ME = 9.1093837015e-31  # mₑ, mass of an electron, kg

########

"""
    Species

Ionosphere constituent `Species`.

# Fields

- `charge::Float64`: signed species charged in Coulombs.
- `mass::Float64`: species mass in kilograms.
- `numberdensity`: a callable that returns number density in number per cubic meter as a
    function of height in meters.
- `collisionfrequency`: a callable that returns the collision frequency in collisions per
    second as a function of height in meters.

!!! note
    `numberdensity` and `collisionfrequency` will be converted to `FunctionerWrapper` types
    that tell the compiler that these functions will always return values of type `Float64`.
    A limited test is run to check if this is true, but otherwise it is up to the user
    to ensure these functions return only values of type `Float64`.
"""
struct Species
    charge::Float64  # C
    mass::Float64  # kg
    numberdensity::FunctionWrapper{Float64,Tuple{Float64}}  # m⁻³
    collisionfrequency::FunctionWrapper{Float64,Tuple{Float64}}  # s⁻¹

    function Species(charge, mass, numberdensity, collisionfrequency)
        _checkFloat64(numberdensity, 0:1e3:110e3)
        _checkFloat64(collisionfrequency, 0:1e3:110e3)
        new(charge, mass,
            FunctionWrapper{Float64,Tuple{Float64}}(numberdensity),
            FunctionWrapper{Float64,Tuple{Float64}}(collisionfrequency))
    end
end
Base.eachindex(s::Species) = 1
Base.length(s::Species) = 1
Base.getindex(s::Species, i) = i == 1 ? s : throw(BoundsError)

function _checkFloat64(f, x)
    for i in eachindex(x)
        v = f(x[i])
        v isa Float64 || throw(TypeError(Symbol(f), Float64, v))
    end
end

########

"""
    BField

Background magnetic field vector of strength `B` in Tesla with direction cosines `dcl`,
`dcm`, and `dcn` corresponding to ``x``, ``y``, and ``z`` directions parallel,
perpendicular, and up to the waveguide.
"""
struct BField
    B::Float64
    dcl::Float64
    dcm::Float64
    dcn::Float64

    function BField(B, dcl, dcm, dcn)
        if iszero(B)
            @warn "BField magnitude of exactly 0 is not supported. Setting B = 1e-16."
            B = 1e-16
        end
        new(B, dcl, dcm, dcn)
    end
end

"""
    BField(B, dip, azimuth)

Return a `BField` of strength `B` in Tesla, `dip` angle in radians from the horizontal, and
`azimuth` angle in radians from the propagation direction, positive towards ``y``.

Dip (inclination) is the angle made with the horizontal by the background magnetic field
lines. Positive dip corresponds to when the magnetic field vector is directed downward into
Earth. Therefore, it ranges between ``+π`` in Earth's geographic northern hemisphere and
``-π`` in the southern hemisphere.
"""
function BField(B, dip, azimuth)
    # BUG: How can we avoid this?
    # Look at the Booker quartic roots for dip angles -1,+1. They are way
    # outside the normal. The roots aren't accurately found, although I'm not
    # sure where the issue is.
    abs(dip) <= deg2rad(1) &&
        @warn "magnetic dip angles between ±1° have known numerical issues."
    abs(dip) > π &&
        @warn "magnetic dip angle should be in radians."
    abs(azimuth) > 2π &&
        @warn "magnetic azimuth angle should be in radians."

    Sdip, Cdip = sincos(dip)
    Saz, Caz = sincos(azimuth)
    BField(B, Cdip*Caz, Cdip*Saz, -Sdip)
end

"""
    dip(b::BField)

Return the dip angle in radians from `b`.
"""
dip(b::BField) = -asin(b.dcn)

"""
    azimuth(b::BField)

Return the azimuth angle in radians from `b`.
"""
azimuth(b::BField) = atan(b.dcm,b.dcl)

function isisotropic(b::BField)
    b.B <= 1e-12 && return true

    tolerance = deg2rad(0.15)

    bdip = dip(b)
    abs(bdip) < tolerance && return true

    baz = azimuth(b)
    (abs(baz - π/2) < tolerance || abs(baz - 3π/2) < tolerance) && return true

    return false
end

########

"""
    Ground(ϵᵣ, σ)

Isotropic earth ground characterized by a relative permittivity `ϵᵣ` and a conductivity `σ`
in Siemens per meter.
"""
struct Ground
    ϵᵣ::Int
    σ::Float64
end

"Default ground conductivity indices from LWPC."
const GROUND = Dict(
    1=>Ground(5, 1e-5),
    2=>Ground(5, 3e-5),
    3=>Ground(10, 1e-4),
    4=>Ground(10, 3e-4),
    5=>Ground(15, 1e-3),
    6=>Ground(15, 3e-3),
    7=>Ground(15, 1e-2),
    8=>Ground(15, 3e-2),
    9=>Ground(15, 1e-1),
    10=>Ground(81, 4.0))

########

@doc raw"""
    waitprofile(z, h′, β; cutoff_low=0, threshold=1e12)

Compute the electron number density in electrons per cubic meter at altitude `z` in meters
using Wait's exponential profile [[Wait1964]; [Thomson1993]] with parameters
`h′` in kilometers and `β` in inverse kilometers.

The profile is:
```math
Nₑ = 1.43 × 10¹³ \exp(-0.15 h') \exp[(β - 0.15)(z/1000 - h')]
```

Optional arguments:
- `cutoff_low=0`: when `z` is below `cutoff_low`, return zero.
- `threshold=1e12`: when density is greater than `threshold`, return `threshold`.

See also: [`electroncollisionfrequency`](@ref), [`ioncollisionfrequency`](@ref)

# References

[Wait1964]: J. R. Wait and K. P. Spies, “Characteristics of the earth-ionosphere waveguide
    for VLF radio waves,” U.S. National Bureau of Standards, Boulder, CO, Technical Note
    300, Dec. 1964.

[Thomson1993]: N. R. Thomson, “Experimental daytime VLF ionospheric parameters,” Journal of
    Atmospheric and Terrestrial Physics, vol. 55, no. 2, pp. 173–184, Feb. 1993,
    doi: 10.1016/0021-9169(93)90122-F.
"""
function waitprofile(z, hp, β; cutoff_low=0, threshold=1e12)
    if z > cutoff_low
        # Using form with single `exp` for speed
        Ne = 1.43e13*exp(-0.15*hp - (β - 0.15)*(hp - z*0.001))  # converting `z` to km
        if Ne > threshold
            return oftype(Ne, threshold)
        else
            return Ne
        end
    else
        zero(promote_type(Float64, typeof(z)))
    end
end

@doc raw"""
    electroncollisionfrequency(z)

Compute the electron-neutral collision frequency in collisions per second at height `z` in
meters based on Wait's conductivity profile [[Wait1964]; [Thomson1993]].

The profile is:
```math
νₑ(z) = 1.816 × 10¹¹ \exp(-0.15⋅z/1000)
```

See also: [`waitprofile`](@ref), [`ioncollisionfrequency`](@ref)

# References

[Wait1964]: J. R. Wait and K. P. Spies, “Characteristics of the earth-ionosphere waveguide
    for VLF radio waves,” U.S. National Bureau of Standards, Boulder, CO, Technical Note
    300, Dec. 1964.

[Thomson1993]: N. R. Thomson, “Experimental daytime VLF ionospheric parameters,” Journal of
    Atmospheric and Terrestrial Physics, vol. 55, no. 2, pp. 173–184, Feb. 1993,
    doi: 10.1016/0021-9169(93)90122-F.
"""
function electroncollisionfrequency(z)
    return 1.816e11*exp(-0.15e-3*z)  # e-3 converts `z` to km
end

@doc raw"""
    ioncollisionfrequency(z)

Compute the ion-neutral collision frequency in collisions per second at height `z` in meters
from [Morfitt1976].

The profile is:
```math
νᵢ(z) = 4.54 × 10⁹ \exp(-0.15⋅z/1000)
```

See also: [`waitprofile`](@ref), [`electroncollisionfrequency`](@ref)

# References

[Morfitt1976]: D. G. Morfitt and C. H. Shellman, “‘MODESRCH’, an improved computer program
    for obtaining ELF/VLF/LF mode constants in an Earth-ionosphere waveguide,” Naval
    Electronics Laboratory Center, San Diego, CA, NELC/IR-77T, Oct. 1976.
"""
function ioncollisionfrequency(z)
    return 4.54e9*exp(-0.15e-3*z)  # e-3 converts `z` to km
end

########
# Useful plasma functions

"""
    plasmafrequency(n, q, m)
    plasmafrequency(z, s::Species)

Return the plasma frequency ``ωₚ = √(nq²/(ϵ₀m))`` for a ``cold'' plasma.
"""
plasmafrequency(n, q, m) = sqrt(n*q^2/(E0*m))
plasmafrequency(z, s::Species) = plasmafrequency(s.numberdensity(z), s.charge, s.mass)

"""
    gyrofrequency(q, m, B)
    gyrofrequency(s::Species, b::BField)

Return the _signed_ gyrofrequency ``Ω = qB/m``.
"""
gyrofrequency(q, m, B) = q*B/m
gyrofrequency(s::Species, b::BField) = gyrofrequency(s.charge, s.mass, b.B)

"""
    magnetoionicparameters(z, frequency::Frequency, bfield::BField, species::Species)

Compute the magnetoionic parameters `X`, `Y`, and `Z` for height `z`.

```math
X = N e² / (ϵ₀ m ω²)
Y = e B / (m ω)
Z = ν / ω
```

# References

[Budden1955a]: K. G. Budden, “The numerical solution of differential equations governing
    reflexion of long radio waves from the ionosphere,” Proc. R. Soc. Lond. A, vol. 227,
    no. 1171, pp. 516–537, Feb. 1955, pp. 517.

[Budden1988]: K. G. Budden, “The propagation of radio waves: the theory of radio
    waves of low power in the ionosphere and magnetosphere,” First paperback edition.
    New York: Cambridge University Press, 1988, pp. 39.

[Ratcliffe1959]: J. A. Ratcliffe, "The magneto-ionic theory & its applications to the
    ionosphere," Cambridge University Press, 1959.
"""
function magnetoionicparameters(z, frequency::Frequency, bfield::BField, species::Species)
    m = species.mass
    ω = frequency.ω

    invω = inv(ω)
    invE0ω = invω/E0

    X, Y, Z = _magnetoionicparameters(z, invω, invE0ω, bfield, species)

    return X, Y, Z
end

# Specialized for performance of `susceptibility`
@inline function _magnetoionicparameters(z, invω, invE0ω, bfield, species)
    N, nu = species.numberdensity, species.collisionfrequency
    e, m = species.charge, species.mass
    B = bfield.B

    invmω = invω/m

    X = N(z)*e^2*invE0ω*invmω  # = N(z)*e^2/(E0*m*ω^2)
    Y = e*B*invmω  # = e*B/(m*ω)  # [Ratcliffe1959] pg. 182 specifies that sign of `e` is included here
    Z = nu(z)*invω  # = nu(z)/ω

    return X, Y, Z
end

"""
    waitsparameter(z, frequency::Frequency, bfield::BField, species::Species)

Compute Wait's conductivity parameter `Wr`.

```math
ωᵣ = Xω²/ν = ωₚ²/ν
```
"""
function waitsparameter(z, frequency::Frequency, bfield::BField, species::Species)
    X, Y, Z = magnetoionicparameters(z, frequency, bfield, species)

    ω = frequency.ω
    nu = species.collisionfrequency

    return X*ω^2/nu(z)
end


text:
"""
    RegisterAction(;
        file="$(contractuser(default_file("github", "workflows", "register.yml")))",
        destination="register.yml",
        prompt="Version to register or component to bump",
    )

Add a GitHub Actions workflow for registering a package with the General registry via workflow dispatch.
See [here](https://github.com/julia-actions/RegisterAction) for more information.

## Keyword Arguments
- `file::AbstractString`: Template file for the workflow file.
- `destination::AbstractString`: Destination of the workflow file,
  relative to `.github/workflows`.
- `prompt::AbstractString`: Prompt for workflow dispatch.
"""
@plugin struct RegisterAction <: FilePlugin
    file::String = default_file("github", "workflows", "register.yml")
    destination::String = "register.yml"
    prompt::String = "Version to register or component to bump"
end

source(p::RegisterAction) = p.file
destination(p::RegisterAction) = joinpath(".github", "workflows", p.destination)
tags(::RegisterAction) = "<<", ">>"

view(p::RegisterAction, ::Template, ::AbstractString) = Dict("PROMPT" => p.prompt)


text:
module Visualizing
  # using Graphs
  # using NetworkLayout
  # using GLMakie
  # using GraphMakie

  using ..Sources
  using ..DependencyTree
  using ..Converting

  function visualize(source::AbstractSource, file_title::String)
    sentences = Converting.convert(source) |> (text -> split(text, r"SENTENCE_END")) .|> strip

    for (index, sentence) in enumerate(sentences)
      filename = "$(file_title)_$(index).pdf"
      draw_sentence_graph(sentence, filename)
    end
  end

  function draw_sentence_graph(sentence, filename::String)
    g = Graphs.SimpleDiGraph()
    add_vertex!(g)
    nodelabel = ["ROOT-0"]
    edgelabel = Vector{String}()
    lines = split(sentence, "\n")

    foreach(lines) do line
      _, _, token = split(line)
      add_vertex!(g)
      push!(nodelabel, token)
    end

    for (index, line) in enumerate(lines)
      dependency_label, head_token, = split(line)
      head_index = findfirst(token -> head_token == token, nodelabel)

      push!(edgelabel, dependency_label)
      add_edge!(g, head_index , index + 1)
    end

    GraphMakie.graphplot(
      g,
      nlabels=nodelabel,
      nlabels_distance = 15,
      elabels = edgelabel,
      layout=NetworkLayout.Buchheim(),
      elabels_rotation=0,
      elabels_distance=25
    )
  end
end


text:
# This file is used for test purposes and not recommended for reading
using Knet
include("model.jl")
include("preprocess.jl")

function test_initparams(layerconfig, encembed, decembed, invocab, outvocab, winit, atype)
    # Desired encoder settings
    d_enc_embedding = (invocab, encembed)
    d_enc_weight = (encembed + layerconfig[1], 4layerconfig[1])
    d_enc_bias = (1, 4layerconfig[1])

    # Desired decoder settings
    d_dec_embedding = (outvocab, decembed)
    d_dec_weight = (decembed + layerconfig[1], 4layerconfig[1])
    d_dec_bias = (1, 4layerconfig[1])

    params = initparams(atype, layerconfig, encembed, decembed, invocab, outvocab, winit)
    # for a single lstm layer we need 8 parameters and each additional layer brings extra 4 parameters
    @assert length(params) == 4 + 4*length(layerconfig)
    @assert size(params[end]) == d_dec_embedding
    @assert size(params[end-1]) == d_enc_embedding
    @assert size(params[end-2]) == (layerconfig[end], outvocab)
    @assert size(params[end-3]) == (1, outvocab)

    @assert size(params[1]) == d_enc_weight
    @assert size(params[2]) == d_enc_bias

    @assert size(params[3]) ==d_dec_weight
    @assert size(params[4]) ==d_dec_bias

    info("Parameter initialization passed")
end

function test_initstates(atype, layerconfig, batchsize)
    # There has to be 4length(layerconfig) many states (half encoder half decoder)
    # The first hidden of the first layer of the encoder has to be 0
    len = length(layerconfig)
    states = initstates(atype, layerconfig, batchsize)
    @assert length(states) == 4len
    @assert states[1+2len] == nothing

    info("State initialization is passed")
end

function test_loss(parameters, states, inputseq, outputseq)
    loss(parameters, states, inputseq, outputseq)
    info("loss implementation is passed")
end

function test_lossgradient(parameters, states, inputseq, outputseq)
    N = 100 # check random gradients of the function paramaters
    gradcheck(loss, parameters, states, inputseq, outputseq;gcheck=10)
end

function test_loss_test(param, state, data; perp=true)
    totloss = 0.0
    numofbatch = 0
    for sequence in data
        val = loss(param, state, sequence, sequence)
        totloss += (perp? exp(val) : val)
        numofbatch += 1
    end
    return totloss / numofbatch
end

function update!(param, state, inputseq, outputseq; lr=1.0, gclip=0.0)
    gloss = lossgradient(param, state, inputseq, outputseq)
    gscale = lr

    for k=1:length(param)
        axpy!(-gscale, gloss[k], param[k])
    end

end

function train!(param, state, data)
    for sequence in data
        update!(param, state, sequence, sequence; lr=1, gclip=0)
    end
end

function tokenize1(vocabulary, prediction)
    index = findmax(prediction)[2]
    return filter(x->vocabulary[x] == index, collect(keys(vocabulary)))[1]
end

function test_train(params, states, data)
    inloss = test_loss_test(params, states, data; perp=true)
    println("initial loss is $inloss")
    for epoch=1:20
        train!(params, states, data)
        y = test_loss_test(params, states, data; perp=true)
        println("loss after epoch $epoch: $y")
    end
    info("training test worked")
end

function test_tokenize(data)
    correct = ["<s>", "The", "dog", "ran", "out", "of", "memory", "</s>"]
    indexes = Any[]
    for item in correct
        push!(indexes, data.word_to_index[item])
    end

    # generate fake numeric sequence
    t_sequence = Any[]
    for i=1:length(indexes)
        k = rand(1, length(data.word_to_index))
        k[indexes[i]] = 1
        push!(t_sequence, k)
    end

    # correct version of the sequence is known in advance for the test purposes
    for i=1:length(t_sequence)
        word = tokenize1(data.word_to_index, t_sequence[i])
        @assert(word == correct[i])
    end
    info("tokenize test is passed")
end

function test()
    # create simple test data
    batchsize = 3
    data = Data("readme_data.txt";batchsize=batchsize)
    seqs=Any[]; for item in data;push!(seqs, item);end;
    sequence = seqs[2] # <s> The dog ran out of memory </s>

    vsize = length(data.word_to_index)

    # test for single lstm layers
    global layerconfig = [64]
    encembed = 128
    decembed = 64
    invocab = vsize
    outvocab = vsize
    winit = 0.1
    atype = eval(parse("Array{Float32}"))

    # initialization test
    test_initparams(layerconfig, encembed, decembed, invocab, outvocab, winit, atype)
    test_initstates(atype, layerconfig, batchsize)

    # model initialization
    params = initparams(atype, layerconfig, encembed, decembed, invocab, outvocab, winit)
    states = initstates(atype, layerconfig, batchsize)

    # loss implementation test
    test_loss(params, states, sequence, sequence)
    #test_loss_test(params, states, data)

    # gradient calculation test
    # test_lossgradient(params, states, sequence, sequence)

    # tokenization test
    test_tokenize(data)
end
test()


text:
using Base.Test

#Better Julia practice:
#Write your tests as runnable scripts - you should be able to play all the file and know what went wrong
# write your own tests here
@test 42 == 42


text:
########################################
# Color utils
########################################
struct Color
    c::Union{Symbol, Int}
    isbold::Bool
end
Color() = Color(-1, false)
Color(c::Union{Symbol, Int}) = Color(c, false)
iscolored(c::Color) = c.c != -1 || c.isbold

function colorparse(s)
    i = prevind(s, ncodeunits(s) + 1)
    i1 = i
    i0 = i
    col = Color()
    modifiers = 0
    while i > 0
        c = s[i]
        if c == ':'
            if i0 <= i1
                cs = s[i0:i1]
                if cs == "bold"
                    col = Color(col.c, true)
                    modifiers += 1
                    i1 = prevind(s, i)
                else
                    intcode = tryparse(Int, cs)
                    code = intcode === nothing ? Symbol(cs) : intcode
                    if haskey(Base.text_colors, code) || code == :func
                        col = Color(code, col.isbold)
                        modifiers += 1
                        i1 = prevind(s, i)
                    end
                end
            end
        elseif c == '}' || c == '{'
            break
        end
        i0 = i
        i = prevind(s, i)
        modifiers == 2 && break
    end

    if iscolored(col) && col.c == -1
        col = Color(:normal, col.isbold)
    end

    return i1, col
end

########################################
# Token utils
########################################
struct Token
    val::String
    c::Color
end
iscolored(t::Token) = iscolored(t.c)
Base.isempty(t::Token) = isempty(t.val)

function Token(val)
    i, col = colorparse(val)
    Token(val[1:i], col)
end

########################################
# tokenizer and auxilary functions
########################################
function addtoken!(out, s, i0, i1)
    i0 > i1 && return
    i0 > ncodeunits(s) && return
    token = Token(s[i0:i1])
    isempty(token) && return
    push!(out, token)

    return
end

function recolor!(out, s, i0, i1, tind)
    i0 > ncodeunits(s) && return
    addtoken!(out, s, i0, i1)
    _, col = colorparse(s[i0:i1])
    if iscolored(col)
        for j in tind:length(out)
            iscolored(out[j]) || (out[j] = Token(out[j].val, col))
        end
    end

    return
end

function tokenize!(out, s, level = 1, i0 = 1)
    i00 = i0
    i = i0
    i1 = prevind(s, i0)
    tind = length(out) + 1
    while i <= ncodeunits(s)
        c = s[i]
        if c == '{'
            addtoken!(out, s, i0, i1)
            i = tokenize!(out, s, level + 1, nextind(s, i))
            i1 = prevind(s, i)
            i0 = i
            continue
        elseif c == '}'
            level == 1 && error("Unexpected } at position $i")
            recolor!(out, s, i0, i1, tind)
            return nextind(s, i)
        end
        i1 = i
        i = nextind(s, i)
    end
    level == 1 || error("Unbalanced { at position $i00")
    recolor!(out, s, i0, i1, tind)

    return i
end

tokenize(v::Vector{Token}) = v
function tokenize(s)
    out = Token[]
    tokenize!(out, s)

    return out
end


text:
using DifferentialEquations
using ProgressMeter
using PyCall
using MATLAB
using Printf

# pygui(:qt)
using PyPlot
include("PSS_Definitions_Dimer.jl")
# pygui(true);
using Roots
function Q1_find_dimer(Q2,P2,H_dimer)
    Q1_to_find(Q1)=Hamiltonian_Dimer([Q1,Q2],[0,P2],1)-H_dimer
    try
        Q1=find_zero(Q1_to_find,.01,maxeval=100,maxfnevals=300,tol=1e-15)
     catch
        Q1=zeros(0)
    end
end

H=range(.1,stop=0.25,length=14)

@showprogress 1 "Computing..." for Energy in H
location="/mnt/bdd38f66-9ece-451a-b915-952523c139d2/"
h=replace(@sprintf("%.15f",Energy),"."=>"_")
h_title=@sprintf("%.8f",Energy)
file_name=location*"DIMER_PSS"*h*".fig"
file_name=h*".fig"
mat"figure();set(gcf, 'Position',  [0, 0, 1500, 1500]); hold on;"
mat"xlabel('Q')"
mat"ylabel('P')"
mat"title($h_title)"
mat"axis([-3, 3, -4, 4])"
# Q1=Q1_find_dimer(0,.4,H)
# t_end=1000.0;
# q0,p0=[zeros(2) for i in 1:2]
# q0[1]=Q1;
# q0[2]=.2;
# p0[1]=0;
# p0[2]=.4;
# Hamiltonian_Dimer(q0,p0,1)
# prob= HamiltonianProblem{true}(Hamiltonian_Dimer, q0, p0, (0., t_end));
# @time t,A=solve(prob, Vern9(),reltol=1e-13,abstol=1e-15,maxiters=1e15);

# q1=A[1,:]
# q2=A[2,:]
# p1=A[3,:]
# p2=A[4,:];

# Q1=(q1+q2)/sqrt(2)
# Q2=(q1-q2)/sqrt(2)
# P1=(p1+p2)/sqrt(2)
# P2=(p1-p2)/sqrt(2);

# figure()
# plot(q1,p1,"b",linewidth=.1)
# plot(q2,p2,"r",linewidth=.1);

# figure()
# plot(Q1,P1,"b",linewidth=.1)
# plot(Q2,P2,"r",linewidth=.1);

t_end=1e5;
# figure()
# @time Q_PSS,P_PSS=PSS_function( q0[2], p0[2], H, t_end,count);
# plot(Q_PSS,P_PSS,".", markersize=4, c
#     ="k");

N_iter_Q=200;#50
Q_start=-2
Q_end=2
N_iter_P=200
P_start=-1.4
P_end=1.4
ArrP=range(P_start,stop=P_end,length=N_iter_P)
ArrQ=range(Q_start,stop=Q_end,length=N_iter_Q)
Brandons_Colors=["#393b79" ,"#5254a3","#6b6ecf","#9c9ede" ,"#637939","#8ca252" ,"#b5cf6b" ,"#cedb9c" ,"#8c6d31","#bd9e39" ,"#e7ba52","#e7cb94","#843c39","#ad494a" ,"#d6616b","#e7969c" ,"#7b4173" ,"#a55194","#ce6dbd" ,"#de9ed6"];

# figure()
@showprogress for k=1:N_iter_Q
    Q_PSS,P_PSS=PSS_function(ArrQ[k], 0, Energy, t_end);
    current_color=Brandons_Colors[mod(k,length(Brandons_Colors))+1]
    if Q_PSS!=0
        mat"plot($Q_PSS,$P_PSS,'.','MarkerSize',1,'color',$current_color); hold on;"
        # plot(Q_PSS,P_PSS,".", markersize=1,c=current_color);
    end
end
@showprogress for k=1:N_iter_P
    Q_PSS,P_PSS=PSS_function(0,ArrP[k], Energy, t_end);
    current_color=Brandons_Colors[mod(k,length(Brandons_Colors))+1]
    if Q_PSS!=0
        mat"plot($Q_PSS,$P_PSS,'.','MarkerSize',1,'color',$current_color); hold on;"
        #  plot(Q_PSS,P_PSS,".", markersize=1,c=current_color);
    end
end

mat"savefig($file_name)"

end


text:
using FdeSolver
using Plots

## inputs
tSpan = [0, 25]        # [intial time, final time]
y0 = [34, 6]           # initial values
β = [0.98, 0.99]       # order of derivatives

## ODE model

par = [0.55, 0.028, 0.84, 0.026] # model parameters

function F(t, y, par)

    α1 = par[1]      # growth rate of the prey population
    β1 = par[2]      # rate of shrinkage relative to the product of the population sizes
    γ = par[3]       # shrinkage rate of the predator population
    δ = par[4]       # growth rate of the predator population as a factor of the product
                     # of the population sizes

    u = y[1]         # population size of the prey species at time t[n]
    v = y[2]         # population size of the predator species at time t[n]

    F1 = α1 .* u .- β1 .* u .* v
    F2 = - γ .* v .+ δ .* u .* v

    [F1, F2]

end

## Solution
t, Yapp = FDEsolver(F, tSpan, y0, β, par)

# plotting
plot(t, Yapp, linewidth = 5, title = "Solution to LV model with 2 FDEs",
     xaxis = "Time (t)", yaxis = "y(t)", label = ["Prey" "Predator"])
plot!(legendtitle = "Population of")


text:
using Portfolio
using Test

@testset "Portfolio.jl" begin

    @test

end


text:
# Autogenerated wrapper script for SISL_jll for armv6l-linux-musleabihf
export libsisl

JLLWrappers.@generate_wrapper_header("SISL")
JLLWrappers.@declare_library_product(libsisl, "libsisl.so")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libsisl,
        "lib/libsisl.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
abstract AtomicUnit{T}

atoms = [
    #     Name   Symbol  Valence
    (:Hydrogen,  :H,     1),
    (  :Carbon,  :C,     4)
]

for (atom, sym, valence) in atoms
    @eval begin
        type $atom{T<:Number} <: AtomicUnit{T}
            value::Complex{T}    # atomic value
            valence::Int         # number of valence electrons
            degree::Int          # number of shared valence electons
            free::Int            # number of free valence electrons

            $atom(value) = new(value, $valence, 0, $valence)
        end

        $atom{T<:Real}(value::Complex{T}) = $atom{T}(value)
        $atom{T<:Real}(value::T) = $atom(complex(value))
        $atom() = $atom(0.0)

        typealias $sym $atom
    end
end

function Base.show(io::IO, a::AtomicUnit)
    print(io, "$(typeof(a))(value = $(a.value), valence = $(a.valence), degree = $(a.degree), free = $(a.free))")
end

Base.(:(==))(a₁::AtomicUnit, a₂::AtomicUnit) = a₁.valence == a₂.valence


text:
default_wrapper_specs = Dict{String, Vector{Pair{Type{T} where T <: AbstractGymWrapper, Tuple}}}(
    "Pong"=>[
        ActionSetWrapper=>([0, 2, 3],),
        RageQuitWrapper=>(7.0,),
                   # reward_thresh, reward_max, min_acceptable_ent_ratio, ignore_actions
        ActionEntropyWrapper=>(1.0, 0.5, 0.2, []),
        SurvivalRewardWrapper=>(2.0, 1.0, 99, 500)
    ],
    "Breakout"=>[
        ActionSetWrapper=>([0, 2, 3],),
                   # reward_thresh, reward_max, min_acceptable_ent_ratio, ignore_actions
        ActionEntropyWrapper=>(1.0, 0.5, 0.2, [1]),
    ],
    "deepmind_defaults"=>[
        EpisodicLifeWrapper=>(),
        RageQuitWrapper=>(7,),
        GreyMeanWrapper=>(),
        # GreyChanWrapper=>(2,),
        DownsizeWrapper=>(84, 84)
        MultiFrameWrapper=>(4,), # stack last 4 frames
        MaxAndSkipWrapper=>(4,), # repeat action for 4 frames
    ]
)


text:
@testset "regression" begin
    y = randn(10)
    l = GaussianLikelihood()
    @test AGP.treat_labels!(y, l) == y
    @test predict_y(l, y) == y
    @test predict_y(l, (y,)) == y
end


text:
# Test escaping of $
@test "query(\$var:ID!){country(var:\$var){name}}" == gql"query($var:ID!){country(var:$var){name}}"

# Test turning off validation errors
str = @gql_str """
{countries{name}}
query{countries{name}}
""" false
@test_throws GraphQLParser.ValidationException GraphQLParser.is_valid_executable_document(str; throw_on_error=true)


text:
using BenchmarkTools

using Jecco
using Jecco.AdS5_3_1

ABCS  = zeros(4)
AA    = zeros(2,2)
BB    = zeros(2,2)
CC    = zeros(2,2)
SS    = zeros(2)
abccS = 0

svars    = tuple(ones(11)...)
fvars    = tuple(ones(37)...)
sdvars   = tuple([ZeroPotential(); ones(70)]...)
bdgvars  = tuple(ones(71)...)
phidvars = tuple([ZeroPotential(); ones(72)]...)
avars    = tuple([ZeroPotential(); ones(76)]...)
xivars   = tuple(ones(84)...)

# inner grid

println("S equation inner grid coefficients:")
@btime Jecco.AdS5_3_1.S_eq_coeff!($ABCS, $svars, AdS5_3_1.Inner())

println("Fxy equations inner grid coefficients:")
@btime Jecco.AdS5_3_1.Fxy_eq_coeff!($AA, $BB, $CC, $SS, $fvars, AdS5_3_1.Inner())

println("Sd equation inner grid coefficients:")
@btime Jecco.AdS5_3_1.Sd_eq_coeff!($ABCS, $sdvars, AdS5_3_1.Inner())

println("B2d equation inner grid coefficients:")
@btime Jecco.AdS5_3_1.B2d_eq_coeff!($ABCS, $bdgvars, AdS5_3_1.Inner())

println("B1dGd equations inner grid coefficients:")
@btime Jecco.AdS5_3_1.B1dGd_eq_coeff!($AA, $BB, $CC, $SS, $bdgvars, AdS5_3_1.Inner())

println("phid equations inner grid coefficients:")
@btime Jecco.AdS5_3_1.phid_eq_coeff!($ABCS, $phidvars, AdS5_3_1.Inner())

println("A equations inner grid coefficients:")
@btime Jecco.AdS5_3_1.A_eq_coeff!($ABCS, $avars, AdS5_3_1.Inner())

# outer grid

println("S equation outer grid coefficients:")
@btime Jecco.AdS5_3_1.S_eq_coeff!($ABCS, $svars, AdS5_3_1.Outer())

println("Fxy equations outer grid coefficients:")
@btime Jecco.AdS5_3_1.Fxy_eq_coeff!($AA, $BB, $CC, $SS, $fvars, AdS5_3_1.Outer())

println("Sd equation outer grid coefficients:")
@btime Jecco.AdS5_3_1.Sd_eq_coeff!($ABCS, $sdvars, AdS5_3_1.Outer())

println("B2d equation outer grid coefficients:")
@btime Jecco.AdS5_3_1.B2d_eq_coeff!($ABCS, $bdgvars, AdS5_3_1.Outer())

println("B1dGd equations outer grid coefficients:")
@btime Jecco.AdS5_3_1.B1dGd_eq_coeff!($AA, $BB, $CC, $SS, $bdgvars, AdS5_3_1.Outer())

println("phid equations outer grid coefficients:")
@btime Jecco.AdS5_3_1.phid_eq_coeff!($ABCS, $phidvars, AdS5_3_1.Outer())

println("A equations outer grid coefficients:")
@btime Jecco.AdS5_3_1.A_eq_coeff!($ABCS, $avars, AdS5_3_1.Outer())

#xi equation

println("xi_t equations coefficients:")
@btime $abccS = Jecco.AdS5_3_1.xi_t_eq_coeff($xivars, AdS5_3_1.Outer())


text:
using Documenter, InfrastructureSystems

if haskey(ENV, "GITHUB_ACTIONS")
    ENV["JULIA_DEBUG"] = "Documenter"
end

makedocs(
    modules = [InfrastructureSystems],
    format = Documenter.HTML(prettyurls = haskey(ENV, "GITHUB_ACTIONS"),),
    sitename = "InfrastructureSystems.jl",
    pages = Any[ # Compat: `Any` for 0.4 compat
        "Home" => "index.md",
        "User Guide" => "man/guide.md",
        "API" => Any[
            "InfrastructureSystems" => "api/InfrastructureSystems.md"
        ]
    ]
)

deploydocs(
    repo = "github.com/NREL-SIIP/InfrastructureSystems.jl.git",
    target = "build",
    branch = "gh-pages",
    devbranch = "master",
    devurl = "dev",
    versions = ["stable" => "v^", "v#.#"],
    push_preview = true,
)


text:
"""
    Cnk(n, k, σ, prec)

Calculation of the Cnk using the direct formula with arbitrary precision, using a specified number of bits (prec). This is unstable for large values of k, more so for small values of σ.

# Examples
```julia-repl
julia> GibbsTypePriors.Cnk(6, 5, 0.5, 200)
[0.234375000000000000000000000000000000000000000000000000000000 +/- 5.85e-61]
```
"""
function Cnk(n, k, σ, prec)
    RF = RealField(prec)
    σ_arb = RF(σ)
    return 1//fac(RF(k)) * sum([ (-1)^i * binom(k, i, RF) * risingfac(-i * σ_arb,n) for i in 0:k])
end

function noncentral_generalised_factorial_coefficient_prec(n, k, s, r, prec)::arb
    RF = RealField(prec)
    arb_1 = RF(1)
    arb_0 = RF(0)
    arb_r = RF(r)
    arb_s = RF(s)
    return noncentral_generalised_factorial_coefficient(n, k, arb_s, arb_r; arb_0=arb_0, arb_1=arb_1)
end

function Cnk_rec(n, k, σ, prec)
    # factor_k = factorial(k)
#   sum((-1)^(1:k) * gmp::chooseZ(n = k, k = 1:k) * Rmpfr::pochMpfr(-(1:k) * Gama, n) / factor_k)
  (-1)^(n - k) * noncentral_generalised_factorial_coefficient_prec(n, k, σ, 0, prec)
end

max_stars_repo_name: konkam/GibbsTypePriors

max_issues_repo_name: konkam/GibbsTypePriors

max_forks_repo_name: konkam/GibbsTypePriors


text:
#!/Applications/Julia-1.4.app/Contents/Resources/julia/bin/julia

function helloworld()
    println("Hello, World!") # Bye bye, MATLAB!
end

helloworld()

max_stars_repo_name: faicaltoubali/ENSEEIHT

max_issues_repo_name: faicaltoubali/ENSEEIHT

max_forks_repo_name: faicaltoubali/ENSEEIHT


text:
if !isdefined(:__EXPRESSION_HASHES__)
    __EXPRESSION_HASHES__ = Set{UInt64}()
end

macro runonce(expr)
    h = hash(expr)
    return esc(quote
        if !in($h, __EXPRESSION_HASHES__)
            push!(__EXPRESSION_HASHES__, $h)
            $expr
        end
    end)
end

function ccv_type{T}(::Type{T})
    if haskey(JULIA2CCV_TYPE, T)
        return JULIA2CCV_TYPE[T]
    else
        throw(ArgumentError("CCV can't handle Julia's $T type"))
    end
end

function julia_type(ccv_typ::Cuint)
    if haskey(CCV2JULIA_TYPE, ccv_typ)
        return CCV2JULIA_TYPE[ccv_typ]
    else
        throw(ArgumentError("Unknown CCV type with cide $ccv_typ"))
    end
end


text:
function ObservationData(arg0::ObservationType, arg1::jdouble, arg2::jint, arg3::jint)
    return ObservationData((ObservationType, jdouble, jint, jint), arg0, arg1, arg2, arg3)
end

function get_loss_of_lock_indicator(obj::ObservationData)
    return jcall(obj, "getLossOfLockIndicator", jint, ())
end

function get_observation_type(obj::ObservationData)
    return jcall(obj, "getObservationType", ObservationType, ())
end

function get_signal_strength(obj::ObservationData)
    return jcall(obj, "getSignalStrength", jint, ())
end

function get_value(obj::ObservationData)
    return jcall(obj, "getValue", jdouble, ())
end


text:
r = 10;
A = pi*r^2; # pi has a default value!
println(A);

max_stars_repo_name: fahimfarhan/legendary-coding-odyssey

max_issues_repo_name: fahimfarhan/legendary-coding-odyssey

max_forks_repo_name: fahimfarhan/legendary-coding-odyssey


text:
module TestCodeInfoTools

using CodeInfoTools
using Test

include("canvas.jl")
include("builder.jl")
include("evaluation.jl")
include("misc.jl")

end # module


text:
module Uprating

#
# This module uprates the main model data using values from CPI, Nominal GDP, average wages, etc. that
# I've hacked together from SFC,UPRATING and Bank Of England Data.
# See the worksheets under `data/prices/`.
#
using DataFrames
using CSV

using ScottishTaxBenefitModel
using .RunSettings
using .Definitions
using .Utils

"""
Semi-complete indexing routine using UPRATING quarterly data.
"""

export uprate, UPRATE_MAPPINGS

Uprate_Map = Dict(
    upr_earnings => :average_earnings,
    upr_housing_rents => :actual_rents_for_housing,
    upr_housing_oo => :mortgage_interest_payments,
    upr_unearned => :nominal_gdp,
    upr_costs => :cpi,
    upr_cpi => :cpi,
    # not used upr_gdp_deflator => :gdp_deflator,
    upr_nominal_gdp => :nominal_gdp,
    upr_shares => :equity_prices
)

function make_uprate_types() :: Dict
    d = Dict()
    for i in instances( Incomes_Type )
        if i <= odd_jobs || i in [trade_unions_etc,friendly_societies,work_expenses,avcs,other_deductions]
            d[i] = upr_earnings
        elseif i in [
            private_pensions,
            national_savings,
            bank_interest,
            individual_savings_account,
            # dividends,
            property,
            royalties,
            bonds_and_gilts,
            other_investment_income,
            other_income,
            pension_contributions_employee,
            pension_contributions_employer,
            loan_repayments ]
            d[i] = upr_nominal_gdp
        elseif i in [stocks_shares]
            d[i] = upr_shares
        elseif i in [
            alimony_and_child_support_received,
            health_insurance,
            alimony_and_child_support_paid,
            care_insurance,
            student_loan_repayments,

            education_allowances,
            foster_care_payments,
            student_grants,
            student_loans,
            free_school_meals ]
            d[i] = upr_cpi
        else # keep bens as they are: freeze, plus we just need that a payment has been made
            d[i] = upr_no_uprate
        end
    end

    d[A_Current_account] = upr_nominal_gdp
    d[A_NSB_Ordinary_account] = upr_nominal_gdp
    d[A_NSB_Investment_account] = upr_nominal_gdp
    d[A_Not_Used] = upr_nominal_gdp
    d[A_Savings_investments_etc] = upr_nominal_gdp
    d[A_Government_Gilt_Edged_Stock] = upr_nominal_gdp
    d[A_Unit_or_Investment_Trusts] = upr_shares
    d[A_Stocks_Shares_Bonds_etc] = upr_shares
    d[A_PEP] = upr_nominal_gdp
    d[A_National_Savings_capital_bonds] = upr_nominal_gdp
    d[A_Index_Linked_National_Savings_Certificates] = upr_nominal_gdp
    d[A_Fixed_Interest_National_Savings_Certificates] = upr_nominal_gdp
    d[A_Pensioners_Guaranteed_Bonds] = upr_nominal_gdp
    d[A_SAYE] = upr_nominal_gdp
    d[A_Premium_bonds] = upr_nominal_gdp
    d[A_National_Savings_income_bonds] = upr_nominal_gdp
    d[A_National_Savings_deposit_bonds] = upr_nominal_gdp
    d[A_First_Option_bonds] = upr_nominal_gdp
    d[A_Yearly_Plan] = upr_nominal_gdp
    d[A_ISA] = upr_nominal_gdp
    d[A_Fixd_Rate_Svngs_Bonds_or_Grntd_Incm_Bonds_or_Grntd_Growth_Bonds] = upr_nominal_gdp
    d[A_GEB] = upr_nominal_gdp
    d[A_Basic_Account] = upr_nominal_gdp
    d[A_Credit_Unions] = upr_nominal_gdp
    d[A_Endowment_Policy_Not_Linked] = upr_nominal_gdp

    d
end

const UPRATE_MAPPINGS =  make_uprate_types()
#
# FIXME type unstable? use the trick in the
#
UPRATING_DATA = nothing

"""
Load Quarterly UPRATING data into a dataframe, and recast everything relative to the target date (Y,Q).
See docs/notes.md on the data. Dataframe is a private global.
"""
function load_prices( settings :: Settings, reload :: Bool = false )

    global UPRATING_DATA
    if ((UPRATING_DATA !== nothing) && ( ! reload ))
        return
    end

    upr = CSV.File("$(PRICES_DIR)/$(settings.prices_file)"; delim = '\t', comment = "#") |> DataFrame
    nrows = size(upr)[1]
    ncols = size(upr)[2]
    println( "read $nrows rows and $ncols cols ")
    lcnames = Symbol.(basiccensor.(string.(names(upr))))
    rename!(upr, lcnames)

    upr[!,:year] = zeros(Int64, nrows)
    upr[!,:q] = zeros(Int8, nrows) #zeros(Union{Int64,Missing},np)

    # add year, quarter cols parsed from the 'YYYY QQ' field
    dp = r"([0-9]{4}) Q([1-4])"
    for i in 1:nrows
        rc = match(dp, upr[i, :date])
        if (rc !== nothing)
            upr[i, :year] = parse(Int64, rc[1])
            upr[i, :q] = parse(Int8, rc[2])
        end
    end

    # Make all relative to the target y,q
    pnew = findfirst((upr.year.== settings.to_y ) .& (upr.q.== settings.to_q ))
    for col in 1:ncols
        baser = upr[pnew,col]
        println( "on col $col $(lcnames[col]); baser=$baser")
        if ! (lcnames[col] in [:q, :year, :date ]) # got to be a better way
            upr[!,col] .= baser./upr[!,col]
        end
    end
    UPRATING_DATA = upr
end

function uprate( item :: Number, from_y::Integer, from_q::Integer, itype::Uprate_Item_Type)::Number
    # FIXME this is likely much too slow..
    global UPRATING_DATA
    if itype == upr_no_uprate
        return item
    end
    global Uprate_Map
    # global FROM_Y, FROM_Q
    colsym = Uprate_Map[itype]
    p = UPRATING_DATA[((UPRATING_DATA.year.==from_y).&(UPRATING_DATA.q.==from_q)), colsym][1]
    return item * p
end

end

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
@testset "montecarlo" begin
end


text:
# Base.@kwdef struct CSVSchema{T₁ <: AbstractVendor,  T₂ <: AbstractProduct, T₃ <: AbstractData} <: AbstractSchema
#     vendor::T₁ = GenericVendor()
#     product::T₂ = GenericProduct()
#     data::T₃ =
# end

struct CSVSchema{T₁ <: AbstractVendor,  T₂ <: AbstractProduct, T₃ <: AbstractData} <: AbstractSchema
    vendor::T₁
    product::T₂
    data::T₃
end


text:
using Test
using YaoHIR
using YaoLocations
using CompilerPluginTools
using ZXCalculus
using ZXCalculus: Phase, BlockIR
using YaoHIR: Chain, Gate, Ctrl, shift, Rz

ir = @make_ircode begin
    Expr(:call, :+, 1, 1)::Int
    Expr(:call, :+, 3, 3)::Int
end

circ = Chain(
    Gate(shift(1.0), Locations(1)),
    Gate(Rz(Core.SSAValue(2)), Locations(1)),
)

bir = BlockIR(ir, 4, circ)

zxd = convert_to_zxd(bir)
c = clifford_simplification(zxd)
qc_tl = convert_to_chain(phase_teleportation(zxd))
@test length(qc_tl) == 1