text:
using Dates
using Test

using MarketData

using TimeSeries

@testset "split" begin

@testset "find methods" begin
    @testset "find returns correct row numbers array" begin
        @test timestamp(cl[findall(cl .> op)])[1] == Date(2000, 1, 3)
        @test length(findall(cl .> op))           == 244
    end

    @testset "findwhen returns correct Dates array" begin
       @test findwhen(cl .> op)[2]      == Date(2000, 1, 5)
       @test length(findwhen(cl .> op)) == 244
    end

    @testset "findall(f::Function, ta)" begin
        @test findall(cl .> 100) == findall(x -> x > 100, cl)
        @test findall(cl .> 100) == findall(x -> x[4] > 100, ohlc)
    end
end

@testset "split date operations" begin
    @testset "from and to correctly subset non-zero and zero-length time arrays" begin
        @test length(from(cl, Date(2001,12,28))) == 2
        @test length(from(cl, Date(2002,1,1)))   == 0
        @test length(from(from(cl, Date(2002,1,1)), Date(2012,1,1))) == 0

        @test length(to(cl, Date(2000,1,4)))     == 2
        @test length(to(cl, Date(1999,1,4)))     == 0
        @test length(to(to(cl, Date(1999,1,4)), Date(1912,1,1))) == 0
    end

    @testset "when method correctly subset" begin
        @test timestamp(when(cl, day, 4))[1]              == Date(2000,1,4)
        @test timestamp(when(cl, dayname, "Friday"))[1]   == Date(2000,1,7)
        @test timestamp(when(cl, week, 5))[1]             == Date(2000,1,31)
        @test timestamp(when(cl, month, 5))[1]            == Date(2000,5,1)
        @test timestamp(when(cl, monthname, "June"))[1]   == Date(2000,6,1)
        @test timestamp(when(cl, year, 2001))[1]          == Date(2001,1,2)
        @test timestamp(when(cl, dayofweek, 1))[1]        == Date(2000,1,3)
        # all the days in the nth week of each month
        @test timestamp(when(cl, dayofweekofmonth, 5))[1] == Date(2000,1,31)
        @test timestamp(when(cl, dayofyear, 365))[1]      == Date(2001,12,31)
        @test timestamp(when(cl, quarterofyear, 4))[1]    == Date(2000,10,2)
        @test timestamp(when(cl, dayofquarter, 1))[1]     == Date(2001,10,1)
    end
end

@testset "head, tail, first and last methods" begin
    @testset "head, tail, first and last methods work with default n value on single column TimeArray" begin
        @test length(head(cl,6))  == 6
        @test timestamp(head(cl)) == [Date(2000,1,3), Date(2000,1,4), Date(2000,1,5),
                                      Date(2000,1,6), Date(2000,1,7), Date(2000,1,10)]
        @test values(head(cl))    == [111.94, 102.5, 104.0, 95.0, 99.5, 97.75]

        @test length(tail(cl,6))  == 6
        @test timestamp(tail(cl)) == [Date(2001,12,21), Date(2001,12,24), Date(2001,12,26),
                                      Date(2001,12,27), Date(2001,12,28), Date(2001,12,31)]
        @test values(tail(cl))    == [21.0, 21.36, 21.49, 22.07, 22.43, 21.9]

        @test length(first(cl))       == 1
        @test timestamp(first(cl))[1] == Date(2000,1,3)
        @test values(first(cl))[1]    == 111.94
        @test meta(first(cl))         == "AAPL"

        @test length(last(cl))       == 1
        @test timestamp(last(cl))[1] == Date(2001,12,31)
        @test values(last(cl))[1]    == 21.9
        @test meta(last(cl))         == "AAPL"
    end

    @testset "head, tail, first and last methods work with default n value on multi column TimeArray" begin
        @test length(head(ohlc))         == 6
        @test values(head(ohlc, 1))      == [104.88 112.5 101.69 111.94]

        @test length(tail(ohlc))         == 6
        @test values(tail(ohlc, 1))      == [22.51 22.66 21.83 21.9]

        @test length(first(ohlc))       == 1
        @test timestamp(first(ohlc))[1] == Date(2000,1,3)
        @test values(first(ohlc))       == [104.88 112.5 101.69 111.94]
        @test meta(first(ohlc))         == "AAPL"

        @test length(last(ohlc))       == 1
        @test timestamp(last(ohlc))[1] == Date(2001,12,31)
        @test values(last(ohlc))       == [22.51 22.66 21.83 21.9]
        @test meta(last(ohlc))         == "AAPL"
    end

    @testset "head, tail, first and last methods work with custom periods on single column TimeArray" begin
        @test length(head(cl, 2))   == 2
        @test length(head(cl, 500)) == length(cl)
        @test length(tail(cl, 2))   == 2
        @test length(tail(cl, 500)) == length(cl)

        @test length(first(cl))     == 1
        @test length(last(cl))      == 1
    end

    @testset "head, tail, first and last methods work with custom periods on multi column TimeArray" begin
        @test length(head(ohlc, 2))   == 2
        @test length(head(ohlc, 500)) == length(ohlc)
        @test length(tail(ohlc, 2))   == 2
        @test length(tail(ohlc, 500)) == length(ohlc)

        @test length(first(ohlc))     == 1
        @test length(last(ohlc))      == 1
    end
end

end  # @testset "split"


text:
"""
    Household_StochasticDP

Stochastic DP for a standard household life-cycle model with one shock on consumption.
developed based on `Household.jl` and `src.jl`
"""
module Household_StochasticDP
    import src # math dependency
    import Household # for basic methods

# --------------------- CRRA UTILITY
function u_CRRA(c::Real ; γ::Real = 1.5, ψ::Real = 1E-6)
    # assertion
    @assert( 0 <= c < Inf, "c must be non-negative and finite" )
    @assert( 0 <= γ < Inf, "γ must be non-negative and finite" )
    @assert( ψ >= 0, "ψ must be non-negative" )
    # branch
    γ == 1 ? (return log(c + ψ)::Real) : nothing
    return ( ((c+ψ)^(1-γ) - 1) / (1-γ) )::Real
end # u_CRRA

end # Household_StochasticDP

ext: jl

max_stars_repo_name: Clpr/EconToolBox

max_issues_repo_name: Clpr/EconToolBox

max_forks_repo_name: Clpr/EconToolBox


text:
# by lazarusA # HIDE
# using GLMakie # HIDE
using CairoMakie
CairoMakie.activate!() # HIDE
let
    x = LinRange(0,2π,50) # with 50 points in the interval.
    fig = Figure(resolution = (700, 450), font =:sans)
    ax = Axis(fig, xlabel = "x", ylabel = "")
    line1 = lines!(x, sin.(x), color = :red)
    line2 = lines!(x, cos.(x), color = :blue)
    points1 = scatter!(x, cos.(x), color = :black, markersize = 5)
    points2 = scatter!(x, -cos.(x), color = :red, strokecolor = :red,
        markersize = 5, marker = '■')
    leg = Legend(fig, [line1, [line2, points1], points2],
        [" sin(x)", " cos(x)",  "-cos(x)"], "f(x)", markersize = 8,
        tellheight =false, tellwidth = true, valign = :center)
    fig[1, 1] = ax
    fig[1, 2] = leg
    fig
    #save("FigLinesAndScattersLegOut.png", fig, px_per_unit = 2)
    save(joinpath(@__DIR__, "output", "FigLinesAndScattersLegOut.png"), fig, px_per_unit = 2) # HIDE
end

using Pkg # HIDE
Pkg.status("CairoMakie") # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
#!/usr/bin/env julia

include("../src/mesh/trees.jl")

using .Trees: Tree, refine!, length, capacity, leaf_cells, refine_box!,
              minimum_level, maximum_level, count_leaf_cells
using TimerOutputs
using Profile
using HDF5: h5open, attrs

const ndim = 2

# Save current mesh with some context information as an HDF5 file.
function save_mesh_file(filename::String, tree::Tree)
  # Open file (clobber existing content)
  h5open(filename, "w") do file
    # Add context information as attributes
    n_cells = length(tree)
    attrs(file)["ndim"] = ndim
    attrs(file)["n_cells"] = n_cells
    attrs(file)["n_leaf_cells"] = count_leaf_cells(tree)
    attrs(file)["minimum_level"] = minimum_level(tree)
    attrs(file)["maximum_level"] = maximum_level(tree)
    attrs(file)["center_level_0"] = tree.center_level_0
    attrs(file)["length_level_0"] = tree.length_level_0

    # Add tree data
    file["parent_ids"] = @view tree.parent_ids[1:n_cells]
    file["child_ids"] = @view tree.child_ids[:, 1:n_cells]
    file["neighbor_ids"] = @view tree.neighbor_ids[:, 1:n_cells]
    file["levels"] = @view tree.levels[1:n_cells]
    file["coordinates"] = @view tree.coordinates[:, 1:n_cells]
  end

  return filename * ".h5"
end

to = TimerOutput()

coordinates_min = [-16, -16]
coordinates_max = [ 16,  16]
capacity_ = 1000
initial_refinement_level = 2

domain_center = (coordinates_min + coordinates_max) / 2
domain_length = maximum(coordinates_max - coordinates_min)

# Create tree object
@timeit to "create tree" t = Tree(Val{2}(), capacity_, domain_center, domain_length)
println("Initial tree:")
println(t)

# Create initial refinement
@timeit to "initial refinement" for l = 1:initial_refinement_level
  @timeit to "refine!" refine!(t)
  println("After uniform refinement to level $l:")
  println(t)
end

# Add patches
patches = [
           [[0.0, -16.0], [16.0, 0.0]],
           [[0.0, -16.0], [16.0,  0.0]],
           [[0.0, -16.0], [16.0, -8.0]],
           #=[[8.0, -16.0], [16.0, -8.0]],=#
          ]
for (coordinates_min, coordinates_max) in patches
  refine_box!(t, coordinates_min, coordinates_max)
  println("After refinement patch $coordinates_min, $coordinates_max:")
  println(t)
end

save_mesh_file("mesh_test.h5", t)
exit(0)

# Create non-uniform refinement
@timeit to "local refinement" begin
  @timeit to "refine!" refine!(t, 4)
  println("After refining node 4:")
  println(t)
  @timeit to "refine!" refine!(t, 5)
  println("After refining node 5:")
  println(t)
end

# Refine one level further
@timeit to "refine!" refine!(t)
println("After refining everything once more:")
println(t)

refine_box!(t, -8, 8)
println("After refining center nodes:")
println(t)

# Print tree information
@show leaf_nodes(t)
@show t.coordinates[:, leaf_nodes(t)]

print_timer(to)
println()


text:
using PQEcolaPoint
using Dates

# what is the contract?
short_put_contract = PutContractModel()
short_put_contract.ticker = "XYZ"
short_put_contract.expiration_date = Date(2022, 04, 14)
short_put_contract.strike_price = 65.0
short_put_contract.premium = 0.0
short_put_contract.number_of_contracts = 1
short_put_contract.direction = 1

# other version -
# origin_point = encode(JIT_lattice; S = 70.12, DTE = 53.0, IV = 0.4926, K = 65.0)
val_γ = γ(short_put_contract; Sₒ=70.12, number_of_levels=80, σ=0.5129, T=(53.0 / 365), μ=0.0047)


text:
const Simplex{M, N, T} = Union{MVector{M, SVector{N, T}}, SVector{M, SVector{N, T}}}
const SimplexType{M, N, T} = Union{Type{MVector{M, SVector{N, T}}}, Type{SVector{M, SVector{N, T}}}}

dimension(::SimplexType{M, N, T}) where {M,N,T} = Val(N)

any_inside(simplex::Simplex{M, N, T}) where {M,N,T} = Tagged(simplex[1])

function support_vector_max(simplex::Simplex{M, N, T}, direction, initial_guess::Tagged) where {M,N,T}
    best_pt, score = gt.argmax(p -> dot(p, direction), simplex)
    Tagged(best_pt)
end

@generated function simplex_face(simplex::Simplex{M, N, T}, i::Integer) where {N,M,T}
    simplex_face_impl(simplex, i)
end

