text:
abstract type AbstractPricer end

include("analytic.jl")
include("carr_madan.jl")
include("lewis.jl")


text:
import Pkg;
Pkg.add("CSV")
Pkg.add("DataFrames")
Pkg.add("JuMP")
Pkg.add("HiGHS")

using Documenter
using SupplyChainOptimization

makedocs(
    sitename = "SupplyChainOptimization",
    format = Documenter.HTML(),
    modules = [SupplyChainOptimization],
    pages = ["index.md",
            "Examples" => ["optimization flows.md", "optimization locations.md", "multi-period optimization.md", "adding special constraints.md", "inventory movements.md"],
            "Internals" => ["optimization model.md"],
            "API" => ["reference.md"],
            "Sponsor" => ["sponsor.md"]
            ]
)

deploydocs(;
    repo="https://github.com/SupplyChef/SupplyChainOptimization.jl",
    devbranch = "main"
)

ext: jl


text:
using CSV
using InteractiveUtils

function main()
    dataframe = CSV.read("sample.csv", header=false, delim=',')
    @show dataframe
    row,col=size(dataframe)
    for  c =1:col
        dataframe[c] *= c
    end
    dataframe |> CSV.write("output.csv",delim=',',writeheader=false)
    run(pipeline(`cat output.tsv`))
end

main()

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/NodeNumbering.jl/blob/master/LICENSE
"""
    element_adjacencies

The Dict shows some element types and their node adjacencies. For example the element type Quad4 has four nodes/element so that node No. 1 is adjacent to nodes 2 and 4 etc.
"""
const element_adjacencies = Dict(
                            :Quad4 => [[2, 4], [1, 3], [2, 4], [1, 3]],
                            :Tri3 => [[2,3], [1,3], [1,2]],
                            :Tet4 => [[2,3,4], [1,3,4], [1,2,4], [1,2,3]]);

"""
    create_adjacency_graph(elements::Dict{Int, Vector{Int}}, element_types::Dict{Int, Symbol})

Create a Dict that shows all nodes and their adjacencies.

# Examples

```jldoctest
julia> elements = Dict(
                  1 => [9, 1, 8, 4],
                  2 => [1, 3, 2, 8],
                  3 => [8, 2, 7, 5],
                  4 => [2, 6, 7]);

julia> element_types = Dict(
                       1 => :Quad4,
                       2 => :Quad4,
                       3 => :Quad4,
                       4 => :Tri3);

julia> create_adjacency_graph(elements, element_types)
Dict{Int64,Array{Int64,1}} with 9 entries:
  7 => [2, 6, 5]
  9 => [1, 4]
  4 => [9, 8]
  2 => [6, 7, 3, 8]
  3 => [1, 2]
  8 => [1, 2, 5, 4]
  5 => [8, 7]
  6 => [2, 7]
  1 => [3, 8, 9]

```

# References

* Wikipedia contributors. "Adjacency list". Wikipedia, The Free Encyclopedia. Wikipedia, The Free Encyclopedia, 7 Jun. 2017. Web. 17 Jul. 2017. https://en.wikipedia.org/wiki/Adjacency_list
"""
function create_adjacency_graph(elements::Dict{Int, Vector{Int}}, element_types::Dict{Int, Symbol})
    neighbours = Dict{Int, Vector{Int}}()
    for (k, nodes) in elements
        neighbour_indices = element_adjacencies[element_types[k]]
        for (i, node) in enumerate(nodes)
            node_neighbours = nodes[neighbour_indices[i]]
            if haskey(neighbours, node)
                 for nodei in node_neighbours
                    if (nodei in neighbours[node]) == false
                       push!(neighbours[node], nodei)
                    end
                end
            else
                neighbours[node] = node_neighbours
            end
        end
    end
    return neighbours
end


text:
module ImageInTerminal

using Crayons
using ColorTypes
using ImageCore
using ImageTransformations

export

    colorant2ansi,
    imshow,
    imshow256,
    imshow24bit

include("colorant2ansi.jl")
include("encodeimg.jl")
include("imshow.jl")

# -------------------------------------------------------------------
# overload default show in the REPL for colorant (arrays)

const colormode = TermColorDepth[TermColor256()]

"""
    use_256()

Triggers `imshow256` automatically if an array of colorants is to
be displayed in the julia REPL. (This is the default)
"""
use_256() = (colormode[1] = TermColor256())

"""
    use_24bit()

Triggers `imshow24bit` automatically if an array of colorants is to
be displayed in the julia REPL.
Call `ImageInTerminal.use_256()` to restore default behaviour.
"""
use_24bit() = (colormode[1] = TermColor24bit())

# colorant arrays
function Base.show(
        io::IO, ::MIME"text/plain",
        img::AbstractVecOrMat{<:Colorant})
    println(io, summary(img), ":")
    ImageInTerminal.imshow(io, img, colormode[1])
end

# colorant
function Base.show(io::IO, ::MIME"text/plain", color::Colorant)
    fgcol = _colorant2ansi(color, colormode[1])
    chr = _charof(alpha(color))
    print(io, Crayon(foreground = fgcol), chr, chr, " ")
    print(io, Crayon(foreground = :white), color)
    print(io, Crayon(reset = true))
end

end # module


text:
export Term

struct Term{C, T} <: AbstractTerm{T}
    α::T
    x::Monomial{C}
end

iscomm(::Type{Term{C, T}}) where {C, T} = C

Base.convert(::Type{Term{C, T}}, t::Term{C, T}) where {C, T} = t
function Base.convert(::Type{Term{C, T}}, t::Term{C}) where {C, T}
    return Term{C, T}(T(t.α), t.x)
end
Term(t::Term) = t

function Base.convert(::Type{Term{C, T}}, x::Monomial{C}) where {C, T}
    return Term{C, T}(one(T), x)
end
Term{C}(x::Monomial{C}) where C = convert(Term{C, Int}, x)
Term(x::Monomial{C}) where C = Term{C}(x)

function Base.convert(::Type{Term{C, T}}, x::PolyVar{C}) where {C, T}
    return convert(Term{C, T}, convert(Monomial{C}, x))
end
Term{C}(x::PolyVar{C}) where C = Term{C}(convert(Monomial{C}, x))
Term(x::PolyVar{C}) where C = Term{C}(x)

function MP.convertconstant(::Type{Term{C, T}}, α) where {C, T}
    return Term{C}(convert(T, α))
end
Term{C}(α::T) where {C, T} = Term{C, T}(α, Monomial{C}())

Base.broadcastable(t::Term) = Ref(t)
#(::Type{TermContainer{C}}){C}(x::PolyVar{C}) = Term(x)
#(::Type{TermContainer{C}}){C}(x::Monomial{C}) = Term(x)
#(::Type{TermContainer{C}}){C}(t::TermContainer{C}) = t
#TermContainer(x::PolyVar) = Term(x)
#TermContainer(x::Monomial) = Term(x)
#TermContainer(t::TermContainer) = t

#Base.convert{C, T}(::Type{TermContainer{C, T}}, x::Union{Monomial{C},PolyVar{C}}) = Term{C, T}(x)
#Base.convert{C, T}(::Type{TermContainer{C, T}}, α::T) = Term{C, T}(α, Monomial{C}())
#Base.convert{C, S, T}(::Type{TermContainer{C, T}}, α::S) = TermContainer{C, T}(T(α))
#(::Type{TermContainer{C}}){C, T}(α::T) = TermContainer{C, T}(α)

#Base.convert{C, T}(::Type{TermContainer{C, T}}, t::Term{C}) = Term{C, T}(t)

Base.copy(t::T) where {T<:Term} = T(copy(t.α), copy(t.x))

MP.coefficient(t::Term) = t.α
MP.monomial(t::Term) = t.x
_vars(t) = _vars(t.x)


text:
@testset "Sampling" begin
  @testset "RegularSampling" begin
    b = Box(P2(0, 0), P2(2, 2))
    ps = sample(b, RegularSampling(3))
    @test collect(ps) == P2[(0,0),(1,0),(2,0),(0,1),(1,1),(2,1),(0,2),(1,2),(2,2)]
    ps = sample(b, RegularSampling(2, 3))
    @test collect(ps) == P2[(0,0),(2,0),(0,1),(2,1),(0,2),(2,2)]

    s = Sphere(P2(0, 0), T(2))
    ps = sample(s, RegularSampling(4))
    ts = P2[(2,0),(0,2),(-2,0),(0,-2)]
    for (p, t) in zip(ps, ts)
      @test p ≈ t
    end

    s = Sphere(P3(0, 0, 0), T(2))
    ps = sample(s, RegularSampling(2,2))
    ts = P3[(1.7320508075688772, 0.0, 1.0),
            (1.7320508075688772, 0.0, -1.0),
            (-1.7320508075688772, 0.0, 1.0),
            (-1.7320508075688772, 0.0, -1.0)]
    for (p, t) in zip(ps, ts)
      @test p ≈ t
    end

    b = Ball(P2(0, 0), T(2))
    ps = sample(b, RegularSampling(4,3))
    ts = P2[(1.0, 0.0), (0.0, 1.0), (-1.0, 0.0), (0.0, -1.0),
            (1.5, 0.0), (0.0, 1.5), (-1.5, 0.0), (0.0, -1.5),
            (2.0, 0.0), (0.0, 2.0), (-2.0, 0.0), (0.0, -2.0)]
    for (p, t) in zip(ps, ts)
      @test p ≈ t
    end

    b = Ball(P3(0, 0, 0), T(2))
    ps = sample(b, RegularSampling(3,2,3))
    ts = P3[(0.7071067811865475, 0.0, 0.7071067811865476),
            (1.0, 0.0, 6.123233995736766e-17),
            (0.7071067811865476, 0.0, -0.7071067811865475),
            (-0.7071067811865475, 8.659560562354932e-17, 0.7071067811865476),
            (-1.0, 1.2246467991473532e-16, 6.123233995736766e-17),
            (-0.7071067811865476, 8.659560562354934e-17, -0.7071067811865475),
            (1.0606601717798212, 0.0, 1.0606601717798214),
            (1.5, 0.0, 9.184850993605148e-17),
            (1.0606601717798214, 0.0, -1.0606601717798212),
            (-1.0606601717798212, 1.2989340843532398e-16, 1.0606601717798214),
            (-1.5, 1.8369701987210297e-16, 9.184850993605148e-17),
            (-1.0606601717798214, 1.29893408435324e-16, -1.0606601717798212),
            (1.414213562373095, 0.0, 1.4142135623730951),
            (2.0, 0.0, 1.2246467991473532e-16),
            (1.4142135623730951, 0.0, -1.414213562373095),
            (-1.414213562373095, 1.7319121124709863e-16, 1.4142135623730951),
            (-2.0, 2.4492935982947064e-16, 1.2246467991473532e-16),
            (-1.4142135623730951, 1.7319121124709868e-16, -1.414213562373095)]
    for (p, t) in zip(ps, ts)
      @test p ≈ t
    end
  end
end


text:
# basic graph tests for trivial cases

using Test
using Statistics
using IncrementalInference

##

@testset "test basic utility functions" begin

@test incrSuffix(:x45_4) == :x45_5

@test incrSuffix(:x45, +3) == :x48

@test incrSuffix(:x45_4, -1) == :x45_3

end

@testset "Test basic single variable graph with one prior..." begin

fg = initfg()

addVariable!(fg, :x0, ContinuousScalar)
addFactor!(fg, [:x0;], Prior(Normal(0.0,1.0)))

# test solved flag
@test getSolvedCount(fg, :x0) == 0
@test !isSolved(getVariable(fg, :x0))

# run solver once
tree = solveTree!(fg)

@test getSolvedCount(fg, :x0) == 1
@test isSolved(fg, :x0)

tree = solveTree!(fg)

@test getSolvedCount(fg, :x0) == 2
@test isSolved(fg, :x0)

# check mean and covariance
@test (getBelief(fg, :x0) |> getKDEMean .|> abs)[1] < 0.5
pts_ = getPoints(getBelief(fg, :x0))
TensorCast.@cast pts[i,j] := pts_[j][i]
@test 0.3 < Statistics.cov( pts[1,:] ) < 1.9

# test free solvable variables (occurs in fixed-/ clique recycling)
addVariable!(fg, :x1, ContinuousScalar, solvable=1)

solveTree!(fg, storeOld=true)

@test getSolvable(fg, :x1) == 0

end

@testset "Test basic single variable graph with one prior offset by 1000..." begin

fg = initfg()

addVariable!(fg, :x0, ContinuousScalar)
addFactor!(fg, [:x0;], Prior(Normal(1000.0,1.0)))

tree = solveTree!(fg)

# check mean and covariance
@test abs((getBelief(fg, :x0) |> getKDEMean)[1]-1000) < 0.5
pts_ = getPoints(getBelief(fg, :x0))
TensorCast.@cast pts[i,j] := pts_[j][i]
@test 0.4 < Statistics.cov( pts[1,:] ) < 1.8

end

@testset "Test basic single variable graph with two identical priors..." begin

fg = initfg()

addVariable!(fg, :x0, ContinuousScalar)
addFactor!(fg, [:x0;], Prior(Normal(0.0,1.0)))
addFactor!(fg, [:x0;], Prior(Normal(0.0,1.0)))

tree = solveTree!(fg)

# check mean and covariance
@test (getBelief(fg, :x0) |> getKDEMean .|> abs)[1] < 0.4
# should be sqrt(0.5) = 0.7, but lands near 0.6 instead -- computation is too confident.
pts_ = getPoints(getBelief(fg, :x0))
TensorCast.@cast pts[i,j] := pts_[j][i]
@test 0.3 < Statistics.cov( pts[1,:] ) < 1.0

end

@testset "Test basic single variable graph with three identical priors..." begin

fg = initfg()

addVariable!(fg, :x0, ContinuousScalar)
addFactor!(fg, [:x0;], Prior(Normal(0.0,1.0)))
addFactor!(fg, [:x0;], Prior(Normal(0.0,1.0)))
addFactor!(fg, [:x0;], Prior(Normal(0.0,1.0)))

tree = solveTree!(fg)

# check mean and covariance
@test (getBelief(fg, :x0) |> getKDEMean .|> abs)[1] < 0.4
# should be sqrt(1/3) = 0.577, but lands near 0.35 instead -- computation is too confident.
pts_ = getPoints(getBelief(fg, :x0))
TensorCast.@cast pts[i,j] := pts_[j][i]
@test 0.1 < Statistics.cov( pts[1,:] ) < 0.75

end

@testset "Test basic single variable graph with two priors at + and - 1..." begin

fg = initfg()

addVariable!(fg, :x0, ContinuousScalar)
addFactor!(fg, [:x0;], Prior(Normal(-1.0,1.0)))
addFactor!(fg, [:x0;], Prior(Normal(+1.0,1.0)))

tree = solveTree!(fg)

# check mean and covariance -- should be zero
@test (getBelief(fg, :x0) |> getKDEMean .|> abs)[1] < 0.8
# should be sqrt(1/2) = 0.707 -- computation results nearer 0.7.
pts_ = getPoints(getBelief(fg, :x0))
TensorCast.@cast pts[i,j] := pts_[j][i]
@test 0.2 < Statistics.cov( pts[1,:] ) < 1.5

