text:
using Random
using StatsBase
using Distributions

struct NaiveBayes
    μ::Array{Float64,2}
    σ::Array{Float64,2}
    priors::Array{Float64,1}
    classes::Array{Int64,1}
end

function fit_params(X, y)
    n_f, n_obs = size(X)
    classes = unique(y)
    nclasses = length(classes)
    μ = zeros(Float64, n_f, nclasses)
    σ = zeros(Float64, n_f, nclasses)
    obs_per_class = (y .== cl for cl in classes)
    priors = mean.(obs_per_class) #i = [0.5,0.5] #
    for (i,o) in enumerate(obs_per_class)
        μ[:, i:i] .= mean(X[:,o],dims=2)
        σ[:, i:i] .= std(X[:,o],dims=2)
    end
    return μ, σ, priors, classes
end

import StatsBase.predict

function predict(d::NaiveBayes, X::Matrix)
    nclasses = length(d.priors)
    n_f, n_obs = size(X)
    predicted_class = zeros(Int, n_obs)
    ll = zeros(Float64, nclasses)
    for obs in 1:n_obs
        for cl in 1:nclasses
            ll[cl] = log(d.priors[cl])
            ll[cl] += sum(logpdf.(Normal.(d.μ[:,cl], d.σ[:,cl]), X[:, obs].+1e-8))
        end
        predicted_class[obs] = d.classes[findmax(ll)[2]]
    end
    return predicted_class
end

import StatsBase.fit
function fit(model::Type{NaiveBayes}, X::Array{Float64,2}, y::Array{Int64,1})
    μ, σ, priors,classes = fit_params(X, y)
    NaiveBayes(μ, σ, priors, classes)
end

function test()
	Random.seed!(0)
	X = rand(5, 100) #features, observatins
	y = ones(Int, 100)
	y[1:50] .= 2
	X[:, 1:50] .*= 2
	X[:, 51:end] ./= 2

	model = fit(NaiveBayes, X,y)
	prediction = predict(model, X)
	@assert prediction == y
end


text:
using IntrusiveCollections
using Test

@testset "IntrusiveCollections.jl" begin
    include("list_tests.jl")
end

ext: jl


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "cspice"
version = v"66"

# Collection of sources required to build cspice
sources = [
    "http://naif.jpl.nasa.gov/pub/naif/toolkit//C/PC_Linux_GCC_64bit/packages/cspice.tar.Z" =>
    "93cd4fbce5818f8b7fecf3914c5756b8d41fd5bdaaeac1f4037b5a5410bc4768",
    "./extra"
]

# These are the platforms built inside the wizard
platforms = supported_platforms()

script = raw"""
cd $WORKSPACE/srcdir
cp -r cspice/src/cspice/ .
patch -d cspice < dskx02.patch
patch -d cspice < subpnt.patch

# Create generated cmake files outside of source tree
mkdir build
cd build

cmake -DCMAKE_INSTALL_PREFIX=${prefix} -DCMAKE_TOOLCHAIN_FILE=/opt/$target/$target.toolchain ..
make -j${nproc} VERBOSE=1
make install VERBOSE=1
"""

products = prefix -> [
    LibraryProduct(prefix, "libcspice", :libcspice)
]

# Dependencies that must be installed before this package can be built
dependencies = []

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
"""
    ismasked(data, I...)

Check if a cell is masked, using the `mask` array.

Used used internally during simulations to skip masked cells.

If `mask` was not passed to the `Output` constructor or `sim!`
it defaults to `nothing` and `false` is always returned.
"""
ismasked(data::AbstractSimData, I...) = ismasked(mask(data), I...)
ismasked(data::GridData, I...) = ismasked(mask(data), I...)
ismasked(mask::Nothing, I...) = false
ismasked(mask::AbstractArray, I...) = @inbounds !(mask[I...])

"""
    isinferred(output::Output, ruleset::Ruleset)
    isinferred(output::Output, rules::Rule...)

Test if a custom rule is inferred and the return type is correct when
`applyrule` or `applyrule!` is run.

Type-stability can give orders of magnitude improvements in performance.
"""
isinferred(output::Output, rules::Rule...) = isinferred(output, rules)
isinferred(output::Output, rules::Tuple) = isinferred(output, Ruleset(rules...))
function isinferred(output::Output, ruleset::Ruleset)
    simdata = precalcrules(SimData(output, ruleset), rules(ruleset))
    map(rules(simdata)) do rule
        isinferred(simdata, rule)
    end
    return true
end
isinferred(simdata::AbstractSimData, rule::Rule) = _isinferred(simdata, rule)
function isinferred(simdata::AbstractSimData,
    rule::Union{NeighborhoodRule,Chain{<:Any,<:Any,<:Tuple{<:NeighborhoodRule,Vararg}}}
)
    grid = simdata[neighborhoodkey(rule)]
    r = max(1, radius(rule))
    T = eltype(grid)
    S = 2r + 1
    buffer = SArray{Tuple{S,S},T,2,S^2}(Tuple(zero(T) for i in 1:S^2))
    rule = _setbuffer(rule, buffer)
    return _isinferred(simdata, rule)
end
function isinferred(simdata::AbstractSimData, rule::SetCellRule)
    rkeys, rgrids = _getreadgrids(rule, simdata)
    wkeys, wgrids = _getwritegrids(rule, simdata)
    simdata = @set simdata.grids = _combinegrids(rkeys, rgrids, wkeys, wgrids)
    readval = _readgrids(rkeys, rgrids, 1, 1)
    @inferred applyrule!(simdata, rule, readval, (1, 1))
    return true
end

function _isinferred(simdata, rule)
    rkeys, rgrids = _getreadgrids(rule, simdata)
    wkeys, wgrids = _getwritegrids(rule, simdata)
    simdata = @set simdata.grids = _combinegrids(rkeys, rgrids, wkeys, wgrids)
    readval = _readgrids(rkeys, rgrids, 1, 1)
    ex_writeval = Tuple(_example_writeval(wgrids))
    writeval = @inferred applyrule(simdata, rule, readval, (1, 1))
    typeof(Tuple(writeval)) == typeof(ex_writeval) ||
        error("return type `$(typeof(Tuple(writeval)))` doesn't match grids `$(typeof(ex_writeval))`")
    return true
end

_example_writeval(grids::Tuple) = map(_example_writeval, grids)
_example_writeval(grid::WritableGridData) = grid[1, 1]

_zerogrids(initgrid::AbstractArray, nframes) = [zero(initgrid) for f in 1:nframes]
_zerogrids(initgrids::NamedTuple, nframes) =
    [map(grid -> zero(grid), initgrids) for f in 1:nframes]

@inline _asiterable(x::Symbol) = (x,)
@inline _asiterable(x::Type{<:Tuple}) = x.parameters
@inline _asiterable(x::Tuple) = x

@inline _astuple(rule::Rule, state) = _astuple(_readkeys(rule), state)
@inline _astuple(::Tuple, state) = state
@inline _astuple(::Symbol, state) = (state,)

@inline _asnamedtuple(x::NamedTuple) = x
@inline _asnamedtuple(x::AbstractArray) = (_default_=x,)
@inline function _asnamedtuple(e::Extent)
    @set! e.init = _asnamedtuple(init(e))
    @set e.padval = _samenamedtuple(init(e), padval(e))
end

@inline _samenamedtuple(init::NamedTuple{K}, padval::NamedTuple{K}) where K = x
@noinline _samenamedtuple(init::NamedTuple{K}, padval::NamedTuple{J}) where {K,J} =
    error("Keys $K and $J do not match")
@inline _samenamedtuple(init::NamedTuple{K}, x::Tuple) where K = NamedTuple{K}(x)
@inline _samenamedtuple(init::NamedTuple, x) = map(_ -> x, init)