function simplex_face_impl(simplex::SimplexType{M, N, T}, i) where {N,M,T}
    Expr(:call, :(SVector),
        Expr(:tuple, [:(i > $j ? simplex[$j] : simplex[$(j+1)]) for j in 1:(M - 1)]...))
end

"""
    normal(face::StaticVector{N, <:StaticVector{N}})

Return any `N`-vector normal to the face spanned by `N` points. The result is not necessarily normalized,
nor is there any guarantee on winding.
"""
function normal end

function normal(face::StaticVector{2, <:StaticVector{2}})
    Δ = face[2] - face[1]
    similar_type(Δ)(-Δ[2], Δ[1])
end

function normal(face::StaticVector{3, <:StaticVector{3}})
    Δ1 = face[2] - face[1]
    Δ2 = face[3] - face[1]
    Δ1 × Δ2
end

function normal(face::StaticVector{N, <:StaticVector{N}}) where N
    Δ = hcat(ntuple(i -> face[i + 1] - face[1], Val(N - 1))...)
    factorization = svd(Δ, full=Val(true))
    factorization.U[:, end]
end


text:
"""
	HardThresholding!(data::AbstractArray, λ::AbstractFloat)

Hard Thresholding function
λ is threshold parameter
Return the number of non-zero elements in the array
"""
function HardThresholding!(data::AbstractArray{T, N}, λ::AbstractFloat) where {T <:Number, N}
	i_end = length(data)
	nnz = i_end
	@views @inbounds @fastmath for i in 1:i_end
		if abs(data[i]) <= λ
			data[i] = zero(T)
			nnz -= 1
		end
	end
	return nnz
end


text:
# p = Params(Ln = 400, W = 400, Ls = 40, Ws = 20, scale = 40, λ = 5., U = 0, α = 0,
#     Δ = 1, d = 0, τ = 1, B = SA[0,0,0], EZ = SA[0,0,1e-6], μN = 1e-6);

# p = Params(Ln = 50, W = 80, Ls = 10, Ws = 0, scale = 20, λ = 5., U = 0, α = 0,
# Δ = 5, d = 0, τ = 1, B = SA[0,0,0], EZ = SA[0,0,1e-6], μN = 1e-6);

p = Params(Ln = 200, W = 320, Ls = 30, Ws = 0, scale = 20, λ = 5., U = 0, α = 0,
Δ = 5, d = 0, τ = 1, B = SA[0,0,0], EZ = SA[0,0,1e-6], μN = 1e-6);

fΦ = 2π *BLG_SQUID.quantumflux(p)
Blist = collect(0.:fΦ/4:5*fΦ)

# ph = rectangle_weaklink(p, false, true, Δx_mask = 20, Δy_mask = 20)
# ic = maxicϕ_exactdiag(Blist, p, mask=true, Δx_mask = 30, Δy_mask = 30)
ic = icϕ_exactdiag(Blist, p, mask=true, Δx_mask = 20, Δy_mask = 20)

@time sp = spectrum(ph(ϕ = 0), method = ArpackPackage(nev = 204, sigma = 0.0im))

function ldosonlattice(number)
    println(sp.energies[number])
   return ldosonlattice(sp.states[:,number], ph(ϕ = 0))
end

ph = rectangle_weaklink(p, false, false)

@time sp = spectrum(ph(ϕ = 0), method = ArpackPackage(nev = 192, sigma = 0.0im))

##########
# check regime

function Bcriteria(B::Array, scale; threshold = 50.)
    ħc = ustrip(u"meV*nm",ħ*c_0)
    a0 = 0.246    #nm
    Bmax = findmax(abs.(B))[1]
    l_B = ħc/Bmax
    l_scale = scale * a0
    println(l_B/ l_scale)
    println(l_B/a0)
    l_B/ l_scale > threshold ? nothing : @warn("lB ≈ a0*scale, decrease scaling")
    l_B/ a0 > threshold ? nothing : @warn("lB ≈ a0, increase BLG area (we are in t
        he Hofstadter regime)")
end

#######################
#
p = Params(Ln = 10, W = 10, Ls = 10, Ws = 0, scale = 10, λ = 5., U = 0, α = 0,
Δ = 5, d = 0, τ = 1, B = SA[0,0,0], EZ = SA[0,0,1e-6], μN = 1e-6);
Blist = [0.]
ic = icϕ_exactdiag(Blist, p, mask = true, Δx_mask = 30, Δy_mask = 30)


text:
module CLagrangianFESpacesTests

using Test
using Gridap.ReferenceFEs
using Gridap.Geometry
using Gridap.TensorValues
using Gridap.Arrays
using Gridap.FESpaces
using Gridap.CellData

domain = (0,1,0,1)
partition = (2,2)
model = CartesianDiscreteModel(domain,partition)
grid = Triangulation(model)
x = collect1d(get_node_coordinates(grid))

V = CLagrangianFESpace(Float64,grid)
cellmat = [rand(4,4) for cell in 1:num_cells(model)]
cellvec = [rand(4) for cell in 1:num_cells(model)]
cellmatvec = pair_arrays(cellmat,cellvec)
test_single_field_fe_space(V,cellmatvec,cellmat,cellvec,grid)

@test get_cell_dof_ids(V) == collect1d(get_cell_node_ids(grid))
@test V.metadata.free_dof_to_node == collect(1:num_nodes(grid))
@test V.metadata.free_dof_to_comp == ones(Int,num_nodes(grid))
@test V.metadata.node_and_comp_to_dof == V.metadata.free_dof_to_node

u(x) = x[1]+x[2]
uh = interpolate(u,V)
uhx = get_free_dof_values(uh)
ux = u.(x)
@test uhx ≈ ux

V = CLagrangianFESpace(VectorValue{2,Float64},grid)
cellmat = [rand(8,8) for cell in 1:num_cells(model)]
cellvec = [rand(8) for cell in 1:num_cells(model)]
cellmatvec = pair_arrays(cellmat,cellvec)
test_single_field_fe_space(V,cellmatvec,cellmat,cellvec,grid)

@test get_cell_dof_ids(V) == [[1,3,7,9,2,4,8,10], [3,5,9,11,4,6,10,12], [7,9,13,15,8,10,14,16], [9,11,15,17,10,12,16,18]]
@test V.metadata.free_dof_to_node == [1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8, 9, 9]
@test V.metadata.free_dof_to_comp == [1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2]
@test V.metadata.node_and_comp_to_dof == VectorValue{2,Int}[(1,2),(3,4),(5,6),(7,8),(9,10),(11,12),(13,14),(15,16),(17,18)]

u(x) = VectorValue(x[1]+x[2],x[2])
uh = interpolate(u,V)
uhx = get_free_dof_values(uh)
ux = u.(x)
@test uhx ≈ collect1d(reinterpret(ux))

model = CartesianDiscreteModel(domain,partition)
grid = Triangulation(model)
labels = get_face_labeling(model)
tags = [1,2,4,5,8]
node_to_tag = get_face_tag_index(labels,tags,0)
masks1 = Bool[true,true,false,true,true]
V = CLagrangianFESpace(Float64,grid,Vector{Float64},node_to_tag,masks1)
@test get_cell_dof_ids(V) == [[-1, -2, 1, 2], [-2, -3, 2, -4], [1, 2, 3, 4], [2, -4, 4, 5]]
@test V.metadata.free_dof_to_node == [4, 5, 7, 8, 9]
@test V.metadata.free_dof_to_comp == [1, 1, 1, 1, 1]
@test V.metadata.dirichlet_dof_to_node == [1,2,3,6]
@test V.metadata.dirichlet_dof_to_comp == [1,1,1,1]
@test V.metadata.node_and_comp_to_dof == [-1, -2, -3, 1, 2, -4, 3, 4, 5]

masks2 = [(true,true),(true,false),(false,false),(false,true),(true,true)]
V = CLagrangianFESpace(VectorValue{2,Float64},grid,Vector{Float64},node_to_tag,masks2)
@test get_cell_dof_ids(V) == [[-1, 1, 3, 5, -2, -3, 4, 6], [1, -4, 5, -5, -3, 2, 6, -6], [3, 5, 7, 9, 4, 6, 8, 10], [5, -5, 9, 11, 6, -6, 10, 12]]
@test V.metadata.free_dof_to_node == [2, 3, 4, 4, 5, 5, 7, 7, 8, 8, 9, 9]
@test V.metadata.free_dof_to_comp == [1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2]
@test V.metadata.dirichlet_dof_to_node == [1, 1, 2, 3, 6, 6]
@test V.metadata.dirichlet_dof_to_comp == [1, 2, 2, 1, 1, 2]
@test V.metadata.node_and_comp_to_dof == VectorValue{2, Int32}[(-1, -2), (1, -3), (-4, 2), (3, 4), (5, 6), (-5, -6), (7, 8), (9, 10), (11, 12)]

reffe = ReferenceFE(lagrangian,Float64,2)
V = FESpace(model,reffe)
@test V.metadata === nothing

reffe = ReferenceFE(lagrangian,Float64,1)
V = FESpace(model,reffe,conformity=:L2)
@test V.metadata === nothing

# Check that the factory uses clagrangian when possible

V = FESpace(model,reffe)
@test V.metadata.node_and_comp_to_dof == [1, 2, 3, 4, 5, 6, 7, 8, 9]

V = FESpace(model,reffe,dirichlet_tags=tags)
@test V.metadata.node_and_comp_to_dof == [-1, -2, -3, 1, 2, -4, 3, 4, -5]

V = FESpace(model,reffe,dirichlet_tags=tags,dirichlet_masks=masks1)
@test V.metadata.node_and_comp_to_dof == [-1, -2, -3, 1, 2, -4, 3, 4, 5]

reffe = ReferenceFE(lagrangian,VectorValue{2,Float64},1)
V = FESpace(model,reffe)
@test V.metadata.node_and_comp_to_dof == VectorValue{2, Int32}[(1, 2), (3, 4), (5, 6), (7, 8), (9, 10), (11, 12), (13, 14), (15, 16), (17, 18)]

V = FESpace(model,reffe,dirichlet_tags=tags)
@test V.metadata.node_and_comp_to_dof  == VectorValue{2, Int32}[(-1, -2), (-3, -4), (-5, -6), (1, 2), (3, 4), (-7, -8), (5, 6), (7, 8), (-9, -10)]

V = FESpace(model,reffe,dirichlet_tags=tags,dirichlet_masks=masks2)
@test V.metadata.node_and_comp_to_dof == VectorValue{2, Int32}[(-1, -2), (1, -3), (-4, 2), (3, 4), (5, 6), (-5, -6), (7, 8), (9, 10), (11, 12)]

# Do not use CLagrangian for models with periodic boundary conditions

domain = (0,1,0,1)
partition = (3,3)
model = CartesianDiscreteModel(domain,partition,isperiodic=(false,true))
reffe = ReferenceFE(lagrangian,Float64,1)
V = FESpace(model,reffe)
@test V.metadata === nothing

end # module


text:
"""
    N(g, x) -> Int

Given a sequence of guesses `g`, returns the number of guesses needed to
correctly guess `x`. Returns `length(g)+1` if the correct answer is never
guesses.
"""
function N(g, x)
    i = findfirst(==(x), g)
    if i === nothing
        return length(g) + 1
    else
        return i
    end
end

"""
    cartesian_power(J, K) -> iterator

Returns an interator over the set of tuples of length `K` with elements in `1:J`.

# Example

```jldoctest
julia> collect(GuessworkQuantumSideInfo.cartesian_power(3,2))
9-element Array{Tuple{Int64,Int64},1}:
 (1, 1)
 (2, 1)
 (3, 1)
 (1, 2)
 (2, 2)
 (3, 2)
 (1, 3)
 (2, 3)
 (3, 3)
```
"""
function cartesian_power(J, K)
    sz = ntuple(i -> 1:J, K)
    inds = CartesianIndices(sz)
    return (Tuple(inds[y]) for y = 1:J^K)
end

