text:
module SkyTEM1DInversion
import ..AbstractOperator.get_misfit, ..main, ..gradientinv
import ..AbstractOperator.Sounding
import ..AbstractOperator.makeoperator
import ..AbstractOperator.getresidual

using ..AbstractOperator, ..AEM_VMD_HMD, Statistics, Distributed, Printf, Dates, StatsBase,
      PyPlot, LinearAlgebra, ..CommonToAll, Random, DelimitedFiles, LinearMaps, SparseArrays, ..GP

import ..Model, ..Options, ..OptionsStat, ..OptionsNonstat

export dBzdt, plotmodelfield!, addnoise_skytem, plotmodelfield!, plotmodelfield_skytem!

const μ₀ = 4*pi*1e-7

mutable struct dBzdt<:Operator1D
    dlow       :: Array{Float64, 1}
    dhigh      :: Array{Float64, 1}
    useML      :: Bool
    σlow       :: Array{Float64, 1}
    σhigh      :: Array{Float64, 1}
    Flow       :: AEM_VMD_HMD.HField
    Fhigh      :: AEM_VMD_HMD.HField
    z          :: Array{Float64, 1}
    nfixed     :: Int
    ρ          :: Array{Float64, 1}
    selectlow  :: Array{Bool, 1}
    selecthigh :: Array{Bool, 1}
    ndatalow   :: Int
    ndatahigh  :: Int
    J          :: AbstractArray
    W          :: SparseMatrixCSC
    res        :: Vector
end

struct Cinv
    # assumes diagonal data cov
    oneoverσ²
end

function (a::Cinv)(x)
    x.*a.oneoverσ²
end

function dBzdt(       Flow       :: AEM_VMD_HMD.HField,
                      Fhigh      :: AEM_VMD_HMD.HField,
                      dlow       :: Array{Float64, 1},
                      dhigh      :: Array{Float64, 1},
                      σlow       :: Array{Float64, 1},
                      σhigh      :: Array{Float64, 1};
                      useML  = false,
                      z = [-1.],
                      ρ = [-1.],
                      nfixed = 1
                    )
    @assert length(Flow.thickness) >= length(z)
    @assert length(Fhigh.thickness) >= length(z)
    @assert size(σlow)  == size(dlow)
    @assert size(σhigh) == size(dhigh)
    ndatalow   = sum(.!isnan.(dlow))
    ndatahigh  = sum(.!isnan.(dhigh))
    selectlow  = .!isnan.(dlow)
    selecthigh = .!isnan.(dhigh)
    # for Gauss-Newton
    calcjacobian = Flow.calcjacobian & Fhigh.calcjacobian
    res = [dlow[selectlow];dhigh[selecthigh]]
    if calcjacobian
        J = [Flow.dBzdt_J'; Fhigh.dBzdt_J'];
        J = J[[selectlow;selecthigh],nfixed+1:length(ρ)]
        Wdiag = [1 ./σlow[selectlow]; 1 ./σhigh[selecthigh]]
        res = [dlow[selectlow];dhigh[selecthigh]]
    else
        res, J, Wdiag = zeros(0), zeros(0), zeros(0)
    end
    W = sparse(diagm(Wdiag))
    dBzdt(dlow, dhigh, useML, σlow, σhigh,
    Flow, Fhigh, z, nfixed, copy(ρ), selectlow, selecthigh, ndatalow, ndatahigh, J, W, res)
end

