text:
using Plots
pyplot()

sizeplot = (800, 800)

rand = randn()

supTitle="\nPRM with connection radius = $connectRadius\n\n"

costTitle= "numsamples vs pathcost\n"
pPRMcost = scatter(nSamples_list, listCosts,
    color=:black,
    title = supTitle * costTitle, ylabel = "euclidean path cost", xlabel = "numsamples",
    xaxis=((0, 300), 0:50:300),
    yaxis=((0,100), 0:20:100))

successTitle = "\nnumsamples vs pSuccess\n"
pPRMsuccess = scatter(nSamples_list, listpSucc,
    #legend=false, yaxis=((0,1.1), 0:0.1:1), xaxis=((0,maxX+20), 0:50:maxX), color = :orange, markersize= 6)
    color = :orange, markersize= 6,
    ylabel = ("P(success) = numSucc/$nTrials trials"), xlabel = "numsamples",
    xaxis=((0, 300), 0:50:300),
    yaxis=((0,1), 0:0.1:1),
    title = successTitle)

plot(pPRMcost, pPRMsuccess, layout=(2,1), window_title="$rand PRM with connection radius = $connectRadius", size = sizeplot,
     #legend=false, xaxis=((0, maxX+50), 0:50:maxX))
    )


text:
import NiLang: Inv
import NiLang.AD: GVar

struct PVar{T,FLT} <: Bundle{T}
    x::T
    logp::FLT
end

PVar(x) = PVar(x, 0.0)
#PVar{T,FLT}(x) where {T,FLT} = PVar(T(x), zero(FLT))
#PVar{T,FLT}(x::Dup) where {T,FLT} = Dup(PVar{T,FLT}(x.x), PVar{T,FLT}(x.twin))
PVar(x::Dup) where T = Dup(PVar(x.x), PVar(x.twin))
(invg::Type{Inv{PVar}})(x::Dup) = Dup((~PVar)(x.x), (~PVar)(x.twin))

PVar{T,FLT}(x::PVar{T,FLT}) where {T,FLT} = x

GVar(x::PVar) = PVar(GVar(x.x), GVar(x.logp))
(invg::Type{Inv{GVar}})(x::PVar) = PVar(invg(x.x), invg(x.logp))

Base.zero(x::PVar) = PVar(zero(x.x))
Base.zero(x::Type{<:PVar{T}}) where T = PVar(zero(T))
Base.copy(x::PVar) = PVar(copy(x.x), copy(x.logp))

# inv kernel and field accesses
invkernel(x::PVar) = x.x
@fieldview NiLang.value(x::PVar) = x.x

# x += field(y) * dt
# TODO: support documentation
@i function update_field(field, x::T, y; dt) where T
    @anc field_out = zero(T)
    get_field(field, field_out, y)
    x += field_out * dt
    (~get_field)(field, field_out, y)
end

# x += field(y) * dt, and update logp.
# TODO: Prove ancilla gradients are discardable! it does not change reversibility.
# i.e. dg_var/dg_ancilla = 0, Beyesian deduction
@i function update_field(field, x::PVar{T}, y::PVar{T}; dt) where T
    @anc field_out = zero(T)
    get_field(field, field_out, y.x)
    # update x
    x.x += field_out * dt

    # update logp
    @routine grad begin
        GVar.((field, field_out, y))
        grad(field_out) ⊕ 1.0
        (~get_field)(field, field_out, y.x)
    end
    x.logp -= grad(y.x) * dt
    ~@routine grad

    (~get_field)(field, field_out, y.x)
end

abstract type Field end
struct LinearField{T} <: RevType
    θ::T
end
Base.zero(lf::LinearField{T}) where T = LinearField(zero(T))
Base.zero(::Type{LinearField{T}}) where T = LinearField(zero(T))
Base.isapprox(lf1::LinearField, lf2::LinearField; kwargs...) = isapprox(lf1.θ, lf2.θ; kwargs...)

NiLang.invkernel(lf::LinearField) = lf.θ

GVar(lf::LinearField) = LinearField(GVar(lf.θ))
(ig::Type{Inv{GVar}})(lf::LinearField) = LinearField((~GVar)(lf.θ))

@i function get_field(lf::LinearField, field_out, x)
    field_out += x * lf.θ
end

ext: jl


text:
# Autogenerated wrapper script for LittleCMS_jll for x86_64-unknown-freebsd
export liblcms2, linkicc, psicc, transicc

