text:
function maxpool2x2relu!(B, A)
    @avx for i₁ ∈ axes(B,1), i₂ ∈ axes(B,2), i₃ ∈ axes(B,3), i₄ ∈ axes(B,4)
        A₁ = A[2i₁-1,2i₂-1,i₃,i₄]
        A₂ = A[2i₁-1,2i₂  ,i₃,i₄]
        A₃ = A[2i₁  ,2i₂-1,i₃,i₄]
        A₄ = A[2i₁  ,2i₂  ,i₃,i₄]
        B[i₁,i₂,i₃,i₄] = max(max(max(A₁, A₂), max(A₃, A₄)), zero(eltype(A)))
    end
    B
end
function maxpool2x2relureverse!(A, B̄, B)
    @avx unroll=(1,1) for i₁ ∈ axes(B,1), i₂ ∈ axes(B,2), i₃ ∈ axes(B,3), i₄ ∈ axes(B,4)
        Bₘ = B[i₁,i₂,i₃,i₄]
        B̄ᵢ = B̄[i₁,i₂,i₃,i₄]
        A[2i₁-1,2i₂-1,i₃,i₄] = (A[2i₁-1,2i₂-1,i₃,i₄] == Bₘ) * B̄ᵢ
        A[2i₁-1,2i₂  ,i₃,i₄] = (A[2i₁-1,2i₂  ,i₃,i₄] == Bₘ) * B̄ᵢ
        A[2i₁  ,2i₂-1,i₃,i₄] = (A[2i₁  ,2i₂-1,i₃,i₄] == Bₘ) * B̄ᵢ
        A[2i₁  ,2i₂  ,i₃,i₄] = (A[2i₁  ,2i₂  ,i₃,i₄] == Bₘ) * B̄ᵢ
    end
end
function maxpool2x2relureversev2!(A, B̄)
    @avx unroll=(1,1) for i₁ ∈ axes(B̄,1), i₂ ∈ axes(B̄,2), i₃ ∈ axes(B̄,3), i₄ ∈ axes(B̄,4)
        A₁ = A[2i₁-1,2i₂-1,i₃,i₄]
        A₂ = A[2i₁-1,2i₂  ,i₃,i₄]
        A₃ = A[2i₁  ,2i₂-1,i₃,i₄]
        A₄ = A[2i₁  ,2i₂  ,i₃,i₄]
        Bₘ = max(max(max(A₁, A₂), max(A₃, A₄)), zero(eltype(A)))
        B̄ᵢ = B̄[i₁,i₂,i₃,i₄]
        A[2i₁-1,2i₂-1,i₃,i₄] = (A₁ == Bₘ) * B̄ᵢ
        A[2i₁-1,2i₂  ,i₃,i₄] = (A₂ == Bₘ) * B̄ᵢ
        A[2i₁  ,2i₂-1,i₃,i₄] = (A₃ == Bₘ) * B̄ᵢ
        A[2i₁  ,2i₂  ,i₃,i₄] = (A₄ == Bₘ) * B̄ᵢ
    end
end

@generated function halve12(x::Tuple{Vararg{<:Any,N}}) where {N}
    out = Expr(:tuple)
    for n in 1:N
        r = Expr(:ref, :x, n)
        if n ≤ 2
            r = Expr(:call, :>>>, r, Expr(:call, Expr(:curly, :Static, 1)))
        end
        push!(out.args, r)
    end
    out
end
function default_alloc_maxpool(mp::MaxPool2x2Layer, img::AbstractArray{T}) where {T}
    s = halve12(maybestaticsize(img))
    _, p = PaddedMatrices.size_permute_tuples(img)
    B = allocarray(T, s)
    PermutedDimsArray(B, p)
end
function (sp::StatckPointer)(::typeof(default_alloc_maxpool), mp::MaxPool2x2Layer, img::AbstractArray{T}) where {T}
    s = halve12(maybestaticsize(img))
    _, p = PaddedMatrices.size_permute_tuples(img)
    sp, B = PtrArray{T}(sp, s)
    sp, PermutedDimsArray(B, p)
end
function alloc_maxpool_pad12(mp::MaxPool2x2Layer, img::AbstractArray{T,4}) where {T}
    strunc = halve12(maybestaticsize(img))
    s = (vadd(strunc[1], Static{2}()), vadd(strunc[2], Static{2}()), strunc[3], strunc[4])
    _, p = PaddedMatrices.size_permute_tuples(img)
    B = allocarray(T, s)
    PermutedDimsArray(B, p)
end
function (sp::StatckPointer)(::typeof(default_alloc_maxpool), mp::MaxPool2x2Layer, img::AbstractArray{T,4}) where {T}
    s = halve12(maybestaticsize(img))
    _, p = PaddedMatrices.size_permute_tuples(img)
    sp, B = PtrArray{T}(sp, s)
    sp, PermutedDimsArray(B, p)
end

struct MaxPool2x2Layer{F} <: AbstractLayer
    f::F
end
MaxPool2x2Layer() = MaxPool2x2Layer(default_alloc_maxpool)

parameters(::MaxPool2x2Layer) = nothing
grad(::MaxPool2x2Layer) = nothing
returns(mp::MaxPool2x2Layer) = mp.o

function forward(sp::StackPointer, mp::MaxPool2x2Layer, img)

    sp, out = stack_pointer_call(mp.f, sp, img)

    maxpool2x2relu!(out, img)
    out, pb, out
end
forward!(mp::MaxPool2x2Layer, img) = (maxpool2x2relu!(returns(mp), img), img)
reverse_grad!(::MaxPool2x2Layer) = nothing
function reverse_chain!(mp::MaxPool2x2Layer, img, ōūt̄)
    maxpool2x2relureversev2!(img, ōūt̄)
end

# function maxpool2x2relureversev3!(A, B̄)
#     for i₁ ∈ axes(B̄,1), i₂ ∈ axes(B̄,2), i₄ ∈ axes(B̄,4)
#         @inbounds @simd ivdep for i₃ ∈ axes(B̄,3)
#             A₁ = A[2i₁-1,2i₂-1,i₃,i₄]
#             A₂ = A[2i₁-1,2i₂  ,i₃,i₄]
#             A₃ = A[2i₁  ,2i₂-1,i₃,i₄]
#             A₄ = A[2i₁  ,2i₂  ,i₃,i₄]
#             Bₘ = max(max(max(A₁, A₂), max(A₃, A₄)), zero(eltype(A)))
#             B̄ᵢ = B̄[i₁,i₂,i₃,i₄]
#             A[2i₁-1,2i₂-1,i₃,i₄] = (A₁ == Bₘ) * B̄ᵢ
#             A[2i₁-1,2i₂  ,i₃,i₄] = (A₂ == Bₘ) * B̄ᵢ
#             A[2i₁  ,2i₂-1,i₃,i₄] = (A₃ == Bₘ) * B̄ᵢ
#             A[2i₁  ,2i₂  ,i₃,i₄] = (A₄ == Bₘ) * B̄ᵢ
#         end
#     end
# end


text:
module CompatCartesian

export @ngenerate, @nsplat

macro ngenerate(itersym, returntypeexpr, funcexpr)
    if isa(funcexpr, Expr) && funcexpr.head == :macrocall && funcexpr.args[1] == Symbol("@inline")
        funcexpr = Base._inline(funcexpr.args[2])
    end
    isfuncexpr(funcexpr) || error("Requires a function expression")
    esc(_ngenerate(itersym, funcexpr))
end

function _ngenerate(itersym::Symbol, funcexpr::Expr)
    prototype = funcexpr.args[1]
    body = funcexpr.args[2]
    varname, T = get_splatinfo(prototype, itersym)
    ex = Expr(:$, itersym)
    sreplace!(body, itersym, ex)
    if !isempty(varname)
        prototype, body = _nsplat(prototype, body, varname, T, itersym)
    else
        body = Expr(:quote, body)
    end
    Expr(:stagedfunction, prototype, body)
end

macro nsplat(itersym, args...)
    if length(args) == 1
        funcexpr = args[1]
    elseif length(args) == 2
        funcexpr = args[2]
    else
        error("Wrong number of arguments")
    end
    if isa(funcexpr, Expr) && funcexpr.head == :macrocall && funcexpr.args[1] == Symbol("@inline")
        funcexpr = Base._inline(funcexpr.args[2])
    end
    isfuncexpr(funcexpr) || error("Second argument must be a function expression")
    prototype = funcexpr.args[1]
    body = funcexpr.args[2]
    varname, T = get_splatinfo(prototype, itersym)
    isempty(varname) && error("Last argument must be a splat")
    prototype, body = _nsplat(prototype, body, varname, T, itersym)
    esc(Expr(:stagedfunction, prototype, body))
end

function _nsplat(prototype, body, varname, T, itersym)
    varsym = Symbol(varname)
    prototype.args[end] = Expr(:..., Expr(:(::), varsym, T)) # :($varsym::$T...)
    varquot = Expr(:quote, varsym)
    bodyquot = Expr(:quote, body)
    newbody = quote
        $itersym = length($varsym)
        quote
            Base.Cartesian.@nexprs $($itersym) (d->$($(Expr(:quote, Symbol(varsym, "_d")))) = $($varquot)[d])
            $(Compat.CompatCartesian.resolvesplats!($bodyquot, $varquot, $itersym))
        end
    end
    prototype, newbody
end

# If using the syntax that will need "desplatting",
#     myfunction(A::AbstractArray, I::NTuple{N, Int}...)
# return the variable name (as a string) and type
function get_splatinfo(ex::Expr, itersym::Symbol)
    if ex.head == :call
        a = ex.args[end]
        if  isa(a, Expr) && a.head == :... && length(a.args) == 1
            b = a.args[1]
            if isa(b, Expr) && b.head == :(::)
                varname = string(b.args[1])
                c = b.args[2]
                if isa(c, Expr) && c.head == :curly && c.args[1] == :NTuple && c.args[2] == itersym
                    T = c.args[3]
                    return varname, T
                end
            end
        end
    end
    "", Void
end

resolvesplats!(arg, varname, N) = arg
function resolvesplats!(ex::Expr, varname, N::Int)
    if ex.head == :call
        for i = 2:length(ex.args)-1
            resolvesplats!(ex.args[i], varname, N)
        end
        a = ex.args[end]
        if isa(a, Expr) && a.head == :... && a.args[1] == Symbol(varname)
            ex.args[end] = :($varname[1]) # Expr(:ref, varname, 1)
            for i = 2:N
                push!(ex.args, :($varname[$i])) # Expr(:ref, varname, i))
            end
        else
            resolvesplats!(a, varname, N)
        end
    else
        for i = 1:length(ex.args)
            resolvesplats!(ex.args[i], varname, N)
        end
    end
    ex
end

isfuncexpr(ex::Expr) =
    ex.head == :function || (ex.head == :(=) && typeof(ex.args[1]) == Expr && ex.args[1].head == :call)
isfuncexpr(arg) = false

sreplace!(arg, sym, val) = arg
function sreplace!(ex::Expr, sym, val)
    for i = 1:length(ex.args)
        ex.args[i] = sreplace!(ex.args[i], sym, val)
    end
    ex
end
sreplace!(s::Symbol, sym, val) = s == sym ? val : s

end

ext: jl


text:
module Client

using HTTP, JSON3, Base64
using ..DataModel, ..Auth, ..Schema

const SERVER = Ref{String}("http://localhost:8080")
const AUTH_TOKEN = Ref{String}()

# function createUser(username, password)
#     body = (; username, password = base64encode(password))
#     resp = HTTP.post(string(SERVER[], "/user"), [], JSON3.write(body); cookies = true)
#     if HTTP.hasheader(resp, Auth.JWT_TOKEN_COOKIE_NAME)
#         AUTH_TOKEN[] = HTTP.header(resp, Auth.JWT_TOKEN_COOKIE_NAME)
#     end
#     return JSON3.read(resp.body, User)
# end

# function loginUser(username, password)
#     body = (; username, password = base64encode(password))
#     resp = HTTP.post(string(SERVER[], "/user/login"), [], JSON3.write(body); cookies = true)
#     if HTTP.hasheader(resp, Auth.JWT_TOKEN_COOKIE_NAME)
#         AUTH_TOKEN[] = HTTP.header(resp, Auth.JWT_TOKEN_COOKIE_NAME)
#     end
#     return JSON3.read(resp.body, User)
# end

# function createAlbum(name, artist, year, songs)
#     body = (; name, artist, year, songs)
#     resp = HTTP.post(string(SERVER[], "/album"), [Auth.JWT_TOKEN_COOKIE_NAME => AUTH_TOKEN[]], JSON3.write(body); cookies = true)
#     return JSON3.read(resp.body, Album)
# end

# function getAlbum(id)
#     resp = HTTP.get(string(SERVER[], "/album/$id"), [Auth.JWT_TOKEN_COOKIE_NAME => AUTH_TOKEN[]]; cookies = true)
#     return JSON3.read(resp.body, Album)
# end

# function updateAlbum(album)
#     resp = HTTP.put(string(SERVER[], "/album/$(album.id)"), [Auth.JWT_TOKEN_COOKIE_NAME => AUTH_TOKEN[]], JSON3.write(album); cookies = true)
#     return JSON3.read(resp.body, Album)
# end

# function deleteAlbum(id)
#     resp = HTTP.delete(string(SERVER[], "/album/$id"), [Auth.JWT_TOKEN_COOKIE_NAME => AUTH_TOKEN[]]; cookies = true)
#     return
# end

# function pickAlbumToListen()
#     resp = HTTP.get(string(SERVER[], "/"), [Auth.JWT_TOKEN_COOKIE_NAME => AUTH_TOKEN[]]; cookies = true)
#     return JSON3.read(resp.body, Album)
# end

function getSystemAlerts()
    resp = HTTP.get(string(SERVER[], "/sysai"); cookies = false)
    return JSON3.read(resp.body, SystemAlert)
end

# systemalerts or any mutable operations are strictly performed through backend
# router must not expose any mutable routes to end user
# without appropriate authentication

# this is demo app, and below function is only for testing purpose

