text:
function project_path(path)
    if !isabspath(path)
        return joinpath(pwd(), path)
    else
        return path
    end
end

"""
create a project or package.

# Arguments

- `path`: path of the project you want to create

# Flags

- `-i, --interactive`: enable to start interactive configuration interface.
"""
@cast function create(path; interactive::Bool = false)
    fullpath = project_path(path)

    if ispath(fullpath)
        error("$path exists, remove it or use a new path")
    end

    if interactive
        t = Template(; dir = dirname(fullpath), interactive = true)
        t(basename(path))
    end

    # TODO: use .ionrc to save user configuration
    # and reuse it next time
    t = Template(; dir = dirname(fullpath))
    t(basename(path))
    return
end

function default_clone_name(url)
    name, _ = splitext(basename(url)) # rm .git
    _name, ext = splitext(name)
    if ext == ".jl" # preserve other extension
        name = _name
    end
    return name
end

"""
clone a package repo to a local directory.

# Arguments

- `url`: a remote or local url of the git repository.
- `to` : a local position, default to be the repository name (without .jl)

"""
@cast function clone(url, to = default_clone_name(url); credential = "")
    LibGit2.clone(url, to)
    return
end

"""
add package/project to the closest project.

# Arguments

- `url`: package name or url to add.

# Options

- `-v, --version <version number>`: package version, default is the latest available version, or master branch for git repos.
- `--rev <branch/commit>`: git revision, can be branch name or commit hash.
- `-s, --subdir <subdir>`: subdir of the package.

# Flags

- `-g, --glob`: add package to global shared environment.

"""
@cast function add(
    url;
    version::String = "",
    rev::String = "",
    subdir::String = "",
    glob::Bool = false,
)
    kwargs = []
    if isurl(url)
        push!(kwargs, "url=\"$url\"")
    else
        push!(kwargs, "name=\"$url\"")
    end

    !isempty(version) && push!(kwargs, "version=\"$version\"")
    !isempty(rev) && push!(kwargs, "rev=\"$rev\"")
    !isempty(subdir) && push!(kwargs, "subdir=\"$subdir\"")

    kw = join(kwargs, ", ")

    script = "using Pkg;"
    if !glob
        err_msg = "cannot install to global environment, use -g, --glob to install a package to global environment"
        script *= "(dirname(dirname(dirname(Pkg.project().path))) in DEPOT_PATH) && error(\"$err_msg\");"
    end

    script *= "Pkg.add(;$kw);"
    cmd = Cmd(["-e", script])

    if glob
        run(`$(Base.julia_cmd()) $cmd`)
    else
        withenv("JULIA_PROJECT" => "@.") do
            run(`$(Base.julia_cmd()) $cmd`)
        end
    end
end

"""
Make a package available for development. If pkg is an existing local path, that path will be recorded in the manifest and used.
Otherwise, a full git clone of pkg is made. Unlike the `dev/develop` command in Julia REPL pkg mode, `ion` will clone the package
to the dev folder of the current project. You can specify `--shared` flag to use shared `dev` folder under `~/.julia/dev`
(specified by `JULIA_PKG_DEVDIR`).

# Arguments

- `url`: URL or local path to the package.

# Flags

- `-s, --shared`: controls whether to use the shared develop folder.

"""
@cast function dev(url; shared::Bool = false)
    shared_flag = shared ? "--shared" : "--local"
    cmd = Cmd(["-e", "using Pkg; pkg\"dev $shared_flag $url\" "])
    withenv("JULIA_PROJECT" => "@.") do
        run(`$(Base.julia_cmd()) $cmd`)
    end
end

"""
Update a package. If no posistional argument is given, update all packages in current project.

# Arguments

- `pkg`: package name.
"""
@cast function update(pkg = "")
    if isempty(pkg)
        cmd = Cmd(["-e", "using Pkg; pkg\"up\" "])
    else
        cmd = Cmd(["-e", "using Pkg; pkg\"up $pkg\" "])
    end
    withenv("JULIA_PROJECT" => "@.") do
        run(`$(Base.julia_cmd()) $cmd`)
    end
end

@doc Docs.doc(update) @cast up(pkg = "") = update(pkg)

"""
"""
@cast function build(pkg = ""; verbose::Bool = false) end

# @cast function register()
# end


text:
using Impero, Test

@wrapper u=1 σ=1
∂x = Operator(nothing, DerivativeMetaData(nothing, "x"))
∂y = Operator(nothing, DerivativeMetaData(nothing, "y"))
∂z = Operator(nothing, DerivativeMetaData(nothing, "z"))
∂t = Operator(nothing, DerivativeMetaData(nothing, "t"))
@testset "Impero Equation Specification" begin
    equ = @to_equation σ=u
    @test equ.rhs == u
    @test equ.lhs == σ

    @pde_system pde_system = [
        σ = ∂x(u),
        ∂t(u) = -∂x(u * u - ∂x(σ)),
    ]
    @test pde_system[1].lhs == σ
    @test pde_system[2].lhs == ∂t(u)
    @test pde_system[1].rhs == ∂x(u)
    @test pde_system[2].rhs == -∂x(u * u - ∂x(σ))
end

@testset "Impero Equation Modification" begin
    @pde_system pde_system = [
        σ = ∂x(u),
        ∂t(u)= ∂x(∂x(σ)),
    ]
    pde_system[2].rhs += -∂x(u*u)
    rhs = (∂x(∂x(σ))+-(∂x((u*u))))
    @test pde_system[2].rhs == rhs
end

@testset "Impero Equation Augmentation" begin
    @pde_system pde_system = [
        σ = ∂x(u),
        ∂t(u)= ∂x(∂x(σ)),
    ]
    push!(pde_system, @to_equation σ=u)
    @test pde_system[3].lhs == σ
    @test pde_system[3].rhs == u
end

ext: jl


text:
module ImportanceSampler

# Samples p(x), to compute f(x), using the sampling distribution q(x).

using HDF5

# Writes samples to a file.

export ISampler, sampleIS, write_samples

"""Define the important functions for the sampler.
p must be a Probability density function, that takes in a single vector of parameter values.
f must be the function you want to evaluate, that takes in a single vector of parameter values.
q must be the sampling distribution, that produces a vector of parameter values. For best results, it should be tuned to match (p * f), but oversample the tails somewhat (e.g., using a student-t distribution when sampling a Gaussian).
"""
type ISampler
    nsamples::Int
    nparam::Int
    p::Function # lnprobability density
    f::Function # Function we wish to compute over this interval. For our inference problems,
    # this is typically a constant.
    q::Function # The proposal distribution, returns lnprobability
    xx::Matrix{Float64}
    p_x::Vector{Float64}
    f_x::Vector{Float64}
    q_x::Vector{Float64}
end

function ISampler(nsamples::Int, nparam::Int, p::Function, q::Function)
    # A dummy function that returns a constant
    f(x) = 1
    return ISampler(nsamples, nparam, p, f, q, Array(Float64, (nparam, nsamples)), Array(Float64, nsamples), Array(Float64, nsamples), Array(Float64, nsamples))
end

"""Load existing samples from an HDF5 file. Put dummy functions for functions since we won't be using them now."""
function ISampler(fname::AbstractString)
    fid = h5open(fname, "r")

    xx = read(fid["xx"])
    p_x = read(fid["p_x"])
    f_x = read(fid["f_x"])
    q_x = read(fid["q_x"])

    nparam, nsamples = size(xx)
    close(fid)

    f(x) = x

    return ISampler(nsamples, nparam, f, f, f, xx, p_x, f_x, q_x)

end

function sampleIS(sampler::ISampler)
    for i=1:sampler.nsamples
        # Draw a parameter combination from q.
        # Also report the likelihood at q(x)
        x, q_x = sampler.q()

        sampler.xx[:,i] = x
        sampler.q_x[i] = q_x

        # Evaluate p(x) and f(x)
        sampler.p_x[i] = sampler.p(x)
        sampler.f_x[i] = sampler.f(x)
    end

end

"""Write the values of the sampler to an HDF5 file."""
function write_samples(sampler::ISampler, fname::AbstractString)
    fid = h5open(fname, "w")
    fid["xx"] = sampler.xx
    fid["p_x"] = sampler.p_x
    fid["f_x"] = sampler.f_x
    fid["q_x"] = sampler.q_x
    close(fid)
end

# Functions to calculate moments.

end # module

max_stars_repo_name: elnjensen/DiskJockey

max_issues_repo_name: elnjensen/DiskJockey

max_forks_repo_name: elnjensen/DiskJockey


text:
function AngularAzElBuilder(arg0::CorrelatedRandomVectorGenerator, arg1::GroundStation, arg2::Vector{jdouble}, arg3::Vector{jdouble}, arg4::ObservableSatellite)
    return AngularAzElBuilder((CorrelatedRandomVectorGenerator, GroundStation, Vector{jdouble}, Vector{jdouble}, ObservableSatellite), arg0, arg1, arg2, arg3, arg4)
end

function build(obj::AngularAzElBuilder, arg0::Vector{SpacecraftState})
    return jcall(obj, "build", AngularAzEl, (Vector{SpacecraftState},), arg0)
end


text:
"""
This module contains various statistical models in the form of [`NestedModel`](@ref)s. These models can be used for examples and for testing.

* [`Models.GaussianShells`](@ref)
* [`Models.CorrelatedGaussian`](@ref)
"""
module Models

using ..NestedSamplers

using Distributions
using LinearAlgebra
using LogExpFunctions

include("shells.jl")
include("correlated.jl")
include("eggbox.jl")

end # module


text:
###############################################################################
###############################################################################
### Conversion between GAP and polymake objects
###############################################################################
###############################################################################

@doc Markdown.doc"""
    ntv_gap2polymake( v::GapObj )

Convert a `GAP` toric variety `v` into a `Polymake` toric variety.
"""
function ntv_gap2polymake(GapNTV::GapObj)
    ff = GAP.Globals.Fan(GapNTV)
    R = GAP.Globals.RayGenerators(ff)
    MC = GAP.Globals.RaysInMaximalCones(ff)
    rays = Matrix{Int}(R)
    cones = [findall(x->x!=0, vec) for vec in [[e for e in mc] for mc in MC]]
    Incidence = Polymake.IncidenceMatrix(cones)
    pmntv = Polymake.fulton.NormalToricVariety(
        RAYS = Oscar.matrix_for_polymake(rays),
        MAXIMAL_CONES = Incidence,
    )
    return pmntv
end
export ntv_gap2polymake

@doc Markdown.doc"""
    ntv_polymake2gap( v::Polymake.BigObject )

Convert a `Polymake` toric variety `v` into a `GAP` toric variety.
"""
function ntv_polymake2gap(polymakeNTV::Polymake.BigObject)
    rays = Matrix{Int}(Oscar.Polymake.common.primitive(polymakeNTV.RAYS))
    gap_rays = GapObj( rays, recursive = true )
    cones = [findall(x->x!=0, v) for v in eachrow(polymakeNTV.MAXIMAL_CONES)]
    gap_cones = GapObj( cones, recursive = true )
    fan = Oscar.GAP.Globals.Fan( gap_rays, gap_cones )
    variety = Oscar.GAP.Globals.ToricVariety( fan )
    return variety
end
export ntv_polymake2gap

function extract_gap_divisor_coeffs(GapDivisor::GapObj)
    return Vector{Int}(GAP.Globals.UnderlyingListOfRingElements(GAP.Globals.UnderlyingGroupElement(GapDivisor)))
end


text:
###include("Crystal.jl")
###include("DFToutMod.jl")
#using XMLDict

####################### QE specific
"""
    module QE

Module for running Quantum Espresso. Generic DFT version below.
"""
module QE
"""
Scripts to run Quantum Espresso
"""
#export load_xml
using LinearAlgebra

#using ExXML
using XMLDict
using GZip
using ..DFToutMod
using ..CrystalMod:crystal
using ..Utility:arr2str
using ..Utility:str_w_spaces
using ..Utility:parse_str_ARR_float
#using ..TB:get_grid
using ..CrystalMod:get_grid
using ..Atomdata:atoms