end

@testset "Test basic single variable graph with two priors at + and - 1, offset by -1000..." begin

fg = initfg()

addVariable!(fg, :x0, ContinuousScalar)
addFactor!(fg, [:x0;], Prior(Normal(-1.0-1000,1.0)))
addFactor!(fg, [:x0;], Prior(Normal(+1.0-1000,1.0)))

tree = solveTree!(fg)

# check mean and covariance -- should be zero
@test abs((getBelief(fg, :x0) |> getKDEMean)[1] + 1000) < 0.6
# should be sqrt(1/2) = 0.707 -- computation results nearer 0.7.
pts_ = getPoints(getBelief(fg, :x0))
TensorCast.@cast pts[i,j] := pts_[j][i]
@test 0.2 < Statistics.cov( pts[1,:] ) < 1.1

end

@testset "Test basic two variable graph with two identical priors and weak connection..." begin

fg = initfg()

addVariable!(fg, :x0, ContinuousScalar)
addVariable!(fg, :x1, ContinuousScalar)
addFactor!(fg, [:x0;], Prior(Normal(0.0,1.0)))
addFactor!(fg, [:x1;], Prior(Normal(0.0,1.0)))
addFactor!(fg, [:x0;:x1;], LinearRelative(Normal(0.0,10.0)))

tree = solveTree!(fg)

# check mean and covariance -- should be zero
@test (getBelief(fg, :x0) |> getKDEMean .|> abs)[1] < 0.6
@test (getBelief(fg, :x1) |> getKDEMean .|> abs)[1] < 0.6

pts_ = getPoints(getBelief(fg, :x0))
TensorCast.@cast pts[i,j] := pts_[j][i]
@test 0.4 < Statistics.cov( pts[1,:] ) < 2.3
pts_ = getPoints(getBelief(fg, :x1))
TensorCast.@cast pts[i,j] := pts_[j][i]
@test 0.4 < Statistics.cov( pts[1,:] ) < 2.4

end

@testset "Test basic two variable graph with two separated priors and weak connection..." begin

fg = initfg()

addVariable!(fg, :x0, ContinuousScalar)
addVariable!(fg, :x1, ContinuousScalar)
addFactor!(fg, [:x0;], Prior(Normal(-1.0,1.0)))
addFactor!(fg, [:x1;], Prior(Normal(+1.0,1.0)))
addFactor!(fg, [:x0;:x1;], LinearRelative(Normal(0.0,10.0)))

tree = solveTree!(fg)

# check mean and covariance -- should be near each prior
@test abs((getBelief(fg, :x0) |> getKDEMean)[1]+1) < 0.75
@test abs((getBelief(fg, :x1) |> getKDEMean)[1]-1) < 0.75

pts_ = getPoints(getBelief(fg, :x0))
TensorCast.@cast pts[i,j] := pts_[j][i]
@test 0.3 < Statistics.cov( pts[1,:] ) < 2.5

pts_ = getPoints(getBelief(fg, :x1))
TensorCast.@cast pts[i,j] := pts_[j][i]
@test 0.3 < Statistics.cov( pts[1,:] ) < 2.5

end

@testset "Test basic two variable graph with two separated priors and strong connection..." begin

fg = initfg()

addVariable!(fg, :x0, ContinuousScalar)
addVariable!(fg, :x1, ContinuousScalar)
addVariable!(fg, :x2, ContinuousScalar)
addFactor!(fg, [:x0;], Prior(Normal(-1.0,1.0)))
addFactor!(fg, [:x2;], Prior(Normal(+1.0,1.0)))
addFactor!(fg, [:x0;:x1;], LinearRelative(Normal(0.0,1.0)))
addFactor!(fg, [:x1;:x2;], LinearRelative(Normal(0.0,1.0)))

tree = solveTree!(fg)

# check mean and covariance -- should between two priors somewhere
@test abs((getBelief(fg, :x0) |> getKDEMean)[1] + 1) < 0.9
@test abs((getBelief(fg, :x1) |> getKDEMean)[1]) < 0.9
@test abs((getBelief(fg, :x2) |> getKDEMean)[1] - 1) < 0.9

pts_ = getPoints(getBelief(fg, :x0))
TensorCast.@cast pts[i,j] := pts_[j][i]
@test 0.3 < Statistics.cov( pts[1,:] ) < 1.8

pts_ = getPoints(getBelief(fg, :x1))
TensorCast.@cast pts[i,j] := pts_[j][i]
@test 0.3 < Statistics.cov( pts[1,:] ) < 2.0

pts_ = getPoints(getBelief(fg, :x2))
TensorCast.@cast pts[i,j] := pts_[j][i]
@test 0.3 < Statistics.cov( pts[1,:] ) < 2.2

end

@testset "Test basic five variable graph with two separated priors and nominal connection..." begin

fg = initfg()

addVariable!(fg, :x0, ContinuousScalar)
addVariable!(fg, :x1, ContinuousScalar)
addVariable!(fg, :x2, ContinuousScalar)
addVariable!(fg, :x3, ContinuousScalar)
addVariable!(fg, :x4, ContinuousScalar)
addFactor!(fg, [:x0;], Prior(Normal(-3.0,1.0)))
addFactor!(fg, [:x4;], Prior(Normal(+3.0,1.0)))
addFactor!(fg, [:x0;:x1;], LinearRelative(Normal(0.0,1.0)))
addFactor!(fg, [:x1;:x2;], LinearRelative(Normal(0.0,1.0)))
addFactor!(fg, [:x2;:x3;], LinearRelative(Normal(0.0,1.0)))
addFactor!(fg, [:x3;:x4;], LinearRelative(Normal(0.0,1.0)))

# #1196
drawGraph(fg, filepath="testgraphplot/myfg.dot", show=false)

tree = solveTree!(fg, storeOld=true)

# using KernelDensityEstimatePlotting
# plotKDE((x->getBelief(fg,x)).([:x0;:x1;:x2;:x3;:x4]))
# using Gadfly, Cairo, Fontconfig
# drawTree(tree,show=true,imgs=true)

# check mean and covariance -- should be zero
X0 = (getBelief(fg, :x0) |> getKDEMean)[1]
X1 = (getBelief(fg, :x1) |> getKDEMean)[1]
X2 = (getBelief(fg, :x2) |> getKDEMean)[1]
X3 = (getBelief(fg, :x3) |> getKDEMean)[1]
X4 = (getBelief(fg, :x4) |> getKDEMean)[1]

@test X0 < X1 < X2 < X3 < X4

@test abs(X0+X4) < 2.2
@test abs(X1+X3) < 2.2
@test abs(X2) < 2.2

pts_ = getPoints(getBelief(fg, :x0))
TensorCast.@cast pts[i,j] := pts_[j][i]
@test 0.2 < Statistics.cov( pts[1,:] ) < 2.8

pts_ = getPoints(getBelief(fg, :x1))
TensorCast.@cast pts[i,j] := pts_[j][i]
@test 0.2 < Statistics.cov( pts[1,:] ) < 2.9

pts_ = getPoints(getBelief(fg, :x2))
TensorCast.@cast pts[i,j] := pts_[j][i]
@test 0.2 < Statistics.cov( pts[1,:] ) < 3.0

pts_ = getPoints(getBelief(fg, :x3))
TensorCast.@cast pts[i,j] := pts_[j][i]
@test 0.2 < Statistics.cov( pts[1,:] ) < 3.1

pts_ = getPoints(getBelief(fg, :x4))
TensorCast.@cast pts[i,j] := pts_[j][i]
@test 0.2 < Statistics.cov( pts[1,:] ) < 3.2

@testset "Test localProduct on solveKey" begin

localProduct(fg,:x2)

localProduct(fg,:x2, solveKey=:graphinit)

end

end

@testset "Test graph reset to init..." begin

fg = initfg()

addVariable!(fg, :x0, ContinuousScalar)
addFactor!(fg, [:x0;], Prior(Normal(1000.0,1.0)))

initAll!(fg)

# init values before solve
pts_ = getPoints(getBelief(fg, :x0))
TensorCast.@cast pts[i,j] := pts_[j][i]
X0 = pts |> deepcopy

tree = solveTree!(fg)

# values after solve
pts_ = getPoints(getBelief(fg, :x0))
TensorCast.@cast pts[i,j] := pts_[j][i]
X0s = pts

@test 1e-10 < norm(X0 - X0s)

resetInitialValues!(fg)

pts_ = getPoints(getBelief(fg, :x0))
TensorCast.@cast pts[i,j] := pts_[j][i]
X0reset = pts |> deepcopy

@test norm(X0 - X0reset) < 1e-10

end

#


text:
using Test, DiffEqOperators
using OrdinaryDiffEq

@testset "Parabolic Heat Equation with Dirichlet BCs" begin
    x = collect(-pi : 2pi/511 : pi)
    u0 = @. -(x - 0.5).^2 + 1/12
    @. u_analytic(x)  = -(x .- 0.5).^2 + 1/12
    A = CenteredDifference(2,2,2π/511,512);
    bc = DirichletBC(u_analytic(-pi - 2pi/511),u_analytic(pi + 2pi/511))
    step(u,p,t) = A*bc*u
    heat_eqn = ODEProblem(step, u0, (0.,10.))
    soln = solve(heat_eqn,Tsit5(),dense=false,tstops=0:0.01:10)

    for t in 0:0.1:10
        @test soln(t)[1] ≈ u0[1] rtol=0.05
        @test soln(t)[end] ≈ u0[end] rtol=0.05
    end
end

@testset "Parabolic Heat Equation with Neumann BCs" begin
    N = 512
    dx = 2π/(N-1)
    x = collect(-pi : dx : pi)
    u0 = @. -(x - 0.5)^2 + 1/12
    B = CenteredDifference(1,2,dx,N)
    deriv_start, deriv_end = (B*u0)[1], (B*u0)[end]

    A = CenteredDifference(2,2,dx,N)
    bc = NeumannBC((deriv_start,deriv_end),dx,1)

    step(u,p,t) = A*bc*u
    heat_eqn = ODEProblem(step, u0, (0.,10.))
    soln = solve(heat_eqn,Tsit5(),dense=false,tstops=0:0.01:10)

    first_order_coeffs_start = [-11/6, 3.0, -3/2, 1/3] * (1/dx)
    first_order_coeffs_end = -reverse([-11/6, 3.0, -3/2, 1/3] * (1/dx))

    for t in 0:0.1:10
        @test sum(first_order_coeffs_start .* soln(t)[1:4]) ≈ deriv_start atol=1e-1
        @test sum(first_order_coeffs_end .* soln(t)[end-3:end]) ≈ deriv_end atol=1e-1
    end
end

@testset "Parabolic Heat Equation with Robin BCs" begin
    N = 512
    dx = 2π/(N-1)
    x = collect(-pi : dx : pi)
    u0 = @. -(x - 0.5)^2 + 1/12
    B = CenteredDifference(1,2,dx,N);
    deriv_start, deriv_end = (B*u0)[1], (B*u0)[end]
    params = [1.0,0.5]

    left_RBC = params[1]*u0[1] - params[2]*deriv_start
    right_RBC = params[1]*u0[end] + params[2]*deriv_end

    A = CenteredDifference(2,2,dx,N);
    bc = RobinBC((params[1],-params[2],left_RBC), (params[1],params[2],right_RBC),dx,1);
    step(u,p,t) = A*bc*u
    heat_eqn = ODEProblem(step, u0, (0.,10.));
    soln = solve(heat_eqn,Tsit5(),dense=false,tstops=0:0.01:10);

    first_order_coeffs_start = [-11/6, 3.0, -3/2, 1/3] * (1/dx)
    first_order_coeffs_end = -reverse([-11/6, 3.0, -3/2, 1/3] * (1/dx))
    val = []

    for t in 0.2:0.1:9.8
        @test params[1]*soln(t)[1] - params[2]*sum(first_order_coeffs_start .* soln(t)[1:4]) ≈ left_RBC atol=1e-1
        # append!(val,params[1]*soln(t)[1] + -params[2]*sum(first_order_coeffs_start .* soln(t)[1:4]) - left_RBC)
        @test params[1]*soln(t)[end] + params[2]*sum(first_order_coeffs_end .* soln(t)[end-3:end]) ≈ right_RBC atol=1e-1
    end
end


text:
abstract type EventOrMessage end

include("Messages.jl")
include("Events.jl")

StructTypes.StructType(::Type{<:EventOrMessage}) = StructTypes.Struct()
StructTypes.StructType(::Type{EventOrMessage}) = StructTypes.AbstractType()
StructTypes.subtypekey(::Type{EventOrMessage}) = :type
StructTypes.subtypes(::Type{EventOrMessage}) = (; Events.event_types..., Messages.message_types...)


text:
module FVReconstructions
export AbstractReconstruction
using KernelAbstractions.Extras: @unroll
import StaticArrays: SUnitRange, SVector

"""
    AbstractReconstruction

Supertype for FV reconstructions.

Concrete types must provide implementions of
    - `width(recon)`
       returns the width of the reconstruction. Total number of points used in
       reconstruction of top and bottom states is `2width(recon) + 1`
       - (::AbstractReconstruction)(state_top, state_bottom, cell_states::SVector,
                                    cell_weights)
      compute the reconstruction

```
(::AbstractReconstruction)(
        state_top,
        state_bottom,
        cell_states,
        cell_weights,
    )
```

Perform the finite volume reconstruction for the top and bottom states using the
tuple of `cell_states` values using the `cell_weights`.
"""
abstract type AbstractReconstruction end
function (::AbstractReconstruction) end

"""
    AbstractSlopeLimiter

Supertype for FV slope limiter

Given two values `Δ1` and `Δ2`, should return the value of
```
Δ2 * ϕ(Δ1 / Δ2)
```
where `0 ≤ ϕ(r) ≤ 2` is the slope limiter
"""
abstract type AbstractSlopeLimiter end

"""
    width(recon::AbstractReconstruction)

Returns the width of the stencil need for the FV reconstruction `recon`. Total
number of values used in the reconstruction are `2width(recon) + 1`
"""
width(recon::AbstractReconstruction) = throw(MethodError(width, (recon,)))

"""
    FVConstant <: AbstractReconstruction

Reconstruction type for cell centered finite volume methods (e.g., constants)
"""
struct FVConstant <: AbstractReconstruction end

width(::FVConstant) = 0

function (::FVConstant)(state_bot, state_top, cell_states::SVector{1}, _)
    @inbounds state_top .= cell_states[1]
    @inbounds state_bot .= cell_states[1]
end
"""
    FVLinear{W = 1} <: AbstractReconstruction

Reconstruction type for limited linear reconstruction finite volume methods.

!!! note
   The optional type parameter `W` is mainly for debugging purposes and allows
   the stencil to be artificially widened to make sure the kernels work with
   wide stencils.

    FVLinear(limiter = VanLeer())
    FVLinear{W}(limiter = VanLeer())

Construct the `FVLinear` reconstruction type with the given slope `limiter`
"""
struct FVLinear{W, L} <: AbstractReconstruction
    limiter::L

