text:
using Base.Test
using ExecutableSpecifications.Gherkin
using ExecutableSpecifications.Gherkin: ScenarioStep
using ExecutableSpecifications
using ExecutableSpecifications: StepDefinitionContext, StepDefinition, StepDefinitionLocation
using ExecutableSpecifications: Executor, StepExecutionResult, QuietRealTimePresenter, executefeature
import ExecutableSpecifications: present

successful_step_definition(::StepDefinitionContext) = ExecutableSpecifications.SuccessfulStepExecution()
failed_step_definition(::StepDefinitionContext) = ExecutableSpecifications.StepFailed("")
error_step_definition(::StepDefinitionContext) = error("Some error")

struct FakeStepDefinitionMatcher <: ExecutableSpecifications.StepDefinitionMatcher
    steps::Dict{ExecutableSpecifications.Gherkin.ScenarioStep, Function}
end

ExecutableSpecifications.findstepdefinition(s::FakeStepDefinitionMatcher, step::ExecutableSpecifications.Gherkin.ScenarioStep) = StepDefinition("some text", s.steps[step], StepDefinitionLocation("", 0))

struct ThrowingStepDefinitionMatcher <: ExecutableSpecifications.StepDefinitionMatcher
    ex::Exception
end

ExecutableSpecifications.findstepdefinition(matcher::ThrowingStepDefinitionMatcher, ::ExecutableSpecifications.Gherkin.ScenarioStep) = throw(matcher.ex)

@testset "Executor        " begin
    @testset "Execute a one-step scenario; No matching step found; Result is NoStepDefinitionFound" begin
        stepdefmatcher = ThrowingStepDefinitionMatcher(ExecutableSpecifications.NoMatchingStepDefinition())
        executor = ExecutableSpecifications.Executor(stepdefmatcher)
        scenario = Scenario("Description", [], [Given("some precondition")])

        scenarioresult = ExecutableSpecifications.executescenario(executor, scenario)

        @test isa(scenarioresult.steps[1], ExecutableSpecifications.NoStepDefinitionFound)
    end

    @testset "Execute a one-step scenario; The matching step is successful; Result is Successful" begin
        given = Given("Some precondition")
        stepdefmatcher = FakeStepDefinitionMatcher(Dict(given => successful_step_definition))
        executor = ExecutableSpecifications.Executor(stepdefmatcher)
        scenario = Scenario("Description", [], [given])

        scenarioresult = ExecutableSpecifications.executescenario(executor, scenario)

        @test isa(scenarioresult.steps[1], ExecutableSpecifications.SuccessfulStepExecution)
    end

    @testset "Execute a one-step scenario; The matching step fails; Result is Failed" begin
        given = Given("Some precondition")
        stepdefmatcher = FakeStepDefinitionMatcher(Dict(given => failed_step_definition))
        executor = ExecutableSpecifications.Executor(stepdefmatcher)
        scenario = Scenario("Description", [], [given])

        scenarioresult = ExecutableSpecifications.executescenario(executor, scenario)

        @test isa(scenarioresult.steps[1], ExecutableSpecifications.StepFailed)
    end

    @testset "Execute a one-step scenario; The matching step throws an error; Result is Error" begin
        given = Given("Some precondition")
        stepdefmatcher = FakeStepDefinitionMatcher(Dict(given => error_step_definition))
        executor = ExecutableSpecifications.Executor(stepdefmatcher)
        scenario = Scenario("Description", [], [given])

        scenarioresult = ExecutableSpecifications.executescenario(executor, scenario)

        @test isa(scenarioresult.steps[1], ExecutableSpecifications.UnexpectedStepError)
    end

    @testset "Execute a two-step scenario; First step throws an error; Second step is Skipped" begin
        given = Given("Some precondition")
        when = When("some action")
        stepdefmatcher = FakeStepDefinitionMatcher(Dict(given => error_step_definition,
                                                        when => successful_step_definition))
        executor = ExecutableSpecifications.Executor(stepdefmatcher)
        scenario = Scenario("Description", [], [given, when])

        scenarioresult = ExecutableSpecifications.executescenario(executor, scenario)

        @test isa(scenarioresult.steps[2], ExecutableSpecifications.SkippedStep)
    end

    @testset "Execute a two-step scenario; First step fails; Second step is Skipped" begin
        given = Given("Some precondition")
        when = When("some action")
        stepdefmatcher = FakeStepDefinitionMatcher(Dict(given => failed_step_definition,
                                                        when => successful_step_definition))
        executor = ExecutableSpecifications.Executor(stepdefmatcher)
        scenario = Scenario("Description", [], [given, when])

        scenarioresult = ExecutableSpecifications.executescenario(executor, scenario)

        @test isa(scenarioresult.steps[2], ExecutableSpecifications.SkippedStep)
    end

    @testset "Execute a two-step scenario; Both steps succeed; All results are Success" begin
        given = Given("Some precondition")
        when = When("some action")
        stepdefmatcher = FakeStepDefinitionMatcher(Dict(given => successful_step_definition,
                                                        when => successful_step_definition))
        executor = ExecutableSpecifications.Executor(stepdefmatcher)
        scenario = Scenario("Description", [], [given, when])

        scenarioresult = ExecutableSpecifications.executescenario(executor, scenario)

        @test isa(scenarioresult.steps[1], ExecutableSpecifications.SuccessfulStepExecution)
        @test isa(scenarioresult.steps[2], ExecutableSpecifications.SuccessfulStepExecution)
    end

    @testset "Execute a three-step scenario; All steps succeeed; All results are Success" begin
        given = Given("Some precondition")
        when = When("some action")
        then = Then("some postcondition")
        stepdefmatcher = FakeStepDefinitionMatcher(Dict(given => successful_step_definition,
                                                        when => successful_step_definition,
                                                        then => successful_step_definition))
        executor = ExecutableSpecifications.Executor(stepdefmatcher)
        scenario = Scenario("Description", [], [given, when, then])

        scenarioresult = ExecutableSpecifications.executescenario(executor, scenario)

        @test isa(scenarioresult.steps[1], ExecutableSpecifications.SuccessfulStepExecution)
        @test isa(scenarioresult.steps[2], ExecutableSpecifications.SuccessfulStepExecution)
        @test isa(scenarioresult.steps[3], ExecutableSpecifications.SuccessfulStepExecution)
    end

    @testset "Execute a scenario; Scenario is provided; Scenario is returned with the result" begin
        given = Given("Some precondition")
        stepdefmatcher = FakeStepDefinitionMatcher(Dict(given => successful_step_definition))
        executor = ExecutableSpecifications.Executor(stepdefmatcher)
        scenario = Scenario("This is a scenario", [], [given])

        scenarioresult = ExecutableSpecifications.executescenario(executor, scenario)

        @test scenarioresult.scenario == scenario
    end

    @testset "Execute a scenario; No unique step definition found; Result is NonUniqueMatch" begin
        stepdefmatcher = ThrowingStepDefinitionMatcher(ExecutableSpecifications.NonUniqueStepDefinition([]))
        executor = ExecutableSpecifications.Executor(stepdefmatcher)
        scenario = Scenario("Description", [], [Given("some precondition")])

        scenarioresult = ExecutableSpecifications.executescenario(executor, scenario)

        @test isa(scenarioresult.steps[1], ExecutableSpecifications.NonUniqueMatch)
    end

    @testset "Execute a ScenarioOutline; Outline has two examples; Two scenarios are returned" begin
        outline = ScenarioOutline("", [], [Given("step <stepnumber>")], ["stepnumber"], ["1" "2"])
        step1 = Given("step 1")
        step2 = Given("step 2")
        stepdefmatcher = FakeStepDefinitionMatcher(Dict(step1 => successful_step_definition,
                                                        step2 => successful_step_definition))
        executor = ExecutableSpecifications.Executor(stepdefmatcher)

        outlineresult = ExecutableSpecifications.executescenario(executor, outline)

        @test length(outlineresult) == 2
    end

    @testset "Execute a ScenarioOutline; Outline has a successful and a failing example; First is success, second is fail" begin
        outline = ScenarioOutline("", [], [Given("step <stepnumber>")], ["stepnumber"], ["1" "2"])
        step1 = Given("step 1")
        step2 = Given("step 2")
        stepdefmatcher = FakeStepDefinitionMatcher(Dict(step1 => successful_step_definition,
                                                        step2 => failed_step_definition))
        executor = ExecutableSpecifications.Executor(stepdefmatcher)

        outlineresult = ExecutableSpecifications.executescenario(executor, outline)

        @test outlineresult[1].steps[1] isa ExecutableSpecifications.SuccessfulStepExecution
        @test outlineresult[2].steps[1] isa ExecutableSpecifications.StepFailed
    end

    @testset "Execute a ScenarioOutline; Outline has three examples; Three scenarios are returned" begin
        outline = ScenarioOutline("", [], [Given("step <stepnumber>")], ["stepnumber"], ["1" "2" "3"])
        step1 = Given("step 1")
        step2 = Given("step 2")
        step3 = Given("step 3")
        stepdefmatcher = FakeStepDefinitionMatcher(Dict(step1 => successful_step_definition,
                                                        step2 => successful_step_definition,
                                                        step3 => successful_step_definition))
        executor = ExecutableSpecifications.Executor(stepdefmatcher)

        outlineresult = ExecutableSpecifications.executescenario(executor, outline)

        @test length(outlineresult) == 3
    end
end

mutable struct FakeRealTimePresenter <: ExecutableSpecifications.RealTimePresenter
    scenarios::Vector{Scenario}
    steps::Vector{ScenarioStep}
    results::Dict{ScenarioStep, StepExecutionResult}
    features::Vector{Feature}

    FakeRealTimePresenter() = new([], [], Dict(), [])
end

present(p::FakeRealTimePresenter, scenario::Scenario) = push!(p.scenarios, scenario)
present(p::FakeRealTimePresenter, step::ScenarioStep) = push!(p.steps, step)
present(p::FakeRealTimePresenter, step::ScenarioStep, result::StepExecutionResult) = p.results[step] = result
present(p::FakeRealTimePresenter, feature::Feature) = push!(p.features, feature)

stepresult(p::FakeRealTimePresenter, step::ScenarioStep) = p.results[step]

@testset "Executor Presentation" begin
    @testset "Execution presentation; Scenario is executed; Scenario is presented" begin
        presenter = FakeRealTimePresenter()
        matcher = FakeStepDefinitionMatcher(Dict())
        executor = Executor(matcher, presenter)

        scenario = Scenario("Some scenario", [], [])
        ExecutableSpecifications.executescenario(executor, scenario)

        @test presenter.scenarios[1] == scenario
    end

    @testset "Execution presentation; Scenario has on Given; Given is presented" begin
        presenter = FakeRealTimePresenter()
        given = Given("some precondition")
        matcher = FakeStepDefinitionMatcher(Dict(given => successful_step_definition))
        executor = Executor(matcher, presenter)

        scenario = Scenario("Some scenario", [], [given])
        ExecutableSpecifications.executescenario(executor, scenario)

        @test presenter.steps[1] == given
    end

    @testset "Execution presentation; Scenario step is successful; Step is presented as successful" begin
        presenter = FakeRealTimePresenter()
        given = Given("some precondition")
        matcher = FakeStepDefinitionMatcher(Dict(given => successful_step_definition))
        executor = Executor(matcher, presenter)

        scenario = Scenario("Some scenario", [], [given])
        ExecutableSpecifications.executescenario(executor, scenario)

        @test stepresult(presenter, given) == ExecutableSpecifications.SuccessfulStepExecution()
    end

    @testset "Execution presentation; Scenario step fails; Next is also presented" begin
        presenter = FakeRealTimePresenter()
        given = Given("some precondition")
        when = When("some action")
        matcher = FakeStepDefinitionMatcher(Dict(given => failed_step_definition,
                                                 when => successful_step_definition))
        executor = Executor(matcher, presenter)

        scenario = Scenario("Some scenario", [], [given, when])
        ExecutableSpecifications.executescenario(executor, scenario)

        @test presenter.steps[2] == when
    end

    @testset "Execution presentation; Scenario step fails; Next has result Skipped" begin
        presenter = FakeRealTimePresenter()
        given = Given("some precondition")
        when = When("some action")
        matcher = FakeStepDefinitionMatcher(Dict(given => failed_step_definition,
                                                 when => successful_step_definition))
        executor = Executor(matcher, presenter)

        scenario = Scenario("Some scenario", [], [given, when])
        ExecutableSpecifications.executescenario(executor, scenario)

        @test stepresult(presenter, when) == ExecutableSpecifications.SkippedStep()
    end

    @testset "Execution presentation; Feature is executed; Feature is presented" begin
        presenter = FakeRealTimePresenter()
        matcher = FakeStepDefinitionMatcher(Dict())
        executor = Executor(matcher, presenter)

        scenario = Scenario("Some scenario", [], [])
        feature = Feature(FeatureHeader("", [], []), [scenario])
        ExecutableSpecifications.executefeature(executor, feature)

        @test presenter.features[1] == feature
    end
end