function createSysAlert(id, entity, unit, status, last_updated, message)
    sysalert = SystemAlert(id, entity, unit, status, last_updated, message)
    # body = (; id, entity, unit, status, last_updated, message)
    Schema.insert!(sysalert)
    # resp = HTTP.post(string(SERVER[], "/album"), [Auth.JWT_TOKEN_COOKIE_NAME => AUTH_TOKEN[]], JSON3.write(body); cookies = true)
    # return JSON3.read(resp.body, Album)
end

end # module


text:
module IodWrapper

using JavaCall

const AbsoluteDate = @jimport org.orekit.time.AbsoluteDate
const CartesianOrbit = @jimport org.orekit.orbits.CartesianOrbit
const Class = @jimport java.lang.Class
const Frame = @jimport org.orekit.frames.Frame
const IodGibbs = @jimport org.orekit.estimation.iod.IodGibbs
const IodGooding = @jimport org.orekit.estimation.iod.IodGooding
const IodLambert = @jimport org.orekit.estimation.iod.IodLambert
const IodLaplace = @jimport org.orekit.estimation.iod.IodLaplace
const JString = @jimport java.lang.String
const KeplerianOrbit = @jimport org.orekit.orbits.KeplerianOrbit
const Object = @jimport java.lang.Object
const PV = @jimport org.orekit.estimation.measurements.PV
const PVCoordinates = @jimport org.orekit.utils.PVCoordinates
const Vector3D = @jimport org.hipparchus.geometry.euclidean.threed.Vector3D

export IodGibbs
export IodGooding
export IodLambert
export IodLaplace
export equals
export estimate
export get_class
export get_range1
export get_range2
export get_range3
export hash_code
export notify
export notify_all
export to_string
export wait

include("iod_gibbs.jl")
include("iod_gooding.jl")
include("iod_lambert.jl")
include("iod_laplace.jl")

end


text:
using Documenter, ReservoirComputing

makedocs(
    modules=[ReservoirComputing],
    clean=true,doctest=false,
    sitename = "ReservoirComputing.jl",
    format = Documenter.HTML(analytics = "UA-90474609-3",
                         assets = ["assets/favicon.ico"],
                         canonical="https://reservoircomputing.sciml.ai/stable/"),
    pages = [
        "ReservoirComputing.jl" => "index.md",
        "General Settings" => Any[
            "Changing Training Algorithms" => "general/different_training.md",
            "Altering States" => "general/states_variation.md",
            "Generative vs Predictive" => "general/predictive_generative.md",
        ],
        "Echo State Network Tutorials" => Any[
            "Lorenz System Forecasting" => "esn_tutorials/lorenz_basic.md",
            "Mackey-Glass Forecasting on GPU" => "esn_tutorials/mackeyglass_basic.md",
            "Using Different Layers" => "esn_tutorials/change_layers.md",
            "Using Different Reservoir Drivers" => "esn_tutorials/different_drivers.md",
            #"Using Different Training Methods" => "esn_tutorials/different_training.md",
            "Deep Echo State Networks" => "esn_tutorials/deep_esn.md",
            "Hybrid Echo State Networks" => "esn_tutorials/hybrid.md",
            ],
        "Reservoir Computing with Cellular Automata" => "reca_tutorials/reca.md",
        "API Documentation"=>Any[
            "Training Algorithms" => "api/training.md",
            "States Modifications" => "api/states.md",
            "Prediction Types" => "api/predict.md",
            "Echo State Networks" => "api/esn.md",
            "ESN Layers" => "api/esn_layers.md",
            "ESN Drivers" => "api/esn_drivers.md",
            "ReCA" => "api/reca.md"
            ]
        ]
    )

deploydocs(
   repo = "github.com/SciML/ReservoirComputing.jl.git";
   push_preview = true
)

max_stars_repo_name: MartinuzziFrancesco/EchoStateNetwork

max_issues_repo_name: MartinuzziFrancesco/EchoStateNetwork

max_forks_repo_name: MartinuzziFrancesco/EchoStateNetwork


text:
######################################
## Auxilary traffic model functions ##
######################################
"""
`traffic_constants` create maximal traffic densities and speeds matrices

**Input parameters**
* `OSMmap` : OpenStreetMapX MapData object with road network data
* `density_factor` : road length reserved for one vehicle
"""
function traffic_constants(OSMmap::MapData,
                           density_factor::Float64 = 5.0)
    #Create maximal densitites matrix
    max_densities = get_max_densities(OSMmap, density_factor)
    #Create speed limit matrix
    max_speeds = OpenStreetMapX.get_velocities(OSMmap)
    return max_densities, max_speeds
end

"""
`get_max_densities` calculate maximal traffic densities on all edges

**Input parameters**
* `OSMmap` : OpenStreetMapX MapData object with road network data
* `density_factor` : road length reserved for one vehicle
"""
function get_max_densities(OSMmap::MapData,
                           density_factor::Float64 = 5.0)
    roads_lanes = Dict{Int,Int}()
    for r in OSMmap.roadways
        OpenStreetMapX.haslanes(r) ? lanes = OpenStreetMapX.getlanes(r) : lanes = 1
        roads_lanes[r.id] = lanes
    end
    segments = OpenStreetMapX.find_segments(OSMmap.nodes, OSMmap.roadways, OSMmap.intersections)
    segments = [[OSMmap.v[s.node0], OSMmap.v[s.node1], roads_lanes[s.parent]] for s in segments]
    sparse_lanes = SparseArrays.sparse([x[1] for x in segments],
                                        [x[2] for x in segments],
                                        [x[3] for x in segments],
                                        length(OSMmap.v),length(OSMmap.v))
    return OSMmap.w .* sparse_lanes / density_factor
end

"""
`init_traffic_variables` create data used in calculating velocities change during simulation

**Input parameters**
* `OSMmap` : OpenStreetMapX MapData object with road network data
* `Agents` : set of agents created with generate_agents function
"""
function init_traffic_variables(OSMmap::MapData,
                                Agents::Vector{Agent})
    #Dictionary with initital densities
    initial_densities = StatsBase.countmap([a.edge for a in Agents if a.start_time == 0.0])
    initial_speeds = OpenStreetMapX.get_velocities(OSMmap)
    return initial_densities, initial_speeds
end

"""
`update_weights!` change speeds in given speed_matrix for edges listed in new_densities

**Input parameters**
* `speed_matrix` : matrix with average speeds on edges
* `new_densities` : dictionary with edges as keys and new traffic density as value
* `max_densities` : matrix with maximal densitites on edges
* `V_max` : matrix with maximal speeds on edges
* `V_min` : minimal speed on road
"""
function update_weights!(speed_matrix::SparseMatrixCSC{Float64,Int},
                        dens::Dict{Tuple{Int,Int},Int},
                        max_densities::SparseMatrixCSC{Float64,Int},
                        V_max::SparseMatrixCSC{Float64,Int},
                        V_min::Float64 = 1.0)
    for k in keys(dens)
        v1, v2 = k
        speed_matrix[v1,v2]  = (V_max[v1,v2] - V_min)* max((1 - dens[k]/max_densities[v1,v2]), 0.0) + V_min
    end
end

"""
`update_weights_and_events!` update events time and weights in speed matrix

**Input parameters**
* `inAgents` : set of agents created with generate_agents function
* `agents_pos` : vector with agents current position (edges)
* `events` : vector with events time
* `speed_matrix` : matrix with average speeds on edges
* `edges` : one or two edges that took part in an event
* `dens` : dictionary with current density on edges
* `max_densities` : matrix with maximal densities
* `V_max` : matrix with speed limits
* `V_min` : minimal speed allowed on road
"""
function update_weights_and_events!(inAgents::Vector{Agent},
                                    agents_pos::Vector{Tuple{Int,Int}},
                                    events::Vector{Float64},
                                    speed_matrix::SparseMatrixCSC{Float64,Int},
                                    edges::Vector{Tuple{Int,Int}},
                                    dens::Dict{Tuple{Int,Int},Int},
                                    max_densities::SparseMatrixCSC{Float64,Int},
                                    V_max::SparseMatrixCSC{Float64,Int},
                                    V_min::Float64 = 1.0)
    for edge in edges
        v1, v2 = edge
        old_speed = speed_matrix[v1,v2]
        #Update speed matrix
        new_speed = (V_max[v1,v2] - V_min)* max((1 - dens[edge]/max_densities[v1,v2]), 0.0) + V_min
        speed_matrix[v1,v2]  = new_speed
        #Adjust event time for agents in modified edges
        speed_factor = old_speed/new_speed
        for i in findall(x-> x==edge, agents_pos)
            @inbounds events[i] = events[i]*speed_factor
        end
    end
end

"""
`next_edge` returns time required to reach next junction for one agent

**Input parameters**
* `Agent` : single instance of Agent
* `speeds` : current speeds matrix
* `lengths` : road lengths matrix
"""
function next_edge(Agent::Agent,
                    speeds::AbstractMatrix,
                    lengths::AbstractMatrix)
    if !Agent.active
        return Inf
    end
    e = Agent.edge
    event_time = lengths[e[1],e[2]]/speeds[e[1],e[2]]
    return event_time
end

"""
`next_edge` returns time required to reach next junction for set of agents

**Input parameters**
* `Agents` : set of agents created with generate_agents function
* `speeds` : current speeds matrix
* `lengths` : road lengths matrix
"""
function next_edge(Agents::Vector{Agent},
                    speeds::AbstractMatrix,
                    lengths::AbstractMatrix)
    return [next_edge(a,speeds,lengths) for a in Agents]
end

"""
`update_event_agent!` update densities matrix, progress agents to next edge and deactivate agents

**Input parameters**
* `inAgent` : agent connected with occuring event
* `curr_time` : current simulation time
* `densities` : current traffic densitites dictionary
* `vertices_map` : mapping from nodes to vertices
"""
function update_event_agent!(inAgent::Agent,
                            curr_time::Float64,
                            densities::Dict{Tuple{Int,Int}, Int},
                            vertices_map::Dict{Int,Int})
    #Decrease density on previous edge
    p_edge = inAgent.edge
    densities[p_edge] -= 1
    #Update agent  position
    if length(inAgent.route) == 2
        #Disable agent and set travelling time
        inAgent.active = false
        inAgent.travel_time = curr_time - inAgent.start_time
        #Return agent's last edge
        return [p_edge]
    else
        #Update agent route and current edge
        inAgent.route = inAgent.route[2:end]
        c_edge = (inAgent.edge[2], vertices_map[inAgent.route[2]])
        inAgent.edge = c_edge
        #Add density on new edge
        haskey(densities, c_edge) ? densities[c_edge] += 1 : densities[c_edge] = 1
        #Return previous and current agent's edge
        return [p_edge, c_edge]
    end
end


text:
### unitful.jl
#
# Copyright (C) 2019 Mosè Giordano.
#
# Maintainer: Mosè Giordano <mose AT gnu DOT org>
# Keywords: uncertainty, error propagation, physics
#
# This file is a part of Measurements.jl.
#
# License is MIT "Expat".
#
### Commentary:
#
# This file integration with Unitful.jl.
#
### Code:

import .Unitful: AbstractQuantity, unit, ustrip
function Measurements.measurement(a::T, b::T) where {T<:AbstractQuantity}
    u = unit(a)
    return measurement(ustrip(u, a), ustrip(u, b)) * u
end
Measurements.measurement(a::AbstractQuantity{T1,D,U1},
                         b::AbstractQuantity{T2,D,U2}) where {T1,T2,D,U1,U2} =
                             measurement(promote(a, b)...)

function Measurements.value(x::AbstractQuantity)
    u = unit(x)
    return value(ustrip(u, x)) * u
end

function Measurements.uncertainty(x::AbstractQuantity)
    u = unit(x)
    return uncertainty(ustrip(u, x)) * u
end


text:
#
# Protein - TMAO (compare new and old implementations)
#

include("../../src/ComplexMixtures.jl")
using PDBTools

dir="/home/leandro/Drive/Alunos/Ivan/PCCP_revision"

atoms = PDBTools.readPDB("$dir/6Mnative.pdb")

protein = PDBTools.select(atoms,"protein")
solute = ComplexMixtures.Selection(protein,nmols=1 )

water = PDBTools.select(atoms,"water")
solvent = ComplexMixtures.Selection(water,natomspermol=3)

options = ComplexMixtures.Options(n_random_samples=10,lastframe=200)

trajectory = ComplexMixtures.Trajectory("$dir/6Mnative.dcd",solute,solvent)

@time lcP = ComplexMixtures.mddf(trajectory,options)


text:
using Test
using HMatrices
using LinearAlgebra
using Random
using StaticArrays

using HMatrices: RkMatrix

include(joinpath(HMatrices.PROJECT_ROOT,"test","testutils.jl"))

Random.seed!(1)

m = 2000
n = 2000

X    = rand(SVector{3,Float64},m)
Y    = [rand(SVector{3,Float64}) for _  in 1:n]
splitter  = CardinalitySplitter(nmax=40)
Xclt      = ClusterTree(X,splitter)
Yclt      = ClusterTree(Y,splitter)
adm       = StrongAdmissibilityStd(eta=3)
rtol      = 1e-5
comp      = PartialACA(rtol=rtol)
K         = laplace_matrix(X,Y)
H         = assemble_hmat(K,Xclt,Yclt;adm,comp,threads=false,distributed=false)

H_full = Matrix(H;global_index=false)

α = rand() - 0.5
β = rand() - 0.5
@testset "hmul!" begin
    C  = deepcopy(H);
    tmp = β*H_full + α*H_full*H_full
    HMatrices.hmul!(C,H,H,α,β,PartialACA(;atol=1e-6))
    @test Matrix(C;global_index=false) ≈ tmp
end

@testset "gemv" begin

    α = rand()-0.5
    β = rand()-0.5
    T = eltype(H)
    m,n = size(H)
    x = rand(T,n)
    y = rand(T,m)

    @testset "serial" begin
        exact  = β*y + α*H_full*x
        approx = mul!(copy(y),H,x,α,β;threads=false,global_index=false)
        @test exact ≈ approx
        exact  = β*y + α*Matrix(H;global_index=true)*x
        approx = mul!(copy(y),H,x,α,β;threads=false,global_index=true)
        @test exact ≈ approx
    end

    @testset "threads" begin
        exact  = β*y + α*H_full*x
        approx = mul!(copy(y),H,x,α,β;threads=true,global_index=false)
        @test exact ≈ approx
        exact  = β*y + α*Matrix(H;global_index=true)*x
        approx = mul!(copy(y),H,x,α,β;threads=false,global_index=true)
        @test exact ≈ approx
    end