"""
    make_povm_outcomes(J, K, remove_repetition)

Returns an iterator over tuples or vectors each representing a POVM outcome. If
`remove_repetition` is true, it returns ``J!/(J-K)!`` outcomes, vectors of
length `K` with elements in ``[1,...,J]`` such that each vector has no
repetition of elements. If `remove_repetition` is false, it returns all ``J^K``
elements of ``[1, ..., J]^K``.
"""
function make_povm_outcomes(J, K, remove_repetition::Bool)
    if remove_repetition
        return multiset_permutations(1:J, K)
    else
        return cartesian_power(J, K)
    end
end

"""
    compatible_problem(f, objective, constraints, numeric_type) -> Convex.Problem

Returns a `Convex.Problem` object corresponding to `f`, which should be either
`maximize` or `minimize`, the objective and constraints given as arguments, and
the `numeric_type`. Since only Float64 numeric types are supported on Convex
versions 0.12 and below, the argument `numeric_type` is only passed to `f` when
it is different from `Float64`.
"""
function compatible_problem(f, objective, constraints, numeric_type)
    if numeric_type != Float64
        f(objective, constraints; numeric_type = numeric_type)
    else
        f(objective, constraints)
    end
end

"""
    guesswork(
        p::AbstractVector{T},
        ρBs::AbstractVector{<:AbstractMatrix};
        solver,
        K::Integer = length(p),
        c = T[1:K..., 5_000],
        dual::Bool = false,
        remove_repetition::Bool = true,
        povm_outcomes = make_povm_outcomes(length(p), K, remove_repetition),
        verbose::Bool = true,
    )

Computes the guesswork for the c-q state specified by a probability vector `p`,
giving the distribution `X`, and `ρBs`, giving the associated quantum states.

The keyword arguments are as follows:

* `solver` is the only required keyword argument; an SDP solver such as SCS or
  MOSEK must be passed.
* `K` corresponds to the maximum number of allowed guesses. The number of
  variables in the primal SDP (and the number of constraints in the dual SDP)
  scales as `length(p)^K`.
* `c` may be given a custom cost vector. If `K < length(p)`, then `c` should be
  of length `K+1`. The last entry, `c[K+1]`, corresponds to the cost of not
  guessing the correct answer within `K` guesses.
* `dual` is a boolean variable indicating whether the primal or dual
  optimization problem should be solved.
* `remove_repetition` is a boolean variable defaulting to true, indicating
  whether repeated guesses of the same value should be removed; as long as `c`
  is increasing, this decreases the size of the SDP without affecting the
  optimal value.
* `povm_outcomes` should be an iterator (or vector) corresponding to the
  possible guessing orders. This defaults to all subsets of length `K` of
  `1:length(p)` without repetition.
* `verbose` is a boolean which indicates if warnings should be printed when the
  problem is not solved optimally.
"""
function guesswork(
    p::AbstractVector{T},
    ρBs::AbstractVector{<:AbstractMatrix};
    solver,
    K::Integer = length(p),
    c = T[1:K..., 5_000],
    dual::Bool = false,
    remove_repetition::Bool = true,
    povm_outcomes = make_povm_outcomes(length(p), K, remove_repetition),
    verbose::Bool = true,
    debug::Bool = false,
) where {T<:Number}
    length(p) == length(ρBs) ||
    throw(ArgumentError("Length of prior and vector of side information must match J"))

    J = length(p)
    dB = size(ρBs[1], 1)
    all(ρB -> size(ρB) == (dB, dB), ρBs) ||
    throw(ArgumentError("All side-information states must be square matrices of the same size"))

    R(g) = sum(c[N(g, x)] * p[x] * ρBs[x] for x in eachindex(p, ρBs))

    num_outcomes = length(povm_outcomes)

    if dual
        Y = ComplexVariable(dB, dB)

        constraints = Convex.Constraint[R(outcome) ⪰ Y for outcome in povm_outcomes]
        push!(constraints, Y' == Y)

        objective = real(tr(Y))

        problem = compatible_problem(Convex.maximize, objective, constraints, T)
        solve!(problem, solver; verbose = verbose)

        # introduced for debugging a problem with SCS with MathOptInterface
        if debug
            @show length(constraints)
            for outcome in povm_outcomes
                @show outcome
                @show eigmin(Hermitian(R(outcome) - evaluate(Y)))
            end
        end

        Y = evaluate(Y)

        # SDP dual values not implemented in Convex.jl yet
        # Es = [Matrix{ComplexF64}(problem.constraints[y].dual) for y = eachindex(povm_outcomes)]
        Es = missing

    else

        Es = [ComplexVariable(dB, dB) for _ = 1:num_outcomes]

        objective =
            real(sum([tr(Es[y] * R(outcome)) for (y, outcome) in enumerate(povm_outcomes)]))
        constraints = Convex.Constraint[E ⪰ 0 for E in Es]
        push!(constraints, sum(Es) == I(dB))

        problem = compatible_problem(Convex.minimize, objective, constraints, T)
        solve!(problem, solver; verbose = verbose)

        Es = evaluate.(Es)
        # SDP dual values not implemented in Convex.jl yet
        # Y = Matrix{ComplexF64}(problem.constraints[end].dual)
        Y = missing
    end

    input_data = (J = J, K = K, c = c, p = p, ρBs = ρBs, dB = dB)

    return (
        optval = problem.optval,
        status = problem.status,
        Y = Y,
        Es = Es,
        povm_outcomes = povm_outcomes,
        input_data...,
    )
end


text:
using OVERTVerify
using Documenter

makedocs(;
    modules=[OvertVerify],
    authors="Amir Maleki, Chelsea Sidrane",
    repo="https://github.com/sisl/OVERTVerify.jl/blob/{commit}{path}#L{line}",
    sitename="OVERTVerify.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://sisl.github.io/OVERTVerify.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/sisl/OVERTVerify.jl",
)


text:
using Documenter, DispatchedTuples

format = Documenter.HTML(
    prettyurls = !isempty(get(ENV, "CI", "")),
    collapselevel = 1,
)

makedocs(
    sitename = "DispatchedTuples.jl",
    strict = true,
    format = format,
    checkdocs = :exports,
    clean = true,
    doctest = true,
    modules = [DispatchedTuples],
    pages = Any[
        "Home" => "index.md",
        "Performance" => "performance.md",
        "API" => "api.md",
    ],
)

deploydocs(
    repo = "github.com/charleskawczynski/DispatchedTuples.jl.git",
    target = "build",
    push_preview = true,
    devbranch = "main",
    forcepush = true,
)


text:
export cf

import ..CFModel

using mlpack._Internal.io

import mlpack_jll
const cfLibrary = mlpack_jll.libmlpack_julia_cf

# Call the C binding of the mlpack cf binding.
function cf_mlpackMain()
  success = ccall((:cf, cfLibrary), Bool, ())
  if !success
    # Throw an exception---false means there was a C++ exception.
    throw(ErrorException("mlpack binding error; see output"))
  end
end

" Internal module to hold utility functions. "
module cf_internal
  import ..cfLibrary

import ...CFModel

# Get the value of a model pointer parameter of type CFModel.
function IOGetParamCFModel(paramName::String)::CFModel
  CFModel(ccall((:IO_GetParamCFModelPtr, cfLibrary), Ptr{Nothing}, (Cstring,), paramName))
end

# Set the value of a model pointer parameter of type CFModel.
function IOSetParamCFModel(paramName::String, model::CFModel)
  ccall((:IO_SetParamCFModelPtr, cfLibrary), Nothing, (Cstring, Ptr{Nothing}), paramName, model.ptr)
end

# Serialize a model to the given stream.
function serializeCFModel(stream::IO, model::CFModel)
  buf_len = UInt[0]
  buf_ptr = ccall((:SerializeCFModelPtr, cfLibrary), Ptr{UInt8}, (Ptr{Nothing}, Ptr{UInt}), model.ptr, Base.pointer(buf_len))
  buf = Base.unsafe_wrap(Vector{UInt8}, buf_ptr, buf_len[1]; own=true)
  write(stream, buf)
end
# Deserialize a model from the given stream.
function deserializeCFModel(stream::IO)::CFModel
  buffer = read(stream)
  CFModel(ccall((:DeserializeCFModelPtr, cfLibrary), Ptr{Nothing}, (Ptr{UInt8}, UInt), Base.pointer(buffer), length(buffer)))
end
end # module

"""
    cf(; [algorithm, all_user_recommendations, input_model, interpolation, iteration_only_termination, max_iterations, min_residue, neighbor_search, neighborhood, normalization, query, rank, recommendations, seed, test, training, verbose])

This program performs collaborative filtering (CF) on the given dataset. Given a
list of user, item and preferences (the `training` parameter), the program will
perform a matrix decomposition and then can perform a series of actions related
to collaborative filtering.  Alternately, the program can load an existing saved
CF model with the `input_model` parameter and then use that model to provide
recommendations or predict values.

The input matrix should be a 3-dimensional matrix of ratings, where the first
dimension is the user, the second dimension is the item, and the third dimension
is that user's rating of that item.  Both the users and items should be numeric
indices, not names. The indices are assumed to start from 0.

A set of query users for which recommendations can be generated may be specified
with the `query` parameter; alternately, recommendations may be generated for
every user in the dataset by specifying the `all_user_recommendations`
parameter.  In addition, the number of recommendations per user to generate can
be specified with the `recommendations` parameter, and the number of similar
users (the size of the neighborhood) to be considered when generating
recommendations can be specified with the `neighborhood` parameter.

For performing the matrix decomposition, the following optimization algorithms
can be specified via the `algorithm` parameter:
 - 'RegSVD' -- Regularized SVD using a SGD optimizer
 - 'NMF' -- Non-negative matrix factorization with alternating least squares
update rules
 - 'BatchSVD' -- SVD batch learning
 - 'SVDIncompleteIncremental' -- SVD incomplete incremental learning
 - 'SVDCompleteIncremental' -- SVD complete incremental learning
 - 'BiasSVD' -- Bias SVD using a SGD optimizer
 - 'SVDPP' -- SVD++ using a SGD optimizer

The following neighbor search algorithms can be specified via the
`neighbor_search` parameter:
 - 'cosine'  -- Cosine Search Algorithm
 - 'euclidean'  -- Euclidean Search Algorithm
 - 'pearson'  -- Pearson Search Algorithm

The following weight interpolation algorithms can be specified via the
`interpolation` parameter:
 - 'average'  -- Average Interpolation Algorithm
 - 'regression'  -- Regression Interpolation Algorithm
 - 'similarity'  -- Similarity Interpolation Algorithm

The following ranking normalization algorithms can be specified via the
`normalization` parameter:
 - 'none'  -- No Normalization
 - 'item_mean'  -- Item Mean Normalization
 - 'overall_mean'  -- Overall Mean Normalization
 - 'user_mean'  -- User Mean Normalization
 - 'z_score'  -- Z-Score Normalization

A trained model may be saved to with the `output_model` output parameter.

To train a CF model on a dataset `training_set` using NMF for decomposition and
saving the trained model to `model`, one could call:

```julia
julia> using CSV
julia> training_set = CSV.read("training_set.csv")
julia> _, model = cf(algorithm="NMF", training=training_set)
```

Then, to use this model to generate recommendations for the list of users in the
query set `users`, storing 5 recommendations in `recommendations`, one could
call

```julia
julia> using CSV
julia> users = CSV.read("users.csv"; type=Int)
julia> recommendations, _ = cf(input_model=model, query=users,
            recommendations=5)
```

# Arguments

 - `algorithm::String`: Algorithm used for matrix factorization.  Default
      value `NMF`.

 - `all_user_recommendations::Bool`: Generate recommendations for all
      users.  Default value `false`.

 - `input_model::CFModel`: Trained CF model to load.
 - `interpolation::String`: Algorithm used for weight interpolation.
      Default value `average`.

 - `iteration_only_termination::Bool`: Terminate only when the maximum
      number of iterations is reached.  Default value `false`.

 - `max_iterations::Int`: Maximum number of iterations. If set to zero,
      there is no limit on the number of iterations.  Default value `1000`.

 - `min_residue::Float64`: Residue required to terminate the factorization
      (lower values generally mean better fits).  Default value `1e-05`.

 - `neighbor_search::String`: Algorithm used for neighbor search.  Default
      value `euclidean`.

 - `neighborhood::Int`: Size of the neighborhood of similar users to
      consider for each query user.  Default value `5`.

 - `normalization::String`: Normalization performed on the ratings.
      Default value `none`.

 - `query::Array{Int, 2}`: List of query users for which recommendations
      should be generated.
 - `rank::Int`: Rank of decomposed matrices (if 0, a heuristic is used to
      estimate the rank).  Default value `0`.

 - `recommendations::Int`: Number of recommendations to generate for each
      query user.  Default value `5`.

 - `seed::Int`: Set the random seed (0 uses std::time(NULL)).  Default
      value `0`.

 - `test::Array{Float64, 2}`: Test set to calculate RMSE on.
 - `training::Array{Float64, 2}`: Input dataset to perform CF on.
 - `verbose::Bool`: Display informational messages and the full list of
      parameters and timers at the end of execution.  Default value `false`.

# Return values

 - `output::Array{Int, 2}`: Matrix that will store output
      recommendations.
 - `output_model::CFModel`: Output for trained CF model.

"""
function cf(;
            algorithm::Union{String, Missing} = missing,
            all_user_recommendations::Union{Bool, Missing} = missing,
            input_model::Union{CFModel, Missing} = missing,
            interpolation::Union{String, Missing} = missing,
            iteration_only_termination::Union{Bool, Missing} = missing,
            max_iterations::Union{Int, Missing} = missing,
            min_residue::Union{Float64, Missing} = missing,
            neighbor_search::Union{String, Missing} = missing,
            neighborhood::Union{Int, Missing} = missing,
            normalization::Union{String, Missing} = missing,
            query = missing,
            rank::Union{Int, Missing} = missing,
            recommendations::Union{Int, Missing} = missing,
            seed::Union{Int, Missing} = missing,
            test = missing,
            training = missing,
            verbose::Union{Bool, Missing} = missing,
            points_are_rows::Bool = true)
  # Force the symbols to load.
  ccall((:loadSymbols, cfLibrary), Nothing, ());

  IORestoreSettings("Collaborative Filtering")

  # Process each input argument before calling mlpackMain().
  if !ismissing(algorithm)
    IOSetParam("algorithm", convert(String, algorithm))
  end
  if !ismissing(all_user_recommendations)
    IOSetParam("all_user_recommendations", convert(Bool, all_user_recommendations))
  end
  if !ismissing(input_model)
    cf_internal.IOSetParamCFModel("input_model", convert(CFModel, input_model))
  end
  if !ismissing(interpolation)
    IOSetParam("interpolation", convert(String, interpolation))
  end
  if !ismissing(iteration_only_termination)
    IOSetParam("iteration_only_termination", convert(Bool, iteration_only_termination))
  end
  if !ismissing(max_iterations)
    IOSetParam("max_iterations", convert(Int, max_iterations))
  end
  if !ismissing(min_residue)
    IOSetParam("min_residue", convert(Float64, min_residue))
  end
  if !ismissing(neighbor_search)
    IOSetParam("neighbor_search", convert(String, neighbor_search))
  end
  if !ismissing(neighborhood)
    IOSetParam("neighborhood", convert(Int, neighborhood))
  end
  if !ismissing(normalization)
    IOSetParam("normalization", convert(String, normalization))
  end
  if !ismissing(query)
    IOSetParamUMat("query", query, points_are_rows)
  end
  if !ismissing(rank)
    IOSetParam("rank", convert(Int, rank))
  end
  if !ismissing(recommendations)
    IOSetParam("recommendations", convert(Int, recommendations))
  end
  if !ismissing(seed)
    IOSetParam("seed", convert(Int, seed))
  end
  if !ismissing(test)
    IOSetParamMat("test", test, points_are_rows)
  end
  if !ismissing(training)
    IOSetParamMat("training", training, points_are_rows)
  end
  if verbose !== nothing && verbose === true
    IOEnableVerbose()
  else
    IODisableVerbose()
  end

  IOSetPassed("output")
  IOSetPassed("output_model")
  # Call the program.
  cf_mlpackMain()

  return IOGetParamUMat("output", points_are_rows),
         cf_internal.IOGetParamCFModel("output_model")
end


text:
__precompile__(true)

"""
Public holiday database
"""
module Holidays

export Holiday, upcoming, recent, nearest

using Base: Dates, Order

using DataStructures

abstract type Source end

"""Holiday type"""
type Holiday
    name::String
    description::String
    regions::Set{String}
    date::Date
end
Holiday(
    name::AbstractString, description::AbstractString,
    regions::AbstractString, date::AbstractString
) = Holiday(name, description, Set{String}(split(regions, '|')), Date(date))

"""Store for holidays"""
type HolidayStore
    holidays::Dict{String, SortedDict{Date, Set{Holiday}}}
    HolidayStore() = new(Dict{String, SortedDict{Date, Set{Holiday}}}())
end

Base.getindex(store::HolidayStore, region::String) = haskey(store.holidays, region) ? store.holidays[region] : error("Unknown region: $region")
Base.getindex(store::HolidayStore, x::Tuple{String, Date}) = get(store.holidays[first(x)], last(x), Set{Holiday}())
Base.getindex(store::HolidayStore, x::Tuple{String, DateTime}) = getindex(store, (first(x), Date(last(x))))
Base.getindex(store::HolidayStore, x::Date) = Set{Holiday}(v for region in keys(store.holidays) for v in store[(region, x)])

"""Gets next upcoming holiday date for region."""
function upcoming(store::HolidayStore, region::String, date::Date; value::Bool=false)
    holidays = store[region]
    date = try deref_key((holidays, searchsortedfirst(holidays, date))) catch lastdayofyear(today()) + Day(1) end
    value ? holidays[date] : date
end

"""Gets most recent holiday date for region."""
function recent(store::HolidayStore, region::String, date::Date; value::Bool=false)
    holidays = store[region]
    date = try deref_key((holidays, searchsortedlast(holidays, date))) catch firstdayofyear(today()) end
    value ? holidays[date] : date
end

"""Gets nearest holiday date for region. Returns upcoming holiday date if tied."""
function nearest(store::HolidayStore, region::String, date::Date; value::Bool=false)
    x, y = recent(store, region, date), upcoming(store, region, date)
    date = (y - date) <= (date - x) ? y : x
    value ? store[region][date] : date
end

"""Loads holidays from a source into the store."""
function load!(store::HolidayStore, holidays::Vector{Holiday})
    for holiday in holidays
        for region in holiday.regions
            if !haskey(store.holidays, region)
                target = store.holidays[region] = valtype(store.holidays)()
                target[holiday.date] = valtype(target)()
            else
                target = store.holidays[region]
                if !haskey(target, holiday.date)
                    target[holiday.date] = valtype(target)()
                end
            end
            push!(target[holiday.date], holiday)
        end
    end
    store
end

const _store = HolidayStore()

"""Load the prepared holiday database"""
function __init__()
    DB_DIR = joinpath(@__DIR__, "..", "deps", "db")
    for source in readdir(DB_DIR)
        files = [joinpath(DB_DIR, source, x) for x in readdir(joinpath(DB_DIR, source))]
        holidays = [Holiday(split(strip(line), '\t')...) for file in files for line in readlines(file)]
        # pick newest version for region/date/name
        holidays = unique((x) -> (x.regions, x.date, replace(lowercase(x.name), " ", "")), holidays)
        load!(_store, holidays)
    end
end

# bind to internal store
upcoming(region::String, date::Date; value::Bool=false) = upcoming(_store, region, date; value=value)
recent(region::String, date::Date; value::Bool=false) = recent(_store, region, date; value=value)
nearest(region::String, date::Date; value::Bool=false) = nearest(_store, region, date; value=value)

end


text:
using Gtk

function twoentrywindow()
    txt = "Enter Text Here"
    txtchanged = false

    win = GtkWindow("Keypress Test", 500, 100) |> (GtkFrame() |> (vbox = GtkBox(:v)))
    lab = GtkLabel("Enter some text in the first box and 7500 into the second box.")
    txtent = GtkEntry()
    set_gtk_property!(txtent,:text,"Enter Some Text Here")
    nument = GtkEntry()
    set_gtk_property!(nument,:text,"Enter the number seventy-five thousand here")
    push!(vbox, lab, txtent, nument)

    function keycall(w, event)
        strtxt = get_gtk_property(txtent, :text, String)
        numtxt = get_gtk_property(nument, :text, String)
        if strtxt != txt && occursin("75000", numtxt)
            set_gtk_property!(lab, :label, "You have accomplished the task.")
        end
    end
    signal_connect(keycall, win, "key-press-event")

    cond = Condition()
    endit(w) = notify(cond)
    signal_connect(endit, win, :destroy)
    showall(win)
    wait(cond)
end

twoentrywindow()

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Eive

n = 500
x = rand(n)
otherx = randn(n)
betas = [5.0, 5]
e = randn(n)
u = randn(n)

dirtyx = x .+ u

y = 5 .+ 5 .* x .+ 5 .* otherx .+ e

X = hcat(ones(n), x, otherx)
println("Clean Regression:")
println(X \ y)

Xd = hcat(ones(n), dirtyx, otherx)
println("Dirty Regression:")
println(Xd \ y)

println("EIVE Regression")
result = eive(dirtyx, y, otherx)
println(result)


text:
"Check whether geosteiner executables are in PATH"
function has_geosteiner()
    for cmd in [`efst -h`, `bb -h`]
        try
            run(pipeline(cmd, stdout=devnull, stderr=devnull))
        catch exc
            isa(exc, Base.IOError) && return false
        end
    end
    true
end

"""
Compute a Euclidean Steiner Minimal Tree for given nodes.

Is implemented by "shelling out" to GeoSteiner.
"""
function euclidean_steiner_tree(nodes::Vector{Node})
    @assert has_geosteiner() "Need GeoSteiner executables (efts, bb) in PATH!"

    tnodes = copy(nodes)
    tarcs = Arc[]

    process = open(pipeline(`efst`, `bb`), "r+")
    for node in nodes
        write(process.in, "$(node.x) $(node.y)\n")
    end
    close(process.in)

    # extract relevant lines from postscript output of GeoSteiner
    line_tokens = Array{String}[]
    _started = false
    for line in eachline(process.out)
        # are we in the good part yet?
        if !_started && !occursin(" % fs", line)
            continue
        else
            _started = true
        end

        # are we still in the good part?
        if occursin("Euclidean SMT", line)
            break
        end

        # skip the noise
        if occursin(" % fs", line)
            continue
        end

        # extract the values
        push!(line_tokens, split(line))
    end

    wait(process)

    for tokens in line_tokens
        @assert length(tokens) == 5
        @assert tokens[5] == "S"
        tail, head = 0, 0
        if tokens[2] == "T"
            tail = parse(Int, tokens[1]) + 1
        else
            stein = Node(parse(Float64, tokens[1]), parse(Float64, tokens[2]))
            tail = findfirst(isequal(stein), tnodes)
            if tail == nothing
                push!(tnodes, stein)
                tail = length(tnodes)
            end
        end
        if tokens[4] == "T"
            head = parse(Int, tokens[3]) + 1
        else
            stein = Node(parse(Float64, tokens[3]), parse(Float64, tokens[4]))
            head = findfirst(isequal(stein), tnodes)
            if head == nothing
                push!(tnodes, stein)
                head = length(tnodes)
            end
        end

        push!(tarcs, Arc(tail,head))
    end

    @assert(length(tnodes) <= 2*length(nodes) - 2)
    @assert(length(tarcs) == length(tnodes) - 1)

    Topology(tnodes, tarcs)
end


text:
import Base.ArgumentError

type Sequence{Tx,Ty}
    x::Array{Tx,1} # observation sequence
    y::Array{Ty,1} # label sequence
    Y::Array{Ty,1} # alphabet from which labels are drawn

    n::Int # length of training data
    d::Int # length of label alphabet
    k::Int # number of features

    Z::Float64 # normalization constant
    σ::Float64 # regularization constant

    M::Array{Array{Float64, 2}, 1} # transformation matrix
    α::Array{Array{Float64, 1}, 1} # forward vectors
    β::Array{Array{Float64, 1}, 1} # backward vectors

    Θ::Array{Float64, 1} # parameters aka feature weights

    f::Function # function returning feature vector
    F::Array{Float64, 1} # sum of feature vectors 1:n

    function Sequence{Tx,Ty}(x::AbstractArray{Tx},
                             y::AbstractArray{Ty},
                             f::Function;
                             labels::AbstractArray{Ty}=Ty[],
                             Θ::AbstractArray{Float64}=Float64[],
                             σ::Float64=50.0)

        # If x are unlabeled observations, the labels keyword must provide the
        # alphabet from which labels are drawn.
        if isempty(y)
            if isempty(labels)
                throw(ArgumentError("A label alphabet must be specified for unlabeled observations."))
            end
        else
            if length(x) != length(y)
                throw(ArgumentError("Labels don't match observation length."))
            end
        end

        res = new()
        res.x = x
        res.y = y

        if isempty(labels)
            res.Y = unique(y)
        else
            res.Y = labels
        end

        res.f = f
        res.σ = σ
        res.Z = 0.0

        res.n, = size(res.x)
        res.d, = size(res.Y)
        res.k, = size(res.f(res.Y[1], res.Y[1], res.x, 1))

        res.M = [ zeros(res.d, res.d) for t = 1:res.n ]
        res.α = [ zeros(res.d) for t = 1:res.n ]
        res.β = [ zeros(res.d) for t = 1:res.n ]

        if isempty(Θ)
            Θ = rand(res.k)
        else
            if length(Θ) != res.k
                throw(ArgumentError("Length of Θ doesn't match feature length."))
            end
        end

        res.Θ = zeros(res.k)
        res.F = zeros(res.k)
        update(res, Θ=Θ)

        return res
    end
end

Sequence{Tx,Ty}(x::AbstractArray{Tx}, y::AbstractArray{Ty}, f::Function; args...) = Sequence{Tx,Ty}(x, y, f; args...)
Sequence{Tx,Ty}(x::AbstractArray{Tx}, f::Function; labels::AbstractArray{Ty}=Ty[], args...) = Sequence{Tx,Ty}(x, Ty[], f; labels=labels, args...)

function update(crf::Sequence; Θ::AbstractArray{Float64}=Float64[])
    if !isempty(Θ)
        if crf.Θ == Θ
            return
        end
        for i = 1:crf.k
            crf.Θ[i] = Θ[i]
        end
    end

    # Transition Matrix
    for (t, mt) = enumerate(crf.M)
        for (i, yi) = enumerate(crf.Y), (j, yj) = enumerate(crf.Y)
            mt[i,j] = dot(crf.Θ, crf.f(yi, yj, crf.x, t))
        end
    end

    # Forward Vectors
    for t = 1:crf.n
        for (i, yi) = enumerate(crf.Y)
            if t == 1
                crf.α[t][i] = dot(crf.Θ, crf.f(yi, crf.x, t))
            else
                crf.α[t][i] = logsumexp(crf.M[t][:,i][:] + crf.α[t-1])
            end
        end
    end

    # Backward Vectors
    for t = reverse(1:crf.n)
        for (i, yi) = enumerate(crf.Y)
            if t == crf.n
                crf.β[t][i] = 0.0
            else
                crf.β[t][i] = logsumexp(crf.M[t+1][i,:][:] + crf.β[t+1])
            end
        end
    end

    # Constant normalization factor
    crf.Z = logsumexp(crf.α[end])

    # Calculate global feature vector if y is present
    if !isempty(crf.y)
        crf.F[:] = 0.0
        for t = 1:crf.n
            if t == 1
                v = crf.f(crf.y[t], crf.x, t)
            else
                v = crf.f(crf.y[t-1], crf.y[t], crf.x, t)
            end
            for i = 1:crf.k
                crf.F[i] += v[i]
            end
        end
    end
end

function update(crfs::AbstractArray{Sequence}; Θ::AbstractArray{Float64}=Float64[])
    if isempty(Θ)
        Θ = copy(crfs[1].Θ)
    end
    for crf in crfs
        update(crf, Θ=Θ)
    end
end

function loglikelihood(crf::Sequence; Θ::AbstractArray{Float64}=Float64[])
    if !isempty(Θ)
        update(crf, Θ=Θ)
    end
    return dot(crf.Θ, crf.F) - crf.Z - dot(crf.Θ, crf.Θ) / (crf.σ ^ 2)
end

function loglikelihood(crfs::AbstractArray{Sequence}; Θ::AbstractArray{Float64}=Float64[])
    sum([ loglikelihood(crf, Θ=Θ) for crf in crfs ])
end

function loglikelihood_gradient(crf::Sequence; Θ::AbstractArray{Float64}=Float64[])
    if !isempty(Θ) && (crf.Θ != Θ)
        update(crf, Θ=Θ)
    end
    F = copy(crf.F)
    t = 1
    for (i, yi) = enumerate(crf.Y)
        v = crf.f(yi, crf.x, t)
        p = exp(dot(crf.Θ, v) + crf.β[t][i] - crf.Z)

        for u = 1:crf.k
            F[u] -= v[u] * p
        end
    end
    for t = 2:crf.n
        for (i, yi) = enumerate(crf.Y), (j, yj) = enumerate(crf.Y)
            v = crf.f(yi, yj, crf.x, t)
            p = exp(crf.α[t-1][i] + crf.M[t][i,j] + crf.β[t][j] - crf.Z)

            for u = 1:crf.k
                F[u] -= v[u] * p
            end
        end
    end
    r = sum(crf.Θ) / crf.σ
    for i = 1:crf.k
        F[i] -= r
    end
    return F
end

function loglikelihood_gradient(crfs::AbstractArray{Sequence}; Θ::AbstractArray{Float64}=Float64[])
    sum([ loglikelihood_gradient(crf, Θ=Θ) for crf in crfs ])
end

function label(crf::Sequence)
    d = zeros(crf.d)
    Q = zeros(crf.d, crf.d)

    M = AbstractArray{Int}[]
    for t in 1:crf.n
        for i = 1:crf.d, j = 1:crf.d
            Q[i,j] = crf.M[t][i,j] + d[j]
        end
        m = zeros(Int, crf.d)
        for i = 1:crf.d
            j = indmax(Q[i,:])
            m[i] = j
            d[i] = Q[i,j]
        end
        push!(M, m)
    end

    result = Int[ indmax(d) ]
    for p in reverse(M)
        push!(result, p[result[end]])
    end
    return [ crf.Y[i] for i in reverse(result[2:end]) ]
end

function label(crfs::AbstractArray{Sequence})
    [ label(crf) for crf in crfs ]
end


text:
module ModelicaScriptingTools

using Base.Filesystem: isfile
using Test: @test
using CSV: CSV
using OMJulia: OMCSession, sendExpression, Parser
using ZMQ: send, recv # only needed for sendExpressionRaw which is a workaround for OMJulia bugs
import ZMQ # for accessing ZMQ.StateError
using DataFrames: DataFrame
using PyCall: PyNULL, pyimport_conda, pyimport, @py_str
using Markdown: Markdown
import Documenter

export moescape, mounescape, MoSTError, loadModel, getSimulationSettings,
    getVariableFilter, simulate, regressionTest, testmodel,
    setupOMCSession, closeOMCSession, withOMC, ModelicaBlocks, getDocAnnotation,
    getequations, getcode, getvariables, mdescape, variabletable, commonprefix,
    commonhierarchy, findvarnames, deprefix, aliasdict, findidentifiers,
    getfunctions, findfuncnames, functionlist, uniquehierarchy,
    replacefuncnames, uniquefunctions, getVersion, installAndLoad

include("Simulation.jl")
include("Testing.jl")
include("Documentation.jl")

function __init__()
    __init__doc()
end

end

max_stars_repo_name: THM-MoTE/MoST

max_issues_repo_name: THM-MoTE/MoST

max_forks_repo_name: THM-MoTE/MoST


text:
abstract type StaticIRNode end

struct TrainableParameterNode <: StaticIRNode
    name::Symbol
    typ::Union{Symbol,Expr,QuoteNode}
end

struct ArgumentNode <: StaticIRNode
    name::Symbol
    typ::Union{Symbol,Expr,QuoteNode}
    compute_grad::Bool
end

struct JuliaNode <: StaticIRNode
    fn::Function
    inputs::Vector{StaticIRNode}
    name::Symbol
    typ::Union{Symbol,Expr,QuoteNode}
end

struct RandomChoiceNode <: StaticIRNode
    dist::Distribution
    inputs::Vector{StaticIRNode}
    addr::Symbol
    name::Symbol
    typ::Union{Symbol,Expr,QuoteNode}
end

struct GenerativeFunctionCallNode <: StaticIRNode
    generative_function::GenerativeFunction
    inputs::Vector{StaticIRNode}
    addr::Symbol
    name::Symbol
    typ::Union{Symbol,Expr,QuoteNode}
end

struct StaticIR
    nodes::Vector{StaticIRNode}
    trainable_param_nodes::Vector{TrainableParameterNode}
    arg_nodes::Vector{ArgumentNode}
    choice_nodes::Vector{RandomChoiceNode}
    call_nodes::Vector{GenerativeFunctionCallNode}
    julia_nodes::Vector{JuliaNode}
    return_node::StaticIRNode
    accepts_output_grad::Bool
end

mutable struct StaticIRBuilder
    nodes::Vector{StaticIRNode}
    node_set::Set{StaticIRNode}
    trainable_param_nodes::Vector{TrainableParameterNode}
    arg_nodes::Vector{ArgumentNode}
    choice_nodes::Vector{RandomChoiceNode}
    call_nodes::Vector{GenerativeFunctionCallNode}
    julia_nodes::Vector{JuliaNode}
    return_node::Union{Nothing,StaticIRNode}
    vars::Set{Symbol}
    addrs_to_choice_nodes::Dict{Symbol,RandomChoiceNode}
    addrs_to_call_nodes::Dict{Symbol,GenerativeFunctionCallNode}
    accepts_output_grad::Bool
end

function StaticIRBuilder()
    nodes = Vector{StaticIRNode}()
    node_set = Set{StaticIRNode}()
    trainable_param_nodes = Vector{TrainableParameterNode}()
    arg_nodes = Vector{ArgumentNode}()
    choice_nodes = Vector{RandomChoiceNode}()
    call_nodes = Vector{GenerativeFunctionCallNode}()
    julia_nodes = Vector{JuliaNode}()
    return_node = nothing
    vars = Set{Symbol}()
    addrs_to_choice_nodes = Dict{Symbol,RandomChoiceNode}()
    addrs_to_call_nodes = Dict{Symbol,GenerativeFunctionCallNode}()
    accepts_output_grad = false
    StaticIRBuilder(nodes, node_set, trainable_param_nodes, arg_nodes, choice_nodes, call_nodes,
        julia_nodes,
        return_node, vars, addrs_to_choice_nodes, addrs_to_call_nodes,
        accepts_output_grad)
end

function build_ir(builder::StaticIRBuilder)
    if builder.return_node === nothing
        builder.return_node = add_constant_node!(builder, nothing)
    end
    StaticIR(
        builder.nodes,
        builder.trainable_param_nodes,
        builder.arg_nodes,
        builder.choice_nodes,
        builder.call_nodes,
        builder.julia_nodes,
        builder.return_node,
        builder.accepts_output_grad)
end

function check_unique_var(builder::StaticIRBuilder, name::Symbol)
    if name in builder.vars
        error("Variable name $name is not unique")
    end
end

function check_inputs_exist(builder::StaticIRBuilder, input_nodes)
    for input_node in input_nodes
        if !(input_node in builder.node_set)
            error("Node $input_node was not previously added to the IR")
        end
    end
end

function check_addr_unique(builder::StaticIRBuilder, addr::Symbol)
    if haskey(builder.addrs_to_choice_nodes, addr) || haskey(builder.addrs_to_call_nodes, addr)
        error("Address $addr was not unique")
    end
end

function _add_node!(builder::StaticIRBuilder, node::StaticIRNode)
    push!(builder.nodes, node)
    push!(builder.node_set, node)
end

function add_trainable_param_node!(builder::StaticIRBuilder,
                                   name::Symbol;
                                   typ::Union{Symbol,Expr,QuoteNode}=QuoteNode(Any))
    check_unique_var(builder, name)
    node = TrainableParameterNode(name, typ)
    _add_node!(builder, node)
    push!(builder.trainable_param_nodes, node)
    node
end

function add_argument_node!(builder::StaticIRBuilder;
                            name::Symbol=gensym(),
                            typ::Union{Symbol,Expr,QuoteNode}=QuoteNode(Any),
                            compute_grad=false)
    check_unique_var(builder, name)
    node = ArgumentNode(name, typ, compute_grad)
    _add_node!(builder, node)
    push!(builder.arg_nodes, node)
    node
end

function add_julia_node!(builder::StaticIRBuilder, fn::Function;
                         inputs::Vector=[],
                         name::Symbol=gensym(),
                         typ::Union{Symbol,Expr,QuoteNode}=QuoteNode(Any))
    check_unique_var(builder, name)
    check_inputs_exist(builder, inputs)
    node = JuliaNode(fn, inputs, name, typ)
    _add_node!(builder, node)
    push!(builder.julia_nodes, node)
    node
end

function add_constant_node!(builder::StaticIRBuilder, val,
        name::Symbol=gensym(),
        typ::Union{Symbol,Expr,QuoteNode}=QuoteNode(typeof(val)))
    check_unique_var(builder, name)
    # NOTE: not wrapping it in a Diffed means it is interpreted as a constant
    node = JuliaNode(() -> val, [], name, typ)
    _add_node!(builder, node)
    push!(builder.julia_nodes, node)
    node
end

function add_addr_node!(builder::StaticIRBuilder, dist::Distribution;
                        inputs::Vector=[], addr::Symbol=gensym(),
                        name::Symbol=gensym())
    check_unique_var(builder, name)
    check_addr_unique(builder, addr)
    check_inputs_exist(builder, inputs)
    typ = QuoteNode(get_return_type(dist))
    node = RandomChoiceNode(dist, inputs, addr, name, typ)
    _add_node!(builder, node)
    builder.addrs_to_choice_nodes[addr] = node
    push!(builder.choice_nodes, node)
    node
end

function add_addr_node!(builder::StaticIRBuilder, gen_fn::GenerativeFunction;
                        inputs::Vector=[], addr::Symbol=gensym(),
                        name::Symbol=gensym())
    check_unique_var(builder, name)
    check_addr_unique(builder, addr)
    check_inputs_exist(builder, inputs)
    typ = QuoteNode(get_return_type(gen_fn))
    node = GenerativeFunctionCallNode(gen_fn, inputs, addr, name, typ)
    _add_node!(builder, node)
    builder.addrs_to_call_nodes[addr] = node
    push!(builder.call_nodes, node)
    node
end

function set_return_node!(builder::StaticIRBuilder, node::StaticIRNode)
    if builder.return_node !== nothing
        error("Return node already set")
    end
    builder.return_node = node
    nothing
end

function set_accepts_output_grad!(builder::StaticIRBuilder, value::Bool)
    builder.accepts_output_grad = value
end

export StaticIR, StaticIRBuilder, build_ir
export add_trainable_param_node!
export add_argument_node!
export add_julia_node!
export add_constant_node!
export add_addr_node!
export set_return_node!


text:
# using Compose

using Compose
compose(context(), fill("red"), circle(0.5,0.5,0.2))

# the x, y, r of the three circles:
xs = [0.25, 0.5, 0.75];
ys = [0.75, 0.5, 0.25];
rs = [0.1, 0.1, 0.1]; # or simply [0.1] is ok

# the colors to fill for the three circles
using Colors
colors = [LCHab(92, 10, 77), LCHab(68, 74, 192), LCHab(78, 84, 29)];

fig_e = compose(context(),
circle(xs, ys, rs), fill(colors));

function resultandtree(fig_x)
  yh = 0.8;
  cap(x,y,str,clr) = (context(), text(x,y,str,hcenter,vtop),fill(clr));
  c1 = cap(0.5,0,"result","black");
  c2 = cap(0.5,0,"introspect","white");
  compose(context(0,0,1,yh),
          (context(0,0,0.5,yh),c1, fig_x),
    	  (context(0.5,0,0.5,yh),c2, introspect(fig_x)))
end

resultandtree(fig_e)

# the start points of the 9 sub-coordinate systems
x0 = [0,0,0,1,1,1,2,2,2]./3;
y0 = [0,1,2,0,1,2,0,1,2]./3;
# and their x range and y range
xw = yh = 1.0./3;

subcoord(x0,y0) = (context(x0,y0,xw,yh), fig_e);
fig_f = compose(context(),
        subcoord.(x0,y0)... );

resultandtree(fig_f)

subcoord(x0,y0) = (context(x0,y0,xw,yh;units=UnitBox(0,0,2,2)), fig_e);
fig_g = compose(context(),
        subcoord.(x0,y0)... );

resultandtree(fig_g)

h = 0.7;
w = 1/3;
box = (context(), rectangle(), fill(nothing),stroke("black"));
cap(str) = (context(), text(0.75,0.95,str),fill("black"));
compose(context(0.0,0.0,1.0,h),
         (context(0.0,0.0,w,h), cap("fig_e"), fig_e, box),
		 (context(1/3,0.0,w,h), cap("fig_f"), fig_f, box),
		 (context(2/3,0.0,w,h), cap("fig_g"), fig_g, box))

v1 = (0.5,0.5); # vector
absv = hypot(v1...); # length of vector
clr = stroke(RGBA(90/255,39/255,41/255,absv)) # color of arrow

"""
`unitarrow(arrangle, arrlength)` creates a good arrow. Noted that the maximum length of arrow should not exceeds 1.
"""
function unitarrow(arrangle, arrlength; arrwidth=0.1, arrheight=0.15)
	head_x =0.5 + 0.5*arrlength;
	center = (0.5, 0.5);
	headpoint = (head_x, 0.5);
	sidepoint_a = headpoint.+(-arrheight,-arrwidth/2);
	sidepoint_b = headpoint.+(-arrheight,arrwidth/2);
	function lim_x(sidepoint_a)
		sidepoint_a = maximum.(([sidepoint_a[1], center[1]], [sidepoint_a[2]]))
	end
	sidepoint_a = lim_x(sidepoint_a);
	sidepoint_b = lim_x(sidepoint_b);

	eq_triangle = [sidepoint_a,
				   sidepoint_b,
				   headpoint]
	tailwidth = arrwidth/5;
	arrowtail = (context(), polygon(
			[(1-headpoint[1],headpoint[2]-tailwidth/2),
			 (1-headpoint[1],headpoint[2]+tailwidth/2),
			 (sidepoint_a[1], headpoint[2]+tailwidth/2),
			 (sidepoint_a[1], headpoint[2]-tailwidth/2)
				]))
		# rectangle(headpoint[1]-tailwidth/2,0,tailwidth,1-arrheight))
	arrowheadx = (context(units=UnitBox(0,0,1,1)),polygon(eq_triangle));
	arw = compose(
		context(units=UnitBox(0,0,1,1), rotation=Rotation(arrangle, 0.5, 0.5)),
		arrowheadx,
		arrowtail);
end

theta = atan(v1[2], v1[1])-pi/2;
unitarrow(theta, absv)

function mymeshgrid(rangex, rangey)
  xgrid = Float64[];
  ygrid = Float64[];
  for i in rangex
    for j in rangey
      push!(xgrid, i);
		push!(ygrid, j);
    end
  end
  return (xgrid, ygrid)
end
subcoord2(x0,y0,obj) = (context(x0,y0,1,1), obj);
subcoord1(x0,y0) = (context(x0,y0,1,1), rectangle(), stroke("black"), fill(nothing));

n = 5;
(xs, ys) = mymeshgrid(collect(0:n),collect(0:n));
m = length(xs);
thetas = rand(1:360, m) ./ (2pi);
arrlen = rand(1:100, m) ./ 100;
objs = unitarrow.(thetas, arrlen);
compose(context(;units=UnitBox(0,0,n+1,n+1)),
        subcoord2.(xs,ys,objs)..., subcoord1.(xs,ys)...)

function arrowhead(θ)
	eq_triangle = [(0, 1/sqrt(3)),
		           (-1/3, -2/(2 * sqrt(3))),
		           (1/3, -2/(2 * sqrt(3)))]

	compose(context(units=UnitBox(-1,-1,2,2), rotation=Rotation(θ, 0, 0)),
				polygon(eq_triangle))
	# Rotation(θ, x, y): Rotate all forms in context around point (x,y) by angle θ in radians.
end

function quiver(points, vecs)
	xmin = minimum(first.(points))
	ymin = minimum(last.(points))
	xmax = maximum(first.(points))
	ymax = maximum(last.(points))
	hs = map(x->hypot(x...), vecs)
	hs = hs / maximum(hs)

	vector(p, v, h) = all(iszero, v) ? context() :
		(context(),
		    (context((p.+v.*6 .- .2)..., .4,.4),
				arrowhead(atan(v[2], v[1]) - pi/2)),
		stroke(RGBA(90/255,39/255,41/255,h)),
		fill(RGBA(90/255,39/255,41/255,h)),
		line([p, p.+v.*8]))

	compose(context(units=UnitBox(xmin,ymin,xmax,ymax)),
         vector.(points, vecs, hs)...)
end

quiver([(j-1,i-1) for i=0:n, j=0:n],
       [len.*(cos(theta), sin(theta)) for theta in thetas, len in arrlen])

max_stars_repo_name: okatsn/UnpredictableLife

max_issues_repo_name: okatsn/UnpredictableLife

max_forks_repo_name: okatsn/UnpredictableLife


text:
using LightGraphs

input = joinpath(@__DIR__, "input")
caverns = readlines(input)

struct Coord
    x::Int
    y::Int
end

function adjacent(coord::Coord)
    x, y = coord.x, coord.y
    return [
        Coord(x, y + 1),
        Coord(x, y - 1),
        Coord(x + 1, y),
        Coord(x - 1, y)
    ]
end

function map_caverns!(
    cs::Array{String,1},
    origins::Array{Int,1},
    c2v::Dict{Coord,Int},
    c2t::Dict{Coord,Char},
    v2c::Dict{Int,Coord},
    v2t::Dict{Int,Char}
)
    n = 1
    for j = 1:length(cs)
        line = cs[j]
        for i = 1:length(line)
            t = line[i]
            c2v[Coord(i, j)] = n
            c2t[Coord(i, j)] = t
            v2c[n] = Coord(i, j)
            v2t[n] = t
            ## capture any starting points
            if t == '@'
                push!(origins, n)
            end
            n += 1
        end
    end
    return n
end

function paths!(g::SimpleGraph{Int}, c2v::Dict{Coord,Int}, v2t::Dict{Int,Char})
    for loc in collect(keys(c2v))
        adj = adjacent(loc)
        for a in adj
            if haskey(c2v, a)
                src = c2v[loc]
                dst = c2v[a]
                if v2t[src] != '#' && v2t[dst] != '#'
                    add_edge!(g, src, dst)
                end
            end
        end
    end
end

function path_to_key!(cache::Dict{Pair{Int,Int},Array{Int,1}}, g::SimpleGraph{Int}, src::Int, dst::Int)
    k = Pair(src, dst)
    if haskey(cache, k)
        return cache[k]
    else
        path = map(e -> e.dst, a_star(g, src, dst))
        cache[k] = path
        return path
    end
end

function reachable(path::Array{Int,1}, have_keys::Set{Int}, v2t::Dict{Int,Char})
    target = v2t[last(path)]
    for v in path
        c = v2t[v]
        ## another key in path, don't bother exploring
        if c >= 'a' && c <= 'z' && c != target && !(v in have_keys)
            return false
        ## door in path, only continue if unlockable
        elseif c >= 'A' && c <= 'Z'
            ks = [uppercase(v2t[n]) for n in have_keys]
            if !(c in ks)
                return false
            end
        end
    end
    return true
end

function explore(g::SimpleGraph{Int}, origins::Array{Int,1}, v2t::Dict{Int,Char})
    nodes = [(Set{Int}(origins), 0, Set{Int}())]
    seen = Set{Tuple{Set{Int},Int,Set{Int}}}()

    ks = collect(keys(filter(vt -> vt.second >= 'a' && vt.second <= 'z', v2t)))

    # step count of shortest complete journey
    shortest_so_far = Inf

    # cache of paths
    cache = Dict{Pair{Int,Int},Array{Int,1}}()

    while !isempty(nodes)
        node = popfirst!(nodes)
        # skip exploring if we've already started from this node
        if node in seen
            continue
        else
            push!(seen, node)
            vertices, steps, have_keys = node
            ## have all keys, update `shortest_so_far` and stop exploring from
            ## here
            if length(have_keys) == length(ks)
                if steps < shortest_so_far
                    shortest_so_far = steps
                end
                continue
            end
            ## trip is longer than an already completed trip, don't explore
            ## further from here
            if steps > shortest_so_far
                continue
            end
            ## for each robot, create nodes for all possible moves
            for v in vertices
                for k in ks
                    if !(k in have_keys)
                        p = path_to_key!(cache, g, v, k)
                        if length(p) > 0 && reachable(p, have_keys, v2t)
                            new_vertices = copy(vertices)
                            delete!(new_vertices, v)
                            push!(new_vertices, k)
                            new_have_keys = copy(have_keys)
                            push!(new_have_keys, k)
                            new = (new_vertices, steps + length(p), new_have_keys)
                            if !(new in seen)
                                push!(nodes, new)
                            end
                        end
                    end
                end
            end
        end
    end

    return shortest_so_far
end

## starting point(s)
origins = Array{Int}(undef, (0))
## store coords to vertex id, coords to type, vertex id to coords and vertex id to type
c2v = Dict{Coord,Int}()
c2t = Dict{Coord,Char}()
v2c = Dict{Int,Coord}()
v2t = Dict{Int,Char}()
## vs is total number of nodes
vs = map_caverns!(caverns, origins, c2v, c2t, v2c, v2t)

g = SimpleGraph(vs)
paths!(g, c2v, v2t)

p1 = explore(g, origins, v2t)
# @assert p1 == 4762

## get old entrance coords
c = v2c[first(origins)]

## reset everything for p2

## store coords to vertex id, coords to type, vertex id to coords and vertex id to type
c2v = Dict{Coord,Int}()
c2t = Dict{Coord,Char}()
v2c = Dict{Int,Coord}()
v2t = Dict{Int,Char}()

## vs is total number of nodes
vs = map_caverns!(caverns, origins, c2v, c2t, v2c, v2t)

## mess with map for p2
v2t[c2v[c]] = '#'
for coord in adjacent(c)
    v = c2v[coord]
    v2t[v] = '#'
end
origins = Array{Int}(undef, (0))
for coord in [Coord(c.x + 1, c.y + 1), Coord(c.x + 1, c.y - 1), Coord(c.x - 1, c.y + 1), Coord(c.x - 1, c.y - 1)]
    v = c2v[coord]
    v2t[v] = '@'
    push!(origins, v)
end

h = SimpleGraph(vs)
paths!(h, c2v, v2t)

p2 = explore(h, origins, v2t)
@show p2
@assert p2 == 1876

print("-----------------------------------------------------------------------\n")
print("many worlds interpretation -- part one\n    distance : $p1\n")
print("many worlds interpretation -- part two\n    distance : $p2\n")
print("-----------------------------------------------------------------------\n")


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Libgpg_error_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Libgpg_error")
JLLWrappers.@generate_main_file("Libgpg_error", UUID("7add5ba3-2f88-524e-9cd5-f83b8a55f7b8"))
end  # module Libgpg_error_jll


text:
function Kmeans(X,k;MaxIter=20)
    m=size(X,1)
    n=size(X,2)

    #Initial Centroids
    C=X[rand(1:m,k),:]
    Idx=zeros(m,MaxIter)
    for i=1:MaxIter
       Idx[:,i]=FindCentroids(X,C)
        C=ComputeCentroids(X,Idx[:,i])
    end
    Id=round.(Int,Idx[:,end])
    return Id,C
end


text:
#VERSION >= v"0.4" #&& __precompile__()

module JuliaTwitter

#Must remain in this order
include("types/TwContributor.jl")
include("types/TwCoordinate.jl")
include("types/TwEntity.jl")
include("types/TwUser.jl")
include("types/Tweet.jl")

using OAuth, Requests, Compat

export
  TwAuth,
  TwUser,
  Tweet,
  TwCoordinate,
  TwContributor,
  auth,
  get_oauth

type TwAuth
  consumer_key::ASCIIString
  consumer_secret::ASCIIString
  access_token::ASCIIString
  access_secret::ASCIIString
end

# A few testing functions, these will be organized into files later
function get_oauth(auth::TwAuth, url::ASCIIString, options::Dict)
  oauth_request_resource(url, "GET", options, auth.consumer_key, auth.consumer_secret, auth.access_token, auth.access_secret)
end

#Include REST APIs
include("rest/get_account.jl")
include("rest/get_application.jl")
include("rest/get_blocks.jl")
include("rest/get_direct_messages.jl")
include("rest/get_favorites.jl")
include("rest/get_followers.jl")
include("rest/get_friends.jl")
include("rest/get_friendships.jl")
include("rest/get_geo.jl")
include("rest/get_help.jl")
include("rest/get_lists.jl")
include("rest/get_mutes.jl")
include("rest/get_saved_searches.jl")
include("rest/get_search.jl")
include("rest/get_statuses.jl")
include("rest/get_trends.jl")
include("rest/get_users.jl")

#Export REST APIs
export
  get_account_settings,
  get_account_verify_credentials,
  get_application_rate_limit_status,
  get_blocks_ids,
  get_blocks_list,
  get_dm,
  get_dm_sent,
  get_dm_show,
  get_favorites_list,
  get_followers_ids,
  get_followers_list,
  get_friends_ids,
  get_friends_list,
  get_friendships_incoming,
  get_friendships_lookup,
  get_friendships_no_retweets_ids,
  get_friendships_outgoing,
  get_friendships_show,
  get_geo_id_place_id,
  get_geo_reverse_geocode,
  get_geo_search,
  get_help_configuration,
  get_help_languages,
  get_help_privacy,
  get_help_tos,
  get_lists_list,
  get_lists_members,
  get_lists_members_show,
  get_lists_memberships,
  get_lists_ownerships,
  get_lists_show,
  get_lists_statuses,
  get_lists_subscribers,
  get_lists_subscribers_show,
  get_lists_subscriptions,
  get_mutes_users_ids,
  get_mutes_users_list,
  get_saved_searches,
  get_saved_searches_show_id,
  get_search,
  get_searchq,
  get_statuses,
  get_statuses_home_timeline,
  get_statuses_lookup,
  get_statuses_mentions_timeline,
  get_statuses_oembed,
  get_statuses_retweets_id,
  get_statuses_retweets_of_me,
  get_statuses_retweeters_ids,
  get_statuses_show_id,
  get_statuses_user_timeline,
  get_trends_available,
  get_trends_closest,
  get_trends_place,
  get_users,
  get_users_lookup,
  get_users_profile_banner,
  get_users_search,
  get_users_show,
  get_users_suggestions,
  get_users_suggestions_slug,
  get_users_suggestions_slug_members

#Streaming
include("streaming/streaming.jl")

#Include all active dbhandles
#include("../dbhandle/active.jl")

#Include User Config File
include("../user_config.jl")

auth = TwAuth(USER_AUTH_DICT[:consumer_key],
              USER_AUTH_DICT[:consumer_secret],
              USER_AUTH_DICT[:access_token],
              USER_AUTH_DICT[:access_secret])

export TwitterDB, stream_tweets

end

max_stars_repo_name: nick-paul/JuliaTwitter

max_issues_repo_name: nick-paul/JuliaTwitter

max_forks_repo_name: nick-paul/JuliaTwitter


text:
#
# readoutModesAvailability.jl -
#
# This example demonstrates how to list the readout modes available on the
# camera.  The user is then able to able to select one, and grab an image with
# it.
#
# Please refer to the "simpleAcquisition.jl" example if you're looking for a
# simpler example to acquire images.
#

module ReadoutModesAvailability

# In this example, we want to use the low-level interface as it is the closest
# to the C API.
import NuvuCameras: NC

function main()

    # Opens the acquisition channel using the automatic detection and 4 loop
    # buffers (recommanded).
    myCam = NC.open(NC.Cam, NC.AUTO_UNIT, NC.AUTO_CHANNEL, 4)

    # Investigates about the readout modes available on the camera.
    nmodes = NC.getNbrReadoutModes(myCam)
    println("Here are the available readout modes:\n")
    for i in 1:nmodes
        amptyp, ampstr, vfreq, hfreq = NC.getReadoutMode(myCam, i)
        println(" > Readout mode number $i:")
 	println("   Amplifier type: $ampstr")
        println("   Vertical frequency: $(vfreq)Hz")
        println("   Horizontal frequency: $(hfreq)Hz\n")
    end

    # Second method to inquire the available readout modes, doing so per Ampli type.
    for (typ, descr) in ((NC.CONV, "conventionnal"),
                         (NC.EM,   "EM (electron multiplying)"))
        namplis = NC.getAmpliTypeAvail(myCam, typ)
        if namplis > 0
            println("\nThe $descr amplifier type is supported.")
 	    println("The frequencies supported for this amplifier type are:")
 	    for i in 0:namplis-1
                vfreq, hfreq, modenum = NC.getFreqAvail(myCam, NC.CONV, i)
                println(" > Readout mode number $modenum:")
                println("   Vertical frequency: $(vfreq)Hz")
                println("   Horizontal frequency: $(hfreq)Hz\n")
 	    end
        end
    end

    print("\nPlease select the desired mode you want to load on the camera: ")
    flush(STDOUT)
    modenum = parse(Cint, readline(STDIN))

    NC.setReadoutMode(myCam, modenum)

    println("\nThe readout mode has been successfully loaded on the camera.")
    println("We will now do a dummy grab on the camera and save it to validate it.")
    print("Hit the Return key to start...")
    flush(STDOUT)
    readline(STDIN)

    # Open the shutter for the acquisition.
    NC.setShutterMode(myCam, NC.OPEN)

    # Launches an acquisition by the framegrabber and requests an image from
    # the camera (this function does not wait for the acquisition to be
    # complete before returning).
    NC.start(myCam, 1)

    # Reads the image received, if a timeout occurs an error code will be
    # returned.
    myImage = NC.read(myCam)

    # Close the shutter, now that the acquistion is complete.
    NC.setShutterMode(myCam, NC.CLOSE)

    # Saves the image acquired
    NC.saveImage(myCam, myNcImage, "NewReadoutImage", NC.TIF,
                 "This is an image grab with the readout mode selected", true)

    println("You just acquired one image, congratulation!.\n")

    # Closes the acquisition channel no longer in use.
    NC.CamClose(myCam)

    if is_windows()
	print("\nHit the Return key to close the command window.")
        flush(STDOUT)
        readline(STDIN)
    else
	println("The program ended.")
    end
end

end # module


text:
module CommunalHelperCassetteJumpFixController

using ..Ahorn, Maple
using Ahorn.CommunalHelper

@mapdef Entity "CommunalHelper/CassetteJumpFixController" CassetteJumpFixController(
    x::Integer,
    y::Integer,
    persistent::Bool=false,
    off::Bool=false,
)

const placements = Ahorn.PlacementDict(
    "Cassette Jump Fix Controller (Communal Helper)" => Ahorn.EntityPlacement(
        CassetteJumpFixController,
    ),
)

const sprite = "objects/CommunalHelper/cassetteJumpFixController/icon"

function Ahorn.selection(entity::CassetteJumpFixController)
    x, y = Ahorn.position(entity)
    return Ahorn.getSpriteRectangle(sprite, x, y)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::CassetteJumpFixController)
    Ahorn.drawSprite(ctx, sprite, 0, 0)
