text:
function isexpr(ex, symbol)
  return ex.head == symbol
end

find_symbols(x::Number) = []
find_symbols(x::Symbol) = [x]
find_symbols(ex) = reduce(vcat, [], map(find_symbols, ex.args[2:end]))

macro fill_series(ex)
  l, r = ex.args # left- and right-hand sides of equation
  println(string("L: ", l, " R: ", r) )
  var_t =                   # expecting left to be array[expression]
                            # so second arg is expression in brackets
    if isa(l.args[2], Symbol)
      l.args[2]               # it's just a variable name so bind it
      println(string("args: ", l.args, " is a symbol"))
    elseif isexpr(l.args[2], :call) #if it's a call it's something like t+1
      println(string("args: ", l.args, " is a call"))
      #find first symbol that occurs in this call
      find_symbols(l)[1]
    end
  println(string("variable is: ", var_t))
  # find biggest and smallest offsets from left and right
  offsets = extrema(vcat(get_offsets(l), get_offsets(r)))
  println(string("Got offsets: ", offsets))
  #get range with expression of (1-smallest_offset):(legth_of_given_array-biggest_offset)
  loopindex = :($(1 - offsets[1]):(length($(l.args[1])) - $(offsets[2])))
  println(loopindex)
  println(string("for ", var_t, " in ", loopindex))
  println(string("    ", ex))
  quote
    for $var_t in $loopindex
      $ex
    end
  end
end

function get_offsets(ex_::Expr)
  println(string(ex_, " is of type ", typeof(ex_)))
  if isexpr(ex_,:call)
    println(string("ex ", ex_, " is a call"))
    # it's a call, so first arg is operator and the rest are its arguments so
    # we need to find refs in it
     return [[get_offsets(a) for a in ex_.args[2:end]]...]
   end

  if isexpr(ex_,:ref)
     println(string("ex ", ex_, " is a ref with arg 2: ", ex_.args[2]))
     return get_offset_from_ref(ex_.args[2])
   end

  warn("Not expecting to be here")
  return Int64[]
end

get_offsets(x) = Int64[]

get_offset_from_ref(s::Symbol) = 0
get_offset_from_ref(x::Number) = x

function get_offset_from_ref(ex_::Expr)
  if isexpr(ex_,:call)
    if ex_.args[1] == :+
      println(string("ex ", ex_, " arg 1 is + so gonna sum offsets from ", ex_.args[2:end]))
      return sum([get_offset_from_ref(a) for a in ex_.args[2:end]])
    end

    if ex_.args[1] == :-
      println(string("ex ", ex_, " arg 1 is - so gonna sub offsets from ", ex_.args[2:end]))
      #expression is in form of y[t-2] so get offset of t(=0) and sub offset of 2(=2)
      return (get_offset_from_ref(ex_.args[2]) - sum([get_offset_from_ref(a) for a in ex_.args[3:end]]))
    end
  end
  warn("Didn’t expect to get here")
  return(0)
end


text:
mutable struct LightBox{T,Int} <: AbstractEnvironment where {T}
    state::Array{T, 1}
    reward::Float64
    done::Bool
    time_step::Int
	stochastic::Bool
	randchance::T

function LightBox(::Type{T}, stochastic::Bool=true, chance::Float64=0.9) where{T}
        state = zeros(T, 20)
        new{T,Int}(state, 0., false, 0, stochastic, chance)
    end
end

function get_gamma(env::LightBox)
	return 1.0
end

function get_state_desc(env::LightBox{T})::Array{T,2} where {T}
    ranges = zeros(T, (20, 2))
    ranges[:, 2] .= T(1.)
	return ranges
end

function get_action_desc(env::LightBox)::Int
    return 20
end

function step!(env::LightBox, action::Int, rng::AbstractRNG)
    env.time_step += 1
    if action <= 0 || action > 20
        error("Action needs to be an integer in [1, 20]")
    end

    #values must be on for key to be on
    dependencies = Dict{Int,Tuple}(9=>(0, 3, 6), 10=>(1, 4), 11=>(2, 5),
                    12=>(2, 1), 13=>(5, 4), 14=>(8, 7, 6),
                    15=>(9, 10), 16=>(10, 11), 17=>(12, 13),
                    18=>(13, 14), 19=>(16, 17))

    #key must be on for corresponding values to be on
    causal = Dict{Int,Tuple}(0=>(0, 9), 1=>(1, 10, 12), 2=>(2, 11, 12), 3=>(3, 9),
              4=>(4,10,13), 5=>(5, 11, 13), 6=>(6, 9, 14), 7=>(7, 14),
              8=>(8, 14), 9=>(9, 15), 10=>(10, 15, 16), 11=>(11, 16),
	         12=>(12, 17), 13=>(13, 17, 18), 14=>(14, 18), 15=>(15,),
             16=>(16, 19), 17=>(17, 19), 18=>(18,), 19=>(19,))

    if !env.stochastic || rand(rng) < env.randchance
        # update state
        if is_legal_move(env.state, action, dependencies)
            env.state[action] = 1.0 - env.state[action]
            update_lights!(env.state, action, causal)
        else
            env.state .= 0
        end
	end

    env.reward = -1.0
	max_steps = 5000
    env.done = env.time_step ≥ max_steps
    if isapprox(env.state[end], 1.0)
        env.reward = 1.0
        env.done = true
    end
    return env.reward
end

function is_legal_move(state, action, dependencies)
    if isapprox(state[action], 1.) || action < 10
        return true #always legal to turn off a light and toggle a base light
    end

    for i in dependencies[action-1]
        if isapprox(state[i+1], 0.)
            return false
        end
    end
    return true
end

function update_lights!(state, action, causal)
    #if light is off, turn off all dependent lights
    if isapprox(state[action], 0.)
        for i in causal[action-1][2:end]
            state[i+1] = 0.
            update_lights!(state, i+1, causal)
        end
    end
end

function get_state(env::LightBox)
    return env.state
end

function get_reward(env::LightBox)::Float64
    return env.reward
end

function is_terminal(env::LightBox)::Bool
    return env.done
end

function env_reset!(env::LightBox, rng)
    env.time_step = 0
    env.state .= 0.0
    env.done = false
end

function clone(env::LightBox)::LightBox
    env2 = LightBox(env.stochastic, env.randchance)
    env2.time_step = deepcopy(env.time_step)
    env2.state = deepcopy(env.state)
    env2.done = deepcopy(env.done)
    env2.reward = deepcopy(env.reward)
    return env2
end

ext: jl

max_issues_repo_name: ScottJordan/EvaluationOfRLAlgs

max_forks_repo_name: ScottJordan/EvaluationOfRLAlgs


text:
using Random
using Statistics
using Test

# Include testing functions
include("test_dimensions.jl")
include("test_symbolic_functions.jl")
include("test_demo_large.jl")

_testset_output_dim("large")
_testset_symbolic_functions(
    "large", ["amg_rs", "amg_sa", "diagonal", "none"], ["abs", "square"]
)
_testset_demo_large("large", ["amg_rs"], ["abs"])

@testset "Broken inputs" begin
    @test_throws ArgumentError tvdiff([0, 1, 2], 1, 0.1; scale="large", precond="bad_input")
    @test_throws ArgumentError tvdiff(
        [0, 1, 2], 1, 0.1; scale="large", diff_kernel="bad_input"
    )
    @test_throws MethodError tvdiff([0, 1, 2], 1, 0.1; scale="l4rge")
    @test_throws DimensionMismatch tvdiff(
        [0, 1, 2], 1, 0.1; u_0=[1, 2, 3, 4], scale="large"
    )
end


text:
module ArnoldiMethodInterface

using ArnoldiMethod
using ArnoldiMethodTransformations
using ..Common

import ..AbstractEigenproblem
import ..AbstractLinearEigenproblem
import ..AbstractCFEigenproblem
import ..AbstractNonlinearEigenproblem
import ..HelmholtzProblem
import ..MaxwellProblem
import ..DEFAULT_LINEAR_EIGENSOLVER
import ..HelmholtzLEP
import ..MaxwellLEP

################################################################################
#LEP

"""
	partialschur(::AbstractLinearEigenproblem, ω; lupack=$DEFAULT_LUPACK, kwargs...) -> decomp, history

	partialschur(::AbstractLinearEigenproblem, ω, ωs::Array, ψs::ElectricField; lupack=$DEFAULT_LUPACK, kwargs...) -> decomp, history
"""
function ArnoldiMethod.partialschur(
			lep::AbstractLinearEigenproblem,
			ω::Number,
			args...;
			lupack::AbstractLUPACK=DEFAULT_LUPACK,
			kwargs...)

	A, B = lep(ω,args...)
	return partialschur(A, B, ω^2; diag_inv_B=true, lupack=lupack, kwargs...)
end

"""
	partialeigen(lep::AbstractLinearEigenproblem, decomp, ω) -> ωs, electricfields

order of `decomp` and `ω` does not matter
"""
function ArnoldiMethod.partialeigen(
			lep::AbstractLinearEigenproblem,
			ω::Number,
			decomp)

	ωs::Vector{ComplexF64}, ψ::Matrix{ComplexF64} = partialeigen(decomp, ω^2)
	normalize!(lep.sim, ψ, lep.εpχ, size(ψ,1)÷2+INDEX_OFFSET) # Normalize according to (ψ₁,ψ₂)=δ₁₂
	if typeof(lep) <: HelmholtzLEP
		Ψ = ScalarField(lep,ψ)
	elseif typeof(lep) <: MaxwellLEP
		Ψ = ElectricField(lep,ψ)
	end
	return sqrt.(ωs), Ψ