end

max_stars_repo_name: WaveProp/HMatrices

max_forks_repo_name: WaveProp/HMatrices


text:
"""
alex_test()

Test for duplication of Table
"""
function alex_test()
    (historye, condhiste, alphanorme) = vtst()
    u0 = ones(2)
    maxit = 20
    maxm = 2
    vdim = 3 * maxm + 3
    Vstore = zeros(2, vdim)
    VstoreS = zeros(2, 2*maxm + 4)
    m = 2
    aout = aasol(alexfp!, u0, m, Vstore; rtol = 1.e-10)
    alexerr = (
        reldiff(aout.history, historye) +
        reldiff(aout.stats.condhist, condhiste) +
        reldiff(aout.stats.alphanorm, alphanorme)
    )
    aoutS = aasol(alexfp!, u0, m, VstoreS; rtol = 1.e-10)
    alexerrS = (
        reldiff(aoutS.history, aout.history) +
        reldiff(aoutS.stats.condhist, aout.stats.condhist) +
        reldiff(aoutS.stats.alphanorm, aout.stats.alphanorm)
    )
    alexok2=(alexerrS < 1.e-15)
    solerr = reldiff(aout.history, historye)
#   put this back to solerr < 1.e-5  once 1.7 is fixed
    solok = (solerr < 1.e-3)
    solok || println("alex solution error","  ",solerr)
    conderr = reldiff(aout.stats.condhist, condhiste)
#   Something's broken with 1.7 in windoze/linux
#   put this back to conderr < 1.e-5  once 1.7 is fixed
    condok = (conderr < 1.e-2)
    condok || println("alex condition error","  ",conderr)
    aerr = reldiff(aout.stats.alphanorm, alphanorme)
#   put this back to aerr < 1.e-5  once 1.7 is fixed
    aok = (aerr < 1.e-5)
    aok || println("alex coefficient error")
    aout.idid || println("idid is wrong for m=2")
    alexok2 = alexok2 && solok && condok && aok
    aout = aasol(alexfp!, u0, 0, Vstore; rtol = 1.e-10)
    aout.idid && println("idid is wrong for m=0")
    alexok0 = ~aout.idid && (aout.errcode == 10)
    alexok = alexok2 && alexok0
    return alexok
end

function vtst()
    historye = [
        6.50111e-01
        4.48661e-01
        2.61480e-02
        7.25389e-02
        1.53107e-04
        1.18512e-05
        1.82476e-08
        1.04804e-13
    ]
    condhiste = [
        1.00000e+00
        2.01556e+10
        1.37776e+09
        3.61344e+10
        2.54947e+11
        3.67672e+10
    ]
    alphanorme = [
        1.00000e+00
        4.61720e+00
        2.15749e+00
        1.18377e+00
        1.00000e+00
        1.00171e+00
    ]
    return (historye, condhiste, alphanorme)
end

function alexfp!(G, u)
G[1]=cos(.5*(u[1]+u[2]))
G[2]=G[1]+ 1.e-8 * sin(u[1]*u[1])
return G
end


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

name = "Fontconfig"
version = v"2.13.1-0"

# Collection of sources required to build Fontconfig
sources = [
    "https://www.freedesktop.org/software/fontconfig/release/fontconfig-2.13.1.tar.gz" =>
    "9f0d852b39d75fc655f9f53850eb32555394f36104a044bb2b2fc9e66dbbfa7f",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
if [[ $target == *"linux"* ]]; then
    wget "https://github.com/hessammehr/UtilLinuxBuilder/releases/download/v2.33.0/UtilLinuxBuilder.v2.33.0.$target.tar.gz"
    tar xzf UtilLinuxBuilder*.tar.gz --directory $prefix
fi
wget https://github.com/hessammehr/GperfBuilder/releases/download/v3.1.0-0/Gperf.v3.1.0.x86_64-linux-gnu.tar.gz
tar xzf Gperf*.tar.gz --directory $prefix

cd $prefix/lib/pkgconfig/

if [ ! -f zlib.pc ]; then
(cat <<EOF
prefix=/workspace/destdir
exec_prefix=${prefix}
libdir=${exec_prefix}/lib
sharedlibdir=${libdir}
includedir=${prefix}/include

Name: zlib
Description: zlib compression library
Version: 1.2.11

Requires:
Libs: -L${libdir} -L${sharedlibdir} -lz
Cflags: -I${includedir}
EOF
) > zlib.pc
fi

(cat <<EOF
prefix=/workspace/destdir
exec_prefix=${prefix}
libdir=${exec_prefix}/lib
sharedlibdir=${libdir}
includedir=${prefix}/include

Name: bzip2
Description: bzip2 compression library
Version: 1.0.6

Requires:
Libs: -L${libdir} -L${sharedlibdir} -lbz2
Cflags: -I${includedir}
EOF
) > bzip2.pc

sed -i 's/2.9.1/22.1.16/g' freetype2.pc
cd $WORKSPACE/srcdir/fontconfig-2.13.1/
./configure --prefix=$prefix --host=$target
sed -i 's/\/uuid//g' test/Makefile
make
make install

if [ $target = "x86_64-unknown-freebsd11.1" ]; then
cd $WORKSPACE/srcdir
cd $prefix/lib/pkgconfig/
ls
cd ..
ls
cd pkgconfig/

fi

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:i686, libc=:glibc),
    Linux(:x86_64, libc=:glibc),
    Linux(:aarch64, libc=:glibc),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf),
    Linux(:i686, libc=:musl),
    Linux(:x86_64, libc=:musl),
    Linux(:aarch64, libc=:musl),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf),
    MacOS(:x86_64),
    Windows(:i686),
    Windows(:x86_64)
]

# The products that we will ensure are always built
products(prefix) = [
    ExecutableProduct(prefix, "fc-validate", :fc_validate),
    LibraryProduct(prefix, "libfontconfig", :libfontconfig),
    ExecutableProduct(prefix, "fc-cat", :fc_cat),
    ExecutableProduct(prefix, "fc-conflist", :fc_conflist),
    ExecutableProduct(prefix, "fc-match", :fc_match),
    ExecutableProduct(prefix, "fc-pattern", :fc_pattern),
    ExecutableProduct(prefix, "fc-scan", :fc_scan),
    ExecutableProduct(prefix, "fc-list", :fc_list),
    ExecutableProduct(prefix, "fc-cache", :fc_cache),
    ExecutableProduct(prefix, "fc-query", :fc_query)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "https://github.com/hessammehr/ExpatBuilder/releases/download/v2.2.6/build_ExpatBuilder.v2.2.6.jl",
    "https://github.com/JuliaGraphics/FreeTypeBuilder/releases/download/v2.9.1-2/build_FreeType2.v2.9.1.jl",
    "https://github.com/staticfloat/Bzip2Builder/releases/download/v1.0.6-1/build_Bzip2.v1.0.6.jl",
    "https://github.com/bicycle1885/ZlibBuilder/releases/download/v1.0.4/build_Zlib.v1.2.11.jl"
]

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

max_stars_repo_name: hessammehr/FontconfigBuilder

max_issues_repo_name: hessammehr/FontconfigBuilder

max_forks_repo_name: hessammehr/FontconfigBuilder


text:
unpack(x) = (@assert length(x)==1; x[1])
isjl(path) = splitext(path)[2] == ".jl"

function walkpaths(f, root)
    for (dir, _, files) in walkdir(root)
        for file in files
            path = joinpath(dir, file)
            f(path)
        end
    end
end

ismacrocall(x) = false
ismacrocall(ex::Expr) = ex.head == :macrocall

iscall(ex::Expr) = ex.head == :call
iscall(x) = false

iscall_fun_or_macro(x) = iscall(x) | ismacrocall(x)

calle(ex::Expr) = (@assert iscall_fun_or_macro(ex); ex.args[1])
callarg_unique(ex::Expr) = ex |> callargs |> unpack

function callargs(ex::Expr)
    @assert iscall_fun_or_macro(ex)
    ex.args[2:end]
end

const INCOMPLETE = parse("[")
isincomplete(ex) = ex == INCOMPLETE

immutable MacroCall{f} end
immutable Call{f} end
immutable AnyExpr end
immutable NoExpr end
immutable Lambda end

extrait(ex::Any) = NoExpr()
function extrait(ex::Expr)
    if iscall(ex)
        f = calle(ex)
        typeof(f) == Expr || return Call{f}()  # HACK: julia does not like expressions as type parameters
        return AnyExpr()
    elseif ismacrocall(ex)
        return MacroCall{calle(ex)}()
    elseif ex.head == Symbol("->")
        return Lambda()
    else
        return AnyExpr()
    end
end


text:
"""
The `mode` type parameter must be either `:inflate` or `:deflate`.
"""
mutable struct Source{mode,T<:BufferedInputStream}
    input::T
    zstream::ZStream
    state::State
    reset_on_end::Bool
end

# inflate source constructors
# ---------------------------

function InflateSource(input::T, raw::Bool, gzip::Bool,
                       reset_on_end::Bool) where T<:BufferedInputStream
    return Source{:inflate,T}(
        input,
        init_inflate_zstream(raw, gzip),
        initialized,
        reset_on_end)
end

function InflateSource(input::BufferedInputStream, bufsize::Integer, raw::Bool,
                       gzip::Bool, reset_on_end::Bool)
    return InflateSource(input, raw, gzip, reset_on_end)
end

function InflateSource(input::IO, bufsize::Integer, raw::Bool, gzip::Bool,
                       reset_on_end::Bool)
    input_stream = BufferedInputStream(input, bufsize)
    return InflateSource(input_stream, raw, gzip, reset_on_end)
end

function InflateSource(input::Vector{UInt8}, bufsize::Integer, raw::Bool,
                       gzip::Bool, reset_on_end::Bool)
    return InflateSource(BufferedInputStream(input), raw, gzip, reset_on_end)
end

"""
    ZlibInflateInputStream(input[; <keyword arguments>])

Construct a zlib inflate input stream to decompress gzip/zlib data.

# Arguments
* `input`: a byte vector, IO object, or BufferedInputStream containing compressed data to inflate.
* `bufsize::Integer=8192`: input and output buffer size.
* `raw::Bool=falso`: if true, data is raw compress data, without zlib metadata
* `gzip::Bool=true`: if true, data is gzip compressed; if false, zlib compressed.
* `reset_on_end::Bool=true`: on stream end, try to find the start of another stream.
"""
function ZlibInflateInputStream(input; bufsize::Integer=8192, raw::Bool=false,
                                gzip::Bool=true, reset_on_end::Bool=true)
    return BufferedInputStream(
        InflateSource(input, bufsize, raw, gzip, reset_on_end),
        bufsize)
end

# deflate source constructors
# ---------------------------

function DeflateSource(
        input::T, raw::Bool, gzip::Bool, level::Integer, mem_level::Integer,
        strategy) where T<:BufferedInputStream
    return Source{:deflate,T}(
        input,
        init_deflate_zstream(raw, gzip, level, mem_level, strategy),
        initialized,
        false)
end

function DeflateSource(input::BufferedInputStream, bufsize::Integer, raw::Bool,
                       gzip::Bool, level::Integer, mem_level::Integer, strategy)
    return DeflateSource(input, raw, gzip, level, mem_level, strategy)
end

function DeflateSource(input::IO, bufsize::Integer, raw::Bool, gzip::Bool,
                       level::Integer,
                       mem_level::Integer, strategy)
    input_stream = BufferedInputStream(input, bufsize)
    return DeflateSource(input_stream, raw, gzip, level, mem_level, strategy)
end

function DeflateSource(input::Vector{UInt8}, bufsize::Integer, raw::Bool,
                       gzip::Bool, level::Integer, mem_level::Integer, strategy)
    return DeflateSource(BufferedInputStream(input), raw, gzip, level, mem_level, strategy)
end

"""
    ZlibDeflateInputStream(input[; <keyword arguments>])

Construct a zlib deflate input stream to compress gzip/zlib data.

# Arguments
* `input`: a byte vector, IO object, or BufferedInputStream containing data to compress.
* `bufsize::Integer=8192`: input and output buffer size.
* `raw::Bool=false`: if true, data is raw compress data, without zlib metadata
* `gzip::Bool=true`: if true, data is gzip compressed; if false, zlib compressed.
* `level::Integer=6`: compression level in 1-9.
* `mem_level::Integer=8`: memory to use for compression in 1-9.
* `strategy=Z_DEFAULT_STRATEGY`: compression strategy; see zlib documentation.
"""
function ZlibDeflateInputStream(input;
                                bufsize::Integer=8192,
                                raw::Bool=false,
                                gzip::Bool=true,
                                level::Integer=6,
                                mem_level::Integer=8,
                                strategy=Z_DEFAULT_STRATEGY)
    return BufferedInputStream(
        DeflateSource(input, bufsize, raw, gzip, level, mem_level, strategy),
        bufsize)
end

"""
    readbytes!(source, buffer, from, to)

Read bytes from the zlib stream to a buffer. Satisfies the BufferedStreams
source interface.
"""
function BufferedStreams.readbytes!(
        source::Source{mode},
        buffer::Vector{UInt8},
        from::Int, to::Int) where mode
    if source.state == finalized
        return 0
    elseif source.state == finished && source.reset_on_end
        reset!(source)
    end

    @trans source (
        initialized => inprogress,
        inprogress  => inprogress
    )

    fillbuffer!(source.input)
    source.zstream.next_out = pointer(buffer, from)
    source.zstream.avail_out = to - from + 1
    _, n_out = process(
        source,
        mode == :deflate && eof(source.input) ? Z_FINISH : Z_NO_FLUSH)
    return n_out
end

