text:
@testset "Simulator: dynamics" begin
    dynamics_model = Simulator(halfcheetah, 1, h=0.02)
    dynamics_model.ip.opts.r_tol = 1e-7
    dynamics_model.ip.opts.κ_tol = 1e-5
    dynamics_model.ip.opts.undercut = 10.0
    nq = dynamics_model.model.nq
    nx = 2nq
    nu = dynamics_model.model.nu
    nw = dynamics_model.model.nw

    y = zeros(nx)
    dx = zeros(nx,nx)
    du = zeros(nx,nu)
    x = RoboDojo.nominal_state(dynamics_model.model) + 0.1*ones(nx)
    u = ones(nu)
    w = ones(nw)

    status = RoboDojo.dynamics(dynamics_model, y, x, u, w)
    @test status
    status = RoboDojo.dynamics_jacobian_state(dynamics_model, dx, x, u, w)
    @test status
    status = RoboDojo.dynamics_jacobian_input(dynamics_model, du, x, u, w)
    @test status

    function out_of_place_explicit_dynamics(dynamics_model, x, u, w)
        nq = dynamics_model.model.nq
        nx = 2nq
        y = zeros(nx)
        RoboDojo.dynamics(dynamics_model, y, x, u, w)
        return y
    end

    dx0 = FiniteDiff.finite_difference_jacobian(x -> out_of_place_explicit_dynamics(dynamics_model, x, u, w), x)
    du0 = FiniteDiff.finite_difference_jacobian(u -> out_of_place_explicit_dynamics(dynamics_model, x, u, w), u)
    @test norm(dx - dx0, Inf) < 1e-4
    @test norm(du - du0, Inf) < 1e-4
end


text:
using IntervalWavelets
using Test

include("DyadicRationals.jl")
include("Filters.jl")
include("InteriorScalingFunction.jl")
include("BoundaryScalingFunction.jl")
include("Basis.jl")

ext: jl


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    BisectPointPartition(normal, point)

A method for partitioning spatial objects into two half spaces
defined by a `normal` direction and a reference `point`.
"""
struct BisectPointPartition{Dim,T} <: PartitionMethod
  normal::Vec{Dim,T}
  point::Point{Dim,T}

  function BisectPointPartition{Dim,T}(normal, point) where {Dim,T}
    new(normalize(normal), point)
  end
end

BisectPointPartition(normal::Vec{Dim,T}, point::Point{Dim,T}) where {Dim,T} =
  BisectPointPartition{Dim,T}(normal, point)

BisectPointPartition(normal::NTuple{Dim,T}, point::NTuple{Dim,T}) where {Dim,T} =
 BisectPointPartition(Vec(normal), Point(point))

function partition(object, method::BisectPointPartition)
  n = method.normal
  p = method.point

  left, right = Int[], Int[]
  for location in 1:nelements(object)
    pₒ = centroid(object, location)
    if (pₒ - p) ⋅ n < zero(coordtype(object))
      push!(left, location)
    else
      push!(right, location)
    end
  end

  Partition(object, [left,right])
end


text:
using DelimitedFiles

data= readdlm("input_day2", '\n')

valid_pass = 0

for line in data

   i1 = findfirst(isequal('-'),line)
   i2 = findfirst(isequal(' '),line)

   #position which has to have letter
   pos_1 = parse(Int, line[1:i1-1])
   #position which has to not have the letter
   pos_2 = parse(Int, line[i1+1:i2-1])

   letter = line[i2+1]
   password = line[i2+4:end]

   if password[pos_1] == letter
      if password[pos_2] != letter
         global valid_pass += 1
      end
   else
      if password[pos_2] == letter
         global valid_pass += 1
      end
   end

end

println(valid_pass)


text:
# ============================================================================ #
# The SimplicialComplex type
# ============================================================================ #
"""
  SimplicialComplex{T}

The simplicial complex. Holds data of type `T` (not implemented.)

Supported operations:

* push!(sxcx, σ)   [√]
* σ in sxcx        [√]
* delete!(sxcx, σ) [Χ]

"""
mutable struct SimplicialComplex{T}
    nodelist::Vector{DictSet{Int, Int}}

    # TODO: store data
    # reverse dict for maximal simplices?
    lbl_values::IntSet
    lbl_holes::IntSet
end

SimplicialComplex{T}() where T =
    SimplicialComplex{T}(DictSet{Int, Int}[], IntSet(), IntSet())
SimplicialComplex() =
    SimplicialComplex{Float64}()

# ---------------------------------------------------------------------------- #
# Base overloads.

Base.show(io::IO, sxcx::SimplicialComplex) =
    print(io, "SimplicialComplex with $(length(sxcx)) nodes")

# Pretty printing.
function Base.print(io::IO, sxcx::SimplicialComplex)
    println(io, "SimplicialComplex:")
    for i in eachindex(sxcx)
        print(io, "$(lpad(i, 5)): [")
        for k in sort(collect(keys(sxcx.nodelist[i].data)))
            for v in sort(collect(sxcx.nodelist[i].data[k]))
                print(io, "($i, $(iszero(k) ? "φ" : k), $v)")
            end
        end
        println(io, "]")
    end
end

Base.length(sxcx::SimplicialComplex) = length(sxcx.nodelist)
Base.eachindex(sxcx::SimplicialComplex) = 1:length(sxcx.nodelist)
Base.isempty(sxcx::SimplicialComplex) = all(isempty.(sxcx.nodelist))

# ============================================================================ #
# Operation helpers.
# ============================================================================ #
"""
  newlabel!(sxcx)

Create a new maximal simplex label and return it.
"""
function newlabel!(sxcx::SimplicialComplex)
    if isempty(sxcx.lbl_holes)
        lbl = length(sxcx.lbl_values) + 1
    else
        lbl = first(sxcx.lbl_holes)
        delete!(sxcx.lbl_holes, lbl)
    end
    push!(sxcx.lbl_values, lbl)
    lbl
end

"""
  insnode!(sxcx, cur, nxt, lbl)

Insert a node that represents a link between `cur` and `nxt` in the maximal
simplex with label `lbl`.
"""
function insnode!(sxcx::SimplicialComplex, cur::Int, nxt::Int, lbl::Int)
    while length(sxcx) < cur
        push!(sxcx.nodelist, DictSet{Int, Int}())
    end

    insert!(sxcx.nodelist[cur], nxt, lbl)
end

"""
  labels(sxcx, from, to)

Return the labels of maximal simplices that contain the link between `from` and
`to`.
"""
function labels(sxcx::SimplicialComplex, from::Int, to::Int)
    sxcx.nodelist[from][to]
end

"""
  alllabels(sxcx)

Return all labels of maximal simplices in `sxcx`.
"""
function alllabels(sxcx::SimplicialComplex)
    sxcx.lbl_values
end

"""
  maximalsimplices(sxcx)

Return all maximal simplices in `sxcx` in a Dict{Int, Vector{Int}}.

TODO: This is slow, make an iterator or avoid doing this entirely.
"""
function maximalsimplices(sxcx::SimplicialComplex)
    maxsxs = Dict{Int, Vector{Int}}()
    for l in sxcx.lbl_values
        maxsxs[l] = []
    end
    for i in eachindex(sxcx)
        for v in values(sxcx.nodelist[i])
            push!(maxsxs[v], i)
        end
    end
    maxsxs
end

# ============================================================================ #
# Operations
# ============================================================================ #
# A σ is in the complex if it is contained in one of the maximal sx's.
function Base.in(σ, sxcx::SimplicialComplex)
    σ = validatesx(σ)
    !isempty(maximalcofaces(sxcx, σ))
end

"""
  maximalcofaces(sxcx, σ)

List the labels of maximal simplices that are cofaces of `σ`.
"""
function maximalcofaces(sxcx::SimplicialComplex, σ)
    last(σ) > length(sxcx) && return Set{Int}()

    s = labels(sxcx, σ[1], σ[2])
    for l in zip(σ[2:end-1], σ[3:end])
        s = intersect(s, labels(sxcx, l...))
        isempty(s) && return s
    end
    s
end

"""
  maximalfaces(sxcx, σ)

List the labels of maximal simplices that are faces of `σ`.
"""
function maximalfaces(sxcx::SimplicialComplex, σ)
    # Find labels of maximal simplices that have a node that is not in σ and
    # setdiff them from all labels.
    s = Set{Int}()
    for i in setdiff(eachindex(sxcx), σ)
        union!(s, values(sxcx.nodelist[i]))
    end

    maxsxs = Set(sxcx.lbl_values)
    setdiff(maxsxs, s)
end

"""
  deletemaximals!(sxcx, lbls)

Delete maximal simplices with labels in lbls.
"""
function deletemaximals!(sxcx::SimplicialComplex, lbls::Set{Int})
    for i in eachindex(sxcx)
        delete!(sxcx.nodelist[i], lbls)
    end
    union!(sxcx.lbl_holes, lbls)
    setdiff!(sxcx.lbl_values, lbls)
end

# DONE.
function Base.push!(sxcx::SimplicialComplex, σ)
    σ = validatesx(σ)
    σ in sxcx && return σ

    # Remove maximal simplices that are faces of σ.
    deletemaximals!(sxcx, maximalfaces(sxcx, σ))

    lbl = newlabel!(sxcx)
    for (i, li) in enumerate(σ), (j, lj) in enumerate(σ[i+1:end])
        insnode!(sxcx, li, lj, lbl)
    end
    insnode!(sxcx, σ[end], 0, lbl)
    σ
end

"""

"""
function facetswithout(σ, τ)
end

function Base.delete!(sxcx::SimplicialComplex, σ)
    maxlbls = maximalcofaces(sxcx, σ)
    # make dim(σ) copies of each of the maximal simplices and in the i-th copy
    # remove nodes with label (σ[i], x, ℓ)

    # Torej v bistvu iz vsakega daš eno točko
end


text:
module NumericalEFT
using StaticArrays
include("common.jl")
export σx, σy, σz, σ0

include("fastmath.jl")
export FastMath

include("utility/utility.jl")
export Utility

##########################  Compact Lehmann representation  ########################################
include("Lehmann/Lehmann.jl")
using .Lehmann
# main module
export Lehmann
# submodules
export Spectral, Sample
export DLRGrid
export tau2dlr, dlr2tau, matfreq2dlr, dlr2matfreq, tau2matfreq, matfreq2tau, tau2tau, matfreq2matfreq

##########################  Monte Carlo Integration  ###############################################
include("MCIntegration/MCIntegration.jl")
using .MCIntegration
# main module
export MCIntegration
export montecarlo, Configuration, Continuous, Discrete

##########################  CompositeGrids  #######################################################
include("CompositeGrids/CompositeGrids.jl")
using .CompositeGrids
# main module
export CompositeGrids
# submodules
export Grid, Interp
export SimpleG, SimpleGrid, AbstractGrid, OpenGrid, ClosedGrid, denseindex
export CompositeG, CompositeGrid

##########################  Green's function  #######################################################
include("GreenFunc/GreenFunc.jl")
using .GreenFunc
# main module
export GreenFunc
export TimeDomain, ImTime, ReTime, ImFreq, ReFreq, DLRFreq
export Green2DLR, toTau, toMatFreq, toDLR, dynamic, instant

##########################  Diagram Builder  #######################################################
include("FeynmanDiagram/FeynmanDiagram.jl")
using .FeynmanDiagram
# main module
export FeynmanDiagram
export SigmaDiag, PolarDiag, Ver3Diag, Ver4Diag, GreenDiag
export Wirreducible, Girreducible, NoBubble, NoHatree, NoFock, Proper
export Response, ChargeCharge, SpinSpin, UpUp, UpDown
export AnalyticProperty, Instant, Dynamic
export GenericPara, Interaction

# submodules
export uidreset, toDataFrame, mergeby, plot_tree
export DiagTree
export TwoBodyChannel, Alli, PHr, PHEr, PPr, AnyChan
export Permutation, Di, Ex, DiEx
export Diagram, addSubDiagram!, toDataFrame
export evalDiagNode!, evalDiagTree!
export Operator, Sum, Prod
export DiagramId, GenericId, Ver4Id, Ver3Id, GreenId, SigmaId, PolarId
export PropagatorId, BareGreenId, BareInteractionId
export BareGreenNId, BareHoppingId, GreenNId, ConnectedGreenNId
export uidreset, toDataFrame, mergeby, plot_tree

# submodules
export Parquet
export ParquetBlocks

# submodules
export ExprTree
export Component, Diagrams
export addpropagator!, addnode!
export setroot!, addroot!
export evalNaive, showTree

include("Atom/Atom.jl")
using .Atom
export Atom
# export Hilbert

end # module


text:
#=
Main generic API

@author : Spencer Lyon <spencer.lyon@stern.nyu.edu>
@date : 2015-04-13 11:29:12

=#

# ----------- #
# Manager API #
# ----------- #
"""
`verbose(m::IterationManager) -> Bool`

Specifies if the iteration manager should be verbose and print status updates.

Checks the `verbose` field of the manager and returns it. If the field
doesn't exist, the default is false.

"""
verbose(m::IterationManager) = isdefined(m, :verbose) ? m.verbose : false

"""
`prefix(m::IterationManager) -> String`

Defines the prefix that should be applied to all printed messages

Checks the `prefix` field of the manager and returns it. If the field
doesn't exist, the default is `""`.

"""
prefix(m::IterationManager) = isdefined(m, :print_prefix) ? m.print_prefix : ""

"""
`print_now(m::IterationManager, n::Int) -> Bool`

Specifies if the iteration manager should print on the `n`th iteration
"""
function print_now(mgr::IterationManager, n::Int)
    if verbose(mgr)
        n == 1 && return true
        isdefined(mgr, :print_skip) ? n % mgr.print_skip == 0 : true
    else
        false
    end
end

# fallback in case x is not an Int
print_now(mgr::IterationManager, x::Any) = verbose(mgr)

# --------- #
# State API #
# --------- #

# default methods for the `by` argument of `update!`
default_by(x::T, y::T) where T <: AbstractArray = maximum(abs.(x .- y))
default_by(x::T, y::T) where T <: Number = abs(x - y)

function default_by(x::S, y::T) where S where T
    msg = "default_by not implemented for types x::$S, y::$T"
    throw(ArgumentError(msg))
end

num_iter(::IterationState) = nothing
LinearAlgebra.norm(::IterationState) = Inf

"""
`default_by(x, y)`