end

end


text:
fname = "winequality.csv"
data = readdlm(fname, ';')
# 1600x12 Array{Any,2}:
#    "fixed acidity"   "volatile acidity"  …    "alcohol"   "quality"
#   7.4               0.7                      9.4         5.0
#   7.8               0.88                     9.8         5.0
#   7.8               0.76                     9.8         5.0
#  11.2               0.28                     9.8         6.0

# if comma (,) separated:
# data2 = readcsv(fname)
# variant:
data3 = readdlm(fname, ';', Float64, '\n', header=true)
#> (1599x12 Array{Float64,2}, 1x12 Array{String,2})
#  7.4  0.7    0.0   1.9  0.076  11.0  34.0  0.9978   3.51  0.56   9.4  5.0
#  7.8  0.88   0.0   2.6  0.098  25.0  67.0  0.9968   3.2   0.68   9.8  5.0
#  7.8  0.76   0.04  2.3  0.092  15.0  54.0  0.997    3.26  0.65   9.8  5.0
# 11.2  0.28   0.56  1.9  0.075  17.0  60.0  0.998    3.16  0.58   9.8  6.0
# ...
#  "fixed acidity"  "volatile acidity"  …  "sulphates"  "alcohol"  "quality"

# extracting data:
row3 = data[3, :] # measurements for wine xyz
# 1x12 Array{Any,2}:
#  7.8  0.88  0.0  2.6  0.098  25.0  67.0  0.9968  3.2  0.68  9.8  5.0
col3 = data[ :, 3]  # measurements of citric acid
# 1600-element Array{Any,1}:
# "citric acid"
# 0.0 0.0 0.04 0.56 0.0 0.0 ⋮ 0.08 0.08 0.1 0.13 0.12 0.47
x = data[:, 2:4] # measurements for "volatile acidity"   "citric acid"   "residual sugar"
y = data[70:75, 2:4]
# 6x3 Array{Any,2}:
 # 0.32   0.57  2.0
 # 0.705  0.05  1.9
 # 0.63   0.08  1.9
 # 0.67   0.23  2.1
 # 0.69   0.22  1.9
 # 0.675  0.26  2.1
 z = [data[:,3] data[:,6] data[:,11]]
 # 1600x3 Array{Any,2}:
 #  "citric acid"    "free sulfur dioxide"    "alcohol"
 # 0.0             11.0                      9.4
 # 0.0             25.0                      9.8
 # 0.04            15.0                      9.8
 # 0.56            17.0                      9.8
 # ⋮
 # 0.08            32.0                     10.5
 # 0.1             39.0                     11.2
 # 0.13            29.0                     11.0
 # 0.12            32.0                     10.2
 # 0.47            18.0                     11.0
 type Wine
    fixed_acidity::Array{Float64}
    volatile_acidity::Array{Float64}
    citric_acid::Array{Float64}
    #...
    quality::Array{Float64}