function process(source::Source{mode}, flush) where mode
    @assert source.state == inprogress
    # counter of processed input/output bytes
    n_in = n_out = 0
    input = source.input
    zstream = source.zstream

    #println("--- Source{", mode, "} ---")
    @label process
    zstream.next_in = pointer(input)
    zstream.avail_in = BufferedStreams.available_bytes(input)
    old_avail_in = zstream.avail_in
    old_avail_out = zstream.avail_out
    if mode == :inflate
        ret = inflate!(zstream, flush)
    else
        ret = deflate!(zstream, flush)
    end
    n_in += old_avail_in - zstream.avail_in
    n_out += old_avail_out - zstream.avail_out
    input.position += old_avail_in - zstream.avail_in

    if ret == Z_OK
        if zstream.avail_in == 0
            if BufferedStreams.fillbuffer!(input) == 0
                flush = Z_FINISH
            end
            @goto process
        end
    elseif ret == Z_STREAM_END
        @trans source inprogress => finished
    elseif ret == Z_BUF_ERROR
        # could not consume more input or produce more output
    elseif ret < 0
        zerror(zstream, ret)
    else
        @assert false
    end

    return n_in, n_out
end

@inline function Base.eof(source::Source{mode}) where mode
    if source.state == initialized ||
        (mode == :inflate && source.state == finished && source.reset_on_end)
        return eof(source.input)
    end
    return source.state == finished || source.state == finalized
end

function Base.close(source::Source{mode}) where mode
    if source.state == finalized
        isopen(source.input) && close(source.input)
        return
    end
    if mode == :inflate
        @zcheck end_inflate!(source.zstream)
    else
        @zcheck end_deflate!(source.zstream)
    end
    @trans source (
        initialized => finalized,
        inprogress  => finalized,
        finished    => finalized
    )
    close(source.input)
    return
end

function reset!(source::Source{mode}) where mode
    if mode == :inflate
        @zcheck reset_inflate!(source.zstream)
    else
        @zcheck reset_deflate!(source.zstream)
    end
    @trans source (
        initialized => initialized,
        inprogress  => initialized,
        finished    => initialized,
        finalized   => initialized
    )
    return source
end

# For backwards compatibility with 0.2 releases.
InflateSource(input::BufferedInputStream, gzip::Bool, reset_on_end::Bool) =
    InflateSource(input, false, gzip, reset_on_end)
InflateSource(input::BufferedInputStream, bufsize::Integer, gzip::Bool,
              reset_on_end::Bool) =
    InflateSource(input, bufsize, false, gzip, reset_on_end)
InflateSource(input::IO, bufsize::Integer, gzip::Bool, reset_on_end::Bool) =
    InflateSource(input, bufsize, false, gzip, reset_on_end)
InflateSource(input::Vector{UInt8}, bufsize::Integer, gzip::Bool,
              reset_on_end::Bool) =
    InflateSource(input, bufsize, false, gzip, reset_on_end)

DeflateSource(input::BufferedInputStream, gzip::Bool, level::Integer,
              mem_level::Integer, strategy) =
    DeflateSource(input, false, gzip, level, mem_level, strategy)
DeflateSource(input::BufferedInputStream, bufsize::Integer, gzip::Bool,
              level::Integer, mem_level::Integer, strategy) =
    DeflateSource(input, bufsize, false, gzip, level, mem_level, strategy)
DeflateSource(input::IO, bufsize::Integer, gzip::Bool, level::Integer,
              mem_level::Integer, strategy) =
    DeflateSource(input, bufsize, false, gzip, level, mem_level, strategy)
DeflateSource(input::Vector{UInt8}, bufsize::Integer, gzip::Bool,
              level::Integer, mem_level::Integer, strategy) =
    DeflateSource(input, bufsize, false, gzip, level, mem_level, strategy)


text:
"""
Provides virtual interface for the gray-box environment and simulator.
"""
module GrayBox

using Random
using Distributions

export Simulation,
       Sample,
       Environment,
       EnvironmentSample,
       environment,
       transition!

"""
    GrayBox.Simulation

Abstract base type for a gray-box simulation.
"""
abstract type Simulation end

"""
    GrayBox.Sample

Holds sampled value and log-probability.
"""
mutable struct Sample{T}
    value::T
    logprob::Real
end

"""
    GrayBox.Environment

Alias type for a dictionary of gray-box environment distributions.

     e.g., `Environment(:variable_name => Sampleable)`
"""
const Environment = Dict{Symbol, Sampleable}

"""
    GrayBox.EnvironmentSample

Alias type for a single environment sample.

    e.g., `EnvironmentSample(:variable_name => Sample(value, logprob))`
"""
const EnvironmentSample = Dict{Symbol, Sample}

"""
    environment(sim::GrayBox.Simulation)

Return all distributions used in the simulation environment.
"""
function environment(sim::Simulation)::Environment end

"""
  GrayBox.State

Alias type for an explicit gray-box simulation state.
"""
const State = Union{Vector{Real}, Nothing}

"""
  GrayBox.state

Get current state of the gray-box simulation.
"""
function state(sim::Simulation)::State end

"""
    transition!(sim::Union{GrayBox.Simulation, GrayBox.Simulation})::Real

Given an input `sample::EnvironmentSample`, apply transition and return the transition log-probability.
"""
function transition!(sim::Simulation, sample::EnvironmentSample)::Real end

"""
    transition!(sim::GrayBox.Simulation)::Real

Apply a transition step, and return the transition log-probability (used with `ASTSeedAction`).
"""
function transition!(sim::Simulation)::Real end

end # module GrayBox


text:
# Gather results from the Vulcan cluster at HLRS.
# The setup is described in `vulcan.md`. This script postprocesses
# the data generated that way by moving output files to get the
# same structure as obtained by running everything locally as
# described in the `README.md` files in each subdirectory.

using DelimitedFiles

function postprocess_vulcan(directory)
  for (root, dirs, files) in walkdir(directory)
    for dir in dirs
      # directories are named "run_IDENTIFIER"
      if !startswith(dir, "run_")
        continue
      end

      id = dir[5:end]
      jobscript = "jobscript_" * id * ".sh"
      if !(jobscript in files) && !("jobscript.sh" in files)
        @warn "There is no jobscript '$jobscript' associated to '$dir' in $root"
        continue
      end

      for (subroot, subdirs, subfiles) in walkdir(joinpath(root, dir))
        for file in subfiles
          if endswith(file, ".dat") || endswith(file, ".txt")
            src = joinpath(subroot, file)
            dst = joinpath(root, file)
            @info "Copying $src to $dst"
            cp(src, dst, force=true)
          end
        end
      end
    end
  end
end

postprocess_vulcan(joinpath(@__DIR__, "numerical_fluxes"))
postprocess_vulcan(joinpath(@__DIR__, "overintegration"))
postprocess_vulcan(joinpath(@__DIR__, "precompute_logs"))
postprocess_vulcan(joinpath(@__DIR__, "primitive_variables"))


text:
function kde!(points::Array{Float64,2}, ks::Array{Float64,1}, weights::Array{Float64,1})
  Nd, Np = size(points)
  if (length(ks) == 1)
    ks = repmat(ks,Nd)
  end

  ks = ks.^2 # Guassian only at this point, taking covariance
  weights = weights./sum(weights);
  #bwsize = length(ks);

  makeBallTreeDensity(points,weights,ks)

  #if (length())
end

function kde!(points::Array{Float64,2}, ks::Array{Float64,1})
  Nd, Np = size(points)
  weights = ones(Np)
  kde!(points, ks, weights)
end

function kde!(points::Array{Float64,1}, ks::Array{Float64,1})
  Np = length(points)
  pts = zeros(1,Np)
  pts[1,:] = points
  weights = ones(Np)
  kde!(pts, ks, weights)
end

function getPoints(bd::BallTreeDensity)
    pts=zeros(bd.bt.dims, bd.bt.num_points)
    perm = bd.bt.permutation[(bd.bt.num_points + 1):end]
    bd.bt.num_points, bd.bt.dims
    res = reshape(bd.bt.centers[(bd.bt.dims*bd.bt.num_points+1):end], bd.bt.dims, bd.bt.num_points)
    ##for i in 1:bd.bt.dims
      pts[:, perm ] = res[:,:]
    ##end
    pts = pts[:,1:bd.bt.num_points]
    return pts
end

function getType(bd::BallTreeDensity)
    return 0
end

#s = zeros(dens.D,dens.N);
#s(:,double(dens.perm(dens.N + (1:dens.N)))+1) = dens.bandwidth(:,dens.N + (1:dens.N));
function getBW(bd::BallTreeDensity, ind::Array{Int64,1}=zeros(Int64,0))
  if length(ind)==0
    ind=1:bd.bt.num_points
  end
  s = zeros(bd.bt.dims,bd.bt.num_points)
  perm = bd.bt.permutation[(bd.bt.num_points + 1):end]
  res = reshape(bd.bandwidth[(bd.bt.dims*bd.bt.num_points + 1):end],bd.bt.dims,bd.bt.num_points)
  s[:, perm ] = res[:,:]
  s = s[:,ind]
  s = sqrt( s );  # stddev gaussian covariance
  return s
end

function getWeights(bd::BallTreeDensity, ind::Array{Int64,1}=zeros(Int64,0))
  if length(ind)==0
    ind=1:bd.bt.num_points
  end
  perm = bd.bt.permutation[(bd.bt.num_points+1):end]
  wts = zeros(bd.bt.num_points)
  wts[perm] = bd.bt.weights[(bd.bt.num_points+1):end]
  wts = wts[ind]
  return wts
end

function marginal(bd::BallTreeDensity, ind::Array{Int64,1})
  pts = getPoints(bd)
  if size(bd.bandwidth,2) > 2*bd.bt.num_points
    sig = getBW(bd)
  else
    sig = getBW(bd,[1])
  end
  wts = getWeights(bd)
  p = kde!(pts[ind,:],sig[ind], wts)
end

function randKernel(N::Int64, M::Int64, t::Int64)
  return randn(N,M)
end

function sample(npd::BallTreeDensity, Npts::Int64)
  pts = getPoints(npd)
  points = zeros(npd.bt.dims, Npts)
  ind = zeros(Int64,Npts)
  bw  = getBW(npd)
  w = getWeights(npd);
  w = cumsum(w)
  w = w./w[end]
  randnums = randKernel(npd.bt.dims, Npts, getType(npd));
  t = [sort(rand(Npts));10];
  ii = 1
  for i in 1:size(pts,2)
    while (w[i] > t[ii])
      points[:,ii] = pts[:,i] + bw[:,i].*randnums[:,ii];
      ind[ii] = i
      ii += 1
    end
  end
  return points, ind
end

function sample(npd::BallTreeDensity, Npts::Int64, ind::Array{Int64,1})
  pts = getPoints(npd)
  points = pts[:,ind] + getBW(npd,ind).*randKernel(npd.bt.dims, length(ind), getType(npd))
  return points, ind
end

function rand(p::BallTreeDensity, N::Int64=1)
    return KernelDensityEstimate.sample(p,N)[1]
end


text:
import Base.*
import Base.size
import Base.eltype
import Base.issym

immutable SqLinOp{Tv,Ti}
    issym::Bool
    value::Tv
    n::Ti
    multFn::Function #this is bad in terms of types?

    SqLinOp(issym,value,n,multFn) = new(issym,value,n,multFn)
end

SqLinOp{Tv,Ti}(issym,value::Tv,index::Ti,multFn) = SqLinOp{Tv,Ti}(issym,value,index,multFn)

eltype{Tv,Ti}(A::SqLinOp{Tv,Ti}) = Tv

size{Tv,Ti}(A::SqLinOp{Tv,Ti}, d::Ti) = A.n

size{Tv,Ti}(A::SqLinOp{Tv,Ti}) = (A.n,A.n)

issym{Tv,Ti}(A::SqLinOp{Tv,Ti}) = A.issym

function *{Tv,Ti}(A::SqLinOp{Tv,Ti}, b::Array{Tv,1})
    return A.multFn(b)
end

function testId{Ti}(n::Ti)
    return M = SqLinOp{Float64,Int64}(true,1.0,n,x -> x)
end


text:
"""
    TypeU(var_1::TypeY)

Lorem ipsum dolor sit amet.
"""
struct TypeU

  var_1::TypeY

end


text:
module Hyperopt

export Hyperoptimizer, @hyperopt, @phyperopt, @thyperopt, printmin, printmax
export RandomSampler, BlueNoiseSampler, LHSampler, CLHSampler, Continuous, Categorical, GPSampler, Max, Min, Hyperband

using Base.Threads: threadid, nthreads
using LinearAlgebra, Statistics, Random
using Juno
using MacroTools
using MacroTools: postwalk, prewalk
using RecipesBase
using Distributed
using LatinHypercubeSampling
using BayesianOptimization, GaussianProcesses
using ThreadPools

const HO_RNG = [MersenneTwister(rand(1:1000)) for _ in 1:nthreads()]

abstract type Sampler end
Base.@kwdef mutable struct Hyperoptimizer{S<:Sampler, F}
    iterations::Int
    params
    candidates
    history = []
    results = []
    sampler::S = RandomSampler()
    objective::F = nothing
end

function Base.getproperty(ho::Hyperoptimizer, s::Symbol)
    s === :minimum && (return minimum(replace(ho.results, NaN => Inf)))
    s === :minimizer && (return ho.history[argmin(ho.results)])
    s === :maximum && (return maximum(replace(ho.results, NaN => Inf)))
    s === :maximizer && (return ho.history[argmax(ho.results)])
    return getfield(ho,s)
end

Base.propertynames(ho::Hyperoptimizer) = (:minimum, :minimizer, :maximum, :maximizer, fieldnames(Hyperoptimizer)...)

include("samplers.jl")

function Hyperoptimizer(iterations::Int, sampler::Sampler = RandomSampler(); kwargs...)
    params = keys(kwargs)
    candidates = []
    for kw in kwargs
        push!(candidates, kw[2])
    end
    Hyperoptimizer(iterations=iterations, params=params, candidates=candidates, sampler=sampler)
end

Base.length(ho::Hyperoptimizer) = ho.iterations

function Base.iterate(ho::Hyperoptimizer, state=1)
    state > ho.iterations && return nothing
    samples = ho.sampler(ho, state)
    push!(ho.history, samples)
    nt = (; Pair.((:i, ho.params...), (state, samples...))...)
    nt, state+1
end