_unwrap(x) = x
_unwrap(::Val{X}) where X = X
_unwrap(::Type{<:Val{X}}) where X = X


text:
export state_tomo, cal_seqs

function tomo_blocks(qubits::Tuple{Vararg{Qubit}}; num_pulses::Int=4)
    if num_pulses == 4
        tomo_set = [Id, X90, Y90, X]
    elseif num_pulses == 6
        tomo_set = [Id, X90, X90m, Y90, Y90m, X]
    else
        error("Only able to handle num_pulses=4 or 6")
    end
    # TODO: replace with lexproduct when https://github.com/JuliaLang/julia/pull/18825 is merged
    pulse_combos = map(reverse, Base.product( fill(tomo_set, length(qubits))... ))
    return [ reduce(⊗, p(q) for (p,q) in zip(pulses, qubits)) for pulses in vec(pulse_combos) ]
end

"""
  state_tomo(sequence, qubits; num_pulses=4)

Applies state tomography after the `sequence` on the `qubits` using either the
4-pulse [Id, X90, Y90, X] or 6-pulse [Id, X90, X90m, Y90, Y90m, X] set of
tomographic readout pulses.
"""
function state_tomo{T<:QGL.SequenceEntry}(seq::Vector{T}, qubits::Tuple{Vararg{Qubit}}; num_pulses::Int=4)
    meas_block = reduce(⊗, MEAS(q) for q in qubits)
    return [[seq; tomo_block; meas_block] for tomo_block in tomo_blocks(qubits; num_pulses=num_pulses)]
end

"""
  cal_seqs(qubits; num_repeats=2)

Returns an array of "calibration sequences" - all combinations in
lexicographical order of Id and X pulses for the `qubits` followed by the
measurement block. To enhance the visuals or estimate noise the sequences can be
inner repeated by `num_repeats`.

# Example
```julia
julia> cal_seqs((q1,q2); num_repeats=2)
8-element Array{Array{QGL.PulseBlock,1},1}:
 QGL.PulseBlock[(Id(q1))⊗(Id(q2)),(MEAS(M-q1))⊗(MEAS(M-q2))⊗(TRIG(digitizerTrig))]
 QGL.PulseBlock[(Id(q1))⊗(Id(q2)),(MEAS(M-q1))⊗(MEAS(M-q2))⊗(TRIG(digitizerTrig))]
 QGL.PulseBlock[(Id(q1))⊗(X(q2)),(MEAS(M-q1))⊗(MEAS(M-q2))⊗(TRIG(digitizerTrig))]
 QGL.PulseBlock[(Id(q1))⊗(X(q2)),(MEAS(M-q1))⊗(MEAS(M-q2))⊗(TRIG(digitizerTrig))]
 QGL.PulseBlock[(X(q1))⊗(Id(q2)),(MEAS(M-q1))⊗(MEAS(M-q2))⊗(TRIG(digitizerTrig))]
 QGL.PulseBlock[(X(q1))⊗(Id(q2)),(MEAS(M-q1))⊗(MEAS(M-q2))⊗(TRIG(digitizerTrig))]
 QGL.PulseBlock[(X(q1))⊗(X(q2)),(MEAS(M-q1))⊗(MEAS(M-q2))⊗(TRIG(digitizerTrig))]
 QGL.PulseBlock[(X(q1))⊗(X(q2)),(MEAS(M-q1))⊗(MEAS(M-q2))⊗(TRIG(digitizerTrig))]
```
"""
function cal_seqs(qubits::Tuple{Vararg{Qubit}}; num_repeats::Int=2)
    cal_set = [Id, X]
    meas_block = reduce(⊗, MEAS(q) for q in qubits)
    # TODO: replace with lexproduct when https://github.com/JuliaLang/julia/pull/18825 is merged
    pulse_combos = map(reverse, Base.product( fill(cal_set, length(qubits))... ))
    pulse_combos = repeat(vec(pulse_combos), inner=num_repeats)
    return [ [reduce(⊗, p(q) for (p,q) in zip(pulses, qubits)), meas_block] for pulses in pulse_combos]
end

"""
Create and optionally compile a RabiWidth experiment for a single qubit assuming
a tanh pulse shape and a piAmp amplitude

RabiWidth(qubit, pulseSpacings; compile=true)

Parameters
-------
qubit           : qubit in scope
pulseSpacings   : iterable of experiment time steps in seconds
compile         : if true, compile the sequence to .h5 files

Returns
-------
seqs            : Array of sequences

Ex: RabiWidth(q, 1e-9*linspace(10, 2010, 201))
"""
function RabiWidth(qubit, pulseSpacings; compile::Bool=true)
    qubit.shape_params[:shape_function] = QGL.PulseShapes.tanh;
	seqs = [[Uθ(qubit, len, qubit.shape_params[:piAmp], 0., qubit.frequency, qubit.shape_params),
        MEAS(qubit)] for len in pulseSpacings]
    if compile
	   compile_to_hardware(seqs, "Rabi");
    end
    return seqs
end

"""
Create and optionally compile a RabiAmp experiment for a single qubit

RabiWidth(qubit, amps; compile=true)

Parameters
-------
qubit           : qubit in scope
amps            : iterable of experiment amplitudes [-1,1]
compile         : if true, compile the sequence to .h5 files

Returns
-------
seqs            : Array of sequences

Ex: RabiAmp(q, linspace(-1, 1, 101))
"""

function RabiAmp(qubit, amps; compile::Bool=true)
	seqs = [[Uθ(qubit, qubit.shape_params[:length], amp, 0.),
        MEAS(qubit)] for amp in amps]
    if compile
        compile_to_hardware(seqs, "Rabi");
    end
    return seqs
end

"""
Create and optionally compile a T1 experiment for a single qubit

InversionRecovery(qubit, amps; compile=true)

qubit           : qubit in scope
pulseSpacings   : iterable of experiment time steps in seconds
compile         : if true, compile the sequence to .h5 files

returns         : Array of sequences

Ex: InversionRecovery(q, 1e-9*linspace(10, 2010, 201))
"""

function InversionRecovery(qubit, pulseSpacings; num_repeats::Int=2, compile::Bool=true)
	seqs = [[X(qubit), Id(qubit,len), MEAS(qubit)] for len in pulseSpacings];
    cals = cal_seqs((qubit,), num_repeats = num_repeats);
    seqs = append!(seqs, cals);
    if compile
        compile_to_hardware(seqs, "T1");
    end
    return seqs
end

"""
Create and optionally compile a Ramsey experiment for a single qubit

Ramsey(qubit, pulseSpacings; TPPIFreq::Float64=0.0, num_repeats::Int=2, compile=true)

Parameters
-------
qubit           : qubit in scope
pulseSpacings   : iterable of experiment time steps in seconds
TPPIFreq        : time-proportinal phase increment frequency in Hz
num_repeats     : number of calibration repeats
compile         : if true, compile the sequence to .h5 files

Returns
-------
seqs            : Array of sequences

Ex: Ramsey(q, 1e-9*linspace(10, 2010, 201))

"""
function Ramsey(qubit, pulseSpacings; TPPIFreq::Float64=0.0, num_repeats::Int=2, compile::Bool=true)
    phases = TPPIFreq * pulseSpacings

    seqs = [[X90(qubit), Id(qubit, l), Uθ(qubit, qubit.shape_params[:length],
        qubit.shape_params[:pi2Amp], phase), MEAS(qubit)] for (l, phase) in zip(pulseSpacings, phases)];
    cals = cal_seqs((qubit,), num_repeats = num_repeats);
    seqs = append!(seqs, cals);
    if compile
        compile_to_hardware(seqs, "Ramsey");
    end
    return seqs
