text:
#------------gale_shapley.jl---------------------------------------------------#
#
# Purpose: To implement a solution to the stable Marriage problem
#
#------------------------------------------------------------------------------#

using DataStructures

function find_preference(n::Int)
    men_pref = fill(-1, n, n)
    women_pref = fill(-1, n, n)
    for i = 1:n
        men_pref[i,:] = shuffle(Vector(1:n))
        women_pref[i,:] = shuffle(Vector(1:n))
    end

    return (men_pref, women_pref)
end

function propose(wid::Int, mid::Int, woman::Array{Int,1}, pairs::Vector{Int})
    mid_0 = find(x -> x == wid, pairs)
    if(length(mid_0) != 0)
        rank_0 = woman[mid_0]
        rank_1 = woman[mid]
        if (rank_1[1] > rank_0[1])
            pairs[mid] = wid
        end
    else
        pairs[mid] = wid
        println(pairs)
    end
end

function gale_shapley(men::Array{Int, 2}, women::Array{Int, 2})
    n = Int(sqrt(length(men)))
    pairs = Vector{Int}(fill(-1.0, n))

    for i = 1:n
        for j = 1:n
            propose(men[j,i], j, women[j,:], pairs)
        end
    end

    return pairs
end

function main(n::Int)
    (men, women) = find_preference(n)
    pairs = gale_shapley(men, women)

    println(pairs)
end

main(10)

max_issues_repo_name: shiffman/simuleios

max_forks_repo_name: shiffman/simuleios


text:
mutable struct Task # It has all what we need to make predictions
	X
	y
	grp
	debiasmeasures

	function Task(X, y, grp, args...)
		cutoff = 0.05
		# groups with representation less than cutoff won't be appearing
		mask = ones(Bool, length(y))
		for group in levels(X[!, grp])
			if sum(X[!, grp].==group) < cutoff*length(y)
				mask[X[!, grp].==group] .= 0
			end
		end
		X, y = X[mask, :], y[mask]
		X[!, grp] = categorical(convert(Array, X[!, grp]))
		return new(X, y, grp, args...)
	end
end

function Task(X, y; grp, debiasmeasures)
	return Task(X, y, grp, debiasmeasures)
end

mutable struct FairnessProblem #It has all what we need to measure fairness using task
	task::Task
	refGrp # Reference Group for disparity calculation
	measures #Measures to evaluate apart from accuracy
	repls # replications
	nfolds # Number of folds for cross validation
	seed # Random Seed
	name # Name of the file or a unique name to indicate the experiment
end

function FairnessProblem(task::Task; refGrp=nothing,
								   measures=nothing, repls=10, nfolds=10,
								   seed=5, name=nothing)
   name = name==nothing ? randstring(['A':'Z'; '0':'9'], 12) : name
   measures = measures==nothing ? debiasmeasures : measures
   refGrp = refGrp==nothing ? StatsBase.mode(task.X[!, task.grp]) : refGrp
   task.X[!, task.grp] = convert(Array{Any}, task.X[!, task.grp])
   task.X[!, task.grp] = string.(task.X[!, task.grp])
   task.X[task.X[!, task.grp].!=refGrp, task.grp] .= "0"
   task.X[task.X[!, task.grp].==refGrp, task.grp] .= "1"
   task.X[!, task.grp] = categorical(convert(Array{String},
												task.X[!, task.grp]))
   transform!(task.X, task.grp => categorical, renamecols=false)
   return FairnessProblem(task, refGrp, measures, repls, nfolds, seed, name)
end

function FairnessProblem(X, y; grp, debiasmeasures, refGrp=nothing,
							measures=nothing, repls=10, nfolds=10,
							seed=5, name=nothing)
	return FairnessProblem(Task(X, y, grp, debiasmeasures),
							refGrp=refGrp, measures=measures, repls=repls,
							nfolds=nfolds, seed=seed, name=name)
end

ext: jl


text:
# Code for Video 1.6:

t1 = (108, 3.14,"Julia",'α') # (108,3.14,"Julia",'α')
typeof(t1) # Tuple{Int64,Float64,ASCIIString,Char}

a, b, c, d = t1 # (108,3.14,"Julia",'α')

a # 108
b # 3.14
c # "Julia"
d # 'α'

a, b = b, a # (3.14,108)
println("a is now $a, b is now $b")
# a is now 3.14, b is now 108

t1[3] # "Julia"

for e in t1
    print(e, " - ")
end
# 108 - 3.14 - Julia - α -

# t1[4] = 'β'
# ERROR: MethodError: `setindex!` has no method matching setindex!(::Tuple{Int64,F
# loat64,ASCIIString,Char}, ::Char, ::Int64)

s = Set(Int64[11, 14, 13, 7, 14, 11])
# Set([7,14,13,11])
s2 = Set(Any[11, 3.14, "Julia", 'α'])
# Set(Any[3.14,'α',"Julia",11])
intersect(s, s2) # Set([11])

diris = Dict("I. setosa" => 50,
             "I. versicolor" => 50,
             "I. virginica" => 50)
# Dict{ASCIIString,Int64} with 3 entries:
#   "I. virginica"  => 50
#   "I. versicolor" => 50
#   "I. setosa"     => 50

diriss = Dict( :setosa => 50,
               :versicolor => 50,
               :virginica => 50)
# Dict{Symbol,Int64} with 3 entries:
#   :versicolor => 50
#   :virginica  => 50
#   :setosa     => 50

diris["I. setosa"] # 50
diriss[:setosa] # 50

diriss[:setosa] = 53
diriss
# Dict{Symbol,Int64} with 3 entries:
#   :versicolor => 50
#   :virginica  => 50
#   :setosa     => 53

haskey(diris, "I. setosa") # true

keys(diris)
# Base.KeyIterator for a Dict{ASCIIString,Int64} with 3 entries. Keys:
#   "I. virginica"
#   "I. versicolor"
#   "I. setosa"

values(diris)
# Base.ValueIterator for a Dict{ASCIIString,Int64} with 3 entries. Values:
#   50
#   50
#   50

for k in keys(diris)
    print(k, " - ")
end
# I. virginica - I. versicolor - I. setosa -

for (k, v) in diris
    println("$k has $v measurements")
end
# I. virginica has 50 measurements
# I. versicolor has 50 measurements
# I. setosa has 50 measurements

max_stars_repo_path: Section 1/1_6.jl

max_stars_repo_name: lytemar/Julia-for-Data-Science-Video

max_issues_repo_path: Section 1/1_6.jl

max_issues_repo_name: lytemar/Julia-for-Data-Science-Video

max_forks_repo_path: Section 1/1_6.jl

max_forks_repo_name: lytemar/Julia-for-Data-Science-Video


text:
# Construct DiGraph
flow_graph = DiGraph(8)

# Load custom dataset
flow_edges = [
    (1,2,10),(1,3,5),(1,4,15),(2,3,4),(2,5,9),
    (2,6,15),(3,4,4),(3,6,8),(4,7,16),(5,6,15),
    (5,8,10),(6,7,15),(6,8,10),(7,3,6),(7,8,10)
]

capacity_matrix = zeros(Int ,nv(flow_graph) ,nv(flow_graph))

for e in flow_edges
    u,v,f = e
    add_edge!(flow_graph,u,v)
    capacity_matrix[u,v] = f
end

# Construct the residual graph
residual_graph = LightGraphs.residual(flow_graph)

# Test with default distances
@test LightGraphs.dinic_impl(residual_graph, 1, 8, LightGraphs.DefaultCapacity(residual_graph))[1] == 3

# Test with capacity matrix
@test LightGraphs.dinic_impl(residual_graph, 1, 8, capacity_matrix)[1] == 28

# Test on disconnected graphs
function test_blocking_flow(residual_graph, source, target, capacity_matrix, flow_matrix)
    #disconnect source
    h = copy(residual_graph)
    for dst in collect(neighbors(residual_graph, source))
        rem_edge!(h, source, dst)
    end
    @test LightGraphs.blocking_flow!(h, source, target, capacity_matrix, flow_matrix) == 0

    #disconnect target and add unreachable vertex
    h = copy(residual_graph)
    for src in collect(in_neighbors(residual_graph, target))
        rem_edge!(h, src, target)
    end
    @test LightGraphs.blocking_flow!(h, source, target, capacity_matrix, flow_matrix) == 0

    # unreachable vertex (covers the case where a vertex isn't reachable from the source)
    h = copy(residual_graph)
    add_vertex!(h)
    add_edge!(h, nv(residual_graph)+1, target)
    capacity_matrix_ = vcat(hcat(capacity_matrix, zeros(Int, nv(residual_graph))), zeros(Int, 1, nv(residual_graph)+1))
    flow_graph_  = vcat(hcat(flow_matrix, zeros(Int, nv(residual_graph))), zeros(Int, 1, nv(residual_graph)+1))

    @test LightGraphs.blocking_flow!(h, source, target, capacity_matrix_, flow_graph_ ) > 0

    #test with connected graph
    @test LightGraphs.blocking_flow!(residual_graph, source, target, capacity_matrix, flow_matrix) > 0
end

flow_matrix = zeros(Int, nv(residual_graph), nv(residual_graph))
test_blocking_flow(residual_graph, 1, 8, capacity_matrix, flow_matrix)


text:
get_backtrace(ex) = sprint(Base.showerror, ex, catch_backtrace())

function get_log_level(l)
    log_level_str = lowercase(l)

    (log_level_str == "debug") ? Logging.Debug :
    (log_level_str == "info")  ? Logging.Info  :
    (log_level_str == "warn")  ? Logging.Warn  : Logging.Error
end

function status_monitor(
    wait_and_close::Function,
    stop_file::AbstractString,
    running_check::Function,
)
    while running_check()
        sleep(5)
        flush(stdout); flush(stderr);
        # stop server if stop is requested
        if isfile(stop_file)
            @warn "Server stop requested."
            flush(stdout); flush(stderr)

            # stop accepting new requests
            wait_and_close()

            rm(stop_file; force=true)
        end
    end
end

function status_monitor(stop_file::AbstractString, zsock::MessageSocket)
    status_monitor(() -> close(zsock.sock), stop_file, () -> isopen(zsock.sock))
end

function status_monitor(
    stop_file::AbstractString,
    event_queue::Channel,
    httpsock::Ref{Sockets.TCPServer},
)
    status_monitor(stop_file, () -> isopen(event_queue)) do
        # wait for queued requests to be processed and close queue
        while isready(event_queue)
            @info "Waiting for queued jobs to finish"
            yield()
        end
        close(httpsock[])
        close(event_queue)
    end
end

function recover(
    name::AbstractString, keep_running::Function,
    do_action::Function, handle_exception::Function;
    backoff=0, backoffmax=120, backoffincrement=1,
)
    while keep_running()
        try
            do_action()
            backoff = 0
        catch ex
            exception_action = handle_exception(ex)
            if exception_action === :exit
                @warn("Stopping", name)
                return
            else # exception_action == :continue
                @error("Recovering from unknown exception", name, backoff, exception=get_backtrace(ex))
                sleep(backoff)
                backoff = min(backoffmax, backoff+backoffincrement)
            end
        end
    end
end


text:
using ITensors
using ITensorMPOCompression
using Test
using Revise

#ITensors.ITensors.disable_debug_checks()
ITensors.ITensors.enable_debug_checks()
#function run_tests()
@testset "ITensorMPOCompression.jl" begin
    include("qx_unittests.jl")
    include("blocking.jl")
    include("hamiltonians.jl")
    include("orthogonalize.jl")
    include("truncate.jl")
end
Nothing #suppress messy dump from Test

#@time run_tests()


