text:
using AutoHashEquals
using MathProgBase
using Serialization

const module_tempdir = joinpath(Base.tempdir(), "julia", string(nameof(@__MODULE__)))
const model_dir = joinpath(module_tempdir, "models")
if !ispath(model_dir)
    mkpath(model_dir)
end

function remove_cached_models()
    if ispath(model_dir)
        rm(model_dir, recursive = true)
        mkpath(model_dir)
    end
end

"""
Supertype for types encoding the family of perturbations allowed.
"""
abstract type PerturbationFamily end

struct UnrestrictedPerturbationFamily <: PerturbationFamily end
Base.show(io::IO, pp::UnrestrictedPerturbationFamily) = print(io, "unrestricted")
Base.hash(a::UnrestrictedPerturbationFamily, h::UInt) = hash(:UnrestrictedPerturbationFamily, h)

abstract type RestrictedPerturbationFamily <: PerturbationFamily end

"""
For blurring perturbations, we currently allow colors to "bleed" across color channels -
that is, the value of the output of channel 1 can depend on the input to all channels.
(This is something that is worth reconsidering if we are working on color input).
"""
@auto_hash_equals struct BlurringPerturbationFamily <: RestrictedPerturbationFamily
    blur_kernel_size::NTuple{2}
end
Base.show(io::IO, pp::BlurringPerturbationFamily) =
    print(io, filter(x -> !isspace(x), "blur-$(pp.blur_kernel_size)"))

@auto_hash_equals struct LInfNormBoundedPerturbationFamily <: RestrictedPerturbationFamily
    norm_bound::Real

    function LInfNormBoundedPerturbationFamily(norm_bound::Real)
        @assert(norm_bound > 0, "Norm bound $(norm_bound) should be positive")
        return new(norm_bound)
    end
end
Base.show(io::IO, pp::LInfNormBoundedPerturbationFamily) =
    print(io, "linf-norm-bounded-$(pp.norm_bound)")

function get_model(
    nn::NeuralNet,
    input::Array{<:Real},
    pp::UnrestrictedPerturbationFamily,
    tightening_solver::MathProgBase.SolverInterface.AbstractMathProgSolver,
    tightening_algorithm::TighteningAlgorithm,
    rebuild::Bool,
    cache_model::Bool,
)::Dict
    d = get_reusable_model(
        nn,
        input,
        pp,
        tightening_solver,
        tightening_algorithm,
        rebuild,
        cache_model,
    )
    @constraint(d[:Model], d[:Input] .== input)
    delete!(d, :Input)
    # NOTE (vtjeng): It is important to set the solver before attempting to add a
    # constraint, as the saved model may have been saved with a different solver (or
    # different) environment. Flipping the order of the two leads to test failures.
    return d
end

function get_model(
    nn::NeuralNet,
    input::Array{<:Real},
    pp::RestrictedPerturbationFamily,
    tightening_solver::MathProgBase.SolverInterface.AbstractMathProgSolver,
    tightening_algorithm::TighteningAlgorithm,
    rebuild::Bool,
    cache_model::Bool,
)::Dict
    d = get_reusable_model(
        nn,
        input,
        pp,
        tightening_solver,
        tightening_algorithm,
        rebuild,
        cache_model,
    )
    return d
end

"""
$(SIGNATURES)

For `UnrestrictedPerturbationFamily`, the search space is simply [0,1] for each pixel.
The model built can thus be re-used for any input with the same input size.
"""
function model_hash(nn::NeuralNet, input::Array{<:Real}, pp::UnrestrictedPerturbationFamily)::UInt
    input_size = size(input)
    return hash(nn, hash(input_size, hash(pp)))
end

"""
$(SIGNATURES)

For `RestrictedPerturbationFamily`, we take advantage of the restricted input search space
corresponding to each nominal (unperturbed) input by considering only activations to the
non-linear units which are possible for some input in the restricted search space. This
reduces solve times, but also means that the model must be rebuilt for each different
nominal input.
"""
function model_hash(nn::NeuralNet, input::Array{<:Real}, pp::RestrictedPerturbationFamily)::UInt
    return hash(nn, hash(input, hash(pp)))
end

function model_filename(nn::NeuralNet, input::Array{<:Real}, pp::PerturbationFamily)::String
    hash_val = model_hash(nn, input, pp)
    input_size = size(input)
    return "$(nn.UUID).$(input_size).$(string(pp)).$(hash_val).jls"
end

function get_reusable_model(
    nn::NeuralNet,
    input::Array{<:Real},
    pp::PerturbationFamily,
    tightening_solver::MathProgBase.SolverInterface.AbstractMathProgSolver,
    tightening_algorithm::TighteningAlgorithm,
    rebuild::Bool,
    cache_model::Bool,
)::Dict

    filename = model_filename(nn, input, pp)
    model_filepath = joinpath(model_dir, filename)

    if isfile(model_filepath) && !rebuild
        notice(MIPVerify.LOGGER, "Loading model from cache.")
        d = open(model_filepath, "r") do f
            deserialize(f)
            # TODO (vtjeng): Identify situations where the saved model has a different name and throw and error.
        end
        d[:TighteningApproach] = "loaded_from_cache"
        d[:Model].ext[:MIPVerify] = MIPVerifyExt(tightening_algorithm)
    else
        notice(
            MIPVerify.LOGGER,
            """
Rebuilding model from scratch. This may take some time as we determine upper and lower bounds for the input to each non-linear unit.""",
        )
        d = build_reusable_model_uncached(nn, input, pp, tightening_solver, tightening_algorithm)
        if cache_model
            notice(
                MIPVerify.LOGGER,
                """
The model built will be cached and re-used for future solves, unless you explicitly set rebuild=true.""",
            )
            open(model_filepath, "w") do f
                serialize(f, d)
            end
        end
    end
    setsolver(d[:Model], tightening_solver)
    return d
end

function build_reusable_model_uncached(
    nn::NeuralNet,
    input::Array{<:Real},
    pp::UnrestrictedPerturbationFamily,
    tightening_solver::MathProgBase.SolverInterface.AbstractMathProgSolver,
    tightening_algorithm::TighteningAlgorithm,
)::Dict

    m = Model(solver = tightening_solver)
    m.ext[:MIPVerify] = MIPVerifyExt(tightening_algorithm) # TODO: consider writing as seperate function

    input_range = CartesianIndices(size(input))

    # v_input will be constrained to `input` in the caller
    v_input = map(_ -> @variable(m), input_range)
    # v_e is the perturbation added
    v_e = map(_ -> @variable(m, lowerbound = -1, upperbound = 1), input_range)
    # v_x0 is the input with the perturbation added
    v_x0 = map(_ -> @variable(m, lowerbound = 0, upperbound = 1), input_range)
    @constraint(m, v_x0 .== v_input + v_e)

    v_output = v_x0 |> nn

    d = Dict(
        :Model => m,
        :PerturbedInput => v_x0,
        :Perturbation => v_e,
        :Output => v_output,
        :Input => v_input,
        :PerturbationFamily => pp,
        :TighteningApproach => string(tightening_algorithm),
    )

    return d
end

function build_reusable_model_uncached(
    nn::NeuralNet,
    input::Array{<:Real},
    pp::BlurringPerturbationFamily,
    tightening_solver::MathProgBase.SolverInterface.AbstractMathProgSolver,
    tightening_algorithm::TighteningAlgorithm,
)::Dict
    # For blurring perturbations, we build a new model for each input. This enables us to get
    # much better bounds.

    m = Model(solver = tightening_solver)
    m.ext[:MIPVerify] = MIPVerifyExt(tightening_algorithm)

    input_size = size(input)
    num_channels = size(input)[4]
    filter_size = (pp.blur_kernel_size..., num_channels, num_channels)

    v_f = map(_ -> @variable(m, lowerbound = 0, upperbound = 1), CartesianIndices(filter_size))
    @constraint(m, sum(v_f) == num_channels)
    v_x0 = map(_ -> @variable(m, lowerbound = 0, upperbound = 1), CartesianIndices(input_size))
    @constraint(m, v_x0 .== input |> Conv2d(v_f))

    v_output = v_x0 |> nn

    d = Dict(
        :Model => m,
        :PerturbedInput => v_x0,
        :Perturbation => v_x0 - input,
        :Output => v_output,
        :BlurKernel => v_f,
        :PerturbationFamily => pp,
        :TighteningApproach => string(tightening_algorithm),
    )

    return d
end

function build_reusable_model_uncached(
    nn::NeuralNet,
    input::Array{<:Real},
    pp::LInfNormBoundedPerturbationFamily,
    tightening_solver::MathProgBase.SolverInterface.AbstractMathProgSolver,
    tightening_algorithm::TighteningAlgorithm,
)::Dict

    m = Model(solver = tightening_solver)
    m.ext[:MIPVerify] = MIPVerifyExt(tightening_algorithm)

    input_range = CartesianIndices(size(input))
    # v_e is the perturbation added
    v_e =
        map(_ -> @variable(m, lowerbound = -pp.norm_bound, upperbound = pp.norm_bound), input_range)
    # v_x0 is the input with the perturbation added
    v_x0 = map(
        i -> @variable(
            m,
            lowerbound = max(0, input[i] - pp.norm_bound),
            upperbound = min(1, input[i] + pp.norm_bound)
        ),
        input_range,
    )
    @constraint(m, v_x0 .== input + v_e)

    v_output = v_x0 |> nn

    d = Dict(
        :Model => m,
        :PerturbedInput => v_x0,
        :Perturbation => v_e,
        :Output => v_output,
        :PerturbationFamily => pp,
        :TighteningApproach => string(tightening_algorithm),
    )

    return d
end

struct MIPVerifyExt
    tightening_algorithm::TighteningAlgorithm
end


text:
using Test
using MARY_fRG.Basics

@testset "线段相关的内容" begin
    seg1 = Segment(Point2D(1.0, 1.0), Point2D(-2, 7.3))
    println(seg1)
    @test seg1.length == 6.977822009767804
end

ext: jl


text:
##################################################
# Base
##################################################

for op in (:+, :-, :*, :&, :|, :!,
           :length, :size,
           :abs, :abs2, :complex, :conj, :exp, :imag, :inv,
           :log, :one, :oneunit, :real, :sqrt, :transpose, :zero,
           :cos, :sin, :tan, :sec, :csc, :cot,
           :acos, :asin, :atan, :asec, :acsc, :acot,
           :cosh, :sinh, :tanh, :sech, :csch, :coth,
           :acosh, :asinh, :atanh, :asech, :acsch, :acoth,
           :iseven, :isinf, :isnan, :isodd)
    @eval @register $(:(Base.$op)) 1
end
Base.adjoint(x::Sym) = combine(Symbol("'"), x)
Base.zero(::Type{Sym{TAG}}) where TAG = apply(zero, Sym(TAG))
Base.one(::Type{Sym{TAG}}) where TAG = apply(one, Sym(TAG))
Base.oneunit(::Type{Sym{TAG}}) where TAG = apply(oneunit, Sym(TAG))