end


text:
@testset "Demoralization user utils" begin
  @testset "vertexsetsize" begin
    partition_int = [[1, 4], [2, 3, 5], [6], [7, 8]]
    partition_vertex = VertexSet([Vertex(block) for block = partition_int])
    @test (QSWalk.vertexsetsize(partition_vertex)) == 8
    # error test
    @test_throws MethodError QSWalk.vertexsetsize(partition_int)

  end

  @testset "default_nm_loc_ham" begin
    @test QSWalk.default_nm_loc_ham(1) == spzeros(ComplexF64, 1, 1)
    @test QSWalk.default_nm_loc_ham(3) == sparse([0. im 0.; -im 0 im; 0 -im 0])
  end

  @testset "nm_loc_ham" begin
    #default option
    @test nm_loc_ham(VertexSet([[1], [2, 3]])) == sparse([0.0+0im 0 0;0 0 im;0 -im 0])
    @test nm_loc_ham(VertexSet([[1], [2], [3]])) == spzeros(ComplexF64, 3, 3)
    #by size version
    @test nm_loc_ham(VertexSet([[1, 2], [3, 4]]), Dict(2 =>[0 1; 1 0])) == sparse([0. 1 0 0; 1 0 0 0; 0 0 0 1; 0 0 1 0])
    @test nm_loc_ham(VertexSet([[1], [2], [3]]), Dict(1 =>ones(Float64, (1, 1)))) == SparseMatrixCSC{Float64}(I, 3, 3)
    #by index version
    M1 = [1. 2;3 5]
    M2 = ones(Float64, 1, 1)
    dict = Dict(Vertex([1, 3]) =>M1, Vertex([2]) =>M2)
    @test nm_loc_ham(VertexSet([[1, 3], [2]]), dict) ==
            sparse([ 1.0+0im 0 2; 0 1 0; 3 0 5 ])
  end

  @testset "Incidence and reverse incidence lists" begin
    A = [1 2 3; 0 3. 4.; 0 0 5.]
    @test QSWalk.incidence_list(A) == [[1], [1, 2], [1, 2, 3]]
    @test QSWalk.incidence_list(sparse(A)) == QSWalk.incidence_list(A)
    @test QSWalk.incidence_list(A; epsilon = 2.5) == [Int64[], [2], [1, 2, 3]]
    @test QSWalk.incidence_list(sparse(A)) == QSWalk.incidence_list(A)
    @test QSWalk.reversed_incidence_list(A) == [[1, 2, 3], [2, 3], [3]]
    @test QSWalk.reversed_incidence_list(sparse(A)) == QSWalk.reversed_incidence_list(A)
    @test QSWalk.reversed_incidence_list(A; epsilon = 2.5) == [[3], [2, 3], [3]]
    @test QSWalk.reversed_incidence_list(sparse(A); epsilon = 2.5) == QSWalk.reversed_incidence_list(A; epsilon = 2.5)
    #errors
    @test_throws ArgumentError QSWalk.reversed_incidence_list(A, epsilon = -1)
    @test_throws ArgumentError QSWalk.incidence_list(A, epsilon = -1)
    @test_throws ArgumentError QSWalk.reversed_incidence_list(sparse(A), epsilon = -1)
    @test_throws ArgumentError QSWalk.incidence_list(sparse(A), epsilon = -1)
  end

  @testset "vertexset creation" begin
    @test QSWalk.revinc_to_vertexset([[1, 3], [2, 3], Int[], [4, 6, 1]]) == VertexSet([[1, 2], [3, 4], [5], [6, 7, 8]])
    @test make_vertex_set([1 2 3; 0 3. 4.; 0 0 5.]) == VertexSet([[1, 2, 3], [4, 5], [6]])
  end

  @testset "Fourier matrix" begin
    @test QSWalk.fourier_matrix(4) ≈ [1 1 1 1; 1 1im -1 -1im;1 -1 1 -1; 1 -1im -1 1im]
    @test QSWalk.fourier_matrix(2) ≈ [1 1; 1 -1.]
    @test QSWalk.fourier_matrix(1) ≈ ones(Float64, (1, 1))
  end

  @testset "nm_lind" begin
    A = sparse([0.0+0.0im 1 0; 1 0 1; 0 1 0])
    B1, B2 = ones(1, 1), ones(2, 2)
    #default
    #needs to be roughly, since exp computing is inexact
    @test nm_lind(A)[1] ≈ [0 1 1 0;
                           1 0 0 1;
                           1 0 0 -1;
                           0 1 1 0]
    @test nm_lind(A)[2] == make_vertex_set(A)

    A = sparse([0.0+0.0im 0 0 0 1;
                0 0 1 0 1;
                0 0 0 0 0;
                0 1 1 0 0;
                0 0 0 0 0])
    L, vset = nm_lind(A)
    v1, v2, v3, v4, v5 = vlist(vset)
    @test nm_lind(A)[1] ≈ [0 0 0 0 0 0 1;
                                             0 0 0 1 0 0 1;
                                             0 0 0 1 0 0 -1;
                                             0 0 0 0 0 0 0;
                                             0 1 1 1 0 0 0;
                                             0 1 1 -1 0 0 0;
                                             0 0 0 0 0 0 0]
    @test nm_lind(A)[2] == make_vertex_set(A)
    @test nm_lind(A, Dict(v1 => B1, v2 => 2*B2, v3 => 3*B1, v4 => 4*B2, v5 => 5*B1))[1] ≈
             [0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  1.0+0.0im;
              0.0+0.0im  0.0+0.0im  0.0+0.0im  2.0+0.0im  0.0+0.0im  0.0+0.0im  2.0+0.0im;
              0.0+0.0im  0.0+0.0im  0.0+0.0im  2.0+0.0im  0.0+0.0im  0.0+0.0im  2.0+0.0im;
              0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im;
              0.0+0.0im  4.0+0.0im  4.0+0.0im  4.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im;
              0.0+0.0im  4.0+0.0im  4.0+0.0im  4.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im;
              0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im]

    A = [0 0 0; 0 0 0; 1 1 0]
    @test nm_lind(A)[1] ≈ [0 0 0 0;
                                             0 0 0 0;
                                             1 1 0 0;
                                             1 -1 0 0]
    @test nm_lind(A)[2] == make_vertex_set(A)

    @test nm_lind(A, Dict(1  => B1, 2 =>3*B2 ))[1] ≈
             [0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im;
              0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im;
              3.0+0.0im  3.0+0.0im  0.0+0.0im  0.0+0.0im;
              3.0+0.0im  3.0+0.0im  0.0+0.0im  0.0+0.0im]
    @test nm_lind(A, Dict(1  => B1, 2 =>3*B2 ))[2] ==
      make_vertex_set(A)

    @test nm_lind(A, Dict(1  => B1, 2 =>3*B2 ))[1] ≈
             [0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im;
              0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im;
              3.0+0.0im  3.0+0.0im  0.0+0.0im  0.0+0.0im;
              3.0+0.0im  3.0+0.0im  0.0+0.0im  0.0+0.0im]

    L, vset = nm_lind(A)
    v1, v2, v3 = vlist(vset)
    @test nm_lind(A, Dict(v1 => B1, v2 => 2*B1, v3 => 3*B2))[1] ≈
             [0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im;
              0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im;
              3.0+0.0im  3.0+0.0im  0.0+0.0im  0.0+0.0im;
              3.0+0.0im  3.0+0.0im  0.0+0.0im  0.0+0.0im]
    @test nm_lind(A, Dict(v1 => B1, v2 => 2*B1, v3 => 3*B2))[2] ==
      make_vertex_set(A)

      A = [2 0 0 3;
    im 0 3im 0;
    1 0 0 0;
    im -im 0 0]
    @test nm_lind(A,Dict(1=>ones(1,1), 2=> ones(2,2)))[1] ≈
              [2.0+0.0im 2.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 3.0+0.0im 3.0+0.0im;
               2.0+0.0im 2.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 3.0+0.0im 3.0+0.0im;
               0.0+1.0im 0.0+1.0im 0.0+0.0im 0.0+0.0im 0.0+3.0im 0.0+0.0im 0.0+0.0im;
               0.0+1.0im 0.0+1.0im 0.0+0.0im 0.0+0.0im 0.0+3.0im 0.0+0.0im 0.0+0.0im;
               1.0+0.0im 1.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im;
               0.0+1.0im 0.0+1.0im 0.0-1.0im 0.0-1.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im;
               0.0+1.0im 0.0+1.0im 0.0-1.0im 0.0-1.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im]

  end

  @testset "nm_glob_ham" begin
    A = sparse([0 1 0;
               1 0 2;
               0 2 0])
    #default
    #needs to be roughly, since exp computing is inexact
    @test nm_glob_ham(A) ≈ [0 1 1 0;
                                   1 0 0 2;
                                   1 0 0 2;
                                   0 2 2 0]
    @test nm_glob_ham(A, Dict((1, 2) => (2+1im)*ones(1, 2), (2, 1) =>1im*ones(2, 1))) ≈
                                    [0 2+1im 2+1im 0;
                                     2-1im 0 0 2im;
                                     2-1im 0 0 2im;
                                     0 -2im -2im 0]

    v1, v2, v3 = vlist(make_vertex_set(A))
    @test nm_glob_ham(A, Dict{Tuple{Vertex,Vertex},AbstractMatrix{<:Number}}((v1, v2) =>2*ones(1, 2), (v2, v3) => -adjoint([1im 2im]))) ==
                                              [0 2 2 0;
                                               2 0 0 2im;
                                               2 0 0 4im;
                                               0 -2im -4im 0]
   A = [0  0  1;
        0  0  1;
        2  2  0]
  v1, v2, v3 = vlist(make_vertex_set(A))
  @test nm_glob_ham(A, Dict{Tuple{Vertex,Vertex},AbstractMatrix{<:Number}}((v1, v3) =>2*ones(1, 2), (v2, v3) =>[1im 2im;])) ≈
                          [0.0+0.0im  0.0+0.0im  3.0+0.0im  3.0+0.0im;
                           0.0+0.0im  0.0+0.0im  0.0+1.5im  0.0+3.0im;
                           3.0+0.0im  0.0-1.5im  0.0+0.0im  0.0+0.0im;
                           3.0+0.0im  0.0-3.0im  0.0+0.0im  0.0+0.0im]

  @test nm_glob_ham(A, Dict{Tuple{Vertex,Vertex},AbstractMatrix{<:Number}}((v1, v3) =>2*ones(1, 1), (v2, v3) =>ones(1,1)*1im),epsilon=3) ≈
                          [0.0+0.0im  0.0+0.0im  0.0+0.0im;
                           0.0+0.0im  0.0+0.0im  0.0+0.0im;
                           0.0+0.0im  0.0+0.0im  0.0+0.0im]
  a_directed = diagm(1=>ones(3))
  a_directed_transposed = diagm(-1=>ones(3))
  a_undirected = a_directed_transposed + a_directed
  @test nm_glob_ham(a_directed) ≈ a_undirected
  @test nm_glob_ham(a_directed_transposed) ≈ a_undirected
  @test_throws DimensionMismatch nm_glob_ham(a_undirected) ≈ a_undirected

  @test nm_glob_ham(Matrix(a_directed)) ≈ a_undirected
  @test nm_glob_ham(Matrix(a_directed_transposed)) ≈ a_undirected
  @test_throws DimensionMismatch nm_glob_ham(Matrix(a_undirected)) ≈ a_undirected
  end