function preprocess_expression(ex)
    ex.head == :for || error("Wrong syntax, Use For-loop syntax, ex: @hyperopt for i=100, param=LinRange(1,10,100) ...")
    params     = []
    candidates = []
    sampler_ = :(RandomSampler())
    ho_ = nothing
    loop = ex.args[1].args # ex.args[1] = the arguments to the For loop
    i = 1
    while i <= length(loop)
        if @capture(loop[i], sampler = sam_) # A sampler was provided
            sampler_ = sam
            deleteat!(loop,i) # Remove the sampler from the args
            continue
        end
        if @capture(loop[i], ho = h_)
            ho_ = h
            deleteat!(loop,i)
            continue
        end
        @capture(loop[i], param_ = list_) || error("Wrong syntax, Use For-loop syntax, ex: @hyperopt for i=100, param=LinRange(1,10,100) ...")
        push!(params, param)
        push!(candidates, list)
        i += 1
    end
    params = ntuple(i->params[i], length(params))
    state_ = sampler_.args[1] === :Hyperband ? :(state = nothing) : :(nothing)
    fun = quote
            $(Expr(:tuple, esc.(params)...)) -> begin
                $(esc(state_))
                $(esc(ex.args[2]))
            end
        end

    params, candidates, sampler_, ho_, fun
end

function optimize(ho::Hyperoptimizer)
    Juno.progress() do id
        for nt = ho
            res = ho.objective(nt...)
            push!(ho.results, res)
            Base.CoreLogging.@logmsg Base.CoreLogging.BelowMinLevel "Hyperopt" progress=nt.i/ho.iterations  _id=id
        end
    end
    ho
end

function create_ho(params, candidates, sampler, ho_, objective_)
    quote
        objective, iters = $(esc(sampler)) isa Hyperband ? $(objective_) : ($(objective_), $(esc(candidates[1])))
        ho = Hyperoptimizer(iterations = iters, params = $(esc(params[2:end])), candidates = $(Expr(:tuple, esc.(candidates[2:end])...)), sampler=$(esc(sampler)), objective = objective)
        if $(esc(ho_)) isa Hyperoptimizer # get info from existing ho. the objective function might be changed, due to variables are captured into the closure, so the type of ho also changed.
            ho.sampler = $(esc(ho_)).sampler
            ho.history = $(esc(ho_)).history # it's important to use the same array, not copy.
            ho.results = $(esc(ho_)).results
        else
            s = (x->x isa Hyperband ? x.inner : x)(ho.sampler)
            s isa Union{LHSampler,CLHSampler,GPSampler} && init!(s, ho)
        end
        ho
    end
end

macro hyperopt(ex)
    pre = preprocess_expression(ex)
    if pre[3].args[1] === :Hyperband
        costfun_ = hyperband_costfun(ex, pre...)
        ho_ = create_ho(pre[1:4]..., costfun_)
        quote
            ho = $ho_
            hyperband(ho)
            ho
        end
    else
        ho_ = create_ho(pre...)
        :(optimize($ho_))
    end
end

function pmacrobody(ex, params, ho_, pmap=pmap)
    quote
        function workaround_function()
            ho = $(ho_)
            # Getting the history right is tricky when using workers. The approach I've found to work is to
            # save the actual array (not copy) in hist, temporarily use a new array that will later be discarded
            # reassign the original array and then append the new history. If a new array is used, the change will not be visible in the original hyperoptimizer
            hist = ho.history
            ho.history = []
            res = $(pmap)(1:ho.iterations) do i
                $(Expr(:tuple, esc.(params)...)),_ = iterate(ho,i)
                res = $(esc(ex.args[2])) # ex.args[2] = Body of the For loop

                res, $(Expr(:tuple, esc.(params[2:end])...))
            end
            ho.history = hist
            append!(ho.results, getindex.(res,1))
            append!(ho.history, getindex.(res,2)) # history automatically appended by the iteration
            ho
        end
        workaround_function()
    end
end

"""
Same as `@hyperopt` but uses `Distributed.pmap` for parallel evaluation of the cost function.
"""
macro phyperopt(ex)
    pre = preprocess_expression(ex)
    pre[3].args[1] === :GPSampler && error("We currently do not support running the GPSampler in parallel. If this is an issue, open an issue ;)")
    ho_ = create_ho(pre...)
    pmacrobody(ex, pre[1], ho_)
end

"""
Same as `@hyperopt` but uses `ThreadPools.tmap` for multithreaded evaluation of the cost function.
"""
macro thyperopt(ex)
    pre = preprocess_expression(ex)
    pre[3].args[1] === :GPSampler && error("We currently do not support running the GPSampler in parallel. If this is an issue, open an issue ;)")
    ho_ = create_ho(pre...)
    pmacrobody(ex, pre[1], ho_, ThreadPools.tmap)
end

function Base.minimum(ho::Hyperoptimizer)
    m,i = findmin(replace(ho.results, NaN => Inf))
    m
end
function Base.maximum(ho::Hyperoptimizer)
    m,i = findmax(replace(ho.results, NaN => -Inf))
    m
end

function printmin(ho::Hyperoptimizer)
    for (param, value) in zip(ho.params, ho.minimizer)
        println(param, " = ", value)
    end
end

function printmax(ho::Hyperoptimizer)
    for (param, value) in zip(ho.params, ho.maximizer)
        println(param, " = ", value)
    end
end

include("plotting.jl")
end # module


text:
using Genie.Router
using DiariesController
using Dates

route("/") do
  serve_static_file("welcome.html")
end

route("/diaries", DiariesController.index, named = :index_diaries)
route("/diaries/show", DiariesController.show, named = :show_diary)
route("/diaries/new", DiariesController.new, named = :new_diary)
route("/diaries/create", DiariesController.create, method = POST, named = :create_diary)
route("/diaries/edit", DiariesController.edit, named = :edit_diary)
route("/diaries/update", DiariesController.update, method = POST, named = :update_diary)
# route("/diaries/destroy", DiariesController.update, method = POST, named = :destroy_diary)

max_stars_repo_name: ky-s/DiaryApp

max_issues_repo_name: ky-s/DiaryApp

max_forks_repo_name: ky-s/DiaryApp


text:
export Transmon,
       Resonator,
       TunableDuffingTransmon,
       TunableFullTransmon,
       TunableTransmon,
       Qubit,
       Duffing

#Resonator
type Resonator <: QSystem
    label::String
    freq::Float64
    dim::Int
end
hamiltonian(r::Resonator) = r.freq*number(r)

#Fixed frequency transmon
type Transmon <: QSystem
    label::String
    E_C::Float64
    E_J::Float64
    dim::Int
end

function Transmon(label::String, nu::Float64, alpha::Float64; dim=3)
    E_C = abs(alpha)
    E_J = (nu + E_C)^2 / (8 * E_C)
    return Transmon(label, E_C, E_J, dim)
end

hamiltonian(t::Transmon) = sqrt(8*t.E_J*t.E_C)*number(t) - t.E_C/12*(X(t)^4)

#Tunable transmon

abstract TunableTransmon <: QSystem

type TunableFullTransmon <: TunableTransmon
    label::String
    E_C::Float64
    E_J::Float64 #sum of junction E_J's
    d::Float64 #asymmetry parameter
    dim::Int
    fluxBias::Float64 # flux bias in units of Phi_0
    flux::Float64 #total flux in units of Phi_0
end
TunableFullTransmon(label::String, E_C::Float64, E_J::Float64, d::Float64, dim::Int, fluxBias::Float64) =
  TunableFullTransmon(label, E_C, E_J, d, dim, fluxBias, fluxBias)

#Helper function to calculate effective EJ for a transmon
scale_EJ(flux::Float64, d::Float64) = cos(pi*flux)*sqrt(1 + d^2*(tan(pi*flux)^2))

function hamiltonian(tt::TunableFullTransmon, t::Float64=0.0)
    myE_J = tt.E_J*scale_EJ(tt.flux, tt.d)
    return (sqrt(8*tt.E_C*myE_J)*number(tt) - (1.0/12)*tt.E_C*(X(tt)^4))
end

#Tunable Duffing approx. transmon
type TunableDuffingTransmon <: TunableTransmon
    label::String
    E_C::Float64
    E_J::Float64 #sum of junction E_J's
    d::Float64 #asymmetry parameter
    dim::Int
    fluxBias::Float64 # flux bias in units of Phi_0
    flux::Float64 #total flux in units of Phi_0
end
TunableDuffingTransmon(label::String, E_C::Float64, E_J::Float64, d::Float64, dim::Int, fluxBias::Float64) = TunableDuffingTransmon(label, E_C, E_J, d, dim, fluxBias, fluxBias)

#Helper function to calculate effective EJ for a transmon
scale_EJ(flux::Float64, d::Float64) = cos(pi*flux)*sqrt(1 + d^2*(tan(pi*flux)^2))

function hamiltonian(tt::TunableDuffingTransmon, t::Float64=0.0)
    myE_J = tt.E_J*scale_EJ(tt.flux, tt.d)
    omega_p = sqrt(8*tt.E_C*myE_J)
    omega = [(omega_p-tt.E_C/2)*i-tt.E_C/2*i^2 for i in 0:(tt.dim-1)]
    return diagm(omega)
end

#Basic two-level qubit
type Qubit <: QSystem
    label::String
    freq::Float64
end
dim(q::Qubit) = 2
hamiltonian(q::Qubit) = q.freq*number(q)

#Duffing oscillator
type Duffing <: QSystem
    label::String
    freq::Float64
    alpha::Float64
    dim::Int
end
hamiltonian(s::Duffing) = (s.freq - 0.5*s.alpha)*number(s) + 0.5*s.alpha * number(s)^2


text:
mutable struct RunningAverageMeter{T<:AbstractFloat}
    momentum::T
    val::Union{Nothing, T}
    avg::T
end

RunningAverageMeter(;momentum::T = 0.99f0) where {T<:AbstractFloat} =
    RunningAverageMeter(momentum, nothing, T(0))

function reset!(r::RunningAverageMeter{T}) where T
    r.val = nothing
    r.avg = T(0)
    return r
end

Zygote.@nograd reset!

function update(r::RunningAverageMeter, val)
    if isnothing(r.val)
        r.avg = val
    else
        r.avg = r.avg * r.momentum + val * (1 - r.momentum)
    end
    r.val = val
    return nothing
end

Zygote.@nograd update

format_string(s) = replace(lowercase(s), " " => "_")


text:
using Base.Test

using BenchmarkTools

using FluxionArrays

expr_type = FluxionArrays.expr_type

@testset "Expression building" begin
    s = 2
    data = rand(10)
    data2 = rand(10, 20)
    x = unknown_array()
    y = unknown_array()
    z = unknown_array()
    a = lazy_array(data)
    b = lazy_array(data2)

    func = x -> x.func

    @test func(x + y) == :+
    @test func(x .+ y) == :.+
    @test func(x + s) == :+
    @test func(x .+ s) == :.+
    @test func(s + x) == :+
    @test func(s .+ x) == :.+
    @test func(x + data) == :+
    @test func(x .+ data) == :.+
    @test func(x + y + z) == :+
    @test func(x .+ y .+ z) == :.+
    @test func(s .+ data2 .+ z) == :.+

    # The calls that involve a literal will be fused
    # into a broadcast of an anonymous function by default.
    # We check that the library builds a proper expression instead.
    @test func(x + 1) == :+
    @test func(x .+ 1) == :.+
    @test func(x + [1 2 3]) == :+
    @test func(x .+ [1 2 3]) == :.+
    @test func(x + y + 1) == :+
    @test func(x .+ y .+ 1) == :.+
    @test func(x + 1 + y) == :+
    @test func(x .+ 1 .+ y) == :.+
    @test func(1 + x + y) == :+
    @test func(1 .+ x .+ y) == :.+
    @test func(z .+ 1 .+ x .+ y) == :.+
    @test func(1 .+ [1 2 3] .+ x) == :.+

end

@testset "Type propagation" begin

    data = rand(10)
    data2 = rand(10, 20)
    x = unknown_array()
    y = unknown_array()
    a = lazy_array(data)
    b = lazy_array(data2)

    @test expr_type(x .+ y) == FluxionArrays.ArrayType(nothing, nothing)

    x = unknown_array(shape=[10, 20, 30])
    y = unknown_array(shape=[10, 20])

    @test size(permutedims(reshape(x .+ y .+ 1, 100, 60), [2, 1])) == (60, 100)
end

@testset "Timings" begin

    data = rand(1000000)
    data2 = rand(1000000)
    x = unknown_array()
    y = unknown_array()
    a = lazy_array(data)
    b = lazy_array(data2)

    x_data = rand(1000000)
    y_data = rand(1000000)

    c = a .+ 1 .+ x
    d = 2 .+ b .+ y
    e = c .+ d .+ size(x, 1) .+ reshape(x, size(x, 1))

    function ref(data, data2, x, y)
        (data .+ 1 .+ x) .+ (2 .+ data2 .+ y) .+ size(x, 1) .+ reshape(x, size(x, 1))
    end

    println("Timing compilation of an array expression")
    display(@benchmark compile_array_function(array_function([$e], [$x, $y])))
    println()
    f = compile_array_function(array_function([e], [x, y]))

    println("Timing the reference function")
    display(@benchmark $ref($data, $data2, $x_data, $y_data))
    println()

    println("Timing the compiled function")
    display(@benchmark $f($x_data, $y_data))
    println()

    res1 = from_backend(f(x_data, y_data))
    res2 = ref(data, data2, x_data, y_data)

    @test isapprox(res1, res2)

end

@testset "Multiple returns" begin
    data = rand(100)
    data2 = rand(100)
    x = unknown_array()
    y = unknown_array()
    a = lazy_array(data)
    b = lazy_array(data2)

    x_data = rand(100)
    y_data = rand(100)

    f = compile_array_function(array_function([a + 1 + x, 2 + b + y], [x, y]))
    res1, res2 = f(x_data, y_data)

    @test isapprox(from_backend(res1), data .+ 1 .+ x_data)
    @test isapprox(from_backend(res2), 2 .+ data2 .+ y_data)
end

@testset "Mutations" begin

    data = rand(100)
    data2 = rand(100)
    x = unknown_array()
    y = unknown_array()
    a = lazy_array(data)
    b = lazy_array(data2)

    x_data = rand(100)
    x_data_copy = copy(x_data)
    y_data = rand(100)

    x_data_backend = to_backend(x_data)

    f = compile_array_function(array_function([a + 1 + x], [x, y], [2 + b + y => x]))
    res1 = from_backend(f(x_data_backend, y_data))

    @test isapprox(res1, data .+ 1 .+ x_data_copy)
    @test isapprox(from_backend(x_data_backend), 2 .+ data2 .+ y_data)

    # Test that an otherwise unused parameter can be mutated
    f = compile_array_function(array_function([a + 1], [x, y], [3 + b + y => x]))
    res1 = from_backend(f(x_data_backend, y_data))

    @test isapprox(res1, data .+ 1)
    @test isapprox(from_backend(x_data_backend), 3 .+ data2 .+ y_data)