Gives the default comparison `by` argument to `managed_iteration` based on the
types of x and y
""" default_by

display_iter(istate::IterationState, prefix="") =
    display_iter(stdout, istate, prefix)

display_iter(io::IO, istate::T, prefix="") where T <: IterationState =
    error("`display_iter` must be implemented directly by type $T")

"""
    display_iter([io::IO=stdout], st::IterationState, [prefix::String])

A `display` method for an iteration state. Must be implemented by all concrete
subtypes of `IterationState`
"""
display_iter

# --------------------------- #
# Combining Manager and State #
# --------------------------- #

"""
    finished(mgr::IterationManager, istate::IterationState)::Bool

Given a manager and a state, determine if the iterations have finished and
should thus be terminated
"""
finished

# ------------- #
# default hooks #
# ------------- #
"""
    pre_hook(mgr::IterationManager, istate::IterationState)

Called before iterations begin.

For the default manager and state this is used to print the column headers for
verbose output printing
"""
pre_hook(mgr::IterationManager, istate::IterationState) =
    verbose(mgr) && display_iter(istate, prefix(mgr))

print_now(mgr::IterationManager, istate::IterationState) =
    print_now(mgr, num_iter(istate))

"""
    iter_hook(mgr::IterationManager, istate::IterationState)

Called after every iteration.

For the default manager and state this is used to print updates on the
iterations if `verbose(mgr) == true`.
"""
iter_hook(mgr::IterationManager, istate::IterationState) =
    print_now(mgr, istate) && display_iter(istate, prefix(mgr))

"""
    post_hook(mgr::IterationManager, istate::IterationState)

Called after iterations have finished

For the default manager this is used to print a warning if the maximum number
of iterations was exceeded.
"""
function post_hook(mgr::IterationManager, istate::IterationState)
    if !(isdefined(mgr, :maxiter))
        return nothing
    end
    if num_iter(istate) >= mgr.maxiter
        @warn "Maximum iterations exceeded. Algorithm may not have converged"
    end
    nothing
end

function managed_iteration(
        f::Base.Callable,
        mgr::IterationManager,
        istate::IterationState{T};
        by=default_by
    ) where T
    pre_hook(mgr, istate)

    while !(finished(mgr, istate))
        v = f(istate.prev)::T  # if we don't get a T back it should be an error
        update!(istate, v; by=by)
        iter_hook(mgr, istate)
    end

    post_hook(mgr, istate)
    istate
end

"""
    managed_iteration!(
        f!::Base.Callable,
        mgr::IterationManager,
        dest::T,
        istate::IterationState{T};
        by::Base.Callable=default_by
    ) where T <: AbstractArray

Given a function with signature `f!(dest::T, src::T)`, run a non-allocating
version of managed_iteration where each step calls the mutating function `f!`
that fills its first argument based on the value of the second argument. All
other behavior is equivalent to other forms of `managed_iteration`

"""
function managed_iteration!(
        f!::Base.Callable,
        mgr::IterationManager,
        dest::T,
        istate::IterationState{T};
        by::Base.Callable=default_by
    ) where T <: AbstractArray
    pre_hook(mgr, istate)

    while !(finished(mgr, istate))
        f!(dest, istate.prev)
        update!(istate, dest; by=by)
        iter_hook(mgr, istate)
    end

    post_hook(mgr, istate)
    istate
end

# kwarg version to create default manger/state
function managed_iteration(
        f::Base.Callable,
        init;
        tol::Float64=NaN,
        maxiter::Int=typemax(Int),
        by::Base.Callable=default_by,
        verbose::Bool=true,
        print_skip::Int=div(maxiter, 5)
    )
    mgr = DefaultManager(tol, maxiter, verbose, print_skip)
    istate = DefaultState(init)
    managed_iteration(f, mgr, istate; by=by)
end


text:
using RecurrenceAnalysis
using DynamicalSystemsBase, Random, Statistics
using Test

RA = RecurrenceAnalysis

rng = Random.seed!(194)
# Trajectories of 200 points
# Examples of the Hénon map based on:
#     Webber CL & Zbilut JP, "Recurrence Quantification Analysis of Nonlinear
#     Dynamical Systems", in: Riley MA & Van Orden GC, Tutorials in Contemporary
#     Nonlinear Methods for the Behavioral Sciences, 2005, 26-94.
#     https://www.nsf.gov/pubs/2005/nsf05057/nmbs/nmbs.pdf
#
trajectories = Dict(
    "Sine wave" => RA.Dataset(map(x->[sin.(x) cos.(x)], StepRangeLen(0.0,0.2,200))),
    "White noise" => RA.Dataset(randn!(zeros(200,2))),
    "Hénon (chaotic)" => trajectory(Systems.henon(a=1.4, b=0.3), 199, Ttr=1000),
    "Hénon (periodic)" => trajectory(Systems.henon(a=1.054, b=0.3), 199, Ttr=1000)
    )
embed_params = Dict(       #(m, τ)
    "Sine wave"   => (9, 7),
    "White noise" => (1, 1),
    "Hénon (chaotic)" => (3, 1),
    "Hénon (periodic)" => (3, 1))
rqa_threshold = Dict(
    "Sine wave"   => 0.15,
    "White noise" => 0.15,
    "Hénon (chaotic)" => 0.15,
    "Hénon (periodic)" => 0.15)

dict_keys = ["Sine wave","White noise","Hénon (chaotic)","Hénon (periodic)"]
@testset "$k" for k in dict_keys
    data = trajectories[k]
    x = data[:,1]
    y = data[1:100,2]
    xe = embed(x, embed_params[k]...)
    ye = embed(y, embed_params[k]...)

    # Distance and recurrence matrices
    ε = rqa_threshold[k]
    dmat = distancematrix(xe, ye)
    crmat = CrossRecurrenceMatrix(xe, ye, ε)
    @test Matrix(crmat.data) == (dmat .≤ ε)
    rmat = RecurrenceMatrix(xe, ε)
    jrmat = JointRecurrenceMatrix(xe, ye, ε)
    sz = size(jrmat)
    @test (jrmat.data .& rmat.data[1:sz[1],1:sz[2]]) == jrmat.data
    # Compare metrics
    m_euc = rmat
    m_max = RecurrenceMatrix(xe, ε, metric="max")
    m_min = RecurrenceMatrix(xe, ε, metric="manhattan")
    @test (m_max.data .& m_euc.data) == m_euc.data
    @test (m_euc.data .& m_min.data) == m_min.data
    # Compare scales
    m_scalemax = CrossRecurrenceMatrix(xe, ye, 0.1, scale=maximum)
    m_scalemean = CrossRecurrenceMatrix(xe, ye, 0.1, scale=mean)
    @test (m_scalemax.data .& m_scalemean.data) == m_scalemean.data
    # Fixed rate for recurrence matrix
    crmat_fixed = CrossRecurrenceMatrix(xe, ye, 0.05; fixedrate=true)
    @test .04 < recurrencerate(crmat_fixed) < .06

    # Recurrence plot
    crp = grayscale(crmat, width=125)
    szplot = size(crp)
    szmat  = size(crmat)
    @test szplot[1] ≈ szplot[2]*szmat[2]/szmat[1] atol = 1

    # RQA
    rqapar = rqa(rmat, theiler=1, lmin=3, border=20)
    rqadiag = rqa(rmat, theiler=1, lmin=3, border=20, onlydiagonal=true)
    for p in keys(rqadiag)
        @test rqapar[p]==rqadiag[p]
    end

    # Windowed RQA
    rmatw = @windowed RecurrenceMatrix(xe, ε, metric=RecurrenceAnalysis.Chebyshev()) 50
    crmatw = @windowed(CrossRecurrenceMatrix(xe, ye, ε),30)
    @windowed jrmatw = JointRecurrenceMatrix(xe, ye, ε) 30
    @test jrmatw[3 .+ (1:30), 3 .+ (1:30)] == jrmat[3 .+ (1:30), 3 .+ (1:30)]
    @windowed(rrw = recurrencerate(rmatw), width=50, step=40)
    @windowed rqaw = rqa(rmatw) width=50 step=40
    @test rqaw[:RR] == rrw

end


text:
"""
    mcmc(gp::GPBase; kwargs...)