end
wine1 = Wine(7.4, 0.7, 0, ..., 5)
wine1 = Wine(data[1, :]...) # ... = splat-operator

# writing to csv file:
writedlm("partial.dat", z, ';')

max_stars_repo_path: Module 1/Chapter08/Old/csv_files.jl

max_issues_repo_path: Module 1/Chapter08/Old/csv_files.jl

max_forks_repo_path: Module 1/Chapter08/Old/csv_files.jl


text:
mutable struct MixintprogSolution
    status
    objval
    sol
    attrs
end

const SymbolInputVector = Union{Vector{Symbol},Symbol}

function no_mip_solver()
    error("""
          An integer programming solver must be specified as the final argument.

          You can use any solver listed on the solvers table of http://www.juliaopt.org/ that has checkmarks in both the Linear/Quadratic and Integer columns.

          A (free) example is Cbc.jl. Once Cbc is installed and imported via "using Cbc", you can specify CbcSolver() as the solver.
          Solver options are specified by using keyword arguments to CbcSolver().
          """)
end

function mixintprog(c::InputVector, A::AbstractMatrix, rowlb::InputVector, rowub::InputVector, vartypes::SymbolInputVector, lb::InputVector, ub::InputVector, solver::AbstractMathProgSolver)
    m = LinearQuadraticModel(solver)
    nrow,ncol = size(A)

    c = expandvec(c, ncol)
    rowlbtmp = expandvec(rowlb, nrow)
    rowubtmp = expandvec(rowub, nrow)
    lb = expandvec(lb, ncol)
    ub = expandvec(ub, ncol)
    vartypes = expandvec(vartypes, ncol)

    # rowlb is allowed to be vector of senses
    if eltype(rowlbtmp) == Char
        realtype = eltype(rowubtmp)
        sense = rowlbtmp
        rhs = rowubtmp
        @assert realtype <: Real
        rowlb = Array{realtype}(undef, nrow)
        rowub = Array{realtype}(undef, nrow)
        for i in 1:nrow
            if sense[i] == '<'
                rowlb[i] = typemin(realtype)
                rowub[i] = rhs[i]
            elseif sense[i] == '>'
                rowlb[i] = rhs[i]
                rowub[i] = typemax(realtype)
            elseif sense[i] == '='
                rowlb[i] = rhs[i]
                rowub[i] = rhs[i]
            else
                error("Unrecognized sense '$(sense[i])'")
            end
        end
    else
        rowlb = rowlbtmp
        rowub = rowubtmp
    end

    loadproblem!(m, A, lb, ub, c, rowlb, rowub, :Min)
    setvartype!(m, vartypes)
    optimize!(m)
    stat = status(m)
    if stat == :Optimal
        attrs = Dict()
        attrs[:objbound] = getobjbound(m)
        return MixintprogSolution(stat, getobjval(m), getsolution(m), attrs)
    else
        return MixintprogSolution(stat, nothing, [], Dict())
    end