text:
export iauC2i06a
"""
Form the celestial-to-intermediate matrix for a given date using the
IAU 2006 precession and IAU 2000A nutation models.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  support function.

Given:
   date1,date2 double       TT as a 2-part Julian Date (Note 1)

Returned:
   rc2i        double[3][3] celestial-to-intermediate matrix (Note 2)

Notes:

   1. The TT date date1+date2 is a Julian Date, apportioned in any
   convenient way between the two arguments.  For example,
   JD(TT)=2450123.7 could be expressed in any of these ways,
   among others:

         date1          date2

      2450123.7           0.0       (JD method)
      2451545.0       -1421.3       (J2000 method)
      2400000.5       50123.2       (MJD method)
      2450123.5           0.2       (date & time method)

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable.  The J2000 method is best matched to the way
   the argument is handled internally and will deliver the
   optimum resolution.  The MJD method and the date & time methods
   are both good compromises between resolution and convenience.

   2. The matrix rc2i is the first stage in the transformation from
   celestial to terrestrial coordinates:

      [TRS]  =  RPOM * R_3(ERA) * rc2i * [CRS]

            =  RC2T * [CRS]

   where [CRS] is a vector in the Geocentric Celestial Reference
   System and [TRS] is a vector in the International Terrestrial
   Reference System (see IERS Conventions 2003), ERA is the Earth
   Rotation Angle and RPOM is the polar motion matrix.

Called:
   iauPnm06a    classical NPB matrix, IAU 2006/2000A
   iauBpn2xy    extract CIP X,Y coordinates from NPB matrix
   iauS06       the CIO locator s, given X,Y, IAU 2006
   iauC2ixys    celestial-to-intermediate matrix, given X,Y and s

References:

   McCarthy, D. D., Petit, G. (eds.), 2004, IERS Conventions (2003),
   IERS Technical Note No. 32, BKG

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauC2i06a(double date1, double date2, double rc2i[3][3])

function iauC2i06a(date1::Real, date2::Real)
   # Allocate return value
   rc2i = zeros(Float64, 3, 3)

   ccall((:iauC2i06a, libsofa_c), Cvoid,
         (Cdouble, Cdouble,
         Ptr{Cdouble}),
         convert(Float64, date1), convert(Float64, date2),
         rc2i)

   return SMatrix{3,3}(rc2i')
end


text:
#
# Date created: 2022-03-09
# Author: aradclif
#
#
############################################################################################
# Building on lvreduce.jl

# function smul_inplacebody(N::Int)
#     body = Expr(:block)
#     a = Expr(:ref, :A, ntuple(d -> Symbol(:i_, d), N)...)
#     e = Expr(:(=), a, Expr(:call, :*, a, :x))
#     push!(body.args, e)
#     body
# end

function transform_inplacebody(f, N::Int)
    body = Expr(:block)
    a = Expr(:ref, :A, ntuple(d -> Symbol(:i_, d), N)...)
    e = Expr(:(=), a, Expr(:call, Symbol(f), a))
    push!(body.args, e)
    body
end

# function smul_inplace_quote(N::Int)
#     ls = loopgen(N)
#     body = smul_inplacebody(N)
#     push!(ls.args, body)
#     return quote
#         @turbo $ls
#         return A
#     end
# end
# @generated function smul!(A::AbstractArray{T, N}, x::T) where {T, N}
#     smul_inplace_quote(N)
# end

function smul!(A::AbstractArray{T, N}, x::T) where {T, N}
    @turbo for i ∈ eachindex(A)
        A[i] *= x
    end
    A
end

function smul!(B::AbstractArray{Tₒ, N}, A::AbstractArray{T, N}, x::Tₒ) where {Tₒ, T, N}
    @turbo for i ∈ eachindex(A)
        B[i] = A[i] * x
    end
    B
end
smul(A::AbstractArray{T, N}, x::Tₒ) where {Tₒ, T, N} = smul!(similar(A, Tₒ), A, x)

function lvmean(A::AbstractArray{T, N}, dims::NTuple{M, Int}) where {T, N, M}
    if ntuple(identity, Val(N)) ⊆ dims
        return hvncat(ntuple(i -> 1, Val(N)), true, lvmean1(A))
    else
        B = lvsum(A, dims=dims)
        Dᴬ = size(A)
        denom = 1
        for d ∈ eachindex(Dᴬ)
            denom = d ∈ dims ? denom * Dᴬ[d] : denom
        end
        x = inv(denom)
        # smul!(B, x)
        # B
        return eltype(B) <: AbstractFloat ? smul!(B, x) : smul(B, x)
    end
end
lvmean(A::AbstractArray{T, N}, dims::Int) where {T, N} = lvmean(A, (dims,))
lvmean(A::AbstractArray{T, N}; dims=:) where {T, N} = lvmean(A, dims)
lvmean(A::AbstractArray{T, N}) where {T, N} = lvmean1(A)
lvmean(A::AbstractArray{T, N}, ::Colon) where {T, N} = lvmean1(A)

function lvmean1(A::AbstractArray{T, N}) where {T, N}
    s = zero(Base.promote_op(+, T, Int))
    @turbo for i ∈ eachindex(A)
        s += A[i]
    end
    s / length(A)
end

################ the mapreduce version analogous to Julia base Statistics
function lvmean(f, A::AbstractArray{T, N}, dims::NTuple{M, Int}) where {T, N, M}
    if ntuple(identity, Val(N)) ⊆ dims
        return hvncat(ntuple(i -> 1, Val(N)), true, lvmean1(f, A))
    else
        B = lvsum(f, A, dims=dims)
        Dᴬ = size(A)
        denom = 1
        for d ∈ eachindex(Dᴬ)
            denom = d ∈ dims ? denom * Dᴬ[d] : denom
        end
        x = inv(denom)
        # smul!(B, x)
        # B
        return eltype(B) <: AbstractFloat ? smul!(B, x) : smul(B, x)
    end
end
lvmean(f, A::AbstractArray{T, N}, dims::Int) where {T, N} = lvmean(f, A, (dims,))
lvmean(f, A::AbstractArray{T, N}; dims=:) where {T, N} = lvmean(f, A, dims)
lvmean(f, A::AbstractArray{T, N}) where {T, N} = lvmean1(f, A)
lvmean(f, A::AbstractArray{T, N}, ::Colon) where {T, N} = lvmean1(f, A)

@generated function lvmean1(f::F, A::AbstractArray{T, N}) where {F, T, N}
    f = F.instance
    Tₒ = Base.promote_op(+, Base.promote_op(f, T), Int)
    quote
        s = zero($Tₒ)
        @turbo for i ∈ eachindex(A)
            s += $f(A[i])
        end
        s / length(A)
    end
end

################ threaded version

function tsmul!(A::AbstractArray{T, N}, x::T) where {T, N}
    @tturbo for i ∈ eachindex(A)
        A[i] *= x
    end
    A
end

function tsmul!(B::AbstractArray{Tₒ, N}, A::AbstractArray{T, N}, x::Tₒ) where {Tₒ, T, N}
    @tturbo for i ∈ eachindex(A)
        B[i] = A[i] * x
    end
    B
end
tsmul(A::AbstractArray{T, N}, x::Tₒ) where {Tₒ, T, N} = tsmul!(similar(A, Tₒ), A, x)

function lvtmean(A::AbstractArray{T, N}, dims::NTuple{M, Int}) where {T, N, M}
    if ntuple(identity, Val(N)) ⊆ dims
        return hvncat(ntuple(i -> 1, Val(N)), true, lvtmean1(A))
    else
        B = lvtsum(A, dims=dims)
        Dᴬ = size(A)
        denom = 1
        for d ∈ eachindex(Dᴬ)
            denom = d ∈ dims ? denom * Dᴬ[d] : denom
        end
        x = inv(denom)
        return eltype(B) <: AbstractFloat ? tsmul!(B, x) : tsmul(B, x)
    end
end
lvtmean(A::AbstractArray{T, N}, dims::Int) where {T, N} = lvtmean(A, (dims,))
lvtmean(A::AbstractArray{T, N}; dims=:) where {T, N} = lvtmean(A, dims)
lvtmean(A::AbstractArray{T, N}) where {T, N} = lvtmean1(A)
lvtmean(A::AbstractArray{T, N}, ::Colon) where {T, N} = lvtmean1(A)

function lvtmean1(A::AbstractArray{T, N}) where {T, N}
    s = zero(Base.promote_op(+, T, Int))
    @tturbo for i ∈ eachindex(A)
        s += A[i]
    end
    s / length(A)
end

################
function lvtmean(f, A::AbstractArray{T, N}, dims::NTuple{M, Int}) where {T, N, M}
    if ntuple(identity, Val(N)) ⊆ dims
        return hvncat(ntuple(i -> 1, Val(N)), true, lvtmean1(f, A))
    else
        B = lvtsum(f, A, dims=dims)
        Dᴬ = size(A)
        denom = 1
        for d ∈ eachindex(Dᴬ)
            denom = d ∈ dims ? denom * Dᴬ[d] : denom
        end
        x = inv(denom)
        return eltype(B) <: AbstractFloat ? tsmul!(B, x) : tsmul(B, x)
    end
end
lvtmean(f, A::AbstractArray{T, N}, dims::Int) where {T, N} = lvtmean(f, A, (dims,))
lvtmean(f, A::AbstractArray{T, N}; dims=:) where {T, N} = lvtmean(f, A, dims)
lvtmean(f, A::AbstractArray{T, N}) where {T, N} = lvtmean1(f, A)
lvtmean(f, A::AbstractArray{T, N}, ::Colon) where {T, N} = lvtmean1(f, A)

@generated function lvtmean1(f::F, A::AbstractArray{T, N}) where {F, T, N}
    f = F.instance
    Tₒ = Base.promote_op(+, Base.promote_op(f, T), Int)
    quote
        s = zero($Tₒ)
        @tturbo for i ∈ eachindex(A)
            s += $f(A[i])
        end
        s / length(A)
    end
end


text:
"""
bicgstab_test.jl

Tests the linear BiCGSTAB code, kl_bicgstab.

This is for CI only. Nothing to see here. Move along.
"""

function bicgstab_test()
    pass3 = test3x3()
    passint = test_integop()
    passr1 = testR1()
    passbicgs = pass3 && passint && passr1
    return passbicgs
end

function test3x3()
    A = [0.001 0 0; 0 0.0011 0; 0 0 1.e4];
    V = zeros(3); b = [1.0; 1.0; 1.0]; x0 = zeros(3);
    eta = 1.e-10
    gout = kl_bicgstab(x0, b, atv, V, 1.e-10; pdata = A)
    pass = (length(gout.reshist)== 5) &&
            (norm(A*gout.sol - b,Inf) < 1.e-12) && gout.idid && (gout.lits==4)
#    return (gout=gout, pass=pass)
pass || println("3x3 test fails")
    return pass
end

function testR1()
    A = Float64.([1 2 3 4 5])
    E = A' * A
    A = I + E
    b = ones(5)
    x0 = zeros(5)
    V = zeros(5)
    gout = kl_bicgstab(x0, b, atv, V, 1.e-7; pdata = A)
    pass = (length(gout.reshist)== 3) && (norm(A*gout.sol - b,Inf) < 1.e-12)
#    return (gout=gout, pass=pass)
pass || println("R1 test fails")
    return pass
end

function test_integop(n=100)
    pdata = integopinit(n)
    f = pdata.f
    ue = pdata.xe
    u0 = zeros(size(f)); V = zeros(size(f))
    gout = kl_bicgstab(u0, f, integop, V, 1.e-10; pdata = pdata)
    realres = (I - pdata.K)*gout.sol - f
pass = ( (norm(realres,Inf) < 1.e-12) && (length(gout.reshist)==4))
#    return (gout=gout, pass=pass)
pass || println("integop test fails")
     return pass
end

function atv(x, A)
    return A * x
end

function integop(u, pdata)
    K = pdata.K
#    f = pdata.f
    return u - K * u
end

function integopinit(n)
    h = 1 / n
    X = collect(0.5*h:h:1.0-0.5*h)
    K = [ker(x, y) for x in X, y in X]
    #    K = zeros(n, n)
    #    for j = 1:n
    #        for i = 1:n
    #            K[i, j] = ker(x[i], x[j])
    #        end
    #    end
    K .*= h
    #    sol = exp.(x) .* log.(2.0 * x .+ 1.0)
    #    sol = usol.(X)
    sol = [usol(x) for x in X]
    f = sol - K * sol
    pdata = (K = K, xe = sol, f = f)
    return pdata
end

function usol(x)
    return exp.(x) .* log.(2.0 * x .+ 1.0)
end

function ker(x, y)
    ker = 0.1 * sin(x + exp(y))
end


text:
# downsample
SUITE["downsample"] = BenchmarkGroup()
let T=Float32, sz=(100, 200, 40)
    x = ones(T, sz)
    SUITE["downsample"]["perm"] = @benchmarkable MIRT.downsample3_perm($x, (2, 2, 2))
    SUITE["downsample"]["loop"] = @benchmarkable MIRT.downsample3_loop($x, [2, 2, 2])
end


text:
#!/usr/bin/env julia

using DataKnots
using NarrativeTest

subs = NarrativeTest.common_subs()

# Ignore the difference in the output of `print(Int)` between 32-bit and 64-bit platforms.
push!(subs, r"Int64" => s"Int(32|64)")

# Normalize printing of `Vector{Bool}`.
Base.show(io::IO, b::Bool) = print(io, get(io, :typeinfo, Any) === Bool ? (b ? "1" : "0") : (b ? "true" : "false"))

# Normalize `LoadError` output under 1.2.
if VERSION >= v"1.2.0-DEV"
    function Base.showerror(io::IO, ex::LoadError, bt; backtrace=true)
        print(io, "LoadError: ")
        showerror(io, ex.error, bt, backtrace=backtrace)
        print(io, "\nin expression starting at $(ex.file):$(ex.line)")
    end
end

# Normalize printing of Enum values.
if VERSION < v"1.2.0-DEV"
    Base.show(io::IO, c::DataKnots.Cardinality) =
        print(io, Symbol(c))
end

# Remove type prefix for Vector{Symbol}.
if VERSION < v"1.4.0-DEV"
    Base.typeinfo_prefix(::IO, ::Vector{Symbol}) = ""
end

# Fix output of character ranges.
if VERSION < v"1.5.0-DEV"
    Base.step(r::StepRangeLen) = r.step
    Base.step(r::StepRangeLen{T}) where {T<:AbstractFloat} = T(r.step)
end

# Normalize printing of vector types.
if VERSION < v"1.6.0-DEV"
    Base.show_datatype(io::IO, x::Type{Vector{T}}) where {T} = print(io, "Vector{$T}")
end

# Normalize printing of type parameters.
if VERSION < v"1.6.0-DEV"
    function Base.show_datatype(io::IO, x::DataType)
        istuple = x.name === Tuple.name
        if (!isempty(x.parameters) || istuple) && x !== Tuple
            n = length(x.parameters)::Int
            if istuple && n > 3 && all(i -> (x.parameters[1] === i), x.parameters)
                print(io, "NTuple{", n, ", ", x.parameters[1], "}")
            else
                Base.show_type_name(io, x.name)
                print(io, '{')
                for (i, p) in enumerate(x.parameters)
                    show(io, p)
                    i < n && print(io, ", ")
                end
                print(io, '}')
            end
        else
            Base.show_type_name(io, x.name)
        end
    end
end

# Set the width to 72 so that MD->PDF via pandoc fits the page.
ENV["COLUMNS"] = "72"

package_path(x) = relpath(joinpath(dirname(abspath(PROGRAM_FILE)), "..", x))
default = package_path.(["doc/src", "README.md"])
runtests(; default=default, subs=subs)


text:
"""
    localtest.jl -- Feel free to modify anything you like in this file!

You can call this file as a script at the command line, or `include` it from within a julia session.

Example command-line usage:

    julia --color=yes localtest.jl                  # Evaluate all simple functions the default number of times (500):
    julia --color=yes localtest.jl 2000             # Evaluate all simple functions a custom number of times (2000 in this case):
    julia --color=yes localtest.jl simple1          # Evaluate only simple1 the default number of times:
    julia --color=yes localtest.jl simple3 2000     # Evaluate only simple3 a custom number of times:
"""

# Include the other relevant files:
include(joinpath("project1_jl", "helpers.jl"))
include(joinpath("project1_jl", "project1.jl"))

if length(ARGS) == 0
    # Default
    K = 500
    probnames = sort(collect(keys(PROBS)))
elseif length(ARGS) == 1
    # One input argument can be either evaluation count or function name, so
    # have to test for both. If it's not an integer it's got to be a function name
    K = tryparse(Int, ARGS[1])
    if K === nothing
        K = 500
        probnames = [ARGS[1]]
    end
elseif length(ARGS) == 2
    # Two inputs can come in at any order, so try both
    K = tryparse(Int, ARGS[1])
    if K === nothing
        K = tryparse(Int, ARGS[2])
        if K === nothing
            throw(ArgumentError("Can't detect which input is intended to be the evaluation count. Make sure it's an integer"))
        end
        probnames = [ARGS[1]]
    else
        probnames = [ARGS[2]]
    end
else
    throw(ArgumentError("Too many command-line inputs to localtest.jl"))
end

# This is the same random search implemented in the autograder.
# We take normally distributed random steps away from x0 and
# pick the best one.
function autograder_random_search(f, x0, n)
    xs = [x0 + randn(length(x0)) for i in 1:n]
    i = argmin(f.(xs))
    return xs[i]
end

function K_random_searches(f, x0, n, K, seed = 42)
    scores = zeros(K)
    for i in 1:K
        Random.seed!(seed+i)
        x_star = autograder_random_search(f, x0(), n)
        scores[i] = f(x_star)

        seed += 1
    end
    return scores
end

function compare_optimize_to_random(probname, K)
    f, g, x0, n = PROBS[probname]

    # Use a consistent random number stream for the comparison
    seed = 42

    random_scores = K_random_searches(f, x0, n, 500, seed)
    scores, n_evals = main(probname, K, optimize, seed)

    if maximum(n_evals) > n
        @warn "Problem $probname saw $(maximum(n_evals)) out of $n allowed calls, which will be an automatic 0 on Gradescope"
    end

    # lower is better
    better = scores .< random_scores

    return mean(better)
end

printstyled("\nTesting $K times\n\n", bold = true)
for nm in probnames
    try
        better_percent = 100 * compare_optimize_to_random(nm, K)

        if better_percent > 55
            printstyled("Pass: optimize does better than random search $(better_percent)% of the time on $(nm)!\n", color = :green)
        else
            printstyled("Fail: optimize only beat random search $(better_percent)% of the time on $nm.\n", color = :red)
        end

    catch e
        println("ERROR IN PROBLEM $nm")
        showerror(stdout, e, catch_backtrace())
    end
end


text:
for T in [
    :Chain, :Parallel, :SkipConnection, :Recur, :Maxout  # container types
  ]
  @eval function Base.show(io::IO, m::MIME"text/plain", x::$T)
    if get(io, :typeinfo, nothing) === nothing  # e.g. top level in REPL
      _big_show(io, x)
    elseif !get(io, :compact, false)  # e.g. printed inside a Vector, but not a Matrix
      _layer_show(io, x)
    else
      show(io, x)
    end
  end
end

function _big_show(io::IO, obj, indent::Int=0, name=nothing)
  pre, post = obj isa Chain{<:AbstractVector} ? ("([", "])") : ("(", ")")
  children = _show_children(obj)
  if all(_show_leaflike, children)
    _layer_show(io, obj, indent, name)
  else
    println(io, " "^indent, isnothing(name) ? "" : "$name = ", nameof(typeof(obj)), pre)
    if obj isa Chain{<:NamedTuple} && children == getfield(obj, :layers)
      # then we insert names -- can this be done more generically?
      for k in Base.keys(obj)
        _big_show(io, obj[k], indent+2, k)
      end
    elseif obj isa Parallel{<:Any, <:NamedTuple}
      _big_show(io, obj.connection, indent+2)
      for k in Base.keys(obj)
        _big_show(io, obj[k], indent+2, k)
      end
    else
      for c in children
        _big_show(io, c, indent+2)
      end
    end
    if indent == 0  # i.e. this is the outermost container
      print(io, rpad(post, 2))
      _big_finale(io, obj)
    else
      println(io, " "^indent, post, ",")
    end
  end
end

_show_leaflike(x) = isleaf(x)  # mostly follow Functors, except for:
_show_leaflike(::Tuple{Vararg{<:Number}}) = true         # e.g. stride of Conv
_show_leaflike(::Tuple{Vararg{<:AbstractArray}}) = true  # e.g. parameters of LSTMcell
_show_leaflike(::Diagonal) = true                        # appears inside LayerNorm

_show_children(x) = trainable(x)  # except for layers which hide their Tuple:
_show_children(c::Chain) = c.layers
_show_children(m::Maxout) = m.layers
_show_children(p::Parallel) = (p.connection, p.layers...)

for T in [
    :Conv, :ConvTranspose, :CrossCor, :DepthwiseConv, :Dense,
    :BatchNorm, :LayerNorm, :InstanceNorm, :GroupNorm,
  ]
  @eval function Base.show(io::IO, m::MIME"text/plain", x::$T)
    if !get(io, :compact, false)
      _layer_show(io, x)
    else
      show(io, x)
    end
  end
end

function _layer_show(io::IO, layer, indent::Int=0, name=nothing)
  _str = isnothing(name) ? "" : "$name = "
  str = _str * sprint(show, layer, context=io)
  print(io, " "^indent, str, indent==0 ? "" : ",")
  if !isempty(params(layer))
    print(io, " "^max(2, (indent==0 ? 20 : 39) - indent - length(str)))
    printstyled(io, "# ", underscorise(sum(length, params(layer))), " parameters"; color=:light_black)
    nonparam = _childarray_sum(length, layer) - sum(length, params(layer))
    if nonparam > 0
      printstyled(io, ", plus ", underscorise(nonparam), indent==0 ? " non-trainable" : ""; color=:light_black)
    end
    _nan_show(io, params(layer))
  end
  indent==0 || println(io)
end

function _big_finale(io::IO, m)
  ps = params(m)
  if length(ps) > 2
    pars = underscorise(sum(length, ps))
    bytes = Base.format_bytes(Base.summarysize(m))
    noncnt = _childarray_sum(_->1, m) - length(ps)
    if noncnt > 0
      nonparam = underscorise(_childarray_sum(length, m) - sum(length, ps))
      printstyled(io, " "^08, "# Total: ", length(ps), " trainable arrays, "; color=:light_black)
      println(io, pars, " parameters,")
      printstyled(io, " "^10, "# plus ", noncnt, " non-trainable, ", nonparam, " parameters, summarysize "; color=:light_black)
      print(io, bytes, ".")
    else
      printstyled(io, " "^18, "# Total: ", length(ps), " arrays, "; color=:light_black)
      print(io, pars, " parameters, ", bytes, ".")
    end
  end
end

_childarray_sum(f, x::AbstractArray{<:Number}) = f(x)
_childarray_sum(f, x) = isleaf(x) ? 0 : sum(y -> _childarray_sum(f, y), Functors.children(x))

# utility functions

underscorise(n::Integer) =
  join(reverse(join.(reverse.(Iterators.partition(digits(n), 3)))), '_')

function _nan_show(io::IO, x)
  if !isempty(x) && _all(iszero, x)
    printstyled(io, "  (all zero)", color=:cyan)
  elseif _any(isnan, x)
    printstyled(io, "  (some NaN)", color=:red)
  elseif _any(isinf, x)
    printstyled(io, "  (some Inf)", color=:red)
  end
end

_any(f, xs::AbstractArray{<:Number}) = any(f, xs)
# _any(f, xs::Union{Tuple,NamedTuple,Zygote.Params}) = any(x -> _any(f, x), xs)
_any(f, xs) = any(x -> _any(f, x), xs)
_any(f, x::Number) = f(x)
# _any(f, x) = false

_all(f, xs) = !_any(!f, xs)


text:
export rho,
    Cp,
    visc,
    therm_cond,
    therm_expand,
    compress,
    std_props

## Std props

function std_props(fluid::Fluid)
    ρ = fluid.ρ
end

## rho

function rho(P, T, fluid::Fluid)
    @error "Property not yet implemented for the fluid. - ρ"
    throw(NotImplementedError())
    return missing
end

function rho(P, T, fluid::Water)
    a = 8.79278807E+02
    b = 1.11820207E+00
    c = -2.40938185E-03
    d = 1.29027705E-06
    e = -5.46746114E-09
    f = 8.06188264E-12

    return a + b * T + c * T^2 + d * P + e * P * T + f * P * T^2
end

function rho(P, T, fluid::BrineNaCl)

    x = fluid.x

    a, b, c, d, g, h, m, n, p, q, r, s = 8.23073018e+02, 1.38381011e+00, -2.85110020e-03, -9.11569454e-18, 4.04861229e-08, 2.38570427e-15, -2.65736884e-09, 9.52410475e-12, 1.27087646e+03, -7.83633987e-01, 6.26293103e-09, -2.11977575e-03

    d1 = a + b * T + c * T^2
    d2 = g * P + h * P^2
    d3 = m * P * T + n * P * T^2 + d * P^2 * T
    d4 = p * x + q * x * T + r * x * P * T + s * x * T^2

    return d1 + d2 + d3 + d4
end

## Cp

function Cp(P, T, fluid::Fluid)
    @error "Property not yet implemented for the fluid. - Cₚ"
    throw(NotImplementedError())
    return missing
end

function Cp(P, T, fluid::Water)
    b = 3.70641927E+03
    c = 7.13116317E-01
    d = 1.05153132E-03
    e = 8.59719978E-07
    f = -4.37409377E-09
    g = -4.26280988E-15
    h = 1.14112829E-17

    return b + 2 * c * T + e * P + 2 * f * P * T + g * P * 2 + 2 * h * P^2 * T
end

function Cp(P, T, fluid::BrineNaCl)

    x = fluid.x

    b, c, d, e, f, g, i, j, k, l, m = (2.48218846e+03, 3.80939153e-06, -3.45541918e-15, -2.57191794e+03, 3.85762099e+03, 1.93733545e-07, 2.27236796e+00, -9.07624069e-09, 1.27273598e-17, 1.38313250e+00, -1.96832796e+00)

    h_1 = b + c * P + d * P^2 + e * x + f * x^2 + g * x * P
    h_2 = 2 * (i + j * P + k * P^2 + l * x + m * x^2) * T
    return h_1 + h_2

end

## k

function therm_cond(P, T, fluid::Fluid)
    @error "Property not yet implemented for the fluid. - k"
    throw(NotImplementedError())
    return missing
end

function therm_cond(P, T, fluid::Water)
    a = -3.98200943E-01
    b = 5.32706765E-03
    c = -6.55323715E-06
    d = 1.19857114E-09
    e = -5.13478793E-12
    f = 9.20668097E-15
    g = -2.77585855E-24

    return a + b * T + c * T^2 + d * P + e * P * T + f * P * T^2 + g * P^2 * T^2
end

function therm_cond(P, T, fluid::BrineNaCl)

    x = fluid.x

    a, b, c, d, g, h, m, n, p, q, r, s, t = (-9.60221679e-01, 1.06170815e-02, -2.25969876e-05, 1.57190072e-08, -8.34381967e-10, -2.63526539e-19, 4.75581304e-12, -3.22756205e-15, 4.47716476e-01, -4.57135183e-03, 1.55078050e-09, -2.82355304e-12, 8.68547458e-06)

    k_1 = a + b * T + c * T^2 + d * T^3
    k_2 = g * P + h * P^2
    k_3 = m * P * T + n * P * T^2
    k_4 = p * x + q * x * T + r * x * P + s * x * P * T + t * x * T^2

    return k_1 + k_2 + k_3 + k_4
end

# mu

function visc(P, T, fluid::Fluid)
    @error "Property not yet implemented for the fluid. - μ"
    throw(NotImplementedError())
    return missing
end

function visc(P, T, fluid::Water)
    a = 2.65387609E+00
    b = 2.37561266E+02
    c = -2.37488893E+02
    d = -6.43059846E+08
    e = 6.43115832E+08
    f = 7.34303718E-07
    g = -9.76896404E-01

    return a * (exp((b - T) / (c + T) * ((P + d) / (P - e)))) + f * T + g
end

function visc(P, T, fluid::BrineNaCl)

    x = fluid.x

    a, b, c, d, g, h, m, p, q, r = (-9.29821669e+00, -5.62181704e+02, -2.60785572e-07, 2.24278114e+03, -1.75402843e-16, 3.70625444e+05, -4.33694111e+05, 4.47607787e-12, -1.44053809e-20, -1.33905655e-05)

    v_1 = a + (b + c * P + d * x + g * P^2) / T
    v_2 = (h + m * x) / T^2
    v_3 = p * P * T + q * P^2 * T + r * T^2 * x^2
    return exp(v_1 + v_2 + v_3)
end

# therm expand

function therm_expand(P, T, fluid::Fluid)
    return -gradient(x -> rho(P, x, fluid), T)[1] / rho(P, T, fluid)
end

## compressibility

function compress(P, T, fluid::Fluid)
    return gradient(x -> rho(x, T, fluid), P)[1] / rho(P, T, fluid)
end

# Solid

## density

function rho(P, T, fluid::Solid)
    @error "Property not yet implemented for solid - ρ"
    throw(NotImplementedError())
    return missing
end

function rho(P, T, solid::Barite)
    return solid.ρ
end

function rho(P, T, solid::CaCO3)
    return solid.ρ
end

function rho(P,T,solid::Hematite)
    return solid.ρ
end

## Cp

function Cp(P, T, solid::Solid)
    @error "Property not yet implemented for solid - Cₚ"
    throw(NotImplementedError())
    return missing
end

function Cp(P, T, solid::Barite)
    return 4600.0
end

function Cp(P, T, solid::CaCO3)
    return solid.cₚ
end

function Cp(P,T,solid::Hematite)
    A=93.43834
    B=108.3577
    C=-50.86447
    D=25.58683
    E=-1.611330
    t=T/1000.0
    CpM=A+B*t+C*t^2+D*t^3+E/(t^2)
    return CpM/159.688*1000.0
end

## therm cond

function therm_cond(P, T, solid::Solid)
    @error "Property not yet implemented for solid - k"
    throw(NotImplementedError())
    return missing
end

function therm_cond(P, T, solid::Barite)
    Tc = T - 273.15
    return -9.0E-5 * Tc + 0.0974
end

function therm_cond(P, T, solid::CaCO3)
    return solid.k
end

function therm_cond(P, T, solid::Hematite)
    return (0.0839-6.63E-5*T)*100.0
end


text:
using OpenCL

const sum_kernel_src = "
   __kernel void sum(__global const float *a,
                     __global const float *b,
                     __global float *c)
    {
      int gid = get_global_id(0);
      c[gid] = a[gid] + b[gid];
    }
"
a = rand(Float32, 50_000)
b = rand(Float32, 50_000)

device, ctx, queue = cl.create_compute_context()

# create opencl buffer objects
# copies to the device initiated when the kernel function is called
a_buff = cl.Buffer(Float32, ctx, (:r, :copy), hostbuf=a)
b_buff = cl.Buffer(Float32, ctx, (:r, :copy), hostbuf=b)
c_buff = cl.Buffer(Float32, ctx, :w, length(a))

# build the program and construct a kernel object
p = cl.Program(ctx, source=sum_kernel_src) |> cl.build!
sum_kernel = cl.Kernel(p, "sum")

# call the kernel object with global size set to the size our arrays
sum_kernel[queue, size(a)](a_buff, b_buff, c_buff)

# perform a blocking read of the result from the device
r = cl.read(queue, c_buff)

# check to see if our result is what we expect!
if isapprox(norm(r - (a+b)), zero(Float32))
    info("Success!")
else
    error("Norm should be 0.0f")
end


text:
export UpdateMassey

"""
    UpdateMassey

 Massey update rule, e.g., see "Whos's #1", Langville and Meyer, p.9