Run MCMC algorithms provided by the Klara package for estimating the hyperparameters of
Gaussian process `gp`.
"""
function mcmc(gp::GPBase; nIter::Int=1000, burn::Int=1, thin::Int=1, ε::Float64=0.1,
              Lmin::Int=5, Lmax::Int=15, lik::Bool=true, noise::Bool=true,
              domean::Bool=true, kern::Bool=true)
    Kgrad = Array{Float64}(undef, gp.nobs, gp.nobs)
    L_bar = Array{Float64}(undef, gp.nobs, gp.nobs)
    params_kwargs = get_params_kwargs(gp; domean=domean, kern=kern, noise=noise, lik=lik)
    count = 0
    function calc_target(gp::GPBase, θ::VecF64) #log-target and its gradient
        count += 1
        try
            set_params!(gp, θ; params_kwargs...)
            update_target_and_dtarget!(gp, Kgrad, L_bar; params_kwargs...)
            return true
        catch err
            if !all(isfinite.(θ))
                return false
            elseif isa(err, ArgumentError)
                return false
            elseif isa(err, LinearAlgebra.PosDefException)
                return false
            else
                throw(err)
            end
        end
    end

    θ_cur = get_params(gp; params_kwargs...)
    D = length(θ_cur)
    leapSteps = 0                   #accumulator to track number of leap-frog steps
    post = Array{Float64}(undef, nIter, D)     #posterior samples
    post[1,:] = θ_cur

    @assert calc_target(gp, θ_cur)
    target_cur, grad_cur = gp.target, gp.dtarget

    num_acceptances = 0
    for t in 1:nIter
        θ, target, grad = θ_cur, target_cur, grad_cur

        ν_cur = randn(D)
        ν = ν_cur + 0.5 * ε * grad

        reject = false
        L = rand(Lmin:Lmax)
        leapSteps +=L
        for l in 1:L
            θ += ε * ν
            if  !calc_target(gp,θ)
                reject=true
                break
            end
            target, grad = gp.target, gp.dtarget
            ν += ε * grad
        end
        ν -= 0.5*ε * grad

        if reject
            post[t,:] = θ_cur
        else
            α = target - 0.5 * ν'ν - target_cur + 0.5 * ν_cur'ν_cur
            u = log(rand())

            if u < α
                num_acceptances += 1
                θ_cur = θ
                target_cur = target
                grad_cur = grad
            end
            post[t,:] = θ_cur
        end
    end
    post = post[burn:thin:end,:]
    set_params!(gp, θ_cur; params_kwargs...)
    @printf("Number of iterations = %d, Thinning = %d, Burn-in = %d \n", nIter,thin,burn)
    @printf("Step size = %f, Average number of leapfrog steps = %f \n", ε,leapSteps/nIter)
    println("Number of function calls: ", count)
    @printf("Acceptance rate: %f \n", num_acceptances/nIter)
    return post'
end


text:
using ForwardDiff
using FixedPointAcceleration
func(x, params) = [params[1]*sqrt(abs(x[1] + x[2])), 1.5*x[1] + params[2]*x[2]]

function ML_estimation(params)
    Inputs = [0.3,900.0]
    fp_simple   = fixed_point(inp -> func(inp, params), Inputs; Algorithm = :Simple)
    return fp_simple.FixedPoint_[1]
end

params = [1.5, 0.5]
ML_estimation(params)
grads = ForwardDiff.gradient(ML_estimation, params)
isa(grads, Array{Float64,1})


text:
# Adapted from:
# SOSDEMO1 --- Sum of Squares Test
# Section 3.1 of SOSTOOLS User's Manual

@testset "SOSDEMO1 with $(factory.optimizer_constructor)" for factory in sdp_factories
    @polyvar x y

    m = SOSModel(factory)

    # Example 2.4 of
    # P. Parrilo and A. Jadbabaie
    # "Approximation of the joint spectral radius using sum of squares."
    # Linear Algebra and its Applications, Elsevier, 2008, 428, 2385-2402

    p = 2*x^4 + 2*x^3*y - x^2*y^2 + 5*y^4

    soscon = @constraint m p >= 0

    JuMP.optimize!(m)

    @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT

    q = gram_matrix(soscon)
    Q = getmat(q)
    @test issymmetric(Q)
    @test isapprox(Q[1, 1], 2, rtol=1e-5)
    @test isapprox(Q[1, 2], 1, rtol=1e-5)
    @test isapprox(Q[3, 3], 5, rtol=1e-5)
    @test abs(Q[2, 3]) < 1e-5
    @test isapprox(Q[2, 2] + 2Q[1, 3], -1, rtol=1e-5)
    sosdec = SOSDecomposition(q)
    @test isapprox(sosdec, sos_decomposition(soscon))
    @test isapprox(sum(sosdec.ps.^2), p; rtol=1e-4, ztol=1e-6)

    M = SOSModel(factory)

    p = 4*x^4*y^6 + x^2 - x*y^2 + y^2

    soscon = @constraint M p >= 0

    JuMP.optimize!(M)

    @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT

    # p should be GramMatrix([1, 0, -1/2, 0, -1, 1, 0, -2/3, 0, 4/3, 0, 0, 2, 0, 4], [y, x, x*y, x*y^2, x^2*y^3])
end


text:
struct Simulator{T}
    model::AbstractModel{T}
    z::Vector{T}
    x::Vector{T}
end

struct SimulatorState{T}
    z::Vector{T}
    a::Vector{T}
end

function simulate(model::AbstractModel{T},  z::AbstractVector{T},  x::AbstractVector{T}) where T
    return Simulator(model,  z,  x)
end

simulate(model::AbstractModel{T},  z::AbstractVector{T}) where T = Simulator(model,  z,  T[])
simulate(model::AbstractModel{T}) where T = Simulator(model, T[], T[])

function Base.iterate(simulator::Simulator{T}) where T
    α = residuals(simulator.model, simulator.z)
    z = forecast(simulator.model, simulator.z, α)
    a = √simulator.model.σ2 * randn(T)
    push!(simulator.z, z + a)
    push!(α, a)
    return (z + a, SimulatorState(simulator.z, α))
end

function Base.iterate(simulator::Simulator{T}, state::SimulatorState{T}) where T
    z = forecast(simulator.model, state.z, state.a)
    a = √simulator.model.σ2 * randn(T)
    push!(state.z, z + a)
    push!(state.a, a)
    return (z + a, state)
end


text:
using DrWatson
@quickactivate "AccuracyAtTop_DeepTopPush"

include(srcdir("utilities.jl"))

# ------------------------------------------------------------------------------------------
# Basic datasets
# ------------------------------------------------------------------------------------------
Dataset_Settings = Dict(
    :dataset => [FashionMNIST, CIFAR100, SVHN2Full, Molecules],
    :posclass => [0, 1, 1, 1],
)

Train_Settings = Dict(
    :batchsize => 32,
    :epochs => 200,
    :optimiser => Descent,
    :steplength => 0.0001,
    :seed => collect(1:10),
)

Model_Settings = Dict(
    :type => [BaseLine, DeepTopPush, PatMatNP],
    :arg => [missing, missing, 0.01],
    :surrogate => [missing, quadratic, quadratic],
    :reg => 1e-3,
    :buffer => [missing, true, missing],
)

run_simulations(Dataset_Settings, Train_Settings, Model_Settings; nepochs_save = 50)
run_evaluation(Dataset_Settings)

# ------------------------------------------------------------------------------------------
# APPerf
# ------------------------------------------------------------------------------------------
Dataset_Settings = Dict(
    :dataset => [FashionMNIST, CIFAR100, Molecules],
    :posclass => [0, 1, 1],
)

Train_Settings = Dict(
    :batchsize => 32,
    :epochs => 200,
    :optimiser => Descent,
    :steplength => 0.0001,
    :seed => 1,
)

Model_Settings = Dict(
    :type => [APPerf],
    :arg => [0.01],
    :surrogate => missing,
    :reg => 1e-3,
    :buffer => missing,
)

run_simulations(Dataset_Settings, Train_Settings, Model_Settings)
run_evaluation(Dataset_Settings)

# ------------------------------------------------------------------------------------------
# Imagenet
# ------------------------------------------------------------------------------------------
Dataset_Settings = Dict(
    :dataset => [ImageNet],
    :posclass => [33:37], # turtles
)

Train_Settings = Dict(
    :batchsize => 32,
    :epochs => 10,
    :optimiser => Descent,
    :steplength => 0.0001,
    :seed => 1,
)

Model_Settings = Dict(
    :type => [DeepTopPush, PatMatNP, PatMatNP],
    :arg => [missing, 0.01, 0.05],
    :surrogate => [hinge, hinge, hinge],
    :reg => 1e-3,
    :buffer => [false, missing, missing],
)

run_simulations(Dataset_Settings, Train_Settings, Model_Settings; nepochs_save = 10)
run_evaluation(Dataset_Settings)

@info "All experiments finished"


text:
using SmartShockFinder

datatyp = "NNPPh"
data_input_size = 3
n_meshes = 4
polydeg = collect(2:9)

generate_validdataset1d(datatyp, n_meshes, polydeg, data_input_size)


text:
#!/usr/bin/env julia

using ArgParse

s = ArgParseSettings()
@add_arg_table s begin
    "-M"
    help = "Print out the list of files this program will create."
    action = :store_true
    "config"
    help = "a YAML configuration file"
    default = "config.yaml"
end

parsed_args = parse_args(ARGS, s)

import YAML
config = YAML.load(open(parsed_args["config"]))

using BenchmarkTools
using DiskTracer.model
using DiskTracer.constants
using DiskTracer.geometry
using DiskTracer.trace

species = config["species"]
transition = config["transition"]
lam0 = lam0s[species*transition]
model = config["model"]
pars = convert_dict(config["parameters"], config["model"])

mol = mols[species*transition]

nr=128
rmax = 500 * AU
nz=128
zmax = 200 * AU

# Get the interpolator object for this particular setup.
interp = DiskTracer.model.get_grids(pars, mol, mol.nu_0, nr, rmax, nz, zmax)

# pick a pixel (xprime, yprime), and trace a full ray through the disk

v0 = 1.0 # km/s

I = trace_pixel(100*AU, -100*AU, v0, mol, pars, interp)

# @btime trace_pixel(100*AU, -100*AU, v0, mol, pars, interp)
println(I)

quit()
# @btime trace_pixel(100*AU, 80*AU, v0, mol, pars, grids)
# @profile trace_pixel(100*AU, 80*AU, v0, mol, pars, grids)
# Profile.print()

rmax = 600*AU #
xprime = 100*AU
yprime = 80*AU

rcyl_min = abs(xprime)
DV2_max, junk, junk2 = interp(rcyl_min, 0.0)
DeltaVmax = sqrt(DV2_max) * 1e-5 # km/s

z1start, z1end = geometry.get_bounding_zps(xprime, yprime, pars, v0, DeltaVmax, rmax)
gpre = precalc_geo(xprime, yprime, pars)
args = (z1start, v0, pars, interp, gpre)

zprime = 2*AU

rcyl = DiskTracer.geometry.get_r_cyl(xprime, yprime, zprime, pars.incl)
z = DiskTracer.geometry.get_z(xprime, yprime, zprime, pars.incl)

# Try evaluating f
println("Get vlos")
@btime DiskTracer.geometry.get_vlos(xprime, yprime, zprime, pars)

println("Get r_cyl")
@btime DiskTracer.geometry.get_r_cyl(xprime, yprime, zprime, pars.incl)

println("Get z")
@btime DiskTracer.geometry.get_z(xprime, yprime, zprime, pars.incl)

println("Interp")
@btime interp(rcyl, z)

println("f")
@btime DiskTracer.trace.f(Float64[0.0, 0.0], Float64[0.0, 0.0], args, 2*AU)

# DiskTracer.trace.f(Float64[0.0, 0.0], Float64[0.0, 0.0], args, 2*AU)
# # Profile.init(delay=1e-16)
# @profile [DiskTracer.trace.f(Float64[0.0, 0.0], Float64[0.0, 0.0], args, 2*AU) for i=1:100]
# Profile.print()


text:
#===============================================================================
        Additional utility functions, not tied to the package
===============================================================================#

# these appear in multiple places of this suite, TODO define them only once
ismutable(el) = ismutable(typeof(el))
ismutable(::Type) = Val(false)
ismutable(::Type{<:Array}) = Val(true)

"""
    _assure_scalar(v)

Accept a scalar or a vector of length one and return it as a scalar. Raise
assertion error if other objects are passed.
"""
function _assure_scalar(v)
    typeof(v) <: Number && return v
    @assert length(v) == 1
    @assert typeof(v[1]) <: Number
    first(v)
end

"""
    _upgrade_to_vec(v, N)

Receive a scalar or a vector of length 1 and return a length `N` vector of
repeats of entry v.
"""
function _upgrade_to_vec(v, N)
    typeof(v) <: Vector && length(v) == N && return v
    @assert length(v) == 1
    repeat([first(v)], N)
end

"""
    _upgrade_to_svec(v, ::Val{N}) where N

Receive a scalar or a vector of length 1 and return a length `N` static vector
of repeats of entry v.
"""
function _upgrade_to_svec(v, ::Val{N}) where N
    typeof(v) <: SVector && length(v) == N && return v
    length(v) == N && return SVector{N}(v)
    @assert length(v) == 1
    first(v) .+ zero(SVector{N})
end

"""
    sigmoid(x, a=1.0)

Sigmoid function. (Inverse of logit).
"""
sigmoid(x, a=1.0) = 1.0 / (1.0 + exp(-a*x))

"""
    logit(x, a=1.0)

Logit function. (Inverse of sigmoid).
"""
logit(x, a=1.0) = (log(x) - log(1-x))/a

"""
    custom_zero(x::T, ::Type{elT}) where {T,elT}

Create a `zero` with eltype `elT`, that is of the same structure as the
collection `x`.

# Examples
```julia-repl
julia> custom_zero(3.0, Bool)
false
julia> custom_zero([3.0, 4.0], Bool)
2-element Array{Bool,1}:
 0
 0
julia> custom_zero(SVector{2}(1, 2), ComplexF64)
2-element SArray{Tuple{2},Complex{Float64},1,2} with indices SOneTo(2):
 0.0 + 0.0im
 0.0 + 0.0im
```
"""
function custom_zero(x::T, ::Type{elT}) where {T,elT}
    T <: Array && return zeros(elT, size(x))
    T <: Number && return zero(elT)
    T <: SArray && return zero(similar_type(T, elT))
end


text:
using Pkg
using SafeTestsets

@time @safetestset "Plasma tests" begin include("plasmatests.jl") end


text:
g = smallgraph("diamondgraph")
@test nv(g) == 4 && ne(g) == 5

g = smallgraph("diamond")
@test nv(g) == 4 && ne(g) == 5

g = smallgraph(:diamond)
@test nv(g) == 4 && ne(g) == 5

g = smallgraph(:bull)
@test nv(g) == 5 && ne(g) == 5

g = smallgraph(:chvatal)
@test nv(g) == 12 && ne(g) == 24

g = smallgraph(:cubical)
@test nv(g) == 8 && ne(g) == 12

g = smallgraph(:desargues)
@test nv(g) == 20 && ne(g) == 30

g = smallgraph(:dodecahedral)
@test nv(g) == 20 && ne(g) == 30

g = smallgraph(:frucht)
@test nv(g) == 20 && ne(g) == 18

g = smallgraph(:heawood)
@test nv(g) == 14 && ne(g) == 21

g = smallgraph(:house)
@test nv(g) == 5 && ne(g) == 6

g = smallgraph(:housex)
@test nv(g) == 5 && ne(g) == 8

g = smallgraph(:icosahedral)
@test nv(g) == 12 && ne(g) == 30

g = smallgraph(:krackhardtkite)
@test nv(g) == 10 && ne(g) == 18

g = smallgraph(:moebiuskantor)
@test nv(g) == 16 && ne(g) == 24

g = smallgraph(:octahedral)
@test nv(g) == 6 && ne(g) == 12

g = smallgraph(:pappus)
@test nv(g) == 18 && ne(g) == 27

g = smallgraph(:petersen)
@test nv(g) == 10 && ne(g) == 15

g = smallgraph(:sedgewickmaze)
@test nv(g) == 8 && ne(g) == 10

g = smallgraph(:tetrahedral)
@test nv(g) == 4 && ne(g) == 6

g = smallgraph(:truncatedcube)
@test nv(g) == 24 && ne(g) == 36

g = smallgraph(:truncatedtetrahedron)
@test nv(g) == 12 && ne(g) == 18 && !is_directed(g)

g = smallgraph(:truncatedtetrahedron_dir)
@test nv(g) == 12 && ne(g) == 18 && is_directed(g)

g = smallgraph(:tutte)
@test nv(g) == 46 && ne(g) == 69

@test_throws ErrorException g = smallgraph(:nonexistent)


text:
import Weave
import Images
import Glob

weaving = true
# filepath = "ch0100.jl"
# img = Images.load("data/array2d.png")
# img = Images.imresize(img, ratio=0.15)
# Images.save("data/array2d_s.png", img)
#
# Weave.weave(filepath; doctype = "md2html", out_path = :pwd)
#
#
# filepath = "ch0100quiz.jl"
# height = 200
# img = Images.load("data/array_cartoon.png")
# ratio = height/size(img)[1]
# img = Images.imresize(img, ratio=ratio)
# Images.save("data/array_cartoon_s.png", img)
# img = Images.load("data/without_arrays.png")
# ratio = height/size(img)[1]
# img = Images.imresize(img, ratio=ratio)
# Images.save("data/without_arrays_s.png", img)
#
# Weave.weave(filepath; doctype = "md2html", out_path = :pwd)
#
#
# filepath = "ch0300.jl"
# Weave.weave(filepath; doctype = "md2html", out_path = :pwd)

# # Use modified tpl to fix a Plotly bug
# filepath = "ch0300quiz.jl"
# Weave.weave(filepath; doctype = "md2html", out_path = :pwd,
#     template="templates/julia_html.tpl",
#     css="templates/skeleton_css.css")

# filepath = "ch0500.jl"
# Weave.weave(filepath; doctype = "md2html", out_path = :pwd,
#     template="templates/julia_html.tpl",
#     css="templates/skeleton_css.css")

# filepath = "ch0500quiz.jl"
# Weave.weave(filepath; doctype = "md2html", out_path = :pwd,
#     template="templates/julia_html.tpl",
#     css="templates/skeleton_css.css")

# filepath = "ch0700.jl"
# Weave.weave(filepath; doctype = "md2html", out_path = :pwd,
#     template="templates/julia_html.tpl",
#     css="templates/skeleton_css.css")

# filepath = "ch0700quiz.jl"
# Weave.weave(filepath; doctype = "md2html", out_path = :pwd,
#     template="templates/julia_html.tpl",
#     css="templates/skeleton_css.css")

# filepath = "ch0900.jl"
# Weave.weave(filepath; doctype = "md2html", out_path = :pwd,
#     template="templates/julia_html.tpl",
#     css="templates/skeleton_css.css")

filepath = "ch0900quiz.jl"
Weave.weave(filepath; doctype = "md2html", out_path = :pwd,
    template="templates/julia_html.tpl",
    css="templates/skeleton_css.css",
    args = Dict("weaving" => true))

# Remove tmp folder
tmpfiles = Glob.glob("jl_??????")
for f in tmpfiles
    rm(f, force=true, recursive=true)
end


text:
using DelimitedFiles

mutable struct Node
    value::Int64
    previous::Node
    next::Node
    Node(v) =(x=new(v); x.previous=x; x.next=x; x)
end

function solve(input)
    raw = readdlm(input,'\n')
    parties = map(parseParty,raw)
    for party in parties
        @show party
        play(party[1],party[2])
        play(party[1],party[2]*100)
    end
end

function parseParty(s)
    m = match(r"(?<players>\d+) players; last marble is worth (?<points>\d+) points", s)
    players = tryparse(Int64,m[:players])
    points = tryparse(Int64,m[:points])
    return [players points]
end

function play(num_players, num_points)
    circle=Node(0)
    current = circle
    player = 1
    scores = zeros(Int64,num_players)
    for marble in 1:num_points
        if marble % 23 == 0
            scores[player]+= marble
            for i in 1:7
                current = current.previous
            end
            scores[player]+= current.value
            current.previous.next = current.next
            current.next.previous = current.previous
            current=current.next
        else
            current = current.next
            new_node = Node(marble)
            new_node.next = current.next
            new_node.previous = current
            current.next.previous = new_node
            current.next = new_node
            current = new_node
        end
        player = (player + 1 ) % num_players + 1
    end
    result = maximum(scores)
    @show result
end

solve("input.txt")


text:
using Optim

"""
    sample_from_cdf(cdf, support_point=0.5; n_retries=10_000, tol=1e-7)