end

"""
    FVLinear(limiter = VanLeer())
    FVLinear{W}(limiter = VanLeer())

Construct the `FVLinear` reconstruction type with the given slope `limiter` and
optional width `W`.
"""
function FVLinear{W}(limiter = VanLeer()) where {W}
    @assert W > 0
    FVLinear{W, typeof(limiter)}(limiter)
end
FVLinear(limiter = VanLeer()) = FVLinear{1}(limiter)

width(::FVLinear{W}) where {W} = W

function (fvrecon::FVLinear)(
    state_bot,
    state_top,
    cell_states::SVector{3},
    cell_weights,
)
    @inbounds wi_top = 1 / (cell_weights[3] + cell_weights[2])
    @inbounds wi_bot = 1 / (cell_weights[2] + cell_weights[1])
    @inbounds @unroll for s in 1:length(state_top)
        # Compute the edge gradient approximations
        Δ_top = wi_top * (cell_states[3][s] - cell_states[2][s])
        Δ_bot = wi_bot * (cell_states[2][s] - cell_states[1][s])

        # Compute the limited slope
        Δ = fvrecon.limiter(Δ_top, Δ_bot)

        # Compute the reconstructions at the cell edges
        state_top[s] = cell_states[2][s] + Δ * cell_weights[2]
        state_bot[s] = cell_states[2][s] - Δ * cell_weights[2]
    end
end

function (fvrecon::FVLinear)(
    state_bot,
    state_top,
    cell_states::SVector{1},
    cell_weights,
)
    FVConstant()(state_bot, state_top, cell_states, cell_weights)
end

function (fvrecon::FVLinear)(
    state_bot,
    state_top,
    cell_states::SVector{D},
    cell_weights,
) where {D}
    W = div(D - 1, 2)
    rng = SUnitRange(W, W + 2)
    @inbounds fvrecon(state_bot, state_top, cell_states[rng], cell_weights[rng])
end

"""
    VanLeer <: AbstractSlopeLimiter

Classic Van-Leer limiter
```
   ϕ(r) = (r + |r|) / (1 + r)
```

### References
    @article{van1974towards,
      title = {Towards the ultimate conservative difference scheme. II.
               Monotonicity and conservation combined in a second-order scheme},
      author = {Van Leer, Bram},
      journal = {Journal of computational physics},
      volume = {14},
      number = {4},
      pages = {361--370},
      year = {1974},
      doi = {10.1016/0021-9991(74)90019-9}
    }
"""
struct VanLeer <: AbstractSlopeLimiter end

function (::VanLeer)(Δ_top, Δ_bot)
    FT = eltype(Δ_top)
    if Δ_top * Δ_bot > 0
        return 2Δ_top * Δ_bot / (Δ_top + Δ_bot)
    else
        return FT(0)
    end
end

"""
    NoLimiter <: AbstractSlopeLimiter

No Limiter:
1) smooth case, no dissipation is added
2) linear case, which guarantee the discreted equation is linear
```
   ϕ(r) = r
```
"""
struct NoLimiter <: AbstractSlopeLimiter end

function (::NoLimiter)(Δ_top, Δ_bot)
    return (Δ_top + Δ_bot) / 2
end

end


text:
#MDDatasets: Support broadcast of DataHR types
#-------------------------------------------------------------------------------

#TODO: centralize findytypes/promote_type with eltype/result_type functions

#==Type definitions
===============================================================================#

mutable struct SubscriptMap
	outidx::Vector{Int} #List of indices (of output array)
	outlen::Int #Number of indices in output subscript
end
SubscriptMap(inlen::Int, outlen::Int) =
	SubscriptMap(Vector{Int}(undef, inlen), outlen)

#==Error generators
===============================================================================#

function error_mismatchedsweep(basesweep::Vector{PSweep}, subsweep::Vector{PSweep})
	msg = "Mismatched sweeps:\n\nSweep1:\n$basesweep\n\nSweep2:\n$subsweep"
	return ArgumentError(msg)
end

#==Helper functions
===============================================================================#

function ensureuniqueids(s::Vector{PSweep})
	n = length(s)
	for i in 1:n
		for j in (i+1):n
			if s[i] == s[j]
				throw(ArgumentError("Sweep id not unique: \"$(s[i])\""))
			end
		end
	end
end

#Find "base" sweep (most complex data configuration to broadcast up to)
#-------------------------------------------------------------------------------
function basesweep(s1::Vector{PSweep}, s2::Vector{PSweep})
	return length(s1)>length(s2) ? s1 : s2
end
basesweep(s::Vector{PSweep}, d::DataHR) = basesweep(s,d.sweeps)
basesweep(s::Vector{PSweep}, d::DataF1) = s
basesweep(s::Vector{PSweep}, d::Number) = s
basesweep(s::Vector{PSweep}, v::Vector{T}) where T<:Number = s
basesweep(d1::DataHR, d2::DataHR) = basesweep(d1.sweeps,d2.sweeps)
basesweep(d1::DataHR, d2) = basesweep(d1.sweeps,d2)
basesweep(d1, d2::DataHR) = basesweep(d2.sweeps,d2)

#Functions to map array dimensions when broadcasting up a DataHR dataset
#-------------------------------------------------------------------------------
function getmap(basesweep::Vector{PSweep}, subsweep::Vector{PSweep})
	ensureuniqueids(basesweep)
	result = SubscriptMap(length(basesweep), length(subsweep))
	found = zeros(Bool, length(subsweep))
	for i in 1:length(basesweep)
		idx = findfirst((x)->(x.id==basesweep[i].id), subsweep)
		if nothing == idx; idx = 0; end
		result.outidx[i] = idx
		if idx>1
			if basesweep[i].v != subsweep[idx].v
				msg = "Mismatched sweeps:\n$basesweep\n$subsweep"
				throw(error_mismatchedsweep(basesweep, subsweep))
			end
			found[idx] = true
		end
	end
	if !all(found); throw(error_mismatchedsweep(basesweep, subsweep)); end
	return result
end
function remap(_map::SubscriptMap, inds::Vector{Int})
	result = Vector{Int}(_map.outlen)
	for i in 1:length(inds)
		idx = _map.outidx[i]
		if idx > 0; result[idx] = inds[i]; end
	end
	return result
end

#==Broadcasting data up-to a given sweep dimension
===============================================================================#
#WARN: Renamed broadcastMDSweep due to strange dispatch issues as "broadcastMD"
function broadcastMDSweep(s::Vector{PSweep}, d::T) where T<:Number
	result = DataHR{T}(s)
	for i in 1:length(result.elem)
		result.elem[i] = d
	end
	return result
end
function broadcastMDSweep(s::Vector{PSweep}, d::DataF1)
	result = DataHR{DataF1}(s)
	for i in 1:length(result.elem)
		result.elem[i] = d
	end
	return result
end
function broadcastMDSweep(s::Vector{PSweep}, d::DataHR{T}) where T
	if s == d.sweeps; return d; end
	_map = getmap(s, d.sweeps)
	result = DataHR{T}(s)
	for inds in subscripts(result)
		result.elem[inds...] = d.elem[remap(_map, inds)...]
	end
	return result
end

#==Broadcast function call on multi-dimensional data
===============================================================================#
#Broadcast data up to base sweep of two first arguments, then call fn
function _broadcast(::Type{T}, s::Vector{PSweep}, fn::Function, args...; kwargs...) where T
	if length(s) < 1
		return fn(args...; kwargs...)
	end
	bargs = Vector{Any}(undef, length(args)) #Broadcasted version of args
	for i in 1:length(args)
		if typeof(args[i])<:DataMD
			bargs[i] = broadcastMDSweep(s, args[i])
		else
			bargs[i] = args[i]
		end
	end
	bkwargs = Vector{Pair}(undef, length(kwargs)) #Broadcasted version of kwargs
	for (i, kv) in enumerate(kwargs)
		(k,v) = kv
		if typeof(v)<:DataMD
			bkwargs[i] = (k => broadcastMDSweep(s, v))
		else
			bkwargs[i] = kv
		end
	end
	result = DataHR{T}(s) #Create empty result
	for i in 1:length(result.elem)
		curargs = Vector{Any}(undef, length(bargs))
		for j in 1:length(bargs)
			if typeof(bargs[j]) <: DataHR
				curargs[j] = bargs[j].elem[i]
			else
				curargs[j] = bargs[j]
			end
		end
		curkwargs = Vector{Pair}(undef, length(bkwargs))
		for j in 1:length(bkwargs)
			(k,v) = bkwargs[j]
			if typeof(v) <: DataHR
				curkwargs[j] = (k => v.elem[i])
			else
				curkwargs[j] = bkwargs[j]
			end
		end
		result.elem[i] = fn(curargs...; curkwargs...)
	end
	return result
end

#Find base sweep for a 1-argument broadcastMD
#-------------------------------------------------------------------------------
fnbasesweep(fn::Function, d) = PSweep[]
fnbasesweep(fn::Function, d::DataHR{T}) where T = d.sweeps

#Ensure collection is composed of DataF1 (ex: DataHR{DataF1}):
#Collapses outer-most dimension of DataHR{Number} to a DataHR{DataF1} value, if necessary
#-------------------------------------------------------------------------------
ensure_coll_DataF1(fn::Function, d) = d #Plain data is ok.
ensure_coll_DataF1(fn::Function, d::DataHR) = DataHR{DataF1}(d)

#Broadcast functions capable of operating directly on 1 base type (Number):
#-------------------------------------------------------------------------------
#DataHR{DataF1/Number}
broadcastMD(::CastType1{Number,1}, fn::Function, d::DataHR{T}, args...; kwargs...) where T =
	_broadcast(T, fnbasesweep(fn, d), fn, d, args...; kwargs...)
#Data reducing (DataHR{DataF1/Number})
function broadcastMD(::CastTypeRed1{Number,1}, fn::Function, d::DataHR{T}, args...; kwargs...) where T<:Number
	d = ensure_coll_DataF1(fn, d) #Collapse DataHR{Number}  => DataHR{DataF1}
	_broadcast(T, fnbasesweep(fn, d), fn, d, args...; kwargs...)
end
function broadcastMD(::CastTypeRed1{Number,1}, fn::Function, d::DataHR{DataF1}, args...; kwargs...)
	TR = promote_type(findytypes(d.elem)...) #TODO: Better way?
	_broadcast(TR, fnbasesweep(fn, d), fn, d, args...; kwargs...)
end

#Broadcast functions capable of operating only on a dataF1 value:
#-------------------------------------------------------------------------------
#TODO: These signatures come in conflict with what is needed for DataRS
#DataF1
function broadcastMD(::CastType1{DataF1,1}, fn::Function, d, args...; kwargs...)
	d = ensure_coll_DataF1(fn, d) #Collapse DataHR{Number}  => DataHR{DataF1}
	_broadcast(DataF1, fnbasesweep(fn, d), fn, d, args...; kwargs...)
end
#Expects DataF1 @ arg #2:
function broadcastMD(::CastType1{DataF1,2}, fn::Function, dany1, d, args...; kwargs...)
	d = ensure_coll_DataF1(fn, d) #Collapse DataHR{Number}  => DataHR{DataF1}
	_broadcast(DataF1, fnbasesweep(fn, d), fn, dany1, d, args...; kwargs...)
end
#Data reducing (DataF1)
function broadcastMD(::CastTypeRed1{DataF1,1}, fn::Function, d, args...; kwargs...)
	d = ensure_coll_DataF1(fn, d) #Collapse DataHR{Number}  => DataHR{DataF1}
	TR = promote_type(findytypes(d.elem)...) #TODO: Better way?
	_broadcast(TR, fnbasesweep(fn, d), fn, d, args...; kwargs...)
end

#Find base sweep for a 2-argument broadcast
#-------------------------------------------------------------------------------
function fnbasesweep(fn::Function, d1, d2)
	local s
	try
		s = basesweep(d1,d2)
	catch
		t1 = typeof(d1); t2 = typeof(d2)
		msg = "No signature found for $fn($t1, $t2, ...)"
		throw(ArgumentError(msg))
	end
end

#Ensure collection is composed of DataF1 (ex: DataHR{DataF1}):
#Collapses outer-most dimension of DataHR{Number} to a DataHR{DataF1} value, if necessary
#-------------------------------------------------------------------------------
function ensure_coll_DataF1(fn::Function, d1, d2)
	try
		d1 = ensure_coll_DataF1(fn, d1)
		d2 = ensure_coll_DataF1(fn, d2)
	catch
		t1 = typeof(d1); t2 = typeof(d2)
		msg = "No signature found for $fn($t1, $t2, ...)"
		throw(ArgumentError(msg))
	end
	return tuple(d1, d2)
end

#Broadcast functions capable of operating directly on base types (Number, Number):
#-------------------------------------------------------------------------------
#DataHR{DataF1/Number} & DataHR{DataF1/Number}:
function broadcastMD(::CastType2{Number,1,Number,2}, fn::Function,
	d1::DataHR{T1}, d2::DataHR{T2}, args...; kwargs...) where {T1,T2}
	_broadcast(promote_type(T1,T2), fnbasesweep(fn, d1, d2), fn, d1, d2, args...; kwargs...)
end
#DataHR{DataF1/Number} & DataF1/Number:
function broadcastMD(::CastType2{Number,1,Number,2}, fn::Function,
	d1::DataHR{T1}, d2::T2, args...; kwargs...) where {T1,T2<:DF1_Num}
	_broadcast(promote_type(T1,T2), fnbasesweep(fn, d1, d2), fn, d1, d2, args...; kwargs...)
end
#DataF1/Number & DataHR{DataF1/Number}:
function broadcastMD(::CastType2{Number,1,Number,2}, fn::Function,
	d1::T1, d2::DataHR{T2}, args...; kwargs...) where {T1<:DF1_Num,T2}
	_broadcast(promote_type(T1,T2), fnbasesweep(fn, d1, d2), fn, d1, d2, args...; kwargs...)
end

#Broadcast functions capable of operating only on a dataF1 value:
#-------------------------------------------------------------------------------
#DataF1, DataF1
function broadcastMD(::CastType2{DataF1,1,DataF1,2}, fn::Function, d1, d2, args...; kwargs...)
	(d1, d2) = ensure_coll_DataF1(fn, d1, d2) #Collapse DataHR{Number}  => DataHR{DataF1}
	_broadcast(DataF1, fnbasesweep(fn, d1, d2), fn, d1, d2, args...; kwargs...)
end
#DataF1, DataF1 @ arg 2/3:
function broadcastMD(::CastType2{DataF1,2,DataF1,3}, fn::Function, dany1, d1, d2, args...; kwargs...)
	(d1, d2) = ensure_coll_DataF1(fn, d1, d2) #Collapse DataHR{Number}  => DataHR{DataF1}
	_broadcast(DataF1, fnbasesweep(fn, d1, d2), fn, dany1, d1, d2, args...; kwargs...)