end

@testset "Randoms" begin

    r1 = unknown_random_state()
    r2 = unknown_random_state()

    shape1 = (10, 100)
    shape2 = (20, 50)

    n1, r1_new = random_normal(r1, 1, 2, shape1...)
    n2, r2_new = random_uniform(r2, 0, 10, shape1...)
    n3, r2_new = random_uniform(r2_new, 0, 11, shape2...)
    n4, r1_new = random_normal(r1_new, 1, 3, shape2...)

    # Also useful: mutating random
    # n3 = random!(r2_new, shape2..., distribution=:uniform, min=0, max=10)
    # r2_new is replaced by a new RandomState() object, while the old one is preserved.
    # can we actually do that?

    expr1 = n1 + n2
    expr2 = n3 + n4

    # No updates
    f = compile_array_function(array_function([expr1, expr2, r1_new, r2_new], [r1, r2]))
    rs1 = random_state(seed=123)
    rs2 = random_state(seed=456)
    arr1, arr2, new_rs1, new_rs2 = f(rs1, rs2)

    # With explicit updates (if the random state is very large and hard to copy, e.g. for Random123)
    f = compile_array_function(array_function([expr1, expr2], [r1, r2], [r1_new => r1, r2_new => r2]))
    rs1 = random_state(seed=123)
    rs2 = random_state(seed=456)
    arr1, arr2 = f(rs1, rs2) # requires rng1 and rng2 to be mutable!

end

@testset "Nested functions" begin

    data = rand(100)
    data2 = rand(100)
    x = unknown_array()
    y = unknown_array()
    a = lazy_array(data)
    b = lazy_array(data2)

    x_data = rand(100)
    x_data_copy = copy(x_data)
    y_data = rand(100)

    f1 = array_function([a + 1 + x, 3 + y], [x, y])
    t1, t2 = f1(y, x)
    f2 = array_function([t1 + x, t2 - y], [x, y])

    # reference
    t1 = data .+ 1 .+ y_data
    t2 = 3 .+ x_data
    ref1 = t1 .+ x_data
    ref2 = t2 .- y_data

    f2c = compile_array_function(f2)
    res1, res2 = f2c(x_data, y_data)

    @test isapprox(from_backend(res1), ref1)
    @test isapprox(from_backend(res2), ref2)

end


text:
"""
A simple `AbstractResultTable` that uses a .csv file as a backend.

    CSVRecord(filename::String, headers::Union{Nothing,Vector{String}}=nothing)

Creates a `CSVRecord` table which writes the results to a .csv file at the
location `filename` in the file system. Optionally, pass `headers` specifying
the names of the fields to be recorded and the first line of the created .csv
will include this header information.
"""
mutable struct CSVRecord <: AbstractResultTable
    fp::IOStream

    function CSVRecord(
        filename::String,
        headers::Union{Nothing,Vector{String}} = nothing,
    )
        fp = open(filename, "w+")
        if headers !== nothing
            println(fp, join(headers, ","))
            flush(fp)
        end
        record = new(fp)
        finalizer(r -> close(r.fp), record)
        return record
    end
end


text:
struct SubAxis{G,I,Inds<:AbstractRange{I}} <: AbstractAxis{I,Inds}
    getindices::G
    parent::Inds
end

_getindices(axis::SubAxis) = getfield(axis, :getindices)
Base.first(::SubAxis) = One()
Base.last(axis::SubAxis) = static_length(_getindices(axis))
Base.eachindex(axis::SubAxis) = first(axis):last(axis)

@propagate_inbounds function Base.view(axis::AbstractAxis, arg::AbstractUnitRange{<:Integer})
    return SubAxis(to_index(axis, arg), axis)
end
@propagate_inbounds function Base.view(axis::AbstractAxis, arg::StepRange{<:Integer})
    return SubAxis(to_index(axis, arg), axis)
end
@propagate_inbounds function Base.view(axis::AbstractAxis, arg)
    return SubAxis(to_index(axis, arg), axis)
end

_sub_offset(axis::SubAxis, arg) = @inbounds(_getindices(axis)[arg])


text:
#
# Quasi Submodule for Inference
#

include("exact.jl")
include("likelihood.jl")
include("gibbs.jl")
include("loopy_belief.jl")


text:
include("lzz_p.jl")
import Base: setindex!, getindex, zero, one, copy,>>, <<

## vector of zz_p{T}
const vec_zz_p{T} = Vector{zz_p{T}} where {T}

function convert(::Type{vec_zz_p{T}}, v::Vector{Int}) where {T}
    return map(zz_p{T}, v)
end

## polynomial with coefficients in zz_p{T}
struct NTL_INIT_zz_pX end
ntl_init_zz_pX = NTL_INIT_zz_pX()

mutable struct zz_pX{T}
    _rep::vec_zz_p{T}
    function zz_pX{T}(X::Int) where {T}
        x = zz_p{T}(X)
        new([x])
    end
    function zz_pX{T}() where {T}
        return zz_pX{T}(0)
    end
    function zz_pX{T}(x::zz_p{T})  where {T}
        new([x])
    end
    function zz_pX{T}(v::vec_zz_p{T}) where {T}
        p=new(v)
        normalize!(p)
    end
    function zz_pX{T}(v::Array{Int64,1}) where {T}
        p=new(v)
        normalize!(p)
    end
    function zz_pX{T}(v::Array{Int64,1},::NTL_INIT_zz_pX) where {T}
        new(v)
    end
    function zz_pX{T}(v::vec_zz_p{T},::NTL_INIT_zz_pX) where {T}
        new(v)
    end
end

convert(::Type{zz_pX{T}}, X::Int) where {T} = zz_pX{T}(X)
convert(::Type{zz_pX{T}}, X::zz_p{T}) where {T} = zz_pX{T}(X)

"""
show(io::IO, p::zz_pX{T}, x::Symbol)

show the polynomial p in variable x
julia> p = zz_pX{7}([1,2,3])
julia> show(stdout, p, :x)
"""
function show(io::IO, p::zz_pX{T}, sym::Symbol) where {T}
    coef = p._rep
    print(io, coef[1])
    x = string(sym)
    len = length(p._rep)
    for i = 2:len
        if !iszero(coef[i])
        print(io, "+", coef[i], "*", x, "^", i-1)
        end
    end
    return nothing
end
function normalize!(p::zz_pX{T})::zz_pX{T} where {T}
    r = p._rep
    m= length(r)
    while (m > 1) && iszero(r[m])
        pop!(r)
        m -= 1
    end
    return p
end

@inline function deg(p::zz_pX{T}) where {T}
    return length(p._rep)-1
end
function coeff(p::zz_pX{T}, i::Int) where {T}
    if i < 0
        error("coeff: negative index")
    end
    if i in 0:deg(p)
        return p._rep[i+1]
    end
    return zz_p{T}(0)
end

function zero(::Type{zz_pX{T}}) where {T}
    return zz_pX{T}(0)
end
function one(::Type{zz_pX{T}}) where {T}
    return zz_pX{T}(1)
end
function setX(::Type{zz_pX{T}}) where {T}
    zz_pX{T}([0,1],ntl_init_zz_pX)
end

@inline function isX(x::zz_pX{T}) where {T}
   return (deg(x) == 1) && isone(x._rep[2]) && iszero(x._rep[1]);
end
@inline function iszero(x::zz_pX{T}) where {T}
    return (deg(x)== 0) && iszero(x._rep[1])
end
@inline function isone(x::zz_pX{T}) where {T}
    return (deg(x)== 0) && isone(x._rep[1])
end

function GetCoeff!(x::zz_p{T}, p::zz_pX{T}, i::Int) where {T}
    x._rep = coeff(p,i)._rep
end
function LeadCoeff(p::zz_pX{T}) where {T}
    p._rep[end]
end
function ConstTerm(p::zz_pX{T}) where {T}
    p._rep[1]
end
function getindex(p::zz_pX{T}, i::Int) where {T}
    if i < 0
        error("getindex: negative index")
    end
    if i in 0:deg(p)
        return getindex(p._rep, i+1)
    end
    return zero(zz_pX{T})
end

function setindex!(p::zz_pX{T},  x::zz_p{T}, i::Int) where {T}
    if i < 0
        error("setindex!: nagative index")
    end
    m = deg(p)
    if (i > m) && (iszero(x))
        return x
    end
    if i>m
        z=zeros(zz_p{T},i-m)
        append!(p._rep, z)
    end
    setindex!(p._rep, x, i+1)
    normalize!(p)
end

setindex!(p::zz_pX{T},  x::Int, i::Int) where {T} = setindex!(p,zz_p{T}(x),i)

@inline function swap!(x::zz_pX{T}, y::zz_pX{T}) where {T}
    z = x._rep; x._rep = y._rep; y._rep = z
    nothing
end

"""
rand(::Type, d::Int)

generate a random polynomial of degree <= d

## Example
julia> P7=zz_pX{7}
julia> rand(P7, 3)

"""
function rand(::Type{zz_pX{T}}, d::Int) where {T}
    if d < 0
        error("rand: nagative degree")
    end
    zz_pX{T}(rand(Int, d+1))
end

function trunc(p::zz_pX{T}, d::Int) where {T}
    if d > deg(p)
        return p
    else
        zz_pX{T}(p._rep[1:d])
    end
end
function trunc!(p::zz_pX{T}, d::Int) where {T}
    if d > deg(p)
        return p
    else
        resize!(p._rep,d)
    end
    normalize!(p)
end
"""
RightShift!(p::zz_pX{T}, d::Int)
p /= x^d

"""
@inline function RightShift!(p::zz_pX{T}, d::Int) where {T}
    if d>deg(p)
        p._rep = [ zero(zz_p{T}) ]
    else
        deleteat!(p._rep, 1:d)
    end
    return p
end
@inline function RightShift(p::zz_pX{T}, d::Int) where {T}
    if d>deg(p)
        return zero(zz_pX{T})
    else
        return zz_pX{T}(p._rep[1+d:end],ntl_init_zz_pX)
    end
 end
 >>(p::zz_pX{T}, d::Int) where {T} = RightShift(p,d)

"""
LeftShift!(p::zz_pX{T}, d::Int)
p *= x^d

"""
@inline function LeftShift!(p::zz_pX{T}, d::Int) where {T}
    if (d>0) && !(iszero(p))
        prepend!(p._rep,zeros(zz_p{T}, d))
    end
    return p
end
@inline function LeftShift(p::zz_pX{T}, d::Int) where {T}
    LeftShift!(copy(p),d)
end
<<(p::zz_pX{T}, d::Int) where {T} = LeftShift(p,d)

function copy(p::zz_pX{T}) where {T}
    zz_pX{T}(copy(p._rep), ntl_init_zz_pX)
end

==(a::zz_pX{T}, b::zz_pX{T}) where {T} = a._rep == b._rep

@inline function diff(a::zz_pX{T}) where {T}
    b=a>>1
    b._rep.*=collect(1:length(b._rep))
    normalize!(b)
end

# ***************************************************************
#                          Addition
# ***************************************************************
@inline function add(x::zz_pX{T}, y::zz_pX{T}) where {T}
    dx =length(x._rep); dy = length(y._rep)
    if dx>dy
        _add_large(x,dx,y, dy)
    elseif dx<dy
        _add_large(y,dy,x, dx)
    else
        _add_equal(x,y)
    end
end
@inline function _add_equal(x::zz_pX{T}, y::zz_pX{T}) where{T}
    zz_pX{T}(x._rep + y._rep)
end
@inline function _add_large(x::zz_pX{T} , dx::Int, y::zz_pX{T}, dy::Int ) where {T}
    z=append!(x._rep[1:dy]+y._rep, x._rep[dy+1:dx])
    zz_pX{T}(z, ntl_init_zz_pX)
end
+(x::zz_pX{T},y::zz_pX{T}) where {T} = add(x, y)
add(X::zz_p{T},y::zz_pX{T}) where {T} = add(convert(zz_pX{T},X), y)
add(X::Int,y::zz_pX{T}) where {T} = add(convert(zz_pX{T},X), y)
add(x::zz_pX{T},Y::zz_p{T}) where {T} = add(x,convert(zz_pX{T},Y))
add(x::zz_pX{T},Y::Int) where {T} = add(x,convert(zz_pX{T},Y))


text:
#    Modified Bessel functions of the first kind of order zero and one
#                       besseli0, besseli1
#
#    Scaled modified Bessel functions of the first kind of order zero and one
#                       besseli0x, besselix
#
#    (Scaled) Modified Bessel functions of the first kind of order nu
#                       besseli, besselix
#
#    Calculation of besseli0 is done in two branches using polynomial approximations [1]
#
#    Branch 1: x < 7.75
#              besseli0 = [x/2]^2 P16([x/2]^2)
#    Branch 2: x >= 7.75
#              sqrt(x) * exp(-x) * besseli0(x) = P22(1/x)
#    where P16 and P22 are a 16 and 22 degree polynomial respectively.
#
#    Remez.jl is then used to approximate the polynomial coefficients of
#    P22(y) = sqrt(1/y) * exp(-inv(y)) * besseli0(inv(y))
#    N,D,E,X = ratfn_minimax(g, (1/1e6, 1/7.75), 21, 0)
#
#    A third branch is used for scaled functions for large values
#
#
#    Calculation of besseli1 is done in two branches using polynomial approximations [2]
#
#    Branch 1: x < 7.75
#              besseli1 = x / 2 * (1 + 1/2 * (x/2)^2 + (x/2)^4 * P13([x/2]^2)
#    Branch 2: x >= 7.75
#              sqrt(x) * exp(-x) * besseli1(x) = P22(1/x)
#    where P13 and P22 are a 16 and 22 degree polynomial respectively.
#
#    Remez.jl is then used to approximate the polynomial coefficients of
#    P13(y) = (besseli1(2 * sqrt(y)) / sqrt(y) - 1 - y/2) / y^2
#    N,D,E,X = ratfn_minimax(g, (1/1e6, 1/7.75), 21, 0)
#
#    A third branch is used for scaled functions for large values
#
#    Horner's scheme is then used to evaluate all polynomials.
#    ArbNumerics.jl is used as the reference bessel implementations with 75 digits.
#
#    Calculation of besseli and besselkx can be done with downward recursion starting with
#    besseli_{nu+1} and besseli_{nu}. Higher orders are determined by a uniform asymptotic
#    expansion similar to besselk (see notes there) using Equation 10.41.3 [3].
#
#
# [1] "Rational Approximations for the Modified Bessel Function of the First Kind
#     - I0(x) for Computations with Double Precision" by Pavel Holoborodko
# [2] "Rational Approximations for the Modified Bessel Function of the First Kind
#     - I1(x) for Computations with Double Precision" by Pavel Holoborodko
# [3] https://dlmf.nist.gov/10.41