# Two-arg operators
for op in (:-, :/, :\, ://, :^, :÷, :isless, :<, :&, :|, :(==))
    @eval @register $(:(Base.$op)) 2
end
Base.getindex(x::Sym, val::Symbol) = Base.getindex(promote(x, QuoteNode(val))...)
Base.getindex(x::Sym, vals...) = Base.getindex(promote(x, map(val -> val isa Symbol ? QuoteNode(val) : val, vals)...)...)
Base.getindex(x::Sym, val::Sym) = combine(:ref, x, val)
Base.getindex(x::Sym, vals::Sym...) = combine(:ref, x, vals...)

Base.getproperty(x::Sym, val::Symbol) = Base.getproperty(promote(x, QuoteNode(val))...)
Base.getproperty(x::Sym, val::Sym{Symbol}) = combine(:(.), x, val)

Base.in(x::Sym, y::Sym) = apply(in, x, y)

@register Base.:(:) --> Sym(:(:)) 2
@register Base.:(:) --> Sym(:(:)) 3

# Multi-arg operators
for op in (:+, :*)
    @eval @register $(:(Base.$op)) 2
    @eval $(:(Base.$op))(xs::Sym...) = apply($(:(Base.$op)), xs...)
end

##################################################
# Linear algebra
##################################################

# One-arg operators
for op in (:norm,)
    @eval @register $(:(LinearAlgebra.$op)) 1
end

# Two-arg operators
for op in (:dot, :cross)
    @eval @register $(:(LinearAlgebra.$op)) 2
end


text:
using AutoDiffSource
using BenchmarkTools
using ReverseDiffSource

rosenbrock(x, y) = sum(100*(y-x.^2).^2 + (1.-x).^2)

@δ rosenbrock(x) = rosenbrock(x[1:length(x)-1], x[2:length(x)])

# verify correctness
@assert checkdiff(rosenbrock, δrosenbrock, randn(3))

# handcrafted derivative
function ∇rosenbrock(x)
    x1 = x[1:end-1]
    x2 = x[2:end]
    ∂x = similar(x)
    ∂x[1:end-1] = - 400*x1.*(x2-x1.^2) - 2*(1-x1)
    ∂x[end] = 0
    ∂x[2:end] += 200*(x2-x1.^2)
    ∂x
end

# verify that it is correct
const x0 = randn(3)
@assert checkgrad(rosenbrock, (x0,), ∇rosenbrock(x0))

rds_rosenbrock = rdiff(rosenbrock, (Vector{Float64}, Vector{Float64}))

# benchmark
const x1 = randn(1_000)
trial_manual = @benchmark (rosenbrock(x1), ∇rosenbrock(x1));
trial_rds = @benchmark ((y, ∂x1, ∂x2) = rds_rosenbrock(x1[1:end-1], x1[2:end]);
                        ∂x = similar(x1);
                        ∂x[1:end-1] = ∂x1;
                        ∂x[end] = 0;
                        ∂x[2:end] += ∂x2;
                        ∂x)
trial_auto = @benchmark ((y, ∇r) = δrosenbrock(x1); ∇r())

# auto should be ~15% faster than handcrafted and 60% faster than reverse diff source
@show trial_manual
@show trial_rds
@show trial_auto

""


text:
"""
`module  JAC.MultiPhotonIonization`
    ... a submodel of JAC that contains all methods for computing multi-photon (single-electron) ionization cross sections.
"""
module MultiPhotonIonization

    using  ..AngularMomentum, ..Basics, ..Continuum, ..InteractionStrength, ..ManyElectron, ..Radial

    """
    `struct  MultiPhotonIonization.Settings`
        ... defines a type for the settings in estimating multi-photon (single-electron) ionization cross sections

        + multipoles              ::Array{EmMultipole}           ... Specifies the multipoles of the radiation field that are to be included.
        + gauges                  ::Array{UseGauge}              ... Specifies the gauges to be included into the computations.
        + photonEnergies          ::Array{Float64,1}             ... List of photon energies.
        + printBefore             ::Bool                         ... True, if all energies and lines are printed before their evaluation.
        + lineSelection           ::LineSelection                ... Specifies the selected levels, if any.

    """
    struct Settings
        multipoles                ::Array{EmMultipole}
        gauges                    ::Array{UseGauge}
        photonEnergies            ::Array{Float64,1}
        printBefore               ::Bool
        lineSelection             ::LineSelection
    end

    """
    `MultiPhotonIonization.Settings()`  ... constructor for the default values of multi-photon (single-electron) ionization estimates.
    """
    function Settings()
        Settings(EmMultipole[], UseGauge[], Float64[], false, LineSelection() )
    end

    # `Base.show(io::IO, settings::MultiPhotonIonization.Settings)`
    # 	... prepares a proper printout of the variable settings::MultiPhotonIonization.Settings.
    function Base.show(io::IO, settings::MultiPhotonIonization.Settings)
        println(io, "multipoles:               $(settings.multipoles)  ")
        println(io, "gauges:                   $(settings.gauges)  ")
        println(io, "photonEnergies:           $(settings.photonEnergies)  ")
        println(io, "printBefore:              $(settings.printBefore)  ")
        println(io, "lineSelection:            $(settings.lineSelection)  ")
    end

    """
    `struct  MultiPhotonIonization.Channel`
        ... defines a type for a multi-photon channel to help characterize multi-photon (single-electron)ionization with
            well-defined multipolarities.

        + NoPhotons      ::EmMultipole            ... Number of photons in the multi-photon ionization
        + multipoles     ::Array{EmMultipole,1}   ... Multipoles of all N incoming photons.
        + gauge          ::EmGauge                ... Gauge for dealing with the (coupled) radiation field.
        + kappa          ::Int64                  ... partial-wave of the free electron
        + symmetry       ::LevelSymmetry          ... total angular momentum and parity of the scattering state
        + phase          ::Float64                ... phase of the partial wave
        + amplitude      ::Complex{Float64}       ... multi-photon ionization amplitude associated with the given channel.
   """
    struct  Channel
        NoPhotons        ::EmMultipole
        multipoles       ::Array{EmMultipole,1}
        gauge            ::EmGauge
        kappa            ::Int64
        symmetry         ::LevelSymmetry
        phase            ::Float64
        amplitude        ::Complex{Float64}
    end

    """
    `struct  MultiPhotonIonization.Line`  ... defines a type for a multi-photon ionization line that may include the definition of channels.

        + initialLevel     ::Level                  ... initial-(state) level
        + finalLevel       ::Level                  ... final-(state) level
        + photonEnergy     ::Float64                ... Energy of the incoming photons; all photons are assumed to have equal energy.
        + crossSection     ::EmProperty             ... Cross section for this multi-photon ionization.
        + hasChannels      ::Bool                   ... Determines whether the individual (sub-) channels are defined in terms of their
                                                        multipolarities, etc., or not.
        + channels         ::Array{MultiPhotonIonization.Channel,1}  ... List of MultiPhotonIonization.Channels of this line.
    """
    struct  Line
        initialLevel       ::Level
        finalLevel         ::Level
        photonEnergy       ::Float64
        crossSection       ::EmProperty
        hasChannels        ::Bool
        channels           ::Array{MultiPhotonIonization.Channel,1}
    end

    # `Base.show(io::IO, line::MultiPhotonIonization.Line)`  ... prepares a proper printout of the variable line::MultiPhotonIonization.Line.
    function Base.show(io::IO, line::MultiPhotonIonization.Line)
        println(io, "initialLevel:      $(line.initialLevel)  ")
        println(io, "finalLevel:        $(line.finalLevel)  ")
        println(io, "photonEnergy:      $(line.photonEnergy)  ")
        println(io, "crossSection:      $(line.crossSection)  ")
        println(io, "hasChannels:       $(line.hasChannels)  ")
        println(io, "channels:          $(line.channels)  ")
    end

    """
    `MultiPhotonIonization.computeLines(finalMultiplet::Multiplet, initialMultiplet::Multiplet, grid::Radial.Grid,
                                              settings::MultiPhotonIonization.Settings; output=true)`
        ... to compute the multiphoton transition amplitudes and all properties as requested by the given settings.
            A list of lines::Array{MultiPhotonIonization.Lines} is returned.
    """
    function  computeLines(finalMultiplet::Multiplet, initialMultiplet::Multiplet, grid::Radial.Grid,
                           settings::MultiPhotonIonization.Settings; output=true)
        println("")
        printstyled("MultiPhotonIonization.computeLines(): The computation of multi-photon ionization amplitudes starts now ... \n", color=:light_green)
        printstyled("---------------------------------------------------- ----------------------------------------------------- \n", color=:light_green)
        println("")
        #
        #
        pathways = "Not yet implemented !"
        #
        if    output    return( pathways )
        else            return( nothing )
        end
    end

    # Inlcude one-electron code
    include("module-MultiPhotonIonization-inc-one-electron.jl")

end # module


text:
using DataStructures
using Base.Test

# construction

@test typeof(OrderedSet()) == OrderedSet{Any}
@test typeof(OrderedSet('a')) == OrderedSet{Char}
@test typeof(OrderedSet(1,2,3,4)) == OrderedSet{Int}

# empty set
d = OrderedSet{Char}()
@test length(d) == 0
@test isempty(d)
@test !('c' in d)
push!(d, 'c')
@test !isempty(d)
empty!(d)
@test isempty(d)

# access, modification

for c in 'a':'z'
    push!(d, c)
end

for c in 'a':'z'
    @test c in d
end

@test collect(d) == ['a':'z']


text:
struct MemoryInfo
    n_qubits::Int
    # NOTE:
    # we alloate memory_slots for all measurement
    # at level of codegen since measure cannot be
    # easily removed it can affect the whole state
    memory_offset::Int # init address of classical memory
    returns::Set{Int}
    # allocate slots for returned measure result
    slots::Dict{Int, Vector{Int}}
    # alloate for used measure result
    registers::Dict{Int,Vector{Int}} # ssa value => list of slots
    # allocate qubits on contiguous address
    qubits::Dict{Int,Int} # location => qobj address
end

function MemoryInfo(ci::CodeInfo; offset::Int=0)
    memory_slots = 0
    qubits = Dict{Int, Int}()
    qreg_head = 0
    slots = Dict{Int, Vector{Int}}()
    registers = Dict{Int, Vector{Int}}()
    measure_ssas = measure_ssa(ci)
    measure_uses = measure_ssa_uses(ci, measure_ssas)
    measure_rets = measure_ssa_returns(ci, measure_ssas)
    qmem_ptr = Ref{Int}(0)

    for (v, stmt) in enumerate(ci.code)
        @switch stmt begin
            @case Expr(:invoke, _, GlobalRef(Intrinsics, :apply), _, _, QuoteNode(locs::Locations))
                allocate_qobj_qubits!(qubits, locs, qmem_ptr)
            @case Expr(:invoke, _, GlobalRef(Intrinsics, :apply), _, _, QuoteNode(locs::Locations), QuoteNode(ctrl::CtrlLocations))
                allocate_qobj_qubits!(qubits, locs, qmem_ptr)
                allocate_qobj_qubits!(qubits, ctrl.storage, qmem_ptr)
            @case Expr(:invoke, _, GlobalRef(Intrinsics, :measure), _, QuoteNode(locs::Locations))
                allocate_qobj_qubits!(qubits, locs, qmem_ptr)
                # NOTE: we should always allocate slot for measure
                slots[v] = collect(memory_slots:memory_slots+length(locs)-1)
                memory_slots += length(locs)

                if v in measure_uses
                    registers[v] = collect(qreg_head:qreg_head+length(locs)-1)
                    qreg_head += length(locs)
                end
            # skip other statement
            @case _
                nothing
        end
    end
    MemoryInfo(qmem_ptr[], offset, measure_rets, slots, registers, qubits)
end

function allocate_qobj_qubits!(qubits::Dict{Int, Int}, locs::Locations, ptr::Ref{Int})
    for loc in locs
        plain_loc = plain(loc)
        haskey(qubits, plain_loc) && continue
        qubits[plain_loc] = ptr[]
        ptr[] += 1
    end
    return qubits
end

# pretty print in multiline mode
function Base.show(io::IO, ::MIME"text/plain", mi::MemoryInfo)
    summary(io, mi)
    println(io, "(")
    println(io, "  n_qubits=", mi.n_qubits, ",")
    println(io, "  memory_offset=", mi.memory_offset, ",")
    if !isempty(mi.slots)
        println(io, "  slots=Dict{Int, Vector{Int}}(")
        for (v, slots) in mi.slots
            println(io, "    ", v, "=>", slots, ",")
        end
        println(io, "  ),")
    end

    if !isempty(mi.registers)
        println(io, "  registers=Dict{Int, Vector{Int}}(")
        for (v, slots) in mi.registers
            println(io, "    ", v, "=>", slots, ",")
        end
        println(io, "  ),")
    end

    println(io, "  qubits=Dict{Int,Int}(", )
    for (loc, addr) in mi.qubits
        println(io, "    ", loc, "=>", addr, ",")
    end
    println(io, "  ),")
    println(io, ")")
end

mutable struct CodeGenState
    src::CodeInfo
    pc::Int
    memory_ptr::Int
    stmt::Any
end

CodeGenState(src) = CodeGenState(src, 0, 0, nothing)

function YaoCompiler.compile(target::QobjQASMTarget, f, tt::Type, options::HardwareFreeOptions)
    interp = YaoInterpreter(;target, options)
    mi = method_instance(f, tt)
    src = Core.Compiler.typeinf_ext_toplevel(interp, mi)
    # we need to remove code coverage effect to pass validate
    src = CompilerPluginTools.rm_code_coverage_effect(src)
    validate(target, src)

    st = CodeGenState(src)
    mi = MemoryInfo(src; offset=target.memory_offset)
    return emit_qobj(target, src, mi, st)
end

function emit_qobj(target::QobjQASMTarget, ci::CodeInfo, mi::MemoryInfo, st::CodeGenState)
    instructions = Instruction[]
    cfg = Core.Compiler.compute_basic_blocks(ci.code)
    worklist = Int[1]

    while !isempty(worklist)
        bb = pop!(worklist)
        bb_insts = []
        stmts = cfg.blocks[bb].stmts
        terminator = ci.code[last(stmts)]
        # this flag removes the measure stmt for reset
        contains_reset = false

        if !(terminator isa GotoIfNot)
            # NOTE:
            # terminator can be a gate call in program
            # that contains no control flow so let's try
            # to convert the terminator to Qobj
            st.pc = last(stmts)
            push_stmt!(bb_insts, emit_stmt(target, terminator, mi, st))
        else
            cond = terminator.cond::SSAValue

            @switch ci.code[cond.id] begin
                @case Expr(:invoke, _, GlobalRef(YaoCompiler, :measure_cmp), result::SSAValue, val::Int) ||
                    Expr(:invoke, _, GlobalRef(YaoCompiler, :measure_cmp), val::Int, result::SSAValue)
                    result = result
                    val = val
                @case _
                    error("if condition must be a measurement result compare with a constant integer")
            end

            condition = mi.registers[result.id]
            length(condition) == 1 || error("Qobj only allows single bit condition for if statement")
            # conditional bb
            non_conditional_branch = Core.Compiler.block_for_inst(cfg.index, terminator.dest)
            conditional_bb = findfirst(!isequal(non_conditional_branch), cfg.blocks[bb].succs)
            push!(worklist, non_conditional_branch)
            conditional_stmts = cfg.blocks[conditional_bb].stmts

            # handles multi-stmt if measure_cmp(a, b) ... end
            for v in conditional_stmts
                st.pc = v
                stmt = emit_stmt(target, ci.code[v], mi, st)
                stmt isa Gate || stmt isa Vector || error("cannot use classical control on $stmt in Qobj(type=QASM)")
                # match reset, if measure value is in the return, we will not use reset
                if !(result.id in mi.returns) && length(conditional_stmts) == 1 && is_reset_stmt(ci.code[result.id], stmt, val, mi)
                    push_stmt!(bb_insts, Reset(;qubits=stmt[1].qubits))
                    contains_reset = true
                    continue
                end
                push_stmt!(bb_insts, condition_stmt(stmt, condition[]))
            end
        end

        for v in stmts[1:end-1]
            contains_reset && v == result.id && continue
            st.pc = v
            push_stmt!(instructions, emit_stmt(target, ci.code[v], mi, st))
        end
        append!(instructions, bb_insts)
    end

    memory_slots = 0
    for inst in instructions
        if inst isa Measure
            memory_slots += length(inst.memory)
        end
    end

    return Experiment(;
        header=Dict{String, Any}(
            "n_qubits"=>mi.n_qubits,
            "description"=>target.description
        ),
        config=ExpConfig(;
            shots=target.nshots,
            seed=target.seed,
            max_credits=target.max_credits,
            memory_slots,
        ),
        instructions,
    )
end

# # return Qobj(;
# # qobj_id=target.qobj_id,
# # type="QASM",
# # schema_version=target.schema_version,
# # header=Dict{String, Any}("description"=>string(src.parent)),
# # config=ExpConfig(;
# #     shots=target.shots,
# #     seed=target.seed,
# #     max_credits=target.max_credits,
# #     memory_slots=st.memory_slots,
# # ),
# # experiments=Experiment[],
# # )

function emit_stmt(target::QobjQASMTarget, stmt, mi::MemoryInfo, st::CodeGenState)
    @switch stmt begin
        @case Expr(:invoke, _, GlobalRef(Intrinsics, :apply), _, gate, QuoteNode(locs))
            emit_gate(target, gate, locs, mi, st)
        @case Expr(:invoke, _, GlobalRef(Intrinsics, :apply), _, gate, QuoteNode(locs), QuoteNode(ctrl))
            emit_ctrl(target, gate, locs, ctrl, mi, st)
        @case Expr(:invoke, _, GlobalRef(Intrinsics, :measure), _, QuoteNode(locs))
            emit_measure(target, locs, mi, st)
        @case Expr(:invoke, _, GlobalRef(Intrinsics, :barrier), _, QuoteNode(locs))
            emit_barrier(target, locs, mi, st)
        @case _ # skip other statement
            nothing
    end
end

function emit_gate(target::QobjQASMTarget, @nospecialize(gate), locs::Locations, mi::MemoryInfo, st::CodeGenState)
    # NOTE: intrinsic instructions only supports single qubit
    # operation, and since we treat single qubit intrinsic
    # on multiple location as an IR semantic we will need to
    # implement it in codegen ourselves
    prog = []
    for loc in locs
        qubits = get_qubit_address(mi, loc)
        inst = @switch gate begin
            # parametric gates
            @case QuoteNode(shift(θ))
                Gate(;name="u1", qubits, params=Any[θ])
            @case QuoteNode(Rx(θ))
                Gate(;name="rx", qubits, params=Any[θ])
            @case QuoteNode(Ry(θ))
                Gate(;name="ry", qubits, params=Any[θ])
            @case QuoteNode(Rz(θ))
                Gate(;name="rz", qubits, params=Any[θ])

            @case QuoteNode(&X)
                Gate(;name="x", qubits)
            @case QuoteNode(&Y)
                Gate(;name="y", qubits)
            @case QuoteNode(&Z)
                Gate(;name="z", qubits)
            @case QuoteNode(&H)
                Gate(;name="h", qubits)
            @case QuoteNode(&T)
                Gate(;name="t", qubits)
            @case QuoteNode(&S)
                Gate(;name="s", qubits)
            @case QuoteNode(val::Operation)
                error("fail to inline operationn $val, please report an issue at YaoCompiler")
            @case SSAValue(id) # parametric intrinsic gate
                error("toplevel Qobj do not support parametric routine")
            @case _
                error("unknown gate value: $gate")
        end
        push!(prog, inst)
    end
    return prog
end

function emit_ctrl(::QobjQASMTarget, gate, locs::Locations, ctrl::CtrlLocations, mi::MemoryInfo, st::CodeGenState)
    qubits = get_qubit_address(mi, locs, ctrl)
    @switch (gate, length(ctrl), length(locs)) begin
        @case (QuoteNode(&X), 1, 1)
            Gate(;name="cx", qubits)
        @case (QuoteNode(&Y), 1, 1)
            Gate(;name="cy", qubits)
        @case (QuoteNode(&Z), 1, 1)
            Gate(;name="cz", qubits)
        @case (QuoteNode(&H), 1, 1)
            Gate(;name="ch", qubits)
        @case (QuoteNode(&X), 2, 1)
            Gate(;name="ccx", qubits)

        @case (QuoteNode(Rx(θ)), 1, 1)
            Gate(;name="crx", qubits, params=Any[θ])
        @case (QuoteNode(Ry(θ)), 1, 1)
            Gate(;name="cry", qubits, params=Any[θ])
        @case (QuoteNode(Rz(θ)), 1, 1)
            Gate(;name="crz", qubits, params=Any[θ])
        @case (QuoteNode(shift(θ)), 1, 1)
            Gate(;name="cu1", qubits, params=Any[θ])
        @case (QuoteNode(Intrinsics.UGate(α, β, γ)), 1, 1)
            Gate(;name="cu3", qubits, params=Any[α, β, γ])
        @case _
            error("control statement is not compatible with Qobj(type=QASM)")
    end
end

function emit_measure(::QobjQASMTarget, locs::Locations, mi::MemoryInfo, st::CodeGenState)
    ret = Measure(;
        qubits=get_qubit_address(mi, locs),
        memory=mi.slots[st.pc],
        # check if the measure result is used
        # if it is used, store it in register
        register=haskey(mi.registers, st.pc) ? mi.registers[st.pc] : nothing
    )
    return ret
end

function emit_barrier(target::QobjQASMTarget, locs::Locations, mi::MemoryInfo, st::CodeGenState)
    return Barrier(;qubits=get_qubit_address(mi, locs))
end

function condition_stmt(stmt::Vector, conditional)
    [condition_stmt(each, conditional) for each in stmt]
end

function condition_stmt(stmt::Gate, conditional)
    return Gate(;
        name=stmt.name,
        qubits=stmt.qubits,
        params=stmt.params,
        texparams=stmt.texparams,
        conditional,
    )
end

function is_reset_stmt(measure, stmt, val, mi::MemoryInfo)
    stmt isa Vector && length(stmt) == 1 || return false
    stmt[1] isa Gate && stmt[1].name == "x" && val == 1 || return false

    # check if measure loc is the same as X gate loc
    @switch measure begin
        @case Expr(:invoke, _, GlobalRef(Intrinsics, :measure), _, QuoteNode(locs::Locations))
            return get_qubit_address(mi, locs) == stmt[1].qubits
        @case _
            return false
    end
end

function get_qubit_address(mi::MemoryInfo, locs::Locations)
    qubits = Vector{Int}(undef, length(locs))
    get_qubit_address!(qubits, mi, locs)
    return qubits
end

function get_qubit_address(mi::MemoryInfo, locs::Locations, ctrl::CtrlLocations)
    qubits = Vector{Int}(undef, length(locs) + length(ctrl))
    get_qubit_address!(qubits, mi, locs)
    get_qubit_address!(qubits, mi, ctrl.storage, length(locs))
    return qubits
end

function get_qubit_address!(qubits::Vector{Int}, mi::MemoryInfo, locs::Locations, offset::Int=0)
    for (k, loc) in enumerate(locs)
        qubits[k + offset] = mi.qubits[plain(loc)]
    end
    return qubits
end

function measure_ssa(ci::CodeInfo)
    measure_ssa = Int[]
    for (v, stmt) in enumerate(ci.code)
        @switch stmt begin
            @case Expr(:invoke, _, GlobalRef(Intrinsics, :measure), _, QuoteNode(locs))
                push!(measure_ssa, v)
            @case _
                nothing
        end
    end
    return measure_ssa
end

measure_ssa_uses(ci::CodeInfo, measure_ssa=measure_ssa(ci)) =
    measure_ssa_uses!(Set{Int}(), ci, measure_ssa)

function measure_ssa_uses!(used::Set{Int}, ci::CodeInfo, measure_ssa=measure_ssa(ci))
    for (v, stmt) in enumerate(ci.code)
        @switch stmt begin
            @case GotoIfNot(SSAValue(id), dest)
                cond_stmt = ci.code[id]
                for useref in Core.Compiler.userefs(cond_stmt)
                    val = useref[]
                    if isa(val, SSAValue) && val.id in measure_ssa
                        push!(used, val.id)
                    end
                end
            @case _
                nothing
        end
    end
    return used
end

function measure_ssa_returns(ci::CodeInfo, measure_ssa=measure_ssa(ci))
    measure_ssa_returns!(Set{Int}(), ci, measure_ssa)
end

function measure_ssa_returns!(ret::Set{Int}, ci::CodeInfo, measure_ssa=measure_ssa(ci))
    for (v, stmt) in enumerate(ci.code)
        stmt isa ReturnNode || continue
        find_return_ssa_deps!(ret, ci, stmt, measure_ssa)
    end
    return ret
end

function find_return_ssa_deps!(ret::Set, ci::CodeInfo, stmt, measure_ssa)
    for useref in Core.Compiler.userefs(stmt)
        val = useref[]
        isa(val, SSAValue) || continue
        if val.id in measure_ssa
            push!(ret, val.id)
        else
            find_return_ssa_deps!(ret, ci, ci.code[val.id], measure_ssa)
        end
    end
    return ret
end


text:
module Timeout
export timeout, TimeoutException

struct TimeoutException <: Exception end

function killtask(task)
    try
        Base.throwto(task, InterruptException())
    catch
    end
end

"""
    timeout(f, secs=3) -> Function

Return an anonymous function that calls function `f`. If the execution time of `f` exceeds `secs`, then throw the `TimeoutException`.
"""
function timeout(f, secs=3)
    function fn()
        done = Channel{Symbol}(1)
        result = Channel(1)

        timer_task = @async begin
            sleep(secs)
            put!(done, :timeout)
            put!(result, nothing)
        end

        main_task = @async begin
            put!(result, f())
            put!(done, :ok)
            killtask(timer_task)
        end

        state = fetch(done)
        r = fetch(result)
        close(done); close(result)
        if state == :timeout
            killtask(main_task)
            throw(TimeoutException())
        end

        return r
    end
end

end # module


text:
# Optimization and solution query

function optimize(model::Model)
    @assert model.ptr_model != C_NULL
    ret = @grb_ccall(optimize, Cint, (Ptr{Void},), model)
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    nothing
end

function computeIIS(model::Model)
    @assert model.ptr_model != C_NULL
    ret = @grb_ccall(computeIIS, Cint, (Ptr{Void},), model)
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    nothing
end

#################################################
#
#  solution status and optimization info
#
#################################################

const GRB_LOADED          = 1
const GRB_OPTIMAL         = 2
const GRB_INFEASIBLE      = 3
const GRB_INF_OR_UNBD     = 4
const GRB_BOUNDED         = 5
const GRB_CUTOFF          = 6
const GRB_ITERATION_LIMIT = 7
const GRB_NODE_LIMIT      = 8
const GRB_TIME_LIMIT      = 9
const GRB_SOLUTION_LIMIT  = 10
const GRB_INTERRUPTED     = 11
const GRB_NUMERIC         = 12
const GRB_SUBOPTIMAL      = 13

const status_symbols = [
    :loaded,
    :optimal,
    :infeasible,
    :inf_or_unbd,
    :unbounded,
    :cutoff,
    :iteration_limit,
    :node_limit,
    :time_limit,
    :solution_limit,
    :interrupted,
    :numeric,
    :suboptimal
]

get_status_code(model::Model) = get_intattr(model, "Status")
get_status(model::Model) = status_symbols[get_status_code(model)]::Symbol

@grb_dbl_attr get_runtime      "Runtime"
@grb_dbl_attr get_objval       "ObjVal"
@grb_dbl_attr get_objbound     "ObjBound"

@grb_int_attr get_sol_count    "SolCount"
@grb_int_attr get_barrier_iter_count "BarIterCount"
#@grb_dbl_attr get_node_count   "NodeCount"

get_iter_count(model::Model) = convert(Int, get_dblattr(model, "IterCount"))
get_node_count(model::Model) = convert(Int, get_dblattr(model, "NodeCount"))

type OptimInfo
    status::Symbol
    runtime::Float64

    sol_count::Int
    iter_count::Int
    barrier_iter_count::Int
    node_count::Int
end

function get_optiminfo(model::Model)
    OptimInfo(
        get_status(model),
        get_runtime(model),

        get_sol_count(model),
        get_iter_count(model),
        get_barrier_iter_count(model),
        get_node_count(model)
    )
end

function show(io::IO, s::OptimInfo)
    println(io, "Gurobi Optimization Info")
    println(io, "    status   = $(s.status)")
    println(io, "    runtime  = $(s.runtime)")
    println(io, "    # solutions = $(s.sol_count)")
    println(io, "    # iters     = $(s.iter_count)")
    println(io, "    # bar iters = $(s.barrier_iter_count)")
    println(io, "    # nodes     = $(s.node_count)")
end

#################################################
#
#  solution query
#
#################################################

get_solution(model::Model) = get_dblattrarray(model, "X", 1, num_vars(model))

const varmap = Dict(
    -3 => :Superbasic,
    -2 => :NonbasicAtUpper,
    -1 => :NonbasicAtLower,
     0 => :Basic
)

const conmap = Dict(
     0 => :Basic,
    -1 => :Nonbasic)

function get_basis(model::Model)
    cval = Array(Cint, num_vars(model))
    cbasis = Array(Symbol, num_vars(model))
    get_intattrarray!(cval, model, "VBasis", 1)
    for it in 1:length(cval)
        cbasis[it] = varmap[cval[it]]
    end

    rval = Array(Cint, num_constrs(model))
    rbasis = Array(Symbol, num_constrs(model))
    get_intattrarray!(rval, model, "CBasis", 1)
    rsense = Array(Cchar, num_constrs(model))
    get_charattrarray!(rsense, model, "Sense", 1)
    for it in 1:length(rval)
        rbasis[it] = conmap[rval[it]]
        if rbasis[it] == :Nonbasic
            if rsense[it] == convert(Cchar,'<')
                rbasis[it] = :NonbasicAtUpper
            else
                rbasis[it] = :NonbasicAtLower
            end
        end
    end
    return cbasis, rbasis
end


text:
using OmegaCore, OmegaCore.Proposals
using OmegaTest
using Distributions
using Test

function test_logpdf_ucond()
  a = 1 ~ Normal(0, 1)
  b = 2 ~ Normal(2, 3)
  z(ω) = (a(ω), b(ω))
  ω, ℓ = propose(Random.GLOBAL_RNG, z)
  @test isinferred(propose, Random.GLOBAL_RNG, z)
  a_, b_ = z(ω)
  @test ℓ == logpdf(Normal(0, 1), a_) + logpdf(Normal(2, 3), b_)
end

# Test doesn't double count
function test_logpdf_ucond_nodouble()
  a = 1 ~ Normal(0, 1)
  b = 2 ~ Normal(2, 3)
  z(ω) = (a(ω), b(ω), a(ω))
  ω, ℓ = propose(Random.GLOBAL_RNG, z)
  @test isinferred(propose, Random.GLOBAL_RNG, z)
  a_, b_ = z(ω)
  @test ℓ == logpdf(Normal(0, 1), a_) + logpdf(Normal(2, 3), b_)
end

function test_logpdf_cond()
  μ = 1 ~ Normal(0, 1)
  x = 2 ~ Normalₚ(μ, 1)
  μ_ = 0.1234
  x_ = -0.54321
  μc = μ |ᶜ x == x_
  ω, ℓ = propose(Random.GLOBAL_RNG, μc)
end

@testset "Logpdf" begin
  test_logpdf_ucond()
  test_logpdf_ucond_nodouble()
end


text:
function JB2008(arg0::JB2008InputParameters, arg1::PVCoordinatesProvider, arg2::BodyShape)
    return JB2008((JB2008InputParameters, PVCoordinatesProvider, BodyShape), arg0, arg1, arg2)
end

function JB2008(arg0::JB2008InputParameters, arg1::PVCoordinatesProvider, arg2::BodyShape, arg3::TimeScale)
    return JB2008((JB2008InputParameters, PVCoordinatesProvider, BodyShape, TimeScale), arg0, arg1, arg2, arg3)
end

function get_density(obj::JB2008, arg0::AbsoluteDate, arg1::Vector3D, arg2::Frame)
    return jcall(obj, "getDensity", jdouble, (AbsoluteDate, Vector3D, Frame), arg0, arg1, arg2)
end

function get_density(obj::JB2008, arg0::FieldAbsoluteDate, arg1::FieldVector3D, arg2::Frame)
    return jcall(obj, "getDensity", RealFieldElement, (FieldAbsoluteDate, FieldVector3D, Frame), arg0, arg1, arg2)
end

function get_density(obj::JB2008, arg0::RealFieldElement, arg1::RealFieldElement, arg10::jdouble, arg11::jdouble, arg12::jdouble, arg13::jdouble, arg14::jdouble, arg2::RealFieldElement, arg3::RealFieldElement, arg4::RealFieldElement, arg5::RealFieldElement, arg6::jdouble, arg7::jdouble, arg8::jdouble, arg9::jdouble)
    return jcall(obj, "getDensity", RealFieldElement, (RealFieldElement, RealFieldElement, jdouble, jdouble, jdouble, jdouble, jdouble, RealFieldElement, RealFieldElement, RealFieldElement, RealFieldElement, jdouble, jdouble, jdouble, jdouble), arg0, arg1, arg10, arg11, arg12, arg13, arg14, arg2, arg3, arg4, arg5, arg6, arg7, arg8, arg9)
end

function get_density(obj::JB2008, arg0::jdouble, arg1::jdouble, arg10::jdouble, arg11::jdouble, arg12::jdouble, arg13::jdouble, arg14::jdouble, arg2::jdouble, arg3::jdouble, arg4::jdouble, arg5::jdouble, arg6::jdouble, arg7::jdouble, arg8::jdouble, arg9::jdouble)
    return jcall(obj, "getDensity", jdouble, (jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble), arg0, arg1, arg10, arg11, arg12, arg13, arg14, arg2, arg3, arg4, arg5, arg6, arg7, arg8, arg9)
end

function get_frame(obj::JB2008)
    return jcall(obj, "getFrame", Frame, ())
end


text:
FitFunction(s::Symbol) = FitFunction(Float64, Val(s))
FitFunction(T::Type{<:AbstractFloat}, s::Symbol) = FitFunction(T, Val(s))

function FitFunction(T::Type{<:AbstractFloat}, s::Val{:Gauss})
    ff = FitFunction{T}( Gauss, 1, 3 )
    set_parameter_names!(ff, ["A", "σ", "μ"])
    set_initial_parameters!(ff, [1, 1, 0])
    set_fitranges!(ff, ([-1, 1],))
    return ff
end
function FitFunction(T::Type{<:AbstractFloat}, s::Val{:PDF_Gauss})
    ff = FitFunction{T}( PDF_Gauss, 1, 2 )
    set_parameter_names!(ff, ["σ", "μ"])
    set_initial_parameters!(ff, [1, 0])
    set_fitranges!(ff, ([-1, 1],))
    return ff
end

function FitFunction(T::Type{<:AbstractFloat}, s::Val{:GaussPlusConstantBackground})
    ff = FitFunction{T}( Gauss_plus_const_background, 1, 4 )
    set_parameter_names!(ff, ["A", "σ", "μ", "offset"])
    set_initial_parameters!(ff, [1, 1, 0, 1])
    set_fitranges!(ff, ([-1, 1],))
    return ff
end
function FitFunction(T::Type{<:AbstractFloat}, s::Val{:GaussPlusLinearBackground})
    ff = FitFunction{T}( Gauss_plus_linear_background, 1, 5 )
    set_parameter_names!(ff, ["A", "σ", "μ", "offset", "lin. slope"])
    set_initial_parameters!(ff, [1, 1, 0, 1, 0])
    set_fitranges!(ff, ([-1, 1],))
    return ff
end

function Gauss(x::T, A::T, σ::T, μ::T)::T where {T <: AbstractFloat}
    r::T = A * exp( -(x - μ)^2 / (2 * σ^2)) / sqrt(2π * σ^2)
    if σ <= 0 || A < 0 return T(Inf) end
    if r < 0 r = 0 end
    return r
end
function Gauss(x::Vector{T}, A::T, σ::T, μ::T)::Vector{T} where {T <: AbstractFloat}
    return Gauss.( x, A, σ, μ)
end
Gauss(x::Union{T, Vector{T}}, p::Vector{T}) where {T <: AbstractFloat} = Gauss(x, p... )

function PDF_Gauss(x::T, σ::T, μ::T)::T where {T <: AbstractFloat}
    r::T = exp( -(x - μ)^2 / (2 * σ^2)) / sqrt(2π * σ^2)
    if σ <= 0 return T(Inf) end
    if r < 0 r = 0 end
    return r
end
function PDF_Gauss(x::Vector{T}, σ::T, μ::T)::Vector{T} where {T <: AbstractFloat}
    return PDF_Gauss.( x, σ, μ)
end
PDF_Gauss(x::Union{T, Vector{T}}, p::Vector{T}) where {T <: AbstractFloat} = PDF_Gauss(x, p... )

function Gauss_plus_linear_background(x::T, A::T, σ::T, μ::T, offset::T, lin_slope::T)::T where {T <: AbstractFloat}
    r::T = A * exp( -(x - μ)^2 / (2 * σ^2)) / sqrt(2π * σ^2) + lin_slope * (x - μ) + offset
    if σ <= 0 || A < 0 || offset < 0 return T(Inf) end
    if r < 0 r = 0 end
    return r
end
function Gauss_plus_linear_background(x::Vector{T}, A::T, σ::T, μ::T, offset::T, lin_slope::T)::Vector{T} where {T <: AbstractFloat}
    return Gauss_plus_linear_background.( x, A, σ, μ, offset, lin_slope)
end
Gauss_plus_linear_background(x::Union{T, Vector{T}}, p::Vector{T}) where {T <: AbstractFloat} = Gauss_plus_linear_background(x, p... )

function Gauss_plus_const_background(x::T, A::T, σ::T, μ::T, offset::T)::T where {T <: AbstractFloat}
    r::T = A * exp( -(x - μ)^2 / (2 * σ^2)) / sqrt(2π * σ^2) + offset
    if σ <= 0 || A < 0 || offset < 0 return T(Inf) end
    if r < 0 r = 0 end
    return r
end
function Gauss_plus_const_background(x::Vector{T}, A::T, σ::T, μ::T, offset::T)::Vector{T} where {T <: AbstractFloat}
    return Gauss_plus_const_background.( x, A, σ, μ, offset)
    if σ <= 0 || A < 0 return T(Inf) end
end
Gauss_plus_const_background(x::Union{T, Vector{T}}, p::Vector{T}) where {T <: AbstractFloat} = Gauss_plus_const_background(x, p... )


text:
@testset "evaluation error handling" begin

using Weave: unwrap_load_err

function get_err_str(str::AbstractString)
    try
        include_string(Main, str)
    catch _err
        err = unwrap_load_err(_err)
        return sprint(showerror, err)
    end
end

err_stmt1 = "using NonExisting"
err_stmt2 = "x = "
err_stmt3 = """
plot(x)
y = 10
f(y
"""

str = """
```julia
$err_stmt1
```

```julia
$err_stmt2
```

```julia; term = true
$err_stmt3
```
"""

err_str1 = get_err_str(err_stmt1)
err_str2 = get_err_str(err_stmt2)
err_str3_1 = get_err_str("plot(x)")
err_str3_2 = get_err_str("f(y")

let doc = mock_run(str; doctype = "github")
    get_output(i) = doc.chunks[i].output

    @test occursin(err_str1, get_output(1))
    @test occursin(err_str2, get_output(2))
    @test occursin(err_str3_1, get_output(3))
    @test occursin(err_str3_2, get_output(3))
end

# TODO: move this into chunk option tests
str = """
```julia; error = true
using  # won't be thrown
```

```julia; error = false
using NonExisting # will be thrown
```
"""

@test_throws ArgumentError mock_run(str; doctype = "github")

end  # @testset "evaluation error handling"


text:
using BenchmarkTools
using Pkg, FinancialFFT;
#get(ENV, "CI", nothing) == "true" ? Pkg.instantiate() : nothing;
path1 = joinpath(dirname(pathof(FinancialFFT)), "..", "benchmark")
test_listTmp = readdir(path1);
BlackList = ["Project.toml", "benchmarks.jl", "runner.jl", "cuda", "af", "Manifest.toml", "bench_kou_rev_diff_grad.jl", "bench_black_mp.jl", "bench_black_mn.jl", "bench_black_mt.jl"]
test_list = [test_element for test_element in test_listTmp if !Bool(sum(test_element .== BlackList))]
println("Running tests:\n")
for (current_test, i) in zip(test_list, 1:length(test_list))
    println("------------------------------------------------------------")
    println("  * $(current_test) *")
    include(joinpath(path1, current_test))
    println("------------------------------------------------------------")
    if (i < length(test_list))
        println("")
    end
end


text:
function relatedness_no_boot_fix2(data::PopData, sample_names::Vector{String}; method::F, inbreeding::Bool) where F
    loci_names = Symbol.(loci(data))
    n_samples = length(samples(data))
    sample_pairs = pairwise_pairs(sample_names)
    n_per_loci = DataFrames.combine(groupby(data.loci, :locus), :genotype => nonmissing => :n)[:, :n]
    allele_frequencies = allele_freq(data)
    relate_vecs = map(i -> Vector{Union{Missing,Float64}}(undef, length(sample_pairs)), 1:length(method))
    shared_loci = Vector{Int}(undef, length(sample_pairs))
    p = Progress(length(sample_pairs), dt = 1, color = :blue)
    popdata_idx = groupby(data.loci, :name)
    @inbounds Base.Threads.@threads for i in 1:length(sample_pairs)
        @inbounds sample_1 = sample_pairs[i][1]
        @inbounds sample_2 = sample_pairs[i][2]
        @inbounds geno1 = popdata_idx[(sample_1,)].genotype
        @inbounds geno2 = popdata_idx[(sample_2,)].genotype
        keep_idx = nonmissings(geno1, geno2)
        @inbounds shared_loci[i] = length(keep_idx)
        @inbounds [relate_vecs[j][i] = @inbounds mth(geno1[keep_idx], geno2[keep_idx], loci_names[keep_idx], allele_frequencies, loc_n = n_per_loci[keep_idx], n_samples = n_samples, inbreeding = inbreeding) for (j,mth) in enumerate(method)]
        next!(p)
    end
    method_colnames = [Symbol("$i") for i in method]
    out_df = DataFrame(:sample_1 => getindex.(sample_pairs, 1), :sample_2 => getindex.(sample_pairs, 2), :n_loci => shared_loci)
    [out_df[:, mth] = relate_vecs[i] for (i, mth) in enumerate(method_colnames)]
    return out_df
end

function relatedness_boot_nonmissing(data::PopData, sample_names::Vector{String}; method::F, iterations::Int, interval::Tuple{Float64, Float64} = (0.025, 0.975), inbreeding::Bool) where F
    loci_names = Symbol.(loci(data))
    sample_pairs = pairwise_pairs(sample_names)
    n_samples = length(samples(data))
    n_per_loci = DataFrames.combine(groupby(data.loci, :locus), :genotype => nonmissing => :n)[:, :n]
    allele_frequencies = allele_freq(data)
    relate_vecs = map(i -> Vector{Union{Missing,Float64}}(undef, length(sample_pairs)), 1:length(method))
    boot_means, boot_medians, boot_ses = map(i -> deepcopy(relate_vecs), 1:3)
    boot_CI = map(i -> Vector{Union{Missing,Tuple{Float64,Float64}}}(undef, length(sample_pairs)), 1:length(method))
    shared_loci = Vector{Int}(undef, length(sample_pairs))
    p = Progress(length(sample_pairs), dt = 1, color = :blue)
    popdata_idx = groupby(data.loci, :name)
    @inbounds Base.Threads.@threads for i in 1:length(sample_pairs)
        @inbounds sample_1 = sample_pairs[i][1]
        @inbounds sample_2 = sample_pairs[i][2]
        @inbounds geno1 = popdata_idx[(sample_1,)].genotype
        @inbounds geno2 = popdata_idx[(sample_2,)].genotype
        # get index for genotype appearing missing in at least one individual in the pair
        keep_idx = nonmissings(geno1, geno2)
        # generate nonmissing genotype data
        gen1, gen2, loc, n_per_loc = (i[keep_idx] for i in (geno1, geno2, loci_names, n_per_loci))
        @inbounds shared_loci[i] = length(keep_idx)
        @inbounds for (j, mthd) in enumerate(method)
            @inbounds relate_vecs[j][i] = mthd(gen1, gen2, loc, allele_frequencies, loc_n = n_per_loc, n_samples = n_samples, inbreeding = inbreeding)
            boot_out = bootstrap_genos_nonmissing(gen1, gen2, loc, n_per_loc, allele_frequencies, method = mthd, iterations = iterations, inbreeding = inbreeding)
            @inbounds boot_means[j][i], boot_medians[j][i], boot_ses[j][i], boot_CI[j][i] = bootstrap_summary(boot_out, interval)
        end
        next!(p)
    end
    method_colnames = [Symbol("$i") for i in method]
    boot_mean_colnames = [Symbol("$i"*"_mean") for i in method]
    boot_median_colnames = [Symbol("$i"*"_median") for i in method]
    boot_se_colnames = [Symbol("$i"*"_SE") for i in method]
    CI_percent = convert(Int64, round(interval[2] - interval[1], digits = 2) * 100)
    boot_CI_colnames = [Symbol("$i"*"_CI_"*"$CI_percent") for i in method]

    out_df = DataFrame(:sample_1 => map(i -> i[1], sample_pairs), :sample_2 => map(i -> i[2], sample_pairs), :n_loci => shared_loci)
    @inbounds for (i, mth) in enumerate(method_colnames)
        out_df[:, mth] = relate_vecs[i]
        out_df[:, boot_mean_colnames[i]] = boot_means[i]
        out_df[:, boot_median_colnames[i]] = boot_medians[i]
        out_df[:, boot_se_colnames[i]] = boot_ses[i]
        out_df[:, boot_CI_colnames[i]] = boot_CI[i]
    end

    return out_df
end

function relatedness_boot_all(data::PopData, sample_names::Vector{String}; method::F, iterations::Int = 100, interval::Tuple{Float64, Float64} = (0.025, 0.975), inbreeding::Bool) where F
    loci_names = Symbol.(loci(data))
    sample_pairs = pairwise_pairs(sample_names)
    n_samples = length(samples(data))
    allele_frequencies = allele_freq(data)
    n_per_loci = DataFrames.combine(groupby(data.loci, :locus), :genotype => nonmissing => :n)[:, :n]
    relate_vecs = map(i -> Vector{Union{Missing,Float64}}(undef, length(sample_pairs)), 1:length(method))
    boot_means, boot_medians, boot_ses = map(i -> deepcopy(relate_vecs), 1:3)
    boot_CI = map(i -> Vector{Union{Missing,Tuple{Float64,Float64}}}(undef, length(sample_pairs)), 1:length(method))
    shared_loci = Vector{Int}(undef, length(sample_pairs))
    p = Progress(length(sample_pairs), dt = 1, color = :blue)
    popdata_idx = groupby(data.loci, :name)
    @inbounds Base.Threads.@threads for i in 1:length(sample_pairs)
        @inbounds sample_1 = sample_pairs[i][1]
        @inbounds sample_2 = sample_pairs[i][2]
        @inbounds geno1 = popdata_idx[(sample_1,)].genotype
        @inbounds geno2 = popdata_idx[(sample_2,)].genotype
        # get index for genotype appearing missing in at least one individual in the pair
        keep_idx = nonmissings(geno1, geno2)
        # generate nonmissing genotype data
        gen1, gen2, loc, n_per_loc = (i[keep_idx] for i in (geno1, geno2, loci_names, n_per_loci))
        @inbounds shared_loci[i] = length(keep_idx)
        @inbounds for (j, mthd) in enumerate(method)
            @inbounds relate_vecs[j][i] = mthd(gen1, gen2, loc, allele_frequencies, loc_n = n_per_loci, n_samples = n_samples, inbreeding = inbreeding)
            boot_out = bootstrap_genos_all(geno1, geno2, loci_names, n_per_loci, allele_frequencies, method = mthd, iterations = iterations, inbreeding = inbreeding)
            @inbounds boot_means[j][i], boot_medians[j][i], boot_ses[j][i], boot_CI[j][i] = bootstrap_summary(boot_out, interval)
        end
        next!(p)
    end
    method_colnames = [Symbol("$i") for i in method]
    boot_mean_colnames = [Symbol("$i"*"_mean") for i in method]
    boot_median_colnames = [Symbol("$i"*"_median") for i in method]
    boot_se_colnames = [Symbol("$i"*"_SE") for i in method]
    CI_percent = convert(Int64, round(interval[2] - interval[1], digits = 2) * 100)
    boot_CI_colnames = [Symbol("$i"*"_CI_"*"$CI_percent") for i in method]

    out_df = DataFrame(:sample_1 => map(i -> i[1], sample_pairs), :sample_2 => map(i -> i[2], sample_pairs), :n_loci => shared_loci)
    @inbounds for (i, mth) in enumerate(method_colnames)
        out_df[:, mth] = relate_vecs[i]
        out_df[:, boot_mean_colnames[i]] = boot_means[i]
        out_df[:, boot_median_colnames[i]] = boot_medians[i]
        out_df[:, boot_se_colnames[i]] = boot_ses[i]
        out_df[:, boot_CI_colnames[i]] = boot_CI[i]
    end
    return out_df
end

function relatedness_boot_all(data::PopData, sample_names::Vector{String}; method::F, iterations::Int = 100, interval::Tuple{Float64, Float64} = (0.025, 0.975), inbreeding::Bool) where F
    loci_names = Symbol.(loci(data))
    sample_pairs = pairwise_pairs(sample_names)
    n_samples = length(samples(data))
    allele_frequencies = allele_freq(data)
    n_per_loci = DataFrames.combine(groupby(data.loci, :locus), :genotype => nonmissing => :n)[:, :n]
    relate_vecs = map(i -> Vector{Union{Missing,Float64}}(undef, length(sample_pairs)), 1:length(method))
    boot_means, boot_medians, boot_ses = map(i -> deepcopy(relate_vecs), 1:3)
    boot_CI = map(i -> Vector{Union{Missing,Tuple{Float64,Float64}}}(undef, length(sample_pairs)), 1:length(method))
    shared_loci = Vector{Int}(undef, length(sample_pairs))
    p = Progress(length(sample_pairs), dt = 1, color = :blue)
    popdata_idx = groupby(data.loci, :name)
    idx = 0
    @inbounds for (sample_n, ind1) in enumerate(sample_names[1:end-1])
        geno1 = popdata_idx[(ind1,)].genotype
        @inbounds @sync Base.Threads.@spawn for ind2 in sample_names[sample_n+1:end]
            idx += 1
            geno2 = popdata_idx[(ind2,)].genotype

            # get index for genotype appearing missing in at least one individual in the pair
            keep_idx = nonmissings(geno1, geno2)
            # generate nonmissing genotype data
            gen1, gen2, loc, n_per_loc = (i[keep_idx] for i in (geno1, geno2, loci_names, n_per_loci))

            # populate shared_loci array
            @inbounds shared_loci[idx] = length(loc)

            @inbounds for (i, mthd) in enumerate(method)
                @inbounds relate_vecs[i][idx] = mthd(gen1, gen2, loc, allele_frequencies, loc_n = n_per_loci, n_samples = n_samples, inbreeding = inbreeding)
                boot_out = bootstrap_genos_all(geno1, geno2, loci_names, n_per_loci, allele_frequencies, method = mthd, iterations = iterations, inbreeding = inbreeding)
                @inbounds boot_means[i][idx], boot_medians[i][idx], boot_ses[i][idx], boot_CI[i][idx] = bootstrap_summary(boot_out, interval)
            end
            update!(p, idx)
        end
    end
    method_colnames = [Symbol("$i") for i in method]
    boot_mean_colnames = [Symbol("$i"*"_mean") for i in method]
    boot_median_colnames = [Symbol("$i"*"_median") for i in method]
    boot_se_colnames = [Symbol("$i"*"_SE") for i in method]
    CI_percent = convert(Int64, round(interval[2] - interval[1], digits = 2) * 100)
    boot_CI_colnames = [Symbol("$i"*"_CI_"*"$CI_percent") for i in method]

    out_df = DataFrame(:sample_1 => map(i -> i[1], sample_pairs), :sample_2 => map(i -> i[2], sample_pairs), :n_loci => shared_loci)
    @inbounds for (i, mth) in enumerate(method_colnames)
        out_df[:, mth] = relate_vecs[i]
        out_df[:, boot_mean_colnames[i]] = boot_means[i]
        out_df[:, boot_median_colnames[i]] = boot_medians[i]
        out_df[:, boot_se_colnames[i]] = boot_ses[i]
        out_df[:, boot_CI_colnames[i]] = boot_CI[i]
    end

    return out_df
end


text:
export PeriodicCurve

struct PeriodicCurve{S<:Space,T,VT} <: PeriodicDomain{T}
    curve::Fun{S,T,VT}
end

==(a::PeriodicCurve, b::PeriodicCurve) = a.curve == b.curve
isempty(::PeriodicCurve) = false

points(c::PeriodicCurve, n::Integer) = c.curve.(points(domain(c.curve),n))

checkpoints(d::PeriodicCurve) = fromcanonical.(Ref(d),checkpoints(domain(d.curve)))

for op in (:(leftendpoint),:(rightendpoint),:(rand))
    @eval $op(c::PeriodicCurve) = c.curve($op(domain(c.curve)))
end

canonicaldomain(c::PeriodicCurve) = domain(c.curve)

fromcanonical(c::PeriodicCurve{S,T},x) where {S<:Space,T<:Number} = c.curve(x)
function tocanonical(c::PeriodicCurve,x)
    rts=roots(c.curve-x)
    @assert length(rts)==1
    first(rts)
end

fromcanonicalD(c::PeriodicCurve,x)=differentiate(c.curve)(x)

function indomain(x,c::PeriodicCurve)
    rts=roots(c.curve-x)
    if length(rts) ≠ 1
        false
    else
        in(first(rts),canonicaldomain(c))
    end
end

reverseorientation(d::PeriodicCurve) = PeriodicCurve(reverseorientation(d.curve))

isambiguous(d::PeriodicCurve) = ncoefficients(d.curve)==0 && isambiguous(domain(d.curve))

convert(::Type{PeriodicCurve{S,T}},::AnyDomain) where {S,T}=Fun(S(AnyDomain()),[NaN])

arclength(d::PeriodicCurve) = linesum(ones(d))


text:
using SatelliteToolbox
using Plots
gr()

#=
設定パラメータ
=#
DataNum =250 #シミュレータ反復回数
dt = 30 ##シミュレータの計算間隔 [s]
t0 = DatetoJD(2019, 12, 19, 8, 27, 21)   #シミュレート開始時刻＠ユリウス日
TLEFileName = "ISS_TLE.txt"

#事前準備関数定義
function set_SPCal(Name, t0)
    tles = read_tle(Name) #TLE読み込み
    #println("readTLE_Complete")
    orbp = init_orbit_propagator(Val{:sgp4}, tles[1]) #SGP4での軌道モデルの読み込み
    #println("OrbitModelDownload_Complete")
    eop_IAU2000A = get_iers_eop(:IAU2000A)  #IAU-2000Aでの地球モデルの読み込み
    #println("EarthModelDownload_Complete")
    o, r, v = propagate_to_epoch!(orbp, t0 - dt/(24*60*60))
    #println("1stStepOrbitCalcurate_Complete")
    return orbp, eop_IAU2000A, o, r, v
end

#ループ毎処理関数定義
function SPCal(n, t0, dt, orbp, eop_IAU2000A)
    #=  衛星位置の赤道面座標系での計算　=#
    o, r, v = step!(orbp, dt)
    #=  衛星位置の地心直交座標系への変換  =#
    xg = rECItoECEF(GCRF(), ITRF(), t0 + n*dt/(24*60*60), eop_IAU2000A) * r
    vg = rECItoECEF(GCRF(), ITRF(), t0 + n*dt/(24*60*60), eop_IAU2000A) * v
    return xg, vg
end

x_ecef_log = zeros(1000, 3)    #衛星位置＠地心直交座標系
x_geod_log = zeros(1000, 3)    #衛星位置＠測地座標系
v_ecef_log = zeros(1000, 3)    #衛星速度＠地心直交座標系

#シミュレーションループ
orbp, eop_IAU2000A, o, r, v = set_SPCal(TLEFileName, t0)

for i=1:DataNum
    x_sat, v_sat = SPCal(i, t0, dt, orbp, eop_IAU2000A)
    x_ecef_log[i,:] = x_sat
    v_ecef_log[i,:] = v_sat
    if i != 1
        print(x_ecef_log[i,1])
        print(",")
        print(x_ecef_log[i,2])
        print(",")
        println(x_ecef_log[i,3])
        #print(",")
        #print((x_ecef_log[i,1]-x_ecef_log[i-1,1])./dt)
        #print(",")
        #print((x_ecef_log[i,2]-x_ecef_log[i-1,2])./dt)
        #print(",")
        #println((x_ecef_log[i,3]-x_ecef_log[i-1,3])./dt)
    end
    x_geo = ECEFtoGeodetic([x_sat[1]; x_sat[2]; x_sat[3]])
    x_geod_log[i,1]=x_geo[1]
    if x_geo[2]>0
        x_geod_log[i,2]=x_geo[2]
    else
        x_geod_log[i,2]=x_geo[2]+2π
    end
    x_geod_log[i,2]=x_geo[2]
    x_geod_log[i,3]=x_geo[3]
    #println(x_geod_log[i,:])
end

#=
#=************************************************
地心直交座標系での衛星位置を3Dプロットしたい場合はコレ
************************************************=#
println(x_ecef_log)

plot(x_ecef_log[:,1], x_ecef_log[:, 2], x_ecef_log[:, 3])
xlabel!("x")
ylabel!("y")
#zlabel!("z")

savefig("test2.png")
=#

#=***********************************************
緯度経度に変換してGifでグラフを作成したい場合はコレ
***********************************************=#
#=
anim = @animate for i=1:DataNum
    x_sat, v_sat = SPCal(i, t0, dt, orbp, eop_IAU2000A)
    x_ecef_log[i,:] = x_sat
    v_ecef_log[i,:] = v_sat
    x_geo = ECEFtoGeodetic([x_ecef_log[i,1]; x_ecef_log[i,2]; x_ecef_log[i,3]])
    x_geod_log[i,1]=x_geo[1]
    if x_geo[2]>0
        x_geod_log[i,2]=x_geo[2]
    else
        x_geod_log[i,2]=x_geo[2]+2π
    end
    x_geod_log[i,3]=x_geo[3]
    println(i)
    plot(rad2deg.(x_geod_log[1:i,2]),rad2deg.(x_geod_log[1:i,1]))
end

gif(anim, "test5_anim.gif", fps = 10)
=#

max_stars_repo_name: yui-project/ACSim

max_issues_repo_name: yui-project/ACSim

max_forks_repo_name: yui-project/ACSim


text:
Pkg.add("CSV")
Pkg.add("DataFrames")
Pkg.add("PowerModels")
Pkg.add("Ipopt")
Pkg.update()

max_stars_repo_name: bishnubhattarai/PFO-ADC-DER-Testbed

max_issues_repo_name: bishnubhattarai/PFO-ADC-DER-Testbed

max_forks_repo_name: bishnubhattarai/PFO-ADC-DER-Testbed


text:
module Particles

using OffsetArrays

export Simulation, run!

include("simulation.jl")

end # module

max_stars_repo_name: michakraus/JuliaAtHomeAssignments

max_issues_repo_name: michakraus/JuliaAtHomeAssignments

max_forks_repo_name: michakraus/JuliaAtHomeAssignments


text:
export
    Plane3,
    get_signed_distance

"""
Plane3

A hyperplane is an affine subspace of dimension n-1 in a space of dimension n.
For example, a hyperplane in a plane is a line; a hyperplane in 3-space is a plane.

This class represents an hyperplane in 3D space. It is the zero set of
n⋅x + d = 0 where n is a unit normal vector of the plane (linear part)
and d is the distance (offset) to the origin.
"""
struct Plane3
    normal::VecE3 # should always be normalized
    offset::Float64

    Plane3() = new(VecE3(1.0,0.0,0.0), 0.0)

    """
    Constructs a plane from its normal n and distance to the origin d
    such that the algebraic equation of the plane is n⋅x + d = 0.
    """
    Plane3(normal::VecE3, offset::Real) = new(normalize(normal), convert(Float64, offset))

    """
    Construct a plane from its normal and a point on the plane.
    """
    function Plane3(normal::VecE3, P::VecE3)
        n = normalize(normal)
        offset = -n⋅P
        return new(n, offset)
    end
end

"""
Constructs a plane passing through the three points.
"""
function Plane3(p0::VecE3, p1::VecE3, p2::VecE3)
    v0 = p2 - p0
    v1 = p1 - p0

    normal = v0×v1
    nnorm = norm(normal)
    if nnorm <= norm(v0)*norm(v1)*eps()
        M = hcat(convert(Vector{Float64}, v0),
                 convert(Vector{Float64}, v1))'
        s = svdfact(M, thin=false)
        normal = convert(VecE3, s[:V][:,2])
    else
        normal = normalize(normal)
    end

    offset = -p0⋅normal

    return Plane3(normal, offset)
end

"""
Returns the signed distance between the plane and a point
"""
get_signed_distance(plane::Plane3, A::VecE3) = plane.normal⋅A + plane.offset

"""
Returns the absolute distance between the plane and a point
"""
get_distance(plane::Plane3, A::VecE3) = abs(get_signed_distance(plane, A))

"""
Returns the projection of a point onto the plane
"""
proj(A::VecE3, plane::Plane3) = A - get_signed_distance(plane, A) * plane.normal

"""
What side of the plane you are on
Is -1 if on the negative side, 1 if on the positive side, and 0 if on the plane
"""
function get_side(plane::Plane3, A::VecE3, ε::Float64=DUMMY_PRECISION)
    signed_dist = get_signed_distance(plane, A)
    if abs(signed_dist) < ε
        return 0
    else
        return sign(signed_dist)
    end
end


text:
using Revise
using Bolt
using ForwardDiff
#using PyPlot
using Plots
pyplot()
using BenchmarkTools
using Printf
using QuadGK
#using ThreadPools

#input ingredients
par = CosmoParams()
logqmin,logqmax = -6,-1
n_q = 15
logq_pts = logqmin:(logqmax-logqmin)/(n_q-1):logqmax
bg = Background(par;x_grid=-20.0:0.1:0.0,nq=n_q) #,logq_grid=logq_pts)
ih = IonizationHistory(Peebles(), par, bg)
k_grid = quadratic_k(0.1bg.H₀, 1000bg.H₀, 100) #quadratically spaced k points

# #quickly check bg against nr approx
# Tγ = (15/ π^2 *bg.ρ_crit *par.Ω_r)^(1/4)
# νfac = (90 * 1.2020569 /(11 * π^4)) * (par.Ω_r * par.h^2 / Tγ)#the factor that goes into nr approx to neutrino energy density
# par.Σm_ν*νfac/par.h^2 *((par.N_ν/3)^(3/4)) /(ρP_0(1,par)[1]/bg.ρ_crit)
# (par.N_ν/3)^(3/4)

function f00(q)
    Tν =  (par.N_ν/3)^(1/4) * (4/11)^(1/3) * (15/ π^2 *3.9669896e-11 *5.042e-5)^(1/4) ##assume instant decouple for now
    gs =  2 #should be 2 for EACH neutrino family (mass eigenstate)
    return gs / (2π)^3 / ( exp(q/Tν) +1)
end

function dlnf0dlnq0(q) #this is actually only used in perts
    Tν =  (par.N_ν/3)^(1/4) * (4/11)^(1/3) * (15/ π^2 *3.9669896e-11 *5.042e-5)^(1/4) ##assume instant decouple for now
    return -q / Tν /(1 + exp(-q/Tν))
end

#find correct factor for normalization...
ρν0 = 7*(2/3)*par.N_ν/8 *(4/11)^(4/3) *par.Ω_r * bg.ρ_crit / 2 #used to be div by Neff, now div by floor(Neff-1)
ρν=ρν0*(exp(-20))^(-4)
aaa=4π  * quadgk(q ->  q^2 * -dlnf0dlnq0(q) *q * f00(q),
            1e-6, 1e-1,rtol=1e-6)[1]/4/ρν0
#!
aaa
# Removed the splines from bg so these won't work anymore
# #check the splining error:
# #use both splines - error of ~1.8e-3
# 4π  * quadgk(q ->  q^2 * -bg.df0(log10(q)) *q * bg.f0(log10(q)),
#             1e-6, 1e-1,rtol=1e-6)[1]/4/ρν0
#
# #only use df0 spline - error of ~ 5e-6, ~ rtol
# 4π  * quadgk(q ->  q^2 * -bg.df0(log10(q)) *q * f00(q),
#             1e-6, 1e-1,rtol=1e-6)[1]/4/ρν0
#
# #only use f0 spline - again we get error of ~1.8e-3 - so this is the problem...
# 4π  * quadgk(q ->  q^2 * -dlnf0dlnq0(q) *q * bg.f0(log10(q)),
#             1e-6, 1e-1,rtol=1e-6)[1]/4/ρν0

#test that ρ_σ is the same as bg when passed ones - it is up to quadgk tol...
bgrho,_ =  (exp(-20)^(-4)) .* ρ_σ(ones(n_q) ,
               zeros(n_q),bg,exp(-20),par)
ρP_0(exp(-20),par)
ρν #analytic answer

#@btime @qthreads
ℓᵧ=8 #cutoff
ℓ_ν=10 #not used except for size here, should pass
ℓ_mν=ℓ_ν
reltol=1e-5 #cheaper  rtol
#solve hierarchy at single x to check
x=-8 #just picking a number
a=exp(x)
pertlen = 2(ℓᵧ+1)+(ℓ_ν+1)+(ℓ_mν+1)*n_q+5
println("pert vector length=",pertlen)
results=zeros(pertlen,length(k_grid))
n_q
for (i_k, k) in enumerate(k_grid)
    println(i_k)
    hierarchy = Hierarchy(BasicNewtonian(), par, bg, ih, k, ℓᵧ, n_q)
    perturb = boltsolve(hierarchy; reltol=reltol)
    u = perturb(x)  #z this can be optimized away, save timesteps at the grid!
    results[:,i_k] = u #z should use unpack somehow
end
results

#Integrate the q moments for ℳ0 and ℳ2 for plotting
ℳρ,ℳσ = zeros(length(k_grid)),zeros(length(k_grid))
ℳθ = zeros(length(k_grid))
Ω_ν =  7*(2/3)*par.N_ν/8 *(4/11)^(4/3) *par.Ω_r
norm𝒩 = 1/(4Ω_ν * bg.ρ_crit / 2)
for (i_k, k) in enumerate(k_grid)
    ℳρ[i_k],ℳσ[i_k] = ρ_σ(results[2(ℓᵧ+1)+(ℓ_ν+1)+1:2(ℓᵧ+1)+(ℓ_ν+1)+n_q,i_k],
                            results[2(ℓᵧ+1)+(ℓ_ν+1)+2*n_q+1:2(ℓᵧ+1)+(ℓ_ν+1)+3*n_q,i_k],
                            bg,a,par)#.*norm𝒩
    ℳθ[i_k],_ = ρ_σ(results[2(ℓᵧ+1)+(ℓ_ν+1)+n_q+1:2(ℓᵧ+1)+(ℓ_ν+1)+2*n_q,i_k],
                            zeros(n_q),
                            bg,a,par)

end
labels = [raw"$\Phi$",raw"$\delta$",raw"$v$",raw"$|\delta_{b}|$",raw"$|v_{b}|$"]
plot(legend=:bottomleft)
#title!("ICs")
for i in 1:5
    plot!(log10.(k_grid/ bg.H₀),
          log10.(abs.(results[2(ℓᵧ+1)+(ℓ_ν+1)+(ℓ_mν+1)*n_q+i,:])),
          label=labels[i])
end
plot!(log10.(k_grid/ bg.H₀), log10.(abs.(results[1,:])),label=raw"$|\Theta_{0}|$")

plot!(log10.(k_grid/ bg.H₀), log10.(abs.( results[2(ℓᵧ+1)+1,:])),
      label=raw"$|\mathcal{N}_{0}|$")
plot!(log10.(k_grid/ bg.H₀), log10.(abs.(results[2(ℓᵧ+1)+2,:])),
      label=raw"$|\mathcal{N}_{1}|$")
plot!(log10.(k_grid/ bg.H₀), log10.(abs.(results[2(ℓᵧ+1)+3,:])),
      label=raw"$|\mathcal{N}_{2}|$")

results[2(ℓᵧ+1)+1,:] ./ ℳρ
results[2(ℓᵧ+1)+2,:] ./ ℳθ
results[2(ℓᵧ+1)+3,:] ./ ℳσ

plot!(log10.(k_grid/ bg.H₀),log10.(abs.(ℳρ)),
      label=raw"$|\mathcal{M}_{0}|$",ls=:dash)
plot!(log10.(k_grid/ bg.H₀),log10.(abs.(ℳθ)),
    label=raw"$|\mathcal{M}_{1}|$",ls=:dash)
plot!(log10.(k_grid/ bg.H₀),log10.(abs.(ℳσ)),
      label=raw"$|\mathcal{M}_{2}|$",ls=:dash)

ylabel!(raw"$\delta_{i}(k)$")
xlabel!(raw"$k / H0$")
title!("z=$(@sprintf("%.0f", exp(-x)-1))")
savefig("../compare/bolt_perts_k_z$(@sprintf("%.0f", exp(-x)-1)).png")


text:
""" [Khairoutdinov2009](@cite) """
function GATE_III_z(::Type{FT}) where {FT}
    # TODO: should the last element in here be 22, not 27, km?
    # TODO: unify with GATE_III_z_in
    z_in = FT[ 0.0,   0.5,  1.0,  1.5,  2.0,   2.5,    3.0,   3.5,   4.0,   4.5,   5.0,  5.5,  6.0,  6.5,
               7.0, 7.5, 8.0,  8.5,   9.0,   9.5, 10.0,   10.5,   11.0, 11.5, 12.0, 12.5, 13.0, 13.5, 14.0,
               14.5, 15.0, 15.5, 16.0, 16.5, 17.0, 17.5, 18.0, 27.0] .* 1000 #z is in meters
    return z_in
end

""" [Khairoutdinov2009](@cite) """
function GATE_III_q_tot(::Type{FT}) where {FT}
    z_in = GATE_III_z(FT)
    r_in = FT[16.5,  16.5, 13.5, 12.0, 10.0,   8.7,    7.1,   6.1,   5.2,   4.5,   3.6,  3.0,  2.3, 1.75, 1.3,
               0.9, 0.5, 0.25, 0.125, 0.065, 0.003, 0.0015, 0.0007,  0.0003,  0.0001,  0.0001,  0.0001,  0.0001,
               0.0001,  0.0001,  0.0001,  0.0001,  0.0001,  0.0001,  0.0001,  0.0001,  0.0001, 0.0001] ./ 1000 # mixing ratio should be in kg/kg
    q_tot = r_in ./ (1 .+ r_in) # convert mixing ratio to specific humidity
    profile = Dierckx.Spline1D(z_in, q_tot; k = 1)
    return profile
end

""" [Khairoutdinov2009](@cite) """
function GATE_III_u(::Type{FT}) where {FT}
    z_in = GATE_III_z(FT)
    U_in = FT[  -1, -1.75, -2.5, -3.6, -6.0, -8.75, -11.75, -13.0, -13.1, -12.1, -11.0, -8.5, -5.0, -2.6, 0.0,
               0.5, 0.4,  0.3,   0.0,  -1.0, -2.5,   -3.5,   -4.5, -4.8, -5.0, -3.5, -2.0, -1.0, -1.0, -1.0,
               -1.5, -2.0, -2.5, -2.6, -2.7, -3.0, -3.0, -3.0] # [m/s]
    profile = Dierckx.Spline1D(z_in, U_in; k = 1)
    return profile
end

""" [Khairoutdinov2009](@cite) """
function GATE_III_T(::Type{FT}) where {FT}
    # temperature is taken from a different input plot at different z levels
    z_in = GATE_III_z(FT)
    T_in = FT[299.184, 294.836, 294.261, 288.773, 276.698, 265.004, 253.930, 243.662, 227.674, 214.266, 207.757, 201.973, 198.278, 197.414, 198.110, 198.110]
    z_T_in = FT[0.0, 0.492, 0.700, 1.698, 3.928, 6.039, 7.795, 9.137, 11.055, 12.645, 13.521, 14.486, 15.448, 16.436, 17.293, 22.0] .* 1000 # for km
    profile = Dierckx.Spline1D(z_T_in, T_in; k = 1)
    return profile
end

""" [Khairoutdinov2009](@cite) """
function GATE_III_tke(::Type{FT}) where {FT}
    return z -> if z <= 2500.0
            FT(1) - z / 3000.0
        else
            FT(0)
        end
end

#=

Unused data:

u_in = FT[-1, -1.75, -2.5, -3.6,  -6.0, -8.75, -11.75,
          -12.9, -13.1, -12.1, -11.0,  -8.5,  -5.0, -2.6,
          0.0, 0.5, 0.4, 0.3, 0.0, -1.0, -3.0, -3.5, -4.5,
          -4.6, -5.0, -3.5, -2.0, -1.0, -1.0, -1.0, -1.5,
          -2.0, -2.5, -2.6, -2.7, -3.0, -3.0]
=#
""" [Khairoutdinov2009](@cite) """
function GATE_III_z_in(::Type{FT}) where {FT}
    #LES z is in meters
    z_in = FT[ 0.0,  0.5,  1.0,  1.5,  2.0,  2.5,  3.0,  3.5,  4.0,  4.5, 5.0,  5.5,   6.0,
               6.5,  7.0,  7.5,  8.0,  8.5,  9.0,  9.5, 10.0, 10.5, 11.0, 11.5, 12.0, 12.5,
              13.0, 13.5, 14.0, 14.5, 15.0, 15.5, 16.0, 16.5, 17.0, 17.5, 18.0] .* 1000.0
    return z_in
end

""" [Khairoutdinov2009](@cite) """
function GATE_III_dqtdt(::Type{FT}) where {FT}
    z_in = GATE_III_z_in(FT)
    # Advective qt forcing  for theta [g/kg/d] converted to [kg/kg/sec]
    r_tend_in = FT[ 0.0,   1.2,  2.0,  2.3,   2.2,   2.1,    1.9,   1.7,   1.5,  1.35,  1.22,  1.08,  0.95,
                   0.82,  0.7,  0.6,   0.5,  0.4,   0.3,  0.2,   0.1,  0.05, 0.0025, 0.0012, 0.0006,  0.0,
                   0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0] ./ (24 * 3600) ./ 1000

    Qtend_in = r_tend_in ./ (1 .+ r_tend_in) # convert mixing ratio to specific humidity
    profile = Dierckx.Spline1D(z_in, Qtend_in; k = 1)
    return profile
end

""" [Khairoutdinov2009](@cite) """
function GATE_III_dTdt(::Type{FT}) where {FT}
    z_in = GATE_III_z_in(FT)
    # Radiative forcing for T [K/d] converted to [K/sec]
    RAD_in = FT[-2.9,  -1.1, -0.8, -1.1, -1.25, -1.35,   -1.4,  -1.4, -1.44, -1.52,  -1.6, -1.54, -1.49,
               -1.43, -1.36, -1.3, -1.25, -1.2, -1.15, -1.1, -1.05,  -1.0,  -0.95,   -0.9,  -0.85, -0.8,
               -0.75, -0.7, -0.6, -0.3,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0] ./ (24 * 3600)
    # Advective qt forcing  for theta [g/kg/d] converted to [kg/kg/sec]
    r_tend_in = FT[ 0.0,   1.2,  2.0,  2.3,   2.2,   2.1,    1.9,   1.7,   1.5,  1.35,  1.22,  1.08,  0.95,
                   0.82,  0.7,  0.6,   0.5,  0.4,   0.3,  0.2,   0.1,  0.05, 0.0025, 0.0012, 0.0006,  0.0,
                   0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0] ./ (24 * 3600) ./ 1000
    # Radiative T forcing [K/d] converted to [K/sec]
    Ttend_in = FT[ 0.0,  -1.0, -2.2, -3.0,  -3.5,  -3.8,   -4.0,  -4.1,  -4.2,  -4.2,  -4.1,  -4.0, -3.85,
                  -3.7, -3.5, -3.25, -3.0, -2.8,  -2.5, -2.1,  -1.7,  -1.3,   -1.0,   -0.7,   -0.5, -0.4,
                  -0.3, -0.2, -0.1,-0.05,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0] ./ (24*3600)

    profile_T = Dierckx.Spline1D(z_in, Ttend_in; k = 1)
    profile_R = Dierckx.Spline1D(z_in, RAD_in; k = 1)
    return z -> profile_T(z) + profile_R(z)
end


text:
using Fortran90Namelists
using Test

@testset "Fortran90Namelists.jl" begin
    include("FortranToJuliaTests.jl")
    include("TokenizeTests.jl")
end