given a cdf function, `sample_from_cdf` attempts to find a sample from the distribution it represents.
This is done via numerical inversion.

 - `cdf` the cdf function
 - `support_point` some point from with in the support of the distribution. This is used as a starting point for the numericla inversion. Ideally, but not nesc, it should be inside the support
 - `tol` how close the numeric approximation of the inverse has to be to the true inverse
 - `n_retries` how many time to restart in an attempt to find

### External Resources

 - https://sciencehouse.wordpress.com/2015/06/20/sampling-from-a-probability-distribution/
 - http://blog.quantitations.com/tutorial/2012/11/20/sampling-from-an-arbitrary-density/

"""
function sample_from_cdf(cdf, support_point::AbstractVector=[0.5]; n_retries=10_000, tol=1e-7, show_optimisation_details=false)
    u = rand()

    loss(xs) = (cdf(xs) - u).^2
    # find x such that u=inverse-cdf(x)
    # that x will be our sample
    start_point=copy(support_point)
    local result, best_x
    for _ in 1:n_retries
        result = optimize(loss, start_point, BFGS())

        if result.minimum < tol
            if show_optimisation_details
                println(result)
            end

            best_x = result.minimizer
            return best_x
        end
        start_point=jitter(start_point, support_point)
    end

    error("Failed to find an sample matching to $u in $(n_retries) attempts. Generally the solution to this is to provide good value for `support_point` that is somewhere in the center of the distribution.")
end

function sample_from_cdf(cdf, pdf, support_point::AbstractVector=[0.5]; n_retries=10_000, tol=1e-7, show_optimisation_details=false)
    u = rand()

    loss(xs) = (cdf(xs) - u).^2
    dloss!(store, xs) = store[1] = pdf(xs)*2*(cdf(xs) -u)

    # find x such that u=inverse-cdf(x)
    # that x will be our sample
    start_point=copy(support_point)
    local result, best_x
    for _ in 1:n_retries
        result = optimize(loss, dloss!, start_point, BFGS())

        best_x = result.minimizer
        if result.minimum < tol && pdf(best_x)>0.0
            if show_optimisation_details
                println(result)
            end

            return best_x
        end
        start_point=jitter(start_point, support_point)
    end

    error("Failed to find an sample matching to $u in $(n_retries) attempts. Generally the solution to this is to provide good value for `support_point` that is somewhere in the center of the distribution.")
end

# Start point jittering process
function jitter(cur_point, original_point)
    if rand()<0.1
        # reset it
        copy(original_point)
    else
        cur_point*4*(rand()-0.5)
        # Jitter the start point. This could be done better with a second start point to have a sense of scale
        # new value is somwehre between -2 times and 2x the previous
    end
end


text:
#TODO: Add doc string
@with_kw struct ArrowCar{A<:AbstractArray{Float64}, C<:Colorant} <: Renderable
    pos::A         = SVector(0.0, 0.0)
    angle::Float64 = 0.0
    length::Float64 = 4.8
    width::Float64 = 1.8
    color::C       = COLOR_CAR_OTHER
    text::String   = "" # some debugging text to print by the car
    id::Int        = 0
end

ArrowCar(pos::AbstractArray, angle::Float64=0.0; length = 4.8, width = 1.8,  color=COLOR_CAR_OTHER, text="", id=0) = ArrowCar(pos, angle, length, width, color, text, id)
ArrowCar(x::Real, y::Real, angle::Float64=0.0; length = 4.8, width = 1.8,  color=COLOR_CAR_OTHER, text="", id=0) = ArrowCar(SVector(x, y), angle, length, width, color, text, id)

function render!(rm::RenderModel, c::ArrowCar)
    x = c.pos[1]
    y = c.pos[2]
    add_instruction!(rm, render_vehicle, (x, y, c.angle, c.length, c.width, c.color))
    add_instruction!(rm, render_text, (c.text, x, y-c.width/2 - 2.0, 10, colorant"white"))
    return rm
end

id(ac::ArrowCar) = ac.id
pos(ac::ArrowCar) = ac.pos


text:
using AMQPClient
const VIRTUALHOST = "/"
const HOST = "127.0.0.1"
using JSON

function fib(n)
    if n == 0
        return 0
    elseif n == 1
        return 1
    else
        return fib(n - 1) + fib(n - 2)
    end
end

function receive()
    # 1. Create a connection to the localhost or 127.0.0.1 of virtualhost '/'
    connection(; virtualhost=VIRTUALHOST, host=HOST) do conn
        println("Connection established")
        # 2. Create a channel to send messages
        channel(conn, AMQPClient.UNUSED_CHANNEL, true) do chan
            # 3. Declare a queue
            println("Channel created")
            queue = "rpc_queue"
            success, queue_name, message_count, consumer_count = queue_declare(chan, queue; durable=true)

            # 4. Receiver requests and dispatch response
            on_receive = (msg) -> begin
                data = String(msg.data)
                # 4.1 Deserialize
                number = JSON.parse(data)
                println("Received requets to calculate fib[$number]")
                try
                    output = fib(number)
                    res = convert(Vector{UInt8}, codeunits(JSON.json(output)))
                    # 4.2 Ack the received request
                    basic_ack(chan, msg.delivery_tag)
                    # 4.3 Get the reply details from the message
                    reply_to = convert(String, msg.properties[:reply_to])
                    corr_id = convert(String, msg.properties[:correlation_id])
                    resp = Message(res, content_type="text/plain", delivery_mode=PERSISTENT, correlation_id=corr_id)
                    # Sent the response
                    basic_publish(chan, resp; exchange="", routing_key=reply_to)
                    println("Sent response for $number as $output")
                catch err
                    println(err)
                end
            end

            # 5. Consume requests for-ever
            success, consumer_tag = basic_consume(chan, queue, on_receive)

            while true
                sleep(1)
            end
        end
    end
end
# Don't exit on Ctrl-C
Base.exit_on_sigint(false)
try
    receive()
catch ex
    if ex isa InterruptException
        println("Interrupted")
    else
        println("Exception: $ex")
    end
end

max_forks_repo_name: tomashq/rabbitmq-tutorials


text:
@testset "c\"... x[N];\"" begin
	c``

	@eval c"""
	typedef char CcharArray[2];
	"""
	@test sizeof(c"CcharArray") == sizeof(Cchar)*2

	@eval c"""
	typedef struct {
		int i;
	} CstructArray[10];
	"""
	@test sizeof(c"CstructArray") == sizeof(Cint)*10

	@eval c"""
	struct CStruct {
		char c;
		int i;
	} __attribute__((packed));
	typedef struct CStruct CStructArray[2];
	"""
	@test sizeof(c"CStructArray") == (sizeof(Cchar)+sizeof(Cint))*2

	@eval c"""
	struct Opaque;
	typedef struct Opaque *PtrArray[3];
	"""
	@test sizeof(c"PtrArray") == sizeof(Cptr)*3

	pa = c"PtrArray"()
	@test sizeof(pa) == sizeof(Cptr)*3
	@test length(pa) == 3
	@test eltype(pa) <: c"struct Opaque *"
	@test pa[1] isa c"struct Opaque *"

	for i in eachindex(pa)
		pa = c"PtrArray"(pa, i => C_NULL)
		@test pa[i] == C_NULL
		pa = c"PtrArray"(pa, i => c"struct Opaque *"(1234))
		@test pa[i] != C_NULL
	end

	pa = c"PtrArray"()
	for ptr in pa
		@test ptr == C_NULL
	end

	ca = c"CStructArray"()
	@test length(ca) == 2
	@test eltype(ca) <: c"struct CStruct"
	@test ca[2] isa c"struct CStruct"

	@eval c"""
	typedef int CarrayArray[4][2];
	"""
	@test sizeof(c"CarrayArray") == sizeof(Cint)*4*2
	@test eltype(c"CarrayArray") <: c"int[2]"
	@test eltype(eltype(c"CarrayArray")) <: Cint

	ca = c"CarrayArray"()
	@test length(ca) == 4
	@test length(ca[1]) == 2
	for i in eachindex(ca), j in eachindex(ca[i])
		@test ca[i][j] == 0
	end

	@eval c"""
	struct UnknownLengthArray {
		int len;
		char str[];
	};
	"""
	@test sizeof(c"struct UnknownLengthArray") == sizeof(Cint)

	ula = c"struct UnknownLengthArray"()
	@test ula.len == 0
	@test ula.str isa Carray{Cchar}
	@test length(ula.str) == 0

	len = 10
	ptr = Libc.malloc(sizeof(c"struct UnknownLengthArray") + len*sizeof(Cchar))

	ula = Cptr{c"struct UnknownLengthArray"}(ptr)
	@test ula.len isa Cptr{Cint}
	@test ula.len == ptr

	ula.len[] = len
	@test ula.str isa Cptr{<:Carray{Cchar}}
	@test ula.str == ptr + sizeof(Cint)
	@test ula.str[1] isa Cptr{Cchar}
	@test ula.str[1] == ptr + sizeof(Cint)
	@test ula.str[5] isa Cptr{Cchar}
	@test ula.str[5] == ptr + sizeof(Cint) + 4*sizeof(Cchar)

	ula.str[1][] = 1
	ula.str[5][] = 5
	@test ula.str[1][] == 1
	@test ula.str[5][] == 5

	Libc.free(ptr)

	# https://github.com/analytech-solutions/CBinding.jl/issues/66
	@eval c"""
	struct file_format {
		unsigned char Header[2000];
	};
	"""
	x = c"struct file_format"(Header = "test")
	@test x.Header[2] == Int('e')
	@test x.Header[2000] == 0
end


text:
function stepAdvection!(
    mb :: ModelBlock,
    Δt :: Float64,
)

    ev = mb.ev
    fi = mb.fi
    tmpfi = mb.tmpfi
    co = mb.co
    cfg = ev.config
    for x = 1:2

        _intmx  = view(tmpfi._INTMX_, :, x)
        _xflx_U = view(fi._Xflx_U_, :, x)
        _xflx_V = view(fi._Xflx_V_, :, x)
        _xflx_W = view(fi._Xflx_W_, :, x)

        calDiffAdv_QUICKEST!(
            _xflx_U,
            _xflx_V,
            _xflx_W,

            _intmx,

            fi._u,
            fi._v,
            fi._w,

            co.amo,
            reshape(co.mtx[:Ks_H_U], :),
            reshape(co.mtx[:Ks_H_V], :),
            0.0,
            Δt,

            co.wksp,
        )

        _ADVx_ = view(fi._ADVX_, :, x)

        _ADVx_[:] = - (
              co.amo.T_DIVx_U * _xflx_U
            + co.amo.T_DIVy_V * _xflx_V
            + co.amo.T_DIVz_W * _xflx_W
        )

        @. _intmx += Δt * _ADVx_

    end

end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
struct RobustAdaptiveMetropolis{d, FT<:AbstractFloat, T <:
    Union{MVector{d,FT}, Vector{FT}}, SST<:StepSize,
    CT<:Cholesky{FT}} <: AdaptState
    L::CT
    α_opt::FT
    dx::T
    n_u::T
    step::SST
end

"""
    r = RobustAdaptiveMetropolis(x0, [α_opt, [step]])

Constructor for RobustAdaptiveMetropolis state.

# Arguments
- `x0`: The initial state vector
- `α_opt`: Desired mean accept rate; default `0.234`.
- `step`: Step size object; default `RAMStepSize(0.66,d)` where
  `d` is state dimension.

If `s` is `RWMState`, then proposal samples may be drawn calling
 `draw!(s, r)` and adaptation is performed with `adapt!(r, s, α)`.