@testset "Feature Executor" begin
    @testset "Execute a feature; Feature has one scenario; Feature result has one scenario result" begin
        presenter = QuietRealTimePresenter()
        given = Given("some precondition")
        matcher = FakeStepDefinitionMatcher(Dict(given => successful_step_definition))
        scenario = Scenario("some scenario", [], [given])
        featureheader = FeatureHeader("Some feature", [], [])
        feature = Feature(featureheader, [scenario])
        executor = Executor(matcher, presenter)

        featureresult = executefeature(executor, feature)

        @test length(featureresult.scenarioresults) == 1
    end

    @testset "Execute a feature; Feature has two scenarios; Result has two scenario results" begin
        presenter = QuietRealTimePresenter()
        given = Given("some precondition")
        matcher = FakeStepDefinitionMatcher(Dict(given => successful_step_definition))
        scenario1 = Scenario("some scenario", [], [given])
        scenario2 = Scenario("some other scenario", [], [given])
        featureheader = FeatureHeader("Some feature", [], [])
        feature = Feature(featureheader, [scenario1, scenario2])
        executor = Executor(matcher, presenter)

        featureresult = executefeature(executor, feature)

        @test length(featureresult.scenarioresults) == 2
    end

    @testset "Execute a feature; Feature has three scenarios; Scenarios are executed in order" begin
        presenter = QuietRealTimePresenter()
        given = Given("some precondition")
        matcher = FakeStepDefinitionMatcher(Dict(given => successful_step_definition))
        scenario1 = Scenario("some scenario", [], [given])
        scenario2 = Scenario("some other scenario", [], [given])
        scenario3 = Scenario("some third scenario", [], [given])
        featureheader = FeatureHeader("Some feature", [], [])
        feature = Feature(featureheader, [scenario1, scenario2, scenario3])
        executor = Executor(matcher, presenter)

        featureresult = executefeature(executor, feature)

        @test featureresult.scenarioresults[1].scenario == scenario1
        @test featureresult.scenarioresults[2].scenario == scenario2
        @test featureresult.scenarioresults[3].scenario == scenario3
    end

    @testset "Execute a feature; Feature has one failing scenario; Scenario result has a failing step" begin
        presenter = QuietRealTimePresenter()
        given = Given("some precondition")
        matcher = FakeStepDefinitionMatcher(Dict(given => failed_step_definition))
        scenario = Scenario("some scenario", [], [given])
        featureheader = FeatureHeader("Some feature", [], [])
        feature = Feature(featureheader, [scenario])
        executor = Executor(matcher, presenter)

        featureresult = executefeature(executor, feature)

        @test featureresult.scenarioresults[1].steps[1] isa ExecutableSpecifications.StepFailed
    end

    @testset "Execute a feature; One Scenario and an Outline with two examples; Three results" begin
        presenter = QuietRealTimePresenter()

        step1 = Given("step 1")
        step2 = Given("step 2")
        scenario = Scenario("some scenario", [], [step1])
        outline = ScenarioOutline("", [], [Given("step <stepnumber>")], ["stepnumber"], ["1" "2"])
        matcher = FakeStepDefinitionMatcher(Dict(step1 => successful_step_definition,
                                                 step2 => successful_step_definition))

        featureheader = FeatureHeader("Some feature", [], [])
        feature = Feature(featureheader, [scenario, outline])
        executor = Executor(matcher, presenter)

        featureresult = executefeature(executor, feature)

        @test length(featureresult.scenarioresults) == 3
    end
end


text:
export MontyHallEnv

const REWARD_OF_GOAT = 10.0
const REWARD_OF_CAR = 1_000.0

mutable struct MontyHallEnv <: AbstractEnv
    doors::Vector{Symbol}
    rng::AbstractRNG
    guest_action::Union{Nothing,Int}
    host_action::Union{Nothing,Int}
    reward::Union{Nothing,Float64}
end

"""
    MontyHallEnv(;rng=Random.GLOBAL_RNG)

Quoted from [wiki](https://en.wikipedia.org/wiki/Monty_Hall_problem):

> Suppose you're on a game show, and you're given the choice of three doors:
> Behind one door is a car; behind the others, goats. You pick a door, say No.
> 1, and the host, who knows what's behind the doors, opens another door, say
> No. 3, which has a goat. He then says to you, "Do you want to pick door No.
> 2?" Is it to your advantage to switch your choice?

Here we'll introduce the first environment which is of [`FULL_ACTION_SET`](@ref).
"""
function MontyHallEnv(; rng = Random.GLOBAL_RNG)
    doors = fill(:goat, 3)
    doors[rand(rng, 1:3)] = :car
    MontyHallEnv(doors, rng, nothing, nothing, nothing)
end

Random.seed!(env::MontyHallEnv, s) = Random.seed!(env.rng, s)

RLBase.action_space(::MontyHallEnv) = Base.OneTo(3)

"""
In the first round, the guest has 3 options, in the second round only two
options are valid, those different then the host's action.
"""
function RLBase.legal_action_space(env::MontyHallEnv)
    if isnothing(env.host_action)
        1:3
    else
        findall(!=(env.host_action), 1:3)
    end
end

"""
For environments of [`FULL_ACTION_SET`], this function must be implemented.
"""
function RLBase.legal_action_space_mask(env::MontyHallEnv)
    mask = BitArray(undef, 3)
    fill!(mask, true)
    if !isnothing(env.host_action)
        mask[env.host_action] = false
    end
    mask
end

function RLBase.state(env::MontyHallEnv)
    if isnothing(env.host_action)
        1
    else
        env.host_action + 1
    end
end

RLBase.state_space(env::MontyHallEnv) = 1:4

function (env::MontyHallEnv)(action)
    if isnothing(env.host_action)
        # first round
        env.guest_action = action
        if env.doors[action] == :car
            env.host_action = findall(!=(action), 1:3)[rand(env.rng, 1:2)]
        else
            for i in 1:3
                if i == action
                    continue
                elseif env.doors[i] == :goat
                    env.host_action = i
                end
            end
        end
    else
        # second round
        if action == env.host_action
            @error "Invalid action. Can not select the same door with the host."
        else
            env.guest_action = action
            env.reward = env.doors[action] == :goat ? REWARD_OF_GOAT : REWARD_OF_CAR
        end
    end
end

RLBase.reward(env::MontyHallEnv) = isnothing(env.reward) ? 0.0 : env.reward

RLBase.is_terminated(env::MontyHallEnv) = !isnothing(env.reward)

function RLBase.reset!(env::MontyHallEnv)
    env.doors .= :goat
    env.doors[rand(env.rng, 1:3)] = :car
    env.guest_action = nothing
    env.host_action = nothing
    env.reward = nothing
end

RLBase.NumAgentStyle(::MontyHallEnv) = SINGLE_AGENT
RLBase.DynamicStyle(::MontyHallEnv) = SEQUENTIAL
RLBase.ActionStyle(::MontyHallEnv) = FULL_ACTION_SET
RLBase.InformationStyle(::MontyHallEnv) = IMPERFECT_INFORMATION  # the distribution of noise and original reward is unknown to the agent
RLBase.StateStyle(::MontyHallEnv) = Observation{Int}()
RLBase.RewardStyle(::MontyHallEnv) = TERMINAL_REWARD
RLBase.UtilityStyle(::MontyHallEnv) = GENERAL_SUM
RLBase.ChanceStyle(::MontyHallEnv) = STOCHASTIC  # the same action lead to different reward each time.

ext: jl


text:
# indicator of the L0 norm ball with given (integer) radius

"""
  IndBallL0(r::Int=1)

Returns the function `g = ind{x : countnz(x) ⩽ r}`, for an integer parameter `r > 0`.
"""

immutable IndBallL0{I <: Integer} <: IndicatorNonconvex
  r::I
  function IndBallL0(r::I)
    if r <= 0
      error("parameter r must be a positive integer")
    else
      new(r)
    end
  end
end

IndBallL0{I <: Integer}(r::I) = IndBallL0{I}(r)

function (f::IndBallL0){T <: RealOrComplex}(x::AbstractArray{T})
  if countnz(x) > f.r
    return +Inf
  end
  return 0.0
end

function prox!{T <: RealOrComplex}(f::IndBallL0, x::AbstractArray{T}, y::AbstractArray{T}, gamma::Real=1.0)
  p = []
  if ndims(x) == 1
    p = selectperm(x, 1:f.r, by=abs, rev=true)
  else
    p = selectperm(x[:], 1:f.r, by=abs, rev=true)
  end
  sort!(p)
  idx = 1
  for i = 1:length(p)
    y[idx:p[i]-1] = 0.0
    y[p[i]] = x[p[i]]
    idx = p[i]+1
  end
  y[idx:end] = 0.0
  return 0.0
end

fun_name(f::IndBallL0) = "indicator of an L0 pseudo-norm ball"
fun_dom(f::IndBallL0) = "AbstractArray{Real}, AbstractArray{Complex}"
fun_expr(f::IndBallL0) = "x ↦ 0 if countnz(x) ⩽ r, +∞ otherwise"
fun_params(f::IndBallL0) = "r = $(f.r)"

function prox_naive{T <: RealOrComplex}(f::IndBallL0, x::AbstractArray{T}, gamma::Real=1.0)
  p = sortperm(abs.(x)[:], rev=true)
  y = similar(x)
  y[p[1:f.r]] = x[p[1:f.r]]
  y[p[f.r+1:end]] = 0.0
  return y, 0.0
end


text:
include("distribution.jl")

export MI, PI, CMI

ϵ = 0.0000001

# predictive information
PI(data::Vector{Int64}; base=2, mode="ML") = MI(hcat(data[1:end-1], data[2:end]), base=base, mode=mode)

function MI_ML(data::Matrix{Int64}; base=2)
  max = maximum(data)
  px  = fe1p(data[:,1])
  py  = fe1p(data[:,2])
  pxy = fe2p(hcat(data[:,[1:2]]))

  r = 0
  for x=1:length(px)
    for y=1:length(py)
      if px[x] != 0.0 && py[y] != 0.0 && pxy[x,y] > 0.0
        r = r + pxy[x,y] * (log(base, pxy[x,y]) - log(base, px[x]*py[y]))
      end
    end
  end
  r
end

# mutual information
function MI(data::Matrix{Int64}; base=2, mode="ML", pseudocount=0)

  modes      = ["ML", "Maximum Likelihood"]
  umodes     = map(x->uppercase(x), modes)
  known_mode = uppercase(mode) in umodes
  pmodes     = map(x->", $x",modes)
  pmodes     = foldl(*, pmodes[1][3:end], pmodes[2:end])
  @assert known_mode "Mode may be any of the following: [$pmodes]."

  if uppercase(mode) in umodes[1:2]
    MI_ML(data, base=base)
  end

end

# MI(X;Y) = ∑ p(x,y) * log( p(x,y) / (p(x) * p(y)))
function MI(pxy::Matrix{Float64}; base=2)
  px = sum(pxy,1)
  py = sum(pxy,2)
  r  = 0
  for x = 1:size(pxy)[1]
    for y = 1:size(pxy)[2]
      if abs(px[x]) > ϵ && abs(py[y]) > ϵ && abs(pxy[x,y]) > ϵ
        r = r + pxy[x,y] * (log(base, pxy[x,y]) - log(base, px[x] * py[y]))
      end
    end
  end
  return r
end

# MI(X;Y|Z) = ∑ p(x,y,z) * log( p(x,y|z) / (p(x|z) * p(y|z)))
function CMI(pxyz::Array{Float64,3}; base=2)
  r = 0

  pz  = sum(pxyz, (1,2))
  pxz = sum(pxyz, 2)
  pyz = sum(pxyz, 1)

  for x = 1:size(pxyz)[1]
    for y = 1:size(pxyz)[2]
      for z = 1:size(pxyz)[3]
        if abs(pxyz[x,y,z]) > ϵ && abs(pz[1,1,z]) > ϵ && abs(pxz[x,1,z]) > ϵ && abs(pyz[1,y,z]) > ϵ
           r = r + pxyz[x,y,z] * ( log(base, pxyz[x,y,z] / pz[1,1,z]) - log(base, pxz[x,1,z] / pz[1,1,z] * pyz[1,y,z] / pz[1,1,z]))
         end
      end
    end
  end
  return r
end


text:
using Rayuela
using Printf

# Load utilities for saving experimental results
include("experiment_utils.jl")