"""
    besseli0(x::T) where T <: Union{Float32, Float64}

Modified Bessel function of the first kind of order zero, ``I_0(x)``.
"""
function besseli0(x::T) where T <: Union{Float32, Float64}
    x = abs(x)
    if x < 7.75
        a = x * x / 4
        return muladd(a, evalpoly(a, besseli0_small_coefs(T)), 1)
    else
        a = exp(x / 2)
        s = a * evalpoly(inv(x), besseli0_med_coefs(T)) / sqrt(x)
        return a * s
    end
end

"""
    besseli0x(x::T) where T <: Union{Float32, Float64}

Scaled modified Bessel function of the first kind of order zero, ``I_0(x)*e^{-x}``.
"""
function besseli0x(x::T) where T <: Union{Float32, Float64}
    T == Float32 ? branch = 50 : branch = 500
    x = abs(x)
    if x < 7.75
        a = x * x / 4
        return muladd(a, evalpoly(a, besseli0_small_coefs(T)), 1) * exp(-x)
    elseif x < branch
        return evalpoly(inv(x), besseli0_med_coefs(T)) / sqrt(x)
    else
        return evalpoly(inv(x), besseli0_large_coefs(T)) / sqrt(x)
    end
end

"""
    besseli1(x::T) where T <: Union{Float32, Float64}

Modified Bessel function of the first kind of order one, ``I_1(x)``.
"""
function besseli1(x::T) where T <: Union{Float32, Float64}
    z = abs(x)
    if z < 7.75
        a = z * z / 4
        inner = (one(T), T(0.5), evalpoly(a, besseli1_small_coefs(T)))
        z = z * evalpoly(a, inner) / 2
    else
        a = exp(z / 2)
        s = a * evalpoly(inv(z), besseli1_med_coefs(T)) / sqrt(z)
        z =  a * s
    end
    if x < zero(x)
        z = -z
    end
    return z
end

"""
    besseli1x(x::T) where T <: Union{Float32, Float64}

Scaled modified Bessel function of the first kind of order one, ``I_1(x)*e^{-x}``.
"""
function besseli1x(x::T) where T <: Union{Float32, Float64}
    T == Float32 ? branch = 50 : branch = 500
    z = abs(x)
    if z < 7.75
        a = z * z / 4
        inner = (one(T), T(0.5), evalpoly(a, besseli1_small_coefs(T)))
        z = z * evalpoly(a, inner) / 2 * exp(-z)
    elseif z < branch
        z = evalpoly(inv(z), besseli1_med_coefs(T)) / sqrt(z)
    else
        z = evalpoly(inv(z), besseli1_large_coefs(T)) / sqrt(z)
    end
    if x < zero(x)
        z = -z
    end
    return z
end

"""
    besseli(nu, x::T) where T <: Union{Float32, Float64}

Modified Bessel function of the first kind of order nu, ``I_{nu}(x)``.
Nu must be real.
"""
function besseli(nu, x::T) where T <: Union{Float32, Float64}
    nu == 0 && return besseli0(x)
    nu == 1 && return besseli1(x)

    if x > maximum((T(30), nu^2 / 4))
        return T(besseli_large_argument(nu, x))
    elseif x <= 2 * sqrt(nu + 1)
        return T(besseli_small_arguments(nu, x))
    elseif nu < 100
        return T(_besseli_continued_fractions(nu, x))
    else
        return T(besseli_large_orders(nu, x))
    end
end

"""
    besselix(nu, x::T) where T <: Union{Float32, Float64}

Scaled modified Bessel function of the first kind of order nu, ``I_{nu}(x)*e^{-x}``.
Nu must be real.
"""
function besselix(nu, x::T) where T <: Union{Float32, Float64}
    nu == 0 && return besseli0x(x)
    nu == 1 && return besseli1x(x)

    if x > maximum((T(30), nu^2 / 4))
        return T(besseli_large_argument_scaled(nu, x))
    elseif x <= 2 * sqrt(nu + 1)
        return T(besseli_small_arguments(nu, x)) * exp(-x)
    elseif nu < 100
        return T(_besseli_continued_fractions_scaled(nu, x))
    else
        return besseli_large_orders_scaled(nu, x)
    end
end
function besseli_large_orders(v, x::T) where T <: Union{Float32, Float64}
    S = promote_type(T, Float64)
    x = S(x)
    z = x / v
    zs = hypot(1, z)
    n = zs + log(z) - log1p(zs)
    coef = SQ1O2PI(S) * sqrt(inv(v)) * exp(v*n) / sqrt(zs)
    p = inv(zs)
    p2  = v^2/fma(max(v,x), max(v,x), min(v,x)^2)

    return coef*Uk_poly_In(p, v, p2, T)
end
function besseli_large_orders_scaled(v, x::T) where T <: Union{Float32, Float64}
    S = promote_type(T, Float64)
    x = S(x)
    z = x / v
    zs = hypot(1, z)
    n = zs + log(z) - log1p(zs)
    coef = SQ1O2PI(S) * sqrt(inv(v)) * exp(v*n - x) / sqrt(zs)
    p = inv(zs)
    p2  = v^2/fma(max(v,x), max(v,x), min(v,x)^2)

    return T(coef*Uk_poly_In(p, v, p2, T))
end
function _besseli_continued_fractions(nu, x::T) where T
    S = promote_type(T, Float64)
    xx = S(x)
    knu, knum1 = up_recurrence(xx, besselk0(xx), besselk1(xx), nu)
    # if knu or knum1 is zero then besseli will likely overflow
    (iszero(knu) || iszero(knum1)) && return throw(DomainError(x, "Overflow error"))
    return 1 / (x * (knum1 + knu / steed(nu, x)))
end
function _besseli_continued_fractions_scaled(nu, x::T) where T
    S = promote_type(T, Float64)
    xx = S(x)
    knu, knum1 = up_recurrence(xx, besselk0x(xx), besselk1x(xx), nu)
    # if knu or knum1 is zero then besseli will likely overflow
    (iszero(knu) || iszero(knum1)) && return throw(DomainError(x, "Overflow error"))
    return 1 / (x * (knum1 + knu / steed(nu, x)))
end
function steed(n, x::T) where T
    MaxIter = 1000
    xinv = inv(x)
    xinv2 = 2 * xinv
    d = x / (n + n)
    a = d
    h = a
    b = muladd(2, n, 2) * xinv
    for _ in 1:MaxIter
        d = inv(b + d)
        a *= muladd(b, d, -1)
        h = h + a
        b = b + xinv2
        abs(a / h) <= eps(T) && break
    end
    return h
end
function besseli_large_argument(v, z::T) where T
    MaxIter = 1000
    a = exp(z / 2)
    coef = a / sqrt(2 * T(pi) * z)
    fv2 = 4 * v^2
    term = one(T)
    res = term
    s = -term
    for i in 1:MaxIter
        i = T(i)
        offset = muladd(2, i, -1)
        term *= T(0.125) * muladd(offset, -offset, fv2) / (z * i)
        res = muladd(term, s, res)
        s = -s
        abs(term) <= eps(T) && break
    end
    return res * coef * a
end
function besseli_large_argument_scaled(v, z::T) where T
    MaxIter = 1000
    coef = inv(sqrt(2 * T(pi) * z))
    fv2 = 4 * v^2
    term = one(T)
    res = term
    s = -term
    for i in 1:MaxIter
        i = T(i)
        offset = muladd(2, i, -1)
        term *= T(0.125) * muladd(offset, -offset, fv2) / (z * i)
        res = muladd(term, s, res)
        s = -s
        abs(term) <= eps(T) && break
    end
    return res * coef
end

function besseli_small_arguments(v, z::T) where T
    S = promote_type(T, Float64)
    x = S(z)
    if v < 20
        coef = (x / 2)^v / factorial(v)
    else
        vinv = inv(v)
        coef = sqrt(vinv / (2 * π)) * MathConstants.e^(v * (log(x / (2 * v)) + 1))
        coef *= evalpoly(vinv, (1, -1/12, 1/288,  139/51840, -571/2488320, -163879/209018880, 5246819/75246796800, 534703531/902961561600))
    end

    MaxIter = 1000
    out = one(S)
    zz = x^2 / 4
    a = one(S)
    for k in 1:MaxIter
        a *= zz / (k * (k + v))
        out += a
        a <= eps(T) && break
    end
    return coef * out
end


text:
# convert uploaded graph in DB to MM-iSAM ready form

using Caesar, CloudGraphs

# Uncomment out for command line operation
# cloudGraph, addrdict = standardcloudgraphsetup(nparticles=true, clearslamindb=true)

# interactive operation
include(joinpath(dirname(@__FILE__),"blandauthremote.jl"))
addrdict["session"] = "SESSSHARK_16_11_14"  # "SESSROX"
cloudGraph, addrdict = standardcloudgraphsetup(addrdict=addrdict)
addrdict["clearslamindb"] = "n"

session = addrdict["session"]
segment = ""
@show clearslamindbdata = addrdict["clearslamindb"]=="yes"

info("remember to set n.ready=1")
if !clearslamindbdata
  println("converting front end data")
  fg = initfg(sessionname=session, cloudgraph=cloudGraph)
  Nparticles = parse(Int, addrdict["num particles"])
  updatenewverts!(fg, N=Nparticles)
else
  println("Clearing slamindb data, leaving front-end data, session: $(session)")
  resetentireremotesession(cloudGraph.neo4j.connection,session, segment=segment)
end

println("Done.")

# Juno.breakpoint("/home/dehann/.julia/v0.5/Caesar/src/cloudgraphs/CloudGraphIntegration.jl", 519 )

# Juno.breakpoint("/home/dehann/.julia/v0.5/IncrementalInference/src/FactorGraph01.jl",353)

# dd = getnewvertdict(fg.cg.neo4j.connection, session)

# using Neo4j
#
# conn = fg.cg.neo4j.connection
# loadtx = transaction(conn)
# query = "match (n:$(session))-[:DEPENDENCE]-(f:NEWDATA:$(session):FACTOR) where n.ready=1 or f.ready=1 return distinct n, f"
# cph = loadtx(query, submit=true)


text:
export NonlinearSystem

struct NLEq
    rhs::Expression
end
function Base.convert(::Type{NLEq}, eq::Equation)
    isequal(eq.lhs, Constant(0)) || return NLEq(eq.rhs - eq.lhs)
    return NLEq(eq.rhs)
end
Base.:(==)(a::NLEq, b::NLEq) = a.rhs == b.rhs

"""
$(TYPEDEF)

A nonlinear system of equations.

# Fields
* `eqs` - Vector of equations defining the system.

# Examples

```
@variables x y z
@parameters σ ρ β

eqs = [0 ~ σ*(y-x),
       0 ~ x*(ρ-z)-y,
       0 ~ x*y - β*z]
ns = NonlinearSystem(eqs, [x,y,z])
```
"""
struct NonlinearSystem <: AbstractSystem
    """Vector of equations defining the system."""
    eqs::Vector{NLEq}
    """Unknown variables."""
    vs::Vector{Expression}
    """Parameters."""
    ps::Vector{Variable}
    function NonlinearSystem(eqs, vs)
        rhss = [eq.rhs for eq ∈ eqs]
        ps = reduce(∪, map(_find_params(vs), rhss); init = vnil())
        new(eqs, vs, collect(ps))
    end
end

vnil() = Set{Variable}()
_find_params(vs) = Base.Fix2(_find_params, vs)
function _find_params(O, vs)
    isa(O, Operation) || return vnil()
    any(isequal(O), vs) && return vnil()
    ps = reduce(∪, map(_find_params(vs), O.args); init = vnil())
    isa(O.op, Variable) && push!(ps, O.op)
    return ps
end

independent_variables(sys::NonlinearSystem) = Set{Variable}()
dependent_variables(sys::NonlinearSystem) = Set{Expression}(sys.vs)
parameters(sys::NonlinearSystem) = Set{Variable}(sys.ps)

function calculate_jacobian(sys::NonlinearSystem)
    rhs = [eq.rhs for eq in sys.eqs]
    jac = expand_derivatives.(calculate_jacobian(rhs, sys.vs))
    return jac
end

function generate_jacobian(sys::NonlinearSystem)
    jac = calculate_jacobian(sys)
    return build_function(jac, clean.(sys.vs), sys.ps, (), NLSysToExpr(sys))
end

struct NLSysToExpr
    sys::NonlinearSystem
end
function (f::NLSysToExpr)(O::Operation)
    any(isequal(O), f.sys.vs) && return O.op.name  # variables
    if isa(O.op, Variable)
        isempty(O.args) && return O.op.name  # 0-ary parameters
        return build_expr(:call, Any[O.op.name; f.(O.args)])
    end
    return build_expr(:call, Any[O.op; f.(O.args)])
end
(f::NLSysToExpr)(x) = convert(Expr, x)

function generate_function(sys::NonlinearSystem, vs, ps; version = nothing)
    rhss = [eq.rhs for eq ∈ sys.eqs]
    vs′ = [clean(v) for v ∈ vs]
    ps′ = [clean(p) for p ∈ ps]
    return build_function(rhss, vs′, ps′, (), NLSysToExpr(sys))
end


text:
using Base.Test,
      Cliffords,
      Distributions,
      QuantumInfo,
      QuantumTomography,
      RandomQuantum,
      SchattenNorms,
      SCS