"""

function RobustAdaptiveMetropolis(x0::T, α_opt::FT=FT(0.234),
  step::StepSize=RAMStepSize(FT(0.66),length(x0))) where {
      FT <: AbstractFloat, T <: AbstractVector{FT}}
    d = length(x0)
    L = Cholesky(eye(FT, d), :L, 0)
    dx = zeros(FT, d)
    n_u = zeros(FT, d)
    RobustAdaptiveMetropolis{d,FT,typeof(dx),typeof(step),typeof(L)}(L,
      α_opt, dx, n_u, step)
end

function RobustAdaptiveMetropolis(x0::T, α_opt::FT=FT(0.234),
  step::StepSize=RAMStepSize(FT(0.66),d)) where {d,
      FT <: AbstractFloat, T<:StaticArray{Tuple{d}, FT}}
    L = Cholesky(MMatrix{d,d}(eye(FT, d)), :L, 0)
    dx = zeros(MVector{d,FT})
    n_u = zeros(MVector{d,FT})
    RobustAdaptiveMetropolis{d,FT,typeof(dx),typeof(step),typeof(L)}(L,
      α_opt, dx, n_u, step)
end

@inline function adapt!(s::RobustAdaptiveMetropolis{d}, sr::RWMState{d},
      α::FT, k::Real) :: Nothing where {d, FT <: AbstractFloat}
    γ::FT = get(s.step, k)
    norm_u::FT = BLAS.nrm2(sr.u)
    norm_u = (norm_u == zero(FT)) ? one(FT) : norm_u

    s.n_u .= sr.u ./ norm_u
    #copy!(s.n_u, sr.u)
    #lmul!(one(FT)/norm_u, s.n_u)

    #mul!(s.dx, s.L.L, s.n_u)
    #lmul!(sqrt(γ*abs(α-s.α_opt)), s.dx)
    lowerTriInplaceMultiplyAdd!(s.dx, sqrt(γ*abs(α-s.α_opt)), s.L.factors,
              s.n_u, zero(FT), s.dx)
    if α >= s.α_opt
        lowrankupdate!(s.L, s.dx)
    else
        lowrankdowndate!(s.L, s.dx)
    end
    nothing
end

@inline function draw!(sr::RWMState, s::RobustAdaptiveMetropolis) :: Nothing
    draw!(sr, s.L)
    nothing
end


text:
# Sample a gas fitness DGP and estimate the parameters for different lenghts of
# time series
using StatsFuns
using StatsBase
using JLD2
N_est = 1

N = 50 ; Ngroups = 1; GTV =Ngroups; useStartVal = false
NnodesPerGroup = floor(N/Ngroups);groupsInds=zeros(Int64,N)
#Distribuisci i nodi tra i gruppi in parti uguali
for i =0:Ngroups-2 groupsInds[1 + Int64(NnodesPerGroup*i):Int64(NnodesPerGroup*(i+1))] = i+1 end;
groupsInds[find(groupsInds .==0)] = Ngroups
indTvGroups = falses(Ngroups);indTvGroups[1:GTV] = true
T = 200

#the values for the unconidtional means to be used in the DGP are obtained estimating
#an equally spaced degree sequence
tmpDegs = Array{Real,1}(linspace(2,N-2,N) )'; tmpdegsT = repmat(tmpDegs,2,1)  #ones(N)*N/2
tmpModel =  DynNets.SdErgm1(tmpdegsT,[zeros(N),zeros(GTV),zeros(GTV)],
                                    ones(Int,N))
est_times =  Array{Float64}(zeros( N_est))
est_conv_flag =  falses(N_est)
est_flag =  falses(N_est)

# Set the value of the parameters common to all replicate of DGP+Estimation
tmp,~ = DynNets.estimate(tmpModel )
realUm = tmp[1]
sim_B = Array(linspace(0.4,0.9,GTV+1)[2:end]) #
sim_A = Array(linspace(0.02,0.2,GTV+1)[2:end]) #
sim_W = realUm.*(1-sim_B[groupsInds])
 #  for n = 1:N_est
n = 1
        Y_T,fit_T  =  DynNets.sampleDynNetsBin1(T,"Fit-Multi-gas";N=N,
                                                    groupsInds=groupsInds,
                                                    AgasGroups =sim_A,
                                                    BgasGroups = sim_B,
                                                    WgasNodes = sim_W )
        degsT = squeeze(sum(Y_T,2),2)
        model_data = DynNets.SdErgm1(   degsT,
                                                [zeros(N),zeros(Ngroups),zeros(Ngroups)],
                                                groupsInds)
        #estimate
        useStartVal ?   start_val = [sim_W, sim_B, sim_A] : start_val = zeros(3,3)
        est_times[n] = @elapsed  hat_pars , conv_flag =
                            DynNets.estimate(model_data,start_values = start_val )
        indTvNodes =   .!(groupsInds .==0);
        est_conv_flag[n] = conv_flag
        est_flag[n] = true
#    end
println([N,T,sum(est_times),sum(est_conv_flag)/sum(est_flag)])

 ##


text:
using One_Dim_Chain
using Test

@testset "One_Dim_Chain.jl" begin
    # Write your tests here.
end


text:
module code_structs

export TreeNode
export Transition

mutable struct TreeNode
	parent::Union{TreeNode, Nothing}
	level::Int64 # Nível da raiz deve ser 1
	cost::Float64 # Custo acumulado no algoritmo de Viterbi
	info_bit::Int64 # Bit de informação associado ao nó
	codeword::Array{Int64, 1} # Palavra código do caminho até o nó
	TreeNode(parent, level, cost, info_bit) = new(
		parent,
		level,
		cost,
		info_bit
	)
end

struct Transition
	in_bit::Int64
	output::Array{Int64}
	next_state::Int64
	Transition(in_bit, output, next_state) = new(
		in_bit, output, next_state
	)
end

end


text:
using Test
using Yota
using Yota: Tape, tracked, record!, exec!, play!, grad, grad!, update!, setfield_nested!
using Yota: gettape, getid, getvalue, setvalue!
using Yota: Input, Constant, Call, Bcast, Assign

include("test_tracked.jl")
include("test_tape.jl")
include("test_macro.jl")
include("test_grad.jl")
include("test_update.jl")
include("test_examples.jl")


text:
################################################################################
#  Copyright 2022, Jose Antonio de la O Serna (UANL), Tom Van Acker (BASF)     #
################################################################################
# TFT.jl                                                                       #
# A Julia package for Taylor-Fourier Transform.                                #
# See http://github.com/timmyfaraday/TFT.jl                                    #
################################################################################

"""
TFT.ftft(prob::AbstractDTFTProblem)

Input:
- `prob::AbstractDTFTProblem`   | DTFT problem struct

Output:
- `sol::AbstractDTFTSolution`   | DTFT solution struct
"""
ftft(prob::AbstractDTFTProblem) = DTFTSolution(fast_estimator(prob), prob)

"""
TFT.ftft(s::Vector{<:Number}, t::Vector{<:Number}, D::Int, F::Real, K::Int)

Input:
- `s::Vector{<:Number}` | discrete signal [?]
- `t::Vector{<:Number}` | discrete time [s]
- `D::Int`              | maximum degree of the derivatives [-]
- `F::Real`             | fundamental frequency [Hz]
- `K::Int`              | degree of the o-spline [-]

Output:
- `sol::DTFTSolution`   | DTFT solution struct
"""
ftft(s::Vector{<:Number}, t::Vector{<:Number}, D::Int, F::Number, K::Int) =
    ftft(build_problem(s, t, D, F, K))

"""
TFT.ftft(s::Vector{<:Number}, t::Vector{<:Number}, h::Vector{<:Int}, D::Int, F::Real, K::Int)

Input:
- `s::Vector{<:Number}` | discrete signal [?]
- `t::Vector{<:Number}` | discrete time [s]
- `h::Vector{<:Int}`    | harmonic numbers [-]
- `D::Int`              | maximum degree of the derivatives [-]
- `F::Real`             | fundamental frequency [Hz]
- `K::Int`              | degree of the o-spline [-]

Output:
- `sol::DTFTSolution`   | DTFT solution struct
"""
ftft(s::Vector{<:Number}, t::Vector{<:Number}, h::Vector{<:Int}, D::Int, F::Number, K::Int) =
    ftft(build_problem(s, t, h, D, F, K))


text:
######################################################################
# random2.jl
#
# random policy for bearing only method
######################################################################
type RandomPolicy2 <: Policy end
export RandomPolicy2
function action(m::SearchDomain, v::Vehicle, o::Float64, f::DF, p::RandomPolicy2)
	xp = rand() * m.length
	yp = rand() * m.length
	return xp - v.x, yp - v.y, 0.0
end


text:
# -*- coding: utf-8 -*-
include("gerda_data_processing/init.jl")
include("gerda_data_processing/fct.jl")

sim_file        = "/remote/ceph2/group/gerda/data/simulation/gerda-mage-sim/gedet/intrinsic_bege/0nbb/edep/raw-gedet-intrinsic_bege-0nbb-edep-000.root"
config_dir      = "/res/Impurity_Scan/config-dep/"
plots_base_path = "pulses/sim/raw_0vbb/plots/"
base_path_raw   = "pulses/sim/raw_0vbb/"
base_path       = "pulses/sim/0vbb/"
n_sim_events    = 5000
sampling_time   = 1u"ns"
hits_threshold  = 0.005; # MeV
generate_wf    = true;
generate_cl_wf = true;

output_path     = base_path_raw*basename(sim_file)*"/";

if !isdir(output_path)
    mkpath(output_path)
end

file = TFile(sim_file)
tree = file["fTree"];

step_lim = 1000000;
steps = []
i = [1]
while i[1] <= length(tree)
    if i[1] + step_lim -1 <= length(tree)
        push!(steps, i[1]:i[1] + step_lim - 1)
    else
        push!(steps, i[1]:length(tree))
    end
    i[1] += step_lim
end

for s in steps

#     IJulia.clear_output(true)
    Base.run(`clear`)
    @info(">------------------------------------------------------<")
    @info("File to be processed: "*basename(sim_file))
    @info("File size: "*string(round(stat(sim_file).size/1024/1024, digits=2))*" mb")
    @info(">------------------------------------------------------<")
    step_count = findall(x->x == s, steps)[1]
    step_count_str = lpad(step_count, 3, "0");

    tree = file["fTree"][s];

    tt = Table(eventnumber = tree.eventnumber[:],
        hits_iddet   = tree.hits_iddet[:],
        hits_totnum  = tree.hits_totnum[:],
        hits_edep    = tree.hits_edep[:],
        hits_xpos    = tree.hits_xpos[:] .* 10, # mm
        hits_ypos    = tree.hits_ypos[:] .* 10, # mm
        hits_zpos    = tree.hits_zpos[:] .* 10) # mm

    #########################################
    # Filter events without energy deposition
    tt = tt |> @filter(_.hits_totnum != 0) |> Table;
    tt = tt |> @filter(length(unique(_.hits_iddet)) == 1) |> Table;
#     tt = tt |> @filter(sum(_.hits_edep) > 0.700 && sum(_.hits_edep) < 1.300) |> Table

    i = 1
    i_max = size(tt, 1)

    evt_num      = []
    multiplicity = []
    iddet        = []
    totnum       = []
    edep         = []
    position     = []

    p = Progress(i_max, dt=0.5,
                 barglyphs=BarGlyphs('|','█', ['▁' ,'▂' ,'▃' ,'▄' ,'▅' ,'▆', '▇'],' ','|',),
                 barlen=10)
    @info("File " * step_count_str * " of " * string(length(steps)))
    @info("Number of events = "*string(i_max))
    @info("Separating events by detector id and filter by minimum energy deposition threshold = "*string(hits_threshold)*" MeV")
    #########################################
    # Separate events by detector ###########
    while i <= i_max
        multi = 0
        for det in unique(tt[i].hits_iddet)
            index = findall(x->x==det, tt[i].hits_iddet)
            if sum( tt[i].hits_edep[ index ] ) >= hits_threshold
                multi += 1
                upside_down = 1
                if parameters[sim_to_channel[det][2]]["upside_down"] == true
                    upside_down = -1
                end

                hits_iddet_mapped = []
                for in in index
                    push!(hits_iddet_mapped, sim_to_channel[det][1])
                end
                push!(iddet, Array{Int64,1}(hits_iddet_mapped))
                push!(edep,  tt[i].hits_edep[ index ])
                push!(position, [ SVector{3}(([ tt[i].hits_xpos[k] .- parameters[sim_to_channel[det][2]]["detcenter_x"],
                                                tt[i].hits_ypos[k] .- parameters[sim_to_channel[det][2]]["detcenter_y"],
                                                upside_down .* (tt[i].hits_zpos[k] .- parameters[sim_to_channel[det][2]]["detcenter_z"] .+ upside_down * parameters[sim_to_channel[det][2]]["height"]/2)
                                    ] * u"mm")...) for k in index ])
            end
        end
        j = 1
        while j <= multi
            push!(evt_num, tt[i].eventnumber)
            push!(totnum, tt[i].hits_totnum)
            push!(multiplicity, multi)
            j += 1
        end
        next!(p)
        i += 1
    end
    tt = nothing
    #########################################
    # Create output table ###################
    filter_index = findall(x -> x == 1, multiplicity)

    data = Table(evtno = evt_num[filter_index],
        event_mult       = multiplicity[filter_index],
        detno            = iddet[filter_index],
        hits_totnum      = totnum[filter_index],
        edep             = VectorOfArrays(edep[filter_index] .*1000 *u"keV"),
        pos              = position[filter_index]
    )

    dets = []
    for detno in data.detno
        push!(dets, unique(detno)[1])
    end
    dets = unique(dets)
    p = Progress(length(dets), dt=0.5,
                    barglyphs=BarGlyphs('|','█', ['▁' ,'▂' ,'▃' ,'▄' ,'▅' ,'▆', '▇'],' ','|',),
                    barlen=10)
    @info("Creating files for each detector")
    for detno in dets
        det = channel_to_name[detno]
        t = data |> @select(:evtno, :event_mult, :detno, :hits_totnum, :edep, :pos) |> @filter(unique(_.detno)[1] == detno) |> Table;
        file_path = output_path * string(detno) * "-" * det * "/"
        if !isdir(file_path)
            mkpath(file_path)
        end
        file_path *= step_count_str * ".h5"
        HDF5.h5open(file_path, "w") do h5f
            LegendHDF5IO.writedata( h5f, "data", Table(evtno = t.evtno, event_mult = t.event_mult, detno = VectorOfArrays(t.detno), hits_totnum = t.hits_totnum,
                                                                    edep = VectorOfArrays(t.edep), pos = VectorOfArrays(t.pos)))
        end
        next!(p)
    end
end


text:
"generates variables for both `active` and `reactive` bus deltas"
function variable_bus_delta_abs(pm::AbstractPowerModel; kwargs...)
    variable_active_delta_abs(pm; kwargs...)
    variable_reactive_delta_abs(pm; kwargs...)
end

""
function variable_active_delta_abs(pm::AbstractPowerModel; nw::Int=pm.cnw, bounded::Bool=true, report::Bool=true)
    p_delta_abs = var(pm, nw)[:p_delta_abs] = @variable(pm.model,
        [i in ids(pm, :bus)], base_name="$(nw)_p_delta_abs",
        start = 0.0
    )

    if bounded
        for (i, bus) in ref(pm, nw, :bus)
            JuMP.set_lower_bound(p_delta_abs[i], 0.0)
            JuMP.set_upper_bound(p_delta_abs[i], 0.5)
        end
    end

    report && sol_component_value(pm, nw, :bus, :p_delta_abs, ids(pm, nw, :bus), p_delta_abs)
end

""
function variable_reactive_delta_abs(pm::AbstractPowerModel; nw::Int=pm.cnw, bounded::Bool=true, report::Bool=true)
     q_delta_abs = var(pm, nw)[:q_delta_abs] = @variable(pm.model,
        [i in ids(pm, :bus)], base_name="$(nw)_q_delta_abs",
        start = 0.0
    )

    if bounded
        for (i, bus) in ref(pm, nw, :bus)
            JuMP.set_lower_bound(q_delta_abs[i], 0.0)
            JuMP.set_upper_bound(q_delta_abs[i], 0.5)
        end
    end

    report && sol_component_value(pm, nw, :bus, :q_delta_abs, ids(pm, nw, :bus), q_delta_abs)
end

""
function variable_reactive_shunt(pm::AbstractPowerModel; nw::Int=pm.cnw, bounded::Bool=true, report::Bool=true)
    bs = var(pm, nw)[:bs] = @variable(pm.model,
        [i in ids(pm, nw, :shunt_var)], base_name="$(nw)_bs",
        start = comp_start_value(ref(pm, nw, :shunt, i), "bs_start")
    )

    if bounded
        for i in ids(pm, nw, :shunt_var)
            shunt = ref(pm, nw, :shunt, i)
            JuMP.set_lower_bound(bs[i], shunt["bmin"])
            JuMP.set_upper_bound(bs[i], shunt["bmax"])
        end
    end

    report && sol_component_value(pm, nw, :shunt, :bs, ids(pm, nw, :shunt_var), bs)
end

""
function variable_reactive_shunt(pm::AbstractWModels; nw::Int=pm.cnw, bounded::Bool=true, report::Bool=true)
    bs = var(pm, nw)[:bs] = @variable(pm.model,
        [i in ids(pm, nw, :shunt_var)], base_name="$(nw)_bs",
        start = comp_start_value(ref(pm, nw, :shunt, i), "bs_start")
    )

    wbs = var(pm, nw)[:wbs] = @variable(pm.model,
        [i in ids(pm, nw, :shunt_var)], base_name="$(nw)_wbs",
        start = 0.0
    )

    if bounded
        for i in ids(pm, nw, :shunt_var)
            shunt = ref(pm, nw, :shunt, i)
            JuMP.set_lower_bound(bs[i], shunt["bmin"])
            JuMP.set_upper_bound(bs[i], shunt["bmax"])
        end
    end

    report && sol_component_value(pm, nw, :shunt, :bs, ids(pm, nw, :shunt_var), bs)
    report && sol_component_value(pm, nw, :shunt, :wbs, ids(pm, nw, :shunt_var), wbs)
end

function variable_branch_flow_slack(pm::AbstractPowerModel; nw::Int=pm.cnw, bounded::Bool=true, report::Bool=true)
    sm_slack = var(pm, nw)[:sm_slack] = JuMP.@variable(pm.model,
        [l in ids(pm, nw, :branch_sm_active)], base_name="$(nw)_sm_slack",
        start = comp_start_value(ref(pm, nw, :branch, l), "sm_slack_start")
    )

    if bounded
        for (l,branch) in ref(pm, nw, :branch_sm_active)
            JuMP.set_lower_bound(sm_slack[l], 0.0)
        end
    end

    report && sol_component_value(pm, nw, :branch, :sm_slack, ids(pm, nw, :branch_sm_active), sm_slack)
end

function variable_vvm_delta(pm::AbstractPowerModel; nw::Int=pm.cnw, bounded::Bool=true, report::Bool=true)
    vvm_delta = var(pm, nw)[:vvm_delta] = JuMP.@variable(pm.model,
        [i in ids(pm, nw, :bus)], base_name="$(nw)_vvm_delta",
        start = comp_start_value(ref(pm, nw, :bus, i), "vvm_delta_start", 0.1)
    )

    report && sol_component_value(pm, nw, :bus, :vvm_delta, ids(pm, nw, :bus), vvm_delta)
end

function variable_pg_delta(pm::AbstractPowerModel; nw::Int=pm.cnw, bounded::Bool=true, report::Bool=true)
    pg_delta = var(pm, nw)[:pg_delta] = JuMP.@variable(pm.model,
        [i in ids(pm, nw, :gen)], base_name="$(nw)_pg_delta",
        start = comp_start_value(ref(pm, nw, :gen, i), "pg_delta_start", 0.1)
    )

    report && sol_component_value(pm, nw, :gen, :pg_delta, ids(pm, nw, :gen), pg_delta)
end

function variable_branch_contigency_flow_violation(pm::AbstractPowerModel; nw::Int=pm.cnw, bounded::Bool=true, report::Bool=true)
    branch_cont_flow_vio = var(pm, nw)[:branch_cont_flow_vio] = JuMP.@variable(pm.model,
        [i in 1:length(ref(pm, :branch_flow_cuts))], base_name="$(nw)_branch_cont_flow_vio",
        #start = comp_start_value(ref(pm, nw, :bus, i), "cont_branch_vio_start")
    )

    if bounded
        for i in 1:length(ref(pm, :branch_flow_cuts))
            JuMP.set_lower_bound(branch_cont_flow_vio[i], 0.0)
        end
    end

    #report && sol_component_value(pm, nw, :gen, :pg_delta, ids(pm, nw, :gen), pg_delta)
end

function variable_gen_contigency_flow_violation(pm::AbstractPowerModel; nw::Int=pm.cnw, bounded::Bool=true, report::Bool=true)
    gen_cont_flow_vio = var(pm, nw)[:gen_cont_flow_vio] = JuMP.@variable(pm.model,
        [i in 1:length(ref(pm, :gen_flow_cuts))], base_name="$(nw)_gen_cont_flow_vio",
        #start = comp_start_value(ref(pm, nw, :bus, i), "cont_branch_vio_start")
    )

    if bounded
        for i in 1:length(ref(pm, :gen_flow_cuts))
            JuMP.set_lower_bound(gen_cont_flow_vio[i], 0.0)
        end
    end

    #report && sol_component_value(pm, nw, :gen, :pg_delta, ids(pm, nw, :gen), pg_delta)
end

function variable_gen_contigency_capacity_violation(pm::AbstractPowerModel; nw::Int=pm.cnw, bounded::Bool=true, report::Bool=true)
    gen_cont_cap_vio = var(pm, nw)[:gen_cont_cap_vio] = JuMP.@variable(pm.model,
        [i in 1:length(ref(pm, :gen_contingencies))], base_name="$(nw)_gen_cont_cap_vio",
        #start = comp_start_value(ref(pm, nw, :bus, i), "cont_branch_vio_start")
    )

    if bounded
        for i in 1:length(ref(pm, :gen_contingencies))
            JuMP.set_lower_bound(gen_cont_cap_vio[i], 0.0)
        end
    end

    #report && sol_component_value(pm, nw, :gen, :pg_delta, ids(pm, nw, :gen), pg_delta)
end


text:
"""
    Model