end
#Data reducing (DataF1, DataF1)
function broadcastMD(::CastTypeRed2{DataF1,1,DataF1,2}, fn::Function, d1, d2, args...; kwargs...)
	(d1, d2) = ensure_coll_DataF1(fn, d1, d2) #Collapse DataHR{Number}  => DataHR{DataF1}
	TR = promote_type(findytypes(d1.elem)...,findytypes(d2.elem)...) #TODO: Better way?
	_broadcast(DataF1, fnbasesweep(fn, d1, d2), fn, d1, d2, args...; kwargs...)
end

#More custom broadcastMD functions
#-------------------------------------------------------------------------------
function broadcastMD(::CastType2{DataF1,1,Number,2}, fn::Function, d1, d2, args...; kwargs...)
	d1 = ensure_coll_DataF1(fn, d1) #Collapse DataHR{Number}  => DataHR{DataF1}
	TR = promote_type(findytypes(d1.elem)...,eltype(d2)) #TODO: Better way?
	_broadcast(DataF1, fnbasesweep(fn, d1, d2), fn, d1, d2, args...; kwargs...)
end

#Last Line


text:
include("../covariatefun.jl")
using Optim

"""
cdf(u, s) is a function of location s and value u, gradients and hessian are
w.r.t the vector (u, s).

INPUTS:
cdf - takes in (x, Fx, y) and outputs cdf
cdf_gradient - takes in (x, Fx, y) and outputs Jacobian
cdf_hessian - takes in (x, Fx, y) and outputs Hessian

OUTPUTS:
"""
function optimizeUCB(cdf, cdf_gradient, cdf_hessian, lx, ux)
    #v is augmented vector consisting of [u, s], where u is value and s is location
    cdf_fixed(v) = cdf(v[2:end], linear_polynomial_basis(v[2:end]), v[1])
    cdf_gradient_fixed(v) = cdf_gradient(v[2:end], linear_polynomial_basis(v[2:end]), v[1])
    cdf_hessian_fixed(v) = cdf_hessian(v[2:end], linear_polynomial_basis(v[2:end]), v[1])
    d = length(lx)  #location is d-dimensional, while label y is 1-dimensional, so augmented vector is (d+1)-dimensional
    fun(x) = x[1] #x\in R^d+1
    function fun_grad!(g, x)
        g[1] = 1
        for i = 2:d
            g[i] = 0
        end
    end;
    function fun_hess!(h, x)
        for i = 1:d
            for j = 1:d
                h[i, j] = 0
            end
        end
    end;
    con_c!(c, x) = (c[1] = cdf_fixed(x); c)
    function con_jacobian!(J::Array{T} where T <:Float64, x)
        grad = cdf_gradient_fixed(x)
        for i = 1:d
            J[i] = grad[i]
        end
        J
    end;
    function con_h!(h, x, λ)
        hess = cdf_hessian_fixed(x)
        for i = 1:d
            for j = 1:d
                h[i, j] += λ[1]*hess[i, j]
            end
        end
    end;
    #lx = [-5.0, -5.0]; ux = [5.0, 5.0] #function box constraint
    lc = [0.24]; uc = [0.26] #approximation to equality constraintx
    x0 = init_constrained_pt(cdf_fixed, lx[2:end], ux[2:end]; quantile = 0.25)
    @info "x0 found with init_constr_point", x0
    @info "cdf_fixed(x0)", cdf_fixed(x0)
    df = TwiceDifferentiable(fun, fun_grad!, fun_hess!, x0)
    @info "lx" lx
    @info "ux" ux
    @info "lc" lc
    @info "uc" uc
    dfc = TwiceDifferentiableConstraints(con_c!, con_jacobian!, con_h!, lx, ux, lc, uc)
    @info x0
    Optim.Options(allow_f_increases = true, successive_f_tol = 2)
    res = optimize(df, dfc, x0, IPNewton())
    minimizer = res.minimizer
    @info "does constraint hold?", cdf_fixed(minimizer)
    return (minimizer, res.minimum, res)
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
using Requests
import Requests: get, post, put, delete, options

function safe(name)
	result = replace(name, ":", " ")
	result = replace(result, "\\", " ")
	result = replace(result, "/", " ")
	result = replace(result, "*", " ")
	result = replace(result, "?", " ")
	result = replace(result, "|", " ")
	result = replace(result, "<", " ")
	result = replace(result, ">", " ")
	result = replace(result, "\"", " ")
	result = replace(result, r"[^\x20-\x7e]", "")
	result = strip(result)
	return result
end

function urlsafe(title)
	result = replace(title, r"[^\x20-\x7e]", "")
	result = strip(result)
	result = replace(result, " ", "+")
	result = replace(result, ",", "")
	result = replace(result, "!", "")
	result = replace(result, "'", "")
	result = replace(result, "-", "")
	return result
end

function main(dataPath)
	lines = parse(readall(pipeline(`cat $[dataPath]movies.csv`, `wc -l`)))
	@time movies = readdlm(dataPath * "movies.csv", ',', header=true, dims=(lines, 3))
	movies = movies[1]

	@printf "Finished loading, going into loop\n"

	count = 0
	for i = 1:length(movies[:, 1])
		title = movies[i, 2]

		if !isfile("movies/$(safe(title))/1.jpg") ||
			!isfile("movies/$(safe(title))/2.jpg")	||
			!isfile("movies/$(safe(title))/3.jpg")
			count += 1
		end
	end

	@show count

	@printf "Done :)\n"
end

main("../ml-latest/fixed/")


text:
function testNHInertiaCart()
close all

% Physical parameters
Param=PhysParameters();

% Grid
nx=60;
ny=2;
Param.Lx=300*1.e3;
Param.Ly=2*Param.Lx/nx;
x0=0;
y0=0;
Param.H=10*1.e3;
Param.OrdPoly=4;
%Horizontal grid size
dx=Param.Lx/nx/(Param.OrdPoly+1);
dz=min(dx/2,Param.H/10);
nz=ceil(Param.H/dz);
nz=40;

Boundary.WE='Period';
Boundary.BT='Period';
Param.hS='';
Param.Grid=CartGrid(nx,ny,Param.Lx,Param.Ly,x0,y0,OrientFaceCart,Boundary,Param);

Param.Grid.nz=nz;
Param.Grid.zP=zeros(nz,1);
Param.Grid.z=zeros(nz+1,1);
Param.Grid.dz=Param.H/nz;
Param.Grid.zP(1)=Param.Grid.dz/2;
for i=2:nz
  Param.Grid.zP(i)=Param.Grid.zP(i-1)+Param.Grid.dz;
end
for i=2:nz+1
  Param.Grid.z(i)=Param.Grid.z(i-1)+Param.Grid.dz;
end

% Model
Param.ProfRho='InertiaGravityWave';
Param.ProfTheta='InertiaGravityWave';
Param.ProfVel='Const';
Param.Damping=false;
Param.Coriolis=false;
Param.Buoyancy=true;
Param.Source=false;
Param.xC=Param.Lx/3;
Param.d=5000;
Param.a=5*1.e3;
Param.Th0=300;
Param.DeltaTh=1.e-2;
Param.uMax=0;
Param.vMax=0;
Param.NBr=1.e-2;
Param.Equation='Compressible';
Param.TopoS='';
Param.NumV=5;
Param.RhoPos=1;
Param.uPos=2;
Param.vPos=3;
Param.wPos=4;
Param.ThPos=5;
Param.Thermo='';

% Discretization
OrdPolyZ=1;
[CG,Param]=Discretization(Param.OrdPoly,OrdPolyZ,@JacobiDG3,Param);
Param.HyperVisc=true;
Param.HyperDCurl=2.e5; %1.e14*(dx/LRef)^3.2;
Param.HyperDGrad=2.e5;
Param.HyperDDiv=2.2e5;
Param.Upwind=false;

% Output
Param.Flat=false;
Param.vtkFileName='InertiaGravityWave';
Param.vtk=0;
vtkGrid=vtkCGGrid(CG,@TransCart,@Topo,Param);
Param.cNames(1).s='Rho';
Param.cNames(2).s='u';
Param.cNames(3).s='w';
Param.cNames(4).s='Th';
Param.cNames(5).s='ThetaPrime';
Param.cNames(6).s='Pres';

% Initial conditions
U=zeros(CG.NumG,nz,Param.NumV);
U(:,:,Param.RhoPos)=Project(@fRho,CG,Param);
[U(:,:,Param.uPos),U(:,:,Param.vPos)]=ProjectVec(@fVel,CG,Param);
U(:,:,Param.ThPos)=Project(@fTheta,CG,Param).*U(:,:,Param.RhoPos);
Param.ThetaBGrd=Project(@fThetaBGrd,CG,Param);
pBGrd=Project(@fpBGrd,CG,Param);
RhoBGrd=Project(@fRhoBGrd,CG,Param);
OP=Param.OrdPoly+1;
NF=Param.Grid.NumFaces;
RhoBGrd=reshape(RhoBGrd(reshape(CG.Glob,OP*OP*NF,1),:)...
  ,OP,OP,NF,nz);
Param.RhoBGrdF=0.5*(RhoBGrd(:,:,:,1:nz-1)+RhoBGrd(:,:,:,2:nz));
Param.pBGrd=reshape(pBGrd(reshape(CG.Glob,OP*OP*NF,1),:)...
  ,OP,OP,NF,nz);

% Integration
dtau=.40;
time=0;
%IntMethod='Rosenbrock'; %'RungeKutta';
IntMethod='RungeKutta';
Param.EndTime=1800;
nIter=Param.EndTime/dtau;
PrintTime=100;
PrintInt=PrintTime/dtau;
Param.RK=RungeKuttaMethod('RK4');
Param.ROS=RosenbrockMethod('ROSRK3');
Param.vtk=vtkOutput(U,vtkGrid,CG,Param);
switch IntMethod
  case 'Rosenbrock'
    for i=1:nIter
      U=Rosenbrock(U,dtau,@FcnNHCurlVec,@Jac,CG,Param);
      time=time+dtau;
      if mod(i,PrintInt)==0
        Param.vtk=vtkOutput(U,vtkGrid,CG,Param);
      end
    end
  case 'RungeKutta'
    for i=1:nIter
      i
      U=RungeKuttaExplicit(U,dtau,@FcnNHCurlVec,CG,Param);
      time=time+dtau;
      if mod(i,PrintInt)==0
        Param.vtk=vtkOutput(U,vtkGrid,CG,Param);
      end
    end
end

end


text:
"""
    struct Domain

Represent a physical domain as a union of entities.

# See also: [`AbstractEntity`](@ref), [`ElementaryEntity`](@ref).
"""
struct Domain
    entities::Vector{AbstractEntity}
end

"""
    entities(Ω::Domain)

Return a vector of all elementary entities making up a domain.
"""
entities(Ω::Domain) = Ω.entities

Domain(ω::AbstractEntity) = Domain([ω,])
Domain() = Domain(AbstractEntity[])

function Base.show(io::IO,d::Domain)
    ents = entities(d)
    n = length(entities(d))
    n == 1 ? print(io,"Domain with $n entity:\n") : print(io,"Domain with $n entities:")
    for ent in ents
        print(io,"\n\t $(ent)")
    end
    return io
end

function geometric_dimension(Ω::Domain)
    l,u = extrema(geometric_dimension(ent) for ent in entities(Ω))
    @assert l == u "geometric dimension of entities in a domain not equal"
    return u
end

"""
    length(Ω:::Domain)

The length of a domain corresponds to the number of elementary entities that
make it.
"""
Base.length(Ω::Domain) = length(entities(Ω))

"""
    skeleton(Ω::Domain)

Return all the boundaries of the domain, i.e. the domain's skeleton.
"""
function skeleton(Ω::Domain)
    ents = AbstractEntity[]
    for ent in entities(Ω)
        append!(ents,boundary(ent))
    end
    Domain(unique!(ents))
end

"""
    ===(Ω1::Domain,Ω2::Domain)

Two `Domain`s are equal if all their entities are equal (regardless of order).
"""
function Base.:(==)(Ω1::Domain, Ω2::Domain)
    return issetequal(entities(Ω1),entities(Ω2))
end

"""
    in(ω::ElementaryEntity,Ω::Domain)

Check whether an `ElementaryEntity` belongs to a `Domain` by recursively
checking whether it belongs to its boundary.
"""
function Base.in(ω::ElementaryEntity, Ω::Domain)
    ents = entities(Ω)
    in(ω, ents) && return true
    # TODO: should we really recurse on the boundary of the entities composing
    # the domain for determining if an entity is in a domain.
    for ent in ents
        in(ω,Domain(boundary(ent))) && (return true)
    end
    return false
end

Base.getindex(Ω::Domain, i) = entities(Ω)[i]

Base.lastindex(Ω::Domain) = length(Ω)

"""
    iterate(Ω::Domain)

Iterating over a domain means iterating over its entities.
"""
function Base.iterate(Ω::Domain, state=1)
    # Check if we are done iterating
    if state > length(Ω)
        return nothing
    end
    # Return (result, state)
    return (Ω[state], state+1)
end

Base.isempty(Ω::Domain) = length(entities(Ω)) == 0

function Base.setdiff(Ω1::Domain, Ω2::Domain)
    Domain(setdiff(entities(Ω1), entities(Ω2)))
end

function Base.union(Ω1::Domain,Ωs::Domain...)
    ents = vcat(entities(Ω1),entities.(Ωs)...)
    Domain(unique!(ents))
end
Base.union(Ω::Domain) = Domain(unique(entities(Ω)))

"""
    assertequaldim(Ω1::Domain,Ω2::Domain)

Check that two domains have same dimension.

If one of the domain (or both) are empty, the assertion is assumed to be true.
"""
function assertequaldim(Ω1::Domain, Ω2::Domain)
    if isempty(Ω1) || isempty(Ω2)
        return true
    else
        msg = "The dimension of the first domain should be equal to the dimension
        of the second domain."
        @assert geometric_dimension(Ω1) == geometric_dimension(Ω2) msg
    end
end

function Base.intersect(Ω1::Domain, Ω2::Domain)
    assertequaldim(Ω1, Ω2)
    Ωinter = Domain(intersect(entities(Ω1), entities(Ω2)))
    if isempty(Ωinter)
        if isempty(Ω1) || isempty(Ω2)
            return Domain()
        else
            return intersect(skeleton(Ω1), skeleton(Ω2))
        end
    else
        return Ωinter
    end
end

function Base.push!(Ω::Domain,ent::ElementaryEntity)
    push!(entities(Ω),ent)
end

function Base.issubset(Ω1::Domain, Ω2::Domain)
    assertequaldim(Ω1, Ω2)
    return issubset(entities(Ω1), entities(Ω2))
end

"""
    boundary(Ω)

Return a domain comprising the external boundary of Ω.

See also: [`external_boundary`](@ref)
"""
boundary(Ω::Domain) = external_boundary(Ω::Domain)