include("Commands.jl")

using ..ThreeBodyTB:TEMPLATES
using ..ThreeBodyTB:PSEUDOS

"""
    function run_pwscf(inputstr, outputstr, nprocs=1, directory="./", use_backup=false)

Run the pw.x code from QE on `inputstr`
"""
function run_pwscf(inputstr, outputstr, nprocs=1, directory="./", use_backup=false)
"""
run command
"""
    directory=rstrip(directory, '/')

    #get commandline
    c_dict = make_commands(nprocs)

    if use_backup
        qe = c_dict["qe_backup"]
    else
        qe = c_dict["qe"]
    end

    command = `$qe $directory/$inputstr`
    println("actual command")
    println(command)
    s=""
    try
        println("Running DFT")

        s = read(command, String)

        println("Writing output")
        f = open(directory*"/"*outputstr, "w")
        write(f, s)
        close(f)

        println("Ran DFT")
        return 0
    catch err
        println(err)
        println("Failed to run DFT")
        println(s)
        return -1
    end

end

"""
    function runSCF(crys::crystal, inputstr=missing, prefix=missing, tmpdir="./", directory="./", functional="PBESOL", wannier=0, nprocs=1, skip=false, calculation="scf", dofree="all", tot_charge = 0.0, smearing = 0.01, magnetic=false, cleanup=false, use_backup=false)

Workflow for doing SCF DFT calculation on `crys`

Return `dftout`
"""
function runSCF(crys::crystal, inputstr=missing, prefix="qe", tmpdir="./", directory="./", functional="PBESOL", wannier=0, nprocs=1, skip=false, calculation="scf", dofree="all", tot_charge = 0.0, smearing = 0.01, magnetic=false, cleanup=false, use_backup=false, grid=missing, klines=missing)
"""
Run SCF calculation using QE
"""
    qeout = []
#    println("runSCF 1")
    if skip == true
        try
            savedir="$tmpdir/$prefix.save"
            qeout = loadXML(savedir)
            println("skip: we loaded SCF data instead of rerunning")
            return qeout
        catch
            println("skipping failed, continue scf calculation")
        end
    end

    if ismissing(inputstr)
        if calculation == "nscf" || calculation == "nscf-sym"
            inputstr="qe.nscf.in"
        elseif calculation == "relax"
            inputstr="qe.relax.in"
        else
            inputstr="qe.in"
        end
    end
    outputstr = inputstr*".out"

    if !isdir(directory)
        println("warning, making directory path $directory")
        mkpath(directory)
    end

    savedir="$tmpdir/$prefix.save"

#    println("runSCF 2")

    tmpdir, prefix, inputfile =  makeSCF(crys, directory, prefix, tmpdir, functional, wannier, calculation, dofree, tot_charge, smearing, magnetic, grid=grid, klines=klines)

    f = open(directory*"/"*inputstr, "w")
    write(f, inputfile)
    close(f)

#    println("runSCF 3")

    ret = run_pwscf(inputstr, outputstr, nprocs, directory, use_backup)
#    println("DFT complete with error")
#    println("runSCF 4")

    updated = false
    if calculation == "relax" || calculation == "vc-relax"
        try
            qeout = loadXML(savedir)
            crys = qeout.crys
            println("updated crys")
            updated = true
        catch
            println("couldn't open results")
        end
    end

    if ret != 0 && updated == false
        println("failed DFT, trying with different mixing")
        tmpdir, prefix, inputfile =  makeSCF(crys, directory, prefix, tmpdir, functional, wannier, calculation, dofree, tot_charge, smearing, magnetic, mixing="TF", grid=grid, klines=klines)

        f = open(directory*"/"*inputstr, "w")
        write(f, inputfile)
        close(f)

        ret = run_pwscf(inputstr, outputstr, nprocs, directory, use_backup)

        if ret != 0
            println("warning, run_pwscf threw an error again: $ret")
            error("runSCF")
        end

    end

#    println("runSCF 5")

#    savedir="$tmpdir/$prefix.save"

    println("directory: ", directory)
    println("outputfile: ", outputstr)
    println("savedir: ", savedir)

    try
        qeout = loadXML(savedir)
    catch
        error("runSCF couldn't open output final try")
    end

    if cleanup
        println("call cleanup")
        doclean(savedir)
    end

#    println("runSCF 6")

    return qeout

end

"""
    function doclean(d)

Clean up wavefunctions in directory `d`
"""
function doclean(d)
    tot = 0

    files=read(`ls $d/`, String)
    for f in split(files, "\n")
        if occursin("wfc", f)
            if tot <= 2
                println("removing $d/$f")
            elseif tot==3
                println("removing $d/$f")
                println("etc...")
                println()
            end
            cmd=`\rm $d/$f`
            s = read(cmd, String)
            tot += 1
        end
    end
    println("clean-up done: $tot wfc files deleted")
end

"""
    function makeSCF(crys::crystal, directory="./", prefix=missing, tmpdir=missing, functional="PBESOL", wannier=0, calculation="scf", dofree="all", tot_charge = 0.0, smearing = 0.01, magnetic=false; mixing="local-TF")

Make QE inputfile for SCF DFT calculation.
"""
function makeSCF(crys::crystal, directory="./", prefix=missing, tmpdir=missing, functional="PBESOL", wannier=0, calculation="scf", dofree="all", tot_charge = 0.0, smearing = 0.01, magnetic=false; mixing="local-TF", grid=missing, klines = missing)
"""
Make inputfile for SCF calculation
"""

    c_dict = make_commands(1)

    template_file=open("$TEMPLATES/template_qe.in")
    temp = read(template_file, String)
    close(template_file)

    temp = replace(temp, "PSEUDODIR" => PSEUDOS)

    temp = replace(temp, "JULIANAT" => crys.nat)

    c = calculation
    if calculation == "nscf-sym"
        c = "nscf"
    end
    temp = replace(temp, "SCF" => c)

    temp = replace(temp, "MIXING" => mixing)

    if calculation == "nscf"
        temp = replace(temp, "nosym = false" => "nosym = true")
        temp = replace(temp, "noinv = false" => "noinv = true")
    end

    settypes = []
    for t in crys.types
        if !(t in settypes)
            push!(settypes, t)
        end
    end

#    println("settypes", settypes)

    ntypes = length(settypes)
    temp = replace(temp, "JULIANTYPE" => ntypes)

    brav="ibrav = 0"
    if abs(tot_charge) > 1e-5
        a = crys.A[1,1]
        brav="ibrav = 1, celldm(1) = $a "
    end
    temp = replace(temp, "JULIABRAV" => brav)

    t=""
    for i = 1:ntypes
        if functional=="PBESOL"
            psp=lowercase(settypes[i])*".pbesol.UPF"
        elseif functional=="LDA"
            psp=lowercase(settypes[i])*".lda.UPF"
        else
            psp=lowercase(settypes[i])*".pbe.UPF"
        end

        mass = atoms[settypes[i]].mass
        t = t * str_w_spaces([settypes[i], mass, psp, "\n"])
    end
    temp = replace(temp, "JULIAPSP\n" => t) ###strip(t,"\n"))

    if abs(tot_charge) < 1e-5
        temp = replace(temp, "JULIACELL\n" => arr2str(crys.A))#strip(arr2str(crys.A), "\n"))
    else
        temp = replace(temp, "JULIACELL\n" => "")#strip(arr2str(crys.A), "\n"))
        temp = replace(temp, "CELL_PARAMETERS\n"  => "")
    end

    t=""

    for i = 1:crys.nat
        t=t*crys.types[i]*" "
        t=t*arr2str(crys.coords[i,:])
        if i != crys.nat
            t=t*"\n"
        end

    end
    temp = replace(temp, "JULIACOORDS" => t)

    temp = replace(temp, "JULIASMEARING" => smearing)

#=    B = transpose(inv(crys.A))
    kden = 55.0
    b1 = norm(B[1,:])
    b2 = norm(B[2,:])
    b3 = norm(B[3,:])

    k1 = convert(Int, round(kden * b1))
    k2 = convert(Int, round(kden * b2))
    k3 = convert(Int, round(kden * b3))

    if k1%2 == 1
        k1 = k1 + 1
    end
    if k2%2 == 1
        k2 = k2 + 1
    end
    if k3%2 == 1
        k3 = k3 + 1
    end

    kpoints = [k1, k2, k3]
=#

    if ismissing(klines)
        if ismissing(grid)
            kpoints = get_grid(crys)
        else
            if length(grid) == 3
                kpoints =  grid
            else
                kpoints = get_grid(crys, grid)
            end
        end

#    if calculation == "nscf"
        #        kpoints = get_grid(c, kden=50.0)
        #    end
        temp = replace(temp, "JULIAKTYPE" => "K_POINTS automatic")
        temp = replace(temp, "JULIAKPOINTS" => arr2str(kpoints)*" 0 0 0 ")

    else

#        println("KLINES xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx")
#        println(klines)

        if calculation != "nscf"
            println("warning, you are specifying the klines without using nscf, i hope you know what you are doing")
        end

        nk = size(klines)[1]
        if size(klines)[2] == 3
            klines2 = ones(nk,4)
            klines2[:,1:3] = klines
            klines = deepcopy(klines2)
        end

        kstr = "$nk\n"
        kstr *= arr2str(klines)

#        println("before ", temp)

        temp = replace(temp, "JULIAKTYPE" => "K_POINTS crystal")
        temp = replace(temp, "JULIAKPOINTS" => kstr)

#        println("after", temp)

    end

    if ismissing(tmpdir)
        tmpdir=string(rand(1)[1])
    end

    if ismissing(prefix)
        prefix=string(rand(1)[1])
    end

    temp = replace(temp, "PREFIX" => prefix)
    temp = replace(temp, "TMPDIR" => tmpdir)

    other=""
    if (typeof(wannier) == Bool && wannier==true) || wannier == 2
        nbandsemi = 0
        nbandval = 0

        for t in crys.types
            nbandsemi += atoms[t].nsemicore
            nbandval += atoms[t].nwan
        end
        nbandtot = 4+1+convert(Int64, round(nbandsemi/2.0 + nbandval * 3.0 / 2.0))   #no spin yet
        other *="nbnd = "*string(nbandtot)*"\n"
    elseif wannier == 1 #fewer bands can be more stable
        nbandsemi = 0
        nbandval = 0

        for t in crys.types
            nbandsemi += atoms[t].nsemicore
            nbandval += atoms[t].nwan
        end
        nbandtot = 0+convert(Int64, round(nbandsemi/2.0 + nbandval * 2.0 / 2.0))   #no spin yet
        other *="nbnd = "*string(nbandtot)*"\n"
    elseif wannier == -1
        nbandsemi = 0
        nbandval = 0

        for t in crys.types
            nbandsemi += atoms[t].nsemicore
            nbandval += atoms[t].nwan
        end
        nbandtot = 0+convert(Int64, round(nbandsemi/2.0 + nbandval * 1.5 / 2.0))   #no spin yet
        other *="nbnd = "*string(nbandtot)*"\n"
    end

    if abs(tot_charge) > 1e-5
        other *= "   tot_charge = $tot_charge , assume_isolated = 'mp' \n"
    end

    if magnetic
        other *= "   nspin = 2 \n"
        for c in 1:ntypes
            other *= "  starting_magnetization( $c ) = 0.7 \n"
        end
    end

    temp = replace(temp, "JULIAOTHER" => other)

    st_free = "cell_dofree = 'all'"
    st_free2 = replace(st_free, "all"=> dofree)

    temp = replace(temp, "JULIADOFREE" => st_free2)

    if calculation == "vc-relax"
        temp = replace(temp, "JULIACELLFACTOR" => "cell_factor = 3.0")
    else
        temp = replace(temp, "JULIACELLFACTOR" => "cell_factor = 1.0")
    end

    return tmpdir, prefix, temp

end