Abstract type representing a model.

Child types of this abstract type are recognized as models by JITrench.

For more information, refer to the "Neural Networks" section of the documentation (here).
"""
abstract type Model end

"""
    plot_model(model::Model, inputs...,  to_file="model.png")
plot model by graphviz.
Note that JITrench builds the graph at the time of calculation, so some kind of input is required.
"""
function plot_model(model::Model, inputs...; to_file="model.png")
    y = model(inputs...)
    plot_graph(y, to_file=to_file)
end

#TODO: better impl
function parameters(model::Model)
    params = []
    for layer in layers(model)
        for param in parameters(layer)
            push!(params, param)
        end
    end
    return params
end

"""
    cleargrads!(model::Model; skip_uninit=false)

clear information about grad of the layers which is returned from `layers`.
if skip_uninit is true, parameters which is uninit(i.e nothing) is skipped.
"""
function cleargrads!(model::T, layers::Function; skip_uninit=false) where T <: Model
    for layer in layers(model)
        for param in parameters(layer)
            if param isa Nothing
                if skip_uninit
                    continue
                else
                    throw(DomainError("It is not possible to operate to clear uninit parameters."))
                end
            else
                param.grad = nothing
            end
        end
    end
end


text:
function show(io::IO, d::Distribution)
    @printf io "%s distribution\n" typeof(d)
    for parameter in typeof(d).names
        if isa(d.(parameter), AbstractArray)
            param = string(ucfirst(string(parameter)),
                           ":\n",
                           d.(parameter),
                           "\n")
        else
            param = string(ucfirst(string(parameter)),
                           ": ",
                           d.(parameter),
                           "\n")
        end
        print(io, param)
    end
end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.10.3
#   kernelspec:
#     display_name: Julia 1.6.4
#     language: julia
#     name: julia-1.6
# ---

# %% [markdown]
# # 通常の信頼区間とベイズ版信用区間の比較

# %%
using Distributions
using Roots
using Plots
using StatsFuns: logistic, logit

# \approx TAB → ≈
# \lessapprox TAB → ⪅
x ⪅ y = x < y || x ≈ y

# %%
# Bernoulli分布モデルでの通常の信頼区間

"""P値函数(exact版)"""
function pvalue_exact(n, p, k)
    bin = Binomial(n, p)
    p0 = pdf(bin, k)
    s = sum(pdf(bin, j) for j in support(bin) if pdf(bin, j) ⪅ p0)
    min(1, s)
end

"""P値函数(正規分布近似)"""
function pvalue_normal(n, p, k)
    bin = Binomial(n, p)
    normal = Normal(mean(bin), std(bin))
    min(1, 2cdf(normal, k), 2ccdf(normal, k))
end

"""信頼区間函数"""
function confidence_interval(pvalue_func, n, k; α = 0.05)
    f(t) = pvalue_func(n, logistic(t), k) - α
    CI = logistic.(find_zeros(f, -10.0, 10.0))
    if length(CI) < 2
        return 2k ≤ n ? [0.0, first(CI)] : [first(CI), 1.0]
    else
        return [first(CI), last(CI)]
    end
end

# %%
# Bernoulli分布モデルと共役事前分布に関するベイズ版信用区間

"""共役事後分布"""
posterior_dist(n, k; a = 1, b = a) = Beta(a + k, b + n - k)

"""ベイズ版P値函数"""
function pvalue_bayes(n, p, k; a = 1, b = a)
    posterior = posterior_dist(n, k; a, b)
    v0 = logpdf(posterior, p)
    f(t) = logpdf(posterior, logistic(t)) - v0
    m = params(posterior) |> ((α, β),) -> (α - 1)/(α + β - 2)
    if m ≤ 0
        s = ccdf(posterior, p)
    elseif m ≥ 1
        s = cdf(posterior, p)
    elseif p ≤ m
        q = logistic(find_zero(f, min(50, logit(m) + 1)))
        s = cdf(posterior, p) + ccdf(posterior, q)
    else
        q = logistic(find_zero(f, max(-50, logit(m) - 1)))
        s = cdf(posterior, q) + ccdf(posterior, p)
    end
    min(1, s)
end

"""ベイズ版信用区間函数"""
function credible_interval(n, k; α = 0.05, a = 1, b = a)
    confidence_interval((n, p, k) -> pvalue_bayes(n, p, k; a, b), n, k; α)
end

# %%
# 信頼区間とベイズ版信用区間の比較

function plot_cis(n, k; α = 0.05, a = 1, b = a)
    ci_exact = confidence_interval(pvalue_exact, n, k; α)
    ci_normal = confidence_interval(pvalue_normal, n, k; α)
    ci_bayes = credible_interval(n, k; α, a, b)

    bin = Binomial(n, k/n)
    μ, σ = mean(bin)/n, std(bin)/n
    p = range(max(1e-8, μ - 3.5σ), min(1-1e-8, μ + 4σ); length=400)
    plot(p, p -> pvalue_bayes(n, p, k; a, b); label="Bayesian p-value", c = 1, lw=1.5)
    plot!(p, p -> pvalue_normal(n, p, k); label="p-value (normal approx.)", c = 2, ls=:dash, lw=1.5)
    plot!(p, p -> pvalue_exact(n, p, k); label="p-value (exact)", c = 3, ls=:dashdot, lw=1.5)
    plot!(; ytick=0:0.05:1)
    plot!(ci_bayes, fill(1.2α, 2); label="Bayesian credible interval", c=1, lw=2)
    plot!(ci_normal, fill(1.0α, 2); label="conf. int. (normal approx.)", c=2, ls=:dash, lw=2)
    plot!(ci_exact, fill(0.8α, 2); label="confidence interval (exact)", c=3, ls=:dashdot, lw=2)
    title!("n = $n, k = $k, α = $α, prior = Beta($a, $b)"; titlefontsize=10)
end

# %%
plot_cis(10, 3)

# %%
plot_cis(100, 30)

# %%
plot_cis(1000, 300)

# %% [markdown]
# $n$ を大きくするにつれて通常の信頼区間とベイズ版信用区間は数値的によく一致するようになる。
#
# このようなことはシンプルなモデルでは一般的に広く成立している.

# %%
plot_cis(1000, 30)

# %%
plot_cis(1000, 3)

# %%
plot_cis(1000, 3; a = 0.5)

# %%
plot_cis(1000, 3; a = 0)

# %%

max_stars_repo_path: 0024/confidence and credible intervals.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
################################################################################
#                          Direction Cosine Matrices
################################################################################

export create_rotation_matrix
export ddcm, dcm2angle, dcm2quat

################################################################################
#                                  Functions
################################################################################

"""
    function create_rotation_matrix(angle::Number, axis::Symbol = :X)

Compute a rotation matrix that rotates a coordinate frame about the axis `axis`
by the angle `angle`. The `axis` must be one of the following symbols: `:X`,
`:Y`, or `:Z`.