"""Return the internal boundaries inside a domain."""
function internal_boundary(Ω::Domain)
    Ω1 = Domain(Ω[1])
    γ  = Domain()
    for ω2 in Ω[2:end]
        Ω2 = Domain(ω2)
        γ1 = Domain(vcat(boundary.(entities(Ω1))...))
        γ2 = Domain(vcat(boundary.(entities(Ω2))...))
        γ = union(γ, intersect(γ1, γ2))
        Ω1 = union(Ω1, Ω2)
    end
    return γ
end

"""Return the external boundaries inside a domain."""
function external_boundary(Ω::Domain)
    return setdiff(skeleton(Ω),internal_boundary(Ω))
end

"""
Return all tags of the elementary entities in the domain `Ω` corresponding to the dimension `d`.
"""
function Base.keys(Ω::Domain, d::Integer)
    if isempty(Ω)
        return Tuple{Int64,Int64}[]
    elseif d == geometric_dimension(Ω)
        return Vector{Tuple{Int64,Int64}}(vcat(key.(entities(Ω))))
    elseif d < geometric_dimension(Ω)
        return unique(keys(skeleton(Ω),d))
    else
        error("Asking for tags with dimension > dimension of domain")
    end
end
function Base.keys(Ω::Domain)
    isempty(Ω) ? Tuple{Int64, Int64}[] : keys(Ω, geometric_dimension(Ω))
end

"""
Return all tags of the elementary entities in the domain `Ω` corresponding to the dimensions contained in `dims`.
"""
function Base.keys(Ω::Domain, dims::Vector{T}) where T <: Integer
    tgs = Vector{Tuple{Int64, Int64}}(undef, 0)
    for d in dims push!(tgs, keys(Ω, d)...) end
    return tgs
end

max_stars_repo_name: WaveProp/WavePropBase

max_issues_repo_name: WaveProp/WavePropBase

max_forks_repo_name: WaveProp/WavePropBase


text:
# Similar to JuMP/test/solvers.jl

function try_import(name::Symbol)
    try
        @eval import $name
        return true
    catch e
        return false
    end
end

# Load an available solver
lp_solver = nothing
glp = try_import(:GLPKMathProgInterface)
glp && (lp_solver = GLPKMathProgInterface.GLPKSolverLP())
if lp_solver === nothing
    cbc = try_import(:Cbc)
    if cbc; import Clp; end
    cbc && (lp_solver = Clp.ClpSolver())
end
if lp_solver === nothing
    grb = try_import(:Gurobi) # Gurobi creates BigFloat when the input is BigInt and then cannot handle it
    grb && (lp_solver = Gurobi.GurobiSolver(OutputFlag=0))
end
if lp_solver === nothing
    cpx = try_import(:CPLEX)
end
if lp_solver === nothing
    xpr = try_import(:Xpress)
    xpr && (lp_solver = Xpress.XpressSolver(OUTPUTLOG=0))
end
if lp_solver === nothing
    mos = try_import(:Mosek)
    mos && (lp_solver = Mosek.MosekSolver(LOG=0))
end
if lp_solver === nothing
    ipt = try_import(:Ipopt)
    ipt && (lp_solver = Ipopt.IpoptSolver(print_level=0))
end
if lp_solver === nothing
    eco = try_import(:ECOS)
    eco && (lp_solver = ECOS.ECOSSolver(verbose=false))
end
if lp_solver === nothing
    scs = try_import(:SCS)
    scs && (lp_solver = SCS.SCSSolver(eps=1e-6,verbose=0))
    @assert lp_solver !== nothing
end


text:
using StaticArrays, Test
using OrdinaryDiffEq

u0 = fill(zero(MVector{2,Float64}), 2)
u0[1] = ones(MVector{2,Float64}) .+ 1
f = (du,u,p,t) -> du .= u
ode = ODEProblem(f, u0, (0.,1.))
sol = solve(ode, Euler(), dt=1e-2)
sol = solve(ode, Tsit5())

u0 = fill(zero(SVector{2,Float64}), 2) .+ 1
u0[1] = ones(SVector{2,Float64}) .+ 1
ode = ODEProblem(f, u0, (0.,1.))
sol = solve(ode, Euler(), dt=1e-2)
sol = solve(ode, Tsit5())

sol = solve(ode, SSPRK22(), dt=1e-2)

u0 = ones(MVector{2,Float64})
ode = ODEProblem(f, u0, (0.,1.))
sol = solve(ode, Euler(), dt=1e-2)
sol = solve(ode, Tsit5(), dt=1e-2)

u0 = ones(SVector{2,Float64})
f = (u,p,t) -> u
ode = ODEProblem(f, u0, (0.,1.))
sol = solve(ode, Euler(), dt=1e-2)
integrator = init(ode, ImplicitEuler())
@test OrdinaryDiffEq.get_W(integrator.cache.nlsolver) isa StaticArrays.LU
sol = solve(ode, ImplicitEuler())
sol = solve(ode, ImplicitEuler(nlsolve=NLAnderson()))
sol = solve(ode, Tsit5(), dt=1e-2)


text:
# Generative model w/ Glow architecture from Kingma & Dhariwal (2018)
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: January 2020

using InvertibleNetworks, LinearAlgebra, Test, Random

# Define network
nx = 32
ny = 32
n_in = 2
n_hidden = 4
batchsize = 2
L = 2
K = 2

###################################################################################################
# Invertibility

# Network and input
G = NetworkGlow(n_in, n_hidden, L, K)
X = rand(Float32, nx, ny, n_in, batchsize)

Y = G.forward(X)[1]
X_ = G.backward(Y, Y)[2]

@test isapprox(norm(X - X_)/norm(X), 0f0; atol=1f-5)

###################################################################################################
# Test gradients are set and cleared

P = get_params(G)
gsum = 0
for p in P
    ~isnothing(p.grad) && (global gsum += 1)
end
@test isequal(gsum, L*K*10)

clear_grad!(G)
gsum = 0
for p in P
    ~isnothing(p.grad) && (global gsum += 1)
end
@test isequal(gsum, 0)

###################################################################################################
# Gradient test

function loss(G, X)
    Y, logdet = G.forward(X)
    f = -log_likelihood(Y) - logdet
    ΔY = -∇log_likelihood(Y)
    ΔX, X_ = G.backward(ΔY, Y)
    return f, ΔX, G.CL[1,1].RB.W1.grad, G.CL[1,1].C.v1.grad
end

# Gradient test w.r.t. input
G = NetworkGlow(n_in, n_hidden, L, K)
X = rand(Float32, nx, ny, n_in, batchsize)
X0 = rand(Float32, nx, ny, n_in, batchsize)
dX = X - X0

f0, ΔX = loss(G, X0)[1:2]
h = 0.1f0
maxiter = 4
err1 = zeros(Float32, maxiter)
err2 = zeros(Float32, maxiter)

print("\nGradient test glow: input\n")
for j=1:maxiter
    f = loss(G, X0 + h*dX,)[1]
    err1[j] = abs(f - f0)
    err2[j] = abs(f - f0 - h*dot(dX, ΔX))
    print(err1[j], "; ", err2[j], "\n")
    global h = h/2f0
end

@test isapprox(err1[end] / (err1[1]/2^(maxiter-1)), 1f0; atol=1f1)
@test isapprox(err2[end] / (err2[1]/4^(maxiter-1)), 1f0; atol=1f1)

# Gradient test w.r.t. parameters
X = rand(Float32, nx, ny, n_in, batchsize)
G = NetworkGlow(n_in, n_hidden, L, K)
G0 = NetworkGlow(n_in, n_hidden, L, K)
Gini = deepcopy(G0)

# Test one parameter from residual block and 1x1 conv
dW = G.CL[1,1].RB.W1.data - G0.CL[1,1].RB.W1.data
dv = G.CL[1,1].C.v1.data - G0.CL[1,1].C.v1.data

f0, ΔX, ΔW, Δv = loss(G0, X)
h = 0.1f0
maxiter = 4
err3 = zeros(Float32, maxiter)
err4 = zeros(Float32, maxiter)

print("\nGradient test glow: input\n")
for j=1:maxiter
    G0.CL[1,1].RB.W1.data = Gini.CL[1,1].RB.W1.data + h*dW
    G0.CL[1,1].C.v1.data = Gini.CL[1,1].C.v1.data + h*dv

    f = loss(G0, X)[1]
    err3[j] = abs(f - f0)
    err4[j] = abs(f - f0 - h*dot(dW, ΔW) - h*dot(dv, Δv))
    print(err3[j], "; ", err4[j], "\n")
    global h = h/2f0
end

@test isapprox(err3[end] / (err3[1]/2^(maxiter-1)), 1f0; atol=1f1)
@test isapprox(err4[end] / (err4[1]/4^(maxiter-1)), 1f0; atol=1f1)

###################################################################################################
# Jacobian-related tests

# Gradient test

# Initialization
G = NetworkGlow(n_in, n_hidden, L, K); G.forward(randn(Float32, nx, ny, n_in, batchsize))
θ = deepcopy(get_params(G))
G0 = NetworkGlow(n_in, n_hidden, L, K); G0.forward(randn(Float32, nx, ny, n_in, batchsize))
θ0 = deepcopy(get_params(G0))
X = randn(Float32, nx, ny, n_in, batchsize)

# Perturbation (normalized)
dθ = θ-θ0; dθ .*= norm.(θ0)./(norm.(dθ).+1f-10)
dX = randn(Float32, nx, ny, n_in, batchsize); dX *= norm(X)/norm(dX)

# Jacobian eval
dY, Y, _, _ = G.jacobian(dX, dθ, X)

# Test
print("\nJacobian test\n")
h = 0.1f0
maxiter = 5
err5 = zeros(Float32, maxiter)
err6 = zeros(Float32, maxiter)
for j=1:maxiter
    set_params!(G, θ+h*dθ)
    Y_loc, _ = G.forward(X+h*dX)
    err5[j] = norm(Y_loc - Y)
    err6[j] = norm(Y_loc - Y - h*dY)
    print(err5[j], "; ", err6[j], "\n")
    global h = h/2f0
end

@test isapprox(err5[end] / (err5[1]/2^(maxiter-1)), 1f0; atol=1f1)
@test isapprox(err6[end] / (err6[1]/4^(maxiter-1)), 1f0; atol=1f1)

# Adjoint test

set_params!(G, θ)
dY, Y, _, _ = G.jacobian(dX, dθ, X)
dY_ = randn(Float32, size(dY))
dX_, dθ_, _, _ = G.adjointJacobian(dY_, Y)
a = dot(dY, dY_)
b = dot(dX, dX_)+dot(dθ, dθ_)
@test isapprox(a, b; rtol=1f-3)


text:
# Univariate Normal setup, following Richardson & Green (1997).
# This version uses the general version of the sampler.
module NormalNonoptimized

module NormalModel
export Theta, Data, log_likelihood, log_prior, prior_sample!, new_theta, Theta_clear!, Theta_adjoin!, Theta_remove!,
       Hyperparameters, construct_hyperparameters, update_hyperparameters!, update_parameter!

include("Random.jl")
using .Random

const Data = Float64

type Normal_params
    mu::Float64
    sigma::Float64
    n::Int64        # number of data points assigned to this cluster
    sum_x::Float64  # sum of the data points x assigned to this cluster
    sum_xx::Float64 # sum of x^2 for the data points assigned to this cluster
    Normal_params(mu,sigma) = (p=new(); p.mu=mu; p.sigma=sigma; p.n=0; p.sum_x=0.0; p.sum_xx=0.0; p)
end

Normal_logpdf(x,mu,sigma) = (r=(x-mu)/sigma;  -0.5*r*r - 0.5*log(2.0*pi) - log(sigma))
Normal_clear!(p) = (p.sum_x = 0.; p.sum_xx = 0.; p.n = 0)
Normal_adjoin!(p,x) = (p.sum_x += x; p.sum_xx += x*x; p.n += 1)
Normal_remove!(p,x) = (p.sum_x -= x; p.sum_xx -= x*x; p.n -= 1)

type Hyperparameters
    m::Float64  # prior mean of mu
    s::Float64  # prior stddev of mu
    a::Float64  # prior shape of sigma
    b::Float64  # prior rate of sigma (starting value)
    g::Float64  # (hyper)prior shape of b
    h::Float64  # (hyper)prior rate of b
end

Gamma_logpdf(x,a,b) = (a-1)*log(x) - b*x + a*log(b) - lgamma(a)
# log density of SqrtInvGamma (distn of sigma when 1/sigma^2 ~ Gamma(a,b))
SqrtInvGamma_logpdf(x,a,b) = Gamma_logpdf(1/x^2,a,b) + log(2/x^3)

const Theta = Normal_params
log_likelihood(x,p) = Normal_logpdf(x,p.mu,p.sigma)
# prior: Normal(mu|H.m,H.s) * SqrtInvGamma(sigma|H.a,H.b)
log_prior(p,m,s,a,b) = Normal_logpdf(p.mu,m,s) + SqrtInvGamma_logpdf(p.sigma,a,b)
log_prior(p,H) = log_prior(p,H.m,H.s,H.a,H.b)
prior_sample!(p,H) = (p.mu = randn()*H.s + H.m; p.sigma = 1/sqrt(Random.gamma(H.a,H.b)))
new_theta(H) = Normal_params(0.0,1.0)
Theta_clear!,Theta_adjoin!,Theta_remove! = Normal_clear!,Normal_adjoin!,Normal_remove!

function construct_hyperparameters(options)
    x = options.x
    # Use values in Green & Richardson (1997) to enable comparison.
    m = (minimum(x) + maximum(x))/2
    R = maximum(x) - minimum(x)
    s = R
    a = 2.0
    b = 1.0  # b will be updated in Gibbs sampling
    g = 0.2
    h = 10/R^2
    return Hyperparameters(m,s,a,b,g,h)
end

function update_parameter!(theta_a,theta_b,H,active,density)
    n,sum_x,sum_xx = theta_b.n,theta_b.sum_x,theta_b.sum_xx
    # update mean
    ell = 1/H.s^2
    lambda = 1/(theta_a.sigma^2)
    L = lambda*n + ell
    M = (lambda*sum_x + ell*H.m) / L
    if active; theta_b.mu = randn()/sqrt(L) + M; end

    # update standard deviation
    alpha = H.a + 0.5*n
    beta = H.b + 0.5*(sum_xx - 2*theta_b.mu*sum_x + n*theta_b.mu^2)
    if active; theta_b.sigma = 1/sqrt(Random.gamma(alpha,beta)); end

    return (density? log_prior(theta_b,M,1/sqrt(L),alpha,beta) : NaN)
end

function update_hyperparameters!(H,theta,list,t,x,z)
    sum_lambdas = 0.0
    for k = 1:t; sum_lambdas += 1/theta[list[k]].sigma^2; end
    alpha = H.g + H.a*t
    beta = H.h + sum_lambdas
    H.b = Random.gamma(alpha,beta)
end

end # module NormalModel
using .NormalModel

# Include generic code
include("generic.jl")

# Include core sampler code
include("coreNonconjugate.jl")

end # module NormalNonoptimized


text:
module Dyestuff2Test

println()
@info "Dyestuff2 example"