"""
   function makedict(savedir)

Load a data-file-schema.xml as julia dictionary.
"""
function makedict(savedir)

    if isfile(savedir*"/data-file-schema.xml")
        filename=savedir*"/data-file-schema.xml"
    elseif isfile(savedir*"/data-file-schema.xml.gz")
        filename=savedir*"/data-file-schema.xml.gz"
    else
        println("error warning missing data-file-schema.xml or data-file-schema.xml.gz")
        filename=missing
    end

    f = gzopen(filename, "r")
    fs = read(f, String)
    close(f)

    d = xml_dict(fs)

    return d
end

"""
    function loadXML(savedir)

Load a QE SCF DFT calculation into a `dftout` object.
"""
function loadXML(savedir)

#    println("start loadXML")
    convert_ha_ryd = 2.0

    d= makedict(savedir)

    stin = d["espresso"]["output"]["atomic_structure"]
    stout = d["espresso"]["output"]["atomic_structure"]
    a1 = stout["cell"]["a1"]
    a2 = stout["cell"]["a2"]
    a3 = stout["cell"]["a3"]

#    println("a1 ", a1)
#    println("a2 ", a2)
#    println("a3 ", a3)

    A = zeros((3,3))
    A[1,:] = parse_str_ARR_float(a1)
    A[2,:] = parse_str_ARR_float(a2)
    A[3,:] = parse_str_ARR_float(a3)

#    println(A)

    nat = parse(Int, stin[:nat])
    types = []
    pos = zeros(nat,3)

#    println(d["espresso"]["input"]["k_points_IBZ"])

    if "monkhorst_pack" in keys(d["espresso"]["input"]["k_points_IBZ"])
        nk1 = parse(Int, d["espresso"]["input"]["k_points_IBZ"]["monkhorst_pack"][:nk1])
        nk2 = parse(Int, d["espresso"]["input"]["k_points_IBZ"]["monkhorst_pack"][:nk2])
        nk3 = parse(Int, d["espresso"]["input"]["k_points_IBZ"]["monkhorst_pack"][:nk3])
    else
        nk1 = 0
        nk2 = 0
        nk3 = 0
    end

    prefix = d["espresso"]["input"]["control_variables"]["prefix"]
    outdir = d["espresso"]["input"]["control_variables"]["outdir"]

    tot_charge = parse(Float64, d["espresso"]["input"]["bands"]["tot_charge"])

    kgrid = [nk1, nk2,nk3]
#    println("nat ", nat)
#    println(st["atomic_positions"]["atom"])

    if nat == 1
        atom =stout["atomic_positions"]["atom"]
        push!(types, atom[:name])
        pos[1,:] = parse_str_ARR_float(atom[""])

    else
        for at = 1:nat
#            println(at)
            atom =stout["atomic_positions"]["atom"][at]
            push!(types, atom[:name])
            pos[at,:] = parse_str_ARR_float(atom[""])
        end
    end

    coords_crys = pos * inv(A)

    if "forces" in keys(d["espresso"]["output"])
        out = d["espresso"]["output"]["forces"]
        f_lin = parse_str_ARR_float(out[""])
        forces = transpose(reshape(f_lin , 3,nat)) * convert_ha_ryd
    else
        forces = -99.0*ones(nat, 3)
    end

    if "stress" in keys(d["espresso"]["output"])

        out = d["espresso"]["output"]["stress"]
        s_lin = parse_str_ARR_float(out[""])
        stress = transpose(reshape(s_lin , 3,3)) * convert_ha_ryd
    else
        stress = -99.0*ones(3,3)
    end

    out = d["espresso"]["output"]["total_energy"]
    energy = parse(Float64,out["etot"]) * convert_ha_ryd
    energy_smear = parse(Float64,out["demet"]) * convert_ha_ryd

    bs = loadXML_bs(d)

    bs.kgrid=kgrid

    d = DFToutMod.makedftout(A, coords_crys, types, energy, energy_smear, forces, stress, bs, prefix=prefix, outdir=outdir, tot_charge=tot_charge)

#    println("end loadXML")

    return d

end

"""
    function loadXML_bs(savedir::String)

Load `bandstructure` from QE xml file.
"""
function loadXML_bs(savedir::String)

    d=makedict(savedir)

    return loadXML_bs(d)

end

"""
    function loadXML_bs(savedir::String)

Load `bandstructure` from QE xml file that was converted to a dict already
"""
function loadXML_bs(d)

    convert_ha_ryd = 2.0

    out = d["espresso"]["output"]["band_structure"]
    nbnd = parse(Int,out["nbnd"])
    nelec = parse(Float64,out["nelec"])
    efermi = parse(Float64,out["fermi_energy"]) * convert_ha_ryd

    st = d["espresso"]["output"]

    b1 = st["basis_set"]["reciprocal_lattice"]["b1"]
    b2 = st["basis_set"]["reciprocal_lattice"]["b2"]
    b3 = st["basis_set"]["reciprocal_lattice"]["b3"]

    B = zeros((3,3))
    B[1,:] = parse_str_ARR_float(b1)
    B[2,:] = parse_str_ARR_float(b2)
    B[3,:] = parse_str_ARR_float(b3)

    nks = parse(Int, out["nks"])
    bandstruct = zeros(nks, nbnd)
    kpts = zeros(nks, 3)
    weights = zeros(nks)
    for b = 1:nks
        ks = out["ks_energies"][b]

        weights[b] = parse(Float64, ks["k_point"][:weight] )
        kpts[b,:] = parse_str_ARR_float(ks["k_point"][""])
        bandstruct[b,:] = parse_str_ARR_float(ks["eigenvalues"][""])
    end

    kpts = kpts * inv(B)

    bandstruct *= convert_ha_ryd

    bs = DFToutMod.makebs(nelec,efermi, kpts, weights, [0,0,0], bandstruct)

    return bs

end

end #module

############################## General DFT functions
"""
    module DFT

This is the generic DFT interface. Only QE is currently implemented however.
"""
module DFT
"""
Scripts to run DFT codes
"""
#using ..CrystalMod
using ..QE
using ..CrystalMod:crystal

#include("Atomdata.jl")
#using ..Atomdata:atoms

"""
    function runSCF(crys::crystal; inputstr=missing, prefix=missing, tmpdir="./", directory="./", functional="PBESOL", wannier=0, nprocs=1, code="QE", skip=false, calculation="scf", dofree="all", tot_charge = 0.0, smearing = missing, magnetic=false, cleanup=false, use_backup=false)

Workflow for generic DFT SCF calculation. `code` can only by "QE"
"""
function runSCF(crys::crystal; inputstr=missing, prefix=missing, tmpdir="./", directory="./", functional="PBESOL", wannier=0, nprocs=1, code="QE", skip=false, calculation="scf", dofree="all", tot_charge = 0.0, smearing = missing, magnetic=false, cleanup=false, use_backup=false, grid=missing, klines=missing)

    if ismissing(smearing)
        if calculation=="scf"
            smearing = 0.01
        else
            smearing = 0.02
        end
    end

    println(functional)
    if code == "QE"
        if ismissing(prefix)
            prefix="qe"
        end
        qeout = missing
        try
            qeout = QE.runSCF(crys, inputstr, prefix, tmpdir, directory, functional, wannier, nprocs, skip, calculation, dofree, tot_charge, smearing, magnetic, cleanup, use_backup, grid, klines)
            return qeout
        catch
            println("WARNING failure, restart qe with higher smearing, hope that helps!!")
            qeout = QE.runSCF(crys, inputstr, prefix, tmpdir, directory, functional, wannier, nprocs, skip, calculation, dofree, tot_charge, smearing*5, magnetic, cleanup, true, grid, klines)
            return qeout

        end

    else
        println("code variable not recognized", code)
        println("currently supported: QE")
    end

end

end #end RunDFT


text:
using Test, Batched, LinearAlgebra

@testset "Testing batched_gemm! with $elty" for elty in [Float32, Float64, ComplexF32, ComplexF64]

    alpha, beta = rand(elty), rand(elty)
    A = rand(elty, 10, 10, 100)
    B = rand(elty, 10, 10, 100)
    C = rand(elty, 10, 10, 100)
    test_C = copy(C)

    Batched.batched_gemm!('N', 'N', alpha, A, B, beta, C)

    for k in 1:100
        test_C[:, :, k] = alpha * (A[:, :, k]) * B[:, :, k] + beta * test_C[:, :, k]
    end

    @test C ≈ test_C
end

@testset "Testing batched_gemm with $elty" for elty in [Float32, Float64, ComplexF32, ComplexF64]
    alpha = rand(elty)
    A = rand(elty, 10, 10, 100)
    B = rand(elty, 10, 10, 100)

    C = Batched.batched_gemm('N', 'N', alpha, A, B)

    test_C = zeros(elty, 10, 10, 100)
    for k in 1:100
        test_C[:, :, k] = alpha * (A[:, :, k]) * B[:, :, k]
    end

    @test test_C ≈ C
end

@testset "Testing batched_gemm (uniform scaling) with $elty" for elty in [Float32, Float64, ComplexF32, ComplexF64]
    A = BatchedUniformScaling(rand(elty, 100))
    B = rand(elty, 10, 10, 100)

    C = Batched.batched_gemm(A, B)

    test_C = zeros(elty, 10, 10, 100)
    for k in 1:100
        test_C[:, :, k] = (A[k]) * B[:, :, k]
    end

    @test test_C ≈ C
end

@testset "Testing batched_gemm (uniform scaling & transpose) with $elty" for elty in [Float32, Float64, ComplexF32, ComplexF64]
    A = BatchedUniformScaling(rand(elty, 100))
    vB = rand(elty, 10, 8, 100)
    B = batched_transpose(vB)

    C = Batched.batched_gemm(A, B)

    test_C = zeros(elty, 8, 10, 100)
    for k in 1:100
        test_C[:, :, k] = (A[k]) * transpose(vB[:, :, k])
    end

    @test test_C ≈ C
end


text:
module UIControls

using ..Twinkle
import ..UIVariables
using Parameters
import UUIDs

#------------------------------------------------------------------------------
# EXports
#------------------------------------------------------------------------------
export  Container,
        Slider,
        Button,
        MeshCatViewer,
        Label,
        Image,
        PanZoom,
        Field,
        ButtonToggle,
        RadioGroup,
        CheckBox,
        ExpansionPanel,
        Accordion,
        Tabs,
        Tab,
        Divider,
        Card,
        Markdown,
        MDJulia,

        CodeSnip,
        CodeSnipJulia,

        VContainer,
        HContainer,
        HContainerSpace,
        H1Label,
        H2Label,
        H3Label,
        H4Label,

        DummyExport

# mutable struct Slider <: UIControl
# end

# typedict(x) = Dict(fn=>getfield(x, fn) for fn ∈ fieldnames(typeof(x)))
typedict(x) = Dict(fn=>getfield(x, fn) for fn ∈ filter(x->!startswith(string(x), "_"), fieldnames(typeof(x))))

# function on(func::Function, control::AbstractUIControl)
#     control._func = func
# end

# function asInt(control::AbstractUIControl)::Int
#     try
#         val = value(control)
#         @show val, typeof(val)
#         if (isa(val, String))
#             return Int(floor(parse(Float64, val)))
#         end
#         if (isa(val, Int64))
#             return val
#         end
#         if (isa(val, Float64))
#             return Int(floor(val))
#         end
#     catch
#         return 0
#     end
# end

# function asFloat(control::AbstractUIControl)::Int
#     try
#         val = value(control)
#         if (isa(val, String))
#             return parse(Float64, val)
#         end
#         if (isa(val, Int))
#             return convert(Float64, val)
#         end
#         if (isa(val, Float64))
#             return val
#         end
#         @error "should not reach here"
#     catch
#         return 0
#     end
# end

@with_kw mutable struct Container <: AbstractUIControl
    type::String = "container"
    direction::String = "row"
    gap::String = "10px"
    align::String = "left left"
    children::Vector{UIControls.AbstractUIControl} = []

    _app::Union{Nothing, Twinkle.AbstractUIApp} = nothing