# Example

```jldocstest
julia> create_rotation_matrix(pi/2, :X)
3×3 StaticArrays.SArray{Tuple{3,3},Float64,2,9}:
 1.0   0.0          0.0
 0.0   6.12323e-17  1.0
 0.0  -1.0          6.12323e-17
```

"""
function create_rotation_matrix(angle::Number, axis::Symbol = :X)
    cos_angle = cos(angle)
    sin_angle = sin(angle)

    if axis == :X
        dcm = DCM(1,      0,          0,
                  0, +cos_angle, +sin_angle,
                  0, -sin_angle, +cos_angle)'
        return dcm
    elseif axis == :Y
        dcm = DCM(+cos_angle, 0, -sin_angle,
                       0,     1,      0,
                  +sin_angle, 0, +cos_angle)'
        return dcm
    elseif axis == :Z
        dcm = DCM(+cos_angle, +sin_angle, 0,
                  -sin_angle, +cos_angle, 0,
                       0,          0,     1)'
        return dcm
    else
        error("axis must be :X, :Y, or :Z");
    end
end

################################################################################
#                                 Conversions
################################################################################

# Euler Angles
# ==============================================================================

"""
    function dcm2angle(dcm::DCM, rot_seq::Symbol=:ZYX)

Convert the DCM `dcm` to Euler Angles (see `EulerAngles`) given a rotation
sequence `rot_seq`.

The rotation sequence is defined by a `:Symbol`. The possible values are:
`:XYX`, `XYZ`, `:XZX`, `:XZY`, `:YXY`, `:YXZ`, `:YZX`, `:YZY`, `:ZXY`, `:ZXZ`,
`:ZYX`, and `:ZYZ`. If no value is specified, then it defaults to `:ZYX`.

# Example

```julia-repl
julia> D = DCM([1. 0. 0.; 0. 0. -1; 0. -1 0.]);

julia> dcm2angle(D,:XYZ)
ReferenceFrameRotations.EulerAngles{Float64}(1.5707963267948966, 0.0, -0.0, :XYZ)
```

"""
function dcm2angle(dcm::DCM, rot_seq::Symbol=:ZYX)
    if rot_seq == :ZYX

        EulerAngles(atan(+dcm[1,2],+dcm[1,1]),
                    asin(-dcm[1,3]),
                    atan(+dcm[2,3],+dcm[3,3]),
                    rot_seq)

    elseif rot_seq == :XYX

        EulerAngles(atan(+dcm[1,2],-dcm[1,3]),
                    acos(+dcm[1,1]),
                    atan(+dcm[2,1],+dcm[3,1]),
                    rot_seq)

    elseif rot_seq == :XYZ

        EulerAngles(atan(-dcm[3,2],+dcm[3,3]),
                    asin(+dcm[3,1]),
                    atan(-dcm[2,1],+dcm[1,1]),
                    rot_seq)

    elseif rot_seq == :XZX

        EulerAngles(atan(+dcm[1,3],+dcm[1,2]),
                    acos(+dcm[1,1]),
                    atan(+dcm[3,1],-dcm[2,1]),
                    rot_seq)

    elseif rot_seq == :XZY

        EulerAngles(atan(+dcm[2,3],+dcm[2,2]),
                    asin(-dcm[2,1]),
                    atan(+dcm[3,1],+dcm[1,1]),
                    rot_seq)

    elseif rot_seq == :YXY

        EulerAngles(atan(+dcm[2,1],+dcm[2,3]),
                    acos(+dcm[2,2]),
                    atan(+dcm[1,2],-dcm[3,2]),
                    rot_seq)

    elseif rot_seq == :YXZ

        EulerAngles(atan(+dcm[3,1],+dcm[3,3]),
                    asin(-dcm[3,2]),
                    atan(+dcm[1,2],+dcm[2,2]),
                    rot_seq)

    elseif rot_seq == :YZX

        EulerAngles(atan(-dcm[1,3],+dcm[1,1]),
                    asin(+dcm[1,2]),
                    atan(-dcm[3,2],+dcm[2,2]),
                    rot_seq)

    elseif rot_seq == :YZY

        EulerAngles(atan(+dcm[2,3],-dcm[2,1]),
                    acos(+dcm[2,2]),
                    atan(+dcm[3,2],+dcm[1,2]),
                    rot_seq)

    elseif rot_seq == :ZXY

        EulerAngles(atan(-dcm[2,1],+dcm[2,2]),
                    asin(+dcm[2,3]),
                    atan(-dcm[1,3],+dcm[3,3]),
                    rot_seq)

    elseif rot_seq == :ZXZ

        EulerAngles(atan(+dcm[3,1],-dcm[3,2]),
                    acos(+dcm[3,3]),
                    atan(+dcm[1,3],+dcm[2,3]),
                    rot_seq)

    elseif rot_seq == :ZYZ

        EulerAngles(atan(+dcm[3,2],+dcm[3,1]),
                    acos(+dcm[3,3]),
                    atan(+dcm[2,3],-dcm[1,3]),
                    rot_seq)

    else
        throw(ArgumentError("The rotation sequence :$rot_seq is not valid."))
    end
end

# Quaternions
# ==============================================================================

"""
    function dcm2quat(dcm::DCM)

Convert the DCM `dcm` to a quaternion. The type of the quaternion will be
automatically selected by the constructor `Quaternion` to avoid `InexactError`.

# Remarks

By convention, the real part of the quaternion will always be positive.
Moreover, the function does not check if `dcm` is a valid direction cosine
matrix. This must be handle by the user.

This algorithm was obtained from:

    http://www.euclideanspace.com/maths/geometry/rotations/conversions/matrixToQuaternion/

# Example

```julia-repl
julia> dcm = angle2dcm(pi/2,0.0,0.0,:XYZ);

julia> q   = dcm2quat(dcm)
Quaternion{Float64}:
  + 0.7071067811865476 + 0.7071067811865475.i + 0.0.j + 0.0.k
```

"""
function dcm2quat(dcm::DCM)
    if  tr(dcm) > 0
        # f = 4*q0
        f = 2sqrt(tr(dcm)+1)

        return Quaternion(f/4,
                          (dcm[2,3]-dcm[3,2])/f,
                          (dcm[3,1]-dcm[1,3])/f,
                          (dcm[1,2]-dcm[2,1])/f)

    elseif (dcm[1,1] > dcm[2,2]) && (dcm[1,1] > dcm[3,3])
        # f = 4*q1
        f = 2sqrt(1 + dcm[1,1] - dcm[2,2] - dcm[3,3])

        # Real part.
        q0 = (dcm[2,3]-dcm[3,2])/f

        # Make sure that the real part is always positive.
        s = (q0 > 0) ? +1 : -1

        return Quaternion(s*q0,
                          s*f/4,
                          s*(dcm[1,2]+dcm[2,1])/f,
                          s*(dcm[3,1]+dcm[1,3])/f)

    elseif (dcm[2,2] > dcm[3,3])
        # f = 4*q2
        f = 2sqrt(1 + dcm[2,2] - dcm[1,1] - dcm[3,3])

        # Real part.
        q0 = (dcm[3,1]-dcm[1,3])/f

        # Make sure that the real part is always positive.
        s = (q0 > 0) ? +1 : -1

        return Quaternion(s*q0,
                          s*(dcm[1,2]+dcm[2,1])/f,
                          s*f/4,
                          s*(dcm[3,2]+dcm[2,3])/f)

    else
        # f = 4*q3
        f = 2sqrt(1 + dcm[3,3] - dcm[1,1] - dcm[2,2])

        # Real part.
        q0 = (dcm[1,2]-dcm[2,1])/f

        # Make sure that the real part is always positive.
        s = (q0 > 0) ? +1 : -1

        return Quaternion(s*q0,
                          s*(dcm[1,3]+dcm[3,1])/f,
                          s*(dcm[2,3]+dcm[3,2])/f,
                          s*f/4)
    end
end

################################################################################
#                                  Kinematics
################################################################################

"""
    function ddcm(Dba::DCM, wba_b::AbstractArray)

Compute the time-derivative of the DCM `dcm` that rotates a reference frame `a`
into alignment to the reference frame `b` in which the angular velocity of `b`
with respect to `a`, and represented in `b`, is `wba_b`.

# Returns

The time-derivative of the DCM `Dba` (3x3 matrix of type `SMatrix{3,3}`).

# Example

```julia-repl
julia> D = DCM(Matrix{Float64}(I,3,3));

julia> ddcm(D,[1;0;0])
3×3 StaticArrays.SArray{Tuple{3,3},Float64,2,9}:
 0.0   0.0  0.0
 0.0   0.0  1.0
 0.0  -1.0  0.0
```