using GLMCopula, RDatasets, Test

# Dataframe with columns: Batch (Categorical), Yield (Int32)
dyestuff2 = dataset("lme4", "Dyestuff2")
groups = unique(dyestuff2[!, :Batch])
n, p, m = length(groups), 1, 1
gcs = Vector{GaussianCopulaVCObs{Float64}}(undef, n)
for (i, grp) in enumerate(groups)
    gidx = dyestuff2[!, :Batch] .== grp
    ni = count(gidx)
    y = Float64.(dyestuff2[gidx, :Yield])
    X = ones(ni, 1)
    V = [ones(ni, ni)]
    gcs[i] = GaussianCopulaVCObs(y, X, V)
end
gcm = GaussianCopulaVCModel(gcs);

# initialize β and τ from least square solution
@info "Initial point:"
init_β!(gcm)
@show gcm.β
# update σ2 from β and τ using the MM algorithm
fill!(gcm.Σ, 1)
update_Σ!(gcm)
@show gcm.τ
@show gcm.Σ;
@test loglikelihood!(gcm, true, false) ≈ -81.43652408
@show gcm.∇β
@show gcm.∇τ
@show gcm.∇Σ

# fit model using NLP on profile loglikelihood
@info "MLE:"
@time GLMCopula.fit!(gcm, IpoptSolver(print_level=0))
@show gcm.β
@show gcm.τ
@show gcm.Σ
@test loglikelihood!(gcm, true, false) ≈  -81.43651914
@show gcm.∇β
@show gcm.∇τ
@show gcm.∇Σ

end


text:
cd("/home/graham_s/julia/vw/ScottishTaxBenefitModel" )
using Revise,Pkg
Pkg.activate(".")
# cd( relpath(@__FILE__)*"../")
include( "../test/testutils.jl" )

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
using Random
using CSV
using DataFrames

function dataloader(name; header = true, dir = "learningdatasets")
    hd = homedir()
    dir = joinpath(hd, dir)
    mkpath(dir)
    cd(dir)
    if name == "MNIST"
        if !(isfile("mnist_train.csv"))
            yn = Base.prompt("Can I download the mnist_train.csv(110MB)? (y/n)")
            yn == "y" && download("https://pjreddie.com/media/files/mnist_train.csv", "mnist_train.csv")
        end
        df_tr = CSV.read("mnist_train.csv", header = false, DataFrame)
        if !(isfile("mnist_test.csv"))
            yn = Base.prompt("Can I download the mnist_test.csv(18MB)? (y/n)")
            yn == "y" && download("https://pjreddie.com/media/files/mnist_test.csv", "mnist_test.csv")
        end
        df_te = CSV.read("mnist_test.csv",header = false , DataFrame)
        return df_tr, df_te
    elseif name == "iris"
        url = "https://gist.githubusercontent.com/netj/8836201/raw/6f9306ad21398ea43cba4f7d537619d0e07d5ae3/iris.csv"
        if !(isfile("iris.csv"))
            yn = Base.prompt("Can I download the iris.csv(4KB)? (y/n)")
            yn == "y" && download(url, "iris.csv")
        end
        df = CSV.read("iris.csv", DataFrame)
        return df
    elseif name == "BostonHousing"
        url = "https://raw.githubusercontent.com/selva86/datasets/master/BostonHousing.csv"
        if !(isfile("BostonHousing.csv"))
            yn = Base.prompt("Can I download the BostonHousing.csv(36KB)? (y/n)")
            yn == "y" && download(url, "BostonHousing.csv")
        end
        df = CSV.read("BostonHousing.csv", DataFrame)
        return df
    end
    try
        df = CSV.read(name, header = header, DataFrame)
        return df
    catch
        @warn "$name doesn't exist, so you can't get the data!"
    end
end

struct DataSplitter
    #TODO: make it possible to DataFrame split
    splitted_data::Tuple
    function DataSplitter(xs, narray, ndata, sb::Bool, seed, test_size::Int)
        if sb == true
            seed != nothing && Random.seed!(seed)
            narray = shuffle(narray)
        end
        narray = (narray[1:test_size-1], narray[test_size:end])
        splitted_data = []
        for x in xs
            push!(splitted_data, index2element(x, narray))
        end
        new(Tuple(splitted_data))
    end
end

function index2element(x, index::Tuple)
    tr, te = index
    if ndims(x) == 1
        return ([x[i] for i in tr], [x[i] for i in te])
    else
        return ([x[i, :] for i in tr], [x[i, :] for i in te])
    end
end

#ndata is dims = 1
function split_data(x...; sb = true, train_size = nothing, test_size = nothing, seed = nothing)
    x = collect(x)
    ndata = size(x[1])[1]
    narray = collect(1 : 1 : ndata)
    test_size == train_size == nothing && throw(ArgumentError("test_size and train_size wasn't specified. you must specify either one."))
    if test_size == nothing
        #number of train data is ceiled
        test_size = (0<=train_size<=1) ? ndata-ceil(train_size*ndata) : ndata-train_size
    else
        test_size = (0<=test_size<=1) ? ndata-ceil(train_size*ndata) : ndata-train_size
    end
    return DataSplitter(x, narray, ndata, sb, seed, Int(test_size)).splitted_data
end


text:
push!(LOAD_PATH, "src/")
using ArgParse,FASTX, nanoJulia, DataFrames, CSV

# Arguments parser
function parse_commandline()
	s = ArgParseSettings()
	@add_arg_table! s begin
		"--recursive", "-r"
			help = "Search recursively through folders for fastq files"
			action = :store_true
		"--inputdir", "-i"
			help = "The path to directory of fastq files"
			arg_type = String
			required = true
		"--outputdir", "-o"
			help = "Folder to output QC results"
			required = true
	end

	return parse_args(s)
end

function main()
	parsed_args = parse_commandline()
	recursive = parsed_args["recursive"]
	inputdir = parsed_args["inputdir"]
	outputdir = parsed_args["outputdir"]

	if isdir(outputdir)
		println("\033[1;31m[Warning] The output folder has already existed, so it will be renamed with additional suffix \"_bak\"\033[0m")
		new_outputdir = string(outputdir, "_bak")
		mv(outputdir, new_outputdir)
	end
	mkdir(outputdir)

	length2qualityTextFile = joinpath(outputdir, "readlength_vs_readquality.tsv")

	if recursive
		println("\033[1;32m* Search fast5 files recursively.......\033[0m")
		filelist = []
		for (root, dirs, files) in walkdir(inputdir)
			for file in files
				if occursin(r"fastq$|fq$", file)
					push!(filelist, joinpath(root, file))
				end
			end
		end
	else
		println("\033[1;32m* Search fast5 files.......\033[0m")
		fqfiles = [fq for fq in readdir(inputdir) if occursin(r"fastq$|fq$", fq)]
		filelist = joinpath.(inputdir, readdir(inputdir))
	end

	numFiles = length(filelist)
	println("\033[1m  - There're $numFiles fastq files required to be parsed!\033[0m")
	println("\033[1;32m* Start parsing fastq files...\033[0m")
	output_table = mapreduce(f -> readFastq(f), vcat, filelist)
	println("\033[1;32m* Writing output file...\033[0m")
	output_table |> CSV.write(length2qualityTextFile)
	totalLength = nanoJulia.totalLen(output_table.length)
	n50 = nanoJulia.calculate_N50(output_table.length, totalLength)
	println("\033[1;32m* Generating summry file and figures...\033[0m")
	generateStatSummary(output_table, totalLength, n50, outputdir)
end

main()


text:
using Plots

function midpoint_method(f::Function, x, u, t, h)
    k1, y = f(x, u, t)
    k2, _ = f(x + h/2*k1, u, t + h/2)
    k3, _ = f(x - h*k1 + 2*h*k2, u, t + h)
    x_VPG = x + h*k2
    d_VPG = h/6*(k1 - 2*k2 + k3)
    return x_VPG, y, d_VPG
end

function pt1(x, u)
    Tm = 10
    xdot = -1/Tm * x + 1/Tm * u
    y = x
    return xdot, y
end

function subtractor(u1, u2)
    return u1 - u2
end

function hysteresis(u)
    h_e = 0.085
    h_a = 0.065
    if u >= h_e
        y = 1;
    elseif u <= -h_e
        y = -1;
    else
        y = 0;
    end

    if i > 1
        if hys_state[i - 1] == 1 && u <= h_a
            y = 0;
        elseif hys_state[i - 1] == -1 && u >= -h_a
            y = 0;
        end
    end
    append!(hys_state, y)
    return y
end

function sys_PRM(x, u, t)
    _, y1::Float64 = pt1(x, 0)
    y2::Float64 = subtractor(u, y1)
    y3::Float64 = hysteresis(y2)
    xdot1, _ = pt1(x, y3)
    return xdot1, [y1 y2 y3]
end

t_step = 1;
Tm = 10
t0 = 0
tf = 20
epsilon = 1e-3
h_max = 2*Tm
h_min = 12*epsilon
u_values = []
x_values = []
y_values = Vector{Float64}()
d_values = []
t_values = []
h_values = []
hys_state = []
x = 0
h = 0.2
t = 0
i = 1
while t <= tf
    if t < t_step
        global u = 0
    else
        global u = 0.17
    end

    global x, y, d = midpoint_method(sys_PRM, x, u, t, h);
    global h = 0.9*h*min(max((epsilon/max(abs(d)))^(1/3), 0.3), 2)
    global h = min(max(h, h_min), h_max)

    append!(x_values, x)
    append!(y_values, y)
    append!(d_values, d)
    append!(h_values, h)
    append!(t_values, t)

    global t = t + h
    global i = i + 1
end

y_values = reshape(y_values, (3, :))

p1 = plot(t_values, y_values[1, :])
p2 = plot(t_values, y_values[2, :])
p3 = plot(t_values, y_values[3, :])
p4 = plot(t_values, h_values)
p5 = plot(t_values, d_values)

plot(p1, p2, p3, p4, p5, layout=(5, 1))

max_stars_repo_name: VeroAsghar/PRM

max_issues_repo_name: VeroAsghar/PRM

max_forks_repo_name: VeroAsghar/PRM


text:
module HTTPClient

include("HTTPC.jl")

end


text:
module FLOAT65

import Base: convert, promote_rule, reinterpret,
             show, showcompact, read, write,
             sizeof, hash,
             typemin, typemax, realmin, realmax,
             isinf, isnan, issubnormal, isfinite, isinteger,
             zero, one, sign, signbit, abs, eps,
             (==),(!=),(<),(<=),(>=),(>), isless, isequal,
             (+),(-),(*),(/),(\),(%),(^),sqrt,
             (.+),(.-),(.*),(./),(.\),(.^),
             min, max, minmax,
             div, fld, rem, mod, mod1, rem1, fld1,
             exp, expm1, log, log1p,
             sin, cos, tan, csc, sec, cot,
             asin, acos, atan, acsc, asec, acot,
             sinh, cosh, tanh, csch, sech, coth,
             asinh, acosh, atanh, acsch, asech, acoth

export Float65,
       getstate, setstate, clearstate

include("type/abstract.jl")
include("type/Float65.jl")
include("type/conversion.jl")
include("type/basics.jl")
include("type/predicates.jl")
include("type/compare.jl")
include("type/io.jl")
include("type/wrapfuns.jl")

include("behavior/Reinterpret.jl")
include("behavior/ebit.jl")
include("behavior/tinyhuge.jl")
include("behavior/absorb.jl")

if (isdefined(Main,:NoStatePropagation) && Main.NoStatePropagatione==true)

include("math/assignState/arith.jl")
include("math/assignState/elementary.jl")

else # propagate the set state

include("math/propagateState/arith.jl")
include("math/propagateState/elementary.jl")

end

if isdefined(Main, :Distributions)
    include("interop/Distributions.jl")
end

end # module


text:
module ReactiveMPTestingHelpers

using Test
using ReactiveMP

import ReactiveMP: OneDivNVector
import ReactiveMP: deep_eltype

@testset "Helpers" begin

    @testset "OneDivNVector" begin

        for type in [ Float64, Float32, BigFloat ], len in [ 3, 5, 10 ]
            iter = OneDivNVector(type, len)

            @test eltype(iter) === type
            @test length(iter) === len
            @test size(iter)   === (len, )
            @test collect(iter) == fill(one(type) / len, len)
            @test vec(iter)     == fill(one(type) / len, len)
            @test sizeof(iter) === 0

            sim = similar(iter)

            @test sim === iter
            @test eltype(sim) === type
            @test length(sim) === len
            @test size(sim)   === (len, )
            @test collect(sim) == fill(one(type) / len, len)
            @test vec(sim)     == fill(one(type) / len, len)
            @test sizeof(sim) === 0

            sim = similar(iter, Float16)
            @test eltype(sim) === Float16
            @test length(sim) === len
            @test size(sim)   === (len, )
            @test collect(sim) == fill(one(Float16) / len, len)
            @test vec(sim)     == fill(one(Float16) / len, len)
            @test sizeof(sim) === 0
        end

        @test eltype(OneDivNVector(3)) === Float64
        @test eltype(OneDivNVector(5)) === Float64

        @test_throws AssertionError OneDivNVector(-2)
        @test_throws AssertionError OneDivNVector(-10)
        @test_throws AssertionError OneDivNVector(Vector, 2)
        @test_throws AssertionError OneDivNVector(Matrix, 10)

    end

    @testset "deep_eltype" begin

        for type in [ Float32, Float64, Complex{Float64}, BigFloat ]

            @test deep_eltype(type) === type
            @test deep_eltype(zero(type)) === type

            vector = zeros(type, 10)
            matrix = zeros(type, 10, 10)
            vector_of_vectors  = [ vector, vector ]
            vector_of_matrices = [ matrix, matrix ]
            matrix_of_vector   = [ vector vector; vector vector ]
            matrix_of_matrices = [ matrix matrix; matrix matrix ]

            @test deep_eltype(vector) === type
            @test deep_eltype(matrix) === type
            @test deep_eltype(vector_of_vectors) === type
            @test deep_eltype(vector_of_matrices) === type
            @test deep_eltype(matrix_of_vector) === type
            @test deep_eltype(matrix_of_matrices) === type

        end

    end

end

end


text:
"""
Subroutine dmid

This subroutine computes the projection of x
on the n-dimensional interval [xl,xu].

MINPACK-2 Project. March 1999.
Argonne National Laboratory.
Chih-Jen Lin and Jorge J. More'.
"""
function dmid(n,x,xl,xu)
    @inbounds for i=1:n
        x[i] = max(xl[i],min(x[i],xu[i]))
    end
    return
end

@inline function dmid(n::Int, x::CuDeviceArray{Float64,1},
                      xl::CuDeviceArray{Float64,1}, xu::CuDeviceArray{Float64,1})
    tx = threadIdx().x
    ty = threadIdx().y

    if tx <= n && ty == 1
        @inbounds x[tx] = max(xl[tx], min(x[tx], xu[tx]))
    end
    CUDA.sync_threads()

    return
end


text:
export barrier_and!, barrier_or!