function getresidual(aem::dBzdt, log10σ::Vector{Float64}; computeJ=false)
    # aem.ρ[:] = 10 .^(-log10σ)
    aem.Flow.calcjacobian = computeJ
    aem.Fhigh.calcjacobian = computeJ
    getfield!(-log10σ, aem)
    fl, fh = aem.Flow.dBzdt[aem.selectlow], aem.Fhigh.dBzdt[aem.selecthigh]
    dl, dh = aem.dlow[aem.selectlow], aem.dhigh[aem.selecthigh]
    aem.res[:] = [fl-dl; fh-dh]
    if computeJ
        selectlow, selecthigh, nfixed = aem.selectlow, aem.selecthigh, aem.nfixed
        Flow, Fhigh = aem. Flow, aem.Fhigh
        copy!(aem.J, [Flow.dBzdt_J'[selectlow,nfixed+1:nfixed+length(log10σ)];
                      Fhigh.dBzdt_J'[selecthigh,nfixed+1:nfixed+length(log10σ)]])
    end
    nothing
end

mutable struct SkyTEMsoundingData <: Sounding
    sounding_string :: String
    X :: Float64
    Y :: Float64
    Z :: Float64
    fid :: Float64
    linenum :: Int
    rRx :: Float64
    zRxLM :: Float64
    zTxLM :: Float64
    zRxHM :: Float64
    zTxHM :: Float64
    rTx :: Float64
    lowpassfcs :: Array{Float64, 1}
    LM_times :: Array{Float64, 1}
    LM_ramp :: Array{Float64, 2}
    HM_times :: Array{Float64, 1}
    HM_ramp :: Array{Float64, 2}
    LM_noise :: Array{Float64, 1}
    HM_noise :: Array{Float64, 1}
    LM_data :: Array{Float64, 1}
    HM_data :: Array{Float64, 1}
end

function SkyTEMsoundingData(;rRx=-12., zRxLM=12., zTxLM=12.,
                            zRxHM=12., zTxHM=12., rTx=-12.,lowpassfcs=[-1, -2.],
                            LM_times=[1., 2.], LM_ramp=[1 2; 3 4],
                            HM_times=[1., 2.], HM_ramp=[1 2; 3 4],
                            LM_noise=[1.], HM_noise=[1.], LM_data=[1.], HM_data=[1.],
                            sounding_string="sounding", X=nothing, Y=nothing, Z=nothing,
                            linenum=nothing, fid=nothing)
    @assert rRx > 0 && rTx > 0
    @assert zRxLM <0 && zTxLM <0
    @assert zRxHM <0 && zTxHM <0
    @assert all(lowpassfcs .> 0)
    @assert all(diff(LM_times) .>0 )
    @assert all(diff(HM_times) .>0 )
    @assert all(diff(LM_ramp[:,1]) .>0 )
    @assert all(diff(HM_ramp[:,1]) .>0 )
    @assert all((LM_noise .>0) .| isnan.(LM_noise))
    @assert all((HM_noise .>0) .| isnan.(HM_noise))
    @assert length(LM_data) == length(LM_noise)
    @assert length(HM_data) == length(HM_noise)
    SkyTEMsoundingData(sounding_string, X, Y, Z, fid, linenum, rRx, zRxLM, zTxLM, zRxHM, zTxHM, rTx,
    lowpassfcs, LM_times, LM_ramp, HM_times, HM_ramp, LM_noise, HM_noise,
    LM_data, HM_data)
    # @show (sounding_string, rRx, zRxLM, zTxLM, zRxHM, zTxHM, rTx,
    # lowpassfcs, LM_times, LM_ramp, HM_times, HM_ramp, LM_noise, HM_noise,
    # LM_data, HM_data)
end

function getlocationinhdr(str, rows)
    for (irow, row) in enumerate(rows)
        lasttabpos = findlast('\t', row)
        (str == row[lasttabpos+1:end]) && return irow
    end
    nothing
end

function getcolNo(str)
	index = findfirst('\t',str)
	value = str[1:index-1]
    return parse(Int,value)
end

function getcolNo_except(str)
    index1 = findall("\t",str)[1][1]
    index2 = findall("\t",str)[2][1]
    value1 = str[1:index1-1]
    value2 = str[index1+1:index2-1]
    return [parse(Int,value1), parse(Int,value2)]
end

function read_survey_files(dfnfile::String;
    fname_specs_halt="",
    frame_height = "",
    frame_dz = "",
    frame_dx = "",
    frame_dy = "",
    LM_Z = "",
    HM_Z = "",
    LM_σ = "",
    HM_σ = "",
    X = "",
    Y = "",
    Z = "",
    fid = "",
    linenum = "",
    LM_drop = nothing,
    HM_drop = nothing,
    relerror = false,
    units=1e-12,
    figsize = (9,7),
    makesounding = false,
    dotillsounding = nothing,
    startfrom = 1,
    skipevery = 1,
    multnoise = 0.03,
    )
    #throw an error if LM_drop or HM_drop are not provided
    (isnothing(LM_drop) || isnothing(HM_drop)) &&
        throw(ArgumentError("user must specify drop gates for LM and HM"))
    prefix = getgdfprefix(dfnfile)
    hdrfile = prefix*".hdr"
    if !isfile(hdrfile)
        dfn2hdr(dfnfile)
    else
        @warn "using existing "*hdrfile
    end
    dfn = readlines(hdrfile)
    frame_height, frame_dz, frame_dx,
    frame_dy, LM_Z, HM_Z, HM_σ,LM_σ,
    X, Y, Z, fid, linenum = map(x->getlocationinhdr(x, dfn), [frame_height, frame_dz, frame_dx,
    frame_dy, LM_Z, HM_Z,HM_σ, LM_σ, X , Y, Z, fid, linenum])

    # use function and map to get column number
    frame_height, frame_dz, frame_dx,
    frame_dy,X, Y, Z, fid, linenum = map(x -> getcolNo(dfn[x]), [frame_height, frame_dz, frame_dx,
    frame_dy, X, Y, Z, fid, linenum])
    LM_Z, HM_Z, LM_σ, HM_σ= map(x -> getcolNo_except(dfn[x]),[LM_Z,HM_Z,LM_σ, HM_σ])
    LM_Z[1] = LM_Z[1] + LM_drop
    HM_Z[1] = HM_Z[1] + HM_drop
    fname_dat = prefix*".dat" # the name of DAT file associated with DFN
    s_array = read_survey_files(fname_dat = fname_dat,
                        fname_specs_halt = fname_specs_halt,
                        LM_Z             = LM_Z,
                        HM_Z             = HM_Z,
                        frame_height     = frame_height,
                        frame_dz         = frame_dz,
                        frame_dy         = frame_dy,
                        frame_dx         = frame_dx,
                        X                = X,
                        Y                = Y,
                        Z                = Z,
                        fid              = fid,
                        units            = units,
                        relerror         = relerror,
                        LM_σ             = LM_σ,
                        HM_σ             = HM_σ,
                        figsize          = figsize,
                        linenum          = linenum,
                        startfrom        = startfrom,
                        skipevery        = skipevery,
                        dotillsounding   = dotillsounding,
                        multnoise        = multnoise,
                        makesounding     = makesounding)
    s_array # return sounding array
end

function read_survey_files(;
    fname_dat="",
    fname_specs_halt="",
    frame_height = -2,
    frame_dz = -2,
    frame_dx = -2,
    frame_dy = -2,
    LM_Z = [-2, -2],
    HM_Z = [-2, -2],
    LM_σ = [-2, 2],
    HM_σ = [-2, 2],
    relerror = false,
    units=1e-12,
    figsize = (9,7),
    makesounding = false,
    dotillsounding = nothing,
    startfrom = 1,
    skipevery = 1,
    multnoise = 0.03,
    X = -1,
    Y = -1,
    Z = -1,
    fid = -1,
    linenum = -1,
    nanchar = "*")
    @assert frame_height > 0
    @assert frame_dz > 0
    @assert frame_dx > 0
    @assert frame_dy > 0
    @assert all(LM_Z .> 0)
    @assert all(HM_Z .> 0)
    if relerror
        @assert all(LM_σ .> 0)
        @assert all(HM_σ .> 0)
    end
    @assert X > 0
    @assert Y > 0
    @assert Z > 0
    @assert linenum > 0
    @assert fid > 0
    @info "reading $fname_dat"
    if dotillsounding!= nothing
        soundings = readdlm(fname_dat)[startfrom:skipevery:dotillsounding,:]
    else
        soundings = readdlm(fname_dat)[startfrom:skipevery:end,:]
    end
    soundings[soundings .== nanchar] .= NaN
    easting = soundings[:,X]
    northing = soundings[:,Y]
    topo = soundings[:,Z]
    fiducial = soundings[:,fid]
    whichline = soundings[:,linenum]
    d_LM = soundings[:,LM_Z[1]:LM_Z[2]]
    d_HM = soundings[:,HM_Z[1]:HM_Z[2]]
    if relerror
        σ_LM = soundings[:,LM_σ[1]:LM_σ[2]]
        σ_HM = soundings[:,HM_σ[1]:HM_σ[2]]
    end
    zTx = soundings[:,frame_height]
    zRx = -(zTx + soundings[:,frame_dz])
    zTx = -zTx
    rRx = sqrt.(soundings[:,frame_dx].^2 + soundings[:,frame_dy].^2)

    @info "reading $fname_specs_halt"
    include(fname_specs_halt)
    @assert size(d_LM, 2) == length(LM_times)
    @assert size(d_HM, 2) == length(HM_times)
    if !relerror
        @assert size(d_LM, 2) == length(LM_noise)
        @assert size(d_HM, 2) == length(HM_noise)
        LM_noise[:] .*= units
        HM_noise[:] .*= units
    else
        σ_LM[:] .*= units
        σ_HM[:] .*= units
    end
    d_LM[:]     .*= units
    d_HM[:]     .*= units
    f = figure(figsize=figsize)
    ax = Array{Any, 1}(undef, 4)
    ax[1] = subplot(2,2,1)
    nsoundings = size(soundings, 1)
    plot_dLM = permutedims(d_LM)
    plot_dLM[plot_dLM .<0] .= NaN
    pcolormesh(1:nsoundings, LM_times, log10.(plot_dLM), shading="nearest")
    xlabel("sounding #")
    cbLM = colorbar()
    cbLM.set_label("log₁₀ d_LM")
    ylabel("LM time s")
    axLM = ax[1].twiny()
    axLM.semilogy(LM_noise, LM_times)
    axLM.set_xlabel("high alt noise")
    ax[2] = subplot(2,2,3,sharex=ax[1], sharey=ax[1])
    plot_dHM = permutedims(d_HM)
    plot_dHM[plot_dHM .<0] .= NaN
    pcolormesh(1:nsoundings, HM_times, log10.(plot_dHM), shading="nearest")
    xlabel("sounding #")
    cbHM = colorbar()
    cbHM.set_label("log₁₀ d_HM")
    ylabel("HM time s")
    ax[2].invert_yaxis()
    axHM = ax[2].twiny()
    axHM.semilogy(HM_noise, HM_times)
    axHM.set_xlabel("high alt noise")
    ax[3] = subplot(2,2,2, sharex=ax[1])
    plot(1:nsoundings, zRx, label="Rx")
    plot(1:nsoundings, zTx, label="Tx")
    legend()
    xlabel("sounding #")
    ylabel("height m")
    ax[3].invert_yaxis()
    ax[4] = subplot(2,2,4, sharex=ax[1])
    ax[4].plot(1:nsoundings, rRx)
    xlabel("sounding #")
    ylabel("rRx m")
    plt.tight_layout()
    if makesounding
        s_array = Array{SkyTEMsoundingData, 1}(undef, nsoundings)
        for is in 1:nsoundings
            l, f = Int(whichline[is]), fiducial[is]
            @info "read $is out of $nsoundings"
            dlow, dhigh = vec(d_LM[is,:]), vec(d_HM[is,:])
            if !relerror
                σLM = sqrt.((multnoise*dlow).^2 + LM_noise.^2)
                σHM = sqrt.((multnoise*dhigh).^2 + HM_noise.^2)
            else
                σLM, σHM = vec(σ_LM[is,:]), vec(σ_HM[is,:])
            end
            s_array[is] = SkyTEMsoundingData(rRx=rRx[is], zRxLM=zRx[is], zTxLM=zTx[is],
                zRxHM=zRx[is], zTxHM=zTx[is], rTx=rTx, lowpassfcs=lowpassfcs,
                LM_times=LM_times, LM_ramp=LM_ramp,
                HM_times=HM_times, HM_ramp=HM_ramp,
                LM_noise=σLM, HM_noise=σHM, LM_data=dlow, HM_data=dhigh,
                sounding_string="sounding_$(l)_$f",
                X=easting[is], Y=northing[is], Z=topo[is], fid=f,
                linenum=l)
        end
        return s_array
    end
end

function getfield!(m::Model, aem::dBzdt)
    getfield!(m.fstar, aem)
    nothing
end

function getfield!(m::Array{Float64}, aem::dBzdt)
    copyto!(aem.ρ, aem.nfixed+1:length(aem.ρ), 10 .^m, 1:length(m))
    aem.ndatalow>0  && AEM_VMD_HMD.getfieldTD!(aem.Flow,  aem.z, aem.ρ)
    aem.ndatahigh>0 && AEM_VMD_HMD.getfieldTD!(aem.Fhigh, aem.z, aem.ρ)
    nothing
end

function get_misfit(m::Model, opt::Options, aem::dBzdt)
    chi2by2 = 0.0
    if !opt.debug
        getfield!(m, aem)
        aem.ndatalow>0 && (chi2by2 += getchi2by2(aem.Flow.dBzdt, aem.dlow,
                    aem.σlow, aem.selectlow, aem.useML, aem.ndatalow))
        aem.ndatahigh>0 && (chi2by2 += getchi2by2(aem.Fhigh.dBzdt, aem.dhigh,
                    aem.σhigh, aem.selecthigh, aem.useML, aem.ndatahigh))
    end
    return chi2by2
end

function getchi2by2(dBzdt, d, σ, select, useML, ndata)
    r, s, idx = dBzdt, σ, select
    r .= (r - d)./s
    if useML
        chi2by2 = 0.5*ndata*log(norm(r[idx])^2)
    else
        chi2by2 = 0.5*norm(r[idx])^2
    end
end

function computeMLfactor(dBzdt, d, σ, select, ndata)
    if ndata > 0
        r, s, idx = dBzdt, σ, select
        r .= (r - d)./s
        r[idx]'*r[idx]/ndata
    else
        NaN
    end
end

function computeMLfactor(aem)
    mlfact_low = computeMLfactor(aem.Flow.dBzdt, aem.dlow, aem.σlow, aem.selectlow, aem.ndatalow)
    mlfact_high = computeMLfactor(aem.Fhigh.dBzdt, aem.dhigh, aem.σhigh, aem.selecthigh, aem.ndatahigh)
    sqrt(mlfact_low), sqrt(mlfact_high)
end

function plotmodelfield_skytem!(Flow::AEM_VMD_HMD.HField, Fhigh::AEM_VMD_HMD.HField,
                         z::Array{Float64, 1}, ρ::Array{Float64, 1}
                        ;figsize=(10,5))
    plotwaveformgates(timesLM=Flow.times,  rampLM=Flow.ramp,
                      timesHM=Fhigh.times, rampHM=Fhigh.ramp)
    f, ax = plt.subplots(1, 2, figsize=figsize)
    ax[1].step(log10.(ρ[2:end]), z[2:end])
    ax[1].set_xlabel("log₁₀ρ")
    ax[2].set_xlabel("time s")
    ax[1].set_ylabel("depth m")
    ax[2].set_ylabel("V/Am⁴")
    AEM_VMD_HMD.getfieldTD!(Flow, z, ρ)
    AEM_VMD_HMD.getfieldTD!(Fhigh, z, ρ)
    ax[2].loglog(Flow.times,μ₀*Flow.dBzdt, label="low moment")
    ax[2].loglog(Fhigh.times,μ₀*Fhigh.dBzdt, label="high moment")
    ax[1].grid()
    ax[1].invert_yaxis()
    ax[2].grid()
    nicenup(f)
end

function addnoise_skytem(Flow::AEM_VMD_HMD.HField, Fhigh::AEM_VMD_HMD.HField,
                  z::Array{Float64, 1}, ρ::Array{Float64, 1};
                  halt_LM = nothing,
                  halt_HM = nothing,
                  noisefrac  = 0.03,
                  noisefloorlow = μ₀*1e-14,
                  noisefloorhigh = μ₀*1e-14,
                  dz = -1.,
                  extendfrac = -1.,
                  nfixed = -1,
                  figsize=(10,5),
                  rseed=42
                  )
    if halt_LM != nothing
        @assert length(halt_LM) == length(Flow.times)
    else
        halt_LM = zeros(length(Flow.times))
    end
    if halt_HM != nothing
        @assert length(halt_HM) == length(Fhigh.times)
    else
        halt_HM = zeros(length(Fhigh.times))
    end
    @assert all((nfixed, dz, extendfrac) .> 0)
    Random.seed!(rseed)
    AEM_VMD_HMD.getfieldTD!(Flow, z, ρ)
    AEM_VMD_HMD.getfieldTD!(Fhigh, z, ρ)
    dlow  = Flow.dBzdt + sqrt.((noisefrac*abs.(Flow.dBzdt)).^2 + (halt_LM/μ₀).^2).*randn(size(Flow.dBzdt))
    dhigh = Fhigh.dBzdt + sqrt.((noisefrac*abs.(Fhigh.dBzdt)).^2 + (halt_HM/μ₀).^2).*randn(size(Fhigh.dBzdt))
    dlow[abs.(dlow).<noisefloorlow] .= NaN
    dhigh[abs.(dhigh).<noisefloorhigh] .= NaN
    σlow = sqrt.((noisefrac*abs.(dlow)).^2 + (halt_LM/μ₀).^2)
    σhigh = sqrt.((noisefrac*abs.(dhigh)).^2 + (halt_HM/μ₀).^2)
    plotmodelfield!(Flow, Fhigh, z, ρ, dlow, dhigh, σlow, σhigh,
                                figsize=figsize, nfixed=nfixed,
                                dz=dz, extendfrac=extendfrac)
    # returned data is dBzdt not H if there is a μ multiplied
    return μ₀.*(dlow, dhigh, σlow, σhigh)
end

function plotmodelfield!(Flow::AEM_VMD_HMD.HField, Fhigh::AEM_VMD_HMD.HField,
                        z, ρ, dlow, dhigh, σlow, σhigh;
                        figsize=(10,5), nfixed=-1, dz=-1., extendfrac=-1.)
    # expects data and noise in units of H, i.e. B/μ
    @assert all((nfixed, dz, extendfrac) .> 0)
    f, ax = plt.subplots(1, 2, figsize=figsize)
    ax[1].step(log10.(ρ[2:end]), z[2:end])
    if dz > 0
        axn = ax[1].twinx()
        ax[1].get_shared_y_axes().join(ax[1],axn)
        axn.step(log10.(ρ[2:end]), z[2:end])
        yt = ax[1].get_yticks()[ax[1].get_yticks().>=0]
        axn.set_yticks(yt)
        axn.set_ylim(ax[1].get_ylim()[end:-1:1])
        axn.set_yticklabels(string.(Int.(round.(getn.(yt .- z[nfixed+1], dz, extendfrac)))))
    end
    ndatalow   = sum(.!isnan.(dlow))
    ndatahigh  = sum(.!isnan.(dhigh))
    if ndatalow>0
        AEM_VMD_HMD.getfieldTD!(Flow, z, ρ)
        ax[2].loglog(Flow.times,μ₀*Flow.dBzdt, label="low moment")
        ax[2].errorbar(Flow.times, μ₀*vec(dlow), yerr = μ₀*2abs.(vec(σlow)),
                            linestyle="none", marker=".", elinewidth=0, capsize=3)
    end
    if ndatahigh>0
        AEM_VMD_HMD.getfieldTD!(Fhigh, z, ρ)
        ax[2].loglog(Fhigh.times,μ₀*Fhigh.dBzdt, label="high moment")
        ax[2].errorbar(Fhigh.times, μ₀*vec(dhigh), yerr = μ₀*2abs.(vec(σhigh)),
                        linestyle="none", marker=".", elinewidth=0, capsize=3)
    end
    ax[1].grid()
    ax[2].grid()
        ax[1].step(log10.(ρ[2:end]), z[2:end])
    ax[1].set_xlabel("log₁₀ρ")
    ax[2].set_xlabel("time s")
    ax[1].set_ylabel("depth m")
    ax[2].set_ylabel("V/Am⁴")
    axn.set_ylabel("index no.")
    nicenup(f)
end

function plotmodelfield!(aem::dBzdt, Ρ::Vector{T};
                        figsize=(8,5), dz=-1., onesigma=true, onlygetMLsampled=false,
                        extendfrac=-1., fsize=10, alpha=0.1) where T<:AbstractArray
    @assert all((dz, extendfrac) .> 0)
    nfixed, z = aem.nfixed, aem.z
    if !onlygetMLsampled
        sigma = onesigma ? 1.0 : 2.0
        f = figure(figsize=figsize)
        ax = Vector{PyPlot.PyObject}(undef, 3)
        ax[1] = subplot(121)
        ρmin, ρmax = extrema(vcat(Ρ...))
        delρ = ρmax - ρmin
        ax[1].set_xlim(ρmin-0.1delρ,ρmax+0.1delρ)
        ax[1].plot([ρmin-0.1delρ,ρmax+0.1delρ], z[nfixed+1]*[1., 1], color="b")
        ax[2] = subplot(122)
    end
    Flow = aem.Flow
    dlow, σlow = aem.dlow, aem.σlow
    Fhigh = aem.Fhigh
    dhigh, σhigh = aem.dhigh, aem.σhigh
    if !onlygetMLsampled
        aem.ndatalow>0 && ax[2].errorbar(Flow.times, μ₀*dlow, yerr = μ₀*sigma*abs.(σlow),
                            linestyle="none", marker=".", elinewidth=0, capsize=3, label="low moment")
        aem.ndatahigh>0 && ax[2].errorbar(Fhigh.times, μ₀*dhigh, yerr = μ₀*sigma*abs.(σhigh),
                            linestyle="none", marker=".", elinewidth=0, capsize=3, label="high moment")
    else
        errorfact_low, errorfact_high = map(x->zeros(length(Ρ)), 1:2)
    end
    for (i, ρ) in enumerate(Ρ)
        getfield!(ρ,  aem)
        if !onlygetMLsampled
            if aem.ndatalow>0
                Flow.dBzdt[.!aem.selectlow] .= NaN
                ax[2].loglog(Flow.times,μ₀*Flow.dBzdt, "k", alpha=alpha, markersize=2)
            end
            if aem.ndatahigh>0
                Fhigh.dBzdt[.!aem.selecthigh] .= NaN
                ax[2].loglog(Fhigh.times,μ₀*Fhigh.dBzdt, "k", alpha=alpha, markersize=2)
            end
            ax[1].step(log10.(aem.ρ[2:end]), aem.z[2:end], "-k", alpha=alpha)
        else
            errorfact_low[i], errorfact_high[i] = computeMLfactor(aem)
        end
    end
    if !onlygetMLsampled
        ax[1].grid()
        ax[1].set_ylabel("Depth m")
        ax[1].plot(xlim(), z[nfixed+1]*[1, 1], "--k")
        if dz > 0
            axn = ax[1].twinx()
            ax[1].get_shared_y_axes().join(ax[1],axn)
            yt = ax[1].get_yticks()[ax[1].get_yticks().>=z[nfixed+1]]
            axn.set_yticks(yt)
            axn.set_ylim(ax[1].get_ylim()[end:-1:1])
            axn.set_yticklabels(string.(Int.(round.(getn.(yt .- z[nfixed+1], dz, extendfrac)))))
        end
        axn.set_ylabel("Depth index", rotation=-90, labelpad=10)
        ax[1].set_xlabel("Log₁₀ρ")
        ax[1].set_title("Model")
        ax[2].set_ylabel(L"dBzdt \; V/(A.m^4)")
        ax[2].set_xlabel("Time (s)")
        ax[2].set_title("Transient response")
        ax[2].legend()
        ax[2].grid()
        ax[1].invert_xaxis()
        nicenup(f, fsize=fsize)
        nothing, nothing, nothing, nothing
    else
        mean(errorfact_low), std(errorfact_low),
        mean(errorfact_high), std(errorfact_high)
    end
end

function makeoperator(sounding::SkyTEMsoundingData;
                       zfixed   = [-1e5],
                       ρfixed   = [1e12],
                       zstart = 0.0,
                       extendfrac = 1.06,
                       dz = 2.,
                       ρbg = 10,
                       nlayers = 40,
                       ntimesperdecade = 10,
                       nfreqsperdecade = 5,
                       showgeomplot = false,
                       useML = false,
                       calcjacobian = false,
                       modelprimary = false,
                       plotfield = false)
    @assert extendfrac > 1.0
    @assert dz > 0.0
    @assert ρbg > 0.0
    @assert nlayers > 1
    nmax = nlayers+1

    zall, znall, zboundaries = setupz(zstart, extendfrac, dz=dz, n=nlayers, showplot=showgeomplot)
    z, ρ, nfixed = makezρ(zboundaries; zfixed=zfixed, ρfixed=ρfixed)
    ρ[z.>=zstart] .= ρbg
    ## LM operator
    Flm = AEM_VMD_HMD.HFieldDHT(
                          ntimesperdecade = ntimesperdecade,
                          nfreqsperdecade = nfreqsperdecade,
                          lowpassfcs = sounding.lowpassfcs,
                          times  = sounding.LM_times,
                          ramp   = sounding.LM_ramp,
                          nmax   = nmax,
                          zTx    = sounding.zTxLM,
                          rRx    = sounding.rRx,
                          rTx    = sounding.rTx,
                          zRx    = sounding.zRxLM,
                          modelprimary = modelprimary,
                          calcjacobian = calcjacobian
                          )
    ## HM operator
    Fhm = AEM_VMD_HMD.HFieldDHT(
                          ntimesperdecade = ntimesperdecade,
                          nfreqsperdecade = nfreqsperdecade,
                          lowpassfcs = sounding.lowpassfcs,
                          times  = sounding.HM_times,
                          ramp   = sounding.HM_ramp,
                          nmax   = nmax,
                          zTx    = sounding.zTxHM,
                          rRx    = sounding.rRx,
                          rTx    = sounding.rTx,
                          zRx    = sounding.zRxHM,
                          modelprimary = modelprimary,
                          calcjacobian = calcjacobian
                          )
    ## create operator
    dlow, dhigh, σlow, σhigh = (sounding.LM_data, sounding.HM_data, sounding.LM_noise, sounding.HM_noise)./μ₀
    aem = dBzdt(Flm, Fhm, vec(dlow), vec(dhigh), useML=useML,
                vec(σlow), vec(σhigh), z=z, ρ=ρ, nfixed=nfixed)
    if plotfield
        plotwaveformgates(aem)
        plotmodelfield!(Flm, Fhm, z, ρ, dlow, dhigh, σlow, σhigh;
                          figsize=(12,4), nfixed=nfixed, dz=dz, extendfrac=extendfrac)
    end
    aem, znall
end

function plotwaveformgates(aem::dBzdt; figsize=(10,5))
    plotwaveformgates(timesLM=aem.Flow.times,  rampLM=aem.Flow.ramp,
                      timesHM=aem.Fhigh.times, rampHM=aem.Fhigh.ramp)
end

function plotwaveformgates(;timesLM=nothing, rampLM=nothing,
                            timesHM=nothing, rampHM=nothing,
                            figsize=(10,5))
    @assert timesLM != nothing
    @assert timesHM != nothing
    @assert rampLM != nothing
    @assert rampHM != nothing

    f = figure(figsize=figsize)
    s1 = subplot(121)
    plot(rampLM[:,1]*1e6, rampLM[:,2], "-or")
    stem(timesLM*1e6, ones(length(timesLM)))
    ylabel("Amplitude")
    xlabel("time μs")
    title("LM linear time")
    s2 = subplot(122, sharey=s1)
    plot(rampHM[:,1]*1e3, rampHM[:,2], "-or")
    stem(timesHM*1e3, ones(length(timesHM)))
    ylabel("Amplitude")
    xlabel("time ms")
    title("HM linear time")
    # s3 = subplot(222)
    # loglog(rampLM[:,1], rampLM[:,2], "-or")
    # stem(timesLM, ones(length(timesLM)))
    # title("LM log time")
    # s4 = subplot(224, sharex=s3, sharey=s3)
    # loglog(rampHM[:,1], rampHM[:,2], "-or")
    # stem(timesHM, ones(length(timesHM)))
    # xlabel("time s")
    # title("HM log time")
    nicenup(f, fsize=12)
end

function make_tdgp_opt(;
                    rseed = nothing,
                    znall = znall,
                    fileprefix = "sounding",
                    nmin = 2,
                    nmax = 40,
                    K = GP.OrstUhn(),
                    demean = true,
                    sdpos = 0.05,
                    sdprop = 0.05,
                    sddc = 0.008,
                    sampledc = false,
                    fbounds = [-0.5 2.5],
                    λ = [2],
                    δ = 0.1,
                    pnorm = 2,
                    save_freq = 25,
                    restart = false
                    )
    sdev_pos = [sdpos*abs(diff([extrema(znall)...])[1])]
    sdev_prop = sdprop*diff(fbounds, dims=2)[:]
    sdev_dc = sddc*diff(fbounds, dims=2)[:]
    xall = permutedims(collect(znall))
    xbounds = permutedims([extrema(znall)...])

    history_mode = "w"
	restart && (history_mode = "a")

    updatenonstat = false
    needλ²fromlog = false
    if rseed != nothing
        Random.seed!(rseed)
    end
    opt = OptionsStat(fdataname = fileprefix*"_",
                            nmin = nmin,
                            nmax = nmax,
                            xbounds = xbounds,
                            fbounds = fbounds,
                            xall = xall,
                            λ = λ,
                            δ = δ,
                            demean = demean,
                            sdev_prop = sdev_prop,
                            sdev_pos = sdev_pos,
                            sdev_dc = sdev_dc,
                            sampledc = sampledc,
                            pnorm = pnorm,
                            quasimultid = false,
                            K = K,
                            save_freq = save_freq,
                            needλ²fromlog = needλ²fromlog,
                            updatenonstat = updatenonstat,
                            dispstatstoscreen = false,
                            history_mode = history_mode
                            )
    opt
end

function makeoperatorandoptions(soundings::Array{SkyTEMsoundingData, 1};
                        rseed = nothing,
                        fileprefix = "sounding",
                        nmin = 2,
                        nmax = 40,
                        K = GP.OrstUhn(),
                        demean = true,
                        sdpos = 0.05,
                        sdprop = 0.05,
                        sddc = 0.008,
                        sampledc = false,
                        fbounds = [-0.5 2.5],
                        λ = [2],
                        δ = 0.1,
                        pnorm = 2,
                        save_freq = 25,
                        restart = false,
                        nplot = 2,
                        zfixed   = [-1e5],
                        ρfixed   = [1e12],
                        zstart = 0.0,
                        extendfrac = 1.06,
                        dz = 2.,
                        ρbg = 10,
                        nlayers = 40,
                        ntimesperdecade = 10,
                        nfreqsperdecade = 5,
                        showgeomplot = false,
                        useML = false,
                        modelprimary = false
                        )
    aem, opt = [], []
    if rseed != nothing
        Random.seed!(rseed)
    end
    for idx in randperm(length(soundings))[1:nplot]
            aem, znall = makeoperator(soundings[idx],
                        zfixed = zfixed,
                        ρfixed = ρfixed,
                        zstart = zstart,
                        extendfrac = extendfrac,
                        dz = dz,
                        ρbg = ρbg,
                        nlayers = nlayers,
                        ntimesperdecade = ntimesperdecade,
                        nfreqsperdecade = nfreqsperdecade,
                        useML = useML,
                        showgeomplot = showgeomplot,
                        modelprimary = modelprimary,
                        plotfield = true)

                    opt = make_tdgp_opt(znall = znall,
                        fileprefix = soundings[idx].sounding_string,
                        nmin = nmin,
                        nmax = nmax,
                        K = K,
                        demean = demean,
                        sdpos = sdpos,
                        sdprop = sdprop,
                        sddc = sddc,
                        sampledc = sampledc,
                        fbounds = fbounds,
                        save_freq = save_freq,
                        λ = λ,
                        δ = δ,
                        restart = restart
                        )
    end
    # returns generically all the required elements in operator
    # for plotting except sounding_string
    opt
end

function summarypost(soundings::Array{SkyTEMsoundingData, 1}, opt::Options;
            qp1=0.05,
            qp2=0.95,
            burninfrac=0.5,
            zstart = 0.0,
            extendfrac = -1,
            dz = -1,
            nlayers = -1,
            useML=false)

    @assert extendfrac > 1.0
    @assert dz > 0.0
    @assert nlayers > 1
    zall, = setupz(zstart, extendfrac, dz=dz, n=nlayers)

    linename = "_line_$(soundings[1].linenum)_summary.txt"
    fnames = ["rho_low", "rho_mid", "rho_hi", "rho_avg",
              "ddz_mean", "ddz_sdev", "phid_mean", "phid_sdev"].*linename
    if isfile(fnames[1])
        @warn fnames[1]*" exists, reading stored values"
        pl, pm, ph, ρmean,
        vdmean, vddev, χ²mean, χ²sd, = map(x->readdlm(x), fnames)
    else
        # this is a dummy operator for plotting
        aem, = makeoperator(soundings[1])
        # now get the posterior marginals
        opt.xall[:] .= zall
        pl, pm, ph, ρmean, vdmean, vddev = map(x->zeros(length(zall), length(soundings)), 1:6)
        χ²mean, χ²sd = zeros(length(soundings)), zeros(length(soundings))
        for idx = 1:length(soundings)
            @info "$idx out of $(length(soundings))\n"
            opt.fdataname = soundings[idx].sounding_string*"_"
            opt.xall[:] .= zall
            pl[:,idx], pm[:,idx], ph[:,idx], ρmean[:,idx],
            vdmean[:,idx], vddev[:,idx] = CommonToAll.plot_posterior(aem, opt, burninfrac=burninfrac,
                                                    qp1=qp1, qp2=qp2,
                                                    doplot=false)
            χ² = 2*CommonToAll.assembleTat1(opt, :U, temperaturenum=1, burninfrac=burninfrac)
            ndata = sum(.!isnan.(soundings[idx].LM_data)) +
                    sum(.!isnan.(soundings[idx].HM_data))
            χ²mean[idx] = mean(χ²)/ndata
            χ²sd[idx]   = std(χ²)/ndata
            useML && (χ²mean[idx] -= log(ndata))
        end
        # write in grid format
        for (fname, vals) in Dict(zip(fnames, [pl, pm, ph, ρmean, vdmean, vddev, χ²mean, χ²sd]))
            writedlm(fname, vals)
        end
        # write in x, y, z, rho format
        for (i, d) in enumerate([pl, pm, ph, ρmean])
            xyzrho = makearray(soundings, d, zall)
            writedlm(fnames[i][1:end-4]*"_xyzrho.txt", xyzrho)
        end
    end
    pl, pm, ph, ρmean, vdmean, vddev, χ²mean, χ²sd, zall
end

function writetabdelim(fname, opt::Options, soundings::Array{SkyTEMsoundingData, 1};
                        nbins=50, qp1=0.05, qp2=0.95, burninfrac=0.5,
                        zstart = 0.0,
                        extendfrac = -1,
                        dz = -1,
                        nlayers = -1,)
    @assert extendfrac > 1.0
    @assert dz > 0.0
    @assert nlayers > 1
    zall, = setupz(zstart, extendfrac, dz=dz, n=nlayers)
    io = open(fname, "w")
    for (idx, sounding) in enumerate(soundings)
        @info "Sounding number: $idx out of $(length(soundings))"
        opt.fdataname = sounding.sounding_string*"_"
        himage, edges, CI, meanimage, meandiffimage, sdslope = make1Dhist(opt, burninfrac=burninfrac, nbins = nbins, qp1=qp1, qp2=qp2, temperaturenum=1)
        χ² = 2*assembleTat1(opt, :U, temperaturenum=1, burninfrac=0.5)
        ndata = sum(.!isnan.(sounding.LM_data)) + sum(.!isnan.(sounding.HM_data))
        χ²mean, χ²sd = mean(χ²)/ndata, std(χ²)/ndata
        towrite = [:X, :Y, :Z, :fid, :linenum, :rRx, :zRxLM, :zTxLM, :zRxHM, :zTxHM, :rTx]
        for tw in towrite # writes sounding parts
            msg = @sprintf("%e\t", getfield(sounding, tw))
            write(io, msg)
        end
        write(io, "\t")
        # zall never was missing
        for el in [vec(zall); -CI[:,3]; -CI[:,2]; -vec(meanimage); -CI[:,1]; χ²mean; χ²sd; -vec(meandiffimage); -vec(sdslope)] # concatenate everything else
            msg = @sprintf("%e\t", el)
            write(io, msg)
        end
        write(io, "\n")
        flush(io) # slower but ensures write is complete
    end
    close(io)
end

function plotindividualsoundings(soundings::Array{SkyTEMsoundingData, 1}, opt::Options;
    burninfrac=0.5,
    nbins=100,
    figsize  = (12,6),
    zfixed   = [-1e5],
    ρfixed   = [1e12],
    zstart = 0.0,
    extendfrac = 1.06,
    dz = 2.,
    ρbg = 10,
    omittemp = false,
    nlayers = 40,
    ntimesperdecade = 10,
    nfreqsperdecade = 5,
    computeforwards = false,
    nforwards = 100,
    rseed = 11,
    modelprimary = false,
    idxcompute = [1],
    onlygetMLsampled = false)

    if onlygetMLsampled
        computeforwards = true
        errorfac_low, errorfacσ_low, errorfac_high, errorfacσ_high = map(x->zeros(length(soundings)), 1:4)
    end
    zall, = setupz(zstart, extendfrac, dz=dz, n=nlayers)
    opt.xall[:] .= zall
    for idx = 1:length(soundings)
        if in(idx, idxcompute)
            @info "Sounding number: $idx"
            opt.fdataname = soundings[idx].sounding_string*"_"
            aem, znall = makeoperator(soundings[idx],
                zfixed = zfixed,
                ρfixed = ρfixed,
                zstart = zstart,
                extendfrac = extendfrac,
                dz = dz,
                ρbg = ρbg,
                nlayers = nlayers,
                ntimesperdecade = ntimesperdecade,
                nfreqsperdecade = nfreqsperdecade,
                modelprimary = modelprimary)
            if !onlygetMLsampled
                getchi2forall(opt, alpha=0.8, omittemp=omittemp)
                CommonToAll.getstats(opt)
                plot_posterior(aem, opt, burninfrac=burninfrac, nbins=nbins, figsize=figsize)
                ax = gcf().axes
                ax[1].invert_xaxis()
            end
            if computeforwards
                M = assembleTat1(opt, :fstar, temperaturenum=1, burninfrac=burninfrac)
                Random.seed!(rseed)
                if onlygetMLsampled
                    errorfac_low[idx], errorfacσ_low[idx],
                    errorfac_high[idx], errorfacσ_high[idx] = plotmodelfield!(aem, M[randperm(length(M))[1:nforwards]],
                        dz=dz, extendfrac=extendfrac, onesigma=false, alpha=0.2, onlygetMLsampled=true)
                else
                    plotmodelfield!(aem, M[randperm(length(M))[1:nforwards]],
                    dz=dz, extendfrac=extendfrac, onesigma=false, alpha=0.2)
                end
            end
        end
    end
    if onlygetMLsampled
        errorfac_low, errorfacσ_low, errorfac_high, errorfacσ_high
    end
end

function plotsummarygrids1(meangrid, phgrid, plgrid, pmgrid, gridx, gridz, topofine, R, Z, χ²mean, χ²sd, lname; qp1=0.05, qp2=0.95,
                        figsize=(10,10), fontsize=12, cmap="viridis", vmin=-2, vmax=0.5, Eislast=true, Nislast=true,
                        topowidth=2, idx=nothing, omitconvergence=false, useML=false, preferEright=false, preferNright=false,
                        saveplot=false, yl=nothing, botadjust=[0.125, 0.05, 0.75, 0.01], dpi=300, showplot=true)
    f = figure(figsize=figsize)
    dr = diff(gridx)[1]
    f.suptitle(lname*" Δx=$dr m, Fids: $(length(R))", fontsize=fontsize)
    nrows = omitconvergence ? 4 : 5
    icol = 1
    s = Array{Any, 1}(undef, nrows)
    if !omitconvergence
        s[icol] = subplot(nrows, 1, icol)
        if useML
            plot(R, exp.(χ²mean))
            semilogy(R, ones(length(R)), "--k")
            ylabel("variance factor")
            title("Max likelihood variance adjustment")
        else
            plot(R, χ²mean)
            plot(R, ones(length(R)), "--k")
            fill_between(R, vec(χ²mean-χ²sd), vec(χ²mean+χ²sd), alpha=0.5)
            ylabel(L"ϕ_d")
            title("Data misfit")
        end
        icol += 1
    end
    s[icol] = omitconvergence ? subplot(nrows, 1, icol) : subplot(nrows, 1, icol, sharex=s[icol-1])
    imshow(plgrid, cmap=cmap, aspect="auto", vmax=vmax, vmin = vmin,
                extent=[gridx[1], gridx[end], gridz[end], gridz[1]])
    plot(gridx, topofine, linewidth=topowidth, "-k")
    idx == nothing || plotprofile(s[icol], idx, Z, R)
    title("Percentile $(round(Int, 100*qp1)) conductivity")
    ylabel("Height m")
    icol += 1
    s[icol] = subplot(nrows, 1, icol, sharex=s[icol-1], sharey=s[icol-1])
    imshow(pmgrid, cmap=cmap, aspect="auto", vmax=vmax, vmin = vmin,
                extent=[gridx[1], gridx[end], gridz[end], gridz[1]])
    plot(gridx, topofine, linewidth=topowidth, "-k")
    title("Percentile 50 conductivity")
    idx == nothing || plotprofile(s[icol], idx, Z, R)
    ylabel("Height m")
    icol += 1
    s[icol] = subplot(nrows, 1, icol, sharex=s[icol-1], sharey=s[icol-1])
    imshow(meangrid, cmap=cmap, aspect="auto", vmax=vmax, vmin = vmin,
                extent=[gridx[1], gridx[end], gridz[end], gridz[1]])
    plot(gridx, topofine, linewidth=topowidth, "-k")
    title("Mean conductivity")
    ylabel("Height m")
    idx == nothing || plotprofile(s[icol], idx, Z, R)
    icol +=1
    s[icol] = subplot(nrows, 1, icol, sharex=s[icol-1], sharey=s[icol-1])
    imlast = imshow(phgrid, cmap=cmap, aspect="auto", vmax=vmax, vmin = vmin,
                extent=[gridx[1], gridx[end], gridz[end], gridz[1]])
    plot(gridx, topofine, linewidth=topowidth, "-k")
    xlabel("Line distance m")
    title("Percentile $(round(Int, 100*qp2)) conductivity")
    ylabel("Height m")
    idx == nothing || plotprofile(s[icol], idx, Z, R)
    xlim(extrema(gridx))
    map(x->x.tick_params(labelbottom=false), s[1:end-1])
    map(x->x.grid(), s)
    nicenup(f, fsize=fontsize)
    isa(yl, Nothing) || s[end].set_ylim(yl...)
    plotNEWSlabels(Eislast, Nislast, gridx, gridz, s)
    f.subplots_adjust(bottom=0.125)
    cbar_ax = f.add_axes(botadjust)
    cb = f.colorbar(imlast, cax=cbar_ax, orientation="horizontal",)
    cb.ax.set_xlabel("Log₁₀ S/m", fontsize=fontsize)
    cb.ax.tick_params(labelsize=fontsize)
    (preferNright && !Nislast) && s[end].invert_xaxis()
    (preferEright && !Eislast) && s[end].invert_xaxis()
    saveplot && savefig(lname*".png", dpi=dpi)
    showplot || close(f)
end

function plotsummarygrids2(σmeangrid, ∇zmeangrid, ∇zsdgrid, cigrid, gridx, gridz, topofine, lname;
        qp1=0.05, qp2=0.95, Eislast=true, Nislast=true,
        figsize=(10,10), fontsize=12, cmap="viridis", vmin=-2, vmax=0.5, topowidth=2)
    f = figure(figsize=figsize)
    f.suptitle(lname, fontsize=fontsize)
    s1 = subplot(411)
    imshow(σmeangrid, cmap=cmap, aspect="auto", vmax=vmax, vmin = vmin,
                extent=[gridx[1], gridx[end], gridz[end], gridz[1]])
    title("Mean conductivity")
    ylabel("Height m")
    colorbar()
    s2 = subplot(412, sharex=s1)
    imshow(abs.(cigrid), cmap=cmap, aspect="auto",
                extent=[gridx[1], gridx[end], gridz[end], gridz[1]])
    plot(gridx, topofine, linewidth=topowidth, "-k")
    title("CI: $(round(Int, 100*(qp2-qp1))) conductivity")
    ylabel("Height m")
    colorbar()
    s3 = subplot(413, sharex=s2, sharey=s2)
    imshow(∇zmeangrid, cmap=cmap, aspect="auto",
                extent=[gridx[1], gridx[end], gridz[end], gridz[1]])
    plot(gridx, topofine, linewidth=topowidth, "-k")
    title("Mean conductivity vertical derivative")
    ylabel("Height m")
    colorbar()
    s4 = subplot(414, sharex=s2, sharey=s2)
    imlast = imshow(∇zsdgrid, cmap=cmap, aspect="auto",
                extent=[gridx[1], gridx[end], gridz[end], gridz[1]])
    plot(gridx, topofine, linewidth=topowidth, "-k")
    xlabel("Line distance m")
    title("Std dev of conductivity vertical derivative")
    ylabel("Height m")
    colorbar()
    xlim(extrema(gridx))
    nicenup(f, fsize=fontsize)
end

function splitlinesummaryimages(soundings::Array{SkyTEMsoundingData, 1}, opt::Options;
                        qp1=0.05,
                        qp2=0.95,
                        burninfrac=0.5,
                        zstart = 0.0,
                        extendfrac = -1,
                        dz = -1,
                        dr = 10,
                        nlayers = -1,
                        fontsize = 10,
                        vmin = -2,
                        vmax = 0.5,
                        cmap="viridis",
                        figsize=(6,10),
                        topowidth=2,
                        idx = nothing,
                        showderivs = false,
                        omitconvergence = false,
                        useML = false,
                        preferEright = false,
                        preferNright = false,
                        saveplot = true,
                        yl = nothing,
                        showplot = true,
                        botadjust=[0.125, 0.05, 0.75, 0.01],
                        dpi = 300)

    linestartidx = splitsoundingsbyline(soundings)
    nlines = length(linestartidx)
    for i in 1:nlines
        a = linestartidx[i]
        b = i != nlines ?  linestartidx[i+1]-1 : length(soundings)
        summaryimages(soundings[a:b], opt, qp1=qp1, qp2=qp2, burninfrac=burninfrac, zstart=zstart,
            extendfrac=extendfrac, dz=dz, dr=dr, nlayers=nlayers, fontsize=fontsize, vmin=vmin,
            vmax = vmax, cmap=cmap, figsize=figsize, topowidth=topowidth, idx=idx, showderivs=showderivs,
            omitconvergence=omitconvergence, useML=useML, preferEright=preferEright, showplot=showplot,
            preferNright=preferNright, saveplot=saveplot, yl=yl, botadjust=botadjust, dpi=dpi)
    end
    nothing
end

function summaryimages(soundings::Array{SkyTEMsoundingData, 1}, opt::Options;
                        qp1=0.05,
                        qp2=0.95,
                        burninfrac=0.5,
                        zstart = 0.0,
                        extendfrac = -1,
                        dz = -1,
                        dr = 10,
                        nlayers = -1,
                        fontsize = 10,
                        vmin = -2,
                        vmax = 0.5,
                        cmap="viridis",
                        figsize=(6,10),
                        topowidth=2,
                        idx = nothing,
                        showderivs = false,
                        omitconvergence = false,
                        useML = false,
                        preferEright = false,
                        preferNright = false,
                        saveplot = false,
                        yl = nothing,
                        showplot = true,
                        botadjust=[0.125, 0.05, 0.75, 0.01],
                        dpi = 300)
    @assert !(preferNright && preferEright) # can't prefer both labels to the right
    pl, pm, ph, ρmean, vdmean, vddev, χ²mean, χ²sd, zall = summarypost(soundings, opt,
                                                                    zstart=zstart,
                                                                    qp1=qp1,
                                                                    qp2=qp2,
                                                                    dz=dz,
                                                                    extendfrac=extendfrac,
                                                                    nlayers=nlayers,
                                                                    burninfrac=burninfrac, useML=useML)

    phgrid, plgrid, pmgrid, σmeangrid, ∇zmeangrid,
    ∇zsdgrid, gridx, gridz, topofine, R, Z = makesummarygrid(soundings, pl, pm, ph, ρmean,
                                                            vdmean, vddev, zall, dz, dr=dr)

    lname = "Line $(soundings[1].linenum)"
    Eislast, Nislast = whichislast(soundings)
    plotsummarygrids1(σmeangrid, phgrid, plgrid, pmgrid, gridx, gridz, topofine, R, Z, χ²mean, χ²sd, lname, qp1=qp1, qp2=qp2,
                        figsize=figsize, fontsize=fontsize, cmap=cmap, vmin=vmin, vmax=vmax, Eislast=Eislast,
                        Nislast=Nislast, topowidth=topowidth, idx=idx, omitconvergence=omitconvergence, useML=useML,
                        preferEright=preferEright, preferNright=preferNright, saveplot=saveplot, showplot=showplot, dpi=dpi,
                        yl=yl, botadjust=botadjust)
    if showderivs
        cigrid = phgrid - plgrid
        plotsummarygrids2(σmeangrid, ∇zmeangrid, ∇zsdgrid, cigrid, gridx, gridz, topofine, lname, qp1=qp1, qp2=qp2,
                        figsize=figsize, fontsize=fontsize, cmap=cmap, vmin=vmin, vmax=vmax, topowidth=topowidth,
                        Eislast=Eislast, Nislast=Nislast)
    end
end

# for deterministic inversions, read in
function compress(soundings, zall; prefix="", rmfile=true, isfirstparalleliteration=false)
    fname = soundings[1].sounding_string*"_gradientinv.dat"
    !isfile(fname) && throw(AssertionError("file does not exist perhaps soundings already zipped?"))
    fout = prefix == "" ? "zipped.dat" : prefix*"_zipped.dat"
    iomode = "w"
    if isfile(fout)
        isfirstparalleliteration && throw(AssertionError("Zipped file "*fout*" exists, will not overwrite!"))
        iomode = "a"
    end
    io = open(fout, iomode)
    for (i, s) in enumerate(soundings)
        fname = s.sounding_string*"_gradientinv.dat"
        A = readdlm(fname)
        ϕd = A[end,2]
        σgrid = vec(A[end,3:end])
        for el in [[s.X, s.Y, s.Z, s.fid,
                    s.linenum, s.rRx, s.zRxLM, s.zTxLM, s.zRxHM,
                    s.zTxHM, s.rTx]; vec(zall); σgrid; ϕd]
            msg = @sprintf("%e\t", el)
            write(io, msg)
        end
        write(io, "\n")
        flush(io) # slower but ensures write is complete
        rmfile && rm(fname)
    end
    close(io)
end

# plot the convergence and the result
function splitlineconvandlast(soundings, delr, delz;
        zstart=-1, extendfrac=-1, dz=-1, nlayers=-1, cmapσ="jet", vmin=-2.5, vmax=0.5, fontsize=12,
        figsize=(20,5),
        topowidth=1,
        preferEright = false,
        preferNright = false,
        saveplot = true,
        showplot = true,
        postfix = "",
        prefix = "",
        markersize = 2,
        logscale = false,
        dpi=400)
    linestartidx = splitsoundingsbyline(soundings)
    nlines = length(linestartidx)
    fnamecheck = soundings[1].sounding_string*"_gradientinv.dat"
    zall, = setupz(zstart, extendfrac, dz=dz, n=nlayers)
    isfile(fnamecheck) && compress(soundings, zall) # write everything in one file if not done yet
    fzipped = prefix == "" ? "zipped.dat" : prefix*"_zipped.dat"
    A = readdlm(fzipped)
    σ = A[:,end-nlayers:end-1]
    ϕd = A[:,end]
    for i in 1:nlines
        a = linestartidx[i]
        b = i != nlines ?  linestartidx[i+1]-1 : length(soundings)
        plotconvandlast(soundings[a:b], view(σ, a:b, :)', view(ϕd, a:b), delr, delz;
            zall = zall,
            cmapσ=cmapσ, vmin=vmin, vmax=vmax, fontsize=fontsize, postfix=postfix, markersize=markersize,
            figsize=figsize, topowidth=topowidth, preferEright=preferEright, logscale=logscale,
            preferNright=preferNright, saveplot=saveplot, showplot=showplot, dpi=dpi)
    end
    getphidhist(ϕd, doplot=true, saveplot=saveplot, prefix=prefix)
    nothing
end

function getphidhist(ϕd; doplot=false, saveplot=false, prefix="", figsize=(8,4), fontsize=11)
    edges=[0,1.1,2, Inf]
    ϕdcounts = fit(Histogram, filter(!isnan, ϕd), edges).weights
    good, bad, ugly = round.(Int, ϕdcounts./sum(ϕdcounts)*100) # % for type of fit in ranges above
    if doplot
        fout = prefix == "" ? "summaryfits.png" : prefix*"_summaryfits.png"
        f = figure(figsize=figsize)
        width=[diff(edges[1:end-1]);1]
        bar(edges[1:end-1], [good, bad, ugly], align="edge", width=width)
        xlabel(L"\phi_d")
        ylabel("% of soundings")
        isgood = sum(.!isnan.(ϕd))
        title(prefix*" total soundings: $isgood NaN: $(length(ϕd)-isgood)")
        grid()
        nicenup(f, fsize=fontsize)
        saveplot && savefig(fout, dpi=400)
    end
    good, bad, ugly
end

function plotconvandlast(soundings, σ, ϕd, delr, delz;
        zall=nothing, cmapσ="jet", vmin=-2.5, vmax=0.5, fontsize=12,
        figsize=(20,5),
        topowidth=1,
        preferEright = false,
        preferNright = false,
        saveplot = true,
        showplot = true,
        logscale = false,
        postfix = "",
        markersize = 2,
        dpi = 400)
    @assert !isnothing(zall)
    Eislast, Nislast = whichislast(soundings)
    # ϕd, σ = readingrid(soundings, zall)
    img, gridr, gridz, topofine, R = makegrid(σ, soundings, zall=zall, dz=delz, dr=delr)
    fig = plt.figure(figsize=figsize)
    axd = fig.subplot_mosaic(
           """
           A
           B
           C
           C
           C
           C
           C
           """
       )
    lname = "Line_$(soundings[1].linenum)"*postfix
    x0, y0 = soundings[1].X, soundings[1].Y
    zTx = [s.zTxLM for s in soundings]
    xend, yend = soundings[end].X, soundings[end].Y
    good, bad, ugly = getphidhist(ϕd)
    fig.suptitle(lname*" Δx=$delr m, Fids: $(length(R)) "*L"\phi_{d_{0-1.1}}:"*" $good "*
    L"\phi_{d_{1.1-2}}:"*" $bad "*L"\phi_{d_{2-\infty}}:"*" $ugly", fontsize=fontsize)
    ax = fig.axes
    ax[1].plot(R, ones(length(R)), "--k")
    ax[1].plot(R, ϕd, ".", markersize=markersize)
    ax[1].set_ylim(0.316, maximum(ax[1].get_ylim()))
    ax[1].set_ylabel(L"\phi_d")
    logscale && ax[1].set_yscale("log")
    ax[2].plot(R, zTx)
    ax[2].set_ylabel("zTx m")
    [a.tick_params(labelbottom=false) for a in ax[1:2]]
    imlast = ax[3].imshow(img, extent=[gridr[1], gridr[end], gridz[end], gridz[1]], cmap=cmapσ, aspect="auto", vmin=vmin, vmax=vmax)
    ax[3].plot(gridr, topofine, linewidth=topowidth, "-k")
    eg = extrema(gridr)
    ax[3].set_ylabel("mAHD")
    ax[3].set_xlabel("Distance m")
    fig.colorbar(imlast, ax=axd["C"], location="bottom", shrink=0.6, label="Log₁₀ S/m")
    nicenup(fig, fsize=fontsize)
    plotNEWSlabels(Eislast, Nislast, gridr, gridz, [ax[3]], x0, y0, xend, yend,
                    preferEright=preferEright, preferNright=preferNright)
    [a.set_xlim(ax[3].get_xlim()) for a in ax[1:2]]
    saveplot && savefig(lname*".png", dpi=dpi)
    showplot || close(fig)
end

# plot multiple grids with supplied labels
function plotgrids()
    f, ax = plt.subplots(size(grids, 1), 1, figsize=figsize,
                        sharex=true, sharey=true, squeeze=false)
    for (i, g) in enumerate(grids)
        ax[i].imshow(g, cmap=cmap, aspect="auto", alpha=α, vmax=vmax, vmin = vmin,
                extent=[gridx[1], gridx[end], gridz[end], gridz[1]])
        colorbar(img, ax=ax[i])
        titles == nothing || ax[i].set_title(titles[i])
        idx  == nothing || plotprofile.(ax[i], idx)
        elev == nothing || plot(gridx, elev, "-k")
        ylabels == nothing || ax[i].set_ylabel(ylabels[i])
    end
    ax[end].set_xlabel("Easting m")
    transD_GP.CommonToAll.nicenup(f, fsize=fsize)
end

# no nuisances e.g., SkyTEM, transD
function loopacrosssoundings(soundings::Array{S, 1}, opt_in::Options;
                            nsequentialiters   =-1,
                            nparallelsoundings =-1,
                            zfixed             = [-1e5],
                            ρfixed             = [1e12],
                            useML              = false,
                            zstart             = 0.0,
                            extendfrac         = 1.06,
                            dz                 = 2.,
                            ρbg                = 10,
                            nlayers            = 50,
                            ntimesperdecade    = 10,
                            nfreqsperdecade    = 5,
                            Tmax               = -1,
                            nsamples           = -1,
                            nchainsatone       = -1,
                            modelprimary       = false,
                            nchainspersounding = -1) where S<:Sounding

    @assert nsequentialiters  != -1
    @assert nparallelsoundings != -1
    @assert nchainspersounding != -1
    @assert nsamples != - 1
    @assert nchainsatone != -1
    @assert Tmax != -1

    nsoundings = length(soundings)
    opt= deepcopy(opt_in)

    for iter = 1:nsequentialiters
        if iter<nsequentialiters
            ss = (iter-1)*nparallelsoundings+1:iter*nparallelsoundings
        else
            ss = (iter-1)*nparallelsoundings+1:nsoundings
        end
        @info "soundings in loop $iter of $nsequentialiters", ss
        r_nothing = Array{Nothing, 1}(undef, length(ss))
        @sync for (i, s) in Iterators.reverse(enumerate(ss))
            pids = (i-1)*nchainspersounding+i:i*(nchainspersounding+1)
            @info "pids in sounding $s:", pids

            aem, = makeoperator(    soundings[s],
                                    zfixed = zfixed,
                                    ρfixed = ρfixed,
                                    zstart = zstart,
                                    extendfrac = extendfrac,
                                    dz = dz,
                                    ρbg = ρbg,
                                    useML = useML,
                                    nlayers = nlayers,
                                    modelprimary = modelprimary,
                                    ntimesperdecade = ntimesperdecade,
                                    nfreqsperdecade = nfreqsperdecade)

            opt = deepcopy(opt_in)
            opt.fdataname = soundings[s].sounding_string*"_"

            @async r_nothing[i] = remotecall_fetch(main, pids[1], opt, aem, collect(pids[2:end]),
                                    Tmax         = Tmax,
                                    nsamples     = nsamples,
                                    nchainsatone = nchainsatone)

        end # @sync
        @info "done $iter out of $nsequentialiters at $(Dates.now())"
    end
end

# for gradient based inversion
function loopacrosssoundings(soundings::Array{S, 1}, σstart, σ0;
                            nsequentialiters   =-1,
                            zfixed             = [-1e5],
                            ρfixed             = [1e12],
                            zstart             = 0.0,
                            extendfrac         = 1.06,
                            dz                 = 2.,
                            ρbg                = 10,
                            nlayers            = 50,
                            ntimesperdecade    = 10,
                            nfreqsperdecade    = 5,
                            modelprimary       = false,
                            regtype            = :R1,
                            nstepsmax          = 10,
                            ntries             = 6,
                            target             = nothing,
                            lo                 = -3.,
                            hi                 = 1.,
                            λ²min              = 0,
                            λ²max              = 8,
                            λ²frac             = 4,
                            β²                 = 0.,
                            ntestdivsλ²        = 50,
                            αmin               = -4,
                            αmax               = 0,
                            αfrac              = 4,
                            ntestdivsα         = 32,
                            regularizeupdate   = false,
                            knownvalue         = 0.7,
                            firstvalue         = :last,
                            κ                  = GP.Mat52(),
                            breakonknown       = true,
                            dobo               = false,
                            compresssoundings  = true,
                            zipsaveprefix      = "",
                            ) where S<:Sounding

    @assert nsequentialiters  != -1
    nparallelsoundings = nworkers()
    nsoundings = length(soundings)
    zall, = setupz(zstart, extendfrac, dz=dz, n=nlayers) # needed for sounding compression
    for iter = 1:nsequentialiters
        if iter<nsequentialiters
            ss = (iter-1)*nparallelsoundings+1:iter*nparallelsoundings
        else
            ss = (iter-1)*nparallelsoundings+1:nsoundings
        end
        @info "soundings in loop $iter of $nsequentialiters", ss
        pids = workers()
        @sync for (i, s) in enumerate(ss)
            aem, = makeoperator(    soundings[s],
                                    zfixed = zfixed,
                                    ρfixed = ρfixed,
                                    zstart = zstart,
                                    extendfrac = extendfrac,
                                    dz = dz,
                                    ρbg = ρbg,
                                    nlayers = nlayers,
                                    modelprimary = modelprimary,
                                    ntimesperdecade = ntimesperdecade,
                                    nfreqsperdecade = nfreqsperdecade,
                                    calcjacobian = true)

            fname = soundings[s].sounding_string*"_gradientinv.dat"
            σstart_, σ0_ = map(x->x*ones(length(aem.ρ)-1), [σstart, σ0])
            @async remotecall_wait(gradientinv, pids[i], σstart_, σ0_, aem,
                                                regtype            = regtype         ,
                                                nstepsmax          = nstepsmax       ,
                                                ntries             = ntries          ,
                                                target             = target          ,
                                                lo                 = lo              ,
                                                hi                 = hi              ,
                                                λ²min              = λ²min           ,
                                                λ²max              = λ²max           ,
                                                λ²frac             = λ²frac          ,
                                                ntestdivsλ²        = ntestdivsλ²     ,
                                                αmin               = αmin            ,
                                                αmax               = αmax            ,
                                                αfrac              = αfrac           ,
                                                β²                 = β²              ,
                                                ntestdivsα         = ntestdivsα      ,
                                                regularizeupdate   = regularizeupdate,
                                                knownvalue         = knownvalue      ,
                                                firstvalue         = firstvalue      ,
                                                κ                  = κ               ,
                                                breakonknown       = breakonknown    ,
                                                dobo               = dobo            ,
                                                fname              = fname           )

        end # @sync
        isfirstparalleliteration = iter == 1 ? true : false
        compresssoundings && compress(soundings[ss[1]:ss[end]], zall,
            isfirstparalleliteration = isfirstparalleliteration, prefix=zipsaveprefix)
        @info "done $iter out of $nsequentialiters at $(Dates.now())"
    end
end

end

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA

max_forks_repo_name: GeoscienceAustralia/HiQGA


text:
# runtests.jl

include("../src/infrastructure/infrastructure.jl")

using Test

using DataFrames

# TEST MODEL
@testset "Orders" begin
    using AppliSales
    orders = AppliSales.process()
    @test length(orders) == 3
    @test orders[1].org.name == "Scrooge Investment Bank"
    @test orders[1].training.name == "Learn Smiling"
end

@testset "UnpaidInvoices" begin
    invnbr = 1000
    path = "./test_invoicing.sqlite"
    db = connect(path)
    using AppliSales
    orders = AppliSales.process()
    invoices_to_save = [create(order, "A" * string(invnbr += 1)) for order in orders]
    archive(db, "UNPAID", invoices_to_save)
    invoices = retrieve_unpaid_invoices(path=path)
    #invoices = SQLite.DBInterface.execute(db, "select * from UNPAID") |> DataFrame

    @test id(invoices[1]) == "A1001"
    @test currency_ratio(meta(invoices[1])) == 1.0
    @test name(header(invoices[1])) == "Scrooge Investment Bank"
    @test price_per_student(body(invoices[1])) == 1000.0
    @test length(students(body(invoices[1]))) == 1
    @test students(body(invoices[1]))[1] == "Scrooge McDuck"
    @test vat_perc(body(invoices[1])) == 0.21

    cmd = `rm test_invoicing.sqlite`
    run(cmd)
end

@testset "Retrieve UnpaidInvoices" begin
    using AppliSales
    orders = AppliSales.process()
    process(orders)
    unpaid_invoices = retrieve_unpaid_invoices()

    @test id(unpaid_invoices[1]) == "A1001"
    @test currency_ratio(meta(unpaid_invoices[1])) == 1.0
    @test name(header(unpaid_invoices[1])) == "Scrooge Investment Bank"
    @test length(students(body(unpaid_invoices[1]))) == 1
    @test price_per_student(body(unpaid_invoices[1])) == 1000.0
    @test first(students(body(unpaid_invoices[1]))) == "Scrooge McDuck"
    @test vat_perc(body(unpaid_invoices[1])) == 0.21

    cmd = `rm test_invoicing.sqlite`
    run(cmd)
end

@testset "Retrieve BankStatement from CSV" begin
    #stms = read_bank_statements("./bank.csv")
    stms = [BankStatement(Date(2020-01-15), "Duck City Chronicals Invoice A1002", "NL93INGB", 2420.0)]
    @test length(stms) == 1
    #@test length(stms) == 2
    @test amount(stms[1]) == 2420.0
    #@test stms[2].amount == 1210.0
end

@testset "JounalEntry's" begin
    stm1 = BankStatement(Date(2020-01-15), "Duck City Chronicals Invoice A1002", "NL93INGB", 2420.0)
    stms = [stm1]

    orders = AppliSales.process()
    process(orders)

    invoices = retrieve_unpaid_invoices()

    potential_paid_invoices = []
    for unpaid_invoice in invoices
      for s in stms # get potential paid invoices
        if occursin(id(unpaid_invoice), descr(s)) # description contains invoice number
          push!(potential_paid_invoices, create(unpaid_invoice, s))
        end
      end
    end

    @test length(potential_paid_invoices) == 1
    @test id(potential_paid_invoices[1]) == "A1002"
    @test amount(stm((potential_paid_invoices[1]))) == 2420.0

    cmd = `rm test_invoicing.sqlite`
    run(cmd)
end

@testset "process(db, orders)" begin
    orders = AppliSales.process()
    entries = process(orders)
    @test length(entries) == 3
    @test entries[1].from == 1300
    @test entries[1].to == 8000
    @test entries[1].debit == 1000.0
    @test entries[1].vat == 210.0
    @test entries[1].descr == "Learn Smiling"

    cmd = `rm test_invoicing.sqlite`
    run(cmd)
end

@testset "retrieve_unpaid_invoices" begin
    path = "./test_invoicing.sqlite"
    orders = AppliSales.process()
    entries = process(orders; path=path)
    unpaid_invoices = retrieve_unpaid_invoices(path=path)

    @test length(unpaid_invoices) == 3
    @test id(unpaid_invoices[1]) == "A1001"

    cmd = `rm test_invoicing.sqlite`
    run(cmd)
end

@testset "process(unpaid_invoices)" begin
    path = "./test_invoicing.sqlite"
    orders = AppliSales.process()
    process(orders, path=path)
    unpaid_invoices = retrieve_unpaid_invoices(path=path)

    stm1 = BankStatement(Date(2020-01-15), "Duck City Chronicals Invoice A1002", "NL93INGB", 2420.0)
    stms = [stm1]
    entries = process(unpaid_invoices, stms; path=path)
    @test length(entries) == 1
    @test entries[1].from == 1150
    @test entries[1].to == 1300
    @test entries[1].debit == 2420.0

    cmd = `rm test_invoicing.sqlite`
    run(cmd)
end

@testset "disconnect(db)" begin
    path = "./test_invoicing.sqlite"
    db = connect(path)
    orders = AppliSales.process()
    process(orders; path=path)
    disconnect(db)
    try
        retrieve(db, "UNPAID")
    catch e
        @test(typeof(e) == SQLite.SQLiteException)
    end

    cmd = `rm test_invoicing.sqlite`
    run(cmd)
end

ext: jl


text:
using Iterators, Base.Test

# count
# -----

i = 0
for j = count(0, 2)
	@test j == i*2
	i += 1
	i <= 10 || break
end

# take
# ----

i = 0
for j = take(0:2:8, 10)
	@test j == i*2
	i += 1
end
@test i == 5

i = 0
for j = take(0:2:100, 10)
	@test j == i*2
	i += 1
end
@test i == 10

# drop
# ----

i = 0
for j = drop(0:2:10, 2)
	@test j == (i+2)*2
	i += 1
end
@test i == 4

# cycle
# -----

i = 0
for j = cycle(0:3)
	@test j == i % 4
	i += 1
	i <= 10 || break
end

# repeated
# --------

i = 0
for j = repeated(1, 10)
	@test j == 1
	i += 1
end
@test i == 10
i = 0
for j = repeated(1)
	@test j == 1
	i += 1
	i <= 10 || break
end

# chain
# -----

@test collect(chain(1:2:5, 0.2:0.1:1.6)) == [1:2:5, 0.2:0.1:1.6]

# product
# -------

x1 = 1:2:10
x2 = 1:5
@test collect(product(x1, x2)) == vec([(y1, y2) for y1 in x1, y2 in x2])

# distinct
# --------

x = [5, 2, 2, 1, 2, 1, 1, 2, 4, 2]
@test collect(distinct(x)) == unique(x)

# partition
# ---------

@test collect(partition(take(count(1), 6), 2)) == [(1,2), (3,4), (5,6)]
@test collect(partition(take(count(1), 4), 2, 1)) == [(1,2), (2,3), (3,4)]
@test collect(partition(take(count(1), 8), 2, 3)) == [(1,2), (4,5), (7,8)]

# imap
# ----

function test_imap(expected, input...)
  result = collect(imap(+, input...))
  @test result == expected
end

# Empty arrays
test_imap(
  {},
  []
)

# Simple operation
test_imap(
  {1,2,3},
  [1,2,3]
)

# Multiple arguments
test_imap(
  {5,7,9},
  [1,2,3],
  [4,5,6]
)

# Different-length arguments
test_imap(
  {2,4,6},
  [1,2,3],
  count(1)
)

# groupby
# -------

function test_groupby(input, expected)
  result = collect(groupby(input, x -> x[1]))
  @test result == expected
end

# Empty arrays
test_groupby(
  [],
  {}
)

# Singletons
test_groupby(
  ["xxx"],
  {["xxx"]}
)

# Typical operation
test_groupby(
  ["face", "foo", "bar", "book", "baz"],
  {["face", "foo"], ["bar", "book", "baz"]}
)

# Trailing singletons
test_groupby(
  ["face", "foo", "bar", "book", "baz", "xxx"],
  {["face", "foo"], ["bar", "book", "baz"], ["xxx"]}
)

# Leading singletons
test_groupby(
  ["xxx", "face", "foo", "bar", "book", "baz"],
  {["xxx"], ["face", "foo"], ["bar", "book", "baz"]}
)

# Middle singletons
test_groupby(
  ["face", "foo", "xxx", "bar", "book", "baz"],
  {["face", "foo"], ["xxx"], ["bar", "book", "baz"]}
)

# subsets

@test collect(subsets({})) == {{}}

@test collect(subsets([:a])) == {Symbol[], Symbol[:a]}

@test collect(subsets([:a, :b, :c])) ==
      {Symbol[], Symbol[:a], Symbol[:b], Symbol[:a, :b], Symbol[:c],
       Symbol[:a, :c], Symbol[:b, :c], Symbol[:a, :b, :c]}

# every n'th element
@test collect(everynth([], 10)) == []
@test_throws ArgumentError everynth([], 0)
@test collect(everynth(10:20, 3)) == [12,15,18]
@test collect(everynth(10:20, 1)) == [10:20]


text:
module TestInferrability

using Test

include("../../../examples/inferrability.jl")

should_test_module() = lowercase(get(ENV, "JULIA_PKGEVAL", "false")) != "true"

function test()
    @test @inferred(UnionTyped.f((111,))) == Some(111)
    @test @inferred(UnionTyped.f(())) === nothing
    @test ConcretelyTyped.f((111,)) == Some(111)
    @test ConcretelyTyped.f(()) === nothing
end

end  # module


text:
int main() {
    if(false) return 0;
}

max_stars_repo_name: carlostome/Javalette-Compiler

max_issues_repo_name: carlostome/Javalette-Compiler

max_forks_repo_name: carlostome/Javalette-Compiler


text:
# use Convex to get access to its types and functions
using Convex

# import reshape to extend it.
import Base.reshape

export partialtransposematrix, partialtranspose, partialtransposenaive

# switch between coordinates of a matrix and its vectorised form
vecCoordToMatrixCoord(k::Integer, n::Integer) = (1 + (k-1)%n, round(Int, 1+floor((k-1)/n)%n))
matrixCoordToVecCoord(i::Integer, j::Integer, n::Integer) = i + (j-1)*n

#Determine the new coordinates `(x, y)` of an entry
#with coordinates `(i, j)` after partial transposition.

function partialtranseposecoord(i::Integer, j::Integer, n::Integer, l::Integer)
  @assert n % l == 0 "l should divide n"
  @assert n ≥ l "n should be greater than l"

  # calculate which block (i, j) is in
  a = floor(Int, (i-1)/l)+1
  b = floor(Int, (j-1)/l)+1

  x = j + (a-b)*l
  y = i + (b-a)*l

  @assert x > 0
  @assert y > 0
  @assert n ≥ x
  @assert n ≥ y

  return (x, y)
end

#Construct a permutation matrix for a `n × n` matrix
#in vector form where we transpose in `l × l` blocks.

function partialtransposematrix(n::Integer, l::Integer)
  @assert n % l == 0 "l should divide n"
  @assert n ≥ l "n should be greater than l"

  # Julia will iterate over I and J simulatiously
  # and set (I[k], J[k]) to 1 in the returned matrix
  # The value in I corresponds to the vectorised coordinates
  # before partial transposition, whereas the value in J
  # corresponds to the coordinates after partial transposition.
  IM = Array(Int, n^2)
  JM = Array(Int, n^2)
  index = 1
  for i in 1:n
    for j in 1:n
      (x, y) = partialtranseposecoord(i, j, n, l)
      IM[index] = matrixCoordToVecCoord(i, j, n)
      JM[index] = matrixCoordToVecCoord(x, y, n)
      index += 1
    end
  end

  return sparse(IM, JM, 1)
end

""" `rhoOut = partialtranspose(rho)` or `rhoOut = partialtranspose(rho, sys)` or `rhoOut = partialtranspose(rho, sys, dim)`

Applies the partial tanspose to the matrix *rho*.

Inputs:
- *rho* input matrix
- *systems* (optional) single system or a vector of systems to be transposed. By default it is equal to 2. If *dim* is not provided, the only two possible values of *systems* are 1 or 2.
- *dim* (optional) vector with dimensions of the systems. By default it assumes two systems of equal dimensions.

Outputs:
- *rhoOut* matrix *rho* with partial transpose applied
"""
function partialtranspose(ρ::Union{AbstractArray, AbstractExpr}, systems::Vector, dim::Vector)
  result = ρ
  for sys in systems
    result = partialtranspose(result, sys, dim)
  end
  return result
end

#Apply the partial tanspose to a single sytem. If no
#dimensions are specified, it will assume there are
#two systems of identical dimension. If no system is
#specified, it will transpose the last system. This
#will always be the second system, because it's not
#possible to leave out the system while providing
#the dimensions.

function partialtransposesystem(ρ::Union{AbstractArray, AbstractExpr}, sys::Integer = 2, dim::Vector = [])
  # do naive partial transposition
  if dim == []
    l = round(Int, sqrt(size(ρ)[1]))
    if sys == 2
      return partialtranspose(ρ, l, naive = true)
    else
      dim = [l; l]
    end
  end

  @assert 1 ≤ sys ≤ length(dim) "System doesnt exist, expected sys to be between 1 and ", length(dim), ". Is your dim vector correct?"

  # permute systems such that the system to apply PT is last
  perm = collect(1:length(dim))
  deleteat!(perm, sys)
  push!(perm, sys)

  l = dim[sys]
  x = permutesystems(ρ, perm, dim)
  y = partialtransposenaive(x, l)
  z = permutesystems(y, perm, dim)
  return z
end

#Apply partial transpose on systems if naive is false.
#Otherwise do the naive partial transpose in blocks
#of `l × l`, where `l` takes the value `systems`.
#This is equivalent to applying the partial
#transepose to the last system if that particular system
#has dimension `l`.

function partialtranspose(ρ::Union{AbstractArray, AbstractExpr}, systems::Integer = 2, dim::Vector = []; naive::Bool = false)
  if naive
    return partialtransposenaive(ρ, systems)
  end

  return partialtransposesystem(ρ, systems, dim)
end

#This does the naive version of the partial transpose
#that I wrote first it simply applies a transpose in blocks
#of `l × l`. This is equivalent to applying the partial
#transepose to the last system if that particular system
#has dimension `l`.

function partialtransposenaive(ρ::Union{AbstractExpr, AbstractArray, Convex.Variable}, l::Integer)
  # if the input is sparse, keep the output sparse as well
  if !isa(ρ, Convex.Variable) && !isa(ρ, AbstractExpr) && issparse(ρ)
    vectorised = sparsevec(ρ)
  else
    vectorised = vec(ρ)
  end
  n = size(ρ)[1]

  @assert n % l == 0 "l should divide n"
  @assert n ≥ l "n should be greater than l"

  return reshape(partialtransposematrix(n, l) * vectorised, n, n)
end

# some aliases for convenience
matrixCoordToVecCoord(ij::Tuple{Int, Int}, n::Integer) = matrixCoordToVecCoord(ij[1], ij[2], n)
partialtranseposecoord(ij::Tuple{Int, Int}, n::Integer, l::Integer) = partialtranseposecoord(ij[1], ij[2], n, l)
reshape(x::AbstractExpr, mn::Tuple{Int, Int}) = reshape(x, mn[1], mn[2])


text:
module Tensor
    using Einsum
    export ⊕,⊗,⋆
    function ⊕(A::T,B::T) where T <: AbstractMatrix
        @einsum C[i,j,k,l] := A[k,j] * B[i,l]
    end

    function oplus!(C::Array{T,4},A::Array{T,2},B::Array{T,2}) where T <: Number
        K,J = size(A)
        I,L = size(B)
        for l=1:L
            for k=1:K
                for j=1:J
                    for i=1:I
                        C[i,j,k,l] = A[k,j]*B[i,l]
                    end
                end
            end
        end
    end

    function otimes!(C::Array{T,4},A::Array{T,2},B::Array{T,2}) where T <: Number
        I,J = size(A)
        K,L = size(B)
        for l=1:L
            for k=1:K
                for j=1:J
                    for i=1:I
                        C[i,j,k,l] = A[i,j]*B[k,l]
                    end
                end
            end
        end
    end

    function ⊗(A::T,B::T) where T <: AbstractMatrix
        @einsum C[i,j,k,l] := A[i,j] * B[k,l]
    end

    function ⋆(A::Array{T1,2},B::Array{T2,4}) where {T1 <: Number,T2 <: Number}
        @einsum C[i,l] := A[j,k] * B[i,j,k,l]
    end

    function star!(C::Array{T3,2},A::Array{T1,2},B::Array{T2,4}) where {T1 <: Number,T2 <: Number,T3 <: Number}
        J,K = size(A)
        I,JB,KB,L = size(B)
        @assert J == JB "dimension mismatch! J of A = $J and J of B = $JB"
        @assert K == KB "dimension mismatch! K of A = $K and K of B = $KB"
        @. C = 0
        for l=1:L
            for k=1:K
                for j=1:J
                    for i=1:I
                        C[i,l] += A[j,k]*B[i,j,k,l]
                    end
                end
            end
        end
    end

    function star!(C::Array{T3,4},A::Array{T1,4},B::Array{T2,4}) where {T1 <: Number,T2 <: Number,T3 <: Number}
        I,M,N,L = size(A)
        MB,J,K,NB = size(B)
        @assert N == NB "dimension mismatch! N of A = $N and N of B = $NB"
        @assert M == MB "dimension mismatch! M of A = $M and M of B = $MB"
        @. C = 0
        for l=1:L
            for k=1:K
                for j=1:J
                    for i=1:I
                        for n=1:N
                            for m=1:M
                                C[i,j,k,l] += A[i,m,n,l]*B[m,j,k,n]
                            end
                        end
                    end
                end
            end
        end
    end

    function ⋆(A::Array{T1,4},B::Array{T2,4}) where {T1 <: Number,T2 <: Number}
        @einsum C[i,j,k,l] := A[i,m,n,l] * B[m,j,k,n]
    end

    const I0 = [1 0;0 1]
    const IoplusI = I0⊕I0
    const IotimesI = I0⊗I0
end

#=
using ..Tensor
A = rand(2,2)
B = rand(2,2)
C1 = A ⊕ B
display(C1)
println("")
Tensor.oplus!(C1,A,B)
display(C1)
println("")
#exit()

A = rand(2,2)
B = rand(2,2)
C = A ⊗ B
display(C)
println("")
C1 = copy(C)
Tensor.otimes!(C1,A,B)
display(C1)
println("C1")
#exit()

println("")
B = rand(2,2)
C = A ⊗ B
D = B ⋆ C
display(D)
println("\tD")
D1 = copy(D)
Tensor.star!(D1,B,C)
display(D1)
println("")

println("")
E = C1 ⋆ C
display(E)

println("\tD")
E1 = copy(E)
Tensor.star!(E1,C1,C)
display(E1)
println("")

Tensor.star!(E1,Tensor.IoplusI,C)
display(E1)
println("")

=#


text:
#=
Source:
  S. Ahmed and R. Garcia. "Dynamic Capacity Acquisition and Assignment under Uncertainty," Annals of Operations Research, vol.124, pp. 267-283, 2003

Input:
  nR: number of resources
  nN: number of tasks
  nT: number of time periods
  nS: number of scenarios

Sets:
  R: resources
  N: tasks
  T: time periods
  S: scenarios

Variables (1st Stage):
  x[i,t]: capacity acquired for resource i at period t
  u[i,t]: 1 if x[i,t] > 0, 0 otherwise

Variables (2nd Stage):
  y[i,j,t]: 1 if resource i is assigned to task j in period t, 0 otherwise

Parameters (general):
  a[i,t]: linear component of expansion cost for resource i at period t
  b[i,t]: fixed component of expansion cost for resource i at period t
  c[i,j,t,s]: cost of assigning resource i to task j in period t
  c0[j,t,s]: penalty incurred if task j in period t is not served

Parameters (scenario):
  d[j,t,s]: capacity required for to perform task j in period t in scenario s
=#

using StructJuMP
using Random
using MathOptInterface

function DCAP(nR::Int, nN::Int, nT::Int, nS::Int, seed::Int=1)::StructuredModel

    # set random seed (default=1)
    Random.seed!(seed)

    # generate & store instance data
    ## sets
    R = 1:nR
    N = 1:nN
    T = 1:nT
    S = 1:nS

    ## parameters
    a = rand(nR, nT) * 5 .+ 5
    b = rand(nR, nT) * 40 .+ 10
    c = rand(nR, nN, nT, nS) * 5 .+ 5
    c0 = rand(nN, nT, nS) * 500 .+ 500
    d = rand(nN, nT, nS) .+ 0.5
    Pr = ones(nS)/nS

    # construct JuMP.Model
    model = StructuredModel(num_scenarios = nS)

    ## 1st stage
    @variable(model, x[i=R,t=T] >= 0)
    @variable(model, u[i=R,t=T], Bin)
    @objective(model, Min, sum(a[i,t]*x[i,t] + b[i,t]*u[i,t] for i in R for t in T))
    @constraint(model, [i=R,t=T], x[i,t] - u[i,t] <= 0)

    ## 2nd stage
    for s in S
        sb = StructuredModel(parent=model, id = s, prob = Pr[s])
        @variable(sb, y[i=R, j=N, t=T], Bin)
        #@variable(sb, z[j=N,t=T] >= 0) # originally implemented variable (continuous)
        @variable(sb, z[j=N,t=T], Bin)  # modify as SIPLIB 1.0
        @objective(sb, Min, sum(c[i,j,t,s]*y[i,j,t] for i in R for j in N for t in T) + sum(c0[j,t,s]*z[j,t] for j in N for t in T))
        @constraint(sb, [i=R, t=T], -sum(x[i,tau] for tau in 1:t) + sum(d[j,t,s]*y[i,j,t] for j in N) <= 0)
        @constraint(sb, [j=N, t=T], sum(y[i,j,t] for i in R) + z[j,t] == 1)
    end

    return model
end

m = DCAP(2,2,3,3);
DSPopt.classifyConstrs(m)

@testset "Parent model" begin
  @test length(m.variables) == 12
  @test length(m.constraints) == 6
  start, index, value, rlbd, rubd, obj, clbd, cubd, ctype, cname, nqrows, linnzcnt, quadnzcnt, rhs, sense, linstart, linind, linval, quadstart, quadrow, quadcol, quadval = DSPopt.get_model_data(m)
  @show clbd
  @show cubd
  @show ctype
  # print(m)
end

@testset "Child model $i" for (i,subm) in m.children
  start, index, value, rlbd, rubd, obj, clbd, cubd, ctype, cname, nqrows, linnzcnt, quadnzcnt, rhs, sense, linstart, linind, linval, quadstart, quadrow, quadcol, quadval = DSPopt.get_model_data(subm, i)
  @test clbd == zeros(Float64, 18)
  @test cubd == ones(Float64, 18)
  @test ctype == "BBBBBBBBBBBBBBBBBB"
end

@testset "optimize!: $j" for j in [DSPopt.ExtensiveForm]
  status = DSPopt.optimize!(m, solve_type = j, is_stochastic = true, param = "params.txt")
  @test status == MathOptInterface.OPTIMAL
end


text:
using HTTP
import ToStruct
import JSON

response = String(
    HTTP.request(
        "POST",
        "http://localhost:8000/jsonpayload",
        [("Content-Type", "application/json")],
        """{"command":"showsection"}""",
    ).body,
)
println(response)
result = ToStruct.tostruct(ContractSection, JSON.parse(response))
println(result)

max_stars_repo_name: michaelfliegner/BitempWebApp

max_issues_repo_name: michaelfliegner/BitempWebApp


text:
# Functions to create interval graphs in Julia

export IntervalGraph, IntervalDigraph1, IntervalDigraph2
export RandomIntervalGraph, RandomIntervalDigraph1, RandomIntervalDigraph2
export UnitIntervalGraph, RandomUnitIntervalGraph
export UnitIntervalGraphEvolution

"""
`IntervalGraph(Jlist)` creates an interval graph from a list of closed
intervals.
"""
function IntervalGraph(Jlist::Array{ClosedInterval{T},1}) where {T}
    d = Dict{Int,ClosedInterval{T}}()
    n = length(Jlist)
    for k = 1:n
        d[k] = Jlist[k]
    end
    return IntervalGraph(d)
end

"""
`IntervalGraph(f::Dict)` creates asn interval graph from a dictionary
whose keys are the names of the vertices and whose values are
intervals.
"""
function IntervalGraph(f::Dict{K,ClosedInterval{T}}) where {K,T}
    klist = collect(keys(f))
    G = SimpleGraph{K}()
    for v in klist
        add!(G, v)
    end

    n = length(klist)

    for i = 1:n-1
        u = klist[i]
        for j = i+1:n
            v = klist[j]
            if !isempty(f[u] * f[v])
                add!(G, u, v)
            end
        end
    end
    cache_save(G, :IntervalRepresentation, f)
    cache_save(G, :name, "Interval graph")
    return G
end

"""
`edge_helper(A,B)` for closed intervals `A` and `B` is a private
function that returns the following values:

+ `0` if no edge
+ `1` if edge from `A` to `B` (but not reverse)
+ `2` if edge from `B` to `A` (but not reverse)
+ `3` if both
"""
function edge_helper(A::ClosedInterval, B::ClosedInterval)
    # no intersection
    if isempty(A * B)
        return 0
    end

    # containment
    if (left(A) <= left(B) <= right(A)) && (left(A) <= right(B) <= right(A))
        return 3
    end

    if (left(B) <= left(A) <= right(B)) && (left(B) <= right(A) <= right(B))
        return 3
    end

    # overlap
    if left(A) <= left(B)
        return 1
    end
    return 2
end

"""
`IntervalDigraph1(Jlist)` creates a type I interval digraph from a
list of closed intervals.
"""
function IntervalDigraph1(Jlist::Array{ClosedInterval{T},1}) where {T}
    n = length(Jlist)
    G = IntDigraph(n)
    for u = 1:n
        J = Jlist[u]
        for v = u:n
            K = Jlist[v]

            choice = edge_helper(J, K)

            if choice == 1 || choice == 3
                add!(G, u, v)
            end

            if choice == 2 || choice == 3
                add!(G, v, u)
            end

        end
    end
    return G
end

"""
`IntervalDigraph2(send_list, rec_list)` creates a Type II interval
graph from two lists of intervals.
"""
function IntervalDigraph2(
    send_list::Array{ClosedInterval{T},1},
    rec_list::Array{ClosedInterval{T},1},
) where {T}
    n = length(send_list)
    if length(rec_list) != n
        error("send_list and rec_list must be same length")
    end

    G = IntDigraph(n)

    for u = 1:n
        A = send_list[u]
        for v = 1:n
            B = rec_list[v]
            if !isempty(A * B)
                add!(G, u, v)
            end
        end
    end
    return G
end

"""
`IntervalDigraph(f::Dict)` creates a type I interval digraph from a
dictionary whose keys are the names of the vertices and whose values
are intervals.
"""
function IntervalDigraph1(
    snd::Dict{K,ClosedInterval{T}},
    rec::Dict{K,ClosedInterval{T}},
) where {K,T}
    if length(snd) != length(rec)
        error("send and receive dictionaries must have same keys")
    end

    klist = keys(snd)
    for k in klist
        if !haskey(rec)
            error("send and receive dictionaries must have same keys")
        end
    end

    n = length(klist)
    G = SimpleDigraph{K}()

    for i = 1:n
        u = klist[i]
        A = snd[u]
        for j = 1:n
            v = klist[j]
            B = rec[v]

            choice = edge_helper(A, B)

            if choice == 1 || choice == 3
                add!(G, u, v)
            end

            if choice == 2 || choice == 3
                add!(G, v, u)
            end
        end
    end
    return G
end

"""
`IntervalDigraph2(send::Dict,rec::Dict)` creates a Type II interval
digraph from two dictionaries mapping vertices to intervals.
"""
function IntervalDigraph2(
    snd::Dict{K,ClosedInterval{T}},
    rec::Dict{K,ClosedInterval{T}},
) where {K,T}
    if length(snd) != length(rec)
        error("send and receive dictionaries must have same keys")
    end

    klist = keys(snd)
    for k in klist
        if !haskey(rec)
            error("send and receive dictionaries must have same keys")
        end
    end

    n = length(klist)
    G = SimpleGraph{K}()

    for i = 1:n
        u = klist[i]
        A = snd[u]
        for j = 1:n
            v = klist[j]
            B = rec[v]
            if !isempty(A * B)
                add!(G, u, v)
            end
        end
    end

    return G
end

"""
`RandomIntervalGraph(n)` generates a random interval graph with `n`
vertices.
"""
function RandomIntervalGraph(n::Int)
    Jlist = [ClosedInterval(rand(), rand()) for _ = 1:n]
    return IntervalGraph(Jlist)
end

"""
`RandomIntervalDigraph1(n)` generates a random type I interval
digraph with `n` vertices.
"""
function RandomIntervalDigraph1(n::Int)
    Jlist = [ClosedInterval(rand(), rand()) for _ = 1:n]
    return IntervalDigraph1(Jlist)
end

"""
`RandomIntervalDigraph2(n::Int)` generates a random type II
directed interval graph with `n` vertices.
"""
function RandomIntervalDigraph2(n::Int)
    snd_list = [ClosedInterval(rand(), rand()) for _ = 1:n]
    rec_list = [ClosedInterval(rand(), rand()) for _ = 1:n]
    return IntervalDigraph2(snd_list, rec_list)
end

"""
`UnitIntervalGraph(x,t=1)` creates a unit interval graph where the
vector `x` specifies the left end points of the intervals. The
optional parameter `t` specifies the lengths of the intervals.
"""
function UnitIntervalGraph(points::Vector{T}, t::Real = 1) where {T<:Real}
    n = length(points)
    f = Dict{Int,T}()
    for k = 1:n
        f[k] = points[k]
    end
    return UnitIntervalGraph(f, t)
end

"""
`UnitIntervalGraph(f,t=1)` creates a unit interval graph from a
dictionary mapping vertex names to the left end points of their
intervals. The optional parameter `t` specifies the length of the
intervals.
"""
function UnitIntervalGraph(f::Dict{S,T}, t::Real = 1) where {S,T<:Real}
    vtcs = collect(keys(f))
    n = length(vtcs)
    G = SimpleGraph{S}()
    for v in vtcs
        add!(G, v)
    end

    for i = 1:n-1
        u = vtcs[i]
        for j = i+1:n
            v = vtcs[j]
            if abs(f[u] - f[v]) <= t
                add!(G, u, v)
            end
        end
    end
    cache_save(G, :UnitIntervalRepresentation, (f, t))
    cache_save(G, :name, "Unit interval graph")
    return G
end

"""
`RandomUnitIntervalGraph(n,t)` creates a unit interval graph with `n`
vertices whose left end points are chosen iid uniformly from [0,1] and
whose lengths are all given by `t` (default value 0.5).
"""
function RandomUnitIntervalGraph(n::Int, t::Real = 0.5)
    x = rand(n)
    return UnitIntervalGraph(x, t)
end

"""
`UnitIntervalGraphEvolution(points)` gives the sequence of edges as
the lengths of the intervals (whose left end points are specified in
`points`) increases.

This returns a pair consisting of the sequence of edges and the
lengths at which those edges appear.
"""
function UnitIntervalGraphEvolution(points::Vector{S}) where {S<:Real}
    n = length(points)
    nC2 = round(Int, n * (n - 1) / 2)

    edges = Vector{Tuple{Int,Int}}(nC2)
    diffs = Vector{S}(nC2)

    idx = 0

    for i = 1:n-1
        for j = i+1:n
            idx += 1
            edges[idx] = (i, j)
            diffs[idx] = abs(points[i] - points[j])
        end
    end

    p = sortperm(diffs)

    return (edges[p], diffs[p])

end

"""
`UnitIntervalGraphEvolution(n::Int)` is equivalent to
`UnitIntervalGraphEvolution(sort(rand(n)))`.
"""
UnitIntervalGraphEvolution(n::Int) = UnitIntervalGraphEvolution(sort(rand(n)))


text:
@testset "Special operators" begin
    @testset "Projection" begin
        s1 = ParameterSpace() × (Taylor(1) ⊗ Fourier(1, 1.0) ⊗ Chebyshev(1))^1
        s2 = ParameterSpace() × (Taylor(2) ⊗ Fourier(1, 1.0) ⊗ Chebyshev(1))^1
        s3 = ParameterSpace() × (Taylor(1) ⊗ Fourier(1, 1.0) ⊗ Chebyshev(0))^1
        s4 = ParameterSpace() × (Taylor(2) ⊗ Fourier(1, 1.0) ⊗ Chebyshev(0))^1
        a = Sequence(s1, [1.0 ; 1.0:12.0])
        @test project(a, s1) == a
        @test project(a, s2) == Sequence(s2, [1.0, 1.0, 2.0, 0.0, 3.0, 4.0, 0.0, 5.0, 6.0, 0.0, 7.0, 8.0, 0.0, 9.0, 10.0, 0.0, 11.0, 12.0, 0.0])
        @test project(a, s3) == Sequence(s3, [1.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0])
        @test project(a, s4) == Sequence(s4, [1.0, 1.0, 2.0, 0.0, 3.0, 4.0, 0.0, 5.0, 6.0, 0.0])
    end

    @testset "Multiplication" begin
        a_𝒯 = Sequence(Taylor(2), [1.0, -1.0, 1.0])
        ℳ_𝒯 = Multiplication(a_𝒯)
        @test (2ℳ_𝒯/1).sequence == (ℳ_𝒯 + ℳ_𝒯).sequence == (ℳ_𝒯 - (-ℳ_𝒯)).sequence == Multiplication(a_𝒯 + a_𝒯).sequence
        @test (ℳ_𝒯^2).sequence == (ℳ_𝒯 * ℳ_𝒯).sequence == Multiplication(a_𝒯 * a_𝒯).sequence
        @test ℳ_𝒯(a_𝒯) == a_𝒯*a_𝒯 == project(ℳ_𝒯, Taylor(2), Taylor(4), Float64)*a_𝒯
        #
        a_ℱ = Sequence(Fourier(1, 1.0), [0.5, 0.0, 0.5])
        ℳ_ℱ = Multiplication(a_ℱ)
        @test (2ℳ_ℱ/1).sequence == (ℳ_ℱ + ℳ_ℱ).sequence == (ℳ_ℱ - (-ℳ_ℱ)).sequence == Multiplication(a_ℱ + a_ℱ).sequence
        @test (ℳ_ℱ^2).sequence == (ℳ_ℱ * ℳ_ℱ).sequence == Multiplication(a_ℱ * a_ℱ).sequence
        @test ℳ_ℱ(a_ℱ) == a_ℱ*a_ℱ == project(ℳ_ℱ, Fourier(1, 1.0), Fourier(2, 1.0), Float64)*a_ℱ
        #
        a_𝒞 = Sequence(Chebyshev(2), [1.0, 0.5, 0.5])
        ℳ_𝒞 = Multiplication(a_𝒞)
        @test (2ℳ_𝒞/1).sequence == (ℳ_𝒞 + ℳ_𝒞).sequence == (ℳ_𝒞 - (-ℳ_𝒞)).sequence == Multiplication(a_𝒞 + a_𝒞).sequence
        @test (ℳ_𝒞^2).sequence == (ℳ_𝒞 * ℳ_𝒞).sequence == Multiplication(a_𝒞 * a_𝒞).sequence
        @test ℳ_𝒞(a_𝒞) == a_𝒞*a_𝒞 == project(ℳ_𝒞, Chebyshev(2), Chebyshev(4), Float64)*a_𝒞
    end

    @testset "Derivative" begin
        ∂¹ = Derivative(1)
        ∂⁴ = Derivative(4)
        @test ∂¹*∂¹ == ∂¹^2 == Derivative(2)
        #
        a_𝒯 = Sequence(Taylor(2), [1.0, -1.0, 1.0])
        @test ∂¹(a_𝒯) == project(∂¹, Taylor(2), Taylor(1), Float64)(a_𝒯) ==
            differentiate!(Sequence(Taylor(1), [Inf, Inf]), a_𝒯) == Sequence(Taylor(1), [-1.0, 2.0])
        @test ∂⁴(a_𝒯) == project(∂⁴, Taylor(2), Taylor(0), Float64)(a_𝒯) ==
            differentiate!(Sequence(Taylor(0), [Inf]), a_𝒯, 4) == Sequence(Taylor(0), [0.0])
        #
        a_ℱ = Sequence(Fourier(1, 1.0), [0.5, 0.0, 0.5])
        @test ∂¹(a_ℱ) == project(∂¹, Fourier(1, 1.0), Fourier(1, 1.0), ComplexF64)(a_ℱ) ==
            differentiate!(Sequence(Fourier(1, 1.0), ComplexF64[Inf, Inf, Inf]), a_ℱ) == Sequence(Fourier(1, 1.0), [-0.5im, 0.0, 0.5im])
        @test ∂⁴(a_ℱ) == project(∂⁴, Fourier(1, 1.0), Fourier(1, 1.0), ComplexF64)(a_ℱ) ==
            differentiate!(Sequence(Fourier(1, 1.0), ComplexF64[Inf, Inf, Inf]), a_ℱ, 4) == a_ℱ
        #
        a_𝒞 = Sequence(Chebyshev(2), [1.0, 0.5, 0.5])
        @test ∂¹(a_𝒞) == project(∂¹, Chebyshev(2), Chebyshev(1), Float64)(a_𝒞) ==
            differentiate!(Sequence(Chebyshev(1), [Inf, Inf]), a_𝒞) == Sequence(Chebyshev(1), [1.0, 2.0])
        #
        a_𝑇 = Sequence(Taylor(2) ⊗ Fourier(1, 1.0) ⊗ Chebyshev(2), collect(1.0:27.0))
        selectdim(a_𝑇, 2, 0) .= 0.0
        @test differentiate(integrate(a_𝑇, (1, 2, 1)), (1, 2, 1)) ≈ a_𝑇
        #
        a = Sequence(Taylor(2)^2 × Fourier(1, 1.0) × Chebyshev(2), collect(1.0:6.0+3.0+3.0))
        @test ∂¹(a) ==
            differentiate!(Sequence(Taylor(1)^2 × Fourier(1, 1.0) × Chebyshev(1), fill(complex(Inf), 2*2+3+2)), a)
    end

    @testset "Integral" begin
        ∫¹ = Integral(1)
        ∫² = Integral(2)
        @test ∫¹*∫¹ == ∫¹^2 == Integral(2)
        #
        a_𝒯 = Sequence(Taylor(2), [1.0, -1.0, 1.0])
        @test ∫¹(a_𝒯) == project(∫¹, Taylor(2), Taylor(3), Float64)(a_𝒯) ==
            integrate!(Sequence(Taylor(3), [Inf, Inf, Inf, Inf]), a_𝒯) == Sequence(Taylor(3), [0.0, 1.0, -1/2, 1/3])
        @test ∫²(a_𝒯) == project(∫², Taylor(2), Taylor(4), Float64)(a_𝒯) ==
            integrate!(Sequence(Taylor(4), [Inf, Inf, Inf, Inf, Inf]), a_𝒯, 2) == Sequence(Taylor(4), [0.0, 0.0, 1/2, -1/6, 1/12])
        #
        a_ℱ = Sequence(Fourier(1, 1.0), [0.5, 0.0, 0.5])
        @test ∫¹(a_ℱ) == project(∫¹, Fourier(1, 1.0), Fourier(1, 1.0), ComplexF64)(a_ℱ) ==
            integrate!(Sequence(Fourier(1, 1.0), ComplexF64[Inf, Inf, Inf]), a_ℱ) == Sequence(Fourier(1, 1.0), [0.5im, 0.0, -0.5im])
        @test ∫²(a_ℱ) == project(∫², Fourier(1, 1.0), Fourier(1, 1.0), ComplexF64)(a_ℱ) ==
            integrate!(Sequence(Fourier(1, 1.0), ComplexF64[Inf, Inf, Inf]), a_ℱ, 2) == -a_ℱ
        #
        a_𝒞 = Sequence(Chebyshev(2), [1.0, 0.5, 0.5])
        @test ∫¹(a_𝒞) ≈ project(∫¹, Chebyshev(2), Chebyshev(3), Float64)(a_𝒞) ≈
            integrate!(Sequence(Chebyshev(3), [Inf, Inf, Inf, Inf]), a_𝒞)
        #
        a = Sequence(Taylor(2)^2 × Fourier(1, 1.0) × Chebyshev(2), collect(1.0:6.0+3.0+3.0))
        component(a, 2)[0] = 0.0
        @test ∫¹(a) ==
            integrate!(Sequence(Taylor(3)^2 × Fourier(1, 1.0) × Chebyshev(3), fill(complex(Inf), 8+3+4)), a)
    end

    @testset "Scale" begin
        𝒮 = Scale(0.5)
        a  = Sequence(Taylor(3), [1.0, 2.0, 4.0, 8.0])
        @test 𝒮(a) == project(𝒮, Taylor(3), Taylor(3), Float64)(a) == Sequence(Taylor(3), [1.0, 1.0, 1.0, 1.0])
        #
        𝒮 = Scale((0.5, 0.5))
        a  = Sequence(Taylor(3) ⊗ Taylor(2), vec([2.0 ^ sum(α) for α ∈ indices(Taylor(3) ⊗ Taylor(2))]))
        @test 𝒮(a) == project(𝒮, Taylor(3) ⊗ Taylor(2), Taylor(3) ⊗ Taylor(2), Float64)(a) == Sequence(Taylor(3) ⊗ Taylor(2), ones(12))
    end

    @testset "Shift" begin
        𝒮 = Shift(π/2)
        a = Sequence(Fourier(1, 1.0), [0.5, 0.0, 0.5])
        @test 𝒮(a) ≈ Sequence(Fourier(1, 1.0), [-0.5im, 0.0, 0.5im])
    end
end


text:
using LinearAlgebra
using JuMP, Clp
include("ontoSimplex.jl")

function distFunction(states1::Vector{Float64}, states2::Vector{Float64})::Array{Float64,2}
    n1= length(states1); n2= length(states2)
    dMatrix= Array{Float64}(undef, (n1,n2))
    for i= 1:n1, j= 1:n2
		dMatrix[i,j]= abs(states2[j]- states1[i])
    end
    return dMatrix
end

#	computes the Wasserstein distance
function Wasserstein(p1::Vector{Float64}, p2::Vector{Float64}, distMatrix::Array{Float64,2}, rWasserstein::Float64=1.)
	ontoSimplex!(p1); ontoSimplex!(p2)
    n1= length(p1);   n2= length(p2)

    A= kron(ones(n2)', Matrix{Float64}(I, n1, n1))
    B= kron(Matrix{Float64}(I, n2, n2), ones(n1)')

	model = Model(Clp.Optimizer)
	@variable(model, x[i=1:n1*n2] >= 0)
	@objective(model, Min, vec(distMatrix.^rWasserstein)' * x)
	@constraint(model, [A;B] * x .== [p1;p2])
	optimize!(model)
    return (distance= (objective_value(model))^(1/ rWasserstein), π= reshape(value.(x), (n1, n2)))
end

#	Sinkhorn-Knopp iteration algorithm
function Sinkhorn(p1::Vector{Float64}, p2::Vector{Float64}, distMatrix::Array{Float64,2}, rWasserstein::Float64= 1., λ::Float64= 1.)
	ontoSimplex!(p1); ontoSimplex!(p2); count= 0
	βr= Array{Float64}(undef, length(p1));
	γc= ones(size(p2))		# guess a starting value
	distMatrix.^= rWasserstein; K= exp.(-λ * distMatrix)
	while count < 1000		# Sinkhorn iteration
		γc= γc./ (p2'*γc)		# rescale
		βr= p1./ (K* γc)		# vector operation
		γc= p2./ (K'* βr)		# vector operation
		count += 1				# iteration count
	end
#	println("r=", βr, p1'*βr); println("c=", γc, p2'*γc)
	π= Diagonal(βr)*K*Diagonal(γc)
	return (distance= (sum(π.* distMatrix)) ^(1/rWasserstein), π= π)
end

# #	Sinkhorn-Knopp algorithm
# function Sinkhorn(p1::Vector{Float64}, p2::Vector{Float64}, distMatrix::Array{Float64,2}, rWasserstein::Float64= 1., λ::Float64= 1.)
# 	ontoSimplex!(p1); ontoSimplex!(p2)
# 	π= exp.(-λ * (distMatrix.^ rWasserstein))
# 	for i= 1:100
# 		π.*= (p2'./ sum(π, dims=1)) # scale to sum of column= 1
# 		π.*= (p1 ./ sum(π, dims=2)) # scale to sum of lines= 1
# 	end
# 	return (distance= sum(π.* distMatrix)^(1/ rWasserstein), π= π)
# end


text:
push!(LOAD_PATH,"../src/")
using BAT
using Documenter

makedocs(
    modules = [BAT],
    clean = true,
    format = :html,
    sitename = "BAT",
    authors = "Oliver Schulz, Frederik Beaujean, and contributors",
    # linkcheck = !("skiplinks" in ARGS),
    pages = ["Home" => "index.md",
             "Manual" => Any[
                 "man/tutorial.md",
                 "man/basics.md",
             ],
             "auto.md",
             ],
    # Use clean URLs, unless built as a "local" build
    html_prettyurls = !("local" in ARGS),
    strict = true,
    checkdocs = :all,
)