"""
function ddcm(Dba::DCM, wba_b::AbstractArray)
    # Auxiliary variable.
    w = wba_b

    # Check the dimensions.
    if length(wba_b) != 3
        throw(ArgumentError("The angular velocity vector must have three components."))
    end

    wx = SMatrix{3,3}(  0  , -w[3], +w[2],
                      +w[3],   0  , -w[1],
                      -w[2], +w[1],   0  )'

    # Return the time-derivative.
    -wx*Dba
end


text:
using DataGenerators

# grammar fasta;
@generator FastaGen begin
    # Since sequence is first non-terminal:
    start() = sequence()

    # sequence
    #    : section +
    #    ;
    sequence() = join(plus(section))

    # section
    #    : descriptionline
    #    | sequencelines
    #    | commentline
    #    ;
    section() = descriptionline()
    section() = sequencelines()
    section() = commentline()

    # sequencelines
    #    : SEQUENCELINE +
    #    ;
    sequencelines() = join(plus(SEQUENCELINE))

    # descriptionline
    #    : DESCRIPTIONLINE
    #    ;
    descriptionline() = DESCRIPTIONLINE__2()

    # commentline
    #    : COMMENTLINE
    #    ;
    commentline() = COMMENTLINE__2()

    # COMMENTLINE
    #    : ';' .*? EOL
    #    ;
    COMMENTLINE__2() = ";" * choose(String, ".*") * EOL()

    # DESCRIPTIONLINE
    #    : '>' TEXT ('|' TEXT)* EOL
    #    ;
    DESCRIPTIONLINE__2() = ">" * TEXT() * join(mult(intermediary_rule_1)) * EOL()
    intermediary_rule_1() = "|" * TEXT()

    # TEXT
    #    : (DIGIT | LETTER | SYMBOL) +
    #    ;
    TEXT() = join(plus(intermediary_rule_2))
    intermediary_rule_2() = DIGIT()
    intermediary_rule_2() = LETTER()
    intermediary_rule_2() = SYMBOL()

    # EOL
    #    : '\r'? '\n'
    #    ;
    EOL() = join(reps("\r", 0, 1)) * "\n"

    # fragment DIGIT
    # : [0-9]
    # ;
    DIGIT() = choose(String, "[0-9]")

    # fragment LETTER
    #    : [A-Za-z]
    #    ;
    LETTER() = choose(String, "[A-Za-z]")

    # fragment SYMBOL
    #    : '.' | '-' | '+' | '_' | ' ' | '[' | ']' | '(' | ')' | ',' | '/' | ':' | '&' | '\''
    #    ;
    SYMBOL() = "."
    SYMBOL() = "-"
    SYMBOL() = "+"
    SYMBOL() = "_"
    SYMBOL() = "["
    SYMBOL() = "]"
    SYMBOL() = "("
    SYMBOL() = ")"
    SYMBOL() = ","
    SYMBOL() = "/"
    SYMBOL() = ":"
    SYMBOL() = "&"
    SYMBOL() = "'"

    # SEQUENCELINE
    #    : LETTER + EOL
    #    ;
    SEQUENCELINE() = join(plus(LETTER)) * EOL()
end

gen = FastaGen()
# Not sure the fasta.g4 grammar really is correct since it does not seem to match what is
# on the FASTA page in Wikipedia: https://en.wikipedia.org/wiki/FASTA_format
choose(gen)


text:
using ArgParse
using HDF5
using KUnet
using CUDArt

function main()
    s = ArgParseSettings()
    @add_arg_table s begin
        "x"
        help = "HDF5 file for input"
        required = true
        "net"
        help = "Comma separated list of HDF5 layer files"
        required = true
        "y"
        help = "HDF5 file for desired output"
        required = true
        "out"
        help = "File prefix for trained output layer files"
        required = true
        "--nogpu"
        help = "Do not use gpu"
        action = :store_true
        "--batch"
        help = "Minibatch size"
        arg_type = Int
        default = 128
    end
    args = parse_args(s)
    batch = args["batch"]
    KUnet.gpu(!args["nogpu"])
    args["nogpu"] && blas_set_num_threads(20)
    x = h5read(args["x"], "/data")
    xx = x[:,1:batch]
    y = h5read(args["y"], "/data")
    yy = y[:,1:batch]
    if !args["nogpu"]
        xx = CudaArray(xx)
        yy = CudaArray(yy)
    end
    net = map(l->Op(l), split(args["net"],','))
    gc()
    @time KUnet.backprop(net, xx, yy)
    h5write("$(args["out"])1.h5", net[1])
    h5write("$(args["out"])2.h5", net[2])
    @time KUnet.backprop(net, xx, yy)
end

main()


text:
function button(a, window)
    robj = visualize(a).children[]
    const m2id = mouse2id(window)
    is_pressed = droprepeats(map(window.inputs[:key_pressed]) do isclicked
        isclicked && value(m2id).id == robj.id
    end)

    robj[:model] = const_lift(is_pressed, robj[:model], boundingbox(robj)) do ip, old, bb
        ip && return old
        scalematrix(Vec3f0(0.95))*old
    end
    robj, is_pressed
end

function toggle_button(a::Signal, b::Composable, window)
    toggle(b, window, signal=a)
    b, a
end

function toggle_button(a, b, window)
    id = Signal(0)
    ab_bool = toggle(id, window)
    a_b = map(ab_bool) do aORb
        aORb ? a : b
    end
    robj = visualize(a_b)
    push!(id, robj.children[].id)
    robj, ab_bool
end

function toggle(ispressed::Signal{Bool}, window, default=true; signal=Signal(default))
    preserve(map(ispressed) do ispressed
        ispressed && push!(signal, !Bool(value(signal)))
        nothing
    end)
    signal
end

function toggle{T<:Union{Int, Tuple, RenderObject}}(id1::Union{Signal{T}, T}, window, default=true; signal=Signal(default))
    is_clicked = droprepeats(map(window.inputs[:mouse_buttons_pressed]) do mbp
        if GLAbstraction.singlepressed(mbp, GLFW.MOUSE_BUTTON_LEFT)
            id2 = value(mouse2id(window))
            return is_same_id(id2, value(id1))
        end
        false
    end)
    toggle(is_clicked, window, default; signal=signal)
end

function toggle(robj::RenderObject, window, default=true; signal=Signal(default))
    toggle(robj, window, default, signal=signal)
end

function toggle(c::Context, window, default=true; signal=Signal(default))
    toggle(map(x->x.id, tuple(GLAbstraction.extract_renderable(c)...)), window, default, signal=signal)
end


text:
module ObjMesh

include("graphics.jl")
include("utils.jl")

using RayTracer: Vec3, parse_mtllib!
using .Graphics

cache = Dict()

export ObjectMesh

mutable struct ObjectMesh
    ##
    #Load and render OBJ model files
    ##
    min_coords::Vector{Float32}
    max_coords::Vector{Float32}
    vlists::Vector{Array{Float32, 3}}
    texclists::Vector{Array{Float32, 3}}
    clists::Vector{Vector{Vec3}}
    textures::Vector{Union{Vec3, Nothing}}
end

function ObjectMesh(file_path::String)
    ##
    #Load an OBJ model file
    #
    #Limitations:
    #- only one object/group
    #- only triangle faces
    ##

    # Comments
    # mtllib file_name
    # o object_name
    # v x y z
    # vt u v
    # vn x y z
    # usemtl mtl_name
    # f v0/t0/n0 v1/t1/n1 v2/t2/n2

    #logger.debug('loading mesh "%s"' % os.path.basename(file_path))

    # Attempt to load the materials library
    materials = _load_mtl(file_path)
    mesh_file = open(file_path, "r")

    verts = Vector{Vector{Float32}}()
    texs = Vector{Vector{Float32}}()
    normals = Vector{Vector{Float32}}()
    faces = []

    cur_mtl = ""

    # For each line of the input file
    for line in eachline(mesh_file)
        line = rstrip(line, [' ', '\r', '\n'])

        # Skip comments
        (startswith(line, '#') || line == "") && continue

        tokens = split(line, " ")
        tokens = map(t -> strip(t, ' '), tokens)
        tokens = collect(filter(t -> t != "", tokens))

        prefix = tokens[1]
        tokens = tokens[2:end]

        if prefix == "v"
            vert = map(v->parse(Float32, v), tokens)
            push!(verts, vert)
        end

        if prefix == "vt"
            tc = map(v -> parse(Float32, v), tokens)
            push!(texs, tc)
        end

        if prefix == "vn"
            normal = map(v -> parse(Float32, v), tokens)
            push!(normals, normal)
        end

        if prefix == "usemtl"
            mtl_name = tokens[1]
            if mtl_name in keys(materials)
                cur_mtl = mtl_name
            else
                cur_mtl = ""
            end
        end

        if prefix == "f"
            @assert length(tokens) == 3 "only triangle faces are supported"

            face = Vector{Vector{Int}}()
            for token in tokens
                indices = filter(t -> t != "", split(token, '/'))
                indices = map(i -> parse(Int, i), indices)
                @assert length(indices) == 2 || length(indices) == 3
                push!(face, indices)
            end

            push!(faces, [face, cur_mtl])
        end
    end

    # Sort the faces by material name
    sort!(faces, by=f->f[2])
    # Compute the start and end faces for each chunk in the model
    cur_mtl = nothing
    chunks = []
    for (idx, face) in enumerate(faces)
        face, mtl_name = face
        if mtl_name != cur_mtl
            if length(chunks) > 0
                chunks[end]["end_idx"] = idx-1
            end
            push!(chunks, Dict([
                "mtl"=>materials[mtl_name],
                "start_idx"=>idx,
                "end_idx"=>nothing
            ]))
            cur_mtl = mtl_name
        end
    end

    chunks[end]["end_idx"] = length(faces)
    num_faces = length(faces)
    # logger.debug('num verts=%d' % len(verts))
    # logger.debug('num faces=%d' % num_faces)
    # logger.debug('num chunks=%d' % len(chunks))

    # Create arrays to store the vertex data
    list_verts = zeros(Float32, 3, 3, num_faces)
    list_norms = zeros(Float32, 3, 3, num_faces)
    list_texcs = zeros(Float32, 3, 2, num_faces)
    list_colors = Vector{Vec3}()

    # For each triangle
    for (f_idx, face) in enumerate(faces)
        face, mtl_name = face

        # Get the color for this face
        f_mtl = materials[mtl_name]
        #NOTE: May get some failure/error here because of else condition
        #f_color = !isempty(f_mtl) ? f_mtl.color_diffuse : rgb(1f0)
        f_color = f_mtl.color_diffuse
        # For each tuple of indices
        for (l_idx, indices) in enumerate(face)
            # Note: OBJ uses 1-based indexing
            # and texture coordinates are optional
            if length(indices) == 3
                v_idx, t_idx, n_idx = indices
                vert = verts[v_idx]
                texc = texs[t_idx]
                normal = normals[n_idx]
            else
                v_idx, n_idx = indices
                vert = verts[v_idx]
                normal = normals[n_idx]
                texc = [0f0, 0f0]
            end

            list_verts[l_idx, :, f_idx] .= vert
            list_texcs[l_idx, :, f_idx] .= texc
            list_norms[l_idx, :, f_idx] .= normal
            push!(list_colors, f_color)
        end
    end

    # Re-center the object so that the base is at y=0
    # and the object is centered in x and z
    min_coords = minimum(minimum(list_verts, dims=3), dims=1)
    max_coords = minimum(maximum(list_verts, dims=3), dims=1)

    mean_coords = (min_coords .+ max_coords) / 2f0
    #NOTE: Why is mean coords correct but min coords has error margin of ~0.02?
    #@show mean_coords
    min_y = min_coords[2]
    mean_x = mean_coords[1]
    mean_z = mean_coords[3]
    list_verts[:, 2, :] .= list_verts[:, 2, :] .- min_y
    list_verts[:, 1, :] .= list_verts[:, 1, :] .- mean_x
    list_verts[:, 3, :] .= list_verts[:, 3, :] .- mean_z

    # Recompute the object extents after centering
    min_coords = minimum(minimum(list_verts, dims=3), dims=1)[1, :, 1]
    max_coords = minimum(maximum(list_verts, dims=3), dims=1)[1, :, 1]

    # Vertex list, one per chunk
    vlists = Vector{Array{Float32, 3}}()

    # Texture Coordinates list corresponding to vlist
    texclists = Vector{Array{Float32, 3}}()

    # Color list
    clists = Vector{Vector{Vec3}}()

    # Textures, one per chunk
    textures = Vector{Union{Vec3, Nothing}}()

    # For each chunk
    for chunk in chunks
        start_idx = chunk["start_idx"]
        end_idx = chunk["end_idx"]
        num_faces_chunk = end_idx - start_idx + 1

        # Create a vertex list to be used for rendering

        #TODO
        #vlist = pyglet.graphics.vertex_list(
            #3num_faces_chunk,
            #("v3f", reshape(list_verts[:, :, start_idx:end_idx], :)),
            #("t2f", reshape(list_texcs[:, :, start_idx:end_idx], :)),
            #("n3f", reshape(list_norms[:, :, start_idx:end_idx], :)),
            #("c3f", reshape(list_color[:, :, start_idx:end_idx], :)))

        mtl = chunk["mtl"]

        texture = mtl.texture_diffuse

        push!(vlists, list_verts[:, :, start_idx:end_idx])
        push!(texclists, list_texcs[:, :, start_idx:end_idx])
        push!(clists, list_colors[start_idx:end_idx])
        push!(textures, texture)
    end

    ObjectMesh(min_coords, max_coords, vlists, texclists, clists, textures)
end

function get(mesh_name::String)
    ##
    #Load a mesh or used a cached version
    ##

    # Assemble the absolute path to the mesh file
    file_path = get_file_path("src/meshes", mesh_name, "obj")

    file_path ∈ keys(cache) && (return cache[file_path])

    mesh = ObjectMesh(file_path)
    cache[file_path] = mesh

    return mesh
end

function _load_mtl(model_file::String)
    texture_diffuse = nothing
    model_dir, file_name = splitdir(model_file) .* ""

    # Determine the default texture path for the default material
    tex_name = split(file_name, '.')[1] * ""
    tex_path = get_file_path("src/textures", tex_name, "png")

    if ispath(tex_path)
        texture_diffuse = load_texture(tex_path)
    end

    materials = Dict{String, Union{NamedTuple,Nothing}}([
        # This is default material
        ""=> (color_diffuse = Vec3(1.0f0),
              color_ambient = Vec3(1.0f0),
              color_specular = Vec3(1.0f0),
              specular_exponent = 50.0f0,
              texture_ambient = nothing,
              texture_diffuse = texture_diffuse,
              texture_specular = nothing)
        ])

    mtl_path = model_file[1:end-4] * ".mtl"

    if ispath(mtl_path)
        parse_mtllib!(mtl_path, materials, Float32)
    end

    return materials
end

function render(obj_mesh::ObjectMesh)
    for (idx, vlist) in enumerate(obj_mesh.vlists)
        texture = obj_mesh.textures[idx]

        if !isnothing(texture)
            #gl.glEnable(gl.GL_TEXTURE_2D)
            #gl.glBindTexture(texture.target, texture.id)
        else
            #gl.glDisable(gl.GL_TEXTURE_2D)
        end

        #vlist.draw(gl.GL_TRIANGLES)
    end

    #gl.glDisable(gl.GL_TEXTURE_2D)
end

end #module


text:
using PastaQ
using ITensors
using Printf

# 1. Prepation of a thermal state
#
# In this example, we show how to prepare the finite-temperature state
# of a many-body system:
#
# ρ̂(β) = exp(-β Ĥ)
#
# where Ĥ is the Hamiltonian and β is the inverse temperature.
# We specificallty consider the one-dimensional Ising model
#
#   H = - ∑ᵢ σᶻ(i) σᶻ(i+1) - B ∑ᵢ σˣ(i)
#
# where B a the transverse magnetic field.

# 1a. Custom gates
#
# In order to build the thermal density operator, we implement the
# simplest flavor of imaginary-time evolution, breaking the operator
# exp(-βĤ) into a set of two-qubit and single-qubit gates, corresponding
# to the Ising interactions and the transverse field respetively. The
# time evolution to inverse temperature β is broken into elementary steps
# of size τ, where a gate is applied for each term appearing in the Hamiltonian.
#
# In this example, the quantum gates are not contained in the gate set of PastaQ.
# In order to extend, it is ony required to define the gate matrices using a
# format analogous to standard gates defined in gates.jl.

import PastaQ: gate

gate(::GateName"expτZZ"; τ::Float64) = exp(τ * kron(gate("Z"), gate("Z")))

gate(::GateName"expτX"; τ::Float64, B::Float64) = exp(τ * B * gate("X"))

# 1b. Generating the thermal state

N = 10    # Number of spins
B = 1.0   # Transverse magnetic field
β = 1.0   # Inverse temperature
τ = 0.005 # Trotter step

# Depth of the circuit
depth = β ÷ τ

# Ising interactions
zz_layer = [("expτZZ", (j, j + 1), (τ=τ,)) for j in 1:(N - 1)]
# Transverse field
x_layer = [("expτX", j, (τ=τ, B=B)) for j in 1:N]

# Build the gate structure
circuit = []
for d in 1:depth
  append!(circuit, zz_layer)
  append!(circuit, x_layer)
end

#
# 2. Run imaginary-time evolution towards the zero temperature
# ground state.
#

# 2a. Ground state energy with DMRG
#
# We compute the ground state energy by running DMRG
# on the Hamiltonian MPO, whose algorithm is implemented in
# ITensors.jl.

# In order to generate the MPO for the Hamiltonian, we leverage
# the ITensors.jl `AutoMPO()` function, which automatically
# generates the local MPO tensors from a set of pre-definend operators..
sites = siteinds("Qubit", N)
ampo = AutoMPO()
for j in 1:(N - 1)
  # Ising ZZ interactions
  ampo .+= -1, "Z", j, "Z", j + 1
end
for j in 1:N
  # Transverse field X
  ampo .+= -B, "X", j
end
# Generate Hamilotnian MPO
H = MPO(ampo, sites)

# Density-matrix renormalization group
dmrg_iter = 5      # DMRG steps
dmrg_cutoff = 1E-10   # Cutoff
Ψ0 = randomMPS(sites) # Initial state
sweeps = Sweeps(dmrg_iter)
maxdim!(sweeps, 10, 20, 30, 40, 50, 100)
cutoff!(sweeps, dmrg_cutoff)
# Run
println("Running DMRG to get ground state of transverse field Ising model:")
E, Ψ = dmrg(H, Ψ0, sweeps)
@printf("\nGround state energy:  %.8f  \n", E)
println("\n---------------------------------------\n")

#
# 2b. Run the imaginary-time circuit
#

β = 5.0 # Inverse temperature
Δ = 0.5  # Intermediate time-step
depth = Δ ÷ τ # Depth of the circuit
steps = β ÷ Δ # Total number of circuit application

# Initialize the density operator
ρ = PastaQ.identity_mpo(H)

println("Running imaginary time evolution to approximate the density matrix ρ = exp(-βH):")
for b in 1:steps
  # Run the circuit
  global ρ = runcircuit(ρ, circuit; cutoff=1E-12)

  # Normalize the density operatorr
  normalize!(ρ)

  # Measure the energy
  E_th = inner(ρ, H)

  @printf("β = %.1f : tr(ρH) = %.8f \n", (Δ * b), E_th)
end