end
ArnoldiMethod.partialeigen(lep::AbstractLinearEigenproblem,decomp,ω::Number) = partialeigen(lep,ω,decomp)

function iris_eigen_arnoldimethod(
            lep::AbstractLinearEigenproblem,
            ω::Number,
			args...;
            verbose::Bool=false,
            lupack::AbstractLUPACK=DEFAULT_LUPACK,
            kwargs...)

	decomp, history = partialschur(lep, ω, args...; diag_inv_B=true, lupack=lupack, kwargs...)
    history.converged || @warn "$(history.nev - history.nconverged) evecs did not converge"
    verbose ? println(history) : nothing
    return partialeigen(lep, ω, decomp)
end

if DEFAULT_LINEAR_EIGENSOLVER == :ArnoldiMethod
@doc """
$doc_lep_h
`nev` Number of eigenvalues (`6`);
`verbose` Print convergence info (`false`)

Advanced keywords:
`lupack` LU solver (`$DEFAULT_LUPACK`);
`which` See `ArnoldiMethod` documentation (`LM()`);
`tol` Tolerance for convergence (`√eps`);
`mindim` Minimum Krylov dimension (`max(10, nev)`);
`maxdim` Maximum Krylov dimension (`max(20, 2nev)`);
`restarts` Maximum number of restarts (`200`)
""" ->
helmholtzeigen(lep::HelmholtzLEP,args...;kwargs...) = iris_eigen_arnoldimethod(lep,args...;kwargs...)

@doc """
$doc_lep_m
`nev` Number of eigenvalues (`6`);
`verbose` Print convergence info (`false`)

Advanced keywords:
`lupack` LU solver (`$DEFAULT_LUPACK`);
`which` See `ArnoldiMethod` documentation (`LM()`);
`tol` Tolerance for convergence (`√eps`);
`mindim` Minimum Krylov dimension (`max(10, nev)`);
`maxdim` Maximum Krylov dimension (`max(20, 2nev)`);
`restarts` Maximum number of restarts (`200`)
""" ->
maxwelleigen(lep::MaxwellLEP,args...;kwargs...) = iris_eigen_arnoldimethod(lep,args...;kwargs...)
end

################################################################################
#CF

"""
	partialschur(::AbstractCFEigenproblem, ω; η=0, lupack=$DEFAULT_LUPACK, kwargs...) -> decomp, history

	partialschur(::AbstractCFEigenproblem, ω, ωs, ψs; η=0, lupack=$DEFAULT_LUPACK, kwargs...) -> decomp, history
"""
function ArnoldiMethod.partialschur(
			cf::AbstractCFEigenproblem,
			ω::Number,
			args...;
			η::Number=0,
			lupack::AbstractLUPACK=DEFAULT_LUPACK,
			kwargs...)

	A, B = cf(ω,args...)
	return partialschur(A, B, η; diag_inv_B=true, lupack=lupack, kwargs...)
end

"""
	partialeigen(cf::AbstractCFEigenproblem, decomp, η) -> ηs, us
"""
function ArnoldiMethod.partialeigen(
			cf::AbstractCFEigenproblem,
			η::Number,
			decomp)

	ηs::Vector{ComplexF64}, us::Matrix{ComplexF64} = partialeigen(decomp, η)
	normalize!(cf.maxwell.sim, us, cf.F)
	if typeof(cf) <: HelmholtzCF
		U = ScalarField(cf,us)
	elseif typeof(cf) <: MaxwellCF
		U = ElectricField(cf,us)
	end
	return sqrt.(ηs), U
end
ArnoldiMethod.partialeigen(cf::AbstractCFEigenproblem,decomp,η::Number) = partialeigen(cf,η,decomp)

function iris_eigen_arnoldimethod(
            cf::AbstractCFEigenproblem,
            ω::Number,
			args...;
            η::Number=0,
            verbose::Bool=false,
            lupack::AbstractLUPACK=DEFAULT_LUPACK,
			kwargs...)

    decomp, history = partialschur(cf, ω, args...; η=η, diag_inv_B=false, lupack=lupack, kwargs...)
    history.converged || @warn "$(history.nev-history.nconverged) eigenvectors failed to converge"
    verbose ? println(history) : nothing
    return partialeigen(cf, η, decomp)
end

if DEFAULT_LINEAR_EIGENSOLVER == :ArnoldiMethod
@doc """
$doc_cf_h
`nev` Number of eigenvalues (`6`);
`verbose` Print convergence info (`false`);

Advanced keywords:
`lupack` LU solver (`$DEFAULT_LUPACK`);
`which` See `ArnoldiMethod` documentation (`LM()`);
`tol` Tolerance for convergence (`√eps`);
`mindim` Minimum Krylov dimension (`max(10, nev)`);
`maxdim` Maximum Krylov dimension (`max(20, 2nev)`);
`restarts` Maximum number of restarts (`200`)
""" ->
helmholtzeigen(cf::HelmholtzCF,args...;kwargs...) = iris_eigen_arnoldimethod(cf,args...;kwargs...)

@doc """
$doc_cf_m
`nev` Number of eigenvalues (`6`);
`verbose` Print convergence info (`false`);

Advanced keywords:
`lupack` LU solver (`$DEFAULT_LUPACK`);
`which` See `ArnoldiMethod` documentation (`LM()`);
`tol` Tolerance for convergence (`√eps`);
`mindim` Minimum Krylov dimension (`max(10, nev)`);
`maxdim` Maximum Krylov dimension (`max(20, 2nev)`);
`restarts` Maximum number of restarts (`200`)
""" ->
maxwelleigen(cf::MaxwellCF,args...;kwargs...) = iris_eigen_arnoldimethod(cf,args...;kwargs...)
end

end

using .ArnoldiMethodInterface


text:
function searchsortedmapping(A::Array, B::Array)
    ## Assume A \in B
    AtoB = Array{Int64,1}()
    temp_index_A = 1
    temp_index_B = 1

    while (temp_index_A <= length(A))
        if A[temp_index_A] == B[temp_index_B]
            push!(AtoB,temp_index_B)
            temp_index_A += 1
            temp_index_B += 1
        else
            temp_index_B += 1
        end
    end

    return AtoB
end

########################################################
# D4M: Dynamic Distributed Dimensional Data Model
# Architect: Dr. Jeremy Kepner (kepner@ll.mit.edu)
# Software Engineer: Alexander Chen (alexc89@mit.edu)
########################################################


text:
# */ --------------------------------------------------------------------------------------------- /* #
# */ --- help functions -------------------------------------------------------------------------- /* #
# */ --------------------------------------------------------------------------------------------- /* #

@everywhere function circ(h::Int64,H::UnitRange{Int64})::Int64

    (h==minimum(H)) ? idx = maximum(H) : idx = h-1

    return idx
end

# */ --------------------------------------------------------------------------------------------- /* #
# */ --- optimization models --------------------------------------------------------------------- /* #
# */ --------------------------------------------------------------------------------------------- /* #

function RMP!(m::JuMP.Model,ms::ms_type,mp::mp_type,unc::u_type)::JuMP.Model

    # problem structure
    # min   f + ∑ᵢ πᵢ βᵢ)
    # s.t.  f = q₀ᵀx₀ + ∑ᵢ qᵢᵀxᵢ
    #      (x₀,x₁,..,xᵢ,..,xₙ) ∈ 𝕏
    #       βᵢ ∈ Θᵢ,  ∀i

    xr=1:unc.nx
    # */ -- variables ---------------------------------------------- /* #
    @variable(m, f, container=Array)
    @variable(m, x0[ms.P,ms.I0] >= .0, container=Array)
    @variable(m, x[xr,ms.I], container=Array)
    @variable(m, β[ms.I] >= .0, container=Array)
    # */ ----------------------------------------------------------- /* #

    # */ -- obj function ------------------------------------------- /* #
    @objective(m, Min, f + mp.κ*sum(mp.π[i]*β[i] for i in ms.I) )
    # */ ----------------------------------------------------------- /* #

    # */ -- constraints -------------------------------------------- /* #
    @constraint(m, c01, f >= exp10(-6)*sum(mp.π0[i0]*sum(mp.ci[p,i0]*x0[p,i0] for p in ms.P) for i0 in ms.I0) + exp10(-6)*mp.κ*sum(mp.π[i]*sum(mp.cf[p]*x[p,i] for p in ms.P) for i in ms.I) )
    @constraint(m, c02[p=ms.P,i=ms.I], x[p,i] == mp.xh[p,i] + sum(x0[p,i0] for i0 in mp.map[i]) )
    @constraint(m, c03[p=ms.P,i=ms.I], x[p,i] >= .0 )
    @constraint(m, c04[p=ms.P,i=ms.I], x[p,i] <= mp.xm[p] )
    for i in ms.I for j in 1:unc.nh
        fix(x[unc.nx0+j,i],unc.h[i,j];force=true)
    end end
    # */ ----------------------------------------------------------- /* #
    optimize!(m)

    return m
end