end

VContainer(args...) = Container(direction = "column", children=[args...])
HContainer(args...) = Container(direction = "row", align="left center",  children=[args...])
HContainerSpace(args...) = Container(direction = "row", align="space-between center",  children=[args...])

@with_kw mutable struct Slider <: AbstractUIControl
    type::String = "slider"
    text::String = ""
    trailing_text::String = ""
    min::Float64 = 0.0
    max::Float64 = 100.0
    value::Float64 = 0
    step::Float64 = 1.0
    variable::Any = nothing

    _app::Union{Nothing, Twinkle.AbstractUIApp} = nothing
end
# value(c::Slider) = c.value

@with_kw mutable struct Button <: AbstractUIControl
    type::String = "button"
    text::String = "default button"
    variable::Any = nothing
    buttonType::String = "raised"     # normal, raised, stroked, flat, icon, fab, mini-fab
    color::String = "primary"         # normal, primary, accent, warn,
    fileTypes::String = ""            # used for an open file dialog, like ".png,.jpg"

    _variable_name_to_create::String = ""
    _app::Union{Nothing, Twinkle.AbstractUIApp} = nothing
end
# value(c::Button) = c.text

function Button(text::String, var_name::String)
    return Button(text=text, _variable_name_to_create=var_name)
end

@with_kw mutable struct MeshCatViewer <: AbstractUIControl
    type::String = "meshcat-viewer"
    url::String = ""
    width::Any = "100%"
    height::Any = "600px"

    _app::Union{Nothing, Twinkle.AbstractUIApp} = nothing
end
# value(c::Button) = c.text

@with_kw mutable struct Label <: AbstractUIControl
    type::String = "label"
    class::String = "normal"
    text::Any = ""
    style::Any = ""
    variable::Any = nothing
    isHTML::Bool = false

    _app::Union{Nothing, Twinkle.AbstractUIApp} = nothing
end

Label(text::String) = Label(text=text)
H1Label(text::String) = Label(text=text, class="h1")
H2Label(text::String) = Label(text=text, class="h2")
H3Label(text::String) = Label(text=text, class="h3")
H4Label(text::String) = Label(text=text, class="h4")

@with_kw mutable struct Image <: AbstractUIControl
    type::String = "image"
    style::Any = ""
    source::Any = nothing

    _app::Union{Nothing, Twinkle.AbstractUIApp} = nothing
end

@with_kw mutable struct PanZoom <: AbstractUIControl
    type::String = "pan-zoom"
    style::Any = "width: 100%; height: 300px;"
    innerStyle::Any = "width: auto; height: calc(100% - 30px);"
    content::Union{Nothing, UIControls.Container} = nothing

    _app::Union{Nothing, Twinkle.AbstractUIApp} = nothing
end

@with_kw mutable struct Field <: AbstractUIControl
    type::String = "field"
    input::String = "number"
    label::String = "default label"
    hint::String = "hint"
    options::Union{Nothing, Vector{Any}} = nothing
    variable::Any = nothing

    _func::Union{Nothing, Function} = nothing
end

@with_kw mutable struct ButtonToggle <: AbstractUIControl
    type::String = "button-toggle"
    options::Union{Nothing, Vector{Any}} = nothing
    variable::Any = nothing

    _func::Union{Nothing, Function} = nothing
end

@with_kw mutable struct RadioGroup <: AbstractUIControl
    type::String = "radio-group"
    direction::String = "horizontal"
    options::Union{Nothing, Vector{Any}} = nothing
    variable::Any = nothing

    _func::Union{Nothing, Function} = nothing
end

@with_kw mutable struct CheckBox <: AbstractUIControl
    type::String = "checkbox"
    label::String = "default label"
    variable::Any = nothing

    _func::Union{Nothing, Function} = nothing
end

@with_kw mutable struct ExpansionPanel <: AbstractUIControl
    type::String = "expansion-panel"
    title::String = "defaut title"
    subtitle::String = ""
    content::Union{Nothing, UIControls.Container} = nothing
    style::Any = ""
    titleStyle::String = ""
    subtitleStyle::String = ""
    headerStyle::Any = ""

    _app::Union{Nothing, Twinkle.AbstractUIApp} = nothing
end

@with_kw mutable struct Accordion <: AbstractUIControl
    type::String = "accordion"
    panels::Union{Nothing, Vector{UIControls.ExpansionPanel}} = nothing

    _app::Union{Nothing, Twinkle.AbstractUIApp} = nothing
end

@with_kw mutable struct Tab <: AbstractUIControl
    type::String = "tab"
    label::String = "default label"
    content::Union{Nothing, UIControls.Container} = nothing

    _app::Union{Nothing, Twinkle.AbstractUIApp} = nothing
end

@with_kw mutable struct Tabs <: AbstractUIControl
    type::String = "tabs"
    style::Any = ""
    tabs::Union{Nothing, Vector{UIControls.Tab}} = nothing

    _app::Union{Nothing, Twinkle.AbstractUIApp} = nothing
end

@with_kw mutable struct Divider <: AbstractUIControl
    type::String = "divider"

    _app::Union{Nothing, Twinkle.AbstractUIApp} = nothing
end

@with_kw mutable struct Card <: AbstractUIControl
    type::String = "card"
    style::Any = ""
    title::String = ""
    titleVariable::Any = ""
    titleStyle::String = ""
    subtitle::String = ""
    subtitleVariable::Any = ""
    subtitleStyle::String = ""
    content::Union{String, UIControls.Container} = ""
    contentStyle::Any = ""
    variable::Any = ""                  # content source
    isHTML::Bool = false

    _app::Union{Nothing, Twinkle.AbstractUIApp} = nothing
end

@with_kw mutable struct Markdown <: AbstractUIControl
    type::String = "markdown"
    style::Any = ""
    content::String = ""
    variable::Any = ""                  # content source
    lineNumbers::Bool = false

    _app::Union{Nothing, Twinkle.AbstractUIApp} = nothing
end

MDJulia(code) = Markdown(content = """```julia\n$(code)\n```""")

@with_kw mutable struct CodeSnip <: AbstractUIControl
    type::String = "codesnip"
    text::String = ""
    language::String = "julia"
    lineNumbers::Bool = false
    commandLine::Bool = false
    prompt::String = "julia>"
    promptStyle::String = "-webkit-text-fill-color: rgb(51, 200, 51); -webkit-text-stroke-width:0.2px; -webkit-text-stroke-color: black;"
    rawLines::String = "2-999"

    _app::Union{Nothing, Twinkle.AbstractUIApp} = nothing
end

CodeSnipJulia(code) = CodeSnip(text=code, commandLine=true, prompt="julia>")

end # module UIControlsget


text:
mutable struct SMTSolver{name} <: NLSolver
    vars::Dict{Symbol,Type}
    cs::ClauseSet

    function SMTSolver{name}() where {name}
        @assert exists(SMTSolver{name})
        @assert name isa Symbol
        new(Dict(), ClauseSet())
    end
end

const Z3Solver = SMTSolver{:z3}
const YicesSolver = SMTSolver{:yices}

SMTSolver(s::String) = SMTSolver{Symbol(s)}

program_name(::Type{SMTSolver{name}}) where {name} = string(name)
program_name(::T) where {T<:SMTSolver} = program_name(T)
program_name(::YicesSolver) = "$(program_name(YicesSolver)) --logic=QF_NRA"

exists(::Type{T}) where {T<:SMTSolver} = !isnothing(Sys.which(program_name(T)))

fileext(::SMTSolver) = ".smt2"
parsefunc(::SMTSolver) = parse_output_smt
parsefunc(::YicesSolver) = parse_output_yices

function variables!(s::SMTSolver, d::Dict{Symbol,Type})
    push!(s.vars, d...)
end

function write_header(io::IO, s::SMTSolver)
    write(io, "(set-option:produce-models true)\n")
    write(io, "(set-logic QF_NRA)\n")
end

function write_header(io::IO, s::Z3Solver) end
function write_header(io::IO, s::YicesSolver) end

function write_vars(io::IO, s::SMTSolver)
    for (v,t) in s.vars
        write(io, "(declare-const $v Real)\n")
    end
end

function write_vars(io::IO, s::YicesSolver)
    for (v,t) in s.vars
        write(io, "(define $v::real)\n")
    end
end

function write_constraints(io::IO, s::SMTSolver)
    write(io, smt(s.cs; expand_pow=true), "\n")
end

function write_constraints(io::IO, s::Union{Z3Solver,YicesSolver})
    write(io, smt(s.cs; expand_pow=false), "\n")
end

function write_footer(io::IO, s::SMTSolver)
    write(io, "(check-sat)\n")
    write(io, "(get-value ($(join(keys(s.vars), " "))))\n")
end

function write_footer(io::IO, s::YicesSolver)
    write(io, "(check)\n")
    write(io, "(show-model)\n")
end

function write_problem(io::IO, s::T) where {T<:NLSolver}
    write_header(io, s)
    write_vars(io, s)
    write_constraints(io, s)
    write_footer(io, s)
end

function solve(s::SMTSolver; timeout::Int=-1)
    path, io = mktemp()
    write_problem(io, s)
    close(io)
    newpath = string(path, fileext(s))
    mv(path, newpath)
    try
        openproc(parsefunc(s), `$(split(program_name(s))) $newpath`, timeout=timeout)
    finally
        rm(newpath)
    end
end

# ------------------------------------------------------------------------------

function openproc(parse::Function, cmd::Cmd; timeout=-1)
    start = time_ns()
    P = open(pipeline(cmd, stderr=devnull))
    if timeout < 0
        wait(P)
    else
        timedwait(()->!process_running(P), float(timeout))
        if process_running(P)
            @debug "Kill"
            kill(P)
            close(P.in)
            return NLSat.timeout, Second(timeout), nothing
        end
    end
    elapsed = Millisecond(round((time_ns()-start)/1e6))

    if P.exitcode >= 0
        lines = readlines(P)
        status = popfirst!(lines)
        if status == "sat"
            d = parse(lines)
            return NLSat.sat, elapsed, d
        elseif status == "unsat"
            return NLSat.unsat, elapsed, nothing
        elseif status == "unknown"
            return NLSat.unknown, elapsed, nothing
        end

        error("Unknown status: $status")
    end
    return NLSat.unknown, elapsed, nothing
end

function parse_output_smt(s::Array{String})
    m = NLModel()
    sexpr = parse_sexpr(join(s))
    for x in sexpr.vec
        @assert length(x.vec) == 2
        sym = x.vec[1]
        val = try_int(eval(to_expr(x.vec[2])))
        push!(m, sym=>val)
    end
    m
end

function parse_output_yices(s::Array{String})
    m = NLModel()
    sexpr = parse_sexpr(string("(", join(s), ")"))
    for x in sexpr.vec
        @assert length(x.vec) == 3 && x.vec[1] == :(=) "$x"
        sym = x.vec[2]
        val = try_int(eval(to_expr(x.vec[3])))
        push!(m, sym=>val)
    end
    m
end

try_int(x) = x === NaN ? x : isinteger(x) ? convert(Int, x) : rationalize(x)

Base.rationalize(x::Rational) = x


text:
################################
# generative function wrappers #
################################

# generative function that takes a factor graph and the result of variable
# elimination and samples from the conditional distribution on the latent
# variables

struct SampleFactorGraphTrace <: Gen.Trace
    fg::FactorGraph
    args::Tuple{FactorGraph,VariableEliminationResult}
    choices::Gen.DynamicChoiceMap
    log_prob::Float64
end

struct SampleFactorGraph <: GenerativeFunction{Nothing,SampleFactorGraphTrace}
end

const sample_factor_graph = SampleFactorGraph()