mutable struct HSAKernelInstance{T<:Tuple}
    agent::HSAAgent
    exe::HSAExecutable
    sym::String
    args::T
    kernel_object::UInt64
    kernarg_segment_size::UInt32
    group_segment_size::UInt32
    private_segment_size::UInt32
    kernarg_address::Ptr{Nothing}
end

# TODO agent can be inferred from the executable
function HSAKernelInstance(agent::HSAAgent, exe::HSAExecutable, symbol::String, args::Tuple)
    agent_func = @cfunction(iterate_exec_agent_syms_cb, HSA.Status,
                           (HSA.Executable, HSA.Agent, HSA.ExecutableSymbol, Ptr{Cvoid}))
    exec_symbol = Ref{HSA.ExecutableSymbol}()
    exec_ptr = Base.unsafe_convert(Ref{Cvoid}, exec_symbol)
    @debug "Agent symbols"
    HSA.executable_iterate_agent_symbols(exe.executable[], agent.agent,
                                         agent_func, exec_ptr) |> check

    # TODO: Conditionally disable once ROCR is fixed
    if !isassigned(exec_symbol)
        agent_ref = Ref(agent.agent)
        GC.@preserve agent_ref begin
            HSA.executable_get_symbol_by_name(exe.executable[], symbol,
                                              agent_ref, exec_symbol) |> check
        end
    end

    kernel_object = Ref{UInt64}(0)
    getinfo(exec_symbol[], HSA.EXECUTABLE_SYMBOL_INFO_KERNEL_OBJECT,
            kernel_object) |> check

    kernarg_segment_size = Ref{UInt32}(0)
    getinfo(exec_symbol[], HSA.EXECUTABLE_SYMBOL_INFO_KERNEL_KERNARG_SEGMENT_SIZE,
            kernarg_segment_size) |> check

    if kernarg_segment_size[] == 0
        # FIXME: Hidden arguments!
        if length(args) > 0
            kernarg_segment_size[] = sum(sizeof.(args))
        else
            # Allocate some memory anyway, #10
            kernarg_segment_size[] = max(kernarg_segment_size[], 8)
        end
    end

    group_segment_size = Ref{UInt32}(0)
    getinfo(exec_symbol[], HSA.EXECUTABLE_SYMBOL_INFO_KERNEL_GROUP_SEGMENT_SIZE,
            group_segment_size) |> check

    private_segment_size = Ref{UInt32}(0)
    getinfo(exec_symbol[], HSA.EXECUTABLE_SYMBOL_INFO_KERNEL_PRIVATE_SEGMENT_SIZE,
            private_segment_size) |> check

    finegrained_region = get_region(agent, :finegrained)
    kernarg_region = get_region(agent, :kernarg)

    # Allocate the kernel argument buffer from the correct region
    kernarg_address = Ref{Ptr{Nothing}}()
    HSA.memory_allocate(kernarg_region[],
                        kernarg_segment_size[],
                        kernarg_address) |> check

    # Fill kernel argument buffer
    # FIXME: Query kernarg segment alignment
    ctr = 0x0
    for arg in args
        rarg = Ref(arg)
        align = Base.datatype_alignment(typeof(arg))
        rem = mod(ctr, align)
        if rem > 0
            ctr += align-rem
        end
        #@info "Storing $(typeof(arg)) at offset $ctr"
        ccall(:memcpy, Cvoid,
            (Ptr{Cvoid}, Ptr{Cvoid}, Csize_t),
            kernarg_address[]+ctr, rarg, sizeof(arg))
        ctr += sizeof(arg)
    end

    kernel = HSAKernelInstance(agent, exe, symbol, args, kernel_object[],
                               kernarg_segment_size[], group_segment_size[],
                               private_segment_size[], kernarg_address[])
    hsaref!()
    finalizer(kernel) do kernel
        HSA.memory_free(kernel.kernarg_address) |> check
        hsaunref!()
    end
    return kernel
end

struct HSASignalSet{T}
    signals::Vector{HSASignal}
end
HSASignalSet{T}() where T = HSASignalSet{T}(HSASignal[])
Base.wait(ss::HSASignalSet{HSA.BarrierAndPacket}) = wait.(ss.signals)
function Base.wait(ss::HSASignalSet{HSA.BarrierOrPacket})
    #= FIXME
    # We need to hack around the fact that barrier OR packets don't handle more
    # than 5 dependencies. We could implement the waiting in software, and emit
    # a barrier that waits on a signal tied to that waiter.
    =#
    @warn "Waiting on OR barriers waits on all partitioned barriers to complete" maxlog=1
    wait.(ss.signals)
end

barrier_and!(signals::Vector) = barrier_and!(default_queue().queue, signals)
barrier_or!(signals::Vector) = barrier_or!(default_queue().queue, signals)
barrier_and!(queue::HSAQueue, signals::Vector{HSA.Signal}) = barrier_and!(queue, map(HSASignal, signals))
barrier_or!(queue::HSAQueue, signals::Vector{HSA.Signal}) = barrier_or!(queue, map(HSASignal, signals))
barrier_and!(queue::HSAQueue, signals::Vector{HSASignal}) =
    barrier!(HSA.BarrierAndPacket, queue, signals)
barrier_or!(queue::HSAQueue, signals::Vector{HSASignal}) =
    barrier!(HSA.BarrierOrPacket, queue, signals)
function barrier!(T, queue::HSAQueue, signals::Vector{HSASignal})
    outset = HSASignalSet{T}()
    if !isempty(signals)
        for signal_set in Iterators.partition(signals, 5)
            comp_signal = HSASignal()
            _launch!(T, queue, comp_signal) do _packet
                @set! _packet.dep_signal = ntuple(i->length(signal_set)>=i ? signal_set[i].signal[] : HSA.Signal(0), 5)
                _packet
            end
            push!(outset.signals, comp_signal)
        end
    end
    return outset
end

# TODO Docstring
function launch!(queue::HSAQueue, kernel::HSAKernelInstance, signal::HSASignal;
                 workgroup_size=nothing, grid_size=nothing)
    @assert workgroup_size !== nothing "Must specify workgroup_size kwarg"
    @assert grid_size !== nothing "Must specify grid_size kwarg"

    _launch!(HSA.KernelDispatchPacket, queue, signal) do _packet
        @set! _packet.setup = 3 << Int(HSA.KERNEL_DISPATCH_PACKET_SETUP_DIMENSIONS)
        @set! _packet.workgroup_size_x = workgroup_size[1]
        @set! _packet.workgroup_size_y = workgroup_size[2]
        @set! _packet.workgroup_size_z = workgroup_size[3]
        @set! _packet.grid_size_x = grid_size[1]
        @set! _packet.grid_size_y = grid_size[2]
        @set! _packet.grid_size_z = grid_size[3]
        @set! _packet.completion_signal = signal.signal[]
        @set! _packet.kernel_object = kernel.kernel_object
        @set! _packet.kernarg_address = kernel.kernarg_address
        @set! _packet.private_segment_size = kernel.private_segment_size
        @set! _packet.group_segment_size = kernel.group_segment_size
        _packet
    end
end

function _launch!(f, T, queue::HSAQueue, signal::HSASignal)
    # Obtain the current queue write index and queue size
    _queue_size = Ref{UInt32}(0)
    getinfo(queue.agent.agent, HSA.AGENT_INFO_QUEUE_MAX_SIZE, _queue_size) |> check
    queue_size = _queue_size[]
    write_index = HSA.queue_add_write_index_scacq_screl(queue.queue[], UInt64(1))

    # Yield until queue has space
    while true
        read_index = HSA.queue_load_read_index_scacquire(queue.queue[])
        if write_index < read_index + queue_size
            break
        end
        yield()
    end

    # TODO: Make this less ugly
    dispatch_packet = Ref{T}()
    ccall(:memset, Cvoid,
          (Ptr{Cvoid}, Cint, Csize_t),
          dispatch_packet, 0, sizeof(T))
    _packet = dispatch_packet[]
    _packet = f(_packet)
    dispatch_packet = Ref{T}(_packet)

    _queue = unsafe_load(queue.queue[])
    queueMask = UInt32(_queue.size - 1)
    baseaddr_ptr = Ptr{HSA.KernelDispatchPacket}(_queue.base_address)
    baseaddr_ptr += sizeof(HSA.KernelDispatchPacket) * (write_index & queueMask)
    dispatch_packet_ptr = convert(Ptr{HSA.KernelDispatchPacket}, Base.unsafe_convert(Ptr{T}, dispatch_packet))
    unsafe_copyto!(baseaddr_ptr, dispatch_packet_ptr, 1)

    packetheadertype(::Type{HSA.KernelDispatchPacket}) = HSA.PACKET_TYPE_KERNEL_DISPATCH
    packetheadertype(::Type{HSA.BarrierAndPacket}) = HSA.PACKET_TYPE_BARRIER_AND
    packetheadertype(::Type{HSA.BarrierOrPacket}) = HSA.PACKET_TYPE_BARRIER_OR

    # Atomically store the header
    header = Ref{UInt16}(0)
    header[] |= Int(HSA.FENCE_SCOPE_SYSTEM) << Int(HSA.PACKET_HEADER_ACQUIRE_FENCE_SCOPE)
    header[] |= Int(HSA.FENCE_SCOPE_SYSTEM) << Int(HSA.PACKET_HEADER_RELEASE_FENCE_SCOPE)
    header[] |= Int(packetheadertype(T)) << Int(HSA.PACKET_HEADER_TYPE)
    atomic_store_n!(Base.unsafe_convert(Ptr{UInt16}, baseaddr_ptr), header[])

    # Ring the doorbell to dispatch the kernel
    HSA.signal_store_relaxed(_queue.doorbell_signal, Int64(write_index))
end


text:
import Base: info, warn, toc
import Memento: debug

export info, warn, debug, tocc,
       configure_logger, add_file_logger

global LOGGER = Memento.getlogger(current_module())

const DEFAULT_LOG_LEVEL = "warn"

"""
    info(msg)

Prints a message on info level.

### Input

- `msg` - message string
"""
@inline function info(msg::String)
    Memento.info(LOGGER, msg)
end

"""
    warn(msg)

Prints a message on warn level.

### Input

- `msg` - message string
"""
@inline function warn(msg::String)
    Memento.warn(LOGGER, msg)
end

"""
    debug(msg)

Prints a message on debug level.

### Input

- `msg` - message string
"""
@inline function debug(msg::String)
    Memento.debug(LOGGER, msg)
end

"""
    tocc(func)

Prints the output of `toc()` using the given log function (default: `info`).

### Input

- `func` - (optional, default: `info`) log function
"""
@inline function tocc(func::Function=info)
    func("elapsed time: $(toq()) seconds")
end

"""
    configure_logger(level)

Configures the global log level. If no log level is passed, we use the log level
that is defined by the constant `DEFAULT_LOG_LEVEL`.

### Input

- `level` - (optional) the log level; can be either an integer between 0 and 2
            or a string that is supported by the
            [Memento.jl](https://invenia.github.io/Memento.jl/latest/man/intro.html#Logging-levels-1)
            package.
"""
function configure_logger(level::Union{String, Int, Void}=DEFAULT_LOG_LEVEL)
    if level isa String
        level_string = level
    elseif level isa Int
        if level == 0
            level_string = "warn"
        elseif level == 1
            level_string = "info"
        elseif level == 2
            level_string = "debug"
        else
            error("Illegal verbosity input $level.")
        end
    else
        error("Illegal verbosity input $level.")
    end
    return Memento.config!(level_string; fmt="[{level}] {msg}")
end

"""
    add_file_logger(filename)

Sets up an additional logger to a file.

### Input

- `filename` -- (optional, default: `tempname()`) the log file name
"""
function add_file_logger(filename::String=tempname())::Void
    formatter = Memento.DefaultFormatter("[{date}|{level}] {msg}")
    handler = Memento.DefaultHandler(filename, formatter)
    push!(LOGGER, handler)
    return nothing
end


text:
# Allows us to use `safesave(filename, chain)` to ensure that we do not overwrite any chains
DrWatson._wsave(filename, chain::Chains) = write(filename, chain)

max_stars_repo_name: andreaskoher/CovidSurvey

max_issues_repo_name: andreaskoher/CovidSurvey

max_forks_repo_name: andreaskoher/CovidSurvey


text:
abstract type RinexBase end
abstract type RinexHeader <: RinexBase end
abstract type RinexData <: RinexBase end

#= =========================== =#

@with_kw struct RinexNavHeader <: RinexHeader
    version::Tuple{Float64,String,String,SatelliteSystemType}
    runBy::Tuple{String,String,String}
    comment::String
    ionoCorr::String
    timeSysCorr::Dict{TimeCorrectionType,TimeSystemCorrection}
    leapSeconds::Tuple{Int64,Int64,Int64,Int64}
    corrSysTime::String
    deltaUTC::String
    dUTC::String
    ionAlpha::String
    ionBeta::String
    eoH::String
end

@with_kw struct EphEpoch
    si::SatID
    dt::DateTime
    d::Array{Float64,1}
end

@with_kw mutable struct EphEpochStore
    si::SatID
    keys::Array{DateTime,1}
    eph::Dict{DateTime,EphEpoch}
end

@with_kw struct RinexNavData <: RinexData
    data::Dict{SatID,EphEpochStore}
end

@with_kw struct RinexNav <: RinexBase
    filename::String
    header::RinexHeader
    data::RinexNavData
end

# TODO multiple data files

#= =========================== =#

D_SystemNumObs = Dict{SatelliteSystemType,Tuple{Int64,Array{ObsID,1}}}

@with_kw struct RinexObsHeader <: RinexHeader
    version::Tuple{Float64,Char,SatelliteSystemType}
    runBy::Tuple{String,String,String}
    comment::String
    markerName::String
    markerNumber::String
    markerType::String
    observer::Tuple{String,String}
    receiver::String
    antennaType::String
    antennaPosition::ECEF{Float64}
    antennaDeltaHEN::String
    antennaDeltaXYZ::String
    antennaPhaseCtr::String
    antennaBsightXYZ::String
    antennaZeroDirAzi::String
    antennaZeroDirXYZ::String
    centerOfMass::String
    numObs::String
    systemNumObs::D_SystemNumObs
    waveFact::String
    sigStrengthUnit::String
    interval::Float64
    firstTime::Tuple{DateTime,Millisecond}
    lastTime::String
    receiverOffset::String
    systemDCBSapplied::String
    systemPCVSapplied::String
    systemScaleFac::String
    systemPhaseShift::String
    glonassSlotFreqNo::String
    glonassCodPhsBias::String
    leapSeconds::String
    numSats::String
    prnObs::String
    eoH::String
    iscrx::Bool = false
    # recommend_length::Int64
    # depended on first obs last obs and interval
    # 2400 => 1 days 15 s interval
end

struct RinexDatum
    obs::Float64 # empty -> NaN
    lli::Char
    ssi::Char
    function RinexDatum(str::Union{String,SubString})
        v = strip(SubString(str, 1, 13))
        if isempty(v)
            obs = NaN
        else
            obs = parse(Float64, v)
        end
        llic = str[14]
        ssic = str[15]
        new(obs, llic, ssic)
    end