using JpegTurbo_jll
using Libtiff_jll
JLLWrappers.@generate_wrapper_header("LittleCMS")
JLLWrappers.@declare_library_product(liblcms2, "liblcms2.so.2")
JLLWrappers.@declare_executable_product(linkicc)
JLLWrappers.@declare_executable_product(psicc)
JLLWrappers.@declare_executable_product(transicc)
function __init__()
    JLLWrappers.@generate_init_header(JpegTurbo_jll, Libtiff_jll)
    JLLWrappers.@init_library_product(
        liblcms2,
        "lib/liblcms2.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        linkicc,
        "bin/linkicc",
    )

    JLLWrappers.@init_executable_product(
        psicc,
        "bin/psicc",
    )

    JLLWrappers.@init_executable_product(
        transicc,
        "bin/transicc",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
    @exception_with_file(
        exception_name::Symbol,
        docstring::Union{Expr, String},
        error_message_bits::Tuple{Vararg{Union{Expr, String}}},
    ) -> Expr

Create an exception with one field: a path to the file.

# Arguments
- `exception_name::Symbol`: name of the exception
- `docstring::Union{Expr, String}`: documentation string
- `error_message_bits::Tuple{Vararg{Union{Expr, String}}}`: strings and expressions which
  will be interpolated in the `showerror` output

# Returns
- `Expr`: an exception definition (struct + `showerror` overload)
"""
@exception exception_with_file file::String context = begin
    @fields(
        "`file::String`: absolute or relative path to a file that is associated with the
        exception",
    )
    @docstring exception_with_file
end


text:
#=
    Support for >= and >
=#
const UnionGT{T} = Union{Strictly{T, MOI.GreaterThan{T}}, MOI.GreaterThan{T}}

function supports_concreteB(concrete_B)
    added_constraints = MOIB.added_constraint_types(concrete_B)
    length(added_constraints) > 1 && return false
    # The inner constraint should not create any variable (might have unexpected consequences)
    return isempty(MOIB.added_constrained_variable_types(concrete_B))
end

function get_num_vars(fct::SAF)
    return length(fct.terms)
end

function get_num_vars(fct::VAF)
    return length(fct.terms)
end

function get_num_vars(fct::MOI.VectorOfVariables)
    return length(fct.variables)
end


text:
using Test; using PyPlot; using Revise; using BenchmarkTools; using AbstractFFTs; using FFTW;

#= test function without discontinuity in derivative at t=0 =#
function funk1(x)
  return (1/2)*x^2 + (1/3)*x^3 - (1/4)*x^4
end

function funk1dot(x)
  return x + x^2 - x^3
end

function funk2(x)
  return x + (1/2)*x^2 - (1/3)*x^3
end

function funk2dot(x)
  return 1 + x - x^2
end

even(x) = x%2==0

function spectralmultiplier(x, L)
  N = length(x)
  out = similar(x, Complex{Float64})
  for k in 0:(N-1)
    if k<(N/2)
      out[k+1] = x[k+1]*2*π*im*k/L

    elseif k>(N/2)
      out[k+1] = x[k+1]*2*π*im*(k-N)/L

    end
  end
  if even(N); out[N÷2]=0.0; end

  return out
end

function specdifftest(f, fdot)
	yff = fft(f)
	yffdot = spectralmultiplier(yff, t[end])
	ydot_approx = real.(ifft(yffdot))
	plot(t, fdot)
	plot(t, ydot_approx, "--", alpha=0.2)
  return
end

dt = 0.01
t = Vector{Float64}(0.0:dt:8*pi)
#y = -cos.(t)
#ydot = sin.(t)

specdifftest(funk1.(t), funk1dot.(t))

max_stars_repo_name: JuliaRheology/RheoBenchAndTest

max_issues_repo_name: JuliaRheology/RheoBenchAndTest

max_forks_repo_name: JuliaRheology/RheoBenchAndTest


text:
using Dates

tomlpath = joinpath(@__DIR__, "sbm_config.toml")
config = Wflow.Config(tomlpath)

model = Wflow.initialize_sbm_model(config)
@unpack network = model

Wflow.load_dynamic_input!(model)
model = Wflow.update(model)

# test if the first timestep was written to the CSV file
flush(model.writer.csv_io)  # ensure the buffer is written fully to disk
@testset "CSV output" begin
    row = csv_first_row(model.writer.csv_path)

    @test row.time == DateTime("2000-01-01T00:00:00")
    @test row.Q ≈ 7.680739298259266f0
    @test row.volume ≈ 2.7783643457819197f7
    @test row.temp_bycoord ≈ 2.3279826641082764f0
    @test row.temp_byindex ≈ 2.3279826641082764f0
    @test row.Q_6336050 ≈ 0.02314152846188506f0
    @test row.Q_6336510 ≈ 0.011829340015725808f0
    @test row.Q_6836100 ≈ 0.004285433507347682f0
    @test row.Q_6336500 ≈ 0.011434655691502343f0
    @test row.Q_6836190 ≈ 0.0004903061800226998f0
    @test row.Q_6336800 ≈ 0.013376156569068251f0
    @test row.Q_6336900 ≈ 0.0033244344482295423f0
    @test row.Q_6336930 ≈ 0.08966505345584352f0
    @test row.Q_6336910 ≈ 0.0020127260170180474f0
    @test row.Q_6336920 ≈ 0.002580440335740291f0
    @test row.Q_6136100 ≈ 0.0008485772628065181f0
    @test row.Q_6136500 ≈ 0.0006749692699499864f0
    @test row.Q_6136520 ≈ 0.0019351916317406491f0
    @test row.Q_6136150 ≈ 0.002126215608742851f0
    @test row.Q_6136151 ≈ 0.0028804236661004205f0
    @test row.Q_6136160 ≈ 3.279985818360295f0
    @test row.Q_6136200 ≈ 1.259358422161174f0
    @test row.Q_6136201 ≈ 5.803861903796768f0
    @test row.Q_6136202 ≈ 1.6624100470042156f0
    @test row.recharge_1 ≈ -0.027398093386017383f0
end

@testset "NetCDF scalar output" begin
    ds = model.writer.dataset_scalar
    @test ds["time"][1] == DateTime("2000-01-01T00:00:00")
    @test ds["Q"][:] ≈ [
        1.66241f0,
        5.803862f0,
        1.2593584f0,
        3.279986f0,
        0.0028804236f0,
        0.0021262157f0,
        0.0019351917f0,
        0.0006749693f0,
        0.00084857724f0,
        0.0025804404f0,
        0.002012726f0,
        0.089665055f0,
        0.0033244344f0,
        0.013376157f0,
        0.0004903062f0,
        0.011434656f0,
        0.0042854333f0,
        0.01182934f0,
        0.023141528f0,
    ]
    @test ds["Q_gauges"].attrib["cf_role"] == "timeseries_id"
    @test ds["temp_index"][:] ≈ [2.3279827f0]
    @test ds["temp_coord"][:] ≈ [2.3279827f0]
    @test keys(ds.dim) == ["time", "Q_gauges", "temp_bycoord", "temp_byindex"]
end

@testset "first timestep" begin
    sbm = model.vertical

    @test sbm.tt[50069] ≈ 1.2999999523162842f0

    @test model.clock.iteration == 2

    @test sbm.θₛ[50069] ≈ 0.48343977332115173f0
    @test sbm.runoff[50069] == 0.0
    @test sbm.soilevap[50069] == 0.0
    @test sbm.snow[50069] ≈ 0.6029989752244306f0
end

# run the second timestep
Wflow.load_dynamic_input!(model)
model = Wflow.update(model)

@testset "second timestep" begin
    sbm = model.vertical
    @test sbm.θₛ[50069] ≈ 0.48343977332115173f0
    @test sbm.runoff[50069] == 0.0
    @test sbm.soilevap[50069] ≈ 0.005865651540305367f0
    @test sbm.snow[50069] ≈ 0.009696763863612956f0
end

@testset "subsurface flow" begin
    ssf = model.lateral.subsurface.ssf
    @test sum(ssf) ≈ 6.368140761295825f7
    @test ssf[network.land.order[1]] ≈ 5.324295565183691f2
    @test ssf[network.land.order[end-100]] ≈ 7.855716879739626f2
    @test ssf[network.land.order[end]] ≈ 2.161246841709492f2
end

@testset "overland flow" begin
    q = model.lateral.land.q_av
    @test sum(q) ≈ 319.69538082148307f0
    @test q[26625] ≈ 0.0
    @test q[39308] ≈ 0.0
    @test q[network.land.order[end]] ≈ 1.0978547520221912f-5
end

@testset "river flow" begin
    q = model.lateral.river.q_av
    @test sum(q) ≈ 2805.266510229795f0
    @test q[1622] ≈ 0.0015498389046670307f0
    @test q[43] ≈ 7.335279907301223f0
    @test q[network.river.order[end]] ≈ 0.006068754125930134f0
end

@testset "reservoir simple" begin
    res = model.lateral.river.reservoir
    @test res.outflow[1] ≈ 0.2174998592483153f0
    @test res.inflow[1] ≈ 49.52026446039592f0
    @test res.volume[1] ≈ 2.776162917050312f7
    @test res.precipitation[1] ≈ 0.1765228509902954f0
    @test res.evaporation[1] ≈ 0.5372688174247742f0
end

# set these variables for comparison in "changed dynamic parameters"
precip = copy(model.vertical.precipitation)
evap = copy(model.vertical.potential_evaporation)
lai = copy(model.vertical.leaf_area_index)
res_evap = copy(model.lateral.river.reservoir.evaporation)

Wflow.close_files(model, delete_output = false)

# test for setting a pit and multithreading multiple basins (by setting 2 extra pits
# resulting in 3 basins)
tomlpath = joinpath(@__DIR__, "sbm_config.toml")
config = Wflow.Config(tomlpath)
config["model"]["pits"] = true
config["input"]["pits"] = "wflow_pits"
config.endtime = DateTime(2000, 1, 9)

model = Wflow.run(config)

@testset "timing" begin
    # clock has been reset
    calendar = get(config, "calendar", "standard")::String
    @test model.clock.time == Wflow.cftime(config.starttime, calendar)
    @test model.clock.iteration == 1
end

@testset "river flow at basin outlets and downstream of one pit" begin
    q = model.lateral.river.q_av
    @test q[3908] ≈ 8.107901948181988f0 # pit/ outlet
    @test q[3920] ≈ 0.008992215117496462f0 # downstream of pit 3908
    @test q[2474] ≈ 156.71758046363857f0 # pit/ outlet
    @test q[5650] ≈ 0.10849551740603759f0 # pit/ outlet
end

# test changing forcing and cyclic LAI parameter
tomlpath = joinpath(@__DIR__, "sbm_config.toml")
config = Wflow.Config(tomlpath)

config.input.vertical.precipitation =
    Dict("scale" => 2.0, "netcdf" => Dict("variable" => Dict("name" => "P")))
config.input.vertical.potential_evaporation = Dict(
    "scale" => 3.0,
    "offset" => 1.50,
    "netcdf" => Dict("variable" => Dict("name" => "PET")),
)
config.input.vertical.leaf_area_index =
    Dict("scale" => 1.6, "netcdf" => Dict("variable" => Dict("name" => "LAI")))

model = Wflow.initialize_sbm_model(config)
Wflow.load_dynamic_input!(model)
model = Wflow.update(model)
Wflow.load_dynamic_input!(model)
model = Wflow.update(model)

@testset "changed dynamic parameters" begin
    res = model.lateral.river.reservoir
    vertical = model.vertical
    @test vertical.precipitation[2] / precip[2] ≈ 2.0f0
    @test (vertical.potential_evaporation[100] - 1.50) / evap[100] ≈ 3.0f0
    @test vertical.leaf_area_index[100] / lai[100] ≈ 1.6f0
    @test (res.evaporation[2] - 1.50) / res_evap[2] ≈ 3.0000012203408635f0
end

# test cyclic river inflow
tomlpath = joinpath(@__DIR__, "sbm_config.toml")
config = Wflow.Config(tomlpath)

config.input.cyclic = ["vertical.leaf_area_index", "lateral.river.inflow"]
config.input.lateral.river.inflow = "inflow"

model = Wflow.initialize_sbm_model(config)
Wflow.load_dynamic_input!(model)
model = Wflow.update(model)
Wflow.load_dynamic_input!(model)
model = Wflow.update(model)

@testset "river inflow (cyclic)" begin
    @test model.lateral.river.inflow[43] ≈ 0.75
    @test model.lateral.river.q_av[43] ≈ 8.08521095942811
end

# test fixed forcing (precipitation = 2.5)
config = Wflow.Config(tomlpath)
config.input.vertical.precipitation = Dict("value" => 2.5)
model = Wflow.initialize_sbm_model(config)
Wflow.load_fixed_forcing(model)

@testset "fixed precipitation forcing (initialize)" begin
    @test maximum(model.vertical.precipitation) ≈ 2.5
    @test minimum(model.vertical.precipitation) ≈ 0.0
    @test all(isapprox.(model.lateral.river.reservoir.precipitation, 2.5))
end

Wflow.load_dynamic_input!(model)
model = Wflow.update(model)

@testset "fixed precipitation forcing (first timestep)" begin
    @test maximum(model.vertical.precipitation) ≈ 2.5
    @test minimum(model.vertical.precipitation) ≈ 0.0
    @test all(isapprox.(model.lateral.river.reservoir.precipitation, 2.5))
end

Wflow.close_files(model, delete_output = false)


text:
@testset "Set: 100 parts" begin
    @info set_100()
end


text:
using Test
using Example2018

@testset "Polynomials form a ring" begin

	z = Polynomial([0])

	for n in 1:200
		deg = rand(0:100)
		p = Polynomial(rand(-100000:100000, deg+1))
		@test p + z == p
		@test z + p == p
		@test p + (-p) == z
	end
end


text:
__precompile__()

module Gnuplot

import Base.writemime

export plot, plot!, xlims!, ylims!, zlims!, write

type GnuplotHandle
    stdin
    stdout
    proc
end

function GnuplotHandle()
    stdout, stdin, proc = readandwrite(`gnuplot`)
    !Base.process_running(proc) && error("could not start gnuplot")
    GnuplotHandle(stdin,stdout,proc)
end

function __init__()
    global gnuplot = GnuplotHandle()
    global figures = Figure[]
    global logfile = tempname()
end

function send_string(s)
    global gnuplot, logfile
    open(logfile,"a") do io; write(io, s * "\n"); end
    w = write(gnuplot.stdin, s * "\n")
    @assert w > 0 "sending \"$s\" to gnuplot failed"
end

set_options(;kwargs...) = send_string(join(["set $(string(k)) \"$v\"" for (k,v) in kwargs], "\n"))
unset_options(args...) = send_string(join(["unset $(string(a))" for a in args], "\n"))

set_term(term) = send_string("set term $term")

function finish()
    finished = tempname()
    send_string("system(\"touch '$finished'\")")
    while !isfile(finished)
        sleep(0.001)
    end
end

abstract Data
abstract Data2D <: Data
abstract Data3D <: Data

colors = (:red, :blue, :green, :orange, :black, :grey, :magenta, :yellow, :purple)

write_data(data::Data) = return

typealias Tuple2{T} Tuple{T,T}
typealias Tuple4{T} Tuple{T,T,T,T}
typealias Tuple6{T} Tuple{T,T,T,T,T,T}
typealias Tuple8{T} Tuple{T,T,T,T,T,T,T,T}
typealias DashType Union{Symbol,Int,AbstractString,Tuple2{Int},Tuple4{Int},Tuple6{Int},Tuple8{Int}}

type DiscreteData2D <: Data2D
    x::AbstractVector
    y::AbstractVector
    label::AbstractString
    style::Symbol  # e.g. lines, dots, points, linespoints/lp, impulses, steps, fsteps, histeps, boxes,
    linecolor::Union{Symbol,Int}  # :auto, 1, 2, :red, :blue
    linewidth::Union{Symbol,Int,Real}  # :auto, 1, 2, 3
    dashtype::DashType  # :auto, 1, :solid, "-._-.", (5,2), (5,3,1,2)
    pointtype::Union{Symbol,AbstractString,Int}  # :auto, 1,
    pointsize::Union{Symbol,Int,Real}  # :auto, 1, 2, 3
    datafile::AbstractString
end

DiscreteData2D(x, y; label="", style=:auto, linecolor=:auto, linewidth=:auto,
    dashtype=:auto, pointtype=:auto, pointsize=:auto) =
        DiscreteData2D(x, y, label, style, linecolor, linewidth, dashtype, pointtype, pointsize, "")

function write_data(data::DiscreteData2D)
    data.datafile = tempname() * ".dat"
    writedlm(data.datafile, [data.x data.y])
end

type DiscreteData3D <: Data3D
    x::AbstractVector
    y::AbstractVector
    z::AbstractVector
    label::AbstractString
    style::Symbol  # e.g. lines, dots, points, linespoints/lp, impulses, steps, fsteps, histeps, boxes,
    linecolor::Union{Symbol,Int}  # :auto, 1, 2, :red, :blue
    linewidth::Union{Symbol,Int,Real}  # :auto, 1, 2, 3
    dashtype::DashType  # :auto, 1, :solid, "-._-.", (5,2), (5,3,1,2)
    pointtype::Union{Symbol,AbstractString,Int}  # :auto, 1,
    pointsize::Union{Symbol,Int,Real}  # :auto, 1, 2, 3
    datafile::AbstractString
end

DiscreteData3D(x, y, z; label="", style=:auto, linecolor=:auto, linewidth=:auto,
    dashtype=:auto, pointtype=:auto, pointsize=:auto) =
        DiscreteData3D(x, y, z, label, style, linecolor, linewidth, dashtype, pointtype, pointsize, "")

function write_data(data::DiscreteData3D)
    data.datafile = tempname() * ".dat"
    writedlm(data.datafile, [data.x data.y data.z])
end

function plot_string(data::Union{DiscreteData2D,DiscreteData3D})
    s = "\"$(data.datafile)\" title \"$(data.label)\""
    if data.style != :auto
        s *= " with $(string(data.style))"
    end
    if data.linecolor != :auto
        if data.linecolor in colors
            s *= " linecolor rgb \"$(string(data.linecolor))\""
        else
            s *= " linecolor $(string(data.linecolor))"
        end
    end
    if data.linewidth != :auto
        s *= " linewidth $(string(data.linewidth))"
    end
    if data.dashtype != :auto
        if isa(data.dashtype, AbstractString)
            s *= " dashtype \"$(data.dashtype)\""
        else
            s *= " dashtype $(string(data.dashtype))"
        end
    end
    if data.pointtype != :auto
        if isa(data.pointtype, AbstractString)
            s *= " pointtype \"$(data.pointtype)\""
        else
            s *= " pointtype $(string(data.pointtype))"
        end
    end
    if data.pointsize != :auto
        s *= " pointsize $(string(data.pointsize))"
    end
    s
end

type FunctionData2D{T} <: Data2D
    f::T
    label::AbstractString
    style::Symbol  # e.g. lines, dots, points, linespoints/lp, impulses, steps, fsteps, histeps, boxes,
    linecolor::Union{Symbol,Int}  # :auto, 1, 2, :red, :blue
    linewidth::Union{Symbol,Int,Real}  # :auto, 1, 2, 3
    dashtype::DashType  # :auto, 1, :solid, "-._-.", (5,2), (5,3,1,2)
    pointtype::Union{Symbol,AbstractString,Int}  # :auto, 1,
    pointsize::Union{Symbol,Int,Real}  # :auto, 1, 2, 3
    datafile::AbstractString
end

FunctionData2D(f; label="", style=:auto, linecolor=:auto, linewidth=:auto,
    dashtype=:auto, pointtype=:auto, pointsize=:auto) =
        FunctionData2D(f, label, style, linecolor, linewidth, dashtype, pointtype, pointsize, "")

function plot_string(data::FunctionData2D)
    s = "\"$(data.datafile)\" title \"$(data.label)\""
    if data.style != :auto
        s *= " with $(string(data.style))"
    end
    if data.linecolor != :auto
        if data.linecolor in colors
            s *= " linecolor rgb \"$(string(data.linecolor))\""
        else
            s *= " linecolor $(string(data.linecolor))"
        end
    end
    if data.linewidth != :auto
        s *= " linewidth $(string(data.linewidth))"
    end
    if data.dashtype != :auto
        if isa(data.dashtype, AbstractString)
            s *= " dashtype \"$(data.dashtype)\""
        else
            s *= " dashtype $(string(data.dashtype))"
        end
    end
    if data.pointtype != :auto
        if isa(data.pointtype, AbstractString)
            s *= " pointtype \"$(data.pointtype)\""
        else
            s *= " pointtype $(string(data.pointtype))"
        end
    end
    if data.pointsize != :auto
        s *= " pointsize $(string(data.pointsize))"
    end
    s
end

abstract Figure

function write_data(fig::Figure)
    for (i,d) in enumerate(fig.data)
        write_data(d)
    end
end

type Figure2D <: Figure
    data::Vector{Data2D}
    title::AbstractString
    xlabel::AbstractString
    ylabel::AbstractString
    xlims::Union{Symbol,Tuple{Real,Real}}
    ylims::Union{Symbol,Tuple{Real,Real}}
    imgfile::AbstractString
end

Figure2D(data; title="", xlabel="", ylabel="", xlims=:auto, ylims=:auto) = Figure2D(data, title, xlabel, ylabel, xlims, ylims, "")
Figure2D() = Figure2D([])

plot_figure(fig::Figure2D) = send_string("plot " * join([plot_string(d) for d in fig.data], ","))

type Figure3D <: Figure
    data::Vector{Data3D}
    title::AbstractString
    xlabel::AbstractString
    ylabel::AbstractString
    xlims::Union{Symbol,Tuple{Real,Real}}
    ylims::Union{Symbol,Tuple{Real,Real}}
    zlims::Union{Symbol,Tuple{Real,Real}}
    imgfile::AbstractString
end

Figure3D(data; title="", xlabel="", ylabel="", xlims=:auto, ylims=:auto, zlims=:auto) =
    Figure3D(data, title, xlabel, ylabel, xlims, ylims, zlims, "")
Figure3D() = Figure3D([])

plot_figure(fig::Figure3D) = send_string("splot " * join([plot_string(d) for d in fig.data], ","))

function discretize(f,xmin,xmax)
	x = collect(linspace(xmin, xmax, 10))

    # add jitter to avoid singularities ?
    for i in 2:length(x)-1
	    x[i] += 1e-6*(rand()-0.5)
	end

	y = [f(x[i]) for i in 1:length(x)]

	δymax = 1e-2*(maximum(y) - minimum(y))
	δxmin = 1e-4*(maximum(x) - minimum(x))
	δxmax = 1e-3*(maximum(x) - minimum(x))

	count = 0
	i = 1
	while i < length(y)
	    δx = x[i+1] - x[i]
	    δy = abs(y[i+1] - y[i])
	    if δx > δxmax || (δx > δxmin && δy > δymax)
	        xmid = (x[i+1] + x[i])/2.
	        insert!(x, i+1, xmid)
	        insert!(y, i+1, f(xmid))
	    else
	        i += 1
	    end
	    count += 1
	    @assert count < 100000
	end

	x,y
end

function plot()
    global figures
    fig = Figure2D()
    push!(figures, fig)
    fig
end

function plot(f,xmin,xmax;
        label="", style=:line, linecolor=:auto, linewidth=:auto,
        dashtype=:auto, pointtype=:auto, pointsize=:auto,
        title="", xlabel="", ylabel="", xlims=:auto, ylims=:auto)
    global figures
    x, y = discretize(f,xmin,xmax)
    data = DiscreteData2D(x, y; label=label, style=style, linecolor=linecolor,
        linewidth=linewidth, dashtype=dashtype, pointtype=pointtype, pointsize=pointsize)
    fig = Figure2D([data], title=title, xlabel=xlabel, ylabel=ylabel,
    	xlims=xlims, ylims=ylims)
    push!(figures, fig)
    fig
end

function plot!(f,xmin,xmax;
        label="", style=:line, linecolor=:auto, linewidth=:auto,
        dashtype=:auto, pointtype=:auto, pointsize=:auto,
        title="", xlabel="", ylabel="", xlims=:auto, ylims=:auto)
    global figures
    fig = figures[end]
    @assert isa(fig, Figure2D)
    x, y = discretize(f,xmin,xmax)
    data = DiscreteData2D(x, y; label=label, style=style, linecolor=linecolor,
        linewidth=linewidth, dashtype=dashtype, pointtype=pointtype, pointsize=pointsize)
    push!(fig.data, data)
    fig
end

plot(x::AbstractVector;
        label="", style=:line, linecolor=:auto, linewidth=:auto,
        dashtype=:auto, pointtype=:auto, pointsize=:auto,
        title="", xlabel="", ylabel="", xlims=:auto, ylims=:auto) =
    plot(1:length(x),x,label=label,style=style,linecolor=linecolor,
        linewidth=linewidth,dashtype=dashtype,pointtype=pointtype,
        pointsize=pointsize,title=title,xlabel=xlabel,ylabel=ylabel,
        xlims=xlims,ylims=ylims)

plot!(x::AbstractVector;
        label="", style=:line, linecolor=:auto, linewidth=:auto,
        dashtype=:auto, pointtype=:auto, pointsize=:auto,
        title="", xlabel="", ylabel="", xlims=:auto, ylims=:auto) =
    plot!(1:length(x),x,label=label,style=style,linecolor=linecolor,
        linewidth=linewidth,dashtype=dashtype,pointtype=pointtype,
        pointsize=pointsize,title=title,xlabel=xlabel,ylabel=ylabel,
        xlims=xlims,ylims=ylims)

function plot(x::AbstractVector, y::AbstractVector;
        label="", style=:line, linecolor=:auto, linewidth=:auto,
        dashtype=:auto, pointtype=:auto, pointsize=:auto,
        title="", xlabel="", ylabel="", xlims=:auto, ylims=:auto)
    global figures
    data = DiscreteData2D(x, y; label=label, style=style, linecolor=linecolor,
        linewidth=linewidth, dashtype=dashtype, pointtype=pointtype, pointsize=pointsize)
    fig = Figure2D([data], title=title, xlabel=xlabel, ylabel=ylabel,
    	xlims=xlims, ylims=ylims)
    push!(figures, fig)
    fig
end

function plot!(x::AbstractVector, y::AbstractVector;
        label="", style=:line, linecolor=:auto, linewidth=:auto,
        dashtype=:auto, pointtype=:auto, pointsize=:auto)
    global figures
    fig = figures[end]
    @assert isa(fig, Figure2D)
    data = DiscreteData2D(x, y; label=label, style=style, linecolor=linecolor,
        linewidth=linewidth, dashtype=dashtype, pointtype=pointtype, pointsize=pointsize)
    push!(fig.data, data)
    fig
end

function plot(x::AbstractVector, y::AbstractVector, z::AbstractVector;
        label="", style=:line, linecolor=:auto, linewidth=:auto,
        dashtype=:auto, pointtype=:auto, pointsize=:auto,
        title="", xlabel="", ylabel="", xlims=:auto, ylims=:auto, zlims=:auto)
    global figures
    data = DiscreteData3D(x, y, z; label=label, style=style, linecolor=linecolor,
        linewidth=linewidth, dashtype=dashtype, pointtype=pointtype, pointsize=pointsize)
    fig = Figure3D([data], title=title, xlabel=xlabel, ylabel=ylabel, xlims=xlims,
        ylims=ylims, zlims=zlims)
    push!(figures, fig)
    fig
end

function plot!(x::AbstractVector, y::AbstractVector, z::AbstractVector;
        label="", style=:line, linecolor=:auto, linewidth=:auto,
        dashtype=:auto, pointtype=:auto, pointsize=:auto)
    global figures
    fig = figures[end]
    @assert isa(fig, Figure3D)
    data = DiscreteData3D(x, y, z; label=label, style=style, linecolor=linecolor,
        linewidth=linewidth, dashtype=dashtype, pointtype=pointtype, pointsize=pointsize)
    push!(fig.data, data)
    fig
end

function xlims!(xmin,xmax)
	global figures
    fig = figures[end]
	fig.xlims = (xmin,xmax)
	fig
end

function ylims!(ymin,ymax)
	global figures
    fig = figures[end]
	fig.ylims = (ymin,ymax)
	fig
end

function zlims!(zmin,zmax)
    global figures
    fig = figures[end]
    @assert typeof(fig) == Figure3D
    fig.zlims = (zmin,zmax)
    fig
end

function set_range(fig::Figure2D)
    if fig.xlims != :auto
        send_string("set xrange [$(fig.xlims[1]):$(fig.xlims[2])]")
    end
    if fig.ylims != :auto
        send_string("set yrange [$(fig.ylims[1]):$(fig.ylims[2])]")
    end
end

function set_range(fig::Figure3D)
    if fig.xlims != :auto
        send_string("set xrange [$(fig.xlims[1]):$(fig.xlims[2])]")
    end
    if fig.ylims != :auto
        send_string("set yrange [$(fig.ylims[1]):$(fig.ylims[2])]")
    end
    if fig.zlims != :auto
        send_string("set zrange [$(fig.zlims[1]):$(fig.zlims[2])]")
    end
end

function writemime(io::IO, ::MIME"image/svg+xml", fig::Figure)
    fig.imgfile = tempname()
    set_term("svg dashed")
    set_options(
        output=fig.imgfile,
        title=fig.title,
        xlabel=fig.xlabel,
        ylabel=fig.ylabel,
    )
    set_range(fig)
    write_data(fig)
    plot_figure(fig)
    send_string("unset output")
    finish()
    img = open(readbytes, fig.imgfile, "r")
    write(io, img)
    send_string("reset")
end

ext_to_term = Dict(
    "pdf"  => "pdfcairo",
    "png"  => "pngcairo",
    "svg"  => "svg",
    "eps"  => "epscairo",
    "html" => "canvas",
)

function Base.write(filename::AbstractString, fig::Figure)
    fig.imgfile = filename
    ext = splitext(filename)[2][2:end]     # file extension (without the .)
    term = ext_to_term[ext]                # choose terminal automatically based on file type
    set_term("$term")
    set_options(
        output=fig.imgfile,
        title=fig.title,
        xlabel=fig.xlabel,
        ylabel=fig.ylabel,
    )
    set_range(fig)
    write_data(fig)
    plot_figure(fig)
    send_string("reset")                    # reset all options
    send_string("unset output")             # flushes output
end

# function writemime(io::IO, ::MIME"image/png", fig::Figure)
#     fig.imgfile = tempname()
#     set_term("pngcairo dashed")
#     set_options(
#         output=fig.imgfile,
#         title=fig.title,
#         xlabel=fig.xlabel,
#         ylabel=fig.ylabel
#     )
#     write_data(fig)
#     plot_figure(fig)
#     finish()
#     img = open(readbytes, fig.imgfile, "r")
#     write(io, img)
# end

end  # module Gnuplot

max_stars_repo_name: johncsnyder/Gnuplot


text:
"""
GenX: An Configurable Capacity Expansion Model
Copyright (C) 2021,  Massachusetts Institute of Technology
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
A complete copy of the GNU General Public License v2 (GPLv2) is available
in LICENSE.txt.  Users uncompressing this from an archive may not have
received this license file.  If not, see <http://www.gnu.org/licenses/>.
"""

function write_time_weights(path::AbstractString, inputs::Dict)
	T = inputs["T"]     # Number of time steps (hours)
	# Save array of weights for each time period (when using time sampling)
	dfTimeWeights = DataFrame(Time=1:T, Weight=inputs["omega"])
	CSV.write(joinpath(path, "time_weights.csv"), dfTimeWeights)
end

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
using Logging
using Random
using ProgressMeter
using PyMNE # not in Unfold.jl currently
using MixedModelsPermutations

function cluster_permutation_test(mres, dat, times, coeffOfInterest::Array; kwargs...)
    df_p = DataFrame()
    for coef in coeffOfInterest
        df_p_single = cluster_permutation_test(mres, dat, times, coef; kwargs...)
        append!(df_p, df_p_single)
    end
    return df_p
end

function cluster_permutation_test(
    mres::UnfoldLinearMixedModel,
    dat::Array,
    times::StepRangeLen,
    coeffOfInterest::Int;
    clusterFormingThreshold = 2,
    adjacency = nothing,
    nPerm = 500,
    test_times = (minimum(times), maximum(times)),
    stat = :β,
)

    # from where to where to actually test (can be subrange)
    @assert(test_times[1] < test_times[2])
    tRange = findfirst(times .>= test_times[1]):findfirst(times .>= test_times[2])

    # extract observed z-values
    stat_obs = [
        getproperty(m, stat) for m in coefpvalues(mres.modelfit) if
        String(m.coefname) == coefnames(formula(designmatrix(mres)))[2][coeffOfInterest]
    ]

    # cluster observed
    obs_cluster = pymne_cluster(
        stat_obs,
        clusterFormingThreshold;
        tRange = tRange,
        adjacency = adjacency,
    )

    # if there is no cluster in the real data, return empty p_df
    if isempty(obs_cluster[2])
        p_df = DataFrame(:from => [], :to => [], :pval => [], :coefname => [])
        return p_df
    end

    # this is the meat of the function.
    permDat = cluster_permutation(mres, dat, tRange, coeffOfInterest, nPerm; stat = stat)

    # cluster permuted
    perm_cluster = [
        pymne_cluster(permDat[:, p], clusterFormingThreshold; adjacency = adjacency) for
        p = 1:nPerm
    ]

    # get maximal cluster value
    perm_H0 = maximum.([abs.(p[2]) for p in perm_cluster if ~isempty(p[2])])
    perm_H0 = vcat(perm_H0, zeros(nPerm - length(perm_H0)))

    # get pvalues
    p_vals = PyMNE.stats.cluster_level._pval_from_histogram(obs_cluster[2], perm_H0, 0)

    fs = 1 ./ Float64(times.step)

    # get p-values, we have to shift by the test_times[1] though
    if clusterFormingThreshold == "tfce"
        # not sure how to handle this nicer...
        fromList = [o.start ./ fs + times[1] for o in obs_cluster[1]]
        toList = [o.stop ./ fs + times[1] for o in obs_cluster[1]]
    else
        fromList = [o[1].start ./ fs + times[1] for o in obs_cluster[1]]
        toList = [o[1].stop ./ fs + times[1] for o in obs_cluster[1]]
    end
    p_df = DataFrame(:from => fromList, :to => toList, :pval => p_vals)

    p_df[!, :coefname] .= coefnames(formula(mres))[2][coeffOfInterest] # yeah this might break :S
    return p_df
    #println(p_vals)
    #h = hist(perm_H0,bins=100)
    #vlines!(h.axis,abs.(obs_cluster[2]),color="red")

end

cluster_permutation(args...; kwargs...) =
    cluster_permutation(MersenneTwister(1), args...; kwargs...)

function cluster_permutation(
    rng::AbstractRNG,
    mres,
    dat,
    tRange,
    coeffOfInterest,
    nPerm;
    stat = :β,
)
    permDat = Matrix{Float64}(undef, length(tRange), nPerm)
    mm_outer = Unfold.LinearMixedModel_wrapper(formula(mres), dat[1, 1, :], modelmatrix(mres))

    chIx = 1 # for now we only support 1 channel anyway
    #
    #p = Progress(length(tRange))
    #Threads.@threads for tIx =1:length(tRange)
    @showprogress "Processing Timepoints" for tIx = 1:length(tRange)
        #println(string(tIx))
        # splice in the correct data for residual calculation
        mm = deepcopy(mm_outer)
        mm.y .= dat[chIx, tRange[tIx], :]

        # set the previous calculated model-fit
        updateL!(setθ!(mm, Vector(mres.modelfit.θ[tRange[tIx]])))

        # get the coefficient
        H0 = coef(mm)
        # set the one of interest to 0
        H0[coeffOfInterest] = 0
        # run the permutation

        perm = undef
        Logging.with_logger(NullLogger()) do   # remove NLopt warnings
         perm = permutation(
            deepcopy(rng), # important here is to set the same seed to keep flip all time-points the same
            nPerm,
            mm;
            β = H0,
            blup_method = olsranef,
            use_threads = false,
            hide_progress = true,
            inflation_method= (x,y,z) ->I(length(H0)),#MixedModelsPermutations.inflation_factor(mm, olsranef(mm), residuals(mm))

        ) # constant rng to keep autocorr & olsranef for singular models
        end

        # extract the test-statistic

        #perm_z = [m.z for m in perm.coefpvalues if String(m.coefname)==coefnames(mm)[coeffOfInterest]]
        permDat[tIx, :] = [
            getproperty(m, stat) for m in perm.coefpvalues if
            String(m.coefname) == coefnames(mm)[coeffOfInterest]
        ]

        #next!(p)
    end
    return permDat
end

# function to call pymne -> _find_cluster function
function pymne_cluster(
    data,
    clusterFormingThreshold;
    tRange = 1:size(data, 1),
    adjacency = nothing,
)

    #return PyMNE.stats.cluster_level._find_clusters(data,clusterFormingThreshold,adjacency=adjacency,include=tRange)
    # this raised an error due to the named arguments, but was working before. not sure whats going on

    # need to convert to boolean
    tRange_bool = zeros(size(data, 1))
    tRange_bool[tRange] .= 1
    return PyMNE.stats.cluster_level._find_clusters(
        data,
        clusterFormingThreshold,
        0,
        adjacency,
        1,
        tRange_bool,
    )
end

pymne_cluster(data, clusterFormingThreshold::String; kwargs...) =
    pymne_cluster(data, Dict(:start => 0, :step => 0.2); kwargs...)

max_stars_repo_name: unfoldtoolbox/jlunfold

max_issues_repo_name: unfoldtoolbox/jlunfold

max_forks_repo_name: unfoldtoolbox/jlunfold


text:
module PSCOPFFixtures

using PSCOPF

using Dates
using DataStructures
using Printf

"""
creates a network with 2 buses
"""
function network_2buses(;limit::Float64=35., ptdf1=0.5, ptdf2=-0.5)
    network = PSCOPF.Networks.Network()
    # Buses
    PSCOPF.Networks.add_new_bus!(network, "bus_1")
    PSCOPF.Networks.add_new_bus!(network, "bus_2")
    # Branches
    PSCOPF.Networks.add_new_branch!(network, "branch_1_2", limit);
    # PTDF
    PSCOPF.Networks.add_ptdf_elt!(network, "branch_1_2", "bus_1", ptdf1)
    PSCOPF.Networks.add_ptdf_elt!(network, "branch_1_2", "bus_2", ptdf2)

    return network
end

#=
    ECH = 10h
    TS: [11h, 11h15]
    S: [S1,S2]
                        bus 1                   bus 2
                        |                      |
    (limitable) wind_1_1|       "1_2"          |
    Pmin=0, Pmax=100    |                      |
    Csta=0, Cprop=1     |                      |
      S1: 20    S1: 15  |----------------------|
      S2: 30    S2: 30  |         35           |
                        |                      |
                        |                      |
    (pilotable) prod_1_1|                      |(pilotable) prod_2_1
    Pmin=10, Pmax=100   |                      | Pmin=10, Pmax=100
    Csta=45k, Cprop=10  |                      | Csta=80k, Cprop=15
INIT: ON                |                      |INIT: ON
                        |                      | DP=DMO=2h
                        |                      |
           load(bus_1)  |                      |load(bus_2)
      S1: 10     S1: 17 |                      | S1: 40  S1: 48
      S2: 10     S2: 13 |                      | S2: 45  S2: 52
=#
"""
Creates a context with preset uncertainties and for a preset 2 buses network :
    - on bus_1 : a limitable and an pilotable
    - on bus_2 : an pilotable
    pilotable unit of bus_2 has a DP=DMO=2h
Arguments :
    TS  : a 2 elements DateTime vector for the target timesteps
    ech : DateTime giving the current horizon timepoint
"""
function context_2buses_2TS_2S(TS, ech;
                                logs=nothing)
    @assert(length(TS) == 2)
    network = PSCOPFFixtures.network_2buses()

    # Limitables
    PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "wind_1_1", PSCOPF.Networks.LIMITABLE,
                                            0., 100.,
                                            0., 1.,
                                            Dates.Second(0), Dates.Second(0))
    # Pilotables
    PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "prod_1_1", PSCOPF.Networks.PILOTABLE,
                                            10., 100.,
                                            45000., 10.,
                                            Dates.Second(0), Dates.Second(0))
    PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_2", "prod_2_1", PSCOPF.Networks.PILOTABLE,
                                            10., 100.,
                                            80000., 15.,
                                            Dates.Second(2*60*60), Dates.Second(2*60*60))
    # initial generators state
    generators_init_state = SortedDict(
                    "prod_1_1" => PSCOPF.ON,
                    "prod_2_1" => PSCOPF.ON
                )
    # Uncertainties
    uncertainties = PSCOPF.Uncertainties()
    PSCOPF.add_uncertainty!(uncertainties, ech, "wind_1_1", TS[1], "S1", 20.)
    PSCOPF.add_uncertainty!(uncertainties, ech, "wind_1_1", TS[1], "S2", 30.)
    PSCOPF.add_uncertainty!(uncertainties, ech, "wind_1_1", TS[2], "S1", 15.)
    PSCOPF.add_uncertainty!(uncertainties, ech, "wind_1_1", TS[2], "S2", 30.)
    PSCOPF.add_uncertainty!(uncertainties, ech, "bus_1", TS[1], "S1", 10.)
    PSCOPF.add_uncertainty!(uncertainties, ech, "bus_1", TS[1], "S2", 10.)
    PSCOPF.add_uncertainty!(uncertainties, ech, "bus_1", TS[2], "S1", 17.)
    PSCOPF.add_uncertainty!(uncertainties, ech, "bus_1", TS[2], "S2", 13.)
    PSCOPF.add_uncertainty!(uncertainties, ech, "bus_2", TS[1], "S1", 40.)
    PSCOPF.add_uncertainty!(uncertainties, ech, "bus_2", TS[1], "S2", 45.)
    PSCOPF.add_uncertainty!(uncertainties, ech, "bus_2", TS[2], "S1", 48.)
    PSCOPF.add_uncertainty!(uncertainties, ech, "bus_2", TS[2], "S2", 52.)

    context = PSCOPF.PSCOPFContext(network, TS, PSCOPF.PSCOPF_MODE_1,
                                generators_init_state,
                                uncertainties, nothing,
                                logs)

    return context
end

end #module PSCOPFFixtures

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
# Unified API Functions

export getActiveBackend, getBackendCount, AF_BACKEND_DEFAULT,
        AF_BACKEND_CPU, AF_BACKEND_CUDA, AF_BACKEND_OPENCL,
        getAvailableBackends, setBackend, getBackendId, sync,
        getActiveBackendId, getDeviceId, getVersionNumber

AF_BACKEND_DEFAULT = UInt32(0)
AF_BACKEND_CPU = UInt32(1)
AF_BACKEND_CUDA = UInt32(2)
AF_BACKEND_OPENCL = UInt32(4)

# Available only v3.3.0 onwards
function getActiveBackend()
    backend = getActiveBackendId()
    if backend == 0
        println("Default Backend")
    elseif backend == 1
        println("CPU Backend")
    elseif backend == 2
        println("CUDA Backend")
    elseif backend == 4
        println("OpenCL Backend")
    end
end

# Available only v3.3.0 onwards
function getActiveBackendId()
    maj, min, pat = getVersionNumber()
    if maj < 3 || (maj == 3 && min < 3)
        throw(ErrorException("This function is supported only on arrayfire v3.3.0 and onward, current version is v$(maj).$(min).$(pat). Please install a newer version of arrayfire."))
    end
    backend = Base.Ref{Cuint}(0)
    af_get_active_backend(backend)
    backend = Int(backend[])
end

function getAvailableBackends()
    backends = Base.Ref{Cint}(0)
    af_get_available_backends(backends)
    b = Int(backends[])
    if b == 0
        println("None")
    elseif b == 1
        println("CPU")
    elseif b == 2
        println("CUDA")
    elseif b == 3
        println("CPU and CUDA")
    elseif b == 4
        println("OpenCL")
    elseif b == 5
        println("CPU and OpenCL")
    elseif b == 6
        println("CUDA and OpenCL")
    elseif b == 7
        println("CPU, CUDA and OpenCL")
    end
end

function getBackendCount()
    b = Base.Ref{Cuint}(0)
    af_get_backend_count(b)
    Int(b[])
end

function setBackend(back::Cuint)
    af_set_backend(back)
end

function getBackendId(a::AFArray)
    backend = Base.Ref{Cuint}(0)
    af_get_backend_id(backend, a)
    backend = backend[]
    if backend == 0
        println("Default Backend")
    elseif backend == 1
        println("CPU Backend")
    elseif backend == 2
        println("CUDA Backend")
    elseif backend == 4
        println("OpenCL Backend")
    end
    backend
end

function sync()
    b = getActiveBackendId()
    af_sync(b)
end

function sync(b::Int)
    af_sync(b)
end

# There is a af_get_device_id only in ArrayFire  > v3.3.x, possibly only in v3.4.0 onwards

function getDeviceId(a::AFArray)
    maj, min, pat = getVersionNumber()
    if maj < 3 || (maj == 3 && min < 4)
        throw("ArrayFire Version is v$maj.$min.$pat, which does not support this function. Install the latest version of ArrayFire.")
    end
    device = Base.Ref{Cint}(0)
    af_get_device_id(device, a)
    Int(device[])
end

function getVersionNumber()
    maj = Base.Ref{Cint}(0)
    min = Base.Ref{Cint}(0)
    pat = Base.Ref{Cint}(0)
    af_get_version(maj, min, pat)
    maj[], min[], pat[]
end


text:
const stdcomment = "// Autogenerated by BitSAD"

function gethandler end

is_hardware_primitive(sig...) = is_trace_primitive(sig...)

function generatehw(io::IO, f, args...;
                    top = _nameof(f),
                    submodules = [],
                    transforms = [constantreduction!])
    # get tape and transform
    # if f itself is a primitive, do a manual tape
    if is_hardware_primitive(Ghost.get_type_parameters(Ghost.call_signature(f, args...))...)
        tape = Ghost.Tape()
        inputs = Ghost.inputs!(tape, f, args...)
        if _isstruct(f)
            tape.result = push!(tape, Ghost.mkcall(inputs...))
        else
            tape.result = push!(tape, Ghost.mkcall(f, inputs[2:end]...))
        end
    else
        tape = trace(f, args...;
                    isprimitive = is_hardware_primitive,
                    submodules = submodules)
    end
    transform!(_unbroadcast, tape)
    transform!(_squash_binary_vararg, tape)
    tape = Ghost.Tape(tape.ops, tape.result, tape.parent, tape.meta, TupleCtx())
    transform!(_record_tuples_and_splats, tape)
    transform!(_reroute_tuple_index, tape)
    transform!(_desplat, tape)
    # transform!(_squash_tuple_index, tape)

    # extract tape into module
    m = Module(fn = f, name = top)
    extracttrace!(m, tape)
    tape = nothing # don't hold onto tape

    # replace constant SBitstreams with calls to SBitstream constructor on floats
    insertrng!(m)

    # apply transformations
    foreach(t! -> t!(m), transforms)

    # replace constants with Verilog strings
    constantreplacement!(m)

    # generate verilog string
    generateverilog(io, m)

    return io, m
end
function generatehw(f, args...; kwargs...)
    buffer = IOBuffer()
    buffer, m = generatehw(buffer, f, args...; kwargs...)

    return String(take!(buffer)), m
end

_getstructname(::T) where T = lowercase(string(nameof(T)))

_handle_parameter(parameter::Number, submodules) = parameter
_handle_parameter(parameter::AbstractArray{<:Number}, submodules) = parameter
_handle_parameter(parameter::Tuple, submodules) = map(parameter) do p
    _handle_parameter(p, submodules)
end
function _handle_parameter(::T, submodules) where T
    (T ∈ submodules) || @warn "Parameter of type $T will be ignored (cannot encode in Verilog)"

    return nothing
end

function _handle_getproperty!(m::Module, call, param_map, const_map)
    if m.fn == _gettapeval(call.args[1])
        val = _handle_parameter(_gettapeval(call), m.submodules)
        prop = string(_gettapeval(call.args[2]))
        encodable = (val isa Union{AbstractArray, Tuple}) ? !all(isnothing, val) : !isnothing(val)
        if encodable
            m.parameters[prop] = val
            param_map[_getid(call)] = prop
        end
    else
        const_map[_getid(call)] = string(_gettapeval(call))
    end

    return m
end

function _get_materialize_origin(x)
    origin = _gettapeop(x)

    return (origin.fn == Base.materialize) ? _get_materialize_origin(origin.args[1]) : x
end

function extracttrace!(m::Module, tape::Ghost.Tape)
    # we use ids instead of Ghost.Variables as keys b/c
    # there are hashing issues with Ghost.Variable
    # see https://github.com/dfdx/Ghost.jl/issues/20
    param_map = Dict{Int, String}()
    const_map = Dict{Int, String}()
    materialize_map = Dict{Int, Tuple{Ghost.Variable, Any}}()
    # skip first call which is the function being compiled
    for call in tape
        if call isa Ghost.Call
            # remap constants on tape
            fn = (_gettapeval(call.fn) isa Ghost.Constant) ? _gettapeval(call.fn).val : call.fn

            # ignore materialize calls
            if fn == Base.materialize
                origin = _get_materialize_origin(call.args[1])
                materialize_map[_getid(call)] = (origin, _gettapeval(call))
                continue
            end

            # handle calls to getproperty
            if fn == Base.getproperty
                _handle_getproperty!(m, call, param_map, const_map)
                continue
            end

            is_hardware_primitive(Ghost.get_type_parameters(Ghost.call_signature(tape, call))...) ||
                continue

            # create Operator for Ghost.Call (handling broadcast)
            # structs are renamed as Foo -> foo_$id
            # plain functions are name ""
            name = _isstruct(fn) ? _getstructname(fn) * "_$(_getid(fn))" : ""
            isbroadcast = _isbcast(fn)
            fn = isbroadcast ? _gettapeval(call.args[1]) : fn
            op = (name = Symbol(name), type = typeof(fn), broadcasted = isbroadcast)

            # map inputs and outputs of Ghost.Call to Nets
            # set args that are Ghost.Input to :input class
            inputs = map(call.args[(1 + isbroadcast):end]) do arg
                arg, val = get(materialize_map, _getid(arg), (arg, _gettapeval(arg)))
                name = haskey(param_map, _getid(arg)) ? param_map[_getid(arg)] :
                       haskey(const_map, _getid(arg)) ? const_map[_getid(arg)] :
                       _isvariable(arg) ? "net_$(_getid(arg))" : string(val)
                net = Net(val; name = name)

                if _isvariable(arg)
                    if _isinput(arg)
                        net = setclass(net, :input)
                    elseif haskey(param_map, _getid(arg))
                        net = setclass(net, :parameter)
                    elseif haskey(const_map, _getid(arg))
                        net = setclass(net, :constant)
                    end
                else # treat all non-variables as constants
                    net = setclass(net, :constant)
                end

                return net
            end
            outval = isbroadcast ? Base.materialize(_gettapeval(call)) : _gettapeval(call)
            output = Net(outval; name = "net_$(_getid(call))")

            if tape.result == Ghost.Variable(call)
                output = setclass(output, :output)
            end

            addnode!(m, inputs, [output], op)
        end
    end

    return m
end


text:
include("../utils/testUtils.jl")
using qccdSimulator.QCCDevControl_Types
using qccdSimulator.QCCDDevControl
using qccdSimulator.QCCDev_Feasible

# ========= JSON tests =========
function readJSONOK(path::String)::Bool
    qccd1 = readJSON(path)
    qccd2 = giveQccDes()
    return checkEqualQCCD(qccd1,qccd2)
end

function readTimeJSONOK(path ::String)
    readTimeJSON(path)
    @assert OperationTimes[:load] == 5
    @assert OperationTimes[:linear_transport] == 78
    @assert OperationTimes[:loadingHole_transport] == 35
    @assert OperationTimes[:swap] == 2
    @assert OperationTimes[:split] == 55
    return true
end

function readTimeJSONfail(paths ::Array{String})
   errormsg1 = "Time values can't be negative"
   errormsg2 = "invalid JSON"
   errorcount = 0

   for path ∈ paths
        try
            readTimeJSON(path)
        catch e
            @assert startswith(e.msg, errormsg1) || startswith(e.msg, errormsg2)
            errorcount += 1
        end
   end
   @assert errorcount == length(paths)
   return true
end

function readTimeJSONnoFile()
    try
        readTimeJSON("foo")
    catch e
        @assert startswith(e.msg, "Input is not a file")
    end

    return true
end
# ========= END JSON tests =========

# ========= Device comparison tests =========
function QCCDevCtrlOKTest()::Bool
    qccd1 = giveQccCtrl()
    qccd2 = QCCDevCtrl(giveQccDes())
    return checkEqualQCCDevCtrl(qccd1,qccd2)
end

function checkEqualQCCD(qccd1::QCCDevDescription, qccd2::QCCDevDescription):: Bool
    # Compare gate zones
    gateZones1 = qccd1.gateZone.gateZones
    gateZones2 = qccd2.gateZone.gateZones
    @assert size(gateZones1) == size(gateZones2)
    for (gateZone1,gateZone2) in zip(gateZones1,gateZones2)
        @assert gateZone1.id == gateZone2.id
        @assert gateZone1.end0 == gateZone2.end0
        @assert gateZone1.end1 == gateZone2.end1
        @assert gateZone1.capacity == gateZone2.capacity
    end
    # Compare auxiliary zones
    auxZones1 = qccd1.auxZone.auxZones
    auxZones2 = qccd2.auxZone.auxZones
    @assert size(auxZones1) == size(auxZones2)
    for (auxZone1,auxZone2) in zip(auxZones1,auxZones2)
        @assert auxZone1.id == auxZone2.id
        @assert auxZone1.end0 == auxZone2.end0
        @assert auxZone1.end1 == auxZone2.end1
        @assert auxZone1.capacity == auxZone2.capacity
    end
    # Compare junctions
    juns1 = qccd1.junction.junctions
    juns2 = qccd2.junction.junctions
    @assert size(juns1) == size(juns2)
    for (jun1,jun2) in zip(juns1,juns2)
        @assert jun1.id == jun2.id
        @assert jun1.type == jun2.type
    end
    # Compare loading zones
    loadZones1 = qccd1.loadZone.loadZones
    loadZones2 = qccd2.loadZone.loadZones
    @assert size(loadZones1) == size(loadZones2)
    for (loadZone1,loadZone2) in zip(loadZones1,loadZones2)
        @assert loadZone1.id == loadZone2.id
        @assert loadZone1.end0 == loadZone2.end0
        @assert loadZone1.end1 == loadZone2.end1
    end
    return true
end

function checkEqualQCCDevCtrl(qccdc1::QCCDevControl,qccdc2::QCCDevControl):: Bool
    @assert qccdc1.t_now == qccdc2.t_now
    @assert qccdc1.simulate == qccdc2.simulate
    @assert qccdc1.qnoise_esimate == qccdc2.qnoise_esimate
    @assert checkEqualQCCD(qccdc1.dev, qccdc2.dev)
    @assert length(qccdc1.gateZones) == length(qccdc2.gateZones)
    @assert length(qccdc1.junctions) == length(qccdc2.junctions)
    @assert length(qccdc1.auxZones) == length(qccdc2.auxZones)
    @assert length(qccdc1.loadingZones) == length(qccdc2.loadingZones)
    for (key,value) in qccdc1.gateZones
        @assert haskey(qccdc2.gateZones, key)
        @assert qccdc2.gateZones[key].id == value.id
        @assert qccdc2.gateZones[key].capacity == value.capacity
        @assert qccdc2.gateZones[key].chain == value.chain
        @assert qccdc2.gateZones[key].end0 == value.end0
        @assert qccdc2.gateZones[key].end1 == value.end1
    end
    for (key,value) in qccdc1.auxZones
        @assert haskey(qccdc2.auxZones, key)
        @assert qccdc2.auxZones[key].id == value.id
        @assert qccdc2.auxZones[key].capacity == value.capacity
        @assert qccdc2.auxZones[key].chain == value.chain
        @assert qccdc2.auxZones[key].end0 == value.end0
        @assert qccdc2.auxZones[key].end1 == value.end1
    end
    for (key,value) in qccdc1.loadingZones
        @assert haskey(qccdc2.loadingZones, key)
        @assert qccdc2.loadingZones[key].id == value.id
        @assert qccdc2.loadingZones[key].end0 == value.end0
        @assert qccdc2.loadingZones[key].end1 == value.end1
        @assert qccdc2.loadingZones[key].hole == value.hole == nothing
    end
    for (key,value) in qccdc1.junctions
        @assert haskey(qccdc2.junctions, key)
        @assert qccdc2.junctions[key].id == value.id
        @assert qccdc2.junctions[key].ends == value.ends
        @assert qccdc2.junctions[key].type == value.type
    end
    return true
end

function QCCDevCtrlTest()
    qdd::QCCDevDescription = giveQccDes()
    return QCCDevCtrl(qdd)
end
# ========= END Device comparison tests =========

# ========= Junction tests =========
function initJunctionsTest()
    nJunctions = rand(5:30)
    juncTypes = [string(type) for type ∈ JunctionType]
    zones, _junctions = giveZonesJunctions(nJunctions, rand(juncTypes, nJunctions))
    #Convert the generic zones to Gate, auxiliary, and loading zones.
    gateZones = ZoneInfoDesc[]
    auxZones = ZoneInfoDesc[]
    loadZones = LoadZoneInfoDesc[]
    for zone ∈ zones
        r = rand()
        if 0 ≤ r ≤ 0.33
            push!(gateZones, zone)
        elseif 0.33 < r ≤ 0.66
            push!(auxZones, zone)
        else
            load = LoadZoneInfoDesc(zone.id,zone.end0,zone.end1)
            push!(loadZones, load)
        end
    end

    junctions = qccdSimulator.QCCDDevControl._initJunctions(gateZones,
                                            auxZones, loadZones, _junctions)
    for (k,junction) in junctions
        @assert k == junction.id
        juncType = filter(x-> Symbol(x.id)==k,_junctions)[1].type #[1] : get first and only element
        juncType = Symbol(juncType)
        @assert junction.type == juncType
        @assert junction.zoneType == :junction
        @assert length(junction.ends) == typesSizes[juncType]
        for zoneId in junction.ends
            zone = filter(x -> Symbol(x.id) == zoneId, zones)[1] #[1] : get first and only element
            @assert Symbol(zone.end0) == k || Symbol(zone.end1) == k
        end
    end
    return true
end

function initJunctionsTestRepId()
    zones, _junctions = giveZonesJunctions(2, ["T","T"];repJunc = true)
    try
        qccdSimulator.QCCDDevControl._initJunctions(zones, nothing, nothing, _junctions)
    catch e
        @assert startswith(e.msg, "Repeated junction ID")
        return true
    end
    return false
end

function initJunctionsTestIsolated()
    zones, _junctions = giveZonesJunctions(2, ["T","T"];isolatedJunc = true)
    try
        qccdSimulator.QCCDDevControl._initJunctions(nothing, zones, nothing, _junctions)
    catch e
        @assert endswith(e.msg, "is not connected to anything.")
        return true
    end
    return false
end

function initJunctionsTestWrongType()
    zones, _junctions = giveZonesJunctions(2, ["T","T"];wrongJuncType = true)
    loadZones = LoadZoneInfoDesc[]
    for zone ∈ zones
        load = LoadZoneInfoDesc(zone.id,zone.end0,zone.end1)
        push!(loadZones, load)
    end
    try
        qccdSimulator.QCCDDevControl._initJunctions(nothing,nothing,loadZones, _junctions)
    catch e
        @assert(e.msg == "Junction with ID 1 of type T has 2 ends. It should have 3 ends.")
        return true
    end
    return false
end
# ========= END Junction tests =========

# ========= Adjacency tests  =========
function initAdjacencyTest()
    device ::QCCDevDescription = giveQccDes()
    adjacency = qccdSimulator.QCCDDevControl._initAdjacency(device)

    @show adjacency

    for (key, value) in adjacency
        auxGate = filter(x-> Symbol(x.id)==key, device.gateZone.gateZones)
        auxLoad = filter(x-> Symbol(x.id)==key, device.loadZone.loadZones)
        auxAux = filter(x-> Symbol(x.id)==key, device.auxZone.auxZones)
        if !isempty(auxGate)
            @assert sort(value) == sort([auxGate[0].end0, auxGate[0].end1])
        elseif !isempty(auxLoad)
            @assert sort(value) == sort([auxLoad[0].end0, auxLoad[0].end1])
        elseif !isempty(auxAux)
            @assert sort(value) == sort([auxAux[0].end0, auxAux[0].end1])
        end
    end
end

# ========= END Adjacency tests  =========

# ========= Loading zones tests =========
function initLoadingZoneTest()
    loadZoneDesc::LoadZoneDesc = giveQccDes().loadZone
    loadZones = qccdSimulator.QCCDDevControl._initLoadingZones(loadZoneDesc)
    for (key, value) in loadZones
        @assert key == value.id
        aux = filter(x-> Symbol(x.id)==key, loadZoneDesc.loadZones)
        @assert length(aux) == 1
        aux = aux[1]
        tmp = aux.end0 == "" ? nothing : Symbol(aux.end0)
        @assert tmp == value.end0
        tmp = aux.end1 == "" ? nothing : Symbol(aux.end1)
        @assert tmp == value.end1
        @assert value.zoneType == :loadingZone
    end
    return true
end

function initLoadingZoneRepeatedIdTest()
    loadZoneDesc::LoadZoneDesc = giveLoadZoneDescRepeatedId()
    return qccdSimulator.QCCDDevControl._initLoadingZones(loadZoneDesc)
end
# ========= END Loading zones tests =========

# ========= Auxiliary and Gate zones tests =========
function initGateZoneTest()
    gateZoneDesc::GateZoneDesc = giveQccDes().gateZone
    gateZones = qccdSimulator.QCCDDevControl._initGateZone(gateZoneDesc)
    for (key, value) in gateZones
        @assert key == value.id
        aux = filter(x-> Symbol(x.id)==key, gateZoneDesc.gateZones)
        @assert length(aux) == 1
        aux = aux[1]
        @assert aux.capacity == value.capacity
        @assert isempty(value.chain)
        tmp = aux.end0 == "" ? nothing : Symbol(aux.end0)
        @assert tmp == value.end0
        tmp = aux.end1 == "" ? nothing : Symbol(aux.end1)
        @assert tmp == value.end1
        @assert value.zoneType == :gateZone
    end
    return true
end

function initGateZoneRepeatedIdTest()
    gateZoneDesc::GateZoneDesc = giveGateZoneDescRepeatedId()
    return qccdSimulator.QCCDDevControl._initGateZone(gateZoneDesc)
end

function checkInitErrorsTest()
    qdd::QCCDevControl = giveQccCtrl()
    qccdSimulator.QCCDDevControl._checkInitErrors(qdd.junctions,qdd.auxZones,
                                                 qdd.gateZones, qdd.loadingZones)
    return true
end

function initAuxGateZonesTestRepId()
    zones, _ = giveZonesJunctions(2, ["T","T"];repZone = true)
    auxDesc = AuxZoneDesc(zones)
    try
        qccdSimulator.QCCDDevControl._initAuxZones(auxDesc)
    catch e
        @assert startswith(e.msg, "Repeated")
    end
    gateDesc = GateZoneDesc(zones)
    try
        qccdSimulator.QCCDDevControl._initGateZone(gateDesc)
    catch e
        @assert startswith(e.msg, "Repeated")
    end
    return true
end

function initAuxGateZonesTestInvZone()
    try
        qccdSimulator.QCCDDevControl._initAuxZones(AuxZoneDesc(giveZoneInfo(rand(5:10);invZone=true)))
    catch e
        @assert endswith(e.msg, "\"end0\" and \"end1\" must be different")
    end
    try
        qccdSimulator.QCCDDevControl._initGateZone(GateZoneDesc(giveZoneInfo(rand(5:10);invZone=true)))
    catch e
        @assert endswith(e.msg, "\"end0\" and \"end1\" must be different")
    end
    return true
end

function initAuxGateZonesTestWithNothing()
    zones = giveZoneInfo(rand(5:10);giveNothing=true)
    qccdSimulator.QCCDDevControl._initAuxZones(AuxZoneDesc(zones))
    qccdSimulator.QCCDDevControl._initGateZone(GateZoneDesc(zones))
    return true
end

function initAuxZonesTest()
    _auxZones = AuxZoneDesc(giveZoneInfo(rand(10:15)))
    auxZones = qccdSimulator.QCCDDevControl._initAuxZones(_auxZones)
    @assert length(_auxZones.auxZones) == length(auxZones)
    for _auxZone in _auxZones.auxZones
        auxZone = auxZones[Symbol(_auxZone.id)]
        @assert _auxZone.end0 == string(auxZone.end0)
        @assert _auxZone.end1 == string(auxZone.end1)
        @assert _auxZone.capacity == auxZone.capacity
        @assert auxZone.zoneType == :auxZone
    end
    return true
end

# ========= END Auxiliary and Gate zones tests =========

# ========= Init functions check tests =========
function checkInitErrorsTest()
    qdd::QCCDevControl = giveQccCtrl()
    qccdSimulator.QCCDDevControl._checkInitErrors(qdd.junctions,qdd.auxZones,
                                                 qdd.gateZones, qdd.loadingZones)
    return true
end

"""
This function checks the edge cases of __auxCheck.
It repeats the same errors for aux, gate, and loading zones.
"""
function checkInitErrorsTestEdgeCases()
    qdd::QCCDevControl = giveQccCtrl()

    #Aux zone with 'nothing' in both ends should throw error:
    _qdd = deepcopy(qdd)
    k = rand(keys(_qdd.auxZones))
    _qdd.auxZones[k] = AuxZone(k,_qdd.auxZones[k].capacity,nothing,nothing)
    try
        qccdSimulator.QCCDDevControl._checkInitErrors(_qdd.junctions,_qdd.auxZones,
                                                 _qdd.gateZones, _qdd.loadingZones)
    catch e
        @assert endswith(e.msg, "Element cannot be isolated")
    end
    _qdd = deepcopy(qdd)

    # Randomly pick aux zones (by keys) and change their end0 or end1 to be nonsense.
    # This should throw error too.
    for k ∈ rand(keys(_qdd.auxZones), min(length(keys(_qdd.auxZones)), rand(1:10)))
        _qdd = deepcopy(qdd)

        if !isnothing(_qdd.auxZones[k].end0)
            _qdd.auxZones[k] = AuxZone(k,_qdd.auxZones[k].capacity,:nonsense,_qdd.auxZones[k].end1)
        else
            _qdd.auxZones[k] = AuxZone(k,_qdd.auxZones[k].capacity,_qdd.auxZones[k].end0,:nonsense)
        end
        try
            qccdSimulator.QCCDDevControl._checkInitErrors(_qdd.junctions,_qdd.auxZones,
                                                 _qdd.gateZones, _qdd.loadingZones)
        catch e
            @assert endswith(e.msg, "is wrong connected.")
        end
    end

    #Gate zone with 'nothing' in both ends should throw error:
    _qdd = deepcopy(qdd)
    k = rand(keys(_qdd.gateZones))
    _qdd.gateZones[k] = GateZone(k,_qdd.gateZones[k].capacity,nothing,nothing)
    try
        qccdSimulator.QCCDDevControl._checkInitErrors(_qdd.junctions,_qdd.auxZones,
                                                 _qdd.gateZones, _qdd.loadingZones)
    catch e
        @assert endswith(e.msg, "Element cannot be isolated")
    end
    _qdd = deepcopy(qdd)

    # Randomly pick gate zones (by keys) and change their end0 or end1 to be nonsense.
    # This should throw error too.
    for k ∈ rand(keys(_qdd.gateZones), min(length(keys(_qdd.gateZones)), rand(1:10)))
        _qdd = deepcopy(qdd)

        if !isnothing(_qdd.gateZones[k].end0)
            _qdd.gateZones[k] = GateZone(k,_qdd.gateZones[k].capacity,:nonsense,_qdd.gateZones[k].end1)
        else
            _qdd.gateZones[k] = GateZone(k,_qdd.gateZones[k].capacity,_qdd.gateZones[k].end0,:nonsense)
        end
        try
            qccdSimulator.QCCDDevControl._checkInitErrors(_qdd.junctions,_qdd.auxZones,
                                                 _qdd.gateZones, _qdd.loadingZones)
        catch e
            @assert endswith(e.msg, "is wrong connected.")
        end
    end

    #Load zone with 'nothing' in both ends should throw error:
    _qdd = deepcopy(qdd)
    k = rand(keys(_qdd.loadingZones))
    _qdd.loadingZones[k] = LoadingZone(k,nothing,nothing)
    try
        qccdSimulator.QCCDDevControl._checkInitErrors(_qdd.junctions,_qdd.auxZones,
                                                 _qdd.gateZones, _qdd.loadingZones)
    catch e
        @assert endswith(e.msg, "Element cannot be isolated")
    end
    _qdd = deepcopy(qdd)

    # Randomly pick load zones (by keys) and change their end0 or end1 to be nonsense.
    # This should throw error too.
    for k ∈ rand(keys(_qdd.loadingZones), min(length(keys(_qdd.loadingZones)), rand(1:10)))
        _qdd = deepcopy(qdd)

        if !isnothing(_qdd.loadingZones[k].end0)
            _qdd.loadingZones[k] = LoadingZone(k,:nonsense,_qdd.loadingZones[k].end1)
        else
            _qdd.loadingZones[k] = LoadingZone(k,_qdd.gateZones[k].end0,:nonsense)
        end
        try
            qccdSimulator.QCCDDevControl._checkInitErrors(_qdd.junctions,_qdd.auxZones,
                                                 _qdd.gateZones, _qdd.loadingZones)
        catch e
            @assert endswith(e.msg, "is wrong connected.")
        end
    end
    return true
end
# ========= END Init functions check tests =========

# ========= linear_transport tests =========
function linearTransportTestOK()
    qdd::QCCDevControl = giveQccCtrl()
    qdd.qubits[1] = Qubit(1,Symbol(8))
    qdd.loadingZones[Symbol(8)].hole = 1

    t = qccdSimulator.QCCDDevControl.linear_transport(qdd, 10, 1, Symbol(3))
    @assert OperationTimes[:linear_transport] == 78
    @assert t == 88
    @assert qdd.t_now == 88
    @assert qdd.qubits[1].position === Symbol(3)
    @assert isnothing(qdd.loadingZones[Symbol(8)].hole)
    @assert qdd.gateZones[Symbol(3)].chain == [[1]]
    return true
end

function linearTransportTest1()
    qdd::QCCDevControl = giveQccCtrl()
    qdd.qubits[1] = Qubit(1,Symbol(3))
    qdd.qubits[1].status = :inGateZone
    qdd.qubits[1].destination = Symbol(7)

    append!(qdd.gateZones[Symbol(3)].chain, [[1],[2]])
    append!(qdd.auxZones[Symbol(7)].chain, [[3]])
    qccdSimulator.QCCDDevControl.linear_transport(qdd, 10, 1, Symbol(7))
    @assert qdd.qubits[1].position === Symbol(7)
    @assert qdd.gateZones[Symbol(3)].chain == [[2]]
    @assert qdd.auxZones[Symbol(7)].chain == [[3],[1]]
    @assert isnothing(qdd.qubits[1].destination)

    return true
end

function linearTransportTest2()
    qdd::QCCDevControl = giveQccCtrl()
    qdd.qubits[1] = Qubit(1,Symbol(3))
    qdd.qubits[1].status = :inGateZone
    qdd.qubits[1].destination = Symbol(7)

    append!(qdd.gateZones[Symbol(3)].chain, [[2],[1]])
    qccdSimulator.QCCDDevControl.linear_transport(qdd, 10, 1, Symbol(8))
    @assert qdd.qubits[1].position === Symbol(8)
    @assert qdd.gateZones[Symbol(3)].chain == [[2]]
    @assert qdd.loadingZones[Symbol(8)].hole == 1
    @assert qdd.qubits[1].destination == Symbol(7)

    return true
end

function linearTransportTest3()
    qdd::QCCDevControl = giveQccCtrl()
    qdd.qubits[1] = Qubit(1,Symbol(7))
    qdd.qubits[1].status = :inGateZone

    append!(qdd.gateZones[Symbol(3)].chain, [[2]])
    append!(qdd.auxZones[Symbol(7)].chain, [[3],[1]])
    qccdSimulator.QCCDDevControl.linear_transport(qdd, 10, 1, Symbol(3))
    @assert qdd.qubits[1].position === Symbol(3)
    @assert qdd.gateZones[Symbol(3)].chain == [[1],[2]]
    @assert qdd.auxZones[Symbol(7)].chain == [[3]]
    @assert isnothing(qdd.qubits[1].destination)

    return true
end

function isallowedLinearTransportTestTime()
    qdd::QCCDevControl = giveQccCtrl()
    try
        isallowed_linear_transport(qdd, qdd.t_now - 5, 1, :a)
    catch e
        @assert startswith(e.msg,"Time must be higher than")
        return true
    end
    return false
end

function isallowedLinearTransportTestNoIon()
    qdd::QCCDevControl = giveQccCtrl()
    try
        isallowed_linear_transport(qdd, qdd.t_now+1, 1, :a)
    catch e
        @assert startswith(e.msg,"Ion with ID 1 is not in device")
        return true
    end
    return false
end

function isallowedLinearTransportTestNoZone()
    qdd::QCCDevControl = giveQccCtrl()
    qdd.qubits[1] = Qubit(1,Symbol(8))
    qdd.loadingZones[Symbol(8)].hole = 1
    try
        isallowed_linear_transport(qdd, qdd.t_now+1, 1, :nonsense)
    catch e
        @assert startswith(e.msg,"Zone with ID nonsense is not in device")
        return true
    end
    return false
end

function isallowedLinearTransportTestNonAdjacent()
    qdd::QCCDevControl = giveQccCtrl()
    qdd.qubits[1] = Qubit(1,Symbol(8))
    qdd.loadingZones[Symbol(8)].hole = 1
    try
        isallowed_linear_transport(qdd, qdd.t_now+1, 1, Symbol(2))
    catch e
        @assert startswith(e.msg,"Can't do linear transport to a non-adjacent zone.")
        return true
    end
    return false
end

function isallowedLinearTransportTestAllGood()
    qdd::QCCDevControl = giveQccCtrl()
    qdd.qubits[1] = Qubit(1,Symbol(8))
    qdd.loadingZones[Symbol(8)].hole = 1
    isallowed_linear_transport(qdd, qdd.t_now+1, 1, Symbol(3))
    return true
end

function isallowedLinearTransportTestFull()
    dest = Symbol(3)

    # "Load" ion
    qdd::QCCDevControl = giveQccCtrl(;alternateDesc=true)
    qdd.qubits[1] = Qubit(1,Symbol(8))
    qdd.loadingZones[Symbol(8)].hole = 1

    # "Fill up" destination
    halvedCapacity = floor(Int64, qdd.gateZones[dest].capacity / 2)
        # Multiple chains to check if reduce works
    push!(qdd.gateZones[dest].chain,
        rand(Int64, halvedCapacity), rand(Int64, halvedCapacity))

    try
        isallowed_linear_transport(qdd, qdd.t_now+1, 1, dest)
    catch e
        @assert endswith(e.msg,"$dest cannot hold more ions.")
        return true
    end
    return false
end

function isallowedLinearTransportTestBlockedEnd0()
    origin = Symbol(3)

    # "Load" ion
    qdd::QCCDevControl = giveQccCtrl(;alternateDesc=true)
    qdd.qubits[1] = Qubit(1,origin)
    halvedCapacity = floor(Int64, qdd.gateZones[origin].capacity / 2)
    push!(qdd.gateZones[origin].chain, rand(2:100, halvedCapacity), [1])
    halvedCapacity -= 1
    push!(qdd.gateZones[origin].chain, rand(2:100, halvedCapacity))

    try
        isallowed_linear_transport(qdd, qdd.t_now+1, 1, Symbol(7))
    catch e
        @assert endswith(e.msg,"not in the correct end position.")
        return true
    end
    return false
end

function isallowedLinearTransportTestBlockedEnd1()
    origin = Symbol(3)

    # "Load" ion
    qdd::QCCDevControl = giveQccCtrl(;alternateDesc=true)
    qdd.qubits[1] = Qubit(1,origin)
    halvedCapacity = floor(Int64, qdd.gateZones[origin].capacity / 2)
    push!(qdd.gateZones[origin].chain, rand(2:100, halvedCapacity), [1])
    halvedCapacity -= 1
    push!(qdd.gateZones[origin].chain, rand(2:100, halvedCapacity))

    try
        isallowed_linear_transport(qdd, qdd.t_now+1, 1, Symbol(8))
    catch e
        @assert endswith(e.msg,"not in the correct end position.")
        return true
    end
    return false
end

function isallowedLinearTransportTestNotBlockedEnd0()
    origin = Symbol(3)

    # "Load" ion
    qdd::QCCDevControl = giveQccCtrl(;alternateDesc=true)
    qdd.qubits[1] = Qubit(1,origin)
    halvedCapacity = floor(Int64, qdd.gateZones[origin].capacity / 2)
    push!(qdd.gateZones[origin].chain,[1], rand(2:100, halvedCapacity))
    halvedCapacity -= 1
    push!(qdd.gateZones[origin].chain, rand(2:100, halvedCapacity))
    isallowed_linear_transport(qdd, qdd.t_now+1, 1, Symbol(7))
    return true
end

function isallowedLinearTransportTestNotBlockedEnd1()
    origin = Symbol(3)

    # "Load" ion
    qdd::QCCDevControl = giveQccCtrl(;alternateDesc=true)
    qdd.qubits[1] = Qubit(1,origin)
    halvedCapacity = floor(Int64, qdd.gateZones[origin].capacity / 2)
    push!(qdd.gateZones[origin].chain, rand(2:100, halvedCapacity))
    halvedCapacity -= 1
    push!(qdd.gateZones[origin].chain, rand(2:100, halvedCapacity), [1])
    isallowed_linear_transport(qdd, qdd.t_now+1, 1, Symbol(8))
    return true
end
# ========= END linear_transport tests =========

# ========= utils functions tests =========
function giveZoneTest()
    qdd::QCCDevControl = giveQccCtrl()

    zone = giveZone(qdd, :nonsense)
    @assert isnothing(zone)

    k = rand(keys(qdd.gateZones))
    zone = giveZone(qdd, k)
    @assert zone === qdd.gateZones[k]

    k = rand(keys(qdd.junctions))
    zone = giveZone(qdd, k)
    @assert zone === qdd.junctions[k]

    k = rand(keys(qdd.auxZones))
    zone = giveZone(qdd, k)
    @assert zone === qdd.auxZones[k]

    k = rand(keys(qdd.loadingZones))
    zone = giveZone(qdd, k)
    @assert zone === qdd.loadingZones[k]

    return true
end

# ========= END utils functions tests =========

# ========= _time_check functiong test =========
time_check_timeFailsTest() = qccdSimulator.QCCDev_Feasible._time_check( 10, 9, :load)
time_check_modelFailsTest() = qccdSimulator.QCCDev_Feasible._time_check( 10, 12, :test)
time_checkOKTest() = qccdSimulator.QCCDev_Feasible._time_check( 10, 12, :load)
# ========= _time_check functiong test =========

# ========= load functiong test =========
function initQubitTest()
    for i in 1:25
        qubit1 = Qubit(i,:test)
        qubit2 = qccdSimulator.QCCDDevControl.initQubit(:test)
        @assert qubit1.id == qubit2.id == i
        @assert qubit1.status == qubit2.status == :inLoadingZone
        @assert qubit1.position == qubit2.position == :test
        @assert qubit1.destination == qubit2.destination == nothing
    end
    return true
end

function isallowedLoad_zoneNotExistTest()
    qccd = giveQccCtrl()
    isallowed_load(qccd,:test,4)
    return true
end

function isallowedLoad_loadingHoleBusyTest()
    qccd = giveQccCtrl()
    qccd.loadingZones[Symbol(8)].hole = 3
    isallowed_load(qccd,Symbol(8),4)
    return true
end

function isallowedLoad_OK()
    qccd = giveQccCtrl()
    isallowed_load(qccd,Symbol(8),4)
    return true
end

function loadOKTest()
    qccd = giveQccCtrl()
    tmp = qccdSimulator.QCCDDevControl.load(qccd,2,Symbol(8))
    @assert tmp.new_ion_idx == qccdSimulator.QCCDDevControl.NIONS
    @assert tmp.t₀ == 7
    @assert haskey(qccd.qubits, tmp.new_ion_idx)
    @assert qccd.qubits[tmp.new_ion_idx].position == Symbol(8)
    @assert qccd.qubits[tmp.new_ion_idx].destination == nothing
    return true
end
# ========= END load functiong test =========

# ========= swap function test =========
function isallowedSwap_qubitNotExist()
    qccd = giveQccCtrl()
    isallowed_swap(qccd, 99, 1, 10)
end

function isallowedSwap_qubitNotSameZone()
    qccd = giveQccCtrl()
    ion1 = giveQubit(Symbol(1),1)
    ion2 = giveQubit(Symbol(2),2)
    qccd.qubits[ion1.id] = deepcopy(ion1)
    qccd.qubits[ion2.id] = deepcopy(ion2)
    push!(qccd.gateZones[Symbol(1)].chain, [ion1.id])
    push!(qccd.gateZones[Symbol(2)].chain, [ion2.id])
    isallowed_swap(qccd, ion1.id, ion2.id, 10)
end

function isallowedSwap_qubitNotGateZone()
    qccd = giveQccCtrl()
    ion1 = giveQubit(Symbol(4),1)
    ion2 = giveQubit(Symbol(4),2)
    qccd.qubits[ion1.id] = deepcopy(ion1)
    qccd.qubits[ion2.id] = deepcopy(ion2)
    push!(qccd.auxZones[Symbol(4)].chain, [ion1.id])
    push!(qccd.auxZones[Symbol(4)].chain[1], ion2.id)
    isallowed_swap(qccd, ion1.id, ion2.id, 10)
end

function isallowedSwap_qubitsNotAdjacents1()
    qccd = giveQccCtrl()
    ion1 = giveQubit(Symbol(1),1)
    ion2 = giveQubit(Symbol(1),2)
    ion3 = giveQubit(Symbol(1),3)
    qccd.qubits[ion1.id] = deepcopy(ion1)
    qccd.qubits[ion2.id] = deepcopy(ion2)
    qccd.qubits[ion3.id] = deepcopy(ion3)
    push!(qccd.gateZones[Symbol(1)].chain, [ion1.id])
    push!(qccd.gateZones[Symbol(1)].chain[1], ion2.id)
    push!(qccd.gateZones[Symbol(1)].chain[1], ion3.id)
    isallowed_swap(qccd, ion1.id, ion3.id, 10)
end

function isallowedSwap_qubitsNotAdjacents2()
    qccd = giveQccCtrl()
    ion1 = giveQubit(Symbol(1),1)
    ion2 = giveQubit(Symbol(1),2)
    ion3 = giveQubit(Symbol(1),3)
    ion4 = giveQubit(Symbol(1),4)
    ion5 = giveQubit(Symbol(1),5)
    ion6 = giveQubit(Symbol(1),6)
    qccd.qubits[ion1.id] = deepcopy(ion1)
    qccd.qubits[ion2.id] = deepcopy(ion2)
    qccd.qubits[ion3.id] = deepcopy(ion3)
    qccd.qubits[ion4.id] = deepcopy(ion4)
    qccd.qubits[ion5.id] = deepcopy(ion5)
    qccd.qubits[ion6.id] = deepcopy(ion6)
    push!(qccd.gateZones[Symbol(1)].chain, [ion1.id])
    push!(qccd.gateZones[Symbol(1)].chain[1], ion2.id)
    push!(qccd.gateZones[Symbol(1)].chain[1], ion3.id)
    push!(qccd.gateZones[Symbol(1)].chain, [ion4.id])
    push!(qccd.gateZones[Symbol(1)].chain[2], ion5.id)
    push!(qccd.gateZones[Symbol(1)].chain[2], ion6.id)
    isallowed_swap(qccd, ion4.id, ion6.id, 10)
end

function isallowedSwap_OK()
    qccd = giveQccCtrl()
    ion1 = giveQubit(Symbol(1),1)
    ion2 = giveQubit(Symbol(1),2)
    ion3 = giveQubit(Symbol(1),3)
    qccd.qubits[ion1.id] = deepcopy(ion1)
    qccd.qubits[ion2.id] = deepcopy(ion2)
    qccd.qubits[ion3.id] = deepcopy(ion3)
    push!(qccd.gateZones[Symbol(1)].chain, [ion1.id])
    push!(qccd.gateZones[Symbol(1)].chain[1], ion2.id)
    push!(qccd.gateZones[Symbol(1)].chain[1], ion3.id)
    isallowed_swap(qccd, ion2.id, ion3.id, 10)
end

function swap_OK1()
    qccd = giveQccCtrl()
    ion1 = giveQubit(Symbol(1),1)
    ion2 = giveQubit(Symbol(1),2)
    ion3 = giveQubit(Symbol(1),3)
    qccd.qubits[ion1.id] = deepcopy(ion1)
    qccd.qubits[ion2.id] = deepcopy(ion2)
    qccd.qubits[ion3.id] = deepcopy(ion3)
    push!(qccd.gateZones[Symbol(1)].chain, [ion1.id])
    push!(qccd.gateZones[Symbol(1)].chain[1], ion2.id)
    push!(qccd.gateZones[Symbol(1)].chain[1], ion3.id)
    @assert qccd.gateZones[Symbol(1)].chain == [[1,2,3]]
    qccdSimulator.QCCDDevControl.swap(qccd, 2, ion2.id,ion3.id)
    @assert qccd.gateZones[Symbol(1)].chain == [[1,3,2]]
    @assert qccd.t_now == 4
    return true
end

function isallowedSwap_qubitsNotSameChain1()
    qccd = giveQccCtrl()
    ion1 = giveQubit(Symbol(1),1)
    ion2 = giveQubit(Symbol(1),2)
    ion3 = giveQubit(Symbol(1),3)
    ion4 = giveQubit(Symbol(1),4)
    ion5 = giveQubit(Symbol(1),5)
    ion6 = giveQubit(Symbol(1),6)
    qccd.qubits[ion1.id] = deepcopy(ion1)
    qccd.qubits[ion2.id] = deepcopy(ion2)
    qccd.qubits[ion3.id] = deepcopy(ion3)
    qccd.qubits[ion4.id] = deepcopy(ion4)
    qccd.qubits[ion5.id] = deepcopy(ion5)
    qccd.qubits[ion6.id] = deepcopy(ion6)
    push!(qccd.gateZones[Symbol(1)].chain, [ion1.id])
    push!(qccd.gateZones[Symbol(1)].chain[1], ion2.id)
    push!(qccd.gateZones[Symbol(1)].chain[1], ion3.id)
    push!(qccd.gateZones[Symbol(1)].chain, [ion4.id])
    push!(qccd.gateZones[Symbol(1)].chain[2], ion5.id)
    push!(qccd.gateZones[Symbol(1)].chain[2], ion6.id)
    isallowed_swap(qccd, ion3.id, ion5.id, 10)
end

function isallowedSwap_qubitsNotSameChain2()
    qccd = giveQccCtrl()
    ion1 = giveQubit(Symbol(1),1)
    ion2 = giveQubit(Symbol(1),2)
    ion3 = giveQubit(Symbol(1),3)
    ion4 = giveQubit(Symbol(1),4)
    ion5 = giveQubit(Symbol(1),5)
    ion6 = giveQubit(Symbol(1),6)
    qccd.qubits[ion1.id] = deepcopy(ion1)
    qccd.qubits[ion2.id] = deepcopy(ion2)
    qccd.qubits[ion3.id] = deepcopy(ion3)
    qccd.qubits[ion4.id] = deepcopy(ion4)
    qccd.qubits[ion5.id] = deepcopy(ion5)
    qccd.qubits[ion6.id] = deepcopy(ion6)
    push!(qccd.gateZones[Symbol(1)].chain, [ion1.id])
    push!(qccd.gateZones[Symbol(1)].chain[1], ion2.id)
    push!(qccd.gateZones[Symbol(1)].chain[1], ion3.id)
    push!(qccd.gateZones[Symbol(1)].chain, [ion4.id])
    push!(qccd.gateZones[Symbol(1)].chain[2], ion5.id)
    push!(qccd.gateZones[Symbol(1)].chain[2], ion6.id)
    isallowed_swap(qccd, ion3.id, ion4.id, 10)
end

function swap_OK2()
    qccd = giveQccCtrl()
    ion1 = giveQubit(Symbol(1),1)
    ion2 = giveQubit(Symbol(1),2)
    ion3 = giveQubit(Symbol(1),3)
    ion4 = giveQubit(Symbol(1),4)
    ion5 = giveQubit(Symbol(1),5)
    ion6 = giveQubit(Symbol(1),6)
    qccd.qubits[ion1.id] = deepcopy(ion1)
    qccd.qubits[ion2.id] = deepcopy(ion2)
    qccd.qubits[ion3.id] = deepcopy(ion3)
    qccd.qubits[ion4.id] = deepcopy(ion4)
    qccd.qubits[ion5.id] = deepcopy(ion5)
    qccd.qubits[ion6.id] = deepcopy(ion6)
    push!(qccd.gateZones[Symbol(1)].chain, [ion1.id])
    push!(qccd.gateZones[Symbol(1)].chain[1], ion2.id)
    push!(qccd.gateZones[Symbol(1)].chain[1], ion3.id)
    push!(qccd.gateZones[Symbol(1)].chain, [ion4.id])
    push!(qccd.gateZones[Symbol(1)].chain[2], ion5.id)
    push!(qccd.gateZones[Symbol(1)].chain[2], ion6.id)
    @assert qccd.gateZones[Symbol(1)].chain == [[1,2,3],[4,5,6]]
    qccdSimulator.QCCDDevControl.swap(qccd, 2, ion5.id,ion6.id)
    @assert qccd.gateZones[Symbol(1)].chain == [[1,2,3],[4,6,5]]
    @assert qccd.t_now == 4
    return true
end
# ========= END swap function test =========

max_stars_repo_name: TartuQC/trapped-ion-qccd-simulator

max_issues_repo_name: TartuQC/trapped-ion-qccd-simulator

max_forks_repo_name: TartuQC/trapped-ion-qccd-simulator