function Gen.simulate(
        ::SampleFactorGraph,
        args::Tuple{FactorGraph{N},VariableEliminationResult}) where {N}
    (fg, elimination_result) = args
    choices = choicemap()
    values = Vector{Any}(undef, N)
    log_prob = 0.0
    for addr in reverse(elimination_result.elimination_order)
        var_idx = addr_to_idx(fg, addr)
        fg = elimination_result.intermediate_fgs[var_idx]
        dist = conditional_dist(fg, values, addr)
        idx = categorical(dist)
        value = idx_to_value(idx_to_var_node(fg, addr_to_idx(fg, addr)), idx)
        values[var_idx] = value
        choices[addr] = value
        log_prob += log(dist[idx])
    end
    return SampleFactorGraphTrace(fg, args, choices, log_prob)
end

function Gen.generate(
        ::SampleFactorGraph,
        args::Tuple{FactorGraph{N},VariableEliminationResult},
        choices::ChoiceMap) where {N}
    (fg, elimination_result) = args
    values = Vector{Any}(undef, N)
    for addr in elimination_result.elimination_order
        values[addr_to_idx(fg, addr)] = choices[addr]
    end
    log_prob = 0.0
    for addr in reverse(elimination_result.elimination_order)
        fg = elimination_result.intermediate_fgs[addr_to_idx(fg, addr)]
        dist = conditional_dist(fg, values, addr)
        idx = value_to_idx(idx_to_var_node(fg, addr_to_idx(fg, addr)), values[addr_to_idx(fg, addr)])
        log_prob += log(dist[idx])
    end
    trace = SampleFactorGraphTrace(fg, args, DynamicChoiceMap(choices), log_prob)
    return (trace, log_prob)
end

Gen.get_args(trace::SampleFactorGraphTrace) = trace.args
Gen.get_retval(trace::SampleFactorGraphTrace) = nothing
Gen.get_choices(trace::SampleFactorGraphTrace) = trace.choices
Gen.get_score(trace::SampleFactorGraphTrace) = trace.log_prob
Gen.get_gen_fn(trace::SampleFactorGraphTrace) = sample_factor_graph
Gen.project(trace::SampleFactorGraphTrace, ::EmptySelection) = 0.0
Gen.has_argument_grads(::SampleFactorGraph) = (false, false, false)
Gen.accepts_output_grad(::SampleFactorGraph) = false

# TODO: to use this as an MH proposal, it must take its first argument to be a trace
# ideally, we could wrap it in another generatiev function like this:

#@gen function mh_proposal(trace, fg, elimination_result)
    #{*} ~ sample_factor_graph(fg, elimination_result)
#end

# however, this is not supported...
# we could embed it in involutive MH, but unfortunately that requires some code
# in principle, it should be possible to do this without requiring any user code.

# generative function that takes a trace and compiles and samples from the
# conditional distribution on the latent variables

struct CompileAndSampleFactorGraphTrace <: Gen.Trace
    fg::FactorGraph
    args::Tuple{Gen.Trace,Dict{Any,Latent},Dict{Any,Observation},Any}
    choices::Gen.DynamicChoiceMap
    log_prob::Float64
end

struct CompileAndSampleFactorGraph <: GenerativeFunction{Nothing,CompileAndSampleFactorGraphTrace}
end

const compile_and_sample_factor_graph = CompileAndSampleFactorGraph()

function Gen.simulate(
        ::CompileAndSampleFactorGraph,
        args::Tuple{Gen.Trace,Dict{Any,Latent},Dict{Any,Observation},Any})
    (trace, latents, observations, elimination_order) = args
    fg = compile_trace_to_factor_graph(trace, latents, observations)
    elimination_result = variable_elimination(fg, elimination_order)
    trace = Gen.simulate(sample_factor_graph, (fg, elimination_result))
    return CompileAndSampleFactorGraphTrace(fg, args, get_choices(trace), get_score(trace))
end

function Gen.generate(
        ::CompileAndSampleFactorGraph,
        args::Tuple{Gen.Trace,Dict{Any,Latent},Dict{Any,Observation},Any},
        choices::ChoiceMap)
    (trace, latents, observations, elimination_order) = args
    fg = compile_trace_to_factor_graph(trace, latents, observations)
    elimination_result = variable_elimination(fg, elimination_order)
    trace, weight = Gen.generate(sample_factor_graph, (fg, elimination_result), choices)
    return (CompileAndSampleFactorGraphTrace(fg, args, get_choices(trace), get_score(trace)), weight)
end

Gen.get_args(trace::CompileAndSampleFactorGraphTrace) = trace.args
Gen.get_retval(trace::CompileAndSampleFactorGraphTrace) = nothing
Gen.get_choices(trace::CompileAndSampleFactorGraphTrace) = trace.choices
Gen.get_score(trace::CompileAndSampleFactorGraphTrace) = trace.log_prob
Gen.get_gen_fn(trace::CompileAndSampleFactorGraphTrace) = sample_factor_graph
Gen.project(trace::CompileAndSampleFactorGraphTrace, ::EmptySelection) = 0.0
Gen.has_argument_grads(::CompileAndSampleFactorGraph) = (false, false, false)
Gen.accepts_output_grad(::CompileAndSampleFactorGraph) = false


text:
using FTPClient
using FTPServer
using FTPServer: HOMEDIR, uri
using Memento
using Memento.TestUtils: @test_log
using Test

@testset "FTPServer.jl" begin
    FTPServer.init()

    @testset "no-ssl" begin
        FTPServer.serve() do server
            opts = (
                :hostname => FTPServer.hostname(server),
                :port => FTPServer.port(server),
                :username => FTPServer.username(server),
                :password => FTPServer.password(server),
            )

            options = RequestOptions(; opts..., ssl=false)
            ctxt, resp = ftp_connect(options)
            @test resp.code == 226
        end
    end
    @testset "no-ssl debug true" begin
        @test_log(
            FTPServer.LOGGER,
            "info",
            r"^Running Command:.*",
            FTPServer.Server(; debug_command=true),
        )
    end
    @testset "ssl - $mode - $gen_cert" for mode in (:explicit, :implicit), gen_cert in (true, false)
        FTPServer.serve(; security=mode, force_gen_certs=gen_cert) do server
            opts = (
                :hostname => FTPServer.hostname(server),
                :port => FTPServer.port(server),
                :username => FTPServer.username(server),
                :password => FTPServer.password(server),
                :ssl => true,
                :implicit => mode === :implicit,
                :verify_peer => false,
            )

            options = RequestOptions(; opts...)
            # Test implicit/exlicit ftp ssl scheme is set correctly
            @test options.uri.scheme == (mode === :implicit ? "ftps" : "ftpes")
            ctxt, resp = ftp_connect(options)
            @test resp.code == 226
        end
    end

    @testset "openssl 1.1.1e" begin
        FTPServer.serve(".", security=:implicit) do server
            ftp = FTP(uri(server), verify_peer=false)
            # With OpenSSL 1.1.1e this will fail due to a new EOF change that was added.
            # https://github.com/openssl/openssl/issues/11381
            # https://github.com/openssl/openssl/issues/11378
            @test_broken upload(ftp, "$HOMEDIR/test_download.txt", "/FOO")
        end
    end

    FTPServer.cleanup()

    @test !isfile(FTPServer.CERT)
    @test !isfile(FTPServer.KEY)
    @test !isdir(FTPServer.HOMEDIR)

end


text:
using JuMP, GLPK, LinearAlgebra, Random, HiGHS
using MKL
using LinearOptimizationSolver
using BenchmarkTools, BenchmarkPlots, StatsPlots

optimizer = () -> GLPK.Optimizer()