## Parameters (none)
```
"""
struct UpdateMassey <: UpdateRule
# this is a batch calculation that takes no account of past ratings, and
# has no parameters
end

function update_info( rule::UpdateMassey )
    info = Dict(
                :name => "Massey",
                :reference => "\"Whos's #1\", Langville and Meyer, p.9",
                :computation => "simultaneous",
                :state_model => "none",
                :input => "score",
                :output => "deterministic",
                :model => "single",
                :ties => true,
                :factors => false,
                :parameters => [],
                :record => true
                )
    return info
end
update_info( ::Type{UpdateMassey} ) =  update_info( UpdateMassey() )

function update_ratings( rule::UpdateMassey,
                         input_ratings::RatingsList,
                         input_competitions::DataFrame)
    n = size(input_competitions,1)
    m = length( input_ratings.players )
    I = player_indexes( input_ratings.players )

    # construct Massey matrices and vectors
    M = zeros(Int,m,m)
    p = zeros(Int, m)
    d = input_competitions # just an abbreviation
    point_diff = d[:,ScoreA] - d[:,ScoreB]
    for i=1:n
        p[ I[d[i,PlayerA]] ] +=  point_diff[i]
        p[ I[d[i,PlayerB]] ] -=  point_diff[i]
        M[ I[d[i,PlayerA]], I[d[i,PlayerB]] ] -= 1
        M[ I[d[i,PlayerB]], I[d[i,PlayerA]] ] -= 1
        M[ I[d[i,PlayerA]], I[d[i,PlayerA]] ] += 1
        M[ I[d[i,PlayerB]], I[d[i,PlayerB]] ] += 1
    end
    M̅  = copy(M) # "M\u0305" to get Mbar
    M̅[end,:] = ones(size(M[end,:]))
    p̅  = copy(p)
    p̅[end] = 0

    # solve Massey's equation
    r = M̅  \  p̅
    ratings = Dict{String, Float64}()
    for player in input_ratings.players
        ratings[player] = r[ I[player] ]
    end

    # output ratings list
    output_ratings = RatingsList(input_ratings.players, ratings )
    return output_ratings
end


text:
using BenchmarkTools
using ImageDistances
using Distances
using TestImages
using ImageBinarization

const SUITE = BenchmarkGroup()

on_CI = haskey(ENV, "GITHUB_ACTIONS")

function create_distances()
    dists = [
        SqEuclidean(),
        Euclidean(),
        Cityblock(),
        TotalVariation(),
        Minkowski(2),
        Hamming(),
        SumAbsoluteDifference(),
        SumSquaredDifference(),
        MeanAbsoluteError(),
        MeanSquaredError(),
        RootMeanSquaredError(),
        ZNCC(),
        Hausdorff(),
        ModifiedHausdorff(),
        CIEDE2000(),
    ]
    return dists
end

#############
# Distances #
############

SUITE["distances"] = BenchmarkGroup()

function add_distances_benchmarks!(SUITE)

    imgA = testimage("cameraman")
    imgB = testimage("livingroom")

    alg = Otsu()
    imgA = binarize(imgA, alg)
    imgB = binarize(imgA, alg)

    imgC = Float32.(imgA)
    imgD = Float32.(imgB)

    dists = create_distances()

    for dist in (dists)
        Tdist = typeof(dist)
        SUITE["distances"][Tdist] = BenchmarkGroup()
        SUITE["distances"][Tdist]["ImageDistances.jl"] =
            @benchmarkable evaluate($dist, $imgA, $imgB)
        SUITE["distances"][Tdist]["Distances.jl"] =
            @benchmarkable Distances.evaluate($dist, $imgC, $imgD)
    end
end

add_distances_benchmarks!(SUITE)

###########
# Colwise #
###########

SUITE["colwise"] = BenchmarkGroup()

function add_colwise_benchmarks!(SUITE)

    imgA = testimage("cameraman")
    imgB = testimage("livingroom")

    alg = Otsu()
    imgA = binarize(imgA, alg)
    imgB = binarize(imgA, alg)

    imgC = Float32.(imgA)
    imgD = Float32.(imgB)

    # Two lists of images
    imgsA = [imgA, imgB, imgA, imgB]
    imgsB = [imgB, imgA, imgA, imgB]

    imgsC = [imgC, imgD, imgC, imgD]
    imgsD = [imgD, imgC, imgC, imgD]

    dists = create_distances()

    for dist in (dists)
        Tdist = typeof(dist)
        SUITE["colwise"][Tdist] = BenchmarkGroup()
        SUITE["colwise"][Tdist]["ImageDistances.jl"] =
            @benchmarkable colwise($dist, $imgsA, $imgsB)
        SUITE["colwise"][Tdist]["Distances.jl"] =
            @benchmarkable Distances.colwise($dist, $imgsC, $imgsD)
    end
end

add_colwise_benchmarks!(SUITE)

############
# Pairwise #
############

SUITE["pairwise"] = BenchmarkGroup()

function add_pairwise_benchmarks!(SUITE)

    imgA = testimage("cameraman")
    imgB = testimage("livingroom")

    alg = Otsu()
    imgA = binarize(imgA, alg)
    imgB = binarize(imgA, alg)

    imgC = Float32.(imgA)
    imgD = Float32.(imgB)

    # Two lists of images
    imgsA = [imgA, imgB, imgA, imgB]
    imgsB = [imgB, imgA, imgA, imgB]

    imgsC = [imgC, imgD, imgC, imgD]
    imgsD = [imgD, imgC, imgC, imgD]
    dists = create_distances()

    for dist in (dists)
        Tdist = typeof(dist)
        SUITE["pairwise"][Tdist] = BenchmarkGroup()
        SUITE["pairwise"][Tdist]["ImageDistances.jl"] =
            @benchmarkable pairwise($dist, $imgsA, $imgsB)
        SUITE["pairwise"][Tdist]["Distances.jl"] =
            @benchmarkable Distances.pairwise($dist, $imgsC, $imgsD)
    end
end

add_pairwise_benchmarks!(SUITE)


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#    Tests of default tables with Markdown.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

@testset "Markdown" begin
    data = [1 md"[This is a link](https://ronanarraes.com)"
            2 md"**Bold** *italics*"]

    expected = """

    Col. 1
    Col. 2

    1

This is a link

    2

Bold italics

"""

    result = pretty_table(String, data, backend = :html, standalone = false)

    @test result == expected
end


text:
@doc raw"""
    AbstractAffineConnection

Abstract type for affine connections on a manifold.
"""
abstract type AbstractAffineConnection end

"""
    LeviCivitaConnection

The [Levi-Civita connection](https://en.wikipedia.org/wiki/Levi-Civita_connection) of a Riemannian manifold.
"""
struct LeviCivitaConnection <: AbstractAffineConnection end

"""
    IsConnectionManifold <: AbstractTrait

Specify that a certain decorated Manifold is a connection manifold in the sence that it provides
explicit connection properties, extending/changing the default connection properties of a manifold.
"""
struct IsConnectionManifold <: AbstractTrait end

"""
    IsDefaultConnection{G<:AbstractAffineConnection}

Specify that a certain [`AbstractAffineConnection`](@ref) is the default connection for a manifold.
This way the corresponding [`ConnectionManifold`](@ref) falls back to the default methods
of the manifold it decorates.
"""
struct IsDefaultConnection{C<:AbstractAffineConnection} <: AbstractTrait
    connection::C
end
parent_trait(::IsDefaultConnection) = IsConnectionManifold()

"""
    ConnectionManifold{𝔽,,M<:AbstractManifold{𝔽},G<:AbstractAffineConnection} <: AbstractDecoratorManifold{𝔽}

# Constructor

    ConnectionManifold(M, C)

Decorate the [`AbstractManifold`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/types.html#ManifoldsBase.AbstractManifold)  `M` with [`AbstractAffineConnection`](@ref) `C`.
"""
struct ConnectionManifold{𝔽,M<:AbstractManifold{𝔽},C<:AbstractAffineConnection} <:
       AbstractDecoratorManifold{𝔽}
    manifold::M
    connection::C
end

function Base.filter(f, t::TraitList)
    if f(t.head)
        return merge_traits(t.head, filter(f, t.tail))
    else
        return filter(f, t.tail)
    end
end
Base.filter(f, t::EmptyTrait) = t

function active_traits(f, M::ConnectionManifold, args...)
    return merge_traits(
        is_default_connection(M.manifold, M.connection) ?
        IsDefaultConnection(M.connection) : EmptyTrait(),
        IsConnectionManifold(),
        filter(x -> x isa IsGroupManifold, active_traits(f, M.manifold, args...)),
    )
end

@doc raw"""
    christoffel_symbols_first(
        M::AbstractManifold,
        p,
        B::AbstractBasis;
        backend::AbstractDiffBackend = default_differential_backend(),
    )

Compute the Christoffel symbols of the first kind in local coordinates of basis `B`.
The Christoffel symbols are (in Einstein summation convention)

````math
Γ_{ijk} = \frac{1}{2} \Bigl[g_{kj,i} + g_{ik,j} - g_{ij,k}\Bigr],
````

where ``g_{ij,k}=\frac{∂}{∂ p^k} g_{ij}`` is the coordinate
derivative of the local representation of the metric tensor. The dimensions of
the resulting multi-dimensional array are ordered ``(i,j,k)``.
"""
christoffel_symbols_first(::AbstractManifold, ::Any, B::AbstractBasis)
function christoffel_symbols_first(
    M::AbstractManifold,
    p,
    B::AbstractBasis;
    backend::AbstractDiffBackend=default_differential_backend(),
)
    ∂g = local_metric_jacobian(M, p, B; backend=backend)
    n = size(∂g, 1)
    Γ = allocate(∂g, Size(n, n, n))
    @einsum Γ[i, j, k] = 1 / 2 * (∂g[k, j, i] + ∂g[i, k, j] - ∂g[i, j, k])
    return Γ
end
@trait_function christoffel_symbols_first(
    M::AbstractDecoratorManifold,
    p,
    B::AbstractBasis;
    kwargs...,
)

@doc raw"""
    christoffel_symbols_second(
        M::AbstractManifold,
        p,
        B::AbstractBasis;
        backend::AbstractDiffBackend = default_differential_backend(),
    )

Compute the Christoffel symbols of the second kind in local coordinates of basis `B`.
For affine connection manifold the Christoffel symbols need to be explicitly implemented
while, for a [`MetricManifold`](@ref) they are computed as (in Einstein summation convention)

````math
Γ^{l}_{ij} = g^{kl} Γ_{ijk},
````

where ``Γ_{ijk}`` are the Christoffel symbols of the first kind
(see [`christoffel_symbols_first`](@ref)), and ``g^{kl}`` is the inverse of the local
representation of the metric tensor. The dimensions of the resulting multi-dimensional array
are ordered ``(l,i,j)``.
"""
function christoffel_symbols_second(
    M::AbstractManifold,
    p,
    B::AbstractBasis;
    backend::AbstractDiffBackend=default_differential_backend(),
)
    Ginv = inverse_local_metric(M, p, B)
    Γ₁ = christoffel_symbols_first(M, p, B; backend=backend)
    Γ₂ = allocate(Γ₁)
    @einsum Γ₂[l, i, j] = Ginv[k, l] * Γ₁[i, j, k]
    return Γ₂
end

@trait_function christoffel_symbols_second(
    M::AbstractDecoratorManifold,
    p,
    B::AbstractBasis;
    kwargs...,
)

@doc raw"""
    christoffel_symbols_second_jacobian(
        M::AbstractManifold,
        p,
        B::AbstractBasis;
        backend::AbstractDiffBackend = default_differential_backend(),
    )

Get partial derivatives of the Christoffel symbols of the second kind
for manifold `M` at `p` with respect to the coordinates of `B`, i.e.

```math
\frac{∂}{∂ p^l} Γ^{k}_{ij} = Γ^{k}_{ij,l}.
```

The dimensions of the resulting multi-dimensional array are ordered ``(i,j,k,l)``.
"""
christoffel_symbols_second_jacobian(::AbstractManifold, ::Any, B::AbstractBasis)
function christoffel_symbols_second_jacobian(
    M::AbstractManifold,
    p,
    B::AbstractBasis;
    backend::AbstractDiffBackend=default_differential_backend(),
)
    n = size(p, 1)
    ∂Γ = reshape(
        _jacobian(q -> christoffel_symbols_second(M, q, B; backend=backend), p, backend),
        n,
        n,
        n,
        n,
    )
    return ∂Γ
end
@trait_function christoffel_symbols_second_jacobian(
    M::AbstractDecoratorManifold,
    p,
    B::AbstractBasis;
    kwargs...,
)

"""
    connection(M::AbstractManifold)

Get the connection (an object of a subtype of [`AbstractAffineConnection`](@ref))
of [`AbstractManifold`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/types.html#ManifoldsBase.AbstractManifold)  `M`.
"""
connection(::AbstractManifold)

"""
    connection(M::ConnectionManifold)

Return the connection associated with [`ConnectionManifold`](@ref) `M`.
"""
connection(M::ConnectionManifold) = M.connection

decorated_manifold(M::ConnectionManifold) = M.manifold

default_retraction_method(M::ConnectionManifold) = default_retraction_method(M.manifold)

@doc raw"""
    exp(::TraitList{IsConnectionManifold}, M::AbstractDecoratorManifold, p, X)

Compute the exponential map on a manifold that [`IsConnectionManifold`](@ref) `M` equipped with
corresponding affine connection.

If `M` is a [`MetricManifold`](@ref) with a [`IsDefaultMetric`](@ref) trait,
this method falls back to `exp(M, p, X)`.

Otherwise it numerically integrates the underlying ODE, see [`solve_exp_ode`](@ref).
Currently, the numerical integration is only accurate when using a single
coordinate chart that covers the entire manifold. This excludes coordinates
in an embedded space.
"""
exp(::TraitList{IsConnectionManifold}, M::AbstractDecoratorManifold, p, X)

function exp!(::TraitList{IsConnectionManifold}, M::AbstractDecoratorManifold, q, p, X)
    return retract!(
        M,
        q,
        p,
        X,
        ODEExponentialRetraction(ManifoldsBase.default_retraction_method(M)),
    )
end

"""
    gaussian_curvature(M::AbstractManifold, p, B::AbstractBasis; backend::AbstractDiffBackend = default_differential_backend())

Compute the Gaussian curvature of the manifold `M` at the point `p` using basis `B`.
This is equal to half of the scalar Ricci curvature, see [`ricci_curvature`](@ref).
"""
gaussian_curvature(::AbstractManifold, ::Any, ::AbstractBasis)
function gaussian_curvature(M::AbstractManifold, p, B::AbstractBasis; kwargs...)
    return ricci_curvature(M, p, B; kwargs...) / 2
end
@trait_function gaussian_curvature(
    M::AbstractDecoratorManifold,
    p,
    B::AbstractBasis;
    kwargs...,
)

"""
    is_default_connection(M::AbstractManifold, G::AbstractAffineConnection)

returns whether an [`AbstractAffineConnection`](@ref) is the default metric on the manifold `M` or not.
This can be set by defining this function, or setting the [`IsDefaultConnection`](@ref) trait for an
[`AbstractDecoratorManifold`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/decorator.html#ManifoldsBase.AbstractDecoratorManifold).
"""
is_default_connection(M::AbstractManifold, G::AbstractAffineConnection)
@trait_function is_default_connection(
    M::AbstractDecoratorManifold,
    G::AbstractAffineConnection,
)
function is_default_connection(
    ::TraitList{IsDefaultConnection{C}},
    ::AbstractDecoratorManifold,
    ::C,
) where {C<:AbstractAffineConnection}
    return true
end
function is_default_connection(M::ConnectionManifold)
    return is_default_connection(M.manifold, M.connection)
end
is_default_connection(::AbstractManifold, ::AbstractAffineConnection) = false

function retract_exp_ode!(
    M::AbstractManifold,
    q,
    p,
    X,
    ::AbstractRetractionMethod,
    b::AbstractBasis,
)
    sol = solve_exp_ode(M, p, X; basis=b, dense=false)
    copyto!(q, sol)
    return q
end

"""
    ricci_tensor(M::AbstractManifold, p, B::AbstractBasis; backend::AbstractDiffBackend = default_differential_backend())

Compute the Ricci tensor, also known as the Ricci curvature tensor,
of the manifold `M` at the point `p` using basis `B`,
see [`https://en.wikipedia.org/wiki/Ricci_curvature#Introduction_and_local_definition`](https://en.wikipedia.org/wiki/Ricci_curvature#Introduction_and_local_definition).
"""
ricci_tensor(::AbstractManifold, ::Any, ::AbstractBasis)
function ricci_tensor(M::AbstractManifold, p, B::AbstractBasis; kwargs...)
    R = riemann_tensor(M, p, B; kwargs...)
    n = size(R, 1)
    Ric = allocate(R, Size(n, n))
    @einsum Ric[i, j] = R[l, i, l, j]
    return Ric
end
@trait_function ricci_tensor(M::AbstractDecoratorManifold, p, B::AbstractBasis; kwargs...)

@doc raw"""
    riemann_tensor(M::AbstractManifold, p, B::AbstractBasis; backend::AbstractDiffBackend=default_differential_backend())