end

struct EpochHeader
    dt::DateTime
    flag::Char
    numObs::Int64
    RCO::Float64
end

@with_kw struct SatelliteObsData
    satid::SatID
    obsids::Array{ObsID,1}
    dts::Array{DateTime,1}
    obs::Array{Array{RinexDatum,1},1}
end
@with_kw struct RinexObsData <: RinexData
    epochs::Dict{DateTime,EpochHeader}
    comments::Dict{DateTime,String}
    obs::Dict{SatID,SatelliteObsData}
end

struct RinexObs <: RinexBase
    filename::String
    header::RinexObsHeader
    data::RinexObsData
end

# TODO multiple data files => cat?


text:
export class_likelihood_per_instance, class_weights_per_instance

using CUDA
using LoopVectorization: @avx, vifelse

"""
Class Conditional Probability
"""
function class_likelihood_per_instance(lc::LogisticCircuit, nc::Int, data)
    cw = class_weights_per_instance(lc, nc, data)
    one = Float32(1.0)
    isgpu(data) ? (@. one / (one + exp(-cw))) : (@. @avx one / (one + exp(-cw)))
end

function class_likelihood_per_instance(bc, data)
    cw = class_weights_per_instance(bc, data)
    one = Float32(1.0)
    isgpu(data) ? (@. one / (one + exp(-cw))) : (@. @avx one / (one + exp(-cw)))
end

function class_weights_per_instance(lc::LogisticCircuit, nc::Int, data)
    bc = ParamBitCircuit(lc, nc, data)
    class_weights_per_instance(bc, data)
end

function class_weights_per_instance(bc, data)
    if isgpu(data)
        class_weights_per_instance_gpu(to_gpu(bc), data)
    else
        class_weights_per_instance_cpu(bc, data)
    end
end

function class_weights_per_instance_cpu(bc, data)
    ne::Int = num_examples(data)
    nc::Int = size(bc.params, 2)
    cw::Matrix{Float32} = zeros(Float32, ne, nc)
    cw_lock::Threads.ReentrantLock = Threads.ReentrantLock()

    @inline function on_edge_binary(flows, values, prime, sub, element, grandpa, single_child, weights)
        lock(cw_lock) do # TODO: move lock to inner loop?
            for i = 1:size(flows, 1)
                @inbounds edge_flow = values[i, prime] & values[i, sub] & flows[i, grandpa]
                first_true_bit = trailing_zeros(edge_flow) + 1
                last_true_bit = 64 - leading_zeros(edge_flow)
                @simd for j = first_true_bit:last_true_bit
                    ex_id = ((i - 1) << 6) + j
                    if get_bit(edge_flow, j)
                        for class = 1:nc
                            @inbounds cw[ex_id, class] += bc.params[element, class]
                        end
                    end
                end
            end
        end
        nothing
    end

    @inline function on_edge_float(flows, values, prime, sub, element, grandpa, single_child, weights)
        lock(cw_lock) do # TODO: move lock to inner loop?
            @avx for i = 1:size(flows, 1)
                @inbounds edge_flow = values[i, prime] * values[i, sub] / values[i, grandpa] * flows[i, grandpa]
                edge_flow = vifelse(isfinite(edge_flow), edge_flow, zero(eltype(flows)))
                for class = 1:nc
                    @inbounds cw[i, class] += edge_flow * bc.params[element, class]
                end
            end
        end
        nothing
    end

    if isbinarydata(data)
        satisfies_flows(bc.bitcircuit, data; on_edge = on_edge_binary)
    else
        satisfies_flows(bc.bitcircuit, data; on_edge = on_edge_float)
    end

    return cw
end

function class_weights_per_instance_gpu(bc, data)
    ne::Int = num_examples(data)
    nc::Int = size(bc.params, 2)
    cw::CuMatrix{Float32} = CUDA.zeros(Float32, num_examples(data), nc)
    cw_device = CUDA.cudaconvert(cw)
    params_device = CUDA.cudaconvert(bc.params)

    @inline function on_edge_binary(flows, values, prime, sub, element, grandpa, chunk_id, edge_flow, single_child, weights)
        first_true_bit = 1 + trailing_zeros(edge_flow)
        last_true_bit = 64 - leading_zeros(edge_flow)
        for j = first_true_bit:last_true_bit
            ex_id = ((chunk_id - 1) << 6) + j
            if get_bit(edge_flow, j)
                for class = 1:nc
                    CUDA.@atomic cw_device[ex_id, class] += params_device[element, class]
                end
            end
        end
        nothing
    end

    @inline function on_edge_float(flows, values, prime, sub, element, grandpa, ex_id, edge_flow, single_child, weights)
        for class = 1:nc
            CUDA.@atomic cw_device[ex_id, class] += edge_flow * params_device[element, class]
        end
        nothing
    end

    if isbinarydata(data)
        v,f = satisfies_flows(bc.bitcircuit, data; on_edge = on_edge_binary)
    else
        @assert isfpdata(data) "Only floating point and binary data are supported"
        v,f = satisfies_flows(bc.bitcircuit, data; on_edge = on_edge_float)
    end
    CUDA.unsafe_free!(v) # save the GC some effort
    CUDA.unsafe_free!(f) # save the GC some effort

    return cw
end

"""
Class Predictions
"""
function predict_class(lc::LogisticCircuit, nc::Int, data)
    class_likelihoods = class_likelihood_per_instance(lc, nc, data)
    predict_class(class_likelihoods)
end

function predict_class(class_likelihoods)
    _, mxindex = findmax(class_likelihoods; dims=2)
    dropdims(getindex.(mxindex, 2); dims=2)
end

"""
Prediction accuracy
"""
accuracy(lc::LogisticCircuit, nc::Int, data, labels::Vector) =
    accuracy(predict_class(lc, nc, data), labels)

accuracy(predicted_class::Vector, labels::Vector) =
    Float64(sum(@. predicted_class == labels)) / length(labels)

accuracy(class_likelihoods::Matrix, labels::Vector) =
    accuracy(predict_class(class_likelihoods), labels)


text:
export
    earth_pol_to_equ_radius,
    earth_flattening,
    earth_equ_radius,
    minutes_per_day,
    days_per_minute,
    seconds_per_day,
    days_per_second,
    km_per_au,
    standard_rst_altitude,
    sun_rst_altitude

# Don't change these unless you are moving to a new universe.

# Ratio of Earth's polar to equatorial radius.
global earth_pol_to_equ_radius = 0.99664719

# Earth flattening
global earth_flattening = 1/298.257

# Equatorial radius of the Earth in km.
global earth_equ_radius = 6378.14

# How many minutes in a day?
global minutes_per_day = 24 * 60

# How many days in minute?
global days_per_minute = 1 / minutes_per_day

# How many seconds (time) in a day?
global seconds_per_day = 24 * 60 * 60

# How many days in a second?
global days_per_second = 1 / seconds_per_day

# How many kilometres in an astronomical unit?
global km_per_au = 149597870

# For rise-set-transit: altitude deflection caused by refraction
global standard_rst_altitude = -0.00989078087105 # -0.5667 degrees
global sun_rst_altitude      = -0.0145438286569       # -0.8333 degrees

global currentobliquity = 23.5


text:
module TestTreeDP

using Test
import Printf: @sprintf
import TreeLas
import TreeLas.TreeDP
import TreeLas.Dual: dual
import GraphIdx
import GraphIdx.Tree: ChildrenIndex, hierarchy, hierarchy_with, hierarchy_string

include(joinpath(dirname(pathof(TreeLas)), "cxx.jl"))

mean(x) = sum(x) / length(x)

@testset "TreeDP                         " begin

@testset "tree_dp: tree.mini.h5          " begin
    t = TreeDP.Tree([1, 1, 2, 1, 2, 3])
    y = [0.53 0.45 0.02
         0.17 0.96 0.97]
    y = Array(y')
    lam = 0.2
    x = TreeDP.tree_dp(y, t, lam)
    @test size(x) == size(y)
    xt = [0.4  0.4  0.4
          0.37 0.76 0.77]
    xt = Array(xt')
    if false
        hierarchy_with(ChildrenIndex(t.parent)) do i
            @sprintf("%.2f:%.2f", y[i], xt[i])
        end
    end
    @test x ≈ xt
end

@testset "tree_dp: 10                    " begin
    t = TreeDP.Tree(Int[1, 1, 2, 3, 4, 1, 8, 9, 4, 9])
    y = Float64[8.2, 7.0, 9.5, 6.8, 5.8, 6.3, 4.3, 2.2, 1.2, 2.8]
    lam = 1.0

    x = TreeDP.tree_dp(copy(y), t, lam)
    @test size(x) == size(y)
    expect3x = [22.7, 22.7, 22.7, 18.9, 18.9, 21.9,  9.9,  8.2,  8.2,  8.2]
    @test 3x ≈ expect3x
    local cidx = ChildrenIndex(t.parent)
    if true
        let z = expect3x/3 - y
            local alpha = dual(z, t.parent, cidx)
            # @info "alpha=" *string(round.(alpha, digits=4))
            @debug hierarchy_string(
                cidx, i -> @sprintf("%.1f≈%+.1f", -3y[i]+expect3x[i], -3alpha[i])
            )
        end
    end
end

@testset "MGT     : demo 3x7             " begin
    y = [0.62 0.73 0.71 1.5 1.17 0.43 1.08
         0.62 1.73 0.95 1.46 1.6 1.16 0.38
         0.9 0.32 -0.48 0.95 1.08 0.02 0.4]
    local n = length(y)
    local root = 1
    edges = Tuple{Int64,Int64}[
        (1,  2),  (4,  5),  (7, 8), (10, 11), (13, 14), (16, 17),
        (19, 20), (2,  3),  (5, 6), (8,   9), (11, 12), (14, 15),
        (17, 18), (20, 21), (1, 4), (4,   7), (7,  10), (10, 13),
        (13, 16), (16, 19), (2, 5), (5,   8), (8,  11), (11, 14),
        (14, 17), (17, 20), (3, 6), (6,   9), (9,  12), (12, 15),
        (15, 18), (18, 21)
    ]

    γ = [-abs(y[u] - y[v]) for (u, v) in edges]

    parent =
        [1, 5, 6, 1, 4, 5, 8, 9, 6, 7, 8, 9, 10, 15, 18, 13, 16, 17, 20, 17, 18]

    local tree = GraphIdx.Tree.RootedTree(root, parent)
    local dp_mem = TreeDP.TreeDPMem(n)
    local x = fill(NaN, size(y)...)
    local tree_alpha = fill(NaN, n)
    let mu = TreeDP.Ones{Float64}(), z = copy(y)
        TreeDP.tree_dp!(x,
                        z,
                        tree,
                        TreeDP.Ones{Float64}(),
                        mu,
                        dp_mem)
        @test !any(isnan.(x))
    end

    @testset "analyse solution" begin
        @test length(unique(x)) == 2
        @test mean(y) ≈ mean(x)
    end

    @testset "compare to cxx tree_dp" begin
        cxx_x = Cxx.cxx_tree_dp(y, parent, root, 1.0)
        @test cxx_x ≈ x
    end

    @testset "dual!" begin
        @test dp_mem.proc_order[end] == root
        @test dp_mem.proc_order[1:end-1] ==
            [12, 11, 19, 20, 21, 14, 15, 18, 17, 16, 13, 10, 7, 8, 9, 3, 6, 2, 5, 4]
        TreeLas.Dual.dual!(tree_alpha, x, copy(y), dp_mem.proc_order, parent)
        local alpha = vec(x) - vec(y)
        for i in @view dp_mem.proc_order[1:end-1]
            alpha[parent[i]] += alpha[i]
        end
        if !(alpha[2:end] ≈ tree_alpha[2:end])
            @show alpha - tree_alpha
            @test alpha[2:end] ≈ tree_alpha[2:end]
        end
        exceed = findall(abs.(tree_alpha) .> 1 + 1e12)
        if exceed != []
            for ei in exceed
                @show ei
                @show alpha[ei]
            end
        end
    end

    @testset "compare to tree_dual" begin
        cxx_alpha = Cxx.cxx_tree_dual!(x - copy(y), parent, root)
        @test abs.(tree_alpha[2:end]) ≈ abs.(cxx_alpha[2:end])
        sign01(x) = Int((sign(x) + 1) / 2)
        if !(sign01.(tree_alpha[2:end]) == sign01.(cxx_alpha[2:end]))
            hierarchy(ChildrenIndex(parent))
            @test sign01.(tree_alpha[2:end]) == sign01.(cxx_alpha[2:end])
        end
        @test tree_alpha[2:end] ≈ cxx_alpha[2:end]
    end

    local wtree = GraphIdx.Tree.WeightedTree(tree, TreeDP.Ones{Float64}())
    let gam = fill(NaN, n)
        TreeLas.MGT.gap_vec!(gam, x, tree_alpha, wtree)
        @test isnan(gam[root])
        gam[root] = 0.0
        wrong_ids = [i for i in 1:n if abs(gam[i]) > 1e-16]
        if wrong_ids != []
            for i in wrong_ids
                @show i, parent[i]
                @show x[i], x[parent[i]]
                @show tree_alpha[i]
                @show gam[i]
            end
            hierarchy(ChildrenIndex(parent))
            @test wrong_ids == []
        end
    end
end

end
end

max_stars_repo_name: EQt/treelas

max_issues_repo_name: EQt/treelas

max_forks_repo_name: EQt/treelas


text:
# zip of uneven iterators is still not supported in julia... super weird
# see https://discourse.julialang.org/t/collecting-zip/20739
# and https://github.com/JuliaLang/julia/issues/17928
function Base.collect(itr::Base.Iterators.Zip)
  itrsize = Base.IteratorSize(itr)
  itrsize isa Base.HasShape && (itrsize = Base.HasLength())
  Base._collect(1:1, itr, Base.IteratorEltype(itr), itrsize)
end

"""
only continue if true, else return false immediately
"""
macro iftrue(expr)
  quote
    $(esc(expr)) || return false
  end
end

same_expr(e1, e2) = e1 == e2
function same_expr(e1::Expr, e2::Expr)
  # don't differentiate between `f(a) = a` and `function f(a); a; end`
  e1.head == :function && (e1.head = :(=))
  e2.head == :function && (e2.head = :(=))
  @iftrue e1.head == e2.head

  args1 = filter(x -> !isa(x, LineNumberNode), e1.args)
  args2 = filter(x -> !isa(x, LineNumberNode), e2.args)
  @iftrue length(args1) == length(args2)
  # recurse
  all(zip(args1, args2)) do (a, b)
    same_expr(a, b)
  end
end

@test same_expr(:a, :a)
@test same_expr(:(a = 4), :(a = 4))
@test same_expr(quote
  f(a) = a
end, quote
  function f(a)
    a
  end
end)
@test !same_expr(:a, :b)
@test !same_expr(:(a = 4), :(a = 5))
@test !same_expr(quote
  f(a) = a
end, quote
  f(b) = b
end)