end

@testset "Cannonical measurement" begin
  probability = [0.05, 0.1, 0.25, 0.3, 0.01, 0.20, 0.04, 0.05]
  vset = VertexSet([[1, 4], [2, 3, 5], [6], [7, 8]])
  result = [0.35, 0.36, 0.2, 0.09]
  @test nm_measurement(probability, vset) == result

  state = [1/6 1.0 1/6; 1. 2/3 1im; 1/6 -1im 1/6]
  vset = VertexSet([[1, 3], [2]])
  @test nm_measurement(state, vset) ≈ [1.0/3, 2.0/3]
  @test_throws AssertionError nm_measurement(Matrix{Float64}(I,2,2)/2., vset)
  @test_throws AssertionError nm_measurement([1.0/2, 1.0/3, 1.0/6, 0.], vset)

end

@testset "Initial states creation" begin
  vset = VertexSet([[1], [2, 3, 4], [5, 6, 7], [8, 9]])
  @test nm_init(vset[[1, 3, 4]], vset) ≈ spdiagm(0=>[1.0/3, 0, 0, 0, 1.0/9, 1.0/9, 1.0/9, 1.0/6, 1.0/6])
  A1 = ones(ComplexF64, 1, 1)/4
  A2 = [ 1/5+0im 0 1/5; 0 1/10 0 ; 1/5 0 1/5 ]
  A3 = [0.125 -0.125+0im; -0.125 0.125]
  dict = Dict(vset[1] =>A1, vset[3] =>A2, vset[4] =>A3 )
  @test nm_init(dict, vset) ==
        sparse([ 0.25+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im     0.0+0.0im     0.0+0.0im
                0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im     0.0+0.0im     0.0+0.0im
                0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im     0.0+0.0im     0.0+0.0im
                0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im     0.0+0.0im     0.0+0.0im
                0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.2+0.0im  0.0+0.0im  0.2+0.0im     0.0+0.0im     0.0+0.0im
                0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.1+0.0im  0.0+0.0im     0.0+0.0im     0.0+0.0im
                0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.2+0.0im  0.0+0.0im  0.2+0.0im     0.0+0.0im     0.0+0.0im
                0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im   0.125+0.0im  -0.125+0.0im
                0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  -0.125+0.0im   0.125+0.0im])

  dictwrong = Dict(vset[1] =>A2, vset[3] =>A2, vset[4] =>A3 )
  @test_throws AssertionError nm_init(vlist(VertexSet([[1, 2, 3], [4, 5]])), vset)
  @test_throws AssertionError nm_init(dictwrong, vset)
end


text:
@testset "Z4Group, Z4Group0 use" begin
    t1 = PauliTerm("XX", Z4Group0(1))
    t2 = PauliTerm("XX", 1)
    (m1, m2) = Matrix.((t1, t2))
    @test m1 == m2
    @test eltype(m1) == eltype(m2) == Float64
end

@testset "Pauli" begin
    for (n, s) in ((0, :I), (1, :X), (2, :Y), (3, :Z))
        @test Pauli(n) == Pauli(s)
        st = string(s)
        @test Pauli(n) == Pauli(st)
        @test Pauli(n) == Pauli(only(st)) # character
    end

    for (a, b, c) in ((:I, :I, :I), (:I, :X, :X), (:I, :Y, :Y), (:I, :Z, :Z),
                      (:X, :I, :X), (:X, :X, :I), (:X, :Y, :Z), (:X, :Z, :Y),
                      (:Y, :I, :Y), (:Y, :X, :Z), (:Y, :Y, :I), (:Y, :Z, :X),
                      (:Z, :I, :Z), (:Z, :X, :Y), (:Z, :Y, :X), (:Z, :Z, :I))
        @test Pauli(a) * Pauli(b) == Pauli(c)
    end

    m = rand(Pauli, (2, 2))
    @test typeof(m) == Matrix{Pauli}
    @test size(m) == (2, 2)