@everywhere function SP!(m::JuMP.Model,ps::ps_type,pp::pp_type,unc::u_type)::JuMP.Model

    # problem structure
    # min   c₀ + cᵀϕ
    # s.t.  c₀  = C₀y
    #       ϕ   = C y
    #       y ∈ 𝕐
    #       A y ≦ B x
    #       x = xf :(λ)

    # */ -- variables ---------------------------------------------- /* #
    @variable(m, ϕ[1:unc.nc], container=Array)
    @variable(m, c0, container=Array)
    @variable(m, yG[ps.G,ps.S,ps.H] >= 0, container=Array)
    @variable(m, yI[ps.B,ps.S,ps.H] >= 0, container=Array)
    @variable(m, yO[ps.B,ps.S,ps.H] >= 0, container=Array)
    @variable(m, yL[ps.B,ps.S,ps.H] >= 0, container=Array)
    @variable(m, yS[ps.S,ps.H] >= 0, container=Array)
    @variable(m, x0[1:unc.nx], container=Array)
    @variable(m, x[1:unc.nx], container=Array)
    # */ ----------------------------------------------------------- /* #

    # */ -- obj function ------------------------------------------- /* #
    @objective(m, Min, c0 )
    # */ ----------------------------------------------------------- /* #

    # */ -- constraints -------------------------------------------- /* #
    @constraint(m, c01, c0 == exp10(-6)*sum(sum(sum((pp.cvg[g]+pp.cfg[g]/pp.ηg[g])*yG[g,s,h] for g in ps.G0) + pp.cvg[ps.gn]*yG[ps.gn,s,h] + pp.cs*yS[s,h] for h in ps.H)*pp.α for s in ps.S))
    @constraint(m, c02, ϕ[1] == exp10(-6)*sum(sum(sum(pp.eg[g]/pp.ηg[g]*yG[g,s,h] for g in ps.G) for h in ps.H)*pp.α for s in ps.S))
    @constraint(m, c03, ϕ[2] == exp10(-6)*sum(sum(1/pp.ηg[ps.gn]*yG[ps.gn,s,h] for h in ps.H)*pp.α for s in ps.S))
    @constraint(m, c04[b=ps.B,s=ps.S,h=ps.H], yL[b,s,h] - yL[b,s,circ(h,ps.H)] == pp.ηb[b]*yI[b,s,h] - yO[b,s,h])
    @constraint(m, c05[g=ps.G,s=ps.S,h=ps.H], yG[g,s,h] - yG[g,s,circ(h,ps.H)] <= pp.rg[g]*x0[g])
    @constraint(m, c06[g=ps.G,s=ps.S,h=ps.H], yG[g,s,circ(h,ps.H)] - yG[g,s,h] <= pp.rg[g]*x0[g])
    @constraint(m, c07[g=ps.G,s=ps.S,h=ps.H], yG[g,s,h] <= x0[g])
    @constraint(m, c08[b=ps.B,s=ps.S,h=ps.H], yI[b,s,h] <= pp.pb[b]*x0[b+ps.b0])
    @constraint(m, c09[b=ps.B,s=ps.S,h=ps.H], yO[b,s,h] <= pp.pb[b]*x0[b+ps.b0])
    @constraint(m, c10[b=ps.B,s=ps.S,h=ps.H], yL[b,s,h] <= x0[b])
    @constraint(m, c11, sum(yG[g,s,h]*pp.eg[g]/pp.ηg[g] for g in ps.G for s in ps.S for h in ps.H) <= pp.lco*x0[unc.nx0+1])
    @constraint(m, c12[s=ps.S,h=ps.H], sum(yG[g,s,h] for g in ps.G) + sum(yO[b,s,h]-yI[b,s,h] for b in ps.B) + yS[s,h] >= -x0[unc.nx0+2]*pp.pd[s,h] - sum(pp.pr[r,s,h]*x0[r+ps.r0] for r in ps.R))
    @constraint(m, c13[n=1:unc.nx0], x0[n]*exp10(-3) == x[n])
    @constraint(m, c14[n=1:unc.nh ], x0[unc.nx0+n]   == x[unc.nx0+n])
    # */ ----------------------------------------------------------- /* #
    optimize!(m)

    return m
end

@everywhere function LBoracle!(m::JuMP.Model,unc::u_type)::JuMP.Model

    # problem structure
    # min   ϕ + γᵀcᵢ
    # s.t.  ϕ + γᵀcₛ ≧ θₛ + λᵀₛ(x-xₛ),   ∀s
    #       x = xᵢ

    # */ -- variables ---------------------------------------------- /* #
    @variable(m, ϕ, container=Array)
    @variable(m, γ[1:unc.nc] >= .0, container=Array)
    @variable(m, x[1:unc.nx], container=Array)
    # */ ----------------------------------------------------------- /* #

    # */ -- obj function ------------------------------------------- /* #
    @objective(m, Min, 0. )
    # */ ----------------------------------------------------------- /* #

    # */ -- constraints -------------------------------------------- /* #
    # */ ----------------------------------------------------------- /* #
    optimize!(m)

    return m
end

@everywhere function UBoracle!(m::JuMP.Model,unc::u_type)::JuMP.Model

    # problem structure
    # max   ϕ + γᵀxᵢ
    # s.t.  ϕ + γᵀxₛ ≦ θₛ + ϕᵀₛ(c-cₛ),   ∀s
    #       c = cᵢ

    # */ -- variables ---------------------------------------------- /* #
    @variable(m, ϕ, container=Array)
    @variable(m, γ[1:unc.nx] <= .0, container=Array)
    @variable(m, c[1:unc.nc], container=Array)
    # */ ----------------------------------------------------------- /* #

    # */ -- obj function ------------------------------------------- /* #
    @objective(m, Max, 0. )
    # */ ----------------------------------------------------------- /* #

    # */ -- constraints -------------------------------------------- /* #
    # */ ----------------------------------------------------------- /* #
    optimize!(m)

    return m
end

# */ --------------------------------------------------------------------------------------------- /* #
# */ --- functions to generate structs : on main core -------------------------------------------- /* #
# */ --------------------------------------------------------------------------------------------- /* #

function gen_D__type(ms::ms_type,mp::mp_type,ps::ps_type,pp::pp_type,u::u_type,c::Int64,a::Int64,p::Int64,j::Int64,e::Float64)::D__type

    if nworkers()>1
        x = zeros(u.nh+u.nc,u.ni)
        for j in 1:u.nh
            uh = u.h[:,j]
            hmin,hmax = minimum(uh),maximum(uh)
            (hmax > hmin) ? x[j,:] .= (uh.-hmin)./(hmax-hmin) :  x[j,:] .= zeros(u.ni)
        end
        for j in 1:u.nc
            uc = u.c[:,j]
            cmin,cmax = minimum(uc),maximum(uc)
            (cmax > cmin) ? x[j+u.nh,:] .= (uc.-cmin)./(cmax-cmin) : x[j+u.nh,:] .= zeros(u.ni)
        end
        kass = kmeans(x,nworkers()).assignments
        kidx = Array{BitArray{1},1}(undef,0)
        for j in 1:nworkers()
            push!(kidx,kass.==j)
        end
        idx = [i for i in 1:u.ni]
        iidx = Array{Array{Int64,1},1}(undef,0)
        for j in 1:nworkers()
            push!(iidx,idx[kidx[j]])
        end
    else
        kidx = Array{BitArray{1},1}(undef,0)
        push!(kidx,zeros(Int64,u.ni).==0)
        iidx = Array{Array{Int64,1},1}(undef,0)
        push!(iidx,[i for i in 1:u.ni])
    end

    return D__type(ms,mp,ps,pp,u,c,a,p,kidx,iidx,j,e)
end

function gen_R__type(d::D__type)::R__type

    m = JuMP.Model(JuMP.optimizer_with_attributes(()->Gurobi.Optimizer(gurobi_env)))
    JuMP.set_optimizer_attribute(m,"OutputFlag",0)
    RMP!(m,d.ms,d.mp,d.unc)
    v = Array{JuMP.VariableRef,2}(undef,d.unc.nx+1,d.unc.ni)
    for i in 1:d.unc.ni v[:,i] .= vcat([m[:β][i]],m[:x][:,i]) end
    c = Array{JuMP.ConstraintRef,2}(undef,d.J,d.unc.ni)

    return R__type(m,v,c)
end