# === demos query base ===
function run_demos_query_base(
  dataset_name="labelme",     # Dataset to work on
  ntrain::Integer=Int(20e3),  # Number of training/base examples
  m::Integer=8,               # Number of codebooks (-1 will be substracted for non-orthogonal methods)
  h::Integer=256,             # Number of Entries per codebook
  niter::Integer=25,          # Number of training iterations
  ntrials::Integer=10)        # Number of times to repeat the experiment

  nquery, nbase, knn = 0, 0, 0
  if dataset_name == "MNIST"
    nquery, nbase, knn = Int(10e3), Int(60e3), Int(1e3)
  elseif dataset_name == "labelme"
    nquery, nbase, knn = Int(2e3), Int(20019), Int(1e3)
  else
    error("dataset unknown")
  end

  verbose = true

  Xt, Xb, Xq, gt = load_experiment_data(dataset_name, ntrain, nbase, nquery, verbose)
  d, _ = size( Xt )

  # (Semi-)orthogonal methods: PQ, OPQ, ChainQ
  for trial = 1:ntrials
    C, B, train_error, recall = Rayuela.experiment_pq_query_base(Xt, Xq, gt, m, h, niter, knn, verbose)
    save_results_pq_query_base("./results/$(lowercase(dataset_name))/pq_m$(m)_it$(niter).h5", trial, C, B, train_error, recall)
  end
  for trial = 1:ntrials
    init = "natural"
    C, B, R, train_error, recall = Rayuela.experiment_opq_query_base(Xt, Xq, gt, m, h, init, niter, knn, verbose)
    save_results_opq_query_base("./results/$(lowercase(dataset_name))/opq_m$(m)_it$(niter).h5", trial, C, B, R, train_error, recall)
  end

  # Cheap non-orthogonal methods: RVQ, ERVQ
  for trial = 1:ntrials
    C, B, train_error, recall = Rayuela.experiment_rvq_query_base(Xt,        Xq, gt, m-1, h, niter, knn, verbose)
    save_results_pq_query_base("./results/$(lowercase(dataset_name))/rvq_m$(m-1)_it$(niter).h5",  trial, C, B, train_error, recall)
  end
  for trial = 1:ntrials
    C, B, _ = load_rvq("./results/$(lowercase(dataset_name))/rvq_m$(m-1)_it$(niter).h5", m-1, trial)
    C, B, train_error, recall = Rayuela.experiment_ervq_query_base(Xt, B, C, Xq, gt, m-1, h, niter, knn, verbose)
    save_results_pq_query_base("./results/$(lowercase(dataset_name))/ervq_m$(m-1)_it$(niter).h5", trial, C, B, train_error, recall)
  end

  # Precompute init for LSQ/SR
  for trial = 1:ntrials
    C, B, R, train_error = Rayuela.train_opq(Xt, m-1, h, niter, "natural", verbose)
    save_results_opq_query_base("./results/$(lowercase(dataset_name))/opq_m$(m-1)_it$(niter).h5", trial, C, B, R, train_error, [0f0])
  end
  for trial = 1:ntrials
    C, B, R, _ = load_chainq("./results/$(lowercase(dataset_name))/opq_m$(m-1)_it$(niter).h5", m-1, trial)
    C, B, R, chainq_error = Rayuela.train_chainq(Xt, m-1, h, R, B, C, niter, verbose)
    save_results_opq_query_base("./results/$(lowercase(dataset_name))/chainq_m$(m-1)_it$(niter).h5", trial, C, B, R, chainq_error, [0f0])
  end

  ilsiter = 8
  icmiter = 4
  randord = true
  npert   = 4
  nsplits_train = 1

  for trial = 1:ntrials
    C, B, R, chainq_error = load_chainq("./results/$(lowercase(dataset_name))/chainq_m$(m-1)_it$(niter).h5", m-1, trial)
    C, B, R, train_error, recall = Rayuela.experiment_lsq_cuda_query_base(Xt, B, C, R, Xq, gt, m-1, h,
      niter, ilsiter, icmiter, randord, npert, knn, nsplits_train, verbose)
    save_results_lsq_query_base("./results/$(lowercase(dataset_name))/lsq_m$(m-1)_it$(niter).h5", trial, C, B, R, train_error, chainq_error, recall)
  end

  schedule = 1
  p = 0.5

  for trial = 1:ntrials
    sr_method = "SR_D"
    C, B, R, chainq_error = load_chainq("./results/$(lowercase(dataset_name))/chainq_m$(m-1)_it$(niter).h5", m-1, trial)
    C, B, R, train_error, recall = Rayuela.experiment_sr_cuda_query_base( Xt, B, C, R, Xq, gt, m-1, h,
      niter, ilsiter, icmiter, randord, npert, knn, nsplits_train, sr_method, schedule, p, verbose)
    save_results_lsq_query_base("./results/$(lowercase(dataset_name))/srd_m$(m-1)_it$(niter).h5", trial, C, B, R, train_error, chainq_error, recall)
  end

  for trial = 1:ntrials
    sr_method = "SR_C"
    C, B, R, chainq_error = load_chainq("./results/$(lowercase(dataset_name))/chainq_m$(m-1)_it$(niter).h5", m-1, trial)
    C, B, R, train_error, recall = Rayuela.experiment_sr_cuda_query_base( Xt, B, C, R, Xq, gt, m-1, h,
      niter, ilsiter, icmiter, randord, npert, knn, nsplits_train, sr_method, schedule, p, verbose)
    save_results_lsq_query_base("./results/$(lowercase(dataset_name))/src_m$(m-1)_it$(niter).h5", trial, C, B, R, train_error, chainq_error, recall)
  end

end

for niter = [25]
  run_demos_query_base("labelme", Int(20e3), 8, 256, niter)
  # run_demos_query_base("labelme", Int(20e3), 16, 256, niter)
  # run_demos_query_base("MNIST",   Int(60e3), 8, 256, niter)
  # run_demos_query_base("MNIST",   Int(60e3), 16, 256, niter)
end


text:
module OrekitWrapper

using JavaCall

include("AttitudesWrapper/AttitudesWrapper.jl")
include("BodiesWrapper/BodiesWrapper.jl")
include("CompilerWrapper/CompilerWrapper.jl")
include("DataWrapper/DataWrapper.jl")
include("ErrorsWrapper/ErrorsWrapper.jl")
include("EstimationWrapper/EstimationWrapper.jl")
include("FilesWrapper/FilesWrapper.jl")
include("ForcesWrapper/ForcesWrapper.jl")
include("FramesWrapper/FramesWrapper.jl")
include("GeometryWrapper/GeometryWrapper.jl")
include("GnssWrapper/GnssWrapper.jl")
include("ModelsWrapper/ModelsWrapper.jl")
include("OrbitsWrapper/OrbitsWrapper.jl")
include("PropagationWrapper/PropagationWrapper.jl")
include("TimeWrapper/TimeWrapper.jl")
include("UtilsWrapper/UtilsWrapper.jl")

end


text:
using CPLEX
using Base.Test

@testset "Changing variable type" begin
    # linear -> mixed integer -> linear
    m = CPLEX.CplexMathProgModel(CPX_PARAM_SCRIND=0)
    CPLEX.loadproblem!(m, [1 1; 1 0], [0, 0], [Inf, 1], [1, 2], [1.33, -Inf], [Inf, 0.5], :Min)

    @test CPLEX.get_prob_type(m.inner) == :LP
    CPLEX.optimize!(m)
    @test CPLEX.getsolution(m) ≈ [0.5, 0.83]

    CPLEX.setvartype!(m, [:Cont, :Bin])

    @test CPLEX.get_prob_type(m.inner) == :MILP
    CPLEX.optimize!(m)
    @test CPLEX.getsolution(m) ≈ [0.33, 1]

    CPLEX.setvartype!(m, [:Cont, :Cont])

    @test CPLEX.get_prob_type(m.inner) == :LP
    CPLEX.optimize!(m)
    @test CPLEX.getsolution(m) ≈ [0.5, 0.83]
end


text:
export find_set_bits,  gray_to_binary, binary_to_gray

function find_set_bits!{T<:Number}(S::Array{T, 1}, code::Array{Bool, 1})
  ns = 0
  n = length(code)
  @inbounds begin
    for i=1:n
      if code[i] == true
        ns = ns + 1
        S[ns] = i
        #      println("$ns -> $S")
      end
    end
  end
  #  println("\tfind_set_bits!\n\t\tcode: $code\n\t\tS: $S\n\t\tns: $ns")
  return ns
end

function find_set_bits(code::String, set_char::Char ='1')
  n = length(code)
  S = zeros(Integer, n, 1)
  c = 0
  @inbounds begin
    for i=[1:length(code)]
      if code[i] == set_char
        c = c + 1
        S[c] = i
      end
    end
  end
  V = S[1:c]
  return V
end

function binary_to_gray!(gray_code::Array{Bool, 1}, binary_code::Array{Bool, 1})
  n = length(binary_code)
  @inbounds begin
    for i=(n-1):-1:1
      if (binary_code[i] != binary_code[i+1])
        gray_code[i] = true
      else
        gray_code[i] = false
      end
    end
  end
  gray_code[n] = binary_code[n]
  #  println("\tbinary: $binary_code\n\tgray: $gray_code")
end

function int_to_binary_bool!(binary_code::Array{Bool, 1}, d::Integer, n::Integer)
  i = d::Integer # local copy
  for j=[1:n;]
    if i > 0
      binary_code[j] = i % 2
      i = i >> 1
    else
      binary_code[j] = false
    end
  end
  #  println("int_to_binary_bool! i: $d n: $n code: $binary_code")
end

function binary_to_gray(binary_code::String)
  n= length(binary_code)
  gray_code  = ""
  for i=[n:-1:2;]
    # @inbounds c_char = binary_code[(i-1):i]
    # println(n, "-> " , i, " -> ", c_char, " : ", gray_code)
    if binary_code[i-1] != binary_code[i]
      gray_code = "1" * gray_code
    else
      gray_code =  "0" * gray_code
    end
  end
  gray_code =  string(binary_code[1]) * gray_code
  return gray_code
end

function gray_to_binary(gray_code::String)
  n = length(gray_code)
  binary_code = string(gray_code[1])
  for i=2:n
    if gray_code[i] == binary_code[i-1]
      binary_code = binary_code * "0"
    else
      binary_code = binary_code * "1"
    end
  end
  return binary_code
end


text:
using SimpleTraits.BaseTraits
using Compat: view

@test istrait(IsAnything{Any})
@test istrait(IsAnything{Union{}})
@test istrait(IsAnything{Int})

@test !istrait(IsNothing{Any})
@test !istrait(IsNothing{Union{}})
@test !istrait(IsNothing{Int})

@test istrait(IsBits{Int})
@test !istrait(IsBits{Vector{Int}})

@test istrait(IsImmutable{Float64})
@test !istrait(IsImmutable{Vector{Int}})

@test istrait(IsCallable{Function})

if VERSION>v"0.4-" # use @generated functions
    a = collect(1:5)
    b = view(a, 2:3)
    c = view(a, 1:2:5)
    @test istrait(IsContiguous{typeof(b)})
    @test !istrait(IsContiguous{typeof(c)})

    @test istrait(IsFastLinearIndex{Vector})
    @test !istrait(IsFastLinearIndex{AbstractArray})

    if VERSION < v"0.5.0-dev"
        @test istrait(IsCallable{Base.AddFun})
    end
end


text:
"""
    setup_dump_aux_diagnostics(
        ::ClimateMachineConfigType,
        interval::String,
        out_prefix::String;
        writer = NetCDFWriter(),
        interpol = nothing,
    )

Create and return a `DiagnosticsGroup` containing a diagnostic that
simply dumps the auxiliary state variables at the specified
`interval` after being interpolated, into NetCDF files prefixed by
`out_prefix`.
"""
function setup_dump_aux_diagnostics(
    ::ClimateMachineConfigType,
    interval::String,
    out_prefix::String;
    writer = NetCDFWriter(),
    interpol = nothing,
)
    # TODO: remove this
    @assert !isnothing(interpol)

    return DiagnosticsGroup(
        "DumpAux",
        Diagnostics.dump_aux_init,
        Diagnostics.dump_aux_fini,
        Diagnostics.dump_aux_collect,
        interval,
        out_prefix,
        writer,
        interpol,
    )
end

dump_aux_init(dgngrp, currtime) = dump_init(dgngrp, currtime, Auxiliary())

function dump_aux_collect(dgngrp, currtime)
    interpol = dgngrp.interpol
    mpicomm = Settings.mpicomm
    dg = Settings.dg
    Q = Settings.Q
    FT = eltype(Q.data)
    bl = dg.balance_law
    mpirank = MPI.Comm_rank(mpicomm)

    iaux = similar(
        dg.state_auxiliary.data,
        interpol.Npl,
        number_states(bl, Auxiliary()),
    )

    interpolate_local!(interpol, dg.state_auxiliary.data, iaux)

    all_aux_data = accumulate_interpolated_data(mpicomm, interpol, iaux)

    if mpirank == 0
        auxnames = flattenednames(vars_state(bl, Auxiliary(), FT))
        varvals = OrderedDict()
        for (vari, varname) in enumerate(auxnames)
            varvals[varname] = all_aux_data[:, :, :, vari]
        end
        append_data(dgngrp.writer, varvals, currtime)
    end

    MPI.Barrier(mpicomm)
    return nothing
end

function dump_aux_fini(dgngrp, currtime) end

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
# # 205: 2D Nonlinear Poisson equation
# ([source code](SOURCE_URL))

module Example205_NonlinearPoisson2D

using Printf
using VoronoiFVM
using ExtendableGrids
using ExtendableSparse
using GridVisualize

function main(;n=10,Plotter=nothing,verbose=false, unknown_storage=:sparse, max_lureuse=0,
              factorization=LUFactorization())
    h=1.0/convert(Float64,n)
    X=collect(0.0:h:1.0)
    Y=collect(0.0:h:1.0)

    grid=VoronoiFVM.Grid(X,Y)

    eps=1.0e-2

    physics=VoronoiFVM.Physics(
        num_species=1,
        reaction=function(f,u,node)
        f[1]=u[1]^2
        end,

        flux=function(f,u0,edge)
        u=unknowns(edge,u0)
        f[1]=eps*(u[1,1]^2-u[1,2]^2)
        end,

        source=function(f,node)
        x1=node[1]-0.5
        x2=node[2]-0.5
        f[1]=exp(-20.0*(x1^2+x2^2))
        end,

        storage=function(f,u,node)
        f[1]=u[1]
        end)
    sys=VoronoiFVM.System(grid,physics,unknown_storage=unknown_storage)
    enable_species!(sys,1,[1])

    boundary_dirichlet!(sys,1,2,0.1)
    boundary_dirichlet!(sys,1,4,0.1)

    inival=unknowns(sys)
    U=unknowns(sys)
    inival.=0.5

    control=VoronoiFVM.NewtonControl()
    control.verbose=verbose
    control.tol_linear=1.0e-5
    control.max_lureuse=max_lureuse
    control.factorization=factorization
    tstep=0.01
    time=0.0
    u15=0
    p=GridVisualizer(Plotter=Plotter)
    while time<1.0
        time=time+tstep
        solve!(U,inival,sys,control=control,tstep=tstep)
        u15=U[15]
        inival.=U

        if verbose
            @printf("time=%g\n",time)
        end

        scalarplot!(p[1,1],grid,U[1,:],Plotter=Plotter,clear=true,show=true)
        tstep*=1.0
    end
    return u15
end

function test()
    # test at once for iterative solution here
    testval=0.3554284760906605
    main(unknown_storage=:sparse,max_lureuse=0) ≈  testval &&
        main(unknown_storage=:dense,max_lureuse=0) ≈ testval &&
        main(unknown_storage=:sparse,max_lureuse=10) ≈ testval &&
        main(unknown_storage=:dense,max_lureuse=10) ≈ testval &&
        main(unknown_storage=:sparse,max_lureuse=0, factorization=ILU0Preconditioner()) ≈ testval &&
        main(unknown_storage=:dense,max_lureuse=0, factorization=ILU0Preconditioner()) ≈ testval
end
end


text:
using Oceananigans.Grids: total_length, topology

using OffsetArrays: OffsetArray

#####
##### Creating offset arrays for field data by dispatching on architecture.
#####

"""
Return a range of indices for a field located at `Center` centers
`along a grid dimension of length `N` and with halo points `H`.
"""
offset_indices(loc, topo, N, H=0) = 1 - H : N + H

"""
Return a range of indices for a field located at cell `Face`s
`along a grid dimension of length `N` and with halo points `H`.
"""
offset_indices(::Type{Face}, ::Type{Bounded}, N, H=0) = 1 - H : N + H + 1