end

mixintprog(c,A,rowlb,rowub,vartypes,solver::AbstractMathProgSolver) = mixintprog(c,A,rowlb,rowub,vartypes,0,Inf,solver)

# Old versions
mixintprog(c,A,rowlb,rowub,vartypes) = no_mip_solver()
mixintprog(c::InputVector, A::AbstractMatrix, rowlb::InputVector, rowub::InputVector, vartypes::SymbolInputVector, lb::InputVector, ub::InputVector) = no_mip_solver()

export mixintprog


text:
# include various methods for updating nodes
strategies = [
    "conjugacy",
    "opt",
    "advi"
]

for sname in strategies
    include(joinpath("inference", "$(sname).jl"))
end


text:
import Base: show

"""
$(TYPEDEF)
The context for a session, made from a user, robot, and session.
Users can have multiple robots and robots can have multiple sessions.
So this indicates a unique session.
"""
struct Client
    userId::String
    robotId::String
    sessionId::String
end

function show(io::IO, c::Client)
    print(io, "Client: User=$(c.userId), Robot=$(c.robotId), Session=$(c.sessionId)")
end

"""
$(TYPEDEF)
Some calls interact across multiple users, robots, and sessions.
A scope allows you to specify these more complex contexts.
"""
struct Scope
    environmentIds::Vector{String}
    userIds::Vector{String}
    robotIds::Vector{String}
    sessionIds::Vector{String}