function gen_T1_type(d::D__type)::T1_type

    x = zeros(d.unc.nx,d.unc.ni)
    c = convert(Array{Float64,2},d.unc.c')
    β = zeros(d.unc.ni)
    θ = zeros(d.unc.ni)
    λ = zeros(d.unc.nx,d.unc.ni)
    re = Array{Tuple{Float64,Array{Float64,1}},1}(undef,d.unc.ni)

    return T1_type(x,c,β,θ,λ,re,0.)
end

function gen_T2_type(d::D__type)::T2_type

    x  = zeros(d.unc.nx,d.unc.ni)
    c  = convert(Array{Float64,2},d.unc.c')
    β  = zeros(d.unc.ni)
    θl = zeros(d.unc.ni)
    λl = zeros(d.unc.nx,d.unc.ni)
    θu = zeros(d.unc.ni)
    ϕu = zeros(d.unc.nc,d.unc.ni)
    re = Array{Tuple{Float64,Array{Float64,1},Array{Float64,1}},1}(undef,nworkers())
    ro = Array{Tuple{Float64,Float64,Array{Float64,1},Array{Float64,1}},1}(undef,d.unc.ni)
    Ie = zeros(Int64,0)

    return T2_type(x,c,β,θl,λl,θu,ϕu,re,ro,Ie,0.)
end

function gen_H__type(d::D__type)::H__type

    L = zeros(d.J)
    U = zeros(d.J)
    t = zeros(d.J)
    T = zeros(d.J,6)

    return H__type(L,U,t,T,0)
end

function gen_M__type(d::D__type)::M__type

    θ = zeros(nprocs())
    λ = zeros(nprocs(),d.unc.nx)
    ϕ = zeros(nprocs(),d.unc.nc)
    x = zeros(nprocs(),d.unc.nx)
    c = zeros(nprocs(),d.unc.nc)

    return M__type(0,θ,λ,ϕ,x,c)
end

function gen_B1_type(c::Int64,a::Int64,p::Int64,j::Int64,e::Float64)::B1_type

    ms,mp,ps,pp,unc = load_data(c)
    d = gen_D__type(ms,mp,ps,pp,unc,c,a,p,j,e)
    r = gen_R__type(d)
    t = gen_T1_type(d)
    h = gen_H__type(d)

    return B1_type(r,t,h,d)
end

function gen_B2_type(c::Int64,a::Int64,p::Int64,j::Int64,e::Float64)::B2_type

    ms,mp,ps,pp,unc = load_data(c)
    d = gen_D__type(ms,mp,ps,pp,unc,c,a,p,j,e)
    r = gen_R__type(d)
    t = gen_T2_type(d)
    h = gen_H__type(d)
    m = gen_M__type(d)

    return B2_type(r,t,h,m,d)
end

# */ --------------------------------------------------------------------------------------------- /* #
# */ --- functions to generate structs : everywhere ---------------------------------------------- /* #
# */ --------------------------------------------------------------------------------------------- /* #

@everywhere function gen_d__type(d::D__type)::d__type

    s = d.ps
    p = d.pp
    u = d.unc
    w = d.wrks
    j = d.J*nworkers()

    return d__type(s,p,u,w,j)
end

@everywhere function gen_Lb_type(d::d__type)::O__type

    m = JuMP.Model(JuMP.optimizer_with_attributes(()->Gurobi.Optimizer(gurobi_env)))
    JuMP.set_optimizer_attribute(m,"OutputFlag",0)
    JuMP.set_optimizer_attribute(m,"Threads",1)
    LBoracle!(m,d.unc)
    v = vcat(m[:ϕ],m[:γ])
    c = Array{JuMP.ConstraintRef,1}(undef,d.Jw+1)
    o = v'*ones(length(v))
    e = v'*ones(length(v))
    h = zeros(d.Jw+1)

    return O__type(m,c,v,o,e,h,0)
end

@everywhere function gen_Ub_type(d::d__type)::O__type

    m = JuMP.Model(JuMP.optimizer_with_attributes(()->Gurobi.Optimizer(gurobi_env)))
    JuMP.set_optimizer_attribute(m,"OutputFlag",0)
    JuMP.set_optimizer_attribute(m,"Threads",1)
    UBoracle!(m,d.unc)
    v = vcat(m[:ϕ],m[:γ])
    c = Array{JuMP.ConstraintRef,1}(undef,d.Jw+1)
    o = v'*ones(length(v))
    e = v'*ones(length(v))
    h = zeros(d.Jw+1)

    return O__type(m,c,v,o,e,h,0)
end

@everywhere function gen_E__type(d::d__type)::E__type

    m = JuMP.Model(JuMP.optimizer_with_attributes(()->Gurobi.Optimizer(gurobi_env)))
    JuMP.set_optimizer_attribute(m,"OutputFlag",0)
    JuMP.set_optimizer_attribute(m,"Threads",1)
    JuMP.set_optimizer_attribute(m,"Method",2)
    SP!(m,d.ps,d.pp,d.unc)
    v = vcat(m[:c0],m[:ϕ][:])
    o = v'*ones(length(v))

    return E__type(m,v,o)
end

@everywhere function gen_M__type(d::d__type)::M__type

    θ = zeros(d.Jw+1)
    λ = zeros(d.Jw+1,d.unc.nx)
    ϕ = zeros(d.Jw+1,d.unc.nc)
    x = zeros(d.Jw+1,d.unc.nx)
    c = zeros(d.Jw+1,d.unc.nc)

    return M__type(0,θ,λ,ϕ,x,c)
end

@everywhere function gen_t__type(d::d__type)::t__type

    x  = zeros(d.unc.nx,d.unc.ni)
    c  = convert(Array{Float64,2},d.unc.c')

   return t__type(x,c)
end

@everywhere function gen_S1_type(D::D__type)::S1_type

    d = gen_d__type(D)
    e = gen_E__type(d)
    t = gen_t__type(d)

    return S1_type(e,d,t)
end

@everywhere function gen_S2_type(D::D__type)::S2_type

    d = gen_d__type(D)
    e = gen_E__type(d)
    l = gen_Lb_type(d)
    u = gen_Ub_type(d)
    m = gen_M__type(d)
    t = gen_t__type(d)

    return S2_type(e,l,u,m,d,t)
end

# */ --------------------------------------------------------------------------------------------- /* #
# */ --- functions to generate structs : on main core and everywhere ----------------------------- /* #
# */ --------------------------------------------------------------------------------------------- /* #

function generate_structures(c::Int64,a::Int64,p::Int64,j::Int64,e::Float64)

    if a == 1
        B = gen_B1_type(c,a,p,j,e)
        sendto(procs(),d=B.data)
        @everywhere S = gen_S1_type(d)
    end
    if a == 2
        B = gen_B2_type(c,a,p,j,e)
        sendto(procs(),d=B.data)
        @everywhere S = gen_S2_type(d)
    end

    return B,S
end

# */ --------------------------------------------------------------------------------------------- /* #
# */ --- lower level functions : on main core ---------------------------------------------------- /* #
# */ --------------------------------------------------------------------------------------------- /* #

function set_iex!(b::B2_type)::B2_type

    b.temp.Ie = zeros(Int64,0)
    for p in 1:nworkers()
        ie = findmax(b.data.mp.π[b.data.kidx[p]].*(b.temp.θu[b.data.kidx[p]].-b.temp.θl[b.data.kidx[p]]))[2]
        push!(b.temp.Ie,b.data.iidx[p][ie])
    end

    return b
end

# */ --------------------------------------------------------------------------------------------- /* #
# */ --- lower level functions : everywhere ------------------------------------------------------ /* #
# */ --------------------------------------------------------------------------------------------- /* #

@everywhere function solv0_exact(i::Int64)::Tuple{Float64,Array{Float64,1}}

    S.ex.objf = S.ex.vars'*vcat([1.],S.temp.c[:,i])
    @objective(S.ex.m, Min, S.ex.objf)
    fix.(S.ex.m[:x],S.temp.x[:,i];force=true)
    optimize!(S.ex.m)
    θ = objective_value(S.ex.m)
    λ = dual.(FixRef.(S.ex.m[:x]))

    return θ,λ
end

@everywhere function solv1_exact(i::Int64)::Tuple{Float64,Array{Float64,1},Array{Float64,1}}

    S.ex.objf = S.ex.vars'*vcat([1.],S.temp.c[:,i])
    @objective(S.ex.m, Min, S.ex.objf)
    fix.(S.ex.m[:x],S.temp.x[:,i];force=true)
    optimize!(S.ex.m)
    θl = objective_value(S.ex.m)
    λl = dual.(FixRef.(S.ex.m[:x]))
    ϕu = value.(S.ex.m[:ϕ])

    return θl,λl,ϕu
end

@everywhere function update_m_!(s::S2_type)::S2_type

    s.m.θ[s.m.n+1:s.m.n+m0.n]   .= m0.θ[1:m0.n]
    s.m.λ[s.m.n+1:s.m.n+m0.n,:] .= m0.λ[1:m0.n,:]
    s.m.ϕ[s.m.n+1:s.m.n+m0.n,:] .= m0.ϕ[1:m0.n,:]
    s.m.x[s.m.n+1:s.m.n+m0.n,:] .= m0.x[1:m0.n,:]
    s.m.c[s.m.n+1:s.m.n+m0.n,:] .= m0.c[1:m0.n,:]
    s.m.n += m0.n

    return s
end

@everywhere function update_lb!(s::S2_type)::S2_type

    for n in s.olb.n+1:s.m.n
        s.olb.cons[n] = @constraint(s.olb.m, s.olb.m[:ϕ] + s.olb.m[:γ]'*s.m.c[n,:] >= 0. )
    end
    s.olb.help[1:s.m.n] .= s.m.θ[1:s.m.n] .- sum!(ones(s.m.n),s.m.λ[1:s.m.n,:].*s.m.x[1:s.m.n,:])
    s.olb.n = s.m.n

    return s
end

@everywhere function update_ub!(s::S2_type)::S2_type

    for n in s.oub.n+1:s.m.n
        s.oub.cons[n] = @constraint(s.oub.m, s.oub.m[:ϕ] + s.oub.m[:γ]'*s.m.x[n,:] <= 0. )
    end
    s.oub.help[1:s.m.n] .= s.m.θ[1:s.m.n] .- sum!(ones(s.m.n),s.m.ϕ[1:s.m.n,:].*s.m.c[1:s.m.n,:])
    s.oub.n = s.m.n

    return s
end

@everywhere function update_s!(s::S2_type)::S2_type

    update_m_!(s)
    update_lb!(s)
    update_ub!(s)

    return s
end

@everywhere function run_lb(i::Int64)::Tuple{Float64,Array{Float64,1}}

    S.olb.objf = S.olb.vars'*vcat([1.],S.temp.c[:,i])
    @objective(S.olb.m, Min, S.olb.objf)
    set_normalized_rhs.(S.olb.cons[1:S.olb.n], S.olb.help[1:S.olb.n].+ S.m.λ[1:S.olb.n,:]*S.temp.x[:,i])
    optimize!(S.olb.m)
    θl = objective_value(S.olb.m)
    λl = S.m.λ[1:S.m.n,:]'*dual.(S.olb.cons[1:S.olb.n])

    return θl,λl
end

@everywhere function run_ub(i::Int64)::Tuple{Float64,Array{Float64,1}}

    S.oub.objf = S.oub.vars'*vcat([1.],S.temp.x[:,i])
    @objective(S.oub.m, Max, S.oub.objf)
    set_normalized_rhs.(S.oub.cons[1:S.oub.n], S.oub.help[1:S.oub.n].+ S.m.ϕ[1:S.oub.n,:]*S.temp.c[:,i])
    optimize!(S.oub.m)
    θu = objective_value(S.oub.m)
    ϕu = S.m.ϕ[1:S.m.n,:]'*dual.(S.oub.cons[1:S.oub.n])

    return θu,ϕu
end

@everywhere function run_oracles(i::Int64)::Tuple{Float64,Float64,Array{Float64,1},Array{Float64,1}}

    θl,λl = run_lb(i)
    θu,ϕu = run_ub(i)

    return θl,θu,λl,ϕu
end

# */ --------------------------------------------------------------------------------------------- /* #

function step_a!(b::B1_type)::B1_type

    optimize!(b.rmp.m)
    b.temp.β .= value.(b.rmp.m[:β])
    b.temp.x .= value.(b.rmp.m[:x])
    for w in workers()
        @spawnat(w, Core.eval(Main, Expr(:(.=), :(S.temp.x),b.temp.x)))
    end

    return b
end

function step_a!(b::B2_type)::B2_type

    optimize!(b.rmp.m)
    b.temp.β .= value.(b.rmp.m[:β])
    b.temp.x .= value.(b.rmp.m[:x])
    for w in workers()
        @spawnat(w, Core.eval(Main, Expr(:(.=), :(S.temp.x),b.temp.x)))
    end

    return b
end

function step_b!(b::B1_type)::B1_type

    b.hist.L[b.hist.k] = exp10(6)*objective_value(b.rmp.m)

    return b
end

function step_b!(b::B2_type)::B2_type

    b.hist.L[b.hist.k] = exp10(6)*objective_value(b.rmp.m)

    return b
end

function step_c!(b::B1_type)::B1_type

    b.temp.re .= pmap(solv0_exact,1:b.data.unc.ni)
    for i in 1:b.data.unc.ni
        b.temp.θ[i] = b.temp.re[i][1]
        b.temp.λ[:,i] .= b.temp.re[i][2]
    end

    return b
end

function step_c!(b::B2_type)::B2_type

    set_iex!(b)
    b.temp.re .= pmap(solv1_exact,b.temp.Ie)
    for w in 1:nworkers()
        b.m0.θ[w]    = b.temp.re[w][1]
        b.m0.λ[w,:] .= b.temp.re[w][2]
        b.m0.ϕ[w,:] .= b.temp.re[w][3]
        b.m0.x[w,:] .= b.temp.x[:,b.temp.Ie[w]]
        b.m0.c[w,:] .= b.temp.c[:,b.temp.Ie[w]]
    end
    sendto(procs(),m0=b.m0)
    @everywhere update_s!(S)

    return b
end

function step_d!(b::B1_type)::B1_type

    return b
end

function step_d!(b::B2_type)::B2_type

    b.temp.ro .= pmap(run_oracles,1:b.data.unc.ni)
    for i in 1:b.data.unc.ni
        b.temp.θl[i]    = b.temp.ro[i][1]
        b.temp.θu[i]    = b.temp.ro[i][2]
        b.temp.λl[:,i] .= b.temp.ro[i][3]
        b.temp.ϕu[:,i] .= b.temp.ro[i][4]
    end

    return b
end

function step_e!(b::B1_type)::B1_type

    if (b.hist.k==1)
        b.hist.U[b.hist.k] = exp10(6)*(value(b.rmp.m[:f])+b.data.mp.κ*b.data.mp.π'*b.temp.θ)
    else
        b.hist.U[b.hist.k] = min(b.hist.U[b.hist.k-1],exp10(6)*(value(b.rmp.m[:f])+b.data.mp.κ*b.data.mp.π'*b.temp.θ))
    end

    return b
end

function step_e!(b::B2_type)::B2_type

    if (b.hist.k==1)
        b.hist.U[b.hist.k] = exp10(6)*(value(b.rmp.m[:f])+b.data.mp.κ*b.data.mp.π'*b.temp.θu)
    else
        b.hist.U[b.hist.k] = min(b.hist.U[b.hist.k-1],exp10(6)*(value(b.rmp.m[:f])+b.data.mp.κ*b.data.mp.π'*b.temp.θu))
    end

    return b
end

function step_f!(b::B1_type)::B1_type

    b.temp.Δ = (b.hist.U[b.hist.k]-b.hist.L[b.hist.k])/b.hist.U[b.hist.k]*100
    if b.temp.Δ > b.data.δ
        for i in 1:b.data.unc.ni
            @constraint(b.rmp.m,b.rmp.m[:β][i] >= b.temp.θ[i]+b.temp.λ[:,i]'*(b.rmp.m[:x][:,i].-b.temp.x[:,i]))
        end
    end

    return b
end

function step_f!(b::B2_type)::B2_type

    b.temp.Δ = (b.hist.U[b.hist.k]-b.hist.L[b.hist.k])/b.hist.U[b.hist.k]*100
    if b.temp.Δ > b.data.δ
        for i in 1:b.data.unc.ni
            @constraint(b.rmp.m,b.rmp.m[:β][i] >= b.temp.θl[i]+b.temp.λl[:,i]'*(b.rmp.m[:x][:,i].-b.temp.x[:,i]))
        end
    end

    return b
end

# */ --------------------------------------------------------------------------------------------- /* #

function step_0!(b::B2_type)::B2_type

    S.temp.x[:,1] .= round.(vcat(minimum(b.data.mp.xh,dims=2)[:]*.99,minimum(b.data.unc.h,dims=1)[:].*[.99,1.01]); digits=4)
    S.temp.c[:,1] .= round.(minimum(b.data.unc.c,dims=1)[:]*.99; digits=4)
    r = solv1_exact(1)
    b.m0.n = 1
    b.m0.θ[b.m0.n] = r[1]
    b.m0.λ[b.m0.n,:] .= r[2]
    b.m0.ϕ[b.m0.n,:] .= r[3]
    b.m0.x[b.m0.n,:] .= S.temp.x[:,1]
    b.m0.c[b.m0.n,:] .= S.temp.c[:,1]
    sendto(procs(),m0=b.m0)
    @everywhere update_s!(S)
    for w in workers()
        @spawnat(w, Core.eval(Main, Expr(:(.=), :(S.temp.x),vcat(b.data.mp.xh,b.data.unc.h'))))
        @spawnat(w, Core.eval(Main, Expr(:(.=), :(S.temp.c),b.data.unc.c')))
    end
    b.m0.n = nworkers()
    b.temp.ro .= pmap(run_oracles,1:b.data.unc.ni)
    for i in 1:b.data.unc.ni
        b.temp.θl[i]    = b.temp.ro[i][1]
        b.temp.θu[i]    = b.temp.ro[i][2]
        b.temp.λl[:,i] .= b.temp.ro[i][3]
        b.temp.ϕu[:,i] .= b.temp.ro[i][4]
    end

    return b
end

function step!(b::B1_type)::B1_type

    b.hist.k += 1
    b.hist.T[b.hist.k,1] = @elapsed step_a!(b)
    b.hist.T[b.hist.k,2] = @elapsed step_b!(b)
    b.hist.T[b.hist.k,3] = @elapsed step_c!(b)
    b.hist.T[b.hist.k,4] = @elapsed step_d!(b)
    b.hist.T[b.hist.k,5] = @elapsed step_e!(b)
    b.hist.T[b.hist.k,6] = @elapsed step_f!(b)

    return b
end

function step!(b::B2_type)::B2_type

    b.hist.k += 1
    b.hist.T[b.hist.k,1] = @elapsed step_a!(b)
    b.hist.T[b.hist.k,2] = @elapsed step_b!(b)
    b.hist.T[b.hist.k,3] = @elapsed step_c!(b)
    b.hist.T[b.hist.k,4] = @elapsed step_d!(b)
    b.hist.T[b.hist.k,5] = @elapsed step_e!(b)
    b.hist.T[b.hist.k,6] = @elapsed step_f!(b)

    return b
end

# */ --------------------------------------------------------------------------------------------- /* #

function solve_Benders!(b::B1_type;print_output::Int64=1)::B1_type

    print_output == 1 ? print0(b) : nothing
    while (b.hist.k < b.data.J)
        step!(b)
        print_output == 1 ? print_info(b) : nothing
        (b.temp.Δ <= b.data.δ) ? break : nothing
    end
    print_output == 1 ? print_summary(b) : nothing

    return b
end

function solve_Benders!(b::B2_type;print_output::Int64=1)::B2_type

    print_output == 1 ? print0(b) : nothing
    step_0!(b)
    while (b.hist.k < b.data.J)
        step!(b)
        print_output == 1 ? print_info(b) : nothing
        (b.temp.Δ <= b.data.δ) ? break : nothing
    end
    print_output == 1 ? print_summary(b) : nothing

    return b
end

# */ --------------------------------------------------------------------------------------------- /* #


text:
using ImageContrastAdjustment
using Test, ImageCore, ColorTypes, FixedPointNumbers, ImageFiltering, TestImages

@testset "ImageContrastAdjustment.jl" begin
    include("core.jl")
    include("histogram_construction.jl")
    include("histogram_matching.jl")
    include("histogram_equalization.jl")
    include("histogram_midway_equalization.jl")
    include("gamma_adjustment.jl")
    include("linear_stretching.jl")
    include("contrast_stretching.jl")
end


text:
function __init__()
    settings.log_output = stdout;
    atexit(_closelog)
end


text:
using ForwardDiff
using Plots

#-----------------------------------------------------------------------------
# ---- Cross Sections
#-----------------------------------------------------------------------------

model = KineticMixing(1e2, 1e3, 1.0, 1e-3)

cross_section_xx_dd(4e4, model)
cross_section_xx_ss(4e4, model)
cross_section_xx_bb(4e4, model)
cross_section_xx_ee(4e4, model)
cross_section_xx_μμ(4e4, model)
cross_section_xx_ττ(4e4, model)
cross_section_xx_νν(4e4, model)
cross_section_xx_uu(4e4, model)
cross_section_xx_cc(4e4, model)
cross_section_xx_tt(4e4, model)
cross_section_xx_vv(4e4, model)

@code_warntype dm_annihilation_cross_section(4e2, model)

model.width_v

begin
    model = KineticMixing(350.0, 1e3, 1.0, 1e-2)
    cmes = 10 .^ LinRange(log10(2.001model.mx), log10(4000), 1000)
    cs = [dm_annihilation_cross_section(cme, model) for cme in cmes]

    plot(cmes, cs, xaxis = :log, yaxis = :log)
end

#-----------------------------------------------------------------------------
# ---- Thermal Cross Sections ------------------------------------------------
#-----------------------------------------------------------------------------

begin
    model = KineticMixing(9e2, 1e3, 1.0, 1e-3)
    xs = 10 .^ LinRange(-1.0, 3.0, 500)
    tcs = [dm_thermal_cross_section(x, model) for x in xs]

    plot(xs, tcs, xaxis = :log, yaxis = :log)
end

#-----------------------------------------------------------------------------
# ---- other ------------------------------------------------
#-----------------------------------------------------------------------------

begin
    Ts = 10 .^ LinRange(-5.0, 5.0, 200)
    hs = sm_heff.(Ts)
    gs = sm_geff.(Ts)
    gstars = sm_sqrt_gstar.(Ts)

    plot(Ts, hs, label = "heff")
    plot!(Ts, gs, label = "geff")
    plot!(Ts, gstars, label = "√g*")
    xaxis!(:log)
end

#-----------------------------------------------------------------------------
# ---- other ------------------------------------------------
#-----------------------------------------------------------------------------

model = KineticMixing(1e2, 1e3, 1.0, 1e-3)

relic_density(model, MPU())
relic_density(model, ODE(reltol = 1e-12))

#-----------------------------------------------------------------------------
# ---- other ------------------------------------------------
#-----------------------------------------------------------------------------

begin
    model = KineticMixing(3e2, 1e3, 1.0, 1e-3)

    xs = LinRange(1e-3, 50.0, 1000)
    fs = [compute_f(x, model) for x in xs]
    qs = [compute_q(x, model) for x in xs]
    ps = [compute_p(x, model) for x in xs]

    plot(xs, qs, label = "λQ(x)")
    plot!(xs, abs.(ps), label = "|P(x)|")
    plot!(xs, fs * 8e-9, label = "λf(x)")
    xlabel!("x")
    yaxis!(:log)
end

max_stars_repo_name: LoganAMorrison/AsymptoticBoltzmann

max_issues_repo_name: LoganAMorrison/AsymptoticBoltzmann

max_forks_repo_name: LoganAMorrison/AsymptoticBoltzmann


text:
function step_through_env(env::NatureEnv, policy::T) where {T <: MultiPPOManager}
    RLBase.reset!(env)

    get_player_poses() = [p.pos for p in env.players if !p.dead]
    function get_food_poses()
        fps = Vector()
        for f in 1:env.food_types
            push!(fps, Vector())
            for pos in CartesianIndices(env.food_frames[f])
                if env.food_frames[f][pos] > 0f0
                    push!(fps[f], (Tuple(pos)..., env.food_frames[f][pos]))
                end
            end
        end
        fps
    end

    player_poses = []
    food_poses = []
    food_color_maps = (
        (:darkblue, :lightblue),
        (:yellow, :orange))
    nframes = env.episode_len
    frame = 0

    while !is_terminated(env) && frame < nframes
        env |> policy |> env
        push!(player_poses, get_player_poses())
        push!(food_poses, get_food_poses())
        frame += 1
    end

    fig = Figure()
    ax  = Axis(fig[1,1])
    limits!(ax, 0, env.world_size[1], 0, env.world_size[2])
    lsgrid = labelslidergrid!(fig, ["step"], [1:1:frame])
    fig[2, 1] = lsgrid.layout
    slider = lsgrid.sliders[1]

    for f in 1:env.food_types
        fps = @lift [(x, y) for (x, y, _) in food_poses[$(slider.value)][f]]
        fvs = @lift [z for (_, _, z) in food_poses[$(slider.value)][f]]
        GLMakie.scatter!(fps,color=fvs, colormap=food_color_maps[f], markersize=40, marker=:rect)
    end
    pps = @lift player_poses[$(slider.value)]
    GLMakie.scatter!(pps, color=:pink, markersize=60)

    display(fig)
    nothing
end


text:
using KalmanAD

using Test

include("test_param.jl")

@testset "statespacemodel" begin
    include("statespacemodel_test.jl")
end

@testset "observer" begin
    include("observer_test.jl")
end

@testset "filter" begin
    include("filter_test.jl")
end


text:
##
# routines to give a rigorous solution of the real linear system Ax=B

"""
    epsilon_inflation(A::AbstractMatrix{T}, b::AbstractArray{S, N};
                      r=0.1, ϵ=1e-20, iter_max=20) where {T<:Real, S<:Real, N}

    epsilon_inflation(A::AbstractMatrix{T}, b::AbstractArray{S, N};
                      r=0.1, ϵ=1e-20, iter_max=20) where {T<:Interval, S<:Interval, N}

Gives an enclosure of the solution of the square linear system ``Ax=b``
using the ϵ-inflation algorithm,  see algorithm 10.7 of [[RUM10]](@ref)

### Input

* `A`        -- square matrix of size n × n
* `b`        -- vector of length n or matrix of size n × m
* `r`        -- relative inflation, default 10%
* `ϵ`        -- absolute inflation, default 1e-20
* `iter_max` -- maximum number of iterations

### Output

* `x`    -- enclosure of the solution of the linear system
* `cert` -- Boolean flag, if `cert==true`, then `x` is *certified* to contain the true
solution of the linear system, if `cert==false`, then the algorithm could not prove that ``x``
actually contains the true solution.

### Algorithm

Given the real system ``Ax=b`` and an approximate solution ``̃x``, we initialize
``x₀ = [̃x, ̃x]``. At each iteration the algorithm computes the inflation

``y = xₖ * [1 - r, 1 + r] .+ [-ϵ, ϵ]``

and the update

``xₖ₊₁ = Z + (I - CA)y``,

where ``Z = C(b - Ax₀)`` and ``C`` is an approximate inverse of ``A``. If the condition
``xₖ₊₁ ⊂ y `` is met, then ``xₖ₊₁`` is a proved enclosure of ``A⁻¹b`` and `cert` is set to
true. If the condition is not met by the maximum number of iterations, the
latest computed enclosure is returned, but ``cert`` is set to false, meaning the algorithm
could not prove that the enclosure contains the true solution. For interval systems,
``̃x`` is obtained considering the midpoint of ``A`` and ``b``.

### Notes

- This algorithm is meant for *real* linear systems, or interval systems with
very tiny intervals. For interval linear systems with wider intervals, see the
[`solve`](@ref) function.

### Examples

```jldoctest
julia> A = [1 2;3 4]
2×2 Matrix{Int64}:
 1  2
 3  4

julia> b = A * ones(2)
2-element Vector{Float64}:
 3.0
 7.0

julia> x, cert = epsilon_inflation(A, b)
(Interval{Float64}[[0.999999, 1.00001], [0.999999, 1.00001]], true)

julia> ones(2) .∈ x
2-element BitVector:
 1
 1

julia> cert
true
```
"""
function epsilon_inflation(A::AbstractMatrix{T}, b::AbstractArray{S, N};
                           r=0.1, ϵ=1e-20, iter_max=20) where {T<:Real, S<:Real, N}

    r1 = Interval(1 - r, 1 + r)
    ϵ1 = Interval(-ϵ, ϵ)
    R = inv(mid.(A))
    C = I - R * A
    xs = R * mid.(b)
    z = R * (b - (A * Interval.(xs)))
    x = z

    for _ in 1:iter_max
        y = r1 * x .+ ϵ1
        x = z + C * y
        if all(x .⊂ y)
            return xs + x, true
        end
    end
    return xs + x, false
end


text:
function parse_kw(ps::ParseState, ::Type{Val{Tokens.FOR}})
    startbyte = ps.t.startbyte
    start_col = ps.t.startpos[2] + 4

    # Parsing
    kw = INSTANCE(ps)
    format_kw(ps)
    @catcherror ps startbyte ranges = @default ps parse_ranges(ps)
    block = EXPR{Block}(EXPR[], 0, Variable[], "")
    @catcherror ps startbyte @default ps parse_block(ps, block, start_col)
    next(ps)
    ret = EXPR{For}(EXPR[kw, ranges, block, INSTANCE(ps)], ps.nt.startbyte - startbyte, Variable[], "")

    return ret
end

function parse_ranges(ps::ParseState)
    startbyte = ps.nt.startbyte
    defs = []
    arg = @closer ps range @closer ps comma @closer ps ws parse_expression(ps)
    _track_range_assignment(ps, arg)
    if ps.nt.kind == Tokens.COMMA
        arg = EXPR{Block}(EXPR[arg], arg.span, Variable[], "")
        while ps.nt.kind == Tokens.COMMA
            next(ps)
            push!(arg.args, INSTANCE(ps))
            format_comma(ps)

            arg.span += last(arg.args).span
            @catcherror ps startbyte nextarg = @closer ps comma @closer ps ws parse_expression(ps)
            _track_range_assignment(ps, nextarg)
            push!(arg.args, nextarg)
            arg.span += last(arg.args).span
        end
    end
    return arg
end

function _track_range_assignment(ps, x) end
function _track_range_assignment(ps, x::EXPR{BinaryOpCall})
    if ((x.args[2] isa EXPR{OPERATOR{ComparisonOp,Tokens.IN,false}} || x.args[2] isa EXPR{OPERATOR{ComparisonOp,Tokens.ELEMENT_OF,false}}))
        append!(x.defs, _track_assignment(ps::ParseState, x.args[1], x.args[3]))
    end
end

function parse_kw(ps::ParseState, ::Type{Val{Tokens.WHILE}})
    startbyte = ps.t.startbyte
    start_col = ps.t.startpos[2] + 4

    # Parsing
    kw = INSTANCE(ps)
    format_kw(ps)
    @catcherror ps startbyte cond = @default ps @closer ps ws parse_expression(ps)
    block = EXPR{Block}(EXPR[], 0, Variable[], "")
    @catcherror ps startbyte @default ps parse_block(ps, block, start_col)
    next(ps)

    ret = EXPR{While}(EXPR[kw, cond, block, INSTANCE(ps)], ps.nt.startbyte - startbyte, Variable[], "")

    return ret
end

function parse_kw(ps::ParseState, ::Type{Val{Tokens.BREAK}})
    return EXPR{Break}(EXPR[INSTANCE(ps)], ps.nt.startbyte - ps.t.startbyte, Variable[], "")
end

function parse_kw(ps::ParseState, ::Type{Val{Tokens.CONTINUE}})
    return EXPR{Continue}(EXPR[INSTANCE(ps)], ps.nt.startbyte - ps.t.startbyte, Variable[], "")
end

"""
parse_generator(ps)

Having hit `for` not at the beginning of an expression return a generator.
Comprehensions are parsed as SQUAREs containing a generator.
"""
function parse_generator(ps::ParseState, ret)
    startbyte = ps.nt.startbyte
    next(ps)
    kw = INSTANCE(ps)
    ret = EXPR{Generator}(EXPR[ret, kw], ret.span - startbyte, Variable[], "")
    @catcherror ps startbyte ranges = @closer ps paren @closer ps square parse_ranges(ps)

    if ps.nt.kind == Tokens.IF
        if ranges isa EXPR{Block}
            ranges = EXPR{Filter}(EXPR[ranges.args...], ranges.span, Variable[], "")
        else
            ranges = EXPR{Filter}(EXPR[ranges], ranges.span, Variable[], "")
        end
        next(ps)
        unshift!(ranges.args, INSTANCE(ps))
        @catcherror ps startbyte cond = @closer ps paren parse_expression(ps)
        unshift!(ranges.args, cond)
        ranges.span = sum(a.span for a in ranges.args)
        push!(ret.args, ranges)
    else
        if ranges isa EXPR{Block}
            append!(ret.args, ranges.args)
        else
            push!(ret.args, ranges)
        end
    end
    ret.span += ps.nt.startbyte

    # This should reverse order of iterators
    if ret.args[1] isa EXPR{Generator} || ret.args[1] isa EXPR{Flatten}
        ret = EXPR{Flatten}([ret], ret.span, Variable[], "")
    end

    return ret
end


text:
module WaveEquation
include("discretization.jl")
using LinearAlgebra
disc=Discretization

names = ["p", "v"]
WAVE_NUMBER = 8 * pi
LEFT = -1
RIGHT = 1

function A(n)
        Aout = zeros(2,2)
        Aout[1,2] = n
        Aout[2,1] = n
        return Aout
end

function A()
        A = zeros(2,2)
        A[1,2] = 1
        A[2,1] = 1
        return A
end

function reflection_matrix(r)
        R1 = zeros(1,1)
        R1[1,1] = r
        return R1
end

function projection(a, dr::Float64)
        dr_ = ones(1,1) * dr
        r = -ones(1,1)
        P = projection_matrix(A(a), r, dr_)
        return P
end

function sat_left(ops::disc.OperatorsType, P::Matrix)
        B = boundary(sizes(ops.nx), [1])
        I2 = Matrix(I, 2, 2)
        sat_term_data = kron(I2, ops.Hxi) * B' * (A(LEFT) * P)
        sat_term = sat_term_data * B
        return sat_term, sat_term_data
end

function sat_right(ops::disc.OperatorsType, P::Matrix)
        B = boundary(sizes(ops.nx), [ops.nx])
        I2 = Matrix(I, 2, 2)
        sat_term_data = kron(I2, ops.Hxi) * B' * (A(RIGHT) * P)
        sat_term = sat_term_data * B
        return sat_term, sat_term_data
end

function sat(A, side, ops, P)
        if side == LEFT
                B = disc.boundary(sizes(ops.nx), [1])
        else
                B = disc.boundary(sizes(ops.nx), [ops.nx])
        end
        I_ = Matrix(I, 2, 2)
        return disc.sat(I_, A, B, ops, P)
end

function zero_p_left()
        return reflection_matrix(-1)
end

function zero_v_right()
        return reflection_matrix(-1)
end

function non_reflective()
        return zeros(1,1)
end

function leakage(alpha)
        return ones(1,1) * (alpha - 1) / (1 + alpha)
end

function naive_leakage(alpha, ops)
        B = disc.boundary(sizes(ops.nx), [ops.nx])
        P = zeros(2, 2)
        P[2,1] = 1
        P[2,2] = -alpha
        I_ = Matrix(I, 2, 2)
        s1, s2 = disc.sat(I_, I_, B, ops, P)
        return P, s1, s2
end

function energy_norm(Hx)
        return kron(Matrix(I, 2, 2), Hx)
end

function manufactured_solution(x, t)
        k = WAVE_NUMBER
        p = t -> cos.(k * t) .* sin.(k * x)
        v = t -> - sin.(k * t) .* cos.(k * x)
        q = Dict("p" => p(t), "v" => v(t))
        return q
end

function functionals(t)
        k = WAVE_NUMBER
        Fp = (1 - cos(k)) / k * cos(k * t)
        Fv = - sin(k) / k * sin(k * t)
        return Dict("p" => Fp, "v" => Fv)
end

function sizes(nx)
        return disc.variables(names, [nx, nx])
end

function restrict_left(nx)
        left = disc.variables(names, [1, 1])
        block_left = disc.block_diag(restrict(vars, left))
end

end

max_stars_repo_name: ooreilly/sbp

max_issues_repo_name: ooreilly/sbp

max_forks_repo_name: ooreilly/sbp


text:
#implementation of 2-opt

function two_opt{T<:Real}(distmat::Matrix{T}, path::Vector{Int})
  return _two_opt(distmat,path,false)
end

function two_opt!{T<:Real}(distmat::Matrix{T}, path::Vector{Int})
  return _two_opt(distmat,path,true)
end

function _two_opt{T<:Real}(distmat::Matrix{T}, path::Vector{Int}, inplace::Bool)
	# size checks
	n = length(path)
	if size(distmat, 1) != size(distmat, 2)
		error("Distance matrix passed to two_opt must be square.")
	end

	# in place or no?
  path = inplace ? path : copy(path)

	# main loop
	# check every possible switch until no 2-swaps reduce objective
	# if the path passed in is a loop (first/last nodes are the same)
	# then we must keep these the endpoints of the path the same
	# ie just keep it a loop, and therefore it doesn't matter which node is at the end
	# if the path is not a cycle, we can do any reversing we like
	isCycle = path[1] == path[end]
	switchLow = isCycle ? 2 : 1
	switchHigh = isCycle ? n - 1 : n
	prevCost = Inf
	curCost = pathcost(distmat, path)
	while prevCost > pathcost(distmat, path)
		prevCost = curCost
		# we can't change the first
		for i in switchLow:(switchHigh-1)
			for j in (i+1):switchHigh
				altCost = pathcost_rev(distmat, path, i, j)
				if altCost < curCost
					curCost = altCost
					reverse!(path, i, j)
				end
			end
		end
	end
	return path
end

function pathcost{T<:Real}(distmat::Matrix{T}, path::Vector{Int}, lb::Int = 1, ub::Int = length(path))
	cost = zero(T)
	for i in lb:(ub - 1)
		@inbounds cost += distmat[path[i], path[i+1]]
	end
	return cost
end

function pathcost_rev{T<:Real}(distmat::Matrix{T}, path::Vector{Int}, revLow::Int, revHigh::Int)
	cost = zero(T)
	# if there's an initial unreversed section
	if revLow > 1
		for i in 1:(revLow - 2)
			@inbounds cost += distmat[path[i], path[i+1]]
		end
		# from end of unreversed section to beginning of reversed section
		@inbounds cost += distmat[path[revLow - 1], path[revHigh]]
	end
	# main reverse section
	for i in revHigh:-1:(revLow + 1)
		@inbounds cost += distmat[path[i], path[i-1]]
	end
	# if there's an unreversed section after the reversed bit
	n = length(path)
	if revHigh < length(path)
		# from end of reversed section back to regular
		@inbounds cost += distmat[path[revLow], path[revHigh + 1]]
		for i in (revHigh + 1):(n-1)
			@inbounds cost += distmat[path[i], path[i+1]]
		end
	end
	return cost
end


text:
using MLStyle

@testset "uncomp map" begin
    @test [2, 3, 5] == @match [(1, 2), (2, 3), (3, 5)] begin
        [(_, snd) for snd in seq] => seq
        _ => nothing
    end
end

@testset "uncomp filter" begin
    @test [3, 5] == @match [(1, 2), (2, 3), (3, 5)] begin
        [(_, snd) for snd in seq if snd > 2] => seq
        _ => nothing
    end
end


text:
struct APIDispatcher
    global_table::Dict{Symbol,Ptr{Cvoid}}
    instance_table::Dict{Symbol,Ptr{Cvoid}}
    device_tables::Dict{Device,Dict{Symbol,Ptr{Cvoid}}}
    lock::Threads.ReentrantLock
end

APIDispatcher() = APIDispatcher(Dict(), Dict(), Dict(), ReentrantLock())

table(disp::APIDispatcher, ::Nothing) = disp.global_table
table(disp::APIDispatcher, ::Instance) = disp.instance_table
table(disp::APIDispatcher, x) = table(disp, handle(x))

function table(disp::APIDispatcher, device::Device)
    get!(disp.device_tables, device, Dict())
end

table(disp::APIDispatcher, handle::Handle) = table(disp, handle)

function get_fptr(disp::APIDispatcher, handle, key::Symbol)::Ptr{Cvoid}
    val = get(table(disp, handle), key, nothing)
    if isnothing(val) # fetch function pointer
        add_fptr!(disp, handle, key)
    else
        val
    end
end

function add_fptr!(disp::APIDispatcher, handle, key)
    fptr = if isnothing(handle)
        function_pointer(string(key))
    else
        function_pointer(handle, string(key))
    end
    if fptr == C_NULL
        error("Could not retrieve function pointer for '$key'. This can be caused by an extension not being enabled for a function that needs it; see the help with `?` or the documentation for more information.")
    end
    Base.@lock disp.lock table(disp, handle)[key] = fptr
    fptr
end


text:
module Gates

using QuantumCircuits
using QuantumCircuits.QML
# using QGen.Quantum
# using QGen.Quantum.Gates
# using QGen.Quantum.QBase
# using QGen.Quantum.Circuit
# using QGen.Quantum.Registers
# using QGen.Quantum.Math
# using QGen.Genetic.CircOpt
# using QGen.Quantum.Execute
# using QGen.Quantum.Graph
# using QGen.Genetic.CircOpt

export ZZ, XX, YY, trotter, trotter2, findU4paramsZZYYXX, trotterU4,
       findU4paramsZZYYXXx2, trotter2U4, trotter2a

function ZZ(qc, q0, q1, t, usePulse=false)
    if usePulse
        qc.h(q1)
        qc.rzx(q0, q1, t)
        qc.x(q0)
        qc.rzx(q0, q1, -t)
        qc.x(q0)
        qc.h(q1)
    else
        qc.cx(q0, q1)
        qc.rz(q1, 2*t)
        qc.cx(q0, q1)
    end
end

function YY(qc, q0, q1, t, usePulse=false)
    if usePulse
        qc.sdg([q0, q1])
        qc.h(q0)
        qc.rzx(q0, q1, t)
        qc.x(q0)
        qc.rzx(q0, q1, -t)
        qc.x(q0)
        qc.h(q0)
        qc.s([q0, q1])
    else
        qc.rx([q0, q1], π/2)
        qc.cx(q0, q1)
        qc.rz(q1, 2*t)
        qc.cx(q0, q1)
        qc.rx([q0, q1], -π/2)
    end
end

function XX(qc, q0, q1, t, usePulse=false)
    if usePulse
        qc.h(q0)
        qc.rzx(q0, q1, t)
        qc.x(q0)
        qc.rzx(q0, q1, -t)
        qc.x(q0)
        qc.h(q0)
    else
        qc.ry([q0, q1], π/2)
        qc.cx(q0, q1)
        qc.rz(q1, 2*t)
        qc.cx(q0, q1)
        qc.ry([q0, q1], -π/2)
    end
end

function trotter(qc, qubits, t)
    for i in 1:(length(qubits)-1)
        ZZ(qc, qubits[i], qubits[i+1], t)
        YY(qc, qubits[i], qubits[i+1], t)
        XX(qc, qubits[i], qubits[i+1], t)
    end
end

function trotter2(qc, qubits, t)
    for i in 1:(length(qubits)-2)
        ZZ(qc, qubits[i], qubits[i+1], t/2)
        YY(qc, qubits[i], qubits[i+1], t/2)
        XX(qc, qubits[i], qubits[i+1], t/2)
    end

    i = length(qubits) - 1
    ZZ(qc, qubits[i], qubits[i+1], t)
    YY(qc, qubits[i], qubits[i+1], t)
    XX(qc, qubits[i], qubits[i+1], t)

    for i in 1:(length(qubits)-2)
        ZZ(qc, qubits[i], qubits[i+1], t/2)
        YY(qc, qubits[i], qubits[i+1], t/2)
        XX(qc, qubits[i], qubits[i+1], t/2)
    end
end
function trotter2a(qc, qubits, t, isFirst, isLast)
    if isFirst
        for i in 1:(length(qubits)-2)
            ZZ(qc, qubits[i], qubits[i+1], t/2)
            YY(qc, qubits[i], qubits[i+1], t/2)
            XX(qc, qubits[i], qubits[i+1], t/2)
        end
    end

    i = length(qubits) - 1
    ZZ(qc, qubits[i], qubits[i+1], t)
    YY(qc, qubits[i], qubits[i+1], t)
    XX(qc, qubits[i], qubits[i+1], t)

    if isLast
        for i in 1:(length(qubits)-2)
            ZZ(qc, qubits[i], qubits[i+1], t/2)
            YY(qc, qubits[i], qubits[i+1], t/2)
            XX(qc, qubits[i], qubits[i+1], t/2)
        end
    else
        for i in 1:(length(qubits)-2)
            ZZ(qc, qubits[i], qubits[i+1], t)
            YY(qc, qubits[i], qubits[i+1], t)
            XX(qc, qubits[i], qubits[i+1], t)
        end
    end
end
# function trotter2a(qc, qubits, t)
#     for i in 1:(length(qubits)-1)
#         ZZ(qc, qubits[i], qubits[i+1], t/2)
#         YY(qc, qubits[i], qubits[i+1], t/2)
#         XX(qc, qubits[i], qubits[i+1], t/2)
#     end
#     for i in (length(qubits)-1):-1:1
#         ZZ(qc, qubits[i], qubits[i+1], t/2)
#         YY(qc, qubits[i], qubits[i+1], t/2)
#         XX(qc, qubits[i], qubits[i+1], t/2)
#     end
# end

function trotterU4(qc, qubits, params)
    for i in 1:(length(qubits)-1)
        qc.u4(qubits[i], qubits[i+1], params)
    end
end

function trotter2U4(qc, qubits, t, isFirst, isLast, params, params2, params3)
    if isFirst
        for i in 1:(length(qubits)-2)
            qc.u4(qubits[i], qubits[i+1], params2)
        end
    end

    i = length(qubits) - 1
    qc.u4(qubits[i], qubits[i+1], params)

    if isLast
        for i in 1:(length(qubits)-2)
            qc.u4(qubits[i], qubits[i+1], params2)
        end
    else
        for i in 1:(length(qubits)-2)
            qc.u4(qubits[i], qubits[i+1], params3)
        end
    end
end
# function trotter2U4(qc, qubits, t, isFirst, isLast, params, params2, params3)
#     if isFirst
#         for i in 1:(length(qubits)-2)
#             qc.u4(qubits[i], qubits[i+1], params2)
#         end
#     end
#
#     i = length(qubits) - 1
#     qc.u4(qubits[i], qubits[i+1], params2)
#     qc.u4(qubits[i], qubits[i+1], params2)
#
#     if isLast
#         for i in 1:(length(qubits)-2)
#             qc.u4(qubits[i], qubits[i+1], params2)
#         end
#     else
#         for i in 1:(length(qubits)-2)
#             qc.u4(qubits[i], qubits[i+1], params3)
#         end
#     end
# end

function findU4paramsZZYYXX(t; debug=false)
    qc = QCircuit(2)
    ZZ(qc, 0, 1, t)
    YY(qc, 0, 1, t)
    XX(qc, 0, 1, t)
    expmat = tomatrix(qc)

    qr = QuantumRegister(2)
    qc = QCircuit(qr)
    qc.u4(qr[0], qr[1])

    params = getRandParameters(qc)
    setparameters!(qc, params)
    qc = decompose(qc)

    params, _, err, _  = findparam(expmat, qc, debug=debug, trystandard=false)

    @assert err < 1e-5 "The error of U gate should be small but it is $err."
    #println(round.(expmat, digits=2))
    #println(round.(tomatrix(qc), digits=2))
    #println(qc)

    return params
end

function findU4paramsZZYYXXx2(t; debug=false)
    qc = QCircuit(2)
    ZZ(qc, 0, 1, t)
    YY(qc, 0, 1, t)
    XX(qc, 0, 1, t)
    ZZ(qc, 0, 1, t)
    YY(qc, 0, 1, t)
    XX(qc, 0, 1, t)
    expmat = tomatrix(qc)

    qr = QuantumRegister(2)
    qc = QCircuit(qr)
    qc.u4(qr[0], qr[1])

    params = getRandParameters(qc)
    setparameters!(qc, params)
    qc = decompose(qc)

    params, _, err, _  = findparam(expmat, qc, debug=debug, trystandard=false)
    @assert err < 1e-5 "The error of U gate should be small but it is $err."

    #println(round.(expmat, digits=2))
    #println(round.(tomatrix(qc), digits=2))
    #println(qc)

    return params
end

end  # module Gate