"""
Return a range of indices for a field along a 'reduced' dimension.
"""
offset_indices(::Type{Nothing}, topo, N, H=0) = 1 : 1

"""
    offset_data(underlying_data, grid::AbstractGrid, loc)

Returns an `OffsetArray` that maps to `underlying_data` in memory,
with offset indices appropriate for the `data` of a field on
a `grid` of `size(grid)` and located at `loc`.
"""
function offset_data(underlying_data, grid::AbstractGrid, loc)
    ii = offset_indices(loc[1], topology(grid, 1), grid.Nx, grid.Hx)
    jj = offset_indices(loc[2], topology(grid, 2), grid.Ny, grid.Hy)
    kk = offset_indices(loc[3], topology(grid, 3), grid.Nz, grid.Hz)

    return OffsetArray(underlying_data, ii, jj, kk)
end

"""
    new_data([FT=Float64], ::CPU, grid, loc)

Returns an `OffsetArray` of zeros of float type `FT`, with
parent data in CPU memory and indices corresponding to a field on a
`grid` of `size(grid)` and located at `loc`.
"""
function new_data(FT, ::AbstractCPUArchitecture, grid, loc)
    underlying_data = zeros(FT, total_length(loc[1], topology(grid, 1), grid.Nx, grid.Hx),
                                total_length(loc[2], topology(grid, 2), grid.Ny, grid.Hy),
                                total_length(loc[3], topology(grid, 3), grid.Nz, grid.Hz))

    return offset_data(underlying_data, grid, loc)
end

"""
    new_data([FT=Float64], ::GPU, grid, loc)

Returns an `OffsetArray` of zeros of float type `FT`, with
parent data in GPU memory and indices corresponding to a field on a `grid`
of `size(grid)` and located at `loc`.
"""
function new_data(FT, ::AbstractGPUArchitecture, grid, loc)
    underlying_data = CuArray{FT}(undef, total_length(loc[1], topology(grid, 1), grid.Nx, grid.Hx),
                                         total_length(loc[2], topology(grid, 2), grid.Ny, grid.Hy),
                                         total_length(loc[3], topology(grid, 3), grid.Nz, grid.Hz))

    underlying_data .= 0 # Ensure data is initially 0.

    return offset_data(underlying_data, grid, loc)
end

# Default to type of Grid
new_data(arch, grid, loc) = new_data(eltype(grid), arch, grid, loc)


text:
#####
##### instantiate
#####

export instantiate

"""
    instantiate(s::T,
                D::Dict,
                ft::AbstractFilterTypes,
                stop_recursion_type=FreeParameter) where {T}

Instantiate instance of struct `s` and replace
using values in `D` where filter types `ft` apply.
"""
function instantiate(s::T,
                     D::Dict,
                     ft::AbstractFilterTypes) where {T}
  is_leaf(s, prop) = __is_leaf(s, prop, ft)
  _basename = first(split(first(keys(D)), "."))
  return _instantiate(s, D, is_leaf, _basename)
end

function _instantiate(s::T, D::Dict, is_leaf::F, _basename::S) where {T,F<:Function,S<:AbstractString}
  cond_s = is_leaf(nothing, s)
  if T <: Function
    return s
  elseif triggers(cond_s) || T<:FPTYPES
    s_active = haskey(D, _basename) ? D[_basename] : s
    cond_D = is_leaf(nothing, s_active)
    cond = haskey(D, _basename) && aux_leaf(cond_D)
    return cond ? get_val(D[_basename]) : s
  else
    props_new = [_instantiate(getproperty(s, fn), D, is_leaf, _basename*"."*string(fn)) for fn in fieldnames(T)]
    return T <: Tuple ? T(props_new) : T(props_new...)
  end
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Test, LinearAlgebra
@test detect_ambiguities(LinearAlgebra; imported=true, recursive=true) == []

max_stars_repo_name: greimel/julia

max_issues_repo_name: greimel/julia

max_forks_repo_name: greimel/julia


text:
using BinaryProvider # requires BinaryProvider 0.3.0 or later

# Parse some basic command-line arguments
const verbose = "--verbose" in ARGS
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))
products = [
    LibraryProduct(prefix, ["libCbcSolver"], :libcbcsolver),
    LibraryProduct(prefix, ["libCbc"], :libCbc),
]