function qst_test_setup()
    obs = Matrix[ (complex(Pauli(i))+eye(2))/2 for i in 1:3 ]
    append!(obs, Matrix[ (-complex(Pauli(i))+eye(2))/2 for i in 1:3 ])

    ρ = .98*projector(rand(FubiniStudyPureState(2)))+.02*rand(HilbertSchmidtMixedState(2))

    return ρ, obs
end

function test_qst_freels(ρ, obs; n=10_000, asymptotic=false)
    tomo = FreeLSStateTomo(obs)

    asymptotic_means = real(predict(tomo,ρ))

    samples = asymptotic ? asymptotic_means : [ rand(Binomial(n,μ))/n for μ in asymptotic_means ]
    sample_mean = samples
    sample_var  = n*(samples - samples.^2)/(n-1)

    ρest, obj, status = fit(tomo, sample_mean, asymptotic ? ones(length(samples)) : sample_var);

    return status, trnorm(ρ-ρest), obj, ρest
end

function test_qst_freels_gls(ρ, obs; n=10_000, asymptotic=false)
    tomo = FreeLSStateTomo(obs)

    asymptotic_means = real(predict(tomo,ρ))

    samples = asymptotic ? asymptotic_means : [ rand(Binomial(n,μ))/n for μ in asymptotic_means ]
    sample_mean = samples
    sample_var  = n*(samples - samples.^2)/(n-1)

    ρest, obj, status = fit(tomo, sample_mean, asymptotic ? ones(length(samples)) : sample_var, algorithm=:GLS);

    return status, trnorm(ρ-ρest), obj, ρest
end

function test_qst_ls(ρ, obs; n=10_000, asymptotic=false)
    tomo = LSStateTomo(obs)

    asymptotic_means = real(predict(tomo,ρ))

    samples = asymptotic ? asymptotic_means : [ rand(Binomial(n,μ))/n for μ in asymptotic_means ]
    sample_mean = samples
    sample_var  = n*(samples - samples.^2)/(n-1)

    ρest, obj, status = fit(tomo, sample_mean, asymptotic ? ones(length(samples)) : sample_var);

    return status, trnorm(ρ-ρest), obj, ρest
end

function test_qst_ml(ρ, obs; n=10_000, β=0.0, asymptotic=false, alt=false, maxiter=5000, ϵ=1000)
    tomo = MLStateTomo(obs,β)

    asymptotic_means = real(predict(tomo,ρ))

    samples = asymptotic ? asymptotic_means[1:3] : Float64[rand(Binomial(n,μ))/n for μ in asymptotic_means[1:3]]
    append!(samples,1-samples)

    ρest, obj, status = fit(tomo, samples, maxiter=maxiter, δ=1/ϵ, λ=β)

    return status, trnorm(ρ-ρest), obj, ρest
end

function test_qst_hml(ρ, obs; n=10_000, β=0.0, asymptotic=false, alt=false, maxiter=5000, ϵ=1000)
    tomo = MLStateTomo(obs,β)

    asymptotic_means = real(predict(tomo,ρ))

    samples = asymptotic ? asymptotic_means[1:3] : Float64[rand(Binomial(n,μ))/n for μ in asymptotic_means[1:3]]
    append!(samples,1-samples)

    ρest, obj, status = fitB(tomo, samples)

    return status, trnorm(ρ-ρest), obj, ρest
end

function test_qpt_ml(n=1000;ρ=zeros(Float64,0,0))

    obs  = Matrix[ [0 1; 1 0], [0 -1im; 1im 0], pauli(3) ]
    prep = map(projector,Vector[ [1;0], [0;1], 1/sqrt(2)*[1;1], 1/sqrt(2)*[1;-1], 1/sqrt(2)*[1;-1im], 1/sqrt(2)*[1;1im] ] )

    exps = Matrix[ choi_liou_involution(vec(o)*vec(p)') for o in obs, p in prep ]

    pred = reduce(vcat, map(m->vec(m)', vec(exps)) )

    if trace(ρ) == 0.0
        ρ = choi_liou_involution(liou(rand_unitary(2)))
    end

    asymptotic_means = real(pred*vec(ρ))

    samples = [ rand(Bernoulli((μ+1)/2),n) for μ in asymptotic_means ]
    sample_mean = 2*map(mean,samples)-1
    sample_var  = 4*map(var,samples)/n

    ρest, obj, status = qpt_ml(pred, sample_mean, sample_var);

    # get the normalization right
    choi_err = ρ-ρest |> choi_liou_involution |> liou2choi

    println("Status                  : $(status)")
    println("Diamond norm lower bound: $(snorm(choi_err,1))")
    println("χ² error                : $(obj)")
    println("Eigvals ρ:")
    for ev in eigvals(ρ)
        println(ev)
    end
    println("Eigvals ρest:")
    for ev in eigvals(ρest)
        println(real(ev))
    end

    return status, snorm(choi_err,1), obj, abs(choi_err), eigvals(ρ), eigvals(ρest)
end

function test_trb(da,db)
    r = rand_mixed_state(da*db)
    norm(trace(r,[da,db],2)-mat(trb_sop(da,db)*vec(r)),1)
end

ρ, obs = qst_test_setup()

kmax = 100
result = zeros(kmax,6)

for k = 1:kmax

    ρ = .98*projector(rand(FubiniStudyPureState(2)))+.02*rand(HilbertSchmidtMixedState(2))

    status, enorm, obj, ρest = test_qst_freels(ρ, obs, asymptotic=true)
    result[k,1] = enorm
    # println("Constrained LS with ∞ counts:")
    # println("   status    : $status")
    # println("   error     : $enorm")
    # println("   true state: $ρ")
    # println("   estimate  : $ρest")
    @test status == :Optimal
    @test enorm < 1e-6

    status, enorm, _, ρest = test_qst_freels(ρ, obs, asymptotic=false)
    result[k,2] = enorm
    # println("Constrained LS with 10_000 counts:")
    # println("   status    : $status")
    # println("   error     : $enorm")
    # println("   true state: $ρ")
    # println("   estimate  : $ρest")
    @test status == :Optimal
    @test enorm < 5e-2

    status, enorm, obj, ρest = test_qst_ls(ρ, obs, asymptotic=true)
    result[k,1] = enorm
    # println("Constrained LS with ∞ counts:")
    # println("   status    : $status")
    # println("   error     : $enorm")
    # println("   true state: $ρ")
    # println("   estimate  : $ρest")
    @test status == :Optimal
    @test enorm < 1e-6

    status, enorm, _, ρest = test_qst_ls(ρ, obs, asymptotic=false)
    result[k,2] = enorm
    # println("Constrained LS with 10_000 counts:")
    # println("   status    : $status")
    # println("   error     : $enorm")
    # println("   true state: $ρ")
    # println("   estimate  : $ρest")
    @test status == :Optimal || ( status == :UnknownError && enorm < 5e-2 )
    @test enorm < 5e-2

    status, enorm, _, ρest = test_qst_ml(ρ, obs, asymptotic=true)
    result[k,3] = enorm
    # println("Strict ML with mean counts:")
    # println("   status    : $status")
    # println("   error     : $enorm")
    # println("   true state: $ρ")
    # println("   estimate  : $ρest")
    @test status == :Optimal
    @test enorm < 5e-2

    status, enorm, _, ρest = test_qst_ml(ρ, obs, asymptotic=false)
    result[k,4] = enorm
    # println("Strict ML with 10_000 counts:")
    # println("   status    : $status")
    # println("   error     : $enorm")
    # println("   true state: $ρ")
    # println("   estimate  : $ρest")
    @test status == :Optimal || status == :MaxIter
    @test enorm < 5e-2

    status, enorm, _, ρest = test_qst_ml(ρ, obs, β=0.001, asymptotic=true)
    result[k,5] = enorm
    # println("Hedged ML with mean counts and β=0.04:")
    # println("   status    : $status")
    # println("   error     : $enorm")
    # println("   true state: $ρ")
    # println("   estimate  : $ρest")
    @test status == :Optimal
    @test enorm < 5e-2

    status, enorm, _, ρest = test_qst_ml(ρ, obs, β=0.001, asymptotic=false)
    result[k,6] = enorm
    # println("Hedged ML with 10_000 counts and β=0.04:")
    # println("   status    : $status")
    # println("   error     : $enorm")
    # println("   true state: $ρ")
    # println("   estimate  : $ρest")
    @test status == :Optimal
    @test enorm < 5e-2

    # status, enorm, _, ρest = test_qst_hml(ρ, obs, β=0.001, asymptotic=false)
    # result[k,6] = enorm
    # println("Hedged ML with 10_000 counts and β=0.04:")
    # println("   status    : $status")
    # println("   error     : $enorm")
    # println("   true state: $ρ")
    # println("   estimate  : $ρest")
    # println(enorm)
    #@test status == :Optimal
    #@test enorm < 5e-2

    #println(result[k,:])
end
@printf "\n"

#println(result)


text:
using BinaryBuilder, Pkg

name = "MrBayes"
version = v"3.2.7"
sources = [
    GitSource("https://github.com/NBISweden/MrBayes.git", "d50016695db24c58bcb36c83c487bd365fe2a566"),
]

script = raw"""
cd ${WORKSPACE}/srcdir/MrBayes
install_license ${WORKSPACE}/srcdir/MrBayes/COPYING

## MPI VERSION
if [[ "$target" == *-mingw* ]]; then
    COMPILER=${target}-gcc
    EXTRALIBS="-lm -lmsmpi -I${includedir} -L${libdir}"
else
    COMPILER=mpicc
    EXTRALIBS=""
fi

CC=$COMPILER \
LIBS="$EXTRALIBS" \
./configure \
    --with-mpi \
    --prefix=$prefix \
    --build=${MACHTYPE} \
    --host=${target} \
    --includedir="${includedir}" \
    --libdir="${libdir}"

make -j${nproc}
make install
install -Dvm 755 $bindir/mb${exeext} $bindir/mb_MPI${exeext}

## SEQUENTIAL VERSION
./configure --prefix=$prefix --build=${MACHTYPE} --host=${target} --includedir="${includedir}" --libdir="${libdir}"
make -j${nproc}
make install
"""

platforms = supported_platforms()

products = [
    ExecutableProduct("mb", :mb),
    ExecutableProduct("mb_MPI", :mb_MPI),
]

dependencies = [
    Dependency(PackageSpec(name="MPICH_jll", uuid="7cb0a576-ebde-5e09-9194-50597f1243b4"); platforms=filter(!Sys.iswindows, platforms)),
    Dependency(PackageSpec(name="MicrosoftMPI_jll", uuid="9237b28f-5490-5468-be7b-bb81f5f5e6cf"); platforms=filter(Sys.iswindows, platforms)),
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae"); platforms=filter(!Sys.isbsd, platforms)),
    Dependency(PackageSpec(name="LLVMOpenMP_jll", uuid="1d63c593-3942-5779-bab2-d838dc0a180e"); platforms=filter(Sys.isbsd, platforms)),
    Dependency(PackageSpec(name="Readline_jll", uuid="05236dd9-4125-5232-aa7c-9ec0c9b2c25a"); platforms=filter(Sys.isbsd, platforms)),
]

build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies, julia_compat="1.6")


text:
using Parameters

abstract type AbstractDiffusivity end
export AbstractDiffusivity

struct EmptyDiffusivity <: AbstractDiffusivity end
export EmptyDiffusivity

@with_kw struct ConstantDiffusivity{N<:Number} <: AbstractDiffusivity
    D::N = nothing
end
(c::ConstantDiffusivity)(;T::AbstractFloat=0.0,P::AbstractFloat=0.0,mu::AbstractFloat=0.0) = c.D
export ConstantDiffusivity

@with_kw struct StokesDiffusivity{N<:Number} <: AbstractDiffusivity
    r::N
end
(sd::StokesDiffusivity)(;T::N,mu::Q,P::R=0.0) where {N,R,Q<:Number} = @fastmath kB*T/(6*Base.pi*mu*sd.r)
export StokesDiffusivity


text:
module MultiLine

# using StaticArrays
using Interpolations
using Roots
using Plots
gr()

# Types
export Line, Point, Envelope

# methods
export interp, splitat,upper_env!, getx, gety, gets, splitLine, getr, bound

import Base.size,
       Base.getindex,
       Base.setindex!,
       Base.eltype,
       Base.prepend!,
       Base.append!,
       Base.insert!,
       Base.delete!,
       Base.sort!,
       Base.length,
       Base.show

# includes
include("line.jl")
include("envelope.jl")
include("plotting.jl")

end # module

max_stars_repo_name: floswald/MultiLine

max_issues_repo_name: floswald/MultiLine

max_forks_repo_name: floswald/MultiLine


text:
add_noise(X, d, N, noiseVar = 1e-10) = X .+ (noiseVar .* rand(d, N))

# Return the kNN nearest point.
function knn_nearest_point(X::Array{Float64, 2}, knnTree, point, k)
  indices, distances = nearest(knnTree, point, k+1)
  sorted = sortperm(distances)
  X[:,indices[sorted[end]]]
end

# Return the kNN nearest point.
function knn_nearest_distance(X::Array{Float64, 2}, knnTree, point, k)
  indices, distances = nearest(knnTree, point, k+1)
  sorted = sortperm(distances)
  distances[sorted[end]]
end


text:
using Fire

@main function main(x; opt1 = 1, opt2::Int = 2, flag = false)
    println("Parsed args:")
    println("flag=>", flag)
    println("arg=>", x)
    println("opt1=>", opt1)
    println("opt2=>", opt2)
end


text:
# This is a temporary module to validate `AbstractImageFilter` idea
# proposed in https://github.com/JuliaImages/ImagesAPI.jl/pull/3
module HistogramAdjustmentAPI

using ImageCore # ColorTypes is sufficient

# TODO Relax this to all image color types
const GenericGrayImage = AbstractArray{<:Union{Number, AbstractGray}}

"""
    AbstractImageAlgorithm

The root of image algorithms type system
"""
abstract type AbstractImageAlgorithm end

"""
    AbstractImageFilter <: AbstractImageAlgorithm

Filters are image algorithms whose input and output are both images
"""
abstract type AbstractImageFilter <: AbstractImageAlgorithm end

include("histogram_adjustment.jl")

# we do not export any symbols since we don't require
# package developers to implemente all the APIs

end  # module HistogramAdjustmentAPI


text:
using XMem
using Test

include("sample.jl")