end


text:
py"""
class DiffBenTal(torch.nn.Module):
    '''
    DiffBenTal

Differentiable BenTal model.

Atributes:
- `D::Int`: dimension (i.e. number of assets)
- `w::PyObject`: differentiable weights
- `c::Vector{PyObject}`: convex constraints on weights
- `μ::Array{Float64,1}`: predicted mean vector
- `sqrtΣ::Array{Float64,2}`: predicted upper triangular Cholesky factor of covariance matrix
- `solver_args::Dict=Dict("use_indirect" => true, "eps" => 1e-5, "max_iters" => 20000}`: solver arguments
    '''
    def __init__(
        self,
        D,
        w,
        c,
        μ,
        sqrtΣ,
        solver_args={'use_indirect': True, 'eps': 1e-5, 'max_iters': 20000}
    ):
        super(DiffBenTal, self).__init__()
        self.D = D
        self.μ = μ
        Σ = sqrtΣ.T @ sqrtΣ
        self.Σ = Σ
        self.solver_args=solver_args

        θ = cp.Variable()

        params = cp.Parameter(1)

        constraints = [
            cp.SOC(θ, sqrtΣ @ w),
        ]

        objective = cp.Maximize(μ.T @ w - θ * params[0])

        problem = cp.Problem(objective, constraints + c)

        self.cvxpylayer = CvxpyLayer(
            problem,
            parameters=[params],
            variables=[w, θ]
        )

    def forward(self, params):
        return self.cvxpylayer(params, solver_args=self.solver_args)[0]

    def get_variables(self, params):
        params_ = torch.tensor(params, dtype=torch.float64, requires_grad=True)
        return self.cvxpylayer(params_, solver_args=self.solver_args)

    def get_w(self, params):
        params_ = torch.tensor(params, dtype=torch.float64, requires_grad=True)
        return self.get_variables(params_)[0]

    def get_lagrangian(self, params):
        params_ = torch.tensor(params, dtype=torch.float64, requires_grad=True)
        D = self.D
        μ_ = torch.tensor(self.μ, dtype=torch.float64)
        Σ_ = torch.tensor(self.Σ, dtype=torch.float64)
        w_, θ_ = self.cvxpylayer(params_, solver_args=self.solver_args)
        return μ_.T @ w_ - θ_ * params_[0]
"""