end

@testset "PauliTerm" begin
    (I, X, Y, Z) = Pauli.(0:3)
    @test X ⊗ Y == PauliTerm("XY")
    @test PauliTerm("XY", 3).coeff == 3

    for (a, b, c, phase) in ((:I, :I, :I, 1), (:I, :X, :X, 1), (:I, :Y, :Y, 1), (:I, :Z, :Z, 1),
                             (:X, :I, :X, 1), (:X, :X, :I, 1), (:X, :Y, :Z, im), (:X, :Z, :Y, -im),
                             (:Y, :I, :Y, 1), (:Y, :X, :Z, -im), (:Y, :Y, :I, 1), (:Y, :Z, :X, im),
                             (:Z, :I, :Z, 1), (:Z, :X, :Y, im), (:Z, :Y, :X, -im), (:Z, :Z, :I, 1))
        (as, bs, cs) = String.((a, b, c))
        @test PauliTerm(as) * PauliTerm(bs) == phase * PauliTerm(cs)
    end
    @test PauliTerm("XYZ") * PauliTerm("ZYX") == PauliTerm("YIY")
    @test PauliTerm("XYZX") * PauliTerm("ZYXZ") == -im * PauliTerm("YIYY")
    @test PauliTerm("ZYXZ") * PauliTerm("XYZX") == im * PauliTerm("YIYY")
    t = PauliTerm("ZYXZ", 3)
    @test t^3 == t.coeff^2 * t
    @test t^4 == t.coeff^4 * one(t)
    @test PauliTerm("IXYZ") == PauliTerm([0, 1, 2, 3])
    @test PauliTerm("IXYZ", 2) == PauliTerm([0, 1, 2, 3], 2)
    @test PauliTerm(X, Y, Z) == PauliTerm("XYZ")
    @test PauliTerm(:XYZ) == PauliTerm("XYZ")

    @test sqrt(PauliTerm(X,Y)) == PauliSum(["II", "XY"], [0.5 + 0.5im, 0.5 - 0.5im])

    @test adjoint(PauliTerm("XYZ")) == PauliTerm("XYZ")
    @test conj(PauliTerm("XYZ")) ==  PauliTerm("XYZ", -1)
    @test transpose(PauliTerm("XYZ")) == PauliTerm("XYZ", -1)
    @test transpose(PauliTerm("XYZY")) == PauliTerm("XYZY")
end

@testset "PauliSum" begin
    ps = PauliSum(["YYY", "XXX"], [2, 3])
    @test length(ps) == 2
    @test [t for t in ps] == [PauliTerm("XXX", 3), PauliTerm("YYY", 2)]
    add!(ps, PauliTerm("III", -1))
    @test ps[1] == PauliTerm("III", -1)
    @test length(ps) == 3
    add!(ps, PauliTerm("XXX", -3))
    @test length(ps) == 2

    # n=2 tests #51
    for T in (Float64, ComplexF64)
        for n in (2, 4, 8, 16, 32)
            for i in 1:10
                m = rand(T, n, n)
                s = PauliSum(m)
                @test m ≈ Matrix(s)
            end
        end
    end

    term = PauliTerm("XYZ")
    tan_term = tan(term)
    @test first(tan_term)  == tan(1) * term
    @test length(tan_term) == 1
    @test numeric_function(PauliTerm("XX"), x -> x^2 + 2.0 * x^3) == PauliSum([PauliTerm("II"), PauliTerm("XX", 2)])

    a = rand_op_sum(Pauli, 10, 3)
    b = rand_op_sum(Pauli, 10, 3)
    sa = SparseArrays.sparse(a)
    sb = SparseArrays.sparse(b)
    @test SparseArrays.sparse(a * b * a) == sa * sb * sa

    @test all(commutes, group_paulis(rand_op_sum(Pauli, 10, 100)))

end

@testset "Z4Group0 Pauli" begin
    ms = Matrix.(Pauli.((0, 1, 2, 3)))
    zms = ((Z4Group0.(m) for m in ms)...,)
    @test zms == ms
    @test kron(zms...) == kron(ms...)
end

@testset "kron" begin
    t1 = PauliTerm("XYZ")
    t2 = PauliTerm("ZYX")
    @test Matrix(t1) ⊗ Matrix(t2) ≈ Matrix(t1 ⊗ t2)

    t1 = 3.0 * Paulis.X ⊗ Paulis.Y
    t2 = Paulis.X ⊗ Paulis.Y * 3.0
    @test t1 == t2
    @test_broken typeof(t1.coeff) == typeof(t2.coeff)
end

import LinearAlgebra
using LinearAlgebra: eigvals

@testset "eigvals" begin
    @test eigvals(PauliTerm("III")) == ones(2^3)
    for t in (PauliTerm("XYZX"), PauliTerm("XYI"))
        @test eigvals(t) ≈ eigvals(Matrix(t))
    end
end

@testset "properties" begin
    @test LinearAlgebra.tr(Paulis.I) == 2
    @test LinearAlgebra.tr(Paulis.X) == 0
    @test LinearAlgebra.tr(Paulis.Y) == 0
    @test LinearAlgebra.tr(Paulis.Z) == 0
end


text:
Nemo.fit!(::fmpq_rel_series, Int) = nothing
Nemo.fit!(::fmpq_abs_series, Int) = nothing

@doc Markdown.doc"""
    integral(f::RelSeriesElem{T}) -> RelSeriesElem

Return the integral of the power series $f$.
"""
function Nemo.integral(f::RelSeriesElem{T}) where T
  g = parent(f)()
  fit!(g, precision(f)+1)
  set_precision!(g, precision(f)+1)
  v = valuation(f)
  Nemo.set_valuation!(g, v+1)
  for i=0:Nemo.pol_length(f)
    c = Nemo.polcoeff(f, i)
    if !iszero(c)
      setcoeff!(g, i, divexact(c, i+v+1))
    end
  end
  Nemo.renormalize!(g)
  return g
end


text:
module lorenz
using ReachabilityAnalysis
@taylorize function lorenz!(dx, x, params, t)
    σ, β, ρ = 10.0, 8.0/3.0, 28.0
    dx[1] = σ * (x[2] - x[1])
    dx[2] = x[1] * (ρ - x[3]) - x[2]
    dx[3] = x[1] * x[2] - β * x[3]
    return dx
end