## START-VERSION-UPDATE-BLOCK
# Download binaries from hosted location
bin_prefix = "https://github.com/JuliaOpt/CbcBuilder/releases/download/v2.10.3-static-mono-small"
# Listing of files generated by BinaryBuilder:
download_info = Dict(
    Linux(:aarch64, libc=:glibc, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/CbcBuilder.v2.10.3.aarch64-linux-gnu-gcc7.tar.gz", "2f4f8470d5ad1ba94c5456e534acc7dc7b39a65bb8fe80c574d57bb7a09d0a33"),
    Linux(:aarch64, libc=:glibc, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/CbcBuilder.v2.10.3.aarch64-linux-gnu-gcc8.tar.gz", "0d5f8aa1bd29ad1948c32aedbb802e0e3c8ad9fdbece857cc5ec554e37d59592"),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/CbcBuilder.v2.10.3.arm-linux-gnueabihf-gcc7.tar.gz", "317b915dc88ae80946a024e43b7005bc4c355fb6bded736b5e3d3905fc1181da"),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/CbcBuilder.v2.10.3.arm-linux-gnueabihf-gcc8.tar.gz", "0deae3b140379380df1eb91f1e0dfdedd523fc5bedaa905d84a566e2c2d7394a"),
    Linux(:i686, libc=:glibc, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/CbcBuilder.v2.10.3.i686-linux-gnu-gcc4.tar.gz", "b3ae1aa3eb0c4e857c85c03d7b5f47918c89251e2919b3a20c8325ee0c218fbe"),
    Linux(:i686, libc=:glibc, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/CbcBuilder.v2.10.3.i686-linux-gnu-gcc7.tar.gz", "efa9f5a1a177d227de847652ace146e739be1f25d7f9be3e8d86f39cf2c61a05"),
    Linux(:i686, libc=:glibc, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/CbcBuilder.v2.10.3.i686-linux-gnu-gcc8.tar.gz", "ca250bbcc505a228f769d2b9ec96e04537c2a49420f50bf744eface5c988fb53"),
    Windows(:i686, compiler_abi=CompilerABI(:gcc6)) => ("$bin_prefix/CbcBuilder.v2.10.3.i686-w64-mingw32-gcc6.tar.gz", "77ced49c1c47fecfe8eb7b6f4e34a4a5550764d7174b03dc8f252578d681d633"),
    Windows(:i686, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/CbcBuilder.v2.10.3.i686-w64-mingw32-gcc7.tar.gz", "49c443e52a02367a2759867c387023ce031b046b063287b927125ea39e6d9b30"),
    Windows(:i686, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/CbcBuilder.v2.10.3.i686-w64-mingw32-gcc8.tar.gz", "f726efe558ec9ac6137c8ab68e8fd4aa1238caa56f7a1bb3aaf6d539b7d145ac"),
    MacOS(:x86_64, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/CbcBuilder.v2.10.3.x86_64-apple-darwin14-gcc4.tar.gz", "c0487ea6c0a6ef80b3b62a94470d943f7eaea0cf7e69ef8532dbca288d182c18"),
    MacOS(:x86_64, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/CbcBuilder.v2.10.3.x86_64-apple-darwin14-gcc7.tar.gz", "800ef428b72a51cd602c45ea5f33c437379369411ebae749941b9a18010a07eb"),
    MacOS(:x86_64, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/CbcBuilder.v2.10.3.x86_64-apple-darwin14-gcc8.tar.gz", "0222b3178cefe53d416e04169185bfce34808a4f8940b3cebcfcaed8b055e594"),
    Linux(:x86_64, libc=:glibc, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/CbcBuilder.v2.10.3.x86_64-linux-gnu-gcc4.tar.gz", "48a97b0109b3b773757fc83e8c32e5a67612fbaf03f1d15f9cc9a2158845d9a8"),
    Linux(:x86_64, libc=:glibc, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/CbcBuilder.v2.10.3.x86_64-linux-gnu-gcc7.tar.gz", "86cddcacbf3133e6249da8c795bc1216ef425e585a446e00fec3856678204367"),
    Linux(:x86_64, libc=:glibc, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/CbcBuilder.v2.10.3.x86_64-linux-gnu-gcc8.tar.gz", "6be3a422be4f6fd4d382df0b82a50c4b54473ea59c66d775c4713e1ed7adc67d"),
    Windows(:x86_64, compiler_abi=CompilerABI(:gcc6)) => ("$bin_prefix/CbcBuilder.v2.10.3.x86_64-w64-mingw32-gcc6.tar.gz", "fdb78659594db6f69b96d11dd48d83973da13cdaf72f5a69afbddc1d535c66b2"),
    Windows(:x86_64, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/CbcBuilder.v2.10.3.x86_64-w64-mingw32-gcc7.tar.gz", "aae75d7fe4f34da22233ca9ac0248f00c192be27a3ed43e7a1f8e921769e32a8"),
    Windows(:x86_64, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/CbcBuilder.v2.10.3.x86_64-w64-mingw32-gcc8.tar.gz", "123091adf7806c9c3748c4ceef0fe3929737340d8ad0e743b3b8c80dcf53b494"),
)
## END-VERSION-UPDATE-BLOCK

# To fix gcc4 bug in Windows
# https://sourceforge.net/p/mingw-w64/bugs/727/
this_platform = platform_key_abi()
if typeof(this_platform)==Windows && this_platform.compiler_abi.gcc_version == :gcc4
   this_platform = Windows(arch(this_platform), libc=libc(this_platform), compiler_abi=CompilerABI(:gcc6))
end

# no dynamic dependencies until Pkg3 support for binaries
dependencies = [
#     "https://github.com/juan-pablo-vielma/CglBuilder/releases/download/v0.59.10-1/build_CglBuilder.v0.59.10.jl",
#     "https://github.com/JuliaOpt/ClpBuilder/releases/download/v1.16.11-1/build_ClpBuilder.v1.16.11.jl",
#     "https://github.com/juan-pablo-vielma/OsiBuilder/releases/download/v0.107.9-1/build_OsiBuilder.v0.107.9.jl",
#     "https://github.com/juan-pablo-vielma/CoinUtilsBuilder/releases/download/v2.10.14-1/build_CoinUtilsBuilder.v2.10.14.jl",
#     "https://github.com/juan-pablo-vielma/COINMumpsBuilder/releases/download/v1.6.0-1/build_COINMumpsBuilder.v1.6.0.jl",
#     "https://github.com/juan-pablo-vielma/COINMetisBuilder/releases/download/v1.3.5-1/build_COINMetisBuilder.v1.3.5.jl",
#     "https://github.com/juan-pablo-vielma/COINLapackBuilder/releases/download/v1.5.6-1/build_COINLapackBuilder.v1.5.6.jl",
#     "https://github.com/juan-pablo-vielma/COINBLASBuilder/releases/download/v1.4.6-1/build_COINBLASBuilder.v1.4.6.jl",
#     "https://github.com/juan-pablo-vielma/ASLBuilder/releases/download/v3.1.0-1/build_ASLBuilder.v3.1.0.jl"
]

custom_library = false
if haskey(ENV,"JULIA_CBC_LIBRARY_PATH")
    custom_products = [LibraryProduct(ENV["JULIA_CBC_LIBRARY_PATH"],product.libnames,product.variable_name) for product in products]
    if all(satisfied(p; verbose=verbose) for p in custom_products)
        products = custom_products
        custom_library = true
    else
        error("Could not install custom libraries from $(ENV["JULIA_CBC_LIBRARY_PATH"]).\nTo fall back to BinaryProvider call delete!(ENV,\"JULIA_CBC_LIBRARY_PATH\") and run build again.")
    end
end

if !custom_library
    # Install unsatisfied or updated dependencies:
    unsatisfied = any(!satisfied(p; verbose=verbose) for p in products)

    dl_info = choose_download(download_info, this_platform)
    if dl_info === nothing && unsatisfied
        # If we don't have a compatible .tar.gz to download, complain.
        # Alternatively, you could attempt to install from a separate provider,
        # build from source or something even more ambitious here.
        error("Your platform (\"$(Sys.MACHINE)\", parsed as \"$(triplet(platform_key_abi()))\") is not supported by this package!")
    end

    # If we have a download, and we are unsatisfied (or the version we're
    # trying to install is not itself installed) then load it up!
    if unsatisfied || !isinstalled(dl_info...; prefix=prefix)
        # Download and install binaries
        # no dynamic dependencies until Pkg3 support for binaries
        # for dependency in reverse(dependencies)          # We do not check for already installed dependencies
        #    download(dependency,basename(dependency))
        #    evalfile(basename(dependency))
        # end
        install(dl_info...; prefix=prefix, force=true, verbose=verbose)
    end
 end

# Write out a deps.jl file that will contain mappings for our products
write_deps_file(joinpath(@__DIR__, "deps.jl"), products, verbose=verbose)


text:
using Plots, StatsBase, Statistics, LaTeXStrings, JLD, StatsPlots

function FindBorder(dim, Network)
    border0 = []
    for i in 1:dim
        for j in 1:dim
            check = false
            if Network[i,j] == 0
                Neighbors = ((i-1,j),(i,j-1),(i+1,j),(i,j+1))
                for nei in Neighbors
                    if dim+1 ∉ nei && 0 ∉ nei && Network[nei...] == 1
                        check = true
                        break
                    end
                end
                if check == true
                    push!(border0,(i,j))
                end
            end
        end
    end
    return border0
end

function FindCorrelationLength(Network,dim)
    FindOnSpots = findall(x->x==1,Network)
    SpotNums = length(FindOnSpots)
    ΣR² = 0.0; iMC = 0.0; jMC = 0.0
    for indx in FindOnSpots
        iMC += indx[1]/SpotNums
        jMC += indx[2]/SpotNums
    end
    for indx in FindOnSpots
        ΣR² += ((indx[1]-iMC)^2 + (indx[2]-jMC)^2)/SpotNums
    end
    return sqrt(ΣR²)
end

function ClusterGrowth(dim,P)
    Network = zeros(Int,dim,dim)
    # Network[rand(1:dim^2)] = 1
    Network[round(Int, dim/2), round(Int, dim/2)] = 1
    Borders = FindBorder(dim,Network)
    while length(Borders) != 0
        Borders = FindBorder(dim,Network)
        for spot in Borders
            if P > rand()
                Network[spot...] = 1
            else
                Network[spot...] = -1
            end
        end
    end
    return Network
end

function Linear_fit(X, Y)
    A = [hcat(X) reshape(ones(length(X)), length(X), 1)]
    b = reshape(hcat(Y), length(Y), 1)
    line = (A \ b)
    return line
end

# Network = ClusterGrowth(200,0.56)
# heatmap(Network, legend = nothing)

# runnum = 10000
# for p in [0.5,0.55,0.59]
#     totstepdata = []
#     for i in 1:runnum
#         Network = ClusterGrowth(200,p)
#         push!(totstepdata,(FindCorrelationLength(Network,200),length(findall(x-> x==1,Network))))
#         print("\r$i")
#     end
#     println('#')
#     save("../../Data/Q3/Q3-$p-CG.jld", "data", totstepdata)
# end

# AllXiData = []
# AllSData = []
# for p in [0.5,0.55,0.59]
#     Data = load("../../Data/Q3/Q3-$p-CG.jld")["data"]
#     xidata = []
#     sdata = []
#     for data in Data
#         push!(xidata, data[1])
#         push!(sdata, data[2])
#     end
#     push!(AllXiData,xidata)
#     push!(AllSData,sdata)
# end
# AllXiData

#Ploting:
#=
violin(["0.59"],AllXiData[3])
violin!(["0.55"],AllXiData[2])
violin!(["0.50"],AllXiData[1])
plot!(legend = nothing, title = L"Correlation\ Length\ for\ 200\times200\ Network,\ 10000 runs",xlabel = L"P",ylabel=L"\xi_{(P)}")
savefig("../../Figs/Q3/Q3-xi-violon.pdf")

histogram(AllXiData[3],bins = 50, label = L"P=0.59")
histogram!(AllXiData[2],bins = 50, label = L"P=0.55")
histogram!(AllXiData[1],bins = 25, label = L"P=0.50")
plot!(title = L"Correlation\ Length\ for\ 200\times200\ Network,\ 10000 runs",xlabel = L"\xi_{(P)}",ylabel=L"Numbers")
savefig("../../Figs/Q3/Q3-xi-hist.pdf")

violin(["0.59"],AllSData[3])
violin!(["0.55"],AllSData[2])
violin!(["0.50"],AllSData[1])
plot!(legend = nothing, title = L"Cluster\ Size\ for\ 200\times200\ Network,\ 10000 runs",xlabel = L"P",ylabel=L"\S_{(P)}")
savefig("../../Figs/Q3/Q3-S-violon.pdf")
=#

histogram(AllXiData[3],bins = 50, label = L"P=0.59")
histogram!(AllXiData[2],bins = 50, label = L"P=0.55")
histogram!(AllXiData[1],bins = 25, label = L"P=0.50")
plot!(title = L"Cluster\ Size\ for\ 200\times200\ Network,\ 10000 runs",xlabel = L"\S_{(P)}",ylabel=L"Numbers")
savefig("../../Figs/Q3/Q3-S-hist.pdf")

runnum = 10000
AllData = []
for p in 0.5:0.005:0.65
    totstepdata = []
    for i in 1:runnum
        Network = ClusterGrowth(100,p)
        push!(totstepdata,(FindCorrelationLength(Network,100),length(findall(x-> x==1,Network))))
        print("\r$i--$p    ")
    end
    push!(AllData, totstepdata)
end
save("../../Data/Q3/Q3-xi-S-CG.jld", "data", AllData)

AllData
AvgSList

AvgSList = []
STDSList = []
AvgXiList = []
STDXiList = []
for data in AllData
    temspdata = []
    temxipdata = []
    for tup in data
        push!(temxipdata,tup[1])
        push!(temspdata,tup[2])
    end
    push!(AvgSList,mean(temspdata))
    push!(STDSList,std(temspdata))
    push!(AvgXiList,mean(temxipdata))
    push!(STDXiList,std(temxipdata))
end

AllData = load("../../Data/Q3/Q3-xi-S-CG.jld")["data"]

XData = []
YData = []
for data in AllData
    for tup in data
        push!(XData,tup[1])
        push!(YData,tup[2])
    end
end

scatter(log.(XData), log.(YData),markersize = 2,
    alpha = 0.01, markerstrokewidth=0)
scatter!(title = L"Total\ Data\ (0.5>P>0.65,\ 10000\ runs\ for\ each\ P)",
    legend = nothing, ylabel = L"\xi_{(P)}", xlabel = L"S_{(P)}")
savefig("../../Figs/Q3/Q3-S-XI.pdf")

Line = Linear_fit(log.(AvgXiList)[1:end-3],log.(AvgSList)[1:end-3])
X = hcat(log.(AvgXiList)[1]:(log.(AvgXiList)[end]-log.(AvgXiList)[1])/20:log.(AvgXiList)[end])
Y = X .* Line[1,1] .+ Line[2,1]

plot(X,Y, label = L"Y = %$(round(Line[1],digits=3))\dot X + %$(round(Line[2],digits=3))", line = :dash, c = :black)
scatter!(log.(AvgXiList), log.(AvgSList),c = :steelblue, label = L"Data\ Point")
scatter!(title = L"Average\ Data\ (0.5>P>0.65,\ 10000\ runs\ for\ each\ P)",
    ylabel = L"\xi_{(P)}", xlabel = L"S_{(P)}")
savefig("../../Figs/Q3/Q3-S-XI-Avg.pdf")


text:
# SPDX-License-Identifier: X11
# 2020-08-31
# Tsundoku (300pt)

function main()
  N, M, K = parse.(Int, split(readline(), " "))
  A = parse.(Int, split(readline(), " "))
  B = parse.(Int, split(readline(), " "))
  As = runsum(A)
  Bs = runsum(B)

  j = M
  ans = 0
  for i ∈ 0:N
    As[i + 1] > K && break

    while j > 0 && As[i + 1] + Bs[j + 1] > K
      j -= 1
    end

    ans = max(ans, i + j)
  end

  println(ans)
end

function runsum(A)
  As = Vector{Int}(undef, 1 + lastindex(A))
  As[1] = 0
  for i ∈ 1:lastindex(A)
    As[i + 1] = As[i] + A[i]
  end
  return As
end

main()

max_stars_repo_name: jthvai/competitive-programming

max_issues_repo_name: jthvai/competitive-programming

max_forks_repo_name: jthvai/competitive-programming


text:
# QuEST.jl/src/QuEST.jl
#
# Authors:
#  - Dirk Oliver Theis, Ketita Labs
#
# Copyright (c) 2020 Ketita Labs oü, Tartu, Estonia
#
# MIT License
#
# Copyright (c) 2020 Ketita Labs oü, Tartu, Estonia
#
# Permission is hereby granted, free of charge, to any person
# obtaining a copy of this software and associated documentation files
# (the "Software"), to deal in the Software without restriction,
# including without limitation the rights to use, copy, modify, merge,
# publish, distribute, sublicense, and/or sell copies of the Software,
# and to permit persons to whom the Software is furnished to do so,
# subject to the following conditions:
#
# The above copyright notice and this permission notice shall be
# included in all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
# BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
module QuEST

include("_QuEST_Internals.jl")
include("_C_QuEST.jl")

import _C_QuEST: numQubits, isDensityMatrix # ...
export numQubits, isDensityMatrix

################################################################################

module QuESTbase32

"""
QuEST real number type with 32 bit precision.
"""
const Qreal = Float32

"""
Handle for quantum register (i.e., pointer to `struct qureg` from QuEST.h) with 32-bit precision
"""
const QuReg = _QuEST_Internals.QuReg32

function __init__()
    if _QuEST_Internals.EXPERT_USER
        dlopen("libQuEST32_expert",RTLD_LAZY|RTLD_DEEPBIND|RTLD_GLOBAL)
    else
        dlopen("libQuEST32_novice",RTLD_LAZY|RTLD_DEEPBIND|RTLD_GLOBAL)
    end
    dlopen(joinpath(@__DIR__,"C","lib_C_QuEST32.so"),
           RTLD_LAZY|RTLD_DEEPBIND|RTLD_GLOBAL)

    if CHECK-THE-SIZE of qreal, see if it's really correct

    if CHECK-THE-SIZE of qureg, see if it's really correct

end

include("QuESTbase")

end

################################################################################

module QuESTbase64

"""
QuEST real number type with 64 bit precision.
"""
const qreal = Float64

"""
Handle for quantum register (i.e., pointer to `struct qureg` from QuEST.h) with 64-bit precision
"""
const QuReg = _QuEST_Internals.QuReg64

function __init__()
    if _QuEST_Internals.EXPERT_USER
        dlopen("libQuEST64_expert",RTLD_LAZY|RTLD_DEEPBIND|RTLD_GLOBAL)
    else
        dlopen("libQuEST64_novice",RTLD_LAZY|RTLD_DEEPBIND|RTLD_GLOBAL)
    end
    dlopen(joinpath(@__DIR__,"C","lib_C_QuEST64.so"),
           RTLD_LAZY|RTLD_DEEPBIND|RTLD_GLOBAL)

    if CHECK-THE-SIZE of qreal, see if it's really correct

    if CHECK-THE-SIZE of qureg, see if it's really correct

end

include("QuESTbase")

end

## Misc #-----------------------------------------------------------------------
#
# QuEST Env
#
#-------------------------------------------------------------------------------

function seedQuESTDefault() ::Nothing
    ccall(:seedQuESTDefault, Cvoid, () )
    return nothing
end

function seedQuEST(seedarray ::Vector{Int}) ::Nothing
    @assert  ! isempty(seedarray)

    ccall(:seedQuESTDefault, Cvoid, (Ptr{Clong},Cint),
          seedarray, length(seedarray))
    nothing;
end

function startRecordingQASM(qureg ::Qureg) ::Nothing
    ccall(:startRecordingQASM, Cvoid, (Qureg,),
          qureg)
    nothing;
end
function stopRecordingQASM(qureg ::Qureg) ::Nothing
    ccall(:stopRecordingQASM, Cvoid, (Qureg,),
          qureg)
    nothing;
end
function clearRecordedQASM(qureg ::Qureg) ::Nothing
    ccall(:clearRecordedQASM, Cvoid, (Qureg,),
          qureg)
    nothing;
end
function printRecordedQASM(qureg ::Qureg) ::Nothing
    ccall(:printRecordedQASM, Cvoid, (Qureg,),
          qureg)
    nothing;
end
function writeRecordedQASMToFile(qureg    ::Qureg,
                                 filename ::String) ::Nothing
    @assert begin
        ios = open(filename, "w")
        close(ios) === nothing
    end

    ccall(:writeRecordedQASMToFile, Cvoid, (Qureg, Cstring),
          qureg, filename)
    nothing;
end

## Env #------------------------------------------------------------------------
#
# QuEST Env
#
#    - create
#    - destroy
#
#    - seed
#-------------------------------------------------------------------------------

function createQuESTEnv() :: QuESTEnv
    return ccall(:createQuESTEnv, QuESTEnv, () )
end

function destroyQuESTEnv(env ::QuESTEnv) ::Nothing
    ccall(:destroyQuESTEnv, Cvoid, (QuESTEnv,),
          env)
    return nothing
end

function reportQuESTEnv(env ::QuESTEnv) ::Nothing
    ccall(:reportQuESTEnv, Cvoid, (QuESTEnv,),
          env)
    return nothing
end

function getEnvironmentString(env ::QuESTEnv, qureg ::Qureg) ::String
    cstr = Vector{Cchar}(undef,232)
    ccall(:getEnvironmentString, Cvoid, (QuESTEnv, Qureg, Ptr{Cchar}),
          env, qureg, cstr)
    return unsafe_string(pointer(cstr))
end

function copyStateToGPU(qureg ::Qureg) ::Nothing
    ccall(:copyStateToGPU, Cvoid, (Qureg,),
          qureg)
    return nothing
end

function copyStateFromGPU(qureg ::Qureg) ::Nothing
    ccall(:copyStateFromGPU, Cvoid, (Qureg,),
          qureg)
    return nothing
end

# function syncQuESTEnv(env ::QuESTEnv) ::Nothing
#     ccall(:syncQuESTEnv, Cvoid, (QuESTEnv,),
#           env)
#     return nothing
# end

## Qureg #----------------------------------------------------------------------
#
# Qureg
#
#    1. Create & destroy
#    2. Complex N-matrix
#    3. Report & get
#    4. Init
#    5. Query
#    6. Measure Copenhagen
#
#-------------------------------------------------------------------------------

#
#  1. Create & Destroy
#

function createQureg(numQubits ::Int, env ::QuESTEnv) ::Qureg
    @assert 1 ≤ numQubits ≤ 50

    return ccall(:createQureg, Qureg, (Cint,QuESTEnv),
                 numQubits, env)
end

function createDensityQureg(numQubits ::Int, env ::QuESTEnv) ::Qureg
    @assert 1 ≤ numQubits ≤ 50

    return ccall(:createDensityQureg, Qureg, (Cint,QuESTEnv),
                 numQubits, env)
end

function createCloneQureg(qureg ::Qureg, env ::QuESTEnv) ::Qureg
    return ccall(:createCloneQureg, Qureg, (Qureg,QuESTEnv),
                 qureg, env)
end

function destroyQureg(qureg ::Qureg, env ::QuESTEnv) ::Nothing
    ccall(:destroyQureg, Cvoid, (Qureg,QuESTEnv),
          qureg, env)
    return nothing
end

#
# 2. Complex N-Matrix
#

function createComplexMatrixN(numQubits ::Int) ::ComplexMatrixN;
    @assert 1 ≤ numQubits ≤ 50

    return ccall(:createComplexMatrixN, ComplexMatrixN, (Cint,),
                 numQubits)
end

@doc raw"
Function `make_QuEST_matrix(M ::Matrix{Qreal}) ::ComplexMatrixN`

Convenience function for creating (`createComplexMatrixN()`) and filling a QuEST matrix data structure (`struct ComplexMatrixN`).

### Input
* `M` must be a 2ⁿ×2ⁿ matrix for 1 ≤ 𝑛 ≤ 50 (haha) [`@assert`]

### Return value
* An 𝑛 qubit matrix.
"
function make_QuEST_matrix(M ::Matrix{Qreal}) ::ComplexMatrixN
    (R,C) = size(M)
    @assert R==C
    @assert R ≥ 2
    numQubits = Int(round( log2(R) ))
    @assert 2^numQubits == R
    @assert 1 ≤ numQubits ≤ 50

    MQ = createComplexMatrixN(numQubits)
    for c = 1:R                       # Julia is column major
        for r = 1:R
            re_MQ_r = unsafe_load(M.real,r)
            im_MQ_r = unsafe_load(M.real,r)
            unsafe_store!(re_MQ_r, real(M[r,c]), c)
            unsafe_store!(im_MQ_r, imag(M[r,c]), c)
        end
    end
    return MQ
end

function destroyComplexMatrixN(M ::ComplexMatrixN) ::Nothing
    ccall(:destroyComplexMatrixN, Cvoid, (ComplexMatrixN,),
          M)
    nothing
end

@doc raw"
Function

     fill_ComplexMatrix!(M ::ComplexMatrixN, F ::Function) ::Nothing

Fills the entries of `M`: ∀(k,ℓ): 𝑀[k,ℓ] = 𝐹(k,ℓ).
(Indices are in 1, ..., 2ⁿ.)

### Input
* `M` must have been ''created'' (`createComplexMatrixN()`)

### Output
* Entries of `M` are overwritten.
"
function fill_ComplexMatrix!(M ::ComplexMatrixN, M_ ::Function) ::Nothing
    N = 2^M.numQubits
    for k = 1:N
        re_M_k = unsafe_load(M.real,k)
        im_M_k = unsafe_load(M.real,k)
        for ℓ = 1:N
            unsafe_store!(re_M_k, real(M_(k,ℓ)), ℓ)
            unsafe_store!(im_M_k, imag(M_(k,ℓ)), ℓ)
        end
    end
    nothing
end

#
# 3. Report & get
#

function reportState(qureg ::Qureg) ::Nothing
    ccall(:reportState, Cvoid, (Qureg,),
          qureg)
    nothing;
end

function reportStateToScreen(qureg ::Qureg, env ::QuESTEnv, reportRank ::Int) ::Nothing
    ccall(:reportStateToScreen, Cvoid, (Qureg, QuESTEnv, Cint),
          qureg, env, reportRank)
    return nothing
end

function reportQuregParams(qureg ::Qureg) ::Nothing
    ccall(:reportQuregParams, Cvoid, (Qureg,),
          qureg)
    return nothing
end

function getNumQubits(qureg ::Qureg) ::Int
    return ccall(:getNumQubits, Cint, (Qureg,),
                 qureg)
end

function getNumAmps(qureg ::Qureg) ::Int
    return ccall(:getNumAmps, Clonglong, (Qureg,),
                 qureg)
end

## WIP: CONTINUE HERE with adding `@assert`s

#
# 4. Init state
#

function initBlankState(qureg ::Qureg) ::Nothing
    ccall(:initBlankState, Cvoid, (Qureg,),
          qureg)
    return nothing
end

function initZeroState(qureg ::Qureg) ::Nothing
    ccall(:initZeroState, Cvoid, (Qureg,),
          qureg)
    return nothing
end

function initPlusState(qureg ::Qureg) ::Nothing
    ccall(:initPlusState, Cvoid, (Qureg,),
          qureg)
    return nothing
end

function initClassicalState(qureg ::Qureg, stateInd ::Int64) ::Nothing
    ccall(:initClassicalState, Cvoid, (Qureg,Clonglong),
          qureg,stateInd)
    return nothing
end

function initPureState(qureg ::Qureg, pure ::Qureg) ::Nothing
    ccall(:initPureState, Cvoid, (Qureg,Qureg),
          qureg,pure)
    return nothing
end

function initDebugState(qureg ::Qureg) ::Nothing
    ccall(:initDebugState, Cvoid, (Qureg,),
          qureg)
    return nothing
end

function initStateFromAmps(qureg ::Qureg,
                           reals ::Vector{Qreal},
                           imags ::Vector{Qreal}  ) ::Nothing
    ccall(:initStateFromAmps, Cvoid, (Qureg, Ptr{Qreal}, Ptr{Qreal}),
          qureg, reals, imags)
    return nothing
end

function setAmps(qureg    ::Qureg,
                 startIdx ::Int64,
                 reals    ::Vector{Qreal},
                 imags    ::Vector{Qreal},
                 numAmps  ::Int64)          :: Nothing
    ccall(:setAmps, Cvoid, (Qureg, Clonglong, Ptr{Qreal}, Ptr{Qreal}, Clonglong),
          qureg,  startIdx,  reals, imags, numAmps)
    return nothing
end

function cloneQureg(targetQureg ::Qureg, copyQureg ::Qureg) ::Nothing
    ccall(:cloneQureg, Cvoid, (Qureg, Qureg),
          targetQureg, copyQureg)
    return nothing
end

function setWeightedQureg(fac1   ::Complex{Float64},   qureg1 ::Qureg,
                          fac2   ::Complex{Float64},   qureg2 ::Qureg,
                          facOut ::Complex{Float64},   out    ::Qureg) ::Nothing
    ccall(:setWeightedQureg,
          Cvoid,
          (QuEST_h.Complex, Qureg,
           QuEST_h.Complex, Qureg,
           QuEST_h.Complex, Qureg),
          QuEST_h.Complex(Qreal(real(fac1)),Qreal(real(fac1))), qureg1,
          QuEST_h.Complex(Qreal(real(fac2)),Qreal(real(fac2))), qureg2,
          QuEST_h.Complex(Qreal(real(out)),Qreal(real(out))), out)
    nothing;
end

#
# 5. Query state
#

function getAmp(qureg ::Qureg,  idx ::Int) ::Complex{Qreal}
    α = ccall(:getAmp, QuEST_h.Complex, (Qureg, Clonglong),
               qureg, idx)
    return Complex{Qreal}(α.real,α.imag)
end

function getProbAmp(qureg ::Qureg, idx ::Int) :: Qreal
    p = ccall(:getProbAmp, Qreal, (Qureg, Clonglong),
               qureg, idx)
    return p
end

function getDensityAmp(qureg ::Qureg, row ::Int, col ::Int) ::Complex{Qreal}
    α = ccall(:getDensityAmp, QuEST_h.Complex, (Qureg, Clonglong, Clonglong),
              qureg, row, col)
    return Complex{Qreal}(α.real, α.imag)
end

function calcTotalProb(qureg ::Qureg) ::Float64
    one = ccall(:calcTotalProb, Qreal, (Qureg,),
                qureg)
    return one
end

function calcProbOfOutcome(qureg        ::Qureg,
                           measureQubit ::Int,
                           outcome      ::Int)   ::Float64
    p = ccall(:calcProbOfOutcome, Qreal, (Qureg, Cint, Cint),
              qureg, measureQubit, outcome)
    return p
end

function calcExpecPauliProd(qureg        ::Qureg,
                            targetQubits ::Vector{Int32},
                            pauliCodes   ::Vector{Int32},
                            workspace    ::Qureg)          ::Float64
    @assert length(targetQubits) == length(pauliCodes)
    @assert all( σ -> 0 ≤ σ ≤ 3,   pauliCodes )

    expval =
           ccall(:calcExpecPauliProd, Qreal, (Qureg, Ptr{Cint}, Ptr{Cint}, Cint, Qureg),
                 qureg, targetQubits, pauliCodes, length(targetQubits),  workspace)

    return expval
end

function calcExpecPauliSum(qureg         ::Qureg,
                           allPauliCodes ::Vector{Int32},
                           termCoeffs    ::Vector{Qreal},
                           workspace     ::Qureg)          ::Float64

    @assert length(allPauliCodes) ==  length(termCoeffs) * getNumQubits(qureg)
    @assert all( σ -> 0 ≤ σ ≤ 3,  allPauliCodes )

    ex = ccall(:calcExpecPauliSum, Qreal, (Qureg, Ptr{Cint}, Ptr{Qreal}, Cint, Qureg),
               qureg, allPauliCodes, termCoeffs, length(termCoeffs),  workspace)

    return ex
end

function calcInnerProduct(bra ::Qureg, ket ::Qureg) ::Complex{Float64}
    w = ccall(:calcInnerProduct, QuEST_h.Complex, (Qureg, Qureg),
          bra, ket)
    return Complex{Float64}(w.real,w.imag)
end

function calcDensityInnerProduct(ϱ1 ::Qureg, ϱ2 ::Qureg) ::Float64
    w = ccall(:calcDensityInnerProduct, Qreal, (Qureg, Qureg),
              ϱ1, ϱ2)
    return Float64(w)
end

function calcPurity(qureg ::Qureg) ::Float64
    pu = ccall(:calcPurity, Qreal, (Qureg,),
               qureg)
    return pu
end

function calcFidelity(qureg ::Qureg, pureState ::Qureg) ::Float64
    fi = ccall(:calcFidelity, Qreal, (Qureg, Qureg),
               qureg, pureState)
    return fi
end

function calcHilbertSchmidtDistance(a ::Qureg, b ::Qureg) ::Float64
    hsd = ccall(:calcHilbertSchmidtDistance, Qreal, (Qureg, Qureg),
                a, b)
    return hsd
end

#
# 6. Copenhagen measure
#

function  collapseToOutcome(qureg        ::Qureg,
                            measureQubit ::Int,
                            outcome      ::Int)   ::Float64
    p = ccall(:collapseToOutcome, Qreal, (Qureg, Cint, Cint),
              qureg, measureQubit, outcome)
    return p
end

function measure(qureg ::Qureg, measureQubit ::Int) ::Int
    i = ccall(:measure, Cint, (Qureg, Cint),
              qureg, measureQubit)
    return i
end

function measureWithStats(qureg        ::Qureg,
                          measureQubit ::Int     ) ::Tuple{Int,Float64}
    p = Qreal(0.)
    i = ccall(:measureWithStats, Cint, (Qureg, Cint, Ref{Qreal}),
          qureg, measureQubit, p)
    return (Int(i),Float64(p))
end

## Apply operation #------------------------------------------------------------
#
# Apply quantum operation
#
#    1. Gates
#    2. Unitary matrices
#    3. Noise
#    4. Other ops
#
#-------------------------------------------------------------------------------

#
# 1. Gates
#

function phaseShift(qureg       ::Qureg,
                    targetQubit ::Int,
                    angle       ::Qreal) ::Nothing
    ccall(:phaseShift, Cvoid, (Qureg, Cint, Qreal),
           qureg, targetQubit, angle)
    return nothing
end

function controlledPhaseShift(qureg    ::Qureg,
                              idQubit1 ::Int,
                              idQubit2 ::Int,
                              angle    ::Qreal)   ::Nothing
    ccall(:controlledPhaseShift, Cvoid, (Qureg, Cint, Cint, Qreal),
          qureg, idQubit1, idQubit2, angle)
    return nothing
end

function multiControlledPhaseShift(qureg         ::Qureg,
                                   controlQubits ::Vector{Int32},
                                   angle         ::Qreal)         ::Nothing
    ccall(:multiControlledPhaseShift, Cvoid, (Qureg, Ptr{Cint}, Cint, Qreal),
          qureg, controlQubits, length(controlQubits),  angle)
    return nothing
end

function controlledPhaseFlip(qureg ::Qureg, idQubit1 ::Int, idQubit2 ::Int) ::Nothing
    ccall(:controlledPhaseFlip, Cvoid, (Qureg, Cint, Cint),
           qureg, idQubit1, idQubit2)
    return nothing
end

function multiControlledPhaseFlip(qureg         ::Qureg,
                                  controlQubits ::Vector{Int32}) ::Nothing
    ccall(:multiControlledPhaseFlip, Cvoid, (Qureg, Ptr{Cint}, Cint),
          qureg, controlQubits, length(controlQubits))
    return nothing
end

function sGate(qureg ::Qureg, targetQubit ::Int) ::Nothing
    ccall(:sGate, Cvoid, (Qureg, Cint),
           qureg,  targetQubit)
    return nothing
end

function tGate(qureg ::Qureg, targetQubit ::Int) ::Nothing
    ccall(:tGate, Cvoid, (Qureg, Cint),
          qureg, targetQubit)
    return nothing
end

function pauliX(qureg ::Qureg, targetQubit ::Int):: Nothing
    ccall(:pauliX, Cvoid, (Qureg, Cint),
          qureg,  targetQubit)
    return nothing
end
function pauliY(qureg ::Qureg, targetQubit ::Int):: Nothing
    ccall(:pauliY, Cvoid, (Qureg, Cint),
          qureg,  targetQubit)
    return nothing
end
function pauliZ(qureg ::Qureg, targetQubit ::Int):: Nothing
    ccall(:pauliZ, Cvoid, (Qureg, Cint),
          qureg,  targetQubit)
    return nothing
end

function hadamard(qureg ::Qureg,  targetQubit ::Int) ::Nothing
    ccall(:hadamard, Cvoid, (Qureg, Cint),
          qureg, targetQubit)
    return nothing
end

function controlledNot(qureg         ::Qureg,
                       controlQubit  ::Int,
                       targetQubit   ::Int)   ::Nothing
    ccall(:controlledNot, Cvoid, (Qureg, Cint, Cint),
           qureg, controlQubit, targetQubit)
    return nothing
end

function controlledPauliY(qureg         ::Qureg,
                       controlQubit  ::Int,
                       targetQubit   ::Int)   ::Nothing
    ccall(:controlledPauliY, Cvoid, (Qureg, Cint, Cint),
           qureg, controlQubit, targetQubit)
    return nothing
end

function controlledPauliZ(qureg         ::Qureg,
                       controlQubit  ::Int,
                       targetQubit   ::Int)   ::Nothing
    ccall(:controlledPauliZ, Cvoid, (Qureg, Cint, Cint),
           qureg, controlQubit, targetQubit)
    return nothing
end

function rotateX(qureg ::Qureg, rotQubit ::Int, angle ::Float64) ::Nothing
    ccall(:rotateX, Cvoid, (Qureg, Cint, Qreal),
          qureg, rotQubit, angle)
    return nothing
end
function rotateY(qureg ::Qureg, rotQubit ::Int, angle ::Float64) ::Nothing
    ccall(:rotateY, Cvoid, (Qureg, Cint, Qreal),
          qureg, rotQubit, angle)
    return nothing
end
function rotateZ(qureg ::Qureg, rotQubit ::Int, angle ::Float64) ::Nothing
    ccall(:rotateZ, Cvoid, (Qureg, Cint, Qreal),
          qureg, rotQubit, angle)
    return nothing
end

function rotateAroundAxis(qureg ::Qureg, rotQubit ::Int, angle ::Float64, axis ::Vector) ::Nothing
    ccall(:rotateAroundAxis, Cvoid, (Qureg, Cint, Qreal, Vector),
          qureg, rotQubit, angle, axis)
    return nothing
end

function controlledRotateX(qureg         ::Qureg,
                           controlQubit  ::Int,
                           targetQubit   ::Int,
                           angle         ::Qreal)  ::Nothing
    ccall(:controlledRotateX, Cvoid, (Qureg, Cint, Cint, Qreal),
          qureg, controlQubit, targetQubit, angle)
    return nothing
end

function controlledRotateY(qureg         ::Qureg,
                           controlQubit  ::Int,
                           targetQubit   ::Int,
                           angle         ::Qreal)  ::Nothing
    ccall(:controlledRotateY, Cvoid, (Qureg, Cint, Cint, Qreal),
          qureg, controlQubit, targetQubit, angle)
    return nothing
end

function controlledRotateZ(qureg         ::Qureg,
                           controlQubit  ::Int,
                           targetQubit   ::Int,
                           angle         ::Qreal)  ::Nothing
    ccall(:controlledRotateZ, Cvoid, (Qureg, Cint, Cint, Qreal),
          qureg, controlQubit, targetQubit, angle)
    return nothing
end

function multiRotateZ(qureg      ::Qureg,
                      qubits     ::Vector{Int32},
                      angle      ::Qreal)         ::Nothing
    ccall(:multiRotateZ, Cvoid, (Qureg, Ptr{Cint}, Cint, Qreal),
          qureg, qubits, length(qubits), angle)
    return nothing
end

function multiRotatePauli(qureg         ::Qureg,
                          targetQubits  ::Vector{Int32},
                          targetPaulis  ::Vector{Int32},
                          angle         ::Qreal)          ::Nothing
    @assert length(targetQubits) == length(targetPaulis)
    @assert all( σ -> 0 ≤ σ ≤ 3,   targetPaulis )

    ccall(:multiRotatePauli, Cvoid, (Qureg, Ptr{Cint}, Ptr{Cint}, Cint, Qreal),
          qureg, targetQubits, targetPaulis, length(targetPaulis), angle)

    return nothing
end

function controlledRotateAroundAxis(qureg        ::Qureg,
                                    controlQubit ::Int,
                                    targetQubit  ::Int,
                                    angle        ::Qreal,
                                    axis         ::Vector)  ::Nothing
    ccall(:controlledRotateAroundAxis, Cvoid, (Qureg, Cint, Cint, Qreal, Vector),
          qureg, controlQubit, targetQubit, angle, axis)
    return nothing
end

function swapGate(qureg ::Qureg, qubit1 ::Int, qubit2 ::Int) ::Nothing
    ccall(:swapGate, Cvoid, (Qureg, Cint, Cint),
          qureg, qubit1, qubit2)
    return nothing
end

function sqrtSwapGate(qureg ::Qureg, qubit1 ::Int, qubit2 ::Int) ::Nothing
    ccall(:sqrtSwapGate, Cvoid, (Qureg, Cint, Cint),
          qureg, qubit1, qubit2)
    return nothing
end

#
# 2. Apply unitary matrices
#

function compactUnitary(qureg       ::Qureg,
                        targetQubit ::Int,
                        α           ::Complex{Qreal},
                        β           ::Complex{Qreal})   ::Nothing
    alpha = QuEST_h.Complex(real(α),imag(α))
    beta  = QuEST_h.Complex(real(β),imag(β))

    ccall(:compactUnitary, Cvoid, (Qureg, Cint, QuEST_h.Complex, QuEST_h.Complex),
          qureg, targetQubit, alpha, beta)

    return nothing
end

function controlledCompactUnitary(qureg          ::Qureg,
                                  controlQubit   ::Int,
                                  targetQubit    ::Int,
                                  α              ::Complex{Qreal},
                                  β              ::Complex{Qreal} ) ::Nothing
    alpha = QuEST_h.Complex(real(α),imag(α))
    beta  = QuEST_h.Complex(real(β),imag(β))

    ccall(:controlledCompactUnitary, Cvoid, (Qureg, Int, Cint, Complex, Complex),
          qureg, controlQubit, targetQubit, alpha, beta)
    return nothing
    end

function _quest_mtx_2(U ::Matrix{Complex{Qreal}}) ::ComplexMatrix2
    @assert size(U) == (2,2)
    u = ComplexMatrix2(
        ( (real(U[1,1]), real(U[1,2])), (real(U[2,1]), real(U[2,2])) ),
        ( (imag(U[1,1]), imag(U[1,2])), (imag(U[2,1]), imag(U[2,2])) )  )
    return u
end

function unitary(qureg           ::Qureg,
                 targetQubit     ::Int,
                 U               ::Matrix{Complex{Qreal}}) ::Nothing

    u = _quest_mtx_2(U)

    ccall(:unitary, Cvoid, (Qureg, Cint, ComplexMatrix2),
          qureg, targetQubit,  u)
    return nothing
end

function controlledUnitary(qureg        ::Qureg,
                           controlQubit ::Int,
                           targetQubit  ::Int,
                           U            ::Matrix{Complex{Qreal}}) ::Nothing

    u = _quest_mtx_2(U)

    ccall(:controlledUnitary, Cvoid, (Qureg, Cint, Cint, ComplexMatrix2),
          qureg, controlQubit, targetQubit, u)
    return nothing
end

function multiControlledUnitary(qureg         ::Qureg,
                                controlQubits ::Vector{Int32},
                                targetQubit   ::Int,
                                U             ::Matrix{Complex{Qreal}}) ::Nothing
    u = _quest_mtx_2(U)

    ccall(:multiControlledUnitary, Cvoid, (Qureg, Ptr{Cint}, Cint, Cint, ComplexMatrix2),
          qureg, controlQubits, length(controlQubits),  targetQubit, u)
    return nothing
end

function multiStateControlledUnitary(qureg           ::Qureg,
                                     controlQubits   ::Vector{Int32},
                                     controlState    ::Vector{Int32},
                                     targetQubit     ::Int,
                                     U               ::Matrix{Complex{Qreal}}) ::Nothing

    @assert length(controlQubits) == length(controlState)
    u = _quest_mtx_2(U)

    ccall(:multiStateControlledUnitary,
          Cvoid,
          (Qureg, Ptr{Cint}, Ptr{Cint}, Cint, Cint, ComplexMatrix2),
          qureg, controlQubits, controlState, length(controlQubits), targetQubit, u)
    return nothing
end

function _quest_mtx_4(U ::Matrix{Complex{Qreal}}) ::ComplexMatrix4
    @assert size(U) = (4,4)
    u = ComplexMatrix4( ( ( real(U[1,1]), real(U[1,2]), real(U[1,3]), real(U[1,4]) ),
                          ( real(U[2,1]), real(U[2,2]), real(U[2,3]), real(U[2,4]) ),
                          ( real(U[3,1]), real(U[3,2]), real(U[3,3]), real(U[3,4]) ),
                          ( real(U[4,1]), real(U[4,2]), real(U[4,3]), real(U[4,4]) ) ),
                        ( ( imag(U[1,1]), imag(U[1,2]), imag(U[1,3]), imag(U[1,4]) ),
                          ( imag(U[2,1]), imag(U[2,2]), imag(U[2,3]), imag(U[2,4]) ),
                          ( imag(U[3,1]), imag(U[3,2]), imag(U[3,3]), imag(U[3,4]) ),
                          ( imag(U[4,1]), imag(U[4,2]), imag(U[4,3]), imag(U[4,4]) ) )  )
end

function twoQubitUnitary(qureg           ::Qureg,
                         targetQubit1    ::Int,
                         targetQubit2    ::Int,
                         U               ::Matrix{Complex{Qreal}}) ::Nothing
    u = _quest_mtx_4(U)
    ccall(:twoQubitUnitary, Cvoid, (Qureg, Cint, Cint, ComplexMatrix4),
          qureg, targetQubit1, targetQubit2, u)

    return nothing
end

function controlledTwoQubitUnitary(qureg         ::Qureg,
                                   controlQubit  ::Int,
                                   targetQubit1  ::Int,
                                   targetQubit2  ::Int,
                                   U             ::Matrix{Complex{Qreal}}) ::Nothing
    u = _quest_mtx_4(U)
    ccall(:controlledTwoQubitUnitary, Cvoid, (Qureg, Cint, Cint, Cint, ComplexMatrix4),
          qureg, controlQubit, targetQubit1, targetQubit2, u)
    return nothing
end

function multiControlledTwoQubitUnitary(qureg         ::Qureg,
                                        controlQubits ::Vector{Int32},
                                        targetQubit1 ::Int,
                                        targetQubit2 ::Int,
                                        U            ::Matrix{Complex{Qreal}}) ::Nothing
    u = _quest_mtx_4(U)
    ccall(:multiControlledTwoQubitUnitary, Cvoid, (Qureg, Ptr{Cint}, Cint, Cint, Cint, ComplexMatrix4),
          qureg, controlQubits, length(controlQubits),  targetQubit1, targetQubit2, u)
    return nothing
end

function multiQubitUnitary(qureg ::Qureg,
                           targs ::Vector{Int32},
                           u     ::ComplexMatrixN) ::Nothing
    ccall(:multiQubitUnitary, Cvoid, (Qureg, Ptr{Cint}, Cint, ComplexMatrixN),
          qureg, targs, length(targs), u)
    return nothing
end

function controlledMultiQubitUnitary(qureg   ::Qureg,
                                     ctrl    ::Int,
                                     targs   ::Vector{Int32},
                                     u       ::ComplexMatrixN) ::Nothing
    ccall(:controlledMultiQubitUnitary, Cvoid, (Qureg, Cint, Ptr{Cint}, Cint, ComplexMatrixN),
          qureg, ctrl, targs, length(targs), u)
    return nothing
end

function multiControlledMultiQubitUnitary(qureg  ::Qureg,
                                          ctrls  ::Vector{Int32},
                                          targs  ::Vector{Int32},
                                          u      ::ComplexMatrixN) ::Nothing
    ccall(:multiControlledMultiQubitUnitary,
          Cvoid,
          (Qureg, Ptr{Cint}, Cint, Ptr{Cint}, Cint, ComplexMatrixN),
          qureg, ctrls, length(ctrls), targs, length(targs), u)
    return nothing
end

#
# 3. Noise
#

function mixDephasing(qureg ::Qureg, targetQubit ::Int, prob ::Float64) ::Nothing
    ccall(:mixDephasing, Cvoid, (Qureg, Cint, Qreal),
          qureg, targetQubit, prob)
    return nothing
end

function mixDepolarising(qureg ::Qureg, targetQubit ::Int, prob ::Float64) ::Nothing
    ccall(:mixDepolarising, Cvoid, (Qureg, Cint, Qreal),
          qureg, targetQubit, prob)
    return nothing
end

function mixDamping(qureg ::Qureg, targetQubit ::Int, prob ::Float64) ::Nothing
    ccall(:mixDamping, Cvoid, (Qureg, Cint, Qreal),
          qureg, targetQubit, prob)
    return nothing
end

function mixTwoQubitDephasing(qureg  ::Qureg,
                              qubit1 ::Int,
                              qubit2 ::Int,
                              prob   ::Float64)   ::Nothing
    ccall(:mixTwoQubitDephasing, Cvoid, (Qureg, Cint, Cint, Qreal),
          qureg, qubit1, qubit2, prob)
    return nothing
end

function mixTwoQubitDepolarising(qureg  ::Qureg,
                              qubit1 ::Int,
                              qubit2 ::Int,
                              prob   ::Float64)   ::Nothing
    ccall(:mixTwoQubitDepolarising, Cvoid, (Qureg, Cint, Cint, Qreal),
          qureg, qubit1, qubit2, prob)
    return nothing
end

function mixPauli(qureg       ::Qureg,
                  targetQubit ::Int,
                  probX       ::Float64,
                  probY       ::Float64,
                  probZ       ::Float64)   ::Nothing
    ccall(:mixPauli, Cvoid, (Qureg, Cint, Qreal, Qreal, Qreal),
          qureg, targetQubit, probX, probY, probZ)
    return nothing
end

function mixDensityMatrix(combineQureg ::Qureg,
                          prob         ::Qreal,
                          otherQureg   ::Qureg)  ::Nothing
    ccall(:mixDensityMatrix, Cvoid, (Qureg, Qreal, Qureg),
          combineQureg, prob, otherQureg)
    return nothing
end

function mixKrausMap(qureg   ::Qureg,
                     target  ::Int,
                     ops     ::Vector{Matrix{Complex{Qreal}}}) ::Nothing

    @assert 1 ≤ length(ops) ≤ 4
    qops = [ _quest_mtx_2(op) for op in ops ]
    ccall(:mixKrausMap, Cvoid, (Qureg, Cint, Ptr{ComplexMatrix2}, Cint),
          qureg, target, qops, length(ops))
    nothing;
end

function mixTwoQubitKrausMap(qureg    ::Qureg,
                             target1  ::Int,
                             target2  ::Int,
                             ops      ::Vector{Matrix{Complex{Qreal}}}) ::Nothing
    @assert 1 ≤ length(ops) ≤ 16
    qops = [ _quest_mtx_4(op) for op in ops ]
    ccall(:mixTwoQubitKrausMap,
          Cvoid,
          (Qureg, Cint, Cint, Ptr{ComplexMatrix4}, Cint),
          qureg, target1, target2, qops, length(ops))
    nothing;
end

function mixMultiQubitKrausMap(qureg    ::Qureg,
                               targets  ::Vector{Int32},
                               ops      ::Vector{ComplexMatrixN}) ::Nothing
    N = length(targets)
    @assert 1 ≤ length(ops) ≤ (2N)^2
    ccall(:mixMultiQubitKrausMap, Cvoid,
          (Qureg, Ptr{Cint}, Cint, Ptr{ComplexMatrixN}, Cint),
           qureg, targets, length(targets), ops, length(ops))
    nothing;
end

#
# 4. Other ops
#

function applyPauliSum(inQureg        ::Qureg,
                       allPauliCodes  ::Vector{Int32},
                       termCoeffs     ::Vector{Qreal},
                       numSumTerms    ::Int,
                       outQureg       ::Qureg)          ::Nothing
    @assert length(termCoeffs) == numSumTerms * getNumQubits(inQureg)

    ccall(:applyPauliSum, Cvoid, (Qureg, Ptr{Cint}, Ptr{Qreal}, Cint, Qureg),
          inQureg, allPauliCodes, termCoeffs, numSumTerms, outQureg)

    nothing;
end

## Julia Module Init #----------------------------------------------------------
#
# Init
#
#-------------------------------------------------------------------------------
using Libdl: dlopen, RTLD_LAZY, RTLD_DEEPBIND, RTLD_GLOBAL

function __init__()
    # nothing to do...?
end

## Includes #-------------------------------------------------------------------
#
# Include Sub-Modules
#
#-------------------------------------------------------------------------------

include("_Tomog.jl")

end # module
# EOF


text:
function map_lookup_elem(map::RTMap{Name,MT,K,V,ME,F}, key::K) where {Name,MT,K,V,ME,F}
    keyref = ZeroInitRef(K, key)
    GC.@preserve keyref begin
        keyref_ptr = Base.unsafe_convert(Ptr{K}, keyref)
        _map_lookup_elem(map, keyref_ptr)
    end
end
@inline function map_update_elem(map::RTMap{Name,MT,K,V,ME,F}, key::K, value::V, flags::UInt64) where {Name,MT,K,V,ME,F}
    keyref = ZeroInitRef(K, key)
    valref = ZeroInitRef(V, value)
    GC.@preserve keyref valref begin
        keyref_ptr = Base.unsafe_convert(Ptr{K}, keyref)
        valref_ptr = Base.unsafe_convert(Ptr{V}, valref)
        _map_update_elem(map, keyref_ptr, valref_ptr, flags)
    end
end
function map_delete_elem(map::RTMap{Name,MT,K,V,ME,F}, key::K) where {Name,MT,K,V,ME,F}
    keyref = ZeroInitRef(K, key)
    GC.@preserve keyref begin
        keyref_ptr = Base.unsafe_convert(Ptr{K}, keyref)
        _map_delete_elem(map, Base.unsafe_convert(Ptr{K}, keyref))
    end
end

# TODO: Use bpfcall
@generated function _map_lookup_elem(map::RTMap{Name,MT,K,V,ME,F}, key::Ptr{K}) where {Name,MT,K,V,ME,F}
    Context() do ctx
        T_keyp = LLVM.PointerType(convert(LLVMType, K; ctx))
        T_valp = LLVM.PointerType(convert(LLVMType, V; ctx))
        T_jlptr = convert(LLVMType, key; ctx)

        llvm_f, _ = create_function(T_jlptr, [T_jlptr])
        mod = LLVM.parent(llvm_f)

        map_gv = _genmap!(mod, map, ctx)
        T_map_gv = llvmtype(map_gv)

        Builder(ctx) do builder
            entry = BasicBlock(llvm_f, "entry"; ctx)
            position!(builder, entry)
            ft = LLVM.FunctionType(T_valp, [T_map_gv, T_keyp])
            ftp = LLVM.PointerType(ft)
            elem = inttoptr!(builder, parameters(llvm_f)[1], T_keyp)
            f = inttoptr!(builder, ConstantInt(Int64(1); ctx), ftp)
            value = call!(builder, f, [map_gv, elem])
            value = ptrtoint!(builder, value, T_jlptr)
            ret!(builder, value)
        end
        call_function(llvm_f, Ptr{V}, Tuple{Ptr{K}}, :key)
    end
end
@generated function _map_update_elem(map::RTMap{Name,MT,K,V,ME,F}, key::Ptr{K}, val::Ptr{V}, flags::UInt64) where {Name,MT,K,V,ME,F}
    Context() do ctx
        T_cint = convert(LLVMType, Cint; ctx)
        T_keyp = LLVM.PointerType(convert(LLVMType, K; ctx))
        T_valp = LLVM.PointerType(convert(LLVMType, V; ctx))
        T_flags = convert(LLVMType, flags; ctx)
        T_jlptr = convert(LLVMType, key; ctx)

        llvm_f, _ = create_function(T_cint, [T_jlptr, T_jlptr, T_flags])
        mod = LLVM.parent(llvm_f)

        map_gv = _genmap!(mod, map, ctx)
        T_map_gv = llvmtype(map_gv)

        Builder(ctx) do builder
            entry = BasicBlock(llvm_f, "entry"; ctx)
            position!(builder, entry)
            ft = LLVM.FunctionType(T_cint, [T_map_gv, T_keyp, T_valp, T_flags])
            ftp = LLVM.PointerType(ft)
            key_elem = inttoptr!(builder, parameters(llvm_f)[1], T_keyp)
            val_elem = inttoptr!(builder, parameters(llvm_f)[2], T_valp)
            f = inttoptr!(builder, ConstantInt(Int64(2); ctx), ftp)
            value = call!(builder, f, [map_gv, key_elem, val_elem, parameters(llvm_f)[3]])
            ret!(builder, value)
        end
        call_function(llvm_f, Cint, Tuple{Ptr{K},Ptr{V},UInt64}, :key, :val, :flags)
    end
end
@generated function _map_delete_elem(map::RTMap{Name,MT,K,V,ME,F}, key::Ptr{K}) where {Name,MT,K,V,ME,F}
    Context() do ctx
        T_cint = convert(LLVMType, Cint; ctx)
        T_keyp = LLVM.PointerType(convert(LLVMType, K; ctx))
        T_jlptr = convert(LLVMType, key; ctx)

        llvm_f, _ = create_function(T_cint, [T_jlptr])
        mod = LLVM.parent(llvm_f)

        map_gv = _genmap!(mod, map, ctx)
        T_map_gv = llvmtype(map_gv)

        Builder(ctx) do builder
            entry = BasicBlock(llvm_f, "entry"; ctx)
            position!(builder, entry)
            ft = LLVM.FunctionType(T_cint, [T_map_gv, T_keyp])
            ftp = LLVM.PointerType(ft)
            elem = inttoptr!(builder, parameters(llvm_f)[1], T_keyp)
            f = inttoptr!(builder, ConstantInt(Int64(3); ctx), ftp)
            value = call!(builder, f, [map_gv, elem])
            ret!(builder, value)
        end
        call_function(llvm_f, Cint, Tuple{Ptr{K}}, :key)
    end
end

# BPF-safe convert alternative
bpfconvert(T, ::Nothing) = nothing
bpfconvert(::Type{T}, v::T) where {T} = v
function bpfconvert(::Type{T}, v::V) where {T<:Unsigned,V<:Unsigned}
    if v > typemax(T)
        nothing
    elseif typemax(V) > typemax(T)
        unsafe_trunc(T, v)
    else
        convert(T, v)
    end::T
end
function bpfconvert(::Type{T}, v::V) where {T,V}
    if (v > typemax(T)) || (v < typemin(T))
        nothing
    elseif (typemax(V) > typemax(T)) || (typemin(V) < typemin(T))
        unsafe_trunc(T, v)
    else
        convert(T, v)
    end::Union{T,Nothing}
end

@inline Base.getindex(map::RTMap{Name,MT,K,V,ME,F}, idx) where {Name,MT,K,V,ME,F} =
    getindex(map, bpfconvert(K, idx))
Base.getindex(map::RTMap, ::Nothing) = nothing
@inline function Base.getindex(map::AbstractHashMap{Name,MT,K,V,ME,F}, idx::K) where {Name,MT,K,V,ME,F}
    ptr = map_lookup_elem(map, idx)
    if reinterpret(UInt64, ptr) > 0
        return unsafe_load(ptr)
    else
        return nothing
    end
end
@inline function Base.getindex(map::AbstractArrayMap{Name,MT,K,V,ME,F}, idx::K) where {Name,MT,K,V,ME,F}
    if idx > 0
        ptr = map_lookup_elem(map, idx-K(1))
        if reinterpret(UInt64, ptr) > 0
            return unsafe_load(ptr)
        else
            return nothing
        end
    else
        return nothing
    end
end

@inline Base.setindex!(map::RTMap{Name,MT,K,V,ME,F}, value, idx) where {Name,MT,K,V,ME,F} =
    setindex!(map, bpfconvert(V, value), bpfconvert(K, idx))
@inline Base.setindex!(map::RTMap, ::Nothing, idx) = nothing
@inline Base.setindex!(map::RTMap, value, ::Nothing) = nothing
@inline Base.setindex!(map::RTMap, ::Nothing, ::Nothing) = nothing
@inline function Base.setindex!(map::AbstractHashMap{Name,MT,K,V,ME,F}, value::V, idx::K) where {Name,MT,K,V,ME,F}
    map_update_elem(map, idx, value, UInt64(0))
    value
end
@inline function Base.setindex!(map::AbstractArrayMap{Name,MT,K,V,ME,F}, value::V, idx::K) where {Name,MT,K,V,ME,F}
    if idx > 0
        map_update_elem(map, idx-K(1), value, UInt64(0))
    end
    value
end

@inline Base.delete!(map::RTMap{Name,MT,K,V,ME,F}, idx) where {Name,MT,K,V,ME,F} =
    delete!(map, bpfconvert(K, idx))
Base.delete!(map::RTMap, ::Nothing) = nothing
@inline function Base.delete!(map::AbstractHashMap{Name,MT,K,V,ME,F}, idx::K) where {Name,MT,K,V,ME,F}
    map_delete_elem(map, idx)
    map
end
@inline function Base.delete!(map::AbstractArrayMap{Name,MT,K,V,ME,F}, idx::K) where {Name,MT,K,V,ME,F}
    if idx > 0
        map_delete_elem(map, idx-K(1))
    end
    map
end

@inline Base.haskey(map::AbstractHashMap{Name,MT,K,V,ME,F}, idx) where {Name,MT,K,V,ME,F} =
    map[bpfconvert(K, idx)] !== nothing
@inline Base.haskey(map::RTMap, ::Nothing) = false
@inline function Base.haskey(map::AbstractArrayMap{Name,MT,K,V,ME,F}, idx) where {Name,MT,K,V,ME,F}
    if idx > 0
        map[bpfconvert(K, idx)-K(1)] !== nothing
    else
        false
    end
end

@inline function Base.get(map::RTMap{Name,MT,K,V,ME,F}, k::K, v::V) where {Name,MT,K,V,ME,F}
    map_v = map[k]
    if map_v !== nothing
        return map_v
    else
        return v
    end
end
@inline Base.get(map::RTMap{Name,MT,K,V,ME,F}, k, v) where {Name,MT,K,V,ME,F} =
    get(map, bpfconvert(K, k), bpfconvert(V, v))

## Perf

#= TODO
export perf_event_output

function perf_event_output(pt_ctx::Ptr{API.pt_regs}, map::RTMap, flags::Union{Int64,UInt64}, data::D, sz::Union{Int64,UInt64}) where D
    data_ref = Ref{D}(data)
    GC.@preserve data_ref begin
        data_ptr = Base.unsafe_convert(Ptr{D}, data_ref)
        _perf_event_output(pt_ctx, map, reinterpret(UInt64, flags), data_ptr, reinterpret(UInt64, sz))
    end
end
@generated function _perf_event_output(pt_ctx::Ptr{API.pt_regs}, map::RTMap{Name,MT,K,V,ME,F}, flags::UInt64, data::Ptr{D}, sz::UInt64) where {Name,MT,K,V,ME,F,D}
    Context() do ctx
        T_cint = convert(LLVMType, Cint, ctx)
        T_ptctx = convert(LLVMType, pt_ctx, ctx)
        T_datap = LLVM.PointerType(convert(LLVMType, data, ctx))
        T_u64 = LLVM.Int64Type(ctx)
        T_jlptr = convert(LLVMType, Ptr{Cvoid}, ctx)

        llvm_f, _ = create_function(T_cint, [T_ptctx, T_u64, T_jlptr, T_u64])
        mod = LLVM.parent(llvm_f)

        map_gv = _genmap!(mod, map, ctx)
        T_map_gv = llvmtype(map_gv)

        Builder(ctx) do builder
            entry = BasicBlock(llvm_f, "entry", ctx)
            position!(builder, entry)
            ft = LLVM.FunctionType(T_cint, [T_ptctx, T_map_gv, T_u64, T_datap, T_u64])
            ftp = LLVM.PointerType(ft)
            data_elem = inttoptr!(builder, parameters(llvm_f)[3], T_datap)
            f = inttoptr!(builder, ConstantInt(Int64(25), ctx), ftp)
            value = call!(builder, f, [parameters(llvm_f)[1], map_gv, parameters(llvm_f)[2], data_elem, parameters(llvm_f)[4]])
            ret!(builder, value)
        end
        call_function(llvm_f, Cint, Tuple{Ptr{API.pt_regs},UInt64,Ptr{D},UInt64}, :((pt_ctx,flags,data,sz)))
    end
end
=#