DiffBenTal = py"DiffBenTal"

max_stars_repo_name: molet/DiffDRO

max_issues_repo_name: molet/DiffDRO

max_forks_repo_name: molet/DiffDRO


text:
using Tensors
using BenchmarkTools

abstract type AbstractMaterialState end

@generated function Base.:+(state::T, dstate::T) where {T <: AbstractMaterialState}
   expr = [:(state.$p+ dstate.$p) for p in fieldnames(T)]
   return :(T($(expr...)))
end

struct SomeState <: AbstractMaterialState
   stress::SymmetricTensor{2,3,Float64}
end

state = SomeState(SymmetricTensor{2,3,Float64}([1.0, 2.0, 3.0, 4.0, 5.0, 6.0]))

N = 1000
function bench_state(N)
   state = SomeState(SymmetricTensor{2,3,Float64}([1.0, 2.0, 3.0, 4.0, 5.0, 6.0]))
   for i in 1:N
      dstate = SomeState(randn(SymmetricTensor{2,3,Float64}))
      state = state + dstate
   end
   return state
end

# println("Benchmark State{SymmetricTensor{2,3,Float64}}")
# @btime bench_state(N)

struct AnotherState <: AbstractMaterialState
   stress::SymmetricTensor{2,3,Float64,6}
   strain::SymmetricTensor{2,3,Float64,6}
   backstress1::SymmetricTensor{2,3,Float64,6}
   backstress2::SymmetricTensor{2,3,Float64,6}
   plastic_strain::SymmetricTensor{2,3,Float64,6}
   cumeq::Float64
   R::Float64
end
 function bench_chaboche_style_state(N)
   stress = zero(SymmetricTensor{2,3})
   strain = zero(SymmetricTensor{2,3})
   backstress1 = zero(SymmetricTensor{2,3})
   backstress2 = zero(SymmetricTensor{2,3})
   plastic_strain = zero(SymmetricTensor{2,3})
   cumeq = 0.0
   R = 0.0
   state = AnotherState(stress, strain, backstress1,
                        backstress2, plastic_strain, cumeq, R)
   for i in 1:N
      dstress = randn(SymmetricTensor{2,3})
      dstrain = randn(SymmetricTensor{2,3})
      dbackstress1 = randn(SymmetricTensor{2,3})
      dbackstress2 = randn(SymmetricTensor{2,3})
      dplastic_strain = randn(SymmetricTensor{2,3})
      dcumeq = norm(dplastic_strain)
      dR = randn()
      dstate = AnotherState(dstress, dstrain, dbackstress1,
                           dbackstress2, dplastic_strain, dcumeq, dR)
      state = state + dstate
   end
   return state
end

println("Benchmark Chaboche State")
@btime bench_chaboche_style_state(N)


text:
"""
"""
struct ConditionalGaussian <: SciMLBase.AbstractLinearAlgorithm end
struct LinearGaussian <: SciMLBase.AbstractLinearAlgorithm end