Compute the Riemann tensor ``R^l_{ijk}``, also known as the Riemann curvature
tensor, at the point `p` in local coordinates defined by `B`. The dimensions of the
resulting multi-dimensional array are ordered ``(l,i,j,k)``.

The function uses the coordinate expression involving the second Christoffel symbol,
see [`https://en.wikipedia.org/wiki/Riemann_curvature_tensor#Coordinate_expression`](https://en.wikipedia.org/wiki/Riemann_curvature_tensor#Coordinate_expression)
for details.

# See also

[`christoffel_symbols_second`](@ref), [`christoffel_symbols_second_jacobian`](@ref)
"""
riemann_tensor(::AbstractManifold, ::Any, ::AbstractBasis)
function riemann_tensor(
    M::AbstractManifold,
    p,
    B::AbstractBasis;
    backend::AbstractDiffBackend=default_differential_backend(),
)
    n = size(p, 1)
    Γ = christoffel_symbols_second(M, p, B; backend=backend)
    ∂Γ = christoffel_symbols_second_jacobian(M, p, B; backend=backend) ./ n
    R = allocate(∂Γ, Size(n, n, n, n))
    @einsum R[l, i, j, k] =
        ∂Γ[l, i, k, j] - ∂Γ[l, i, j, k] + Γ[s, i, k] * Γ[l, s, j] - Γ[s, i, j] * Γ[l, s, k]
    return R
end
@trait_function riemann_tensor(M::AbstractDecoratorManifold, p, B::AbstractBasis; kwargs...)

@doc raw"""
    solve_exp_ode(
        M::AbstractConnectionManifold,
        p,
        X,
        B::AbstractBasis;
        backend::AbstractDiffBackend = default_differential_backend(),
        solver = AutoVern9(Rodas5()),
        kwargs...,
    )

Approximate the exponential map on the manifold by evaluating the ODE descripting the geodesic at 1,
assuming the default connection of the given manifold by solving the ordinary differential
equation

```math
\frac{d^2}{dt^2} p^k + Γ^k_{ij} \frac{d}{dt} p_i \frac{d}{dt} p_j = 0,
```

where ``Γ^k_{ij}`` are the Christoffel symbols of the second kind, and
the Einstein summation convention is assumed. The argument `solver` follows
the `OrdinaryDiffEq` conventions. `kwargs...` specify keyword
arguments that will be passed to `OrdinaryDiffEq.solve`.

Currently, the numerical integration is only accurate when using a single
coordinate chart that covers the entire manifold. This excludes coordinates
in an embedded space.

!!! note
    This function only works when
    [OrdinaryDiffEq.jl](https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl) is loaded with
    ```julia
    using OrdinaryDiffEq
    ```
"""
function solve_exp_ode(M, p, X; kwargs...)
    throw(
        ErrorException(
            """
            solve_exp_ode not implemented on $(typeof(M)) for point $(typeof(p)), vector $(typeof(X)).
            For a suitable default, enter `using OrdinaryDiffEq`.
            """,
        ),
    )
end


text:
"Write [`Record`](@ref)s to `IO` as string."
function Base.write(io::IO, records::Vector{<:Record}) #Note: catching non-isbits arrays and we assume the indexes have been bumpped and the open ends are correct.
    for record in records
        Base.write(io, record, '\n')
    end
    return position(io)
end

"Write [`Record`](@ref) to `IO` as string."
function Base.write(io::IO, record::Record)
    return write(io, convert(String, record))
end


text:
"""
    masked_crc32c(data)
"""
function masked_crc32c(data)
    x = UInt32(CRC32c.crc32c(data))
    return UInt32(((x >> 15) | UInt32(x << 17)) + 0xa282ead8)
end

function serialize_proto(data::Union{ProtoType, ProtoEnum})
    pb = PipeBuffer()
    writeproto(pb, data)
    pb.data
end

function serialize_proto(data::Any)
    pb = PipeBuffer()
    write(pb, data)
    pb.data
end


text:
__precompile__()

module FilePaths

using Compat

import Base: ==
export
    # Types
    AbstractPath,
    Path,
    PosixPath,
    WindowsPath,
    Mode,
    Status,

    # Methods
    anchor,
    cwd,
    drive,
    home,
    parts,
    root,
    hasparent,
    parents,
    filename,
    extension,
    extensions,
    exists,
    isabs,
    mode,
    created,
    modified,
    relative,
    move,
    remove,
    tmpname,
    tmpdir,
    mktmp,
    mktmpdir,
    chown,
    executable,
    readable,
    writable,
    raw,

    # Macros
    @p_str,

    # Constants
    READ,
    WRITE,
    EXEC

@static if VERSION < v"0.6.0-dev.2514"
    import Base: isexecutable
else
    export isexecutable
end

@compat abstract type AbstractPath <: AbstractString end

# Required methods for subtype of AbstractString
Base.endof(p::AbstractPath) = endof(String(p))
Base.next(p::AbstractPath, i::Int) = next(String(p), i)

# The following should be implemented in the concrete types
Base.String(path::AbstractPath) = error("`String not implemented")
parts(path::AbstractPath) = error("`parts` not implemented.")
root(path::AbstractPath) = error("`root` not implemented.")
drive(path::AbstractPath) = error("`drive` not implemented.")

Base.convert(::Type{AbstractPath}, x::AbstractString) = Path(x)

include("constants.jl")
include("utils.jl")
include("libc.jl")
include("mode.jl")
include("status.jl")
include("posix.jl")
include("windows.jl")
include("path.jl")
include("deprecates.jl")

end # end of module


text:
# This is a helper function that assures we have the
# proportion of positive cases that we specify.
function curate_data(y, X, n, pct)
    N = length(y)
    if n > N
        error("The input data must have more records than the size of the desired output")
    elseif sum(y)/n < pct
        println(sum(y)/n)
        error("The input data must have more positive cases than desired output percentage")
    end
    pos_indcs = find(y)
    neg_indcs = setdiff(1:N, pos_indcs)
    num_pos = round(Int, n * pct)
    num_neg = n - num_pos
    indcs = vcat(sample(pos_indcs, num_pos), sample(neg_indcs, num_neg))
    return (y[indcs], X[indcs, :])
end

function gen_varcov(p, ρ1, ρ2, ρ3)
    Σ = eye(p)
    Σ[2, 1], Σ[1, 2] = ρ1, ρ1
    Σ[3, 1], Σ[1, 3] = ρ2, ρ2
    Σ[3, 2], Σ[2, 3] = ρ3, ρ3
    Σ[4, 3], Σ[3, 4] = ρ3, ρ3
    Σ[5, 3], Σ[3, 5] = ρ2, ρ2
    Σ[5, 4], Σ[4, 5] = ρ1, ρ1
    return Σ
end

function fitmodels(y, X, train, ntrees, subsample, ξ)
    n = length(y)
    n_conditions = 5
    test = setdiff(1:n, train)
    perf = zeros(n_conditions, 5)

    for i = 1:n_conditions
        if i == 1
            model, coeffs = build_adaboost_stumps(y[train], X[train, :], ntrees, subsample)
        elseif i == 2
            model, coeffs = build_adaboost_stumps_prefer(y[train], X[train, :], ntrees, subsample)
        elseif i == 3
            model, coeffs = build_adaboost_stumps_weight(y[train], X[train, :], ntrees, subsample, ξ)
        elseif i == 4
            model, coeffs = build_adaboost_stumps_wghtpref(y[train], X[train, :], ntrees, subsample, ξ)
        elseif i == 5
            X_smoted, y_smoted  = smote(X[train, :], y[train], 5)
            model, coeffs = build_adaboost_stumps(y_smoted, X_smoted, ntrees, subsample)
        end
        train_err = adaboost_train_error(y[train], X[train, :], model, coeffs)
        y_hat::Array{Int,1} = apply_adaboost_stumps(model, coeffs, X[test, :])
        y_hat_train::Array{Int,1} = apply_adaboost_stumps(model, coeffs, X[train, :])

        # train_err, test_err, train_auc, test_auc, f1score
        perf[i, :] = vcat(1 - train_err, mean(y[test] .== y_hat), roc_auc_score(y[train], y_hat_train), roc_auc_score(y[test], y_hat), f1score(roc(y[test], y_hat)))
    end
    perf
end

# @code_warntype fitmodels([1, 0, 1], randn(3, 3), [true, false, true], 10, 0.7, 0.1)
function simdata(n, p, μ_err)
    Σ = gen_varcov(p, 0.6, 0.3, 0.1)
    mvn = MvNormal(ones(p), Σ)
    N = n + round(Int, n * 0.5)                      # add 50% to desired n, so we can use curate_data()
    X1 = rand(mvn, N)'

    ϵ = rand(Normal(μ_err, 0.5), N)
    β = vcat([-3, -4, 1, -5, -4], zeros(p - 5))
    η = X1*β + ϵ                                     # linear predictor w/ error
    pr = 1.0 ./ (1.0 + exp(-η))                      # inv-logit
    y1 = map(π -> rand(Binomial(1, π)), pr)          # simulate outcome variable
    return (y1, X1)
end

function runsim(n, p, μ_err, pct; ntrees = 100, subsample = 0.7, seed = round(Int, time()), ξ = 0.1)
    srand(seed)
    y = zeros(n)
    X = zeros(n, p)
    # We use this loop to ensure that data curation works. In rare
    # cases we get an error resulting from too few positive cases.
    retry = true
    while retry
        y1, X1 = simdata(n, p, μ_err)
        retry = false

        try
            y, X = curate_data(y1, X1, n, pct)
        catch
            retry = true
        end
    end

    ntrain = round(Int, n * 0.7)
    train = sample(1:n, ntrain, replace = false)

    perf = fitmodels(y, X, train, ntrees, subsample, ξ)
    #println("Proportion of minority class in training set: $(mean(y[train]))")
    return perf
end

# runsim(200, 10, 2, 0.1, seed = 111, ξ = 0.1)

function run_simulations(nsims, n, p, μ_err, pct, ntrees = 100, subsample = 0.7; ξ = 0.1)
    n_conditions = 5
    n_trials = n_conditions * nsims
    perf = zeros(n_trials, 7)

    step_size = n_conditions - 1                                            # for stepping through output mat.
    simnum = 1
    for i = 1:n_conditions:n_trials
        seed = round(Int, time())

        perf[i:i+step_size, 1] = simnum
        perf[i:i+step_size, 2] = linspace(1.0, n_conditions, n_conditions)  # indicates the AdaBoost variation
        perf[i:i+step_size, 3:end] = runsim(n,
                                            p,
                                            μ_err,
                                            pct,
                                            ntrees = ntrees,
                                            subsample = subsample,
                                            seed = seed,
                                            ξ = ξ)

        simnum += 1
    end
    return perf
end

# simresults = run_simulations(500, 200, 10, 2, 0.1, ξ = 0.1)

function colmeans(X)
    p = size(X, 2)
    res = zeros(p)
    for j = 1:p
        res[j] = mean(X[:, j])
    end
    res
end

function summarise_sims(X, n_conditions = 5)

    means = zeros(n_conditions, 6)
    boost_typ = 1.0
    for i = 1:n_conditions
        keep = find(x -> x == boost_typ, X[:, 2])
        means[i, 1] = boost_typ
        means[i, 2:end] = colmeans(X[keep, 3:end])
        boost_typ += 1.0
    end
    means
end

# summarise_sims(simresults)

max_stars_repo_name: paulstey/boosting

max_issues_repo_name: paulstey/boosting

max_forks_repo_name: paulstey/boosting


text:
using GLMakie
function peaks(; n = 49)
    x = LinRange(-3, 3, n)
    y = LinRange(-3, 3, n)
    a = 3 * (1 .- x').^2 .* exp.(-(x'.^2) .- (y .+ 1).^2)
    b = 10 * (x' / 5 .- x'.^3 .- y.^5) .* exp.(-x'.^2 .- y.^2)
    c = 1 / 3 * exp.(-(x' .+ 1).^2 .- y.^2)
    return (x, y, a .- b .- c)
end
x, y, z = peaks()
set_theme!(theme_dark())
fig = Figure(resolution = (1200,800))
ax1 = Axis(fig[1,1], aspect = 1)
ax2 = Axis3(fig[1,2]; aspect = (1,1,0.7), perspectiveness = 0.5)
axs = [ax1, ax2]
cmap =  :diverging_bkr_55_10_c35_n256
contourf!(axs[1], x, y, z; levels = 0.1:0.08:0.9, mode = :relative,
    colormap = cmap)
# bug, colormap cannot be transparent
contourf!(axs[2], x, y, z; levels = 14, colormap = cmap)
contour3d!(axs[2], x, y, z; levels = 14, colormap = cmap,
    transparency = true, linewidth = 5)
limits!(axs[1], -3,3,-3,3)
[hidedecorations!(ax; grid = false) for ax in axs]
save("contourf.png", fig)
fig


text:
# integration tests of entire models

facts("Conjugate Normal model") do
    dims = (20, 6)

    # note: it won't matter much how we initialize here
    μ[j] ~ Normal(zeros(dims[2]), ones(dims[2]))
    τ[j] ~ Gamma(1.1 * ones(dims[2]), ones(dims[2]))
    μ0[j] ~ Const(zeros(dims[2]))
    τ0[j] ~ Const(2 * ones(dims[2]))
    a0[j] ~ Const(1.1 * ones(dims[2]))
    b0[j] ~ Const(ones(dims[2]))

    y[i, j] ~ Const(rand(dims))

    # make factors
    obs = @factor LogNormalFactor y μ τ
    μ_prior = @factor LogNormalFactor μ μ0 τ0
    τ_prior = @factor LogGammaCanonFactor τ a0 b0

    m = VBModel([μ, τ, μ0, τ0, a0, b0, y], [obs, μ_prior, τ_prior])

    @fact Set([n.name for n in m.nodes]) --> Set([:μ, :τ, :μ0, :τ0, :a0, :b0, :y])
    @fact Set([typeof(f) for f in m.factors]) --> Set([VinDsl.LogNormalFactor{2}, VinDsl.LogNormalFactor{1}, VinDsl.LogGammaCanonFactor{1}])
    @fact length(m.graph) --> 7
    @fact check_conjugate(τ, m) --> true
    @fact check_conjugate(μ, m) --> true

    update!(m)
end

facts("Univariate ⟷ multivariate naturals extraction") do
    context("vector mean, full covariance") do
        d = 5
        μ[i] ~ MvNormalCanon(zeros(d), diagm(ones(d)))
        Λ[i, i] ~ Wishart(float(d), diagm(ones(d)))
        x[i, j] ~ MvNormalCanon([randn(d) for x in 1:20], [diagm(ones(d)) for x in 1:20])
        f = @factor LogMvNormalCanonFactor x μ Λ

        @fact value(f) --> isfinite
        @fact map(size, naturals(f, μ)[1]) --> ((d,), (d, d))
        @fact map(size, naturals(f, Λ)[1]) --> ((), (d, d))
    end

    context("vector mean, diagonal covariance") do
        d = 5
        μ[i] ~ MvNormalCanon(zeros(d), diagm(ones(d)))
        τ[i] ~ Gamma(1.1 * ones(d), ones(d))
        x[i, j] ~ MvNormalCanon([randn(d) for x in 1:20], [diagm(ones(d)) for x in 1:20])
        f = @factor LogMvNormalDiagCanonFactor x μ τ

        @fact value(f) --> isfinite
        @fact map(size, naturals(f, μ)[1]) --> ((d,), (d, d))
        @fact size(naturals(f, τ)) --> (d,)
    end

    context("vector mean, scalar covariance") do
        d = 5
        μ[i] ~ MvNormalCanon(zeros(d), diagm(ones(d)))
        τ ~ Gamma(1.1, 1.)
        x[i, j] ~ MvNormalCanon([randn(d) for x in 1:20], [diagm(ones(d)) for x in 1:20])
        f = @factor LogMvNormalDiagCanonFactor x μ τ

        @fact value(f) --> isfinite
        @fact map(size, naturals(f, μ)[1]) --> ((d,), (d, d))
        @fact map(size, naturals(f, τ)[1]) --> ((), ())
    end

    context("scalar mean, full covariance") do
        d = 5
        μ ~ Normal(0, 1)
        Λ[i, i] ~ Wishart(float(d), diagm(ones(d)))
        x[i, j] ~ MvNormalCanon([randn(d) for x in 1:20], [diagm(ones(d)) for x in 1:20])
        f = @factor LogMvNormalCanonFactor x μ Λ

        @fact value(f) --> isfinite
        @fact map(size, naturals(f, μ)[1]) --> ((), ())
        @fact map(size, naturals(f, Λ)[1]) --> ((), (d, d))
    end

    context("scalar mean, diagonal covariance") do
        d = 5
        μ ~ Normal(0, 1)
        τ[i] ~ Gamma(1.1 * ones(d), ones(d))
        x[i, j] ~ MvNormalCanon([randn(d) for x in 1:20], [diagm(ones(d)) for x in 1:20])
        f = @factor LogMvNormalDiagCanonFactor x μ τ

        @fact value(f) --> isfinite
        @fact map(size, naturals(f, μ)[1]) --> ((), ())
        @fact map(size, naturals(f, τ)[1]) --> ((), ())
    end

    context("vector-of-scalars mean, diagonal covariance") do
        d = 5
        N = 20
        μ[i] ~ Normal(zeros(d), ones(d))
        τ[i] ~ Gamma(1.1 * ones(d), ones(d))
        x[i, j] ~ MvNormalCanon([randn(d) for x in 1:N], [diagm(ones(d)) for x in 1:N])
        f = @factor LogMvNormalDiagCanonFactor x μ τ

        @fact value(f) --> isfinite
        @fact size(naturals(f, μ)) --> (d,)
        @fact size(naturals(f, τ)) --> (d,)
        @fact size(naturals(f, x)) --> (N,)
        @fact map(size, naturals(f, μ)[1]) --> ((), ())
        @fact map(size, naturals(f, τ)[1]) --> ((), ())
        @fact map(size, naturals(f, x)[1]) --> ((d,), (d, d))
    end

    context("scalar mean, scalar covariance") do
        d = 5
        N = 20
        μ ~ Normal(0, 1)
        τ ~ Gamma(1.1, 1)
        x[i, j] ~ MvNormalCanon([randn(d) for x in 1:N], [diagm(ones(d)) for x in 1:N])
        f = @factor LogMvNormalDiagCanonFactor x μ τ

        @fact value(f) --> isfinite
        @fact map(size, naturals(f, μ)[1]) --> ((), ())
        @fact map(size, naturals(f, τ)[1]) --> ((), ())
        @fact size(naturals(f, x)) --> (N,)
        @fact map(size, naturals(f, x)[1]) --> ((d,), (d, d))
    end
end

facts("Basic Hidden Markov Model") do
    d = 5
    T = 100
    pars = [Dirichlet(d, 1) for i in 1:d]

    A_par = rand(MarkovMatrix(pars))
    π0_par = rand(Dirichlet(d, 1))
    ψ_par = rand(d, T)

    z[i, t] ~ HMM(ψ_par, π0_par, A_par)
    A ~ MarkovMatrix(pars)
    A_0 ~ Const(A_par)
    π0 ~ Dirichlet(π0_par)
    π0_0 ~ Const(π0_par)
    ψ[i, t] ~ Const(rand(d, T))

    f = @factor LogMarkovChainFactor z π0 A
    π_prior = @factor LogDirichletFactor π0 π0_0
    A_prior = @factor LogMarkovMatrixFactor A A_0
    π_nats = naturals(f, π0)
    A_nats = naturals(f, A)
    z_nats = naturals(f, z)

    @fact value(f) --> isfinite
    @fact value(π_prior) --> isfinite
    @fact value(A_prior) --> isfinite
    @fact map(size, naturals(f, π0)[1]) --> ((d,), )
    @fact map(size, naturals(f, A)[1]) --> ((d, d), )
    @fact map(size, naturals(f, z)[1]) --> ((d, T), (d,), (d, d))
    @fact map(size, naturals(π_prior, π0)[1]) --> ((d,), )
    @fact map(size, naturals(A_prior, A)[1]) --> ((d, d), )
end

facts("Linear combination expression node") do
    context("Univariate") do
        dims = (20, 6)
        μ[j] ~ Normal(zeros(dims[2]), ones(dims[2]))
        ν[j] ~ Normal(zeros(dims[2]), ones(dims[2]))
        τ[j] ~ Gamma(1.1 * ones(dims[2]), ones(dims[2]))
        @exprnode w (μ + ν)

        y[i, j] ~ Const(rand(dims))

        # make factors
        obs = @factor LogNormalFactor y w τ
        @fact value(obs) --> isfinite
    end

    context("Multivariate") do
        d = 5
        μ[i] ~ MvNormalCanon(zeros(d), diagm(ones(d)))
        ν[i] ~ MvNormalCanon(zeros(d), diagm(ones(d)))
        Λ[i, i] ~ Wishart(float(d), diagm(ones(d)))
        x[i, j] ~ MvNormalCanon([randn(d) for x in 1:20], [diagm(ones(d)) for x in 1:20])
        @exprnode w (μ + ν)
        f = @factor LogMvNormalCanonFactor x w Λ

        @fact value(f) --> isfinite
    end

    context("Matrix linear combination") do
        d = 5
        q = 7
        μ[i, k] ~ MvNormalCanon([zeros(d) for _ in 1:q], [diagm(ones(d)) for _ in 1:q])
        ν[k] ~ MvNormalCanon(zeros(q), diagm(ones(q)))
        Λ[i, i] ~ Wishart(float(d), diagm(ones(d)))
        x[i, j] ~ MvNormalCanon([randn(d) for _ in 1:20], [diagm(ones(d)) for _ in 1:20])
        @exprnode w (dot(μ[k], ν[k]))
        f = @factor LogMvNormalCanonFactor x w Λ

        # @fact value(f) --> isfinite
    end
end

facts("Updating via explicit optimization") do
    dims = (20, 6)

    # note: it won't matter much how we initialize here
    μ[j] ~ Normal(zeros(dims[2]), ones(dims[2]))
    τ[j] ~ Gamma(1.1 * ones(dims[2]), ones(dims[2]))
    μ0[j] ~ Const(zeros(dims[2]))
    τ0[j] ~ Const(2 * ones(dims[2]))
    a0[j] ~ Const(1.1 * ones(dims[2]))
    b0[j] ~ Const(ones(dims[2]))

    y[i, j] ~ Const(rand(dims))

    # make factors
    obs = @factor LogNormalFactor y μ τ
    μ_prior = @factor LogNormalFactor μ μ0 τ0
    τ_prior = @factor LogGammaCanonFactor τ a0 b0

    m = VBModel([μ, τ, μ0, τ0, a0, b0, y], [obs, μ_prior, τ_prior])

    context("L-BFGS") do
        update!(μ, m, Val{:l_bfgs})
    end
end


text:
################################################################################
# Mutability
"""
    Mutability
The Mutability type is an abstract trait type with children Mutable, Immutable,
and MutabilityUnknown.
"""
abstract type Mutability end
struct Mutable <: Mutability end
struct Immutable <: Mutability end
struct UnknownMutability <: Mutability end
const _MUT_TYPE = Mutable()
const _IMM_TYPE = Immutable()
const _UNK_TYPE = UnknownMutability()
"""
    mutability(obj)
Yields an object of type Mutable, Immutable, or UnknownMutability depending on
whether the given type object is mutable, immutable, or unknown.
"""
mutability(T::Type) = begin
    if typeof(T) === Union
        let mA = mutability(T.a)
            return (mA === mutability(T.b) ? mA : _UNK_TYPE)
        end
    else
        return _UNK_TYPE
    end
end
# Implementations for mutability below
mutability(T::DataType) = begin
    if isbitstype(T)
        return _IMM_TYPE
    elseif !T.isconcretetype
        return _UNK_TYPE
    elseif T.mutable
        return _MUT_TYPE
    else
        return _IMM_TYPE
    end
end
mutability(::Core.TypeofBottom) = _UNK_TYPE
mutability(T::UnionAll) = mutability(T{T.var.ub})
# these two come out wrong in the above code, according to the documentation:
mutability(::Type{String}) = _IMM_TYPE
mutability(::Type{Symbol}) = _IMM_TYPE

################################################################################
# #ismut, #isimm, #isimmtype, #ismuttype functions for testing immutabilities
_ismuttype(::Mutable) = true
_ismuttype(::Immutable) = false
_ismuttype(::UnknownMutability) = false
"""
    ismuttype(T)
Yields true if the given type T is known to be mutable and false otherwise.
"""
ismuttype(T::Type) = _ismuttype(mutability(T))
_isimmtype(::Mutable) = false
_isimmtype(::Immutable) = true
_isimmtype(::UnknownMutability) = false
"""
    isimmtype(T)
Yields true if the given type T is known to be immutable and false otherwise.
"""
isimmtype(T::Type) = _isimmtype(mutability(T))
"""
    ismut(x)
Yields true if x is mutable and false if x is immutable or has unknown mutability.
"""
ismut(::T) where {T} = _ismuttype(mutability(T))
"""
    isimm(x)
Yields true if x is immutable and false if x is mutable or has unknown mutability.
"""
isimm(::T) where {T} = _isimmtype(mutability(T))

################################################################################
# #mutable_cotype #immutable_cotype #abstract_cotype

_mutable_cotype(::Mutable, x) = error(
    "Mutable type $(typeof(x)) cannot have a mutable cotype")
_mutable_cotype(::Immutable, x) = error(
    "Immutable type $(typeof(x)) has not overloaded mutable_cotype()")
_mutable_cotype(::UnknownMutability, x) = error(
    "Type $(typeof(x)) with unknown mutability does not have a mutable cotype")
"""
    mutable_cotype(T)
Yields the mutable equivalent of the immutable or abstract type T. If T is
already a mutable type or there is no mutable cotype of T then an error is
raised.
"""
mutable_cotype(T::Type) = _mutable_cotype(mutability(T))

_immutable_cotype(::Immutable, x) = error(
    "Immutable type $(typeof(x)) cannot have an immutable cotype")
_immutable_cotype(::Mutable, x) = error(
    "Mutable type $(typeof(x)) has not overloaded immutable_cotype()")
_immmutable_cotype(::UnknownMutability, x) = error(
    "Type $(typeof(x)) with unknown mutability does not have an immutable cotype")
"""
    immutable_cotype(T)
Yields the immutable equivalent of the mutable or abstract type T. If T is
already an immutable type or there is no immutable cotype of T then an error is
raised.
"""
immutable_cotype(T::Type) = _immutable_cotype(mutability(T), x)

_abstract_cotype(::Mutable, x) = error(
    "Mutable type $(typeof(x)) has not overloaded abstract_cotype()")
_abstract_cotype(::Immutable, x) = error(
    "Immutable type $(typeof(x)) has not overloaded abstract_cotype()")
_abstract_cotype(::UnknownMutability, x) = error(
    "Type $(typeof(x)) with unknown mutability cannot have an abstract cotype")
"""
    abstract_cotype(T)
Yields the abstract equivalent of the immutable or mutable type T. If T is
already an abstract type or there is no mutable cotype of T then an error is
raised.
"""
abstract_cotype(T::Type) = _abstract_cotype(mutability(T), x)

################################################################################
# #aimtypes
"""
    @aimtypes A I M
Declares that the given types are related types; specifically:
 * A is an abstract type,
 * I <: A is an immutable type,
 * M <: A is a mutable type,
 * freeze(x::M) yields an object of type I, and
 * thaw(x::I) yields an object of type M.
"""
macro aimtypes(A::Symbol, I::Symbol, M::Symbol)
    return quote
        $(esc(:abstract_cotype))(x::Type{$I}) = $A
        $(esc(:abstract_cotype))(x::Type{$M}) = $A
        $(esc(:immutable_cotype))(x::Type{$A}) = $I
        $(esc(:immutable_cotype))(x::Type{$M}) = $I
        $(esc(:mutable_cotype))(x::Type{$A}) = $M
        $(esc(:mutable_cotype))(x::Type{$I}) = $M
        $(esc(:mutability))(::Type{$I}) = Immutable()
        $(esc(:mutability))(::Type{$M}) = Mmutable()
        $(esc(:copy))(x::$I) = x

    end
end

################################################################################
# #freeze, #thaw, #deepfreeze, #deepthaw

_freeze(::Immutable, x) = x
_freeze(::UnknownMutability, x) = error(
    "type $(typeof(x)) has not overloaded the freeze() method")
_freeze(::Mutable, x) = error(
    "type $(typeof(x)) has not overloaded the freeze() method")
"""
    freeze(x)
Yields x if x is immutable or a shallow immutable copy of x if x is not
immutable. If x is of a non-immutable type that has not overloaded the freeze()
method, then an error is raised.
"""
freeze(x) = _freeze(mutability(x), x)

"""
    deepfreeze(x)
Yields a deep copy of x in which x and all of its sub-values have been converted
from mutable into immutable objects. Immutable objects are left as-is unless
their children are changed.

Like the deepcopy function, deepfreeze() can be customized; however, unlike
deepcopy() which employs a deepcopy_internal() co-routine, deepfreeze() should
be overloaded directly for the particular type.
"""
deepfreeze(@nospecialize(x)) = begin
    T = typeof(x)::DataType
    isbitstype(T) && return x
    # There are a few possibilities here:
    # (1) the object is mutable, in which case we'll need to freeze() it
    # (2) the object is immutable, in which case it's fine except maybe its
    #     children;
    # (3) the object may have children whose identities change, in which case
    #     we need to make a new object even if it's immutable
    # The easiest case to catch is the mutable case, in which case we freeze it
    # then deepfreeze the result (to ensure that children are also frozen).
    isimmtype(T) || return deepfreeze(freeze(x))
    # For an immutable object, we start by freezing the children/fields:
    n = nfields(x)
    changed = false
    field_defined = Vector{Bool}()
    field_changed = Vector{Bool}()
    new_field_val = Vector{Any}()
    for 1 in 1:n
        if isdefined(x, i)
            old_obj = getfield(x, i)
            new_obj = deepfreeze(old_obj)
            push!(new_field_val, new_obj)
            push!(field_changed, old_obj !== new_obj)
            push!(field_defined, true)
            old_obj === new_obj || changed = true
        else
            push!(new_field_val, nothing)
            push!(field_changed, false)
            push!(field_defined, false)
        end
    end
    # We now have the updated fields stored as these three arrays; if the
    # field was not defined in the old object then it should not be set in any
    # new object.
    if !changed:
        # In this case, none of the children have changed and the object is
        # known to be immutable; we do not even have to call freeze on it!
        return x
    else
        # The object is immutable, but we need to freeze its children:
        y = ccall(:jl_new_struct_uninit, Any, (Any,), T)
        ftypes = fieldtypes(T)
        for i in 1:n
            isdef = field_defined[i]
            isnew = field_changed[i]
            new_val = new_field_val[i]
            if isdef
                # if it's a new object and the types don't match, we should
                # throw an error; otherwise we can set it
                if isnew && !(typeof(new_val) <: ftypes[i])
                    error("field $i of type $T was not frozen into" *
                          " a sub-type of $(ftypes[i])")
                else
                    ccall(:jl_set_nth_field,
                          Cvoid, (Any, Csize_t, Any), y, i-1,
                          new_val)
                end
            end
        end
        return y
    end
end
const TERMINAL_TYPES = Union{Symbol,String,Core.MethodInstance,Method,GlobalRef,
                             DataType,Union,UnionAll,Task,Regex}
deepfreeze(x::TERMINAL_TYPES) = x
deepfreeze(x::Tuple) = ntuple(i->deepfreeze(x[i]), length(x))
deepfreeze(x::Module) = error("deepfreeze of Modules not supported")
deepfreeze(x::Nothing) = nothing

_thaw(::Mutable, x) = x
_thaw(::UnknownMutability, x) = error(
    "type $(typeof(x)) has not overloaded the freeze() method")
_thaw(::Immutable, x) = error(
    "type $(typeof(x)) has not overloaded the freeze() method")
"""
    thaw(x)
Yields a mutable version of the object x if x is an immutable object; otherwise
if x is already mutable yields x unchanged. To get a mutable copy of an object
regardless of whether it is immutable, use thawcopy(). If x cannot be converted
to an immutable object, an error is raised.
"""
thaw(x::T) = _thaw(mutability(T), x)

"""
    deepthaw(x)
Yields a mutable version of the object x after also deepthawing any relevant
children of x. Note that deepthaw(x) should never affect x or any of its
children whether x is mutable or not (i.e., when the children of x change, a
copy of x with the new children is returned even if x is already mutable).

Unlike thaw(), freeze(), and deepfreeze(), deepthaw() does not raise errors when
it encounters immutable objects that it cannot thaw. This is because it is
presumed that immutable terminals (like String and Int64) should be left as they
are.

If you wish to obtain a mutable deepcopy of an object regardless of whether the
original object is mutable or immutable, use the deepthawcopy() function.
"""
deepthaw(@nospecialize(x)) = begin
    T = typeof(x)::DataType
    isbitstype(T) && return x
    # There are a few possibilities here:
    # (1) the object is immutable, in which case we'll need to thaw() it
    # (2) the object is mutable, in which case it's fine except maybe its
    #     children;
    # (3) the object may have children whose identities change when we deepthaw
    #     them, in which case we need to copy the object then update
    # The easiest case to catch is the immutable case, in which case we thaw it
    # then deepthaw the result (to ensure that children are also deepthawed).
    if !ismuttype(T)
        try
            x = thaw(x)
        catch e
            return x
        end
    end
    # We start by deepthaw'ing the children/fields:
    n = nfields(x)
    changed = false
    field_defined = Vector{Bool}()
    field_changed = Vector{Bool}()
    new_field_val = Vector{Any}()
    for 1 in 1:n
        if isdefined(x, i)
            old_obj = getfield(x, i)
            new_obj = deepthaw(old_obj)
            push!(new_field_val, new_obj)
            push!(field_changed, old_obj !== new_obj)
            push!(field_defined, true)
            old_obj === new_obj || changed = true
        else
            push!(new_field_val, nothing)
            push!(field_changed, false)
            push!(field_defined, false)
        end
    end
    # We now have the updated fields stored as these three arrays; if the
    # field was not defined in the old object then it should not be set in any
    # new object.
    if !changed:
        # In this case, none of the children have changed and the object is
        # known to be mutable already
        return x
    else
        # The object is mutable, but we need to update its children:
        y = ccall(:jl_new_struct_uninit, Any, (Any,), T)
        ftypes = fieldtypes(T)
        for i in 1:n
            isdef = field_defined[i]
            isnew = field_changed[i]
            new_val = new_field_val[i]
            if isdef
                # if it's a new object and the types don't match, we should
                # throw an error; otherwise we can set it
                if isnew && !(typeof(new_val) <: ftypes[i])
                    error("field $i of type $T was not thawed into" *
                          " a sub-type of $(ftypes[i])")
                else
                    ccall(:jl_set_nth_field,
                          Cvoid, (Any, Csize_t, Any), y, i-1,
                          new_val)
                end
            end
        end
        return y
    end
end
deepthaw(x::TERMINAL_TYPES) = x
deepthaw(x::Tuple) = ntuple(i->deepthaw(x[i]), length(x))
deepthaw(x::Module) = x
deepthaw(x::Nothing) = nothing

_thawcopy(::Mutable, x) = copy(x)
_thawcopy(::Immutable, x) = thaw(x)
_thawcopy(::UnknownMutability, x) = copy(x)
"""
    thawcopy(x)
Yields thaw(x) if x is immutable and copy(x) if x is mutable.
"""
thawcopy(x::T) where {T} = _thawcopy(mutability(T), x)

"""
    deepthawcopy(x)
Yields deepthaw(deepcopy(x)).
"""
deepthawcopy(x) = deepthaw(deepcopy(x))

max_issues_repo_name: noahbenson/FunctionalUtilities

max_forks_repo_name: noahbenson/FunctionalUtilities


text:
"""
            Functions to read a snapshot and store it in a dictionary.

    Author: Ludwig Böss
    Contact: lboess@usm.lmu.de
    Created: 2018-12-12

"""

max_stars_repo_name: LudwigBoess/GadJet

max_issues_repo_name: LudwigBoess/GadJet

max_forks_repo_name: LudwigBoess/GadJet


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl

see description in native.jl
=#

using SparseArrays

struct MatrixCompletionJuMP{T <: Real} <: ExampleInstanceJuMP{T}
    size_ratio::Int
    num_rows::Int
end

function build(inst::MatrixCompletionJuMP{T}) where {T <: Float64}
    (size_ratio, num_rows) = (inst.size_ratio, inst.num_rows)
    @assert size_ratio >= 1
    num_cols = size_ratio * num_rows

    (rows, cols, Avals) = findnz(sprandn(num_rows, num_cols, 0.8))
    is_known = vec(Matrix(sparse(rows, cols,
        trues(length(Avals)), num_rows, num_cols)))

    model = JuMP.Model()
    JuMP.@variable(model, X[1:length(is_known)])
    JuMP.@variable(model, t)
    JuMP.@objective(model, Min, t)
    JuMP.@constraint(model, vcat(t, X) in
        MOI.NormSpectralCone(num_rows, num_cols))
    X_unknown = X[.!is_known]
    JuMP.@constraint(model, vcat(1, X_unknown) in
        MOI.GeometricMeanCone(1 + length(X_unknown)))
    JuMP.@constraint(model, X[is_known] .== Avals)

    return model
end


text:
## julia-shell-emacstools.jl --- Helper functions for julia-shell.el
#
# Author: Dennis Ogbe <dogbe@purdue.edu>
#
## Usage:
# This file should not be modified, it is loaded from emacs on the startup of a
# julia shell
#
## TODO:
# -  precompile the module on first load?
#
## License:
#
# Copyright (c) 2015 Dennis Ogbe
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.

#
# PART 1: The EmacsTools Module
#
###############################################################################

module EmacsTools

export get_completions, get_latex_symbols

# print completions of SUBSTRING
function get_completions(substring::AbstractString)
  c = Base.REPLCompletions.completions(substring, length(substring));
  for completion in c[1]
    println(completion);
  end
end

# print all available LaTeX symbols
function get_latex_symbols()
  c = Base.REPLCompletions.latex_symbols;
  for (latex, unicode) in sort!(collect(c), by=x->x[2])
    println("$latex $unicode")
  end
end

end # module

#
# PART 2: Load the EmacsTools module and do other customizations
#
###############################################################################

import EmacsTools


text:
using Mimi

@defcomp climatedynamics begin
    TATM = Variable(index=[time]) # Increase in temperature of atmosphere (degrees C from 1900)
    TOCEAN = Variable(index=[time]) # Increase in temperature of lower oceans (degrees C from 1900)

    FORC = Parameter(index=[time]) # Increase in radiative forcing (watts per m2 from 1900)
    fco22x = Parameter() # Forcings of equilibrium CO2 doubling (Wm-2)
    t2xco2 = Parameter() # Equilibrium temp impact (oC per doubling CO2)
    tatm0 = Parameter() # Initial atmospheric temp change (C from 1900)
    tatm1 = Parameter() # Period 2 atmospheric temp change (C from 1900)
    tocean0 = Parameter() # Initial lower stratum temp change (C from 1900)

    # Transient TSC Correction ("Speed of Adjustment Parameter")
    c1 = Parameter() # Climate equation coefficient for upper level
    c3 = Parameter() # Transfer coefficient upper to lower stratum
    c4 = Parameter() # Transfer coefficient for lower level

    function run_timestep(p, v, d, t)
        #Define function for TATM
        if is_first(t)
            v.TATM[t] = p.tatm0
        elseif t.t == 2
            v.TATM[t] = p.tatm1
        else
            v.TATM[t] = v.TATM[t-1] + p.c1 * ((p.FORC[t] - (p.fco22x/p.t2xco2) * v.TATM[t-1]) - (p.c3 * (v.TATM[t-1] - v.TOCEAN[t-1])))
        end

        #Define function for TOCEAN
        if is_first(t)
            v.TOCEAN[t] = p.tocean0
        else
            v.TOCEAN[t] = v.TOCEAN[t-1] + p.c4 * (v.TATM[t-1] - v.TOCEAN[t-1])
        end
    end
end


text:
# discrete MDP 1
ns = 4
na = 6
R = zeros(ns, na)
T = zeros(ns, na, ns) # sp, a, s

# from state 1
T[2,1, 1] = 1.0  # if go, move to 2
T[1,2:na,1] .= 1.0 # else stay

# from state 2
T[1,2,2] = 0.7 # if (s2, safe) end up in 1
T[3,2,2] = 0.3 # if (s2, safe) end up in 3

T[3,3,2] = 0.5
T[4,3,2] = 0.5
T[2,1,2] = 1.
T[2,4:end,2] .= 1.

# from state 3
T[3,:,3] .= 1

# from state 4
T[4,:,4] .= 1.
T[:,4,4] .= 0.
T[1,4,4] = 1.

mdp = TabularMDP(T, R, 0.95)

MEC = maximal_end_components(mdp, verbose=true)
@test MEC == [[3], [4]]

gw = SimpleGridWorld(size=(20,20))
gw_mec = maximal_end_components(gw, verbose=false)

@test length(gw_mec) == 2*gw.size_y + 2*gw.size_x - 3


text:
module Resource

# Point the REST APIs to to actual functions that will do the tasks
# parsing the URI to translate the string to arguments
# it also Responsesible for create a new users and validate them
# and log them in to the system

using Dates, HTTP, JSON3
using ..Model, ..Service, ..Auth, ..Contexts, ..Workers

const ROUTER = HTTP.Router()

HTTP.@register(ROUTER, "POST", "/album", createAlbum)
createAlbum(req) = Service.createAlbum(JSON3.read(req.body))::Album

HTTP.@register(ROUTER, "GET", "/album/*", getAlbum)
getAlbum(req) = Service.getAlbum(parse(Int, HTTP.URIs.splitpath(req.target)[2]))::Album

HTTP.@register(ROUTER, "PUT", "/album/*", updateAlbum)
updateAlbum(req) = Service.updateAlbum(parse(Int, HTTP.URIs.splitpath(req.target)[2]), JSON3.read(req.body, Album))::Album

HTTP.@register(ROUTER, "DELETE", "/album/*", deleteAlbum)
deleteAlbum(req) = Service.deleteAlbum(parse(Int, HTTP.URIs.splitpath(req.target)[2]))

HTTP.@register(ROUTER, "GET", "/", pickAlbumToListen)
pickAlbumToListen(req) = fetch(Workers.@async(Service.pickAlbumToListen()::Album))

function contextHandler(req)
    withcontext(User(req)) do
        HTTP.Response(200, JSON3.write(HTTP.handle(ROUTER, req)))
    end
end

const AUTH_ROUTER = HTTP.Router(contextHandler)

function authenticate(user::User)
    resp = HTTP.Response(200, JSON3.write(user))
    return Auth.addtoken!(resp, user)
end

HTTP.@register(AUTH_ROUTER, "POST", "/user", createUser)
createUser(req) = authenticate(Service.createUser(JSON3.read(req.body))::User)

HTTP.@register(AUTH_ROUTER, "POST", "/user/login", loginUser)
loginUser(req) = authenticate(Service.loginUser(JSON3.read(req.body, User))::User)

function requestHandler(req)
    start = Dates.now(Dates.UTC)
    @info (timestamp=start, event="ServiceRequestBegin", tid=Threads.threadid(),
            method=req.method, target=req.target)
    local resp
    try
        resp = HTTP.handle(AUTH_ROUTER, req)
    catch e
        if e isa Auth.Unauthenticated
            resp = HTTP.Response(401)
        else
            s = IOBuffer()
            showerror(s, e, catch_backtrace(); backtrace=true)
            errormsg = String(resize!(s.data, s.size))
            @error errormsg
            resp = HTTP.Response(500, errormsg)
        end
    end
    stop = Dates.now(Dates.UTC)
    @info (timestamp=stop, event="ServiceRequestEnd", tid=Threads.threadid(),
            method=req.method, target=req.target, duration=Dates.value(stop - start),
            status=resp.status, bodysize=length(resp.body))
    return resp
end

function run()
    HTTP.serve(requestHandler, "0.0.0.0", 8080)
end

end # module

max_stars_repo_name: vanhp/MicroServiceJulia

max_issues_repo_name: vanhp/MicroServiceJulia

max_forks_repo_name: vanhp/MicroServiceJulia


text:
function _projection(projection_idx, E::AbstractMatrix{N}, e::AbstractVector{N},
                             backend=nothing,
                             algorithm=nothing,
                             prune=true) where {N<:Real}
    require(:Polyhedra; fun_name="_projection")
    require(:CDDLib; fun_name="_projection")
    if backend == nothing
        backend = LazySets.default_cddlib_backend(N)
    end

    if algorithm == nothing
        algorithm = Polyhedra.FourierMotzkin()
    elseif !(algorithm <: EliminationAlgorithm)
        error("the algorithm $algorithm is not a valid elimination algorithm;
              choose among any of $(subtypes(Polyhedra.EliminationAlgorithm))")
    end
    eliminiation_idx = setdiff(1:size(E,2), projection_idx)
    PQ = HPolyhedron(E, e)
    PQ_cdd = polyhedron(PQ, backend=backend)
    W = HPolyhedron(Polyhedra.eliminate(PQ_cdd, eliminiation_idx, algorithm))
    if prune
        success = LazySets.remove_redundant_constraints!(W)
        if !success
            error("the constraints corresponding to the minkowski sum of the given " *
                  "sets are infeasible")
        end
    end
    return W
end

max_stars_repo_name: mforets/InvariantSets

max_issues_repo_name: mforets/InvariantSets


text:
# SOCP Primal Structs for Simple SOCP problems

include("../constraints/constraints.jl")

"""
    SOCP_primals(x, s, t)

Struct to store the primal variables of the form [x; s; t]

`x` and `s` are mx1 and nx1 vectors.
`t` is a real number

See also: [`primalVec`](@ref) and [`primalStruct`](@ref)
"""
mutable struct SOCP_primals
    x
    s
    t
end

"""
    primalVec(y::SOCP_primals)

Converts [`SOCP_primals`](@ref) struct to [x; s; t] vector.
"""
function primalVec(y::SOCP_primals)
    return [y.x; y.s; y.t]
end

"""
    primalStruct(v, xSize::Int64, sSize::Int64, tSize::Int64)

Converts [x; s; t] vector to [`SOCP_primals`](@ref) struct
"""
function primalStruct(v, xSize::Int64, sSize::Int64, tSize::Int64)
    return SOCP_primals(v[1:xSize], v[xSize+1:xSize+sSize], v[end])
end

"""
    getXVals(yList::Array{SOCP_primals, 1})

Returns the "x" values of and an array of [`SOCP_primals`](@ref).
Used primarily with plotting.
"""
function getXVals(yList::Array{SOCP_primals, 1})
    xList = []
    for xst in yList
        push!(xList, xst.x)
    end
    return xList
end

@doc raw"""
    getViolation(yList::Array{SOCP_primals, 1}, c::AL_coneSlack)

returns how much each SOCP constraint is violated
(in 3 mulitdimensional arrays).

Recall that an SOCP is generally
```math
||Ax-b|| ≤ c^⊤ x - d
```

Using slack variables, we have

```math
||s|| - t ≤ 0
```

```math
(Ax - b) - s = 0
```

```math
(c^⊤ x - d) - t = 0
```

This function returns the constraint violation for each of the above 3 in 3
separate arrays.
"""
function getViolation(yList::Array{SOCP_primals, 1}, c::AL_coneSlack)
    coneViolation = []
    affineViolation = []
    lastViolation = []
    for y in yList
        vio = getNormToProjVals(c, y.x, y.s, y.t)
        push!(coneViolation, vio[1])
        push!(affineViolation, vio[2:1+size(y.s, 1)])
        push!(lastViolation, vio[end])
    end
    return coneViolation, affineViolation, lastViolation
end

max_stars_repo_name: danineamati/TrajOptSOCPs


text:
mutable struct WishartPriorXMLElement <: MyXMLElement
    el::XMLOrNothing
    df::Int
    scale_mat::AbstractMatrix{Float64}
    param::MatrixParameter
    id::String

    function WishartPriorXMLElement(scale_mat::AbstractMatrix{Float64}, df::Int,
                                    ref_param::MatrixParameter)
        check_posdef(scale_mat)
        p = size(scale_mat, 1)
        if df < p
            error("Degrees of freedom must be at least the matrix dimension.")
        end
        return new(nothing, df, scale_mat, ref_param, ref_param.id * ".prior")
    end
end

function WishartPriorXMLElement(param::MatrixParameter)
    p = size(param, 1)
    return WishartPriorXMLElement(Diagonal(ones(p)), p, param)
end

function WishartPriorXMLElement(mat::AbstractMatrix{Float64},
                                param::MatrixParameter)
    if size(mat) != size(param)
        error("Incompatible dimensions.")
    end
    return WishartPriorXMLElement(mat, size(mat, 1), param)
end

function make_xml(wp::WishartPriorXMLElement)
    make_xml(wp.param)
    el = new_element(bn.MULTIVARIATE_WISHART_PRIOR)
    set_attributes(el, [(bn.ID, wp.id), (bn.DF, string(wp.df))])
    scale_el = new_child(el, bn.SCALE_MATRIX)
    add_matrix_parameter(scale_el, wp.scale_mat)
    data_el = new_child(el, bn.DATA)
    add_ref_el(data_el, wp.param.el)
    wp.el = el
    return el
end

# function make_Wishart_prior(scale::AbstractArray{Float64, 2},
#                             mp_el::XMLElement,
#                             id::String)

#     el = new_element(bn.MULTIVARIATE_WISHART_PRIOR)
#     set_attributes(el, [(bn.ID, id), (bn.DF, string(size(scale, 1)))])
#     scale_el = new_child(el, bn.SCALE_MATRIX)
#     add_matrix_parameter(scale_el, scale)
#     data_el = new_child(el, bn.DATA)
#     add_ref_el(data_el, mp_el)
#     return el
# end


text:
module KernelTools

export @loophoist, @tile, @test_looporder

if VERSION.minor < 3
    using Cartesian
else
    using Base.Cartesian
end

include("edge.jl")
include("macros.jl")

end # module