function model(X0)
    S = @system(x' = lorenz!(x), dim:3)
    return IVP(S, X0)
end
end # module

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary


text:
module CUBLAS

using ..APIUtils

using ..CUDA
using ..CUDA: CUstream, cuComplex, cuDoubleComplex, libraryPropertyType, cudaDataType
using ..CUDA: libcublas, unsafe_free!, @retry_reclaim, isdebug, @sync, @context!

using GPUArrays

using LinearAlgebra

using BFloat16s

using CEnum

using Memoization

using DataStructures

# core library
include("libcublas_common.jl")
include("error.jl")
include("libcublas.jl")
include("libcublas_deprecated.jl")

# low-level wrappers
include("util.jl")
include("wrappers.jl")

# high-level integrations
include("linalg.jl")

function math_mode!(handle, mode)
    flags = 0

    # https://github.com/facebookresearch/faiss/issues/1385
    if version(handle) > v"11"
        flags = CUBLAS_MATH_DISALLOW_REDUCED_PRECISION_REDUCTION
    end

    flags |= if mode == CUDA.PEDANTIC_MATH
        # prevent use of tensor cores
        if version(handle) < v"11"
            CUBLAS_DEFAULT_MATH
        else
            CUBLAS_PEDANTIC_MATH
        end
    elseif mode == CUDA.DEFAULT_MATH
        # use tensor cores, but don't reduce precision
        if version(handle) < v"11"
            CUBLAS_TENSOR_OP_MATH
        else
            CUBLAS_DEFAULT_MATH
        end
    elseif mode == CUDA.FAST_MATH
        # we'll additionally select a compute-mode with reduced precision whenever possible
        if version(handle) < v"11"
            CUBLAS_TENSOR_OP_MATH
        else
            CUBLAS_TF32_TENSOR_OP_MATH
        end
    end

    cublasSetMathMode(handle, flags)

    return
end

# cache for created, but unused handles
const idle_handles = HandleCache{CuContext,cublasHandle_t}()
const idle_xt_handles = HandleCache{Any,cublasXtHandle_t}()

function handle()
    state = CUDA.active_state()
    handle, stream, math_mode = get!(task_local_storage(), (:CUBLAS, state.context)) do
        new_handle = pop!(idle_handles, state.context) do
            cublasCreate()
        end

        finalizer(current_task()) do task
            push!(idle_handles, state.context, new_handle) do
                @context! skip_destroyed=true state.context cublasDestroy_v2(handle)
            end
        end

        cublasSetStream_v2(new_handle, state.stream)

        math_mode!(new_handle, state.math_mode)

        new_handle, state.stream, state.math_mode
    end::Tuple{cublasHandle_t,CuStream,CUDA.MathMode}

    if stream != state.stream
        cublasSetStream_v2(handle, state.stream)
        task_local_storage((:CUBLAS, state.context), (handle, state.stream, math_mode))
        stream = state.stream
    end

    if math_mode != state.math_mode
        math_mode!(handle, state.math_mode)
        task_local_storage((:CUBLAS, state.context), (handle, stream, state.math_mode))
        math_mode = state.math_mode
    end

    return handle
end

function xt_handle()
    ctxs = Tuple(context(dev) for dev in devices())
    get!(task_local_storage(), (:CUBLASxt, ctxs)) do
        handle = pop!(idle_xt_handles, ctxs) do
            cublasXtCreate()
        end

        finalizer(current_task()) do task
            push!(idle_xt_handles, ctxs, handle) do
                # TODO: which context do we need to destroy this on?
                cublasXtDestroy(handle)
            end
        end

        # select all devices
        devs = convert.(Cint, devices())
        cublasXtDeviceSelect(handle, length(devs), devs)

        handle
    end::cublasXtHandle_t
end

function log_message(cstr)
    # NOTE: we can't `@debug` these messages, because the logging function is called for
    #       every line... also not sure what the i!/I! prefixes mean (info?)
    # NOTE: turns out we even can't `print` these messages, as cublasXt callbacks might
    #       happen from a different thread! we could strdup + uv_async_send, but I couldn't
    #       find an easy way to attach data to that
    # TODO: use a pre-allocated lock-free global message buffer?
    len = ccall(:strlen, Csize_t, (Cstring,), cstr)
    ccall(:write, Cint, (Cint, Cstring, Csize_t), 0, cstr, len)
    return
end

function __runtime_init__()
    # enable library logging when launched with JULIA_DEBUG=CUBLAS
    if isdebug(:init, CUBLAS)
        callback = @cfunction(log_message, Nothing, (Cstring,))
        cublasSetLoggerCallback(callback)
    end
end

end


text:
using Documenter
using Latexify
using Markdown
using Markdown: MD, Paragraph, LineBreak
using Plots

using GeometricIntegrators
using GeometricIntegrators.Integrators.VPRK

using GeometricProblems.Diagnostics
using GeometricProblems.PlotRecipes

function Integrators.integrate!(int::Union{IntegratorVPRK{DT,TT,D,S},IntegratorFIRKimplicit{DT,TT,D,S},IntegratorSRKimplicit{DT,TT,D,S}}, sol::Solution{DT,TT}) where {DT,TT,D,S}
    asol = AtomicSolution(int)

    # loop over initial conditions showing progress bar
    for m in eachsample(sol)
        get_initial_conditions!(sol, asol, m, 1)
        initialize!(int, asol)

        n = 0
        try
            for outer n in eachtimestep(sol)
                integrate!(int, sol, asol, m, n)
            end
        catch ex
            tstr = " in time step " * string(n)

            if nsamples(sol) > 1
                tstr *= " for initial condition " * string(m)
            end

            tstr *= "."

            if isa(ex, DomainError)
                show(stdout, "text/markdown", Markdown.parse("DOMAIN ERROR: Simulation crashed" * tstr))
                @warn("DOMAIN ERROR: Simulation crashed" * tstr)
            elseif isa(ex, NonlinearSolverException)
                show(stdout, "text/markdown", Markdown.parse("SOLVER ERROR: Simulation crashed" * tstr))
                show(stdout, "text/markdown", Markdown.parse(ex.msg))
                @warn("SOLVER ERROR: Simulation crashed" * tstr)
                @warn(ex.msg)
            else
                show(stdout, "text/markdown", Markdown.parse("ERROR: Simulation crashed with $(typeof(ex)) " * tstr))
                show(stdout, "text/markdown", Markdown.parse(ex.msg))
                @warn("ERROR: Simulation crashed with $(typeof(ex)) " * tstr)
                @warn(ex.msg)
                # showerror(stdout, ex, catch_backtrace())
            end
        end
    end

    nothing
end

_linebreak(io) = show(io, "text/markdown", MD(Paragraph([LineBreak()])))

function _plot_figure_md(file, name, filename)
    # if isfile(filename)
        show(file, "text/markdown", Markdown.parse("![$name]($filename)"))
        _linebreak(file)
    # else
    #     show(stdout, "text/markdown", Markdown.parse("ERROR: Plot output $filename does not exist!"))
    #     @warn("Plot output $filename does not exist!")
    # end
end

function write_plots(dir, file, name, fig_suff)

    plot_file = file * ".md"

    open(plot_file, "w") do f
        show(f, "text/markdown", Markdown.parse("# $name"))
        _linebreak(f)

        _plot_figure_md(f, name, "$(dir)/$(file)_solution$(fig_suff)")
        _plot_figure_md(f, name, "$(dir)/$(file)_traces$(fig_suff)")

        show(f, "text/markdown", Markdown.parse("## Energy Error"))
        _linebreak(f)

        _plot_figure_md(f, name, "$(dir)/$(file)_energy_error$(fig_suff)")
        _plot_figure_md(f, name, "$(dir)/$(file)_energy_drift$(fig_suff)")

        show(f, "text/markdown", Markdown.parse("## Constraint"))
        _linebreak(f)

        _plot_figure_md(f, name, "$(dir)/$(file)_constraint_error$(fig_suff)")
    end
end

function plot(sol, equ, dir, file, fig_suff)
    if !isdir(dir)
        mkdir(dir)
    end

    try
        plot_lotka_volterra_2d(sol, equ; nt=lastentry(sol), fmt=:png)
        savefig(dir * "/" * file * fig_suff)

        plot_lotka_volterra_2d_solution(sol, equ; nt=lastentry(sol), fmt=:png)
        savefig(dir * "/" * file * "_solution" * fig_suff)

        plot_lotka_volterra_2d_traces(sol, equ; nt=lastentry(sol), fmt=:png)
        savefig(dir * "/" * file * "_traces" * fig_suff)

        H, ΔH = compute_energy_error(sol.t, sol.q, equ.parameters)
        plotenergyerror(sol.t, ΔH; nt=lastentry(sol), fmt=:png)
        savefig(dir * "/" * file * "_energy_error" * fig_suff)

        plotenergydrift(compute_error_drift(sol.t, ΔH, div(ntime(sol), 10))...; nt=lastentry(sol))
        savefig(dir * "/" * file * "_energy_drift" * fig_suff)

        plotconstrainterror(sol.t, compute_momentum_error(sol.t, sol.q, sol.p, (t,q,k)->ϑ(t,q,equ.parameters,k)); nt=lastentry(sol), fmt=:png)
        savefig(dir * "/" * file * "_constraint_error" * fig_suff)
    catch ex
        if isa(ex, DomainError)
            show(stdout, "text/markdown", Markdown.parse("DOMAIN ERROR: Diagnostics crashed."))
            @warn("DOMAIN ERROR: Diagnostics crashed.")
        else
            show(stdout, "text/markdown", Markdown.parse("ERROR: Plot crashed with $(typeof(ex))"))
            show(stdout, "text/markdown", Markdown.parse(ex.msg))
            @warn("ERROR: Plot crashed with $(typeof(ex))")
            @warn(ex.msg)
            showerror(stdout, ex, catch_backtrace())
        end
    end
end

function run_integrator(iode, tab, integrator, nt, dir, file, fig_suff)
    sol = Solution(iode, Δt, nt)
    int = integrator(iode, tab, Δt)

    integrate!(int, sol)

    plot(sol, iode, dir, file, fig_suff)
end

function run_list(iode, name, list, plot_dir = PLOT_DIR, symp_dir = SYMP_DIR;
                    fig_suff = ".png")

    for run in list
        tab, file = run

        if length(run) ≥ 3
            integrator = run[3]
        else
            integrator = Integrator
        end

        write_plots(plot_dir, file, name, fig_suff)

        show(stdout, "text/markdown", Markdown.parse("## $(tab.name)"))
        _linebreak(stdout)

        show(stdout, "text/markdown", Markdown.parse("[Plots]($file.md)"))
        _linebreak(stdout)

        run_integrator(iode, tab, integrator, nt, plot_dir, file, fig_suff)
        show(stdout, "text/markdown", Markdown.parse("![$name]($plot_dir/$file$fig_suff)"))
    end

    nothing
end


text:
using Documenter, GMMInference, DocumenterMarkdown

runweave=true
runnotebook=false
if runweave
  using Weave
  wd = pwd()
  try
    builddir=joinpath(dirname(Base.pathof(GMMInference)),"..","docs","build")
    mkpath(builddir)
    cd(builddir)
    jmdfiles = filter(x->occursin(r".jmd$",x), readdir(joinpath("..","jmd")))
    for f in jmdfiles
      src = joinpath("..","jmd",f)
      target = joinpath("..","build",replace(f, r"jmd$"=>s"md"))
      if stat(src).mtime > stat(target).mtime
        weave(src,out_path=joinpath("..","build"),
              cache=:user,
              cache_path=joinpath("..","weavecache"),
              doctype="github", mod=Main,
              args=Dict("md" => true))
      end
      target = joinpath("..","build",replace(f, r"jmd$"=>s"ipynb"))
      if (runnotebook && stat(src).mtime > stat(target).mtime)
          notebook(src,out_path=joinpath("..","build"),
                   nbconvert_options="--allow-errors")
      end
    end
  finally
    cd(wd)
  end
  if (isfile("build/temp.md"))
    rm("build/temp.md")
  end
end

makedocs(
  modules=[GMMInference],
  format=Markdown(),
  clean=false,
  pages=[
    "Home" => "index.md", # this won't get used anyway; we use mkdocs instead for interoperability with weave's markdown output.
  ],
  repo="https://github.com/schrimpf/GMMInference.jl/blob/{commit}{path}#L{line}",
  sitename="GMMInference.jl",
  authors="Paul Schrimpf <paul.schrimpf@gmail.com>",
)

run(`mkdocs build`)

#deploydocs(;
#    repo="github.com/schrimpf/GMMInference.jl",
#)

deploy=true
if deploy || "deploy" in ARGS
  run(`mkdocs gh-deploy`)
end


text: open(readlines, "path/to/file")[7]

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Clang.Generators

# using ImPlot.LibCImPlot.CImPlot_jll

using CImGui.CImGui_jll

include_dir = joinpath(CImGui_jll.artifact_dir, "include")

cd(@__DIR__)

const CIMPLOT_H = joinpath(@__DIR__, "cimplot_patched.h") |> normpath

options = load_options(joinpath(@__DIR__, "generator.toml"))

args = ["-I$include_dir", "-DCIMGUI_DEFINE_ENUMS_AND_STRUCTS"]

# add definitions
@add_def ImVec2
@add_def ImVec4
@add_def ImGuiMouseButton
@add_def ImGuiKeyModFlags
@add_def ImS8
@add_def ImU8
@add_def ImS16
@add_def ImU16
@add_def ImS32
@add_def ImU32
@add_def ImS64
@add_def ImU64
@add_def ImTextureID
@add_def ImGuiCond
@add_def ImGuiDragDropFlags
@add_def ImDrawList
@add_def ImGuiContext

#=
jltypes = [Float32, Float64, Int8, UInt8, Int16, UInt16, Int32, UInt32, Int64, UInt64]
typedict = Dict(zip(imtypes,jltypes))
=#
#type_names = ["FloatPtr", "doublePtr", "S8Ptr", "U8Ptr", "S16Ptr", "U16Ptr", "S32Ptr", "U32Ptr", "S64Ptr", "U64Ptr"]

imdatatypes = [:Cfloat, :Cdouble, :ImS8, :ImU8, :ImS16, :ImU16, :ImS32, :ImU32, :ImS64, :ImU64]
plot_types = ["Line", "Scatter", "Stairs", "Shaded", "BarsH", "Bars", "ErrorBarsH", "ErrorBars", "Stems", "VLines", "HLines", "PieChart", "Heatmap", "Histogram", "Histogram2D", "Digital"]

ctx = create_context(CIMPLOT_H, args, options)

build!(ctx, BUILDSTAGE_NO_PRINTING)

function revise_function!(e::Expr)

    # Get name of the function
    fname = string(e.args[1].args[1])

    # Skip if it's not a prefix added by cimplot
    startswith(fname,"ImPlot_") || return e

    # Strip off the prefix to match C++ (since we have a namespace)
    fname = fname[8:end] # remove first 7 characters == 'ImPlot_'

    # Plot functions are templated and have a regular structure
    if startswith(fname, "Plot")

        # iterate the tuple of C argument types in the function's ccall
        cargtypes = e.args[2].args[1].args[4].args
        for (i, cargtype) in enumerate(cargtypes)

            # All pointers in plot functions (where eltype ∈ imdatatypes) are array inputs
            if (Meta.isexpr(cargtype, :curly) && cargtype.args[1] == :Ptr) && cargtype.args[2] ∈ imdatatypes

                # Input arrays allocated in Julia should be passed as Ref: https://docs.julialang.org/en/v1/manual/calling-c-and-fortran-code/#When-to-use-T,-Ptr{T}-and-Ref{T}
                # e.args[2].args[1].args[4].args[i].args[1] = :Ref

                # Annotate the Julia method signature
                sym = e.args[1].args[i+1]
                e.args[1].args[i+1] = Expr(:(::), sym, Expr(:curly, :Union, Expr(:curly, :AbstractArray, cargtype.args[2]), Expr(:curly, :Ref, cargtype.args[2]), Expr(:curly, :Ptr, cargtype.args[2])))
                # Used if you want to calculate stride size below
                #= if cargtype.args[2] ∈ keys(typedict)
                    global stridesize = sizeof(typedict[cargtype.args[2]])
                end =#
            elseif cargtype == :Cint
                sym = e.args[1].args[i+1]
                #= it's possible to set default arg values from DAG expressions too. Examples below are unused for now.
                if sym == :stride
                    e.args[1].args[1+i] = Expr(:call, :(=), Expr(:(::), sym, :Integer), Expr(:call, :Int32, stridesize))
                elseif sym == :offset
                    e.args[1].args[1+i] = Expr(:call, :(=), Expr(:(::), sym, :Integer), Expr(:call, :zero, :Int32))
                else =#
                    e.args[1].args[1+i] = Expr(:(::), sym, :Integer)
                # end
            elseif cargtype ∈ [:Cdouble, :Cfloat]
                sym = e.args[1].args[i+1]
                e.args[1].args[i+1] = Expr(:(::), sym, :Real)
            end
        end

        for ptype in plot_types
            fullname = "Plot" * ptype
            if startswith(fname, fullname)
                if fname[length(fullname)+1] == 'G' || fname[length(fullname)+1] == 'H'
                    fname = fname[1:length(fullname)+1]
                else
                    fname = fname[1:length(fullname)]
                end
            end
        end
    end

    e.args[1].args[1] = Symbol(fname)
    return e
end

function rewrite!(dag::ExprDAG)
    for node in get_nodes(dag)
        for expr in get_exprs(node)
            Meta.isexpr(expr, :function) && revise_function!(expr)
        end
    end
end

rewrite!(ctx.dag)
build!(ctx, BUILDSTAGE_PRINTING_ONLY)


text:
function read_mod_file(filepath::AbstractString)

    println("Reading $(filepath) ...")

    if isfile(filepath)
        f = open(filepath, "r")
    else
        error("No .mod file detected.")
    end

    skip = 0
    comment = 0
    parsed = 0
    consregion = false

    mod = oneProblem()

    for l in eachline(f)
        sl = lstrip(l, ' ')
        if isempty(sl) || length(sl) == 1  # Skip empty line
            skip += 1
        elseif l[1] == '#'
            comment += 1
        else
            slt = split(sl, r" |,")
            slt = [slt[i] for i in 1:length(slt) if length(slt[i]) > 0]
            if slt[1] == "var"
                @assert !consregion
                mod_read_var(f, mod, l)
            elseif slt[1] in ["minimize", "maximize"]
                @assert !consregion
                mod_read_objective(f, mod, l, string(slt[1]))
            elseif ismatch(r"e\d+:", slt[1])
                @assert consregion
                mod_read_equation(f, mod, l)
            elseif slt[1] == "subject"
                consregion = true
            end
            parsed += 1
        end
    end
    close(f)

    println("Finish reading mod file PARSED $(parsed) || COMMENT $(comment) || EMPTY $(skip)")

    return mod
end

function mod_read_var(file::IOStream, mod::oneProblem, lInit::AbstractString)

    one_l = get_one_line(file, one_line=lInit)
    one_l = replace(one_l, ",", " ")
    one_l = replace(one_l, ";", "")

    sl = split(one_l, " ")
    var = ""
    i = 1
    while i < length(sl)
        if sl[i] == "var"
            i += 1
            var = string(sl[i])
            push!(mod.cols, var)
        elseif sl[i] == "integer"
            @assert !isempty(var)
            mod.colsType[var] = "Integer"
        elseif sl[i] == "binary"
            @assert !isempty(var)
            mod.colsType[var] = "Binary"
        elseif sl[i] == ":="
            @assert !isempty(var)
            i += 1
            mod.l[var] = parse(sl[i])
        elseif sl[i] == ">="
            @assert !isempty(var)
            i += 1
            mod.lb[var] = parse(sl[i])
        elseif sl[i] == "<="
            @assert !isempty(var)
            i += 1
            mod.ub[var] = parse(sl[i])
        elseif isempty(sl[i])
            non = 0
        else
            error("Detected exception in variable line. $(i) - $(sl)")
        end
        i += 1
    end

    return
end

function mod_read_equation(file::IOStream, mod::oneProblem, lInit::AbstractString)

    #=
        Example Equation
        #Single Line
        e1:    objvar - x145 - x146 - x147 - x148 - x149 - x150 = 3271.22725820856;
        # Multi Line
        e38:    109.0495*x1 + 100.462*x2 + 115.2937*x3 + 2.860271*x4 + 15.1404*x5
               + x44 - x556 >= 101.1304;
    =#

	eN = string(strip(split(lInit, r" |,")[1], [':', '.', '\r', '\n']))
    eS = split(eN, '.')

    push!(mod.rows, "$(eN)")

    one_l = get_one_line(file, one_line=lInit)
    one_l = replace(one_l, "$(eN):", "")

    if contains(one_l, "<=")
        mod.rowsSense[eN] = "L"
    elseif contains(one_l, ">=")
        mod.rowsSense[eN] = "G"
    elseif contains(one_l, "=")
        mod.rowsSense[eN] = "E"
    else
        error("NO sense detected in equation $(eN). $(one_l)")
    end

    symbolexchange = Dict("E"=>"=", "L"=>"<=", "G"=>">=")
    sl = split(one_l, "$(symbolexchange[mod.rowsSense[eN]])", keep=false)
    mod.rowsLHS[eN] = replace(sl[1], " ", "")
    mod.rowsRHS[eN] = Float64(parse(sl[2]))

    return
end

function mod_read_objective(file::IOStream, mod::oneProblem, lInit::AbstractString, sense::AbstractString)

    sense == "maximize" ? mod.objSense = "maximizing" :  mod.objSense = "minimizing"

    one_l = get_one_line(file, one_line=lInit)
    one_l = replace(one_l, "minimize", "")
    one_l = replace(one_l, "obj", "")
    one_l = replace(one_l, ":", "")
    one_l = replace(one_l, " ", "")

    mod.objective = one_l

    return
end

function mod2jump(modPath::AbstractString;
                    mode::AbstractString="index",
                    ending="m=m",
                    quadNL::Bool=false,
                    loopifpossible=true,
                    outdir::AbstractString="")

    probName = replace(splitdir(modPath)[end], ".mod", "")
    mod = read_mod_file(modPath)
    write_julia_script(probName, mod, mode, quadNL=quadNL, outdir=outdir)

    return
end


text:
"""imaginary time evolution solution to maximum cut problem."""
using Yao

@const_gate ZZ = mat(kron(Z,Z))

function maxcut_circuit(W::AbstractMatrix, τ)
    nbit = size(W, 1)
    ab = chain(nbit)
    for i=1:nbit,j=i+1:nbit
        if W[i,j] != 0
            #push!(ab, put(nbit, (i,j)=>0.5*W[i,j]*ZZ))
            push!(ab, put(nbit, (i,j)=>rot(ZZ, -im*W[i,j]*τ)))
        end
    end
	return ab
end

function maxcut_h(W::AbstractMatrix)
    nbit = size(W, 1)
    ab = Add{nbit}()
    for i=1:nbit,j=i+1:nbit
        if W[i,j] != 0
            #push!(ab, put(nbit, (i,j)=>0.5*W[i,j]*ZZ))
            push!(ab, put(nbit, (i,j)=>0.5*W[i,j]*ZZ))
        end
    end
	return ab
end

using Random, Test
@testset "iqaoa circuit" begin
	Random.seed!(2)
	nbit = 5
	# exact solution is [(1,3,4), (2,5)]
	W = [0 5 2 1 0;
		 5 0 3 2 0;
		 2 3 0 0 0;
		 1 2 0 0 4;
		 0 0 0 4 0]

	c = maxcut_circuit(W, 10)
	h = maxcut_h(W)
	reg = uniform_state(nbit) |> c

	# check the correctness of result
	opt_pl = reg |> probs
	config = argmax(opt_pl)-1
	@show BitStr64{nbit}(config), opt_pl[config+1]
	@test config == bit"01101" || config == bit"10010"
	@test isapprox(expect(h, reg |> normalize!), -5.5)
end