function build_problem(p::Integer, k::Integer)
    Random.seed!(123)
    e_p = rand(1:5,p)
    cc_p = rand(1:9,p)
    oc_p = rand(1:2,p)
    ic_k = rand(1:30,k)
    dd_k = rand(1:5,k)
    du_k = rand(1:1,k)
    r_k = [dd_k[i]*du_k[i] for i in 1:k]

    lenvars = [p, p*k, k];
    n = sum(lenvars)

    I_pp = Matrix(I,p,p)
    I_kk = Matrix(I,k,k)

    Z_pk = zeros(p,k)
    Z_kp = zeros(k,p)

    D = zeros(p,p*k)
    for i in 1:p
        D[i,1+(i-1)*k:i*k] .= 1
    end

    E = zeros(k,p*k)
    for i in 1:p
        E[:,1+(i-1)*k:i*k] = I_kk
    end

    A = [-I_pp D Z_pk; Z_kp E I_kk; Z_kp -E -I_kk]

    b = vcat(e_p, r_k, -r_k)

    c_y = zeros(p*k)

    for i in 1:p, j in 1:k
        c_y[(i-1)*k+j] = du_k[j]*oc_p[i]
    end

    c = -[cc_p' c_y' ic_k']

    m = p+2*k

    A = [A Matrix(I,m,m)]
    c = vcat(c',zeros(m))
    return A, b, c, n, m
end

p = 30
k = 10
A, b, c, n, m = build_problem(p, k)

input = create(A, b, c, solver = 1, verbose=false)
bench1 = @benchmark output = solve(input)
input = create(A, b, c, solver = 0, verbose=false)
bench2 = @benchmark output = solve(input)

function GLPK!(A, b, c)
    model_pl = Model(GLPK.Optimizer);
    m = size(A,1);
    n = size(A,2) - length(b);
    @variable(model_pl, X[1:n + m] >= 0);
    @constraint(model_pl, A*X .== b);
    @objective(model_pl, Max, c'X);
    optimize!(model_pl);
end

function HiGHS!(A, b, c)
    model_pl = Model(HiGHS.Optimizer);
    set_silent(model_pl)
    m = size(A,1);
    n = size(A,2) - length(b);
    @variable(model_pl, X[1:n + m] >= 0);
    @constraint(model_pl, A*X .== b);
    @objective(model_pl, Max, c'X);
    optimize!(model_pl);
end

bench3 = @benchmark GLPK!($A, $b, $c)
bench4 = @benchmark HiGHS!($A, $b, $c)

plotd = plot(bench1,yaxis=:log10,st=:violin)
plot!(bench2,yaxis=:log10,st=:violin,xticks=(1:2,["IP" "Simplex"]))
plot!(bench3,yaxis=:log10,st=:violin,xticks=(1:3,["IP" "Simplex" "GLPK"]))
plot!(bench4,yaxis=:log10,st=:violin,xticks=(1:4,["IP" "Simplex" "GLPK" "HiGHS"]))

savefig(plotd,"examples\\p=$(p),k=$(k).png")


text:
using SmartShockFinder

datatyp = "NNPP"
data_input_size = 2
n_meshes = 4
polydeg = collect(2:9)

generate_traindataset1d(datatyp, n_meshes, polydeg, data_input_size)


text:
using Gadfly

plot(layer(x=1:10, y=1:10, Geom.point, Theme(default_color=colorant"orange")),
     layer(x=collect(1:10) .+ 0.1, y=collect(1:10) .+ 0.1, Geom.point,
           Theme(default_color=colorant"blue"), order=1))


text:
function StandardDeviation()
    return StandardDeviation(())
end

function StandardDeviation(arg0::SecondMoment)
    return StandardDeviation((SecondMoment,), arg0)
end

function StandardDeviation(arg0::StandardDeviation)
    return StandardDeviation((StandardDeviation,), arg0)
end

function StandardDeviation(arg0::jboolean)
    return StandardDeviation((jboolean,), arg0)
end

function StandardDeviation(arg0::jboolean, arg1::SecondMoment)
    return StandardDeviation((jboolean, SecondMoment), arg0, arg1)
end

function clear(obj::StandardDeviation)
    return jcall(obj, "clear", void, ())
end

function copy(obj::StandardDeviation)
    return jcall(obj, "copy", StandardDeviation, ())
end

function evaluate(obj::StandardDeviation, arg0::Vector{jdouble}, arg1::jdouble)
    return jcall(obj, "evaluate", jdouble, (Vector{jdouble}, jdouble), arg0, arg1)
end

function evaluate(obj::StandardDeviation, arg0::Vector{jdouble}, arg1::jdouble, arg2::jint, arg3::jint)
    return jcall(obj, "evaluate", jdouble, (Vector{jdouble}, jdouble, jint, jint), arg0, arg1, arg2, arg3)
end

function evaluate(obj::StandardDeviation, arg0::Vector{jdouble}, arg1::jint, arg2::jint)
    return jcall(obj, "evaluate", jdouble, (Vector{jdouble}, jint, jint), arg0, arg1, arg2)
end

function get_n(obj::StandardDeviation)
    return jcall(obj, "getN", jlong, ())
end

function get_result(obj::StandardDeviation)
    return jcall(obj, "getResult", jdouble, ())
end

function increment(obj::StandardDeviation, arg0::jdouble)
    return jcall(obj, "increment", void, (jdouble,), arg0)
end

function is_bias_corrected(obj::StandardDeviation)
    return jcall(obj, "isBiasCorrected", jboolean, ())
end

function with_bias_correction(obj::StandardDeviation, arg0::jboolean)
    return jcall(obj, "withBiasCorrection", StandardDeviation, (jboolean,), arg0)
end


text:
using Pkg
Pkg.activate(@__DIR__)

using SigmaRidgeRegression
using LinearAlgebra
using StatsBase
using Plots
using MLJ
using LaTeXStrings
using Random

using PGFPlotsX
pgfplotsx()
push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\usepackage{amsmath}")
push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\usepackage{amssymb}")
push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\usepackage{bm}")
push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\usepackage[bbgreekl]{mathbbol}")
push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\newcommand{\sigmacv}{\bbsigma}")
push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\newcommand{\loo}{\operatorname{CV}^*}")
push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\newcommand{\blambda}{\bm{\lambda}}")

main_cols = [:purple :green :grey]

Random.seed!(1)
σ = 4.0
grp = GroupedFeatures([25; 25; 25])
n = 400
p = grp.p

X = randn(n, p)

αs = sqrt.([4.0; 8.0; 12.0])
β = random_betas(grp, αs)
#group_summary(grp, β, norm)
#sum(abs2, β)
Y = X * β .+ σ .* randn(n)

loo_sigma_ridge = LooSigmaRidgeRegressor(; groups=grp, center=false, scale=false)

loo_sigmaridge_machine = machine(loo_sigma_ridge, MLJ.table(X), Y)

fit!(loo_sigmaridge_machine)

σs = loo_sigmaridge_machine.report.params

#  \sigmacv <= 6.2
param_subset = loo_sigmaridge_machine.report.params .<= 6.2
λs = Matrix(hcat(loo_sigmaridge_machine.report.λs...)')
λs_subset = λs[param_subset, :]

βs_list = [
    fit!(
        machine(
            SigmaRidgeRegressor(; groups=grp, σ=σ, center=false, scale=false),
            MLJ.table(X),
            Y,
        ),
    ).fitresult.coef for σ in σs
]
βs = Matrix(hcat(βs_list...)')

pl_left = plot(
    loo_sigmaridge_machine.report.params[param_subset],
    λs_subset;
    legend=:topleft,
    linecolor=main_cols,
    linestyle=[:dot :dashdot :dash],
    xlab=L"\sigmacv",
    ylab=L"\widehat{\lambda}_g(\sigmacv)",
    background_color_legend=:transparent,
    foreground_color_legend=:transparent,
    grid=nothing,
    frame=:box,
    label=["Group 1" "Group 2" "Group 3"],
    thickness_scaling=1.8,
    size=(550, 400),
)
savefig(pl_left, "intro_figure_left.tikz")

pl_right = plot(
    loo_sigmaridge_machine.report.params[param_subset],
    loo_sigmaridge_machine.report.loos[param_subset];
    color=:darkblue,
    grid=nothing,
    frame=:box,
    background_color_legend=:transparent,
    foreground_color_legend=:transparent,
    xlab=L"\sigmacv",
    ylab=L"\loo(\sigmacv)",
    label=nothing,
    thickness_scaling=1.8,
    size=(550, 400),
    ylim=(18.2, 19.2),
    yticks=[18.4; 18.6; 18.8; 19.0],
)

savefig(pl_right, "intro_figure_right.tikz")

pl_left_tree = plot(
    loo_sigmaridge_machine.report.params,
    λs;
    ylim=(0, 15),
    legend=:topleft,
    color=main_cols,
    linestyle=[:dot :dashdot :dash],
    background_color_legend=:transparent,
    foreground_color_legend=:transparent,
    grid=nothing,
    frame=:box,
    xlab=L"\sigmacv",
    ylab=L"\widehat{\lambda}_g(\sigmacv)",
    label=["Group 1" "Group 2" "Group 3"],
    thickness_scaling=1.8,
    size=(550, 400),
)

savefig(pl_left_tree, "christmas_tree_left.tikz")

cols_rep = hcat([fill(col, 1, 25) for col in [:purple; :green; :grey]]...)
ltype_rep = hcat([fill(lt, 1, 25) for lt in [:solid; :dash; :dot]]...)

pl_right_tree = plot(
    σs,
    βs;
    alpha=0.6,
    linewidth=0.8,
    label="",
    grid=nothing,
    frame=:box,
    ylab=L"\widehat{w}_j(\widehat{\blambda}(\sigmacv))",
    xlab=L"\sigmacv",
    color=cols_rep,
    ylim=(-2.1, 2.1),
    thickness_scaling=1.8,
    size=(550, 400),
)

savefig(pl_right_tree, "christmas_tree_right.tikz")


text:
import SimpleMolecule.AbstractMolecularContainer

"""
A test for constructor_1 for ``Atom``.

Arguments
----------
props:Dict
    property dictionary

solution:Atom
    answer key
"""
function test_constructor{T<:AbstractMolecularContainer}(
    props::Dict,
    solution::T)
    println("-------------------------------------------------")
    print_with_color(:blue, "Test constructor\n\n")

    answer = T(; props...)

    for field in fieldnames(solution)
        if obtain(answer, field) != obtain(solution, field)
            println("field: ", field)
            println("\tanswer ", obtain(answer, field))
            println("\tsolution ", obtain(solution, field))
            @test false
        end
    end

    print_with_color(:green, "VERIFIED! constructor 1 for Atom\n\n")
    println("-------------------------------------------------")

    return true
end


text:
using QuantumAlgebra
using QuantumAlgebra: δ, QuExpr, QuTerm, BaseOpProduct, BaseOperator, Param, QuIndex, _map_quexpr_ops, TLSx, TLSCreate, is_normal_form
using Test, Documenter

function myδ(i,j)
    iA,iB = QuIndex.((i,j))
    QuExpr(QuTerm([δ(min(iA,iB),max(iA,iB))],BaseOpProduct()))
end
scal(x) = x*one(QuExpr)

macro test_is_normal_form(x)
    x = esc(x)
    quote
        xn = normal_form($x)
        @test is_normal_form($x) == ($x == xn)
        @test is_normal_form(xn)
    end
end

@time @testset "QuantumAlgebra.jl" begin
    doctest(QuantumAlgebra)

    @test isbitstype(QuantumAlgebra.QuIndex)
    @test isbitstype(QuantumAlgebra.QuOpName)
    if isbitstype(QuantumAlgebra.QuIndices)
        @test isbitstype(QuantumAlgebra.BaseOperator)
        @test isbitstype(QuantumAlgebra.Param)
    end
    @test isbitstype(QuantumAlgebra.δ)

    @testset "auto_normal_form($auto_norm)" for auto_norm in (false,true)
        QuantumAlgebra.auto_normal_form(auto_norm)

        @test QuantumAlgebra._auto_normal_form[] == auto_norm

        @boson_ops b
        @test normal_form(b(:i)*bdag(:j)) == bdag(:j)*b(:i) + myδ(:i,:j)
        @test normal_form(b(:i)*adag(:j)) == adag(:j)*b(:i)

        @fermion_ops c
        @test normal_form(c(:i)*cdag(:j) + cdag(:j)*c(:i)) == myδ(:i,:j)
        # different fermionic species commute by default
        @test normal_form(c(:i)*fdag(:j)) == fdag(:j)*c(:i)

        # this creates two fermionic operators belonging to the same species (which anticommute, e.g., for itinerant and localized electrons)
        # NB: sort order is the order in which names are given (not alphabetic like normally)
        @anticommuting_fermion_group e d
        # fermionic operators of the same species anticommute
        @test iszero(normal_form(d()*e() + e()*d()))
        # ...but different fermionic species commute
        @test normal_form(f()*e() + e()*f()) == 2*e()*f()

        # all equal numbers should be equal scalars (ignore type)
        @test a() + 0 == a() + 0.0
        @test a() + 1 == a() + (1//1 + 0im)
        @test scal(1.5) == scal(3//2)

        # test params and their complex conjugation
        @test Pc"ω" == param(:ω,'n')
        @test Pc"ω_i" == param(:ω,'n',:i)
        @test Pc"ω_i,j,2" == param(:ω,'n',:i,:j,2)
        @test_throws ArgumentError param(:ω,'g')
        @test_throws ArgumentError param(:ω,'n',2,:i,"ag")
        @test_throws ArgumentError a("aa")

        @test adjoint(param(:g,'r')) == param(:g,'r')
        @test adjoint(param(:g,'n')) == param(:g,'c')
        @test adjoint(param(:g,'n',:i)) == param(:g,'c',:i)
        @test (Pc"g_i")' == param(:g,'c',:i)
        @test adjoint(param(:g,'c',:i,1,:m_3)) == param(:g,'n',:i,1,:m_3)

        tmp1 = param(:g,'n',2)*param(:g,'n',1)*param(:g,'c',3)*param(:b,'n')*param(:a,'n')*param(:d,'c',3)*param(:f,'r',1,:i)
        tmpc = param(:g,'c',2)*param(:g,'c',1)*param(:g,'n',3)*param(:b,'c')*param(:a,'c')*param(:d,'n',3)*param(:f,'r',1,:i)
        tmp2 = Pc"a"*Pc"b"*Pc"g_1"*Pc"g_2"*(Pc"g_3")' * (Pc"d_3")' * Pr"f_1,i"
        @test normal_form(tmp1) == normal_form(tmp2)
        @test normal_form(tmp1') == normal_form(tmpc)

        @test a()*5 == 5*a()
        @test a()+5 == 5+a()
        @test 5-a() == -a() + 5
        @test a()-5 == -5 + a()

        @test a()' == adag()
        @test normal_form(a(:m)*adag(:m)) == adag(:m)*a(:m) + 1
        @test normal_form((a(:m)*a(1))') == adag(1)*adag(:m)
        @test normal_form((a(1,2,:k)*a(:m))') == normal_form(adag(1,2,:k)*adag(:m))

        @test normal_form(fdag(:a)*f(:b) + f(:b)*fdag(:a)) == myδ(:a,:b)
        @test iszero(normal_form(f(:a)*f(:b) + f(:b)*f(:a)))
        @test iszero(normal_form(fdag(:a)*fdag(:b) + fdag(:b)*fdag(:a)))

        @test f()' == fdag()
        @test fdag()' == f()
        @test (fdag()*f())' == fdag()*f()
        @test comm(fdag(),f()) == -comm(f(),fdag())

        @test ∑(:j,a(:j))*∑(:i,a(:i)) == ∑(:i,∑(:j,a(:i)*a(:j)))

        tmp = ∑(:i,adag(:i)*a(:i))
        @test tmp' == tmp
        @test normal_form(a(:i)*∑(:i,a(:i))) == normal_form(∑(:i_1,a(:i_1)*a(:i)))
        @test normal_form(adag(:n)*tmp) == normal_form(∑(:i,adag(:n)*adag(:i)*a(:i)))
        @test normal_form(a(:n)   *tmp) == normal_form(∑(:i,adag(:i)*a(:n)*a(:i)) + a(:n))
        @test normal_form(param(:g,'n',:i)*∑(:i,a(:i))) == normal_form(∑(:i_1,param(:g,'n',:i)*a(:i_1)))
        @test normal_form(param(:g,'n',:i_1)*a(:i)*∑(:i,a(:i))) == normal_form(∑(:i_2,param(:g,'n',:i_1)*a(:i)*a(:i_2)))
        @test normal_form(param(:g,'n',:n)*∑(:i,a(:i))) == normal_form(∑(:i,param(:g,'n',:n)*a(:i)))

        @test Pr"gz_i,μ" == param(:gz,'r',(:i,:μ))
        @test Pc"gz_i,μ" == param(:gz,'n',(:i,:μ))
        @test (Pc"gz_i,μ")' == param(:gz,'c',(:i,:μ))
        @test Pr"α_i_1,μ_2" == Pr"α_{i_1,μ_2}" == param(:α,'r',:i_1,:μ_2)
        @test Pr"α_◔,◔" == Pr"α_{◔,◔}" == param(:α,'r',:◔,:◔)

        @testset "normal_form with sums" begin
            # bug report from FJ Matute
            t1 = Pc"a_w"*Pc"B_x"*Pc"C_y"*Pc"D_z"*f(:w)*f(:x)'*f(:y)*f(:z)'
            X1 = ∑(:w,∑(:x,∑(:y,∑(:z,t1))))
            X2 = ∑(:w,∑(:x,∑(:y,∑(:z,normal_form(t1)))))
            @test normal_form(X1) == normal_form(X2)
        end

        @testset "use_σpm($with_σpm)" for with_σpm in (false,true)
            QuantumAlgebra.use_σpm(with_σpm)

            @test QuantumAlgebra.using_σpm() == with_σpm

            @test_throws ArgumentError σx(:i,"aa")
            @test σx(:i,:b) == σx((:i,:b))
            @test σx(:i,:j) != σx(:i,:b)

            α = param(:α)
            tmp = α' * a(1) + α * a(1)
            @test iszero(tmp - tmp)

            @test σx() == σp() + σm()
            @test σx(:i) == σp(:i) + σm(:i)
            @test σx(:i,:j) == σp((:i,:j)) + σm((:i,:j))

            @test normal_form(σy(:i)) == normal_form(-1im*(σp(:i) - σm(:i)))
            @test normal_form(σz(:i)) == normal_form(2*σp(:i)*σm(:i) - 1)

            for s=(σx,σy,σz)
                @test isone(normal_form(s(:i)*s(:i)))
                @test s(:i)' == s(:i)
            end
            @test normal_form(σx(:i)*σy(:i)) == 1im*σz(:i)
            @test normal_form(σy(:i)*σx(:i)) == -1im*σz(:i)
            @test normal_form(σy()*σz()) == 1im*σx()
            @test normal_form(σx(:i,:j)*σz(:i,:j)) == -1im*σy(:i,:j)
            @test normal_form(σx(:i,:j)*σz(:i,:j)) != -1im*σy(:i,:k)

            @test normal_form(a(1) * (σy(1) * a(1))') == normal_form(a(1) * (adag(1) * σy(1))) == normal_form(adag(1)*a(1)*σy(1) + σy(1))

            @test normal_form(σy(:i)*σy(:i)*σy(:i)*σy(:i)*σy(:i)*σy(:i)*σx(:i)) == σx(:i)
            @test normal_form(σy(:i)*σy(:i)*σy(:i)*σy(:i)*σy(:i)*σy(:i)*σx(:j)) == σx(:j)

            @test_throws ArgumentError normal_form(QuExpr(TLSx(:σ)) * QuExpr(TLSCreate(:σ)))

            tmp1 = normal_form(3*Pc"ω"*Pc"g"*expval(σp(:k))*σp(:k)*adag(5)*a(5))
            if QuantumAlgebra.using_σpm()
                @test length(tmp1.terms) == 1
                @test first(tmp1.terms)[1].bares == first((adag(5)*σp(:k)*a(5)).terms)[1].bares
            else
                @test length(tmp1.terms) == 4
            end

            x1 = ∑(:i,∑(:j,a(:i)*a(:j))) + f(:j,:k,:l)*fdag(:j,:i,:m) + a(3,8,:i)*expval(a(:i)*adag(:k))*adag(:α,:μ,:ν)
            x2 = σp(:k) + σx(:m)
            @test normal_form(x1 + x2) == normal_form(x1) + normal_form(x2)
            @test normal_form(x1 * x2) == normal_form(normal_form(x1) * normal_form(x2))
            @test normal_form(x1 * x2 + x2) == normal_form(normal_form(x1) * normal_form(x2)) + normal_form(x2)

            @testset "is_normal_form" begin
                y = one(QuExpr)
                for x in (σx(:i)*σy(:i), a(1) * (σy(1) * a(1))', a(:d)*adag(:c), a(:i)*adag(:j)*σz(:α)*σy(:β)*σx(:α), expval(adag(:c)*a(:d)))
                    @test_is_normal_form x
                    y *= x
                    @test_is_normal_form y
                    @test_is_normal_form ∑(:c,∑(:i,∑(:α,y)))
                end
            end

            @testset "exponent" begin
                for x in (σx(:i)*σy(:i), a(1) * (σy(1) * a(1))', a(:d)*adag(:c), a(:i)*adag(:j)*σz(:α)*σy(:β)*σx(:α), expval(adag(:c)*a(:d)))
                    @test isone(x^0)
                    @test x^1 == x
                    @test x^2 == x*x
                    @test x^3 == x*x*x
                    @test x^4 == x*x*x*x
                    @test_throws ArgumentError x^-1
                    @test_throws ArgumentError x^-2
                end

            end

            @testset "Commutation inside ExpVal/Corr" begin
                for (x1,x2,s) in ((a(:d)*adag(:c), a(:i)*adag(:j)*σz(:α)*σy(:β)*σx(:α),  expval(adag(:c)*a(:d))),
                                (f(:d)*fdag(:c), f(:i)*fdag(:j)*σz(:α)*σy(:β)*σx(:α), -expval(fdag(:c)*f(:d))))
                    x = ∑(:c,∑(:i,∑(:α,3*expval(x1)*x2)))
                    @test normal_form(expval(x1)) == myδ(:c,:d) + s
                    @test normal_form(3*expval(x1)) == expval(normal_form(3*x1))
                    @test normal_form(expval(x1)*expval(x1)) == normal_form(normal_form(expval(x1))*normal_form(expval(x1)))
                    @test normal_form(3*expval(x2)) == expval(normal_form(3*x2))
                    @test normal_form(3*expval(x1)*x2) == expval(3*normal_form(x1))*normal_form(x2)
                    @test normal_form(3*expval(x1)*corr(x2)) == normal_form(expval(3*normal_form(x1))*corr(normal_form(x2)))
                    @test normal_form(x) == normal_form(∑(:α,∑(:c,∑(:i,3*expval(normal_form(x1))*normal_form(x2)))))
                    xn = normal_form(x)
                    xsqn = normal_form(x*x)
                    @test xsqn == normal_form(xn*xn)

                    # do normal_form already on the first partial product, otherwise computation time explodes.
                    # NOTE: this test is broken not because of an error as far as I can tell, but because
                    # depending on the order of computation it produces equivalent terms that we cannot
                    # identify as equal yet, such as:
                    # ∑₁₂ σ⁺(#₁) σ⁺(#₂) σ⁺(β) σ⁻(#₁) and
                    # ∑₁₂ σ⁺(#₁) σ⁺(#₂) σ⁺(β) σ⁻(#₂)
                    #
                    # For use_σxyz, this is even worse, since even though non-adjacent operators are contracted,
                    # there are ambiguities in sums that even lead to different numbers of terms
                    if auto_norm
                        @test normal_form(xsqn*x) == normal_form(xn*xn*xn)
                    else
                        @test_broken normal_form(xsqn*x) == normal_form(xn*xn*xn)
                    end
                end
            end

            @test normal_form(myδ(:i,:k)*a(:k)) == normal_form(a(:i)*myδ(:k,:i))
            @test normal_form(myδ(:i,:k)*a(:k,:i)) == normal_form(a(:i,:i)*myδ(:k,:i))
            @test normal_form(myδ(:i,:k)*myδ(:i,:j)) == normal_form(myδ(:k,:i)*myδ(:j,:k))
            # k cannot be equal to 1 and 3 at the same time
            @test iszero(normal_form(myδ(1,:k)*myδ(:k,3)*σx(:k)))
            @test normal_form(myδ(1,:k)*myδ(:k,1)*σx(:k)) == myδ(1,:k)*σx(1)

            @test iszero(normal_form(comm(σx(5),σy(3))))
            @test iszero(normal_form(comm(σx(5),σx(5))))
            @test normal_form(comm(σx(1),σz(1))) == normal_form(-2im*σy(1))
            @test normal_form(comm(σx(:μ),σy(:ν))) == normal_form(2im*myδ(:μ,:ν)*σz(:ν))
            @test normal_form(1//2im * comm(σx(:m),σy(:m))) == normal_form(σz(:m))
            @test normal_form(σx(:a)*σy(:a)*σz(:a)) == 1im*one(QuExpr)

            @test iszero(normal_form(comm(Pc"g",a(5)+a(3))))
            @test iszero(normal_form(comm(Pc"g",a(5)*a(3))))
            @test iszero(normal_form(comm(a(5)+a(3),Pc"g")))
            @test iszero(normal_form(comm(a(5)*a(3),Pc"g")))

            @test normal_form(comm(a(5)+a(3),adag(5))) == normal_form(comm(a(5),adag(5))+ comm(a(3),adag(5)))
            @test normal_form(comm(a(5)+a(3),adag(5)*a(3))) == normal_form(comm(a(5),adag(5)*a(3))+ comm(a(3),adag(5)*a(3)))

            @test normal_form(adag(2)*σy(:i) - 14Pc"ω") == normal_form(-14Pc"ω" + σy(:i)*adag(2))
            @test normal_form(adag(2)*σy(:i)) == normal_form(σy(:i)*adag(2))

            @test normal_form(σp(1) * σm(1)) == normal_form(1//2*σz(1) + 1//2)
            @test normal_form(σm(1)*σp(1)) == normal_form(1 - σp(1)*σm(1))
            @test normal_form(σm(1)*σp(1)) == normal_form(σp(1)*σm(1) + comm(σm(1),σp(1)))
            @test normal_form(σz(1)) == normal_form(σp(1)*σm(1) - σm(1)*σp(1))

            @test iszero(comm(σp(1),σp(1)))
            @test normal_form(comm(σp(:n),σm(:n))) == normal_form(σz(:n))

            @test normal_form(comm(a(1),   adag(1)*a(1))) == a(1)
            @test normal_form(comm(adag(1),adag(1)*a(1))) == -adag(1)

            @test expval(scal(3)) == scal(3)
            @test corr(scal(3)) == scal(3)
            @test corr(3+adag(2)) == 3 + corr(adag(2))

            @test expval_as_corrs(scal(3)) == scal(3)
            @test expval_as_corrs(a(2)) == corr(a(2))
            @test expval_as_corrs(3Pc"g") == 3Pc"g"
            @test expval_as_corrs(3a(2)) == 3corr(a(2))
            @test expval_as_corrs(3+adag(2)) == 3 + corr(adag(2))

            @test expval_as_corrs(a(2)*a(2)) == corr(a(2)*a(2)) + corr(a(2))*corr(a(2))
            @test expval_as_corrs(a(2)*a(:m)) == corr(a(2)*a(:m)) + corr(a(2))*corr(a(:m))
            tmpas = a.(1:3)
            @test *(tmpas...) == a(1)*a(2)*a(3)
            tmpEVs = corr.(tmpas)
            @test expval_as_corrs(*(8,Pr"g",tmpas...)) == 8*Pr"g" * (corr(*(tmpas...)) + *(tmpEVs...) + tmpEVs[1]*corr(tmpas[2]*tmpas[3]) + tmpEVs[2]*corr(tmpas[1]*tmpas[3]) + tmpEVs[3]*corr(tmpas[1]*tmpas[2]))

            if QuantumAlgebra.using_σpm()
                @test expval_as_corrs(-param(:g,'r',1)*σp(1)) == -param(:g,'r',1)*corr(σp(1))
                @test expval_as_corrs(∑(:i,σp(:i)*σm(:n))) == ∑(:i,corr(σp(:i)*σm(:n))) + ∑(:i,corr(σp(:i))*corr(σm(:n)))
            else
                @test expval_as_corrs(-param(:g,'r',1)*σz(1)) == -param(:g,'r',1)*corr(σz(1))
                @test expval_as_corrs(∑(:i,σy(:i)*σy(:n))) == ∑(:i,corr(σy(:i)*σy(:n))) + ∑(:i,corr(σy(:i))*corr(σy(:n))) - corr(σy(:n))*corr(σy(:n))
            end

            H = ∑(:i,param(:ω,'r',:i)*adag(:i)*a(:i))
            @test normal_form(comm(a(:i),H)) == param(:ω,'r',:i)*a(:i)
            @test normal_form(comm(a(:n),H)) == param(:ω,'r',:n)*a(:n)
            @test normal_form(comm(H,a(:n))) == -param(:ω,'r',:n)*a(:n)
            @test normal_form(comm(adag(:n),H)) == -param(:ω,'r',:n)*adag(:n)
            @test normal_form(comm(adag(:n)*a(:m),H)) == (param(:ω,'r',:m)-param(:ω,'r',:n))*adag(:n)*a(:m)

            @test normal_form(a()*H) == normal_form(∑(:i,param(:ω,'r',:i)*adag(:i)*a(:i)*a()))
            @test normal_form(a(:k)*H) == normal_form(param(:ω,'r',:k)*a(:k) + ∑(:i,param(:ω,'r',:i)*adag(:i)*a(:i)*a(:k)))
            HH = ∑(:i,param(:ω,'r',:i,:i)*adag(:i,:i)*a(:i,:i))
            @test normal_form(a(:k,:k)*HH) == normal_form(param(:ω,'r',:k,:k)*a(:k,:k) + ∑(:i,param(:ω,'r',:i,:i)*adag(:i,:i)*a(:i,:i)*a(:k,:k)))

            @test iszero(Avac(H))
            @test iszero(vacA(H))
            @test iszero(vacA(adag(3)*σp(1)*σm(1)))
            @test iszero(vacA(fdag(:n)))
            @test vacA(f(:n)) == f(:n)
            @test vacA(a(:n)) == a(:n)
            @test Avac(fdag(:n)) == fdag(:n)
            @test iszero(Avac(f(:n)))
            @test iszero(Avac(a(3)*σp(1)*σm(1)))
            @test isone(Avac(σm(1)*σp(1)))
            @test iszero(Avac(σp(1)*σm(1)))
            @test Avac(σp(1)) == σp(1)
            @test vacA(σm(1)) == σm(1)
            if QuantumAlgebra.using_σpm()
                @test iszero(Avac(σm(1)))
                @test iszero(vacA(σp(1)))
            else
                @test Avac(σx(1)) == vacA(σx(1)) == σx(1)
            end
            @test iszero(vacExpVal(σx(1)))
            @test iszero(vacExpVal(σp(1)))
            @test iszero(vacExpVal(∑(:i,σp(:i))))
            @test iszero(vacExpVal(σp(:i)*σm(:k)))
            @test vacExpVal(σm(:j)*σp(:l)) == myδ(:j,:l)
            if QuantumAlgebra.using_σpm()
                @test normal_form(σp(:i)*σm(:k)) == σp(:i)*σm(:k)
            end
            @test normal_form(a(:n)*adag(:n)*a(:n)*adag(:n)) == scal(1) + 3adag(:n)*a(:n) + adag(:n)*adag(:n)*a(:n)*a(:n)

            S = (1/√(2*6))*adag(:n)*adag(:n)*adag(:n) + (1/√2)*adag(:m)
            for (A,val) in [(scal(1),scal(1)),
                            (adag(:n)*a(:n),scal(1.5) + 0.5 * myδ(:n,:m)),
                            (adag(:n)*adag(:n)*a(:n)*a(:n),scal(3))]
                @test vacExpVal(A,S) ≈ val
            end

            tmp = ∑(:i,expval_as_corrs(adag(:n)*a(:i)))
            tmplatex = raw"\sum_{\#_{1}}\langle {a}_{n}^\dagger \rangle_{c} \langle {a}_{\#_{1}} \rangle_{c} + \sum_{\#_{1}}\langle {a}_{n}^\dagger {a}_{\#_{1}} \rangle_{c}"
            @test latex(tmp) == tmplatex
            @test expval_as_corrs(tmp) == tmp
            @test sprint(show,"text/latex",tmp) == "\$$(tmplatex)\$"

            @test latex(a(:i_1,:i_2)) == raw"{a}_{i_{1}i_{2}}"

            @test latex(normal_form(e(:i)*ddag(:j))) == "-{d}_{j}^\\dagger {e}_{i}"
            @test latex(normal_form(f(:i)*ddag(:j))) == "{d}_{j}^\\dagger {f}_{i}"

            lσp = latex(σp())
            lσz = latex(σz())
            if QuantumAlgebra.using_σpm()
                @test lσp == "{{\\sigma}}^+"
                @test lσz == "-1 + 2{{\\sigma}}^+ {{\\sigma}}^-"
            else
                @test lσp == "\\frac{1}{2}{{\\sigma}}^x + \\frac{1}{2}i{{\\sigma}}^y"
                @test lσz == "{{\\sigma}}^z"
            end

            @test QuantumAlgebra.symmetric_index_nums(adag(:i)*adag(:j)*a(:k)*a(:l)) == [2,2]

            @test string(normal_form(f(:i) * ddag(:j))) == "d†(j) f(i)"
            @test string(normal_form(e(:i) * ddag(:j))) == "-d†(j) e(i)"

            @test string(normal_form(∑(:i,a(:i)) * adag(:n))) == "1 + ∑₁ a†(n) a(#₁)"
            @test string(5a(:i) + 3adag(:j) - 3f(1)) == "3 a†(j) - 3 f(1) + 5 a(i)"
            @test string(5im*a(:i) + (3+2im)*adag(:j) - 3//2*f(1)) == "(3+2i) a†(j) - 3//2 f(1) + 5i a(i)"
            if QuantumAlgebra.using_σpm()
                @test string(normal_form(a(5)*adag(5)*σp(3)*expval_as_corrs(adag(5,:i)*a(5)))) == "⟨a†(5i)⟩c ⟨a(5)⟩c σ⁺(3) + ⟨a†(5i) a(5)⟩c σ⁺(3) + ⟨a†(5i)⟩c ⟨a(5)⟩c a†(5) σ⁺(3) a(5) + ⟨a†(5i) a(5)⟩c a†(5) σ⁺(3) a(5)"
            else
                @test string(normal_form(a(5)*adag(5)*σz(3)*expval_as_corrs(adag(5,:i)*a(5)))) == "⟨a†(5i)⟩c ⟨a(5)⟩c σᶻ(3) + ⟨a†(5i) a(5)⟩c σᶻ(3) + ⟨a†(5i)⟩c ⟨a(5)⟩c a†(5) σᶻ(3) a(5) + ⟨a†(5i) a(5)⟩c a†(5) σᶻ(3) a(5)"
            end

            @test julia_expression(QuExpr()) == 0
            # an empty QuTerm is the identity operator!
            @test julia_expression(QuTerm()) == 1

            x = ∑(:i,Pc"g_i,k"*a(:i,:j_2,:K)*adag(:i_1,:J,:k)*σp(:i))
            ex = julia_expression(expval(normal_form(x)))
            if QuantumAlgebra.using_σpm()
                @test ex == :(I[J, j₂] * I[K, k] * g[i₁, K] * σ⁺[i₁] + g[s̄₁, k] * aᴴσ⁺a[i₁, J, k, s̄₁, s̄₁, j₂, K])
            else
                # julia_expression interpolates complex numbers directly as complex, not as expressions. so make sure to do the same here
                halfim = 0.5im
                @test ex == :(0.5 * I[J, j₂] * I[K, k] * g[i₁, K] * σˣ[i₁] + $halfim * I[J, j₂] * I[K, k] * g[i₁, K] * σʸ[i₁] + 0.5 * g[s̄₁, k] * aᴴσˣa[i₁, J, k, s̄₁, s̄₁, j₂, K] + $halfim * g[s̄₁, k] * aᴴσʸa[i₁, J, k, s̄₁, s̄₁, j₂, K])
            end
        end
    end
end


text:
using Random: randn

@testset "detrend()" begin
    N = 1000
    X = randn(N, 2)
    Y = copy(X)
    for i in 1:N
        Y[i,:] = X[i,:] .+ 1 .+ (2*i)
    end

    @test X ≈ detrend(Y) rtol=3/sqrt(N);

    detrend!(Y)
    @test X ≈ Y rtol=3/sqrt(N)
end


text:
#    Copyright (c) 2021 J.A. Duffek
#    Copyright (c) 2000 D.M. Spink
#
#    This program is free software; you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation; either version 2 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with this program; if not, see <http://www.gnu.org/licenses/>.

"""
    bspeval(d::I,c::Matrix{F},k::Vector{F},u::Vector{F}
           )::Matrix{F} where {I<:Integer,F<:AbstractFloat}

Evaluate a univariate B-Spline.

# Arguments:
- `d`: degree of the B-Spline
- `c`: control points, matrix of size `(dim,nc)`
- `k`: knot sequence, vector of size `nk`
- `u`: parametric evaluation points, vector of size `nu`

# Output:
- `p`: Evaluated points, matrix of size `(dim,nu)`

# Examples:
```julia
julia> p = bspeval(d,c,k,u)
```

# Reference:
- Algorithm A3.1 from 'The NURBS BOOK' pg82 ([1997Piegl](@cite))
"""
function bspeval(d::I,c::Matrix{F},k::Vector{F},u::Vector{F}
                )::Matrix{F} where {I<:Integer,F<:AbstractFloat}
# TODO: the original C function passed c, k, u and p as reference
# TODO: this function has optimization potential

nu = length(u);
mc, nc = size(c);

# allocate space for evaluated points
p = zeros(F,(mc,nu));

# allocate space for the basis functions
N = zeros(F,d+1);

# for each parametric point
for col in 1:nu
  # find the span of u[col]
  s = findspan(nc-1,d,u[col],k);
  N = basisfun(s,u[col],d,k);

  tmp1 = s - d + 1;
  for row in 1:mc
    tmp2 = 0;
    for i in 0:d
      tmp2 += N[i+1] * c[row,tmp1+i];
    end # for i
    p[row,col] = tmp2;
  end # for row
end # for col
return p;
end # bspeval


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct QueryTroubleshootingParameters <: SwaggerModel
    targetResourceId::Any # spec type: Union{ Nothing, String } # spec name: targetResourceId

    function QueryTroubleshootingParameters(;targetResourceId=nothing)
        o = new()
        validate_property(QueryTroubleshootingParameters, Symbol("targetResourceId"), targetResourceId)
        setfield!(o, Symbol("targetResourceId"), targetResourceId)
        o
    end
end # type QueryTroubleshootingParameters

const _property_map_QueryTroubleshootingParameters = Dict{Symbol,Symbol}(Symbol("targetResourceId")=>Symbol("targetResourceId"))
const _property_types_QueryTroubleshootingParameters = Dict{Symbol,String}(Symbol("targetResourceId")=>"String")
Base.propertynames(::Type{ QueryTroubleshootingParameters }) = collect(keys(_property_map_QueryTroubleshootingParameters))
Swagger.property_type(::Type{ QueryTroubleshootingParameters }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_QueryTroubleshootingParameters[name]))}
Swagger.field_name(::Type{ QueryTroubleshootingParameters }, property_name::Symbol) =  _property_map_QueryTroubleshootingParameters[property_name]

function check_required(o::QueryTroubleshootingParameters)
    (getproperty(o, Symbol("targetResourceId")) === nothing) && (return false)
    true
end

function validate_property(::Type{ QueryTroubleshootingParameters }, name::Symbol, val)
end