text:
@testset "The Hyperbolic Manifold" begin
    M = Hyperbolic(2)
    @test typicalDistance(M) == sqrt(2)
    @test getValue(HnPoint(1.)) == 1.
    @test getValue(HnTVector(0.)) == 0.
    x = HnPoint([0.,0.,1])
    @test getValue(x) == [0., 0., 1.]
    ξ = HnTVector([1.,1.,0,])
    @test getValue(ξ) == [1., 1., 0.]
    @test_throws ErrorException validateMPoint(M, HnPoint([1.,0.]))
    @test_throws ErrorException validateMPoint(M, HnPoint([0., 0., 0.]))
    @test_throws ErrorException validateMPoint(M, HnPoint([2., 0., 3.])) # -5 instead of -1
    @test validateTVector(M, x, ξ)
    @test_throws ErrorException validateTVector(M,x, HnTVector([0.,0.])) # Dimensions don't agree
    @test_throws ErrorException validateTVector(M,x, HnTVector(getValue(x))) # Not orthogonal
    y = exp(M, x, ξ)
    @test distance(M, x, y) ≈ norm(M, x, ξ)
    @test norm(M, x, ξ) ≈ sqrt(dot(M,x,ξ,ξ))
    @test distance(M, exp(M,x,log(M,x,y)), y) ≈ 0
    @test distance(M, exp(M,x,zeroTVector(M,x)), x) ≈ 0
    @test norm(M, x, log(M,x,x)) ≈ 0
    #
    @test manifoldDimension(x) == 2
    @test manifoldDimension(M) == 2
    @test parallelTransport(M,x,x,ξ) == ξ
    @test norm(M, y, parallelTransport(M,x,y,ξ)+log(M,y,x)) ≈ 0 atol=5*10^(-15)
    #
    @test typeofTVector(typeof(x)) == HnTVector{Float64}
    @test typeofMPoint(typeof(ξ)) == HnPoint{Float64}
    @test typeofTVector(x) == HnTVector{Float64}
    @test typeofMPoint(ξ) == HnPoint{Float64}

    @test "$(M)" == "The 2-Hyperbolic Space."
    @test "$(x)" == "Hn([0.0, 0.0, 1.0])"
    @test "$(ξ)" == "HnT([1.0, 1.0, 0.0])"
end


text:
function julia_main()::Cint
    try
        real_main()
    catch
        Base.invokelatest(Base.display_error, Base.catch_stack())
        return 1
    end
    return 0
end

function real_main()
    # TODO: "real" argument handling
    infile = ARGS[1]
    outdir = ARGS[2]
    mesolve_qutip(infile, outdir; backend=:CUDA)
end

ext: jl


text:
#=
    Bayesian Data Analysis of Titanic data set

    This is also a test of the ARFFFiles package for reading ARFF Files
    (the format used in the data mining system Weka, which happens to
    be one of my favorite data mining system).

    The file used is from my Weka page: http://www.hakank.org/weka/titanic.arff

    This is the simple version of the Titanic dataset with 2201 rows and the following
    attributes, their values (and the transformed values):
    - class: 1st, 2nd, 3rd, crew (=> 1, 2, 3, 4)
    - age  : adult, child        (=> 1, 2)
    - sex  : male, female        (=> 1, 2)
    - survived: yes, no          (=> 1, 0)

    Compared to titanic.jl, this program tests the model with a different train and
    test sets. The accuracy for the titanic.jl model is about 77.6% (using the full dataset).

    Here are some results with different training percentages.

    * train pct 0.99
    Num correct: 19.0
    Num incorrect: 3.0
    Accuracy: 0.8636363636363636

    Confusion matrix:
    2×2 Matrix{Float64}:
    17.0  0.0
     3.0  2.0

    num bad: 3
    Distribution of the training categories:
    0 : 1473 (67.59981642955483%)
    1 : 706 (32.40018357044516%)

    * train_pct 0.75
    Num correct: 417.0
    Num incorrect: 133.0
    Accuracy: 0.7581818181818182

    Confusion matrix:
    2×2 Matrix{Float64}:
    328.0  39.0
     94.0  89.0

    num bad: 133
    Distribution of the training categories:
    0 : 1123 (68.01938219261054%)
    1 : 528 (31.98061780738946%)

    * train pct 0.66
    Num correct: 586.0
    Num incorrect: 162.0
    Accuracy: 0.7834224598930482

    Confusion matrix:
    2×2 Matrix{Float64}:
    461.0   50.0
    112.0  125.0

    num bad: 162
    Distribution of the training categories:
    0 : 979 (67.37783895388851%)
    1 : 474 (32.62216104611149%)

    * train pct 0.33
    Num correct: 1128.0
    Num incorrect: 347.0
    Accuracy: 0.7647457627118645

    Confusion matrix:
    2×2 Matrix{Float64}:
    908.0   89.0
    258.0  220.0

    num bad: 347
    Distribution of the training categories:
    0 : 493 (67.90633608815428%)
    1 : 233 (32.09366391184573%)

    * train_pct 0.1
    Num correct: 1525.0
    Num incorrect: 456.0
    Accuracy: 0.7698132256436143

    Confusion matrix:
    2×2 Matrix{Float64}:
    1209.0  114.0
     342.0  316.0

    num bad: 456
    Distribution of the training categories:
    0 : 167 (75.9090909090909%)
    1 : 53 (24.09090909090909%)

    * train pct 0.004543389368468878 (10 training instances (of 2201)
    Num correct: 1679.0
    Num incorrect: 512.0
    Accuracy: 0.7663167503423094

    Confusion matrix:
    2×2 Matrix{Float64}:
    1370.0  114.0
     398.0  309.0

    num bad: 512
    Distribution of the training categories:
    0 : 6 (60.0%)
    1 : 4 (40.0%)

    * train pct 0.0009086778736937755 (2 training instances of 2201), no standardizing
    This is an extreme case with just 2 training instances (but I had to try this).
    Note that there are 67% change of just guessing the correct value, by always
    selecting the majority category.

    When there are two different categories there's often/sometimes a better chance (~ 76%)
    of getting the correct value, almost as good as with the full dataset...

    Num correct: 1632.0
    Num incorrect: 567.0
    Accuracy: 0.7421555252387448

    Confusion matrix:
    2×2 Matrix{Float64}:
    1274.0  215.0
     352.0  358.0

    num bad: 567
    Distribution of the training categories:
    0 : 1 (50.0%)
    1 : 1 (50.0%)

    But sometimes this fails. Here we have one training instance of each category, but
    got an accuracy 30%, which is worse than just guessing the majority category:

    Num correct: 677.0
    Num incorrect: 1522.0
    Accuracy: 0.30786721236925874

    Confusion matrix:
    2×2 Matrix{Float64}:
    398.0  1091.0
    431.0   279.0

    num bad: 1522
    Distribution of the training categories:
    0 : 1 (50.0%)
    1 : 1 (50.0%)

    Here is an example when the two training instance are from the same category:

    Num correct: 1488.0
    Num incorrect: 711.0
    Accuracy: 0.6766712141882674

    Confusion matrix:
    2×2 Matrix{Float64}:
    1488.0  0.0
     711.0  0.0

    num bad: 711
    Distribution of the training categories:
    0 : 2 (100.0%)

   Another example (with <30% accuracy)
    Num correct: 636.0
    Num incorrect: 1563.0
    Accuracy: 0.2892223738062756

    Confusion matrix:
    2×2 Matrix{Float64}:
     65.0  1425.0
    138.0   571.0

    num bad: 1563
    Distribution of the training categories:
    1 : 2 (100.0%)

=#
using Turing, StatsPlots, Distributions
using DataFrames, StatsBase
using CSV
using ARFFFiles
include("jl_utils.jl")

#
# Read dataset and preprocessing
#
df = DataFrame(ARFFFiles.load(download("http://www.hakank.org/weka/titanic.arff")))
display(describe(df))
println()

# Convert categories -> numeric values
transform!(df, [:class]    => ByRow(i -> Dict("1st"=>1,"2nd"=>2,"3rd"=>3,"crew"=>4)[i]) => :class)
transform!(df, [:age]      => ByRow(i -> Dict("child"=>1,"adult"=>2)[i]) => :age)
transform!(df, [:sex]      => ByRow(i -> Dict("male"=>1,"female"=>2)[i]) => :sex)
transform!(df, [:survived] => ByRow(i -> Dict("no"=>0,"yes"=>1)[i]) => :survived)

display(describe(df))
println()

# Convert to matrix and extract x and y.
num_rows,num_cols = size(df)
println("Num rows: $num_rows Num columns: $num_cols")
data = Matrix(df)
x = data[:,1:end-1]
y = data[:,end]

categories = unique(y) |> sort
println("\nMean value for each categories (full dataset): ", categories)
display(categories .|> s -> (s,mean(x[y.==s,:], dims=1)))
println("\n")

println("\nDistribution of the categories (full dataset):")
for (k,v) in make_hash(y)
  println("$k : $v ($(100*v/num_rows)%)")
end

# train_pct = 0.99
train_pct = 0.75
# train_pct = 0.66
# train_pct = 0.33
# train_pct = 0.1
# train_pct = 10/num_rows # 10 instances
# train_pct = 4/num_rows # 4 instances
# train_pct = 2/num_rows # 2 instances

println("\nTrain percentage: $train_pct")

x_train, x_test, y_train, y_test = split_train_test(x,y,train_pct)
num_train_rows = length(y_train)
num_test_rows = length(y_test)

println("Num training instances: $num_train_rows")
println("Num training instances: $num_test_rows")

println("Distribution of the training categories:")
for (k,v) in make_hash(y_train)
  println("$k : $v ($(100*v/num_train_rows)%)")
end

# Standardize the dataset (better and faster)
# Separate train and test set.
# if train_pct >= 0.2
#     means_train = mean(x_train,dims=1)
#     stds_train = std(x_train,dims=1)
#     x_train = (x_train .- means_train) ./ stds_train
#
#     means_test = mean(x_test,dims=1)
#     stds_test = std(x_test,dims=1)
#     x_train = (x_train .- means_test) ./ stds_test
# end

println("\nmean(x_train): ", mean(x_train,dims=1))
println("std(x_train): ", std(x_train,dims=1))
println()

# This is defined in jl_utils.jl
model = linear_regression_model(x_train,y_train)

println("\nModel:")
# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 1_000)
# chns = sample(model, PG(5), 1_000)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

# chns = sample(model, HMC(0.1,6), 1_000)
chns = sample(model, NUTS(), 1_000)

chns_orig = chns
display(chns)

# display(plot(chns))
println("\nPredictions:")

mpred = linear_regression_model(x_test,missing)
pred = predict(mpred,chns)
# display(pred)

cats = Dict(0=>1,1=>2)
confusion, bads = confusion_matrix(y_test,pred,cats,false)
println("\nConfusion matrix:")
display(confusion)
# println("\nBad classfied:")
# display(bads)
println("\nnum bad: $(length(bads))")

println("Distribution of the training categories:")
for (k,v) in make_hash(y_train)
  println("$k : $v ($(100*v/num_train_rows)%)")
end

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
using OdinSon

Polygon([1 2 3; 3 4 5])
render(PolyLine(Array[[1, 2], [3, 4], [5, 6]]))

diagram = Canvas([Circle([0, 0], 1), Circle([1, 1], 0.5, style=Style(fill=NC"violet", stroke_width=2))])
render(diagram)
xlim(-2, 2)
ylim(-2, 2)

diagram = Canvas([PolyLine(Array[[1, 2], [2, 3], [3, 2]], style=Style(stroke_width=5, stroke=HSV(180, 1, 1)))])
render(diagram)

p1 = Polygon(Array[[1, 2], [2, 3], [3, 2]])
render(Canvas([p1]))

# not sure what to do about nesting Canvas
diagram = Canvas([Circle([0, 0], 1), Canvas([Circle([1, 0.5], 0.1), Circle([0.1, 0.1], 0.01)]), Circle([0.1, 0.1], 0.3)])
render(diagram)

Circle((1, 1), 1)
Circle([1, 2], 3)

circ = Circle(1, 1, 1)
render(circ)

(fig, ax) = subplots(figsize=(5, 5))
ax[:add_patch](render(circ))
show(fig)

#Do I assume a default axes/coordinate system?
Canvas([Circle((x, y), r, style=Style(stroke=NC"red", stroke_width=2)), PolyLine()])

Style(axes=AxesStyle(stroke_width=0.8, stroke=NC"red"))

convert(Point, (1, 2))
convert(Point, [1, 2])
convert(Point, [1, 2, 3])

const _style_attrs = Dict{Symbol, Any}(
    :stroke => NC"black",
    :fill => :none,
    :stroke_width => 1.0)

function _check_attrs(attrs::Symbol)
    unknown = Symbol[]
    for attr in attrs
        if !haskey(_style_attrs, attr)
            push!(unknown, attr)
        end
    end
    return unknown
end

type Style2
    attrs::Dict{Symbol, Any}
end

function Style2(;kwargs...)
    kw = Dict{Symbol, Any}(kwargs)
    unknown = _check_attrs(keys(kw))
    if length(unknown) == 0
        return Style2(merge(_style_attrs, kw))
    else
        throw(ArgumentError("Unknown style attributes: $unknown"))
    end
end

Base.getindex(s::Style2, key::Symbol) = s.attrs[key]
function Base.setindex!(s::Style2, value, key::Symbol)
    if haskey(s.attrs, key)
        s.attrs[key] = value
        return value
    else
        throw(ArgumentError("Unknown style attribute: $key"))
    end
end

s = Style2(stroke_width=5)
s[:stroke_width]
s[:james] = NC"red"

# so pretty nice, but a lot of work and I don't have type checking on my arguments. The
# Type solution is looking nicer. I just want to add a merge method.
#
# Soooo. What to do about the Union{RGB, Void} issues. I like storing the RGB since
# this would mean I can use all the nice facilities of Colors, but due to the potential
# for nothing, I need to do a lot of checks, and lose type stability. I could go to
# Nullable{RGB}, but I am not sure that is the correct semantics, as it is not missing,
# as much as "not colored". I think I need a new type that is more descriptive. nothing
# is really for functions that don't return anything.
#
Viewport
GraphicsView # instead of Viewport for naming consistency a Viewport could be the abstract type
AxesView
GridView
# This will work from 0.5+ since there will no longer be implicit concat
GridView([[Viewport, AxesView], [AxesView]]) #(2 rows, 2 columns, second column spans)
# Is it worth having something that does a reshape? like:
GridView([v1, v2, v3, v4], layout=Layout([2, 2])) # rows and columns
Viewport([], inset=Inset(Viewport())
Viewport([], inset=Inset([Viewport(), AxesView()]))

Canvas(Viewport)

max_stars_repo_name: gabrielgellner/OdinSon

max_issues_repo_name: gabrielgellner/OdinSon

max_forks_repo_name: gabrielgellner/OdinSon


text:
### should refactor to RectangleGeometry with Identity as the default Stat.
struct RectangularBinGeometry <: Gadfly.GeometryElement
    default_statistic::Gadfly.StatisticElement
    tag::Symbol
end

function RectangularBinGeometry(
        default_statistic::Gadfly.StatisticElement=Gadfly.Stat.rectbin();
        tag=empty_tag)
    RectangularBinGeometry(default_statistic, tag)
end

"""
    Geom.rectbin

Draw equal sizes rectangles centered at `x` and `y` positions.  Optionally
specify their `color`.
"""
const rectbin = RectangularBinGeometry

"""
    Geom.rect

Draw colored rectangles with the corners specified by the
`xmin`, `xmax`, `ymin` and `ymax` aesthetics.  Optionally
specify their `color`.
"""
rect() = RectangularBinGeometry(Gadfly.Stat.Identity())

"""
    Geom.histogram2d[(; xbincount=nothing, xminbincount=3, xmaxbincount=150,
                        ybincount=nothing, yminbincount=3, ymaxbincount=150)]

Draw a heatmap of the `x` and `y` aesthetics by binning into rectangles and
indicating density with color.  This geometry is equivalent to
[`Geom.rect`](@ref) with [`Stat.histogram2d`](@ref);  see the latter for more
information.
"""
function histogram2d(; xbincount=nothing, xminbincount=3, xmaxbincount=150,
                       ybincount=nothing, yminbincount=3, ymaxbincount=150)
    RectangularBinGeometry(
        Gadfly.Stat.histogram2d(xbincount=xbincount,
                                xminbincount=xminbincount,
                                xmaxbincount=xmaxbincount,
                                ybincount=ybincount,
                                yminbincount=yminbincount,
                                ymaxbincount=ymaxbincount))
end

default_statistic(geom::RectangularBinGeometry) = geom.default_statistic

element_aesthetics(::RectangularBinGeometry) = [:xmin, :xmax, :ymin, :ymax, :color, :alpha]

# Render rectangular bin (e.g., heatmap) geometry.
#
# Args:
#   geom: rectbin geometry
#   theme: the plot's theme
#   aes: some aesthetics
#
# Returns
#   A compose form.
#
function render(geom::RectangularBinGeometry, theme::Gadfly.Theme, aes::Gadfly.Aesthetics)

    default_aes = Gadfly.Aesthetics()
    default_aes.color = RGBA{Float32}[theme.default_color]
    default_aes.alpha = Float64[theme.alphas[1]]
    aes = inherit(aes, default_aes)

    Gadfly.assert_aesthetics_defined("RectangularBinGeometry", aes, :xmin, :xmax, :ymin, :ymax)
    Gadfly.assert_aesthetics_equal_length("RectangularBinGeometry", aes, :xmin, :xmax, :ymin, :ymax)

    nx = length(aes.xmin)
    ny = length(aes.ymin)
    n = nx

    xmin = aes.xmin
    xwidths = [(x1 - x0)*cx - theme.bar_spacing
               for (x0, x1) in zip(aes.xmin, aes.xmax)]

    ymin = aes.ymin
    ywidths = [(y1 - y0)*cy - theme.bar_spacing
               for (y0, y1) in zip(aes.ymin, aes.ymax)]

    AT, CT = eltype(aes.alpha), eltype(aes.color)
    aes_color = Vector{CT}(undef, n)
    aes_alpha = Vector{Float64}(undef, n)
    alphav = AT <: Int ? theme.alphas[aes.alpha] : aes.alpha
    for (i, (_, c, a)) in enumerate(Compose.cyclezip(xmin, aes.color, alphav))
        aes_color[i] = c
        aes_alpha[i] = a
    end

    allvisible = true
    for c in aes_color
        if c == nothing
            allvisible = false
            break
        end
    end

    if !allvisible
        visibility = aes_color .!== nothing
        aes_color = aes_color[visibility]
        aes_alpha = aes_alpha[visibility]
        xmin = xmin[visibility]
        xwidths = xwidths[visibility]
        ymin = ymin[visibility]
        ywidths = ywidths[visibility]
    end

    return compose!(
        context(),
        rectangle(xmin, ymin, xwidths, ywidths, geom.tag),
        fill(aes_color), fillopacity(aes_alpha),
        stroke(nothing),
        svgclass("geometry"),
        svgattribute("shape-rendering", "crispEdges"))
end


text:
using Pkg
Pkg.activate(".")
using PlotlyJS, Zygote, LinearAlgebra
using Dash, DashCoreComponents, DashHtmlComponents
using DashBootstrapComponents

function frenet_frame(x, y, z)
    ẋ(t) = gradient(x, t)[1]
    ẏ(t) = gradient(y, t)[1]
    ż(t) = gradient(z, t)[1]

    ẍ(t) = hessian(x, t)
    ÿ(t) = hessian(y, t)
    z̈(t) = hessian(z, t)

    α̇(t) = [ẋ(t), ẏ(t), ż(t)]
    α̈(t) = [ẍ(t), ÿ(t), z̈(t)]
    function out(t)
        T = α̇(t)
        T = T / norm(T)
        ad = α̈(t)
        N = ad - dot(ad, T)T
        N = N / norm(N)
        B = cross(T, N)
        T, N, B
    end
end

function uv_surface(uv_domain, ϕ; uv_res = [40, 40])
    us = LinRange(uv_domain[1], uv_domain[2], uv_res[1])
    vs = LinRange(uv_domain[3], uv_domain[4], uv_res[2])
    u_res, v_res = uv_res
    xs = zeros(u_res, v_res)
    ys = zeros(u_res, v_res)
    zs = zeros(u_res, v_res)
    for (i, u) in enumerate(us)
        for (j, v) in enumerate(vs)
            xs[i, j], ys[i, j], zs[i, j] = ϕ(u, v)
        end
    end
    surface(x = xs, y = ys, z = zs, showscale = false)
end

function xyz(ps)
    x = [p[1] for p in ps]
    y = [p[2] for p in ps]
    z = [p[3] for p in ps]
    return x, y, z
end

function torus_knot(m::Int64, n::Int64; R = 2, r = 1, uv_res = [220, 220])
    ϕ(u, v) = [(R + r * cos(u)) * cos(v), (R + r * cos(u)) * sin(v), r * sin(u)]
    α(t) = ϕ(m * t, n * t)
    x(t) = α(t)[1]
    y(t) = α(t)[2]
    z(t) = α(t)[3]

    F = frenet_frame(x, y, z)
    function ψ(u, v)
        T, N, B = F(u)
        ρ = r / 3
        α(u) + ρ * cos(v) * N + ρ * sin(v) * B
    end
    uv_surface([-π, 1.02π, -π, π], ψ, uv_res = uv_res)
end
description = dcc_markdown(raw"Torus knots are obtained by winding a curve around a torus
 ($m$ times in parallel direction and $n$ times in the meridian direction). More concretely, if

 $$\varphi(u,v) = ((R+r\cos(u))\cos(v), (R+r\cos(u))\sin(v), r\sin(u))$$

 is the torus parametrization, the curve

 $$\alpha(t) = \varphi(mt,nt)\text{ for }0\leq t \leq 2\pi $$

is said to be knot of the type $(m,n)$. These knots can be made into surfaces by thickening them using the Frenet frame over $\alpha$
 ")
sd = 4.0
scene3d = Dict(
    :xaxis => Dict(:range => [-sd, sd], :autorange => false),
    :yaxis => Dict(:range => [-sd, sd], :autorange => false),
    :zaxis => Dict(:range => [-sd, sd], :autorange => false),
    :aspectratio => Dict(:x => 1, :y => 1, :z => 1),
)

layout = Layout(scene = scene3d, width = 600, height = 600)

mathjax = "https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.4/MathJax.js?config=TeX-MML-AM_CHTML"
app = dash(external_scripts = [mathjax], external_stylesheets = [dbc_themes.BOOTSTRAP])
app.title = "Dash Torus Knots"

knot = torus_knot(2, 3)
graph = dbc_spinner(dcc_graph(id = "graph", figure = plot(knot, layout)))

slider_m = dcc_slider(
    id = "slider_m",
    min = 1,
    max = 7,
    step = 1,
    value = 3,
    marks = Dict([i => ("m=$(i)") for i = 1:10]),
)
slider_n = dcc_slider(
    id = "slider_n",
    min = 1,
    max = 7,
    step = 1,
    value = 2,
    marks = Dict([i => ("n=$(i)") for i = 1:10]),
)

input = dbc_row([dbc_col(slider_m, width = Dict("size" => 3, "order" => "first", "offset" => 3)),
 dbc_col(slider_n, width = Dict("size" => 3, "order" => "last", "offset" => 0))])
app.layout = dbc_container(
    [
        html_h1(html_center("Torus Knots"))
        html_br()
        description
        html_center(graph)
        html_center(input)
        html_br()
        html_br()
    ],
)

callback!(
    app,
    Output("graph", "figure"),
    Input("slider_m", "value"),
    Input("slider_n", "value"),
) do m, n
    knot = torus_knot(m, n)
    return plot(knot, layout)
end
run_server(app, "0.0.0.0", debug = true)


text:
#!/usr/bin/env julia

using Printf

quad(a, squareRoot, b) = (-b + squareRoot) / 2a

# Solves x for a*x^2 + b*x + c = 0.
# Arguments types can be defined in function definitions
#
function QuadraticSolver(a::Float64, b::Float64, c::Float64)
    squareRoot = sqrt(b^2-4a*c)
    r1 = quad(a, squareRoot, b)
    r2 = quad(a, -squareRoot, b)

    r1, r2
end

# Solve for 2x^2 - 2x -12 = 0
#
result1, result2 = QuadraticSolver(2.0, -2.0, -12.0)
println("result 1: ", result1)
println("result 2: ", result2)

max_stars_repo_name: wbeebe/julia

max_issues_repo_name: wbeebe/julia

max_forks_repo_name: wbeebe/julia


text:
# This file includes functions that build subnetwork models
# with extensions from PowerModels.jl

function decompose(
    data::Dict{String, Any},
    N_gs::Vector{Vector{Int64}},
    modeltype::Type{T},
    build_function::Function;
    extra_ref_extensions=[]
    )::NetDecModel where T <: PM.AbstractPowerModel

    return decompose(data, [Set(i) for i in N_gs], modeltype, build_function, extra_ref_extensions = extra_ref_extensions)
end

function decompose(
    data::Dict{String, Any},
    N_gs::Vector{Set{Int64}},
    modeltype::Type{T},
    build_function::Function;
    extra_ref_extensions=[]
    )::NetDecModel where T <: PM.AbstractPowerModel

    models = T[]
    shared_vars_dict = Dict{Int64, Dict}()
    denet_refs = [ref_add_cut_bus!, ref_add_cut_branch!]
    if T <: PM.AbstractWRMModel
        push!(denet_refs, ref_add_global_bus!)
    end
    for i in eachindex(N_gs)
        N_g = N_gs[i]
        sub_data = generate_subnet_data(data, N_g)
        push!(models, PM.instantiate_model(sub_data, modeltype, build_function, ref_extensions = push!(denet_refs, extra_ref_extensions...)))
        shared_vars_dict[i] = collect_split_vars(models[i])
    end

    return NetDecModel(N_gs, models, shared_vars_dict)
end

# modified from original build_opf
# get rid of ref buses and apply power balance to only nodes in partition
function build_opf_mod(pm::PM.AbstractPowerModel)
    # variable_bus_voltage(pm, bounded=false)
    # vr = var(pm, :vr)
    # vi = var(pm, :vi)
    # for (i, bus) in setdiff(ref(pm, :bus), ref(pm, :cut_bus))
    #     JuMP.set_lower_bound(vr[i], -bus["vmax"])
    #     JuMP.set_upper_bound(vr[i],  bus["vmax"])
    #     JuMP.set_lower_bound(vi[i], -bus["vmax"])
    #     JuMP.set_upper_bound(vi[i],  bus["vmax"])
    #     constraint_voltage_magnitude_bounds(pm, i)
    # end
    PM.variable_bus_voltage(pm)
    PM.variable_gen_power(pm)
    PM.variable_branch_power(pm)
    PM.variable_dcline_power(pm)

    objective_min_fuel_and_flow_cost_mod(pm)

    PM.constraint_model_voltage(pm)

    for i in setdiff(PM.ids(pm, :bus), PM.ids(pm, :cut_bus))
        PM.constraint_power_balance(pm, i)
    end

    for i in PM.ids(pm, :branch)
        PM.constraint_ohms_yt_from(pm, i)
        PM.constraint_ohms_yt_to(pm, i)

        PM.constraint_voltage_angle_difference(pm, i)

        PM.constraint_thermal_limit_from(pm, i)
        PM.constraint_thermal_limit_to(pm, i)
    end

    for i in PM.ids(pm, :dcline)
        PM.constraint_dcline_power_losses(pm, i)
    end
end

function build_opf_mod(pm::PM.AbstractWRMModel)
    function variable_bus_voltage(pm::PM.AbstractWRMModel; nw::Int=PM.nw_id_default, bounded::Bool=true, report::Bool=true)
        wr_min, wr_max, wi_min, wi_max = PM.ref_calc_voltage_product_bounds(PM.ref(pm, nw, :all_buspairs))
        bus_ids = PM.ids(pm, nw, :all_bus)

        w_index = 1:length(bus_ids)
        lookup_w_index = Dict((bi,i) for (i,bi) in enumerate(bus_ids))

        WR_start = zeros(length(bus_ids), length(bus_ids)) + I

        WR = PM.var(pm, nw)[:WR] = JuMP.@variable(pm.model,
            [i=1:length(bus_ids), j=1:length(bus_ids)], Symmetric, base_name="$(nw)_WR", start=WR_start[i,j]
        )
        if report
            PM.sol(pm, nw)[:WR] = WR
        end

        WI = PM.var(pm, nw)[:WI] = JuMP.@variable(pm.model,
            [1:length(bus_ids), 1:length(bus_ids)], base_name="$(nw)_WI", start=0.0
        )
        if report
            PM.sol(pm, nw)[:WI] = WI
        end

        # bounds on diagonal
        for (i, bus) in PM.ref(pm, nw, :all_bus)
            w_idx = lookup_w_index[i]
            wr_ii = WR[w_idx,w_idx]
            wi_ii = WR[w_idx,w_idx]

            if bounded
                JuMP.set_lower_bound(wr_ii, (bus["vmin"])^2)
                JuMP.set_upper_bound(wr_ii, (bus["vmax"])^2)
            else
                JuMP.set_lower_bound(wr_ii, 0)
            end
        end

        # bounds on off-diagonal
        for (i,j) in PM.ids(pm, nw, :all_buspairs)
            wi_idx = lookup_w_index[i]
            wj_idx = lookup_w_index[j]

            if bounded
                JuMP.set_upper_bound(WR[wi_idx, wj_idx], wr_max[(i,j)])
                JuMP.set_lower_bound(WR[wi_idx, wj_idx], wr_min[(i,j)])

                JuMP.set_upper_bound(WI[wi_idx, wj_idx], wi_max[(i,j)])
                JuMP.set_lower_bound(WI[wi_idx, wj_idx], wi_min[(i,j)])
            end
        end

        PM.var(pm, nw)[:w] = Dict{Int,Any}()
        for (i, bus) in PM.ref(pm, nw, :all_bus)
            w_idx = lookup_w_index[i]
            PM.var(pm, nw, :w)[i] = WR[w_idx,w_idx]
        end
        report && PM.sol_component_value(pm, nw, :all_bus, :w, PM.ids(pm, nw, :all_bus), PM.var(pm, nw)[:w])

        PM.var(pm, nw)[:wr] = Dict{Tuple{Int,Int},Any}()
        PM.var(pm, nw)[:wi] = Dict{Tuple{Int,Int},Any}()
        for (i,j) in PM.ids(pm, nw, :all_buspairs)
            w_fr_index = lookup_w_index[i]
            w_to_index = lookup_w_index[j]

            PM.var(pm, nw, :wr)[(i,j)] = WR[w_fr_index, w_to_index]
            PM.var(pm, nw, :wi)[(i,j)] = WI[w_fr_index, w_to_index]
        end
    end
    function constraint_voltage_angle_difference_all(pm::PM.AbstractPowerModel, i::Int; nw::Int=PM.nw_id_default)
        branch = PM.ref(pm, nw, :all_branch, i)
        f_bus = branch["f_bus"]
        t_bus = branch["t_bus"]
        f_idx = (i, f_bus, t_bus)
        pair = (f_bus, t_bus)
        buspair = PM.ref(pm, nw, :all_buspairs, pair)

        if buspair["branch"] == i
            PM.constraint_voltage_angle_difference(pm, nw, f_idx, buspair["angmin"], buspair["angmax"])
        end
    end
    variable_bus_voltage(pm) # re-define W matrix
    PM.variable_gen_power(pm)
    PM.variable_branch_power(pm)
    PM.variable_dcline_power(pm)

    objective_min_fuel_and_flow_cost_mod(pm) # allow for no generator case

    PM.constraint_model_voltage(pm)

    for i in setdiff(PM.ids(pm, :bus), PM.ids(pm, :cut_bus))
        PM.constraint_power_balance(pm, i)
    end

    for i in PM.ids(pm, :branch)
        PM.constraint_ohms_yt_from(pm, i)
        PM.constraint_ohms_yt_to(pm, i)

        # PM.constraint_voltage_angle_difference(pm, i)

        PM.constraint_thermal_limit_from(pm, i)
        PM.constraint_thermal_limit_to(pm, i)
    end

    for i in PM.ids(pm, :all_branch)
        constraint_voltage_angle_difference_all(pm, i)
    end

    for i in PM.ids(pm, :dcline)
        PM.constraint_dcline_power_losses(pm, i)
    end
end

#=
function PM.variable_bus_voltage_on_off(pm::PM.AbstractACRModel; kwargs...)
    variable_bus_voltage(pm; kwargs...)
end

function PM.constraint_model_voltage_on_off(pm::PM.AbstractACRModel; kwargs...)
end

function PM.constraint_ohms_yt_from_on_off(pm::AbstractACRModel, n::Int, f_bus, t_bus, f_idx, t_idx, g, b, g_fr, b_fr, tr, ti, tm)
    p_fr = var(pm, n, :p, f_idx)
    q_fr = var(pm, n, :q, f_idx)
    wrr_ff = var(pm, n, :wrr)[f_bus, f_bus]
    wii_ff = var(pm, n, :wii)[f_bus, f_bus]
    wrr_ft = var(pm, n, :wrr)[f_bus, t_bus]
    wii_ft = var(pm, n, :wii)[f_bus, t_bus]
    wri_ft = var(pm, n, :wri)[f_bus, t_bus]
    wri_tf = var(pm, n, :wri)[t_bus, f_bus]
    z = var(pm, n, :z_branch, i)

    JuMP.@constraint(pm.model, p_fr == z * ((g+g_fr)/tm^2*(wrr_ff + wii_ff) + (-g*tr+b*ti)/tm^2*(wrr_ft + wii_ft) + (-b*tr-g*ti)/tm^2*(wri_tf - wri_ft)) )
    JuMP.@constraint(pm.model, q_fr == z * (-(b+b_fr)/tm^2*(wrr_ff + wii_ff) - (-b*tr-g*ti)/tm^2*(wrr_ft + wii_ft) + (-g*tr+b*ti)/tm^2*(wri_tf - wri_ft)) )
end
=#

# modified from original build_ots
# get rid of ref buses and apply power balance to only nodes in partition
function build_ots_mod(pm::PM.AbstractPowerModel)
    PM.variable_branch_indicator(pm)
    PM.variable_bus_voltage_on_off(pm)
    PM.variable_gen_power(pm)
    PM.variable_branch_power(pm)
    PM.variable_dcline_power(pm)

    objective_min_fuel_and_flow_cost_mod(pm)

    PM.constraint_model_voltage_on_off(pm)

    for i in setdiff(PM.ids(pm, :bus), PM.ids(pm, :cut_bus))
        PM.constraint_power_balance(pm, i)
    end

    for i in PM.ids(pm, :branch)
        PM.constraint_ohms_yt_from_on_off(pm, i)
        PM.constraint_ohms_yt_to_on_off(pm, i)

        PM.constraint_voltage_angle_difference_on_off(pm, i)

        PM.constraint_thermal_limit_from_on_off(pm, i)
        PM.constraint_thermal_limit_to_on_off(pm, i)
    end

    for i in PM.ids(pm, :dcline)
        PM.constraint_dcline_power_losses(pm, i)
    end
end

# modified from original build_opf_bf
# get rid of ref buses and apply power balance to only nodes in partition
function build_opf_bf_mod(pm::PM.AbstractPowerModel)
    PM.variable_bus_voltage(pm)
    PM.variable_gen_power(pm)
    PM.variable_branch_power(pm)
    PM.variable_branch_current(pm)
    PM.variable_dcline_power(pm)

    objective_min_fuel_and_flow_cost_mod(pm)

    PM.constraint_model_current(pm)

    for i in setdiff(PM.ids(pm, :bus), PM.ids(pm, :cut_bus))
        PM.constraint_power_balance(pm, i)
    end

    for i in PM.ids(pm, :branch)
        PM.constraint_power_losses(pm, i)
        PM.constraint_voltage_magnitude_difference(pm, i)

        PM.constraint_voltage_angle_difference(pm, i)

        PM.constraint_thermal_limit_from(pm, i)
        PM.constraint_thermal_limit_to(pm, i)
    end

    for i in PM.ids(pm, :dcline)
        PM.constraint_dcline_power_losses(pm, i)
    end
end

# modified from original objective_min_fuel_and_flow_cost to allow for
# the case with no generators
function objective_min_fuel_and_flow_cost_mod(pm::PM.AbstractPowerModel; kwargs...)
    model = PM.check_cost_models(pm)
    if model == 1
        return PM.objective_min_fuel_and_flow_cost_pwl(pm; kwargs...)
    elseif model == 2
        return PM.objective_min_fuel_and_flow_cost_polynomial(pm; kwargs...)
    elseif model === nothing
        return JuMP.@objective(pm.model, Min, 0)
    else
        Memento.error(_LOGGER, "Only cost models of types 1 and 2 are supported at this time, given cost model type of $(model)")
    end
end

function variable_w_matrix(pm::PM.ACRPowerModel; nw::Int=PM.nw_id_default)
    wrr = PM.var(pm, nw)[:wrr] = JuMP.@variable(pm.model,
        [i in PM.ids(pm, nw, :bus), j in PM.ids(pm, nw, :bus)], base_name="$(nw)_wrr",
        start = PM.comp_start_value(PM.ref(pm, nw, :bus, i), "vr_start", 1.0) * PM.comp_start_value(PM.ref(pm, nw, :bus, j), "vr_start", 1.0)
        )
    wri = PM.var(pm, nw)[:wri] = JuMP.@variable(pm.model,
        [i in PM.ids(pm, nw, :bus), j in PM.ids(pm, nw, :bus)], base_name="$(nw)_wri",
        start = PM.comp_start_value(PM.ref(pm, nw, :bus, i), "vr_start", 1.0) * PM.comp_start_value(PM.ref(pm, nw, :bus, j), "vi_start")
        )
    wii = PM.var(pm, nw)[:wii] = JuMP.@variable(pm.model,
        [i in PM.ids(pm, nw, :bus), j in PM.ids(pm, nw, :bus)], base_name="$(nw)_wii",
        start = PM.comp_start_value(PM.ref(pm, nw, :bus, i), "vi_start") * PM.comp_start_value(PM.ref(pm, nw, :bus, j), "vi_start")
        )
end

#=
function build_acopf_with_free_lines(pm::PM.AbstractPowerModel)
    # this is for adding w variables to the model
    build_opf_mod(pm)
    PM.variable_bus_voltage_magnitude_sqr(pm)
    PM.variable_buspair_voltage_product(pm)

    w  = PM.var(pm,  :w)
    wr = PM.var(pm, :wr)
    wi = PM.var(pm, :wi)
    vr = PM.var(pm, :vr)
    vi = PM.var(pm, :vi)

    for (i, bus) in PM.ref(pm, :bus)
        JuMP.@constraint(pm.model, w[i] == vr[i]^2 + vi[i]^2)
    end

    for (_, branch) in PM.ref(pm, :branch)
        fbus = branch["f_bus"]
        tbus = branch["t_bus"]
        JuMP.@constraint(pm.model, wr[(fbus, tbus)] == vr[fbus] * vr[tbus] + vi[fbus] * vi[tbus])
        JuMP.@constraint(pm.model, wi[(fbus, tbus)] == vi[fbus] * vr[tbus] - vr[fbus] * vi[tbus])
    end
end
=#

function build_acopf_with_free_lines(pm::PM.AbstractPowerModel)
    build_opf_mod(pm)
    variable_w_matrix(pm)
    wrr = PM.var(pm, :wrr)
    wri = PM.var(pm, :wri)
    wii = PM.var(pm, :wii)
    vr = PM.var(pm, :vr)
    vi = PM.var(pm, :vi)
    for (i, _) in PM.ref(pm, :bus), (j, _) in PM.ref(pm, :bus)
        JuMP.@constraint(pm.model, wrr[i,j] == vr[i] * vr[j])
        JuMP.@constraint(pm.model, wri[i,j] == vr[i] * vi[j])
        JuMP.@constraint(pm.model, wii[i,j] == vi[i] * vi[j])
    end
end

function build_acopf_with_free_lines(pm::PM.SDPWRMPowerModel)
    build_opf_mod(pm)
end

function build_socwr_with_free_lines(pm::PM.AbstractWRModel)
    build_opf_mod(pm)
end

function build_socbf_with_free_lines(pm::PM.AbstractBFModel)
    build_opf_bf_mod(pm)
end

function build_acots_with_free_lines(pm::PM.AbstractACPModel)
    build_ots_mod(pm)
end

# function build_acots_with_free_lines(pm::PM.AbstractACRModel)
#     variable_w_matrix(pm)
#     build_ots_mod(pm)
# end

function collect_split_vars(pm::PM.ACPPowerModel)
    p = PM.var(pm, :p)
    q = PM.var(pm, :q)
    z = PM.var(pm, :z_branch)
    shared_vars_dict = Dict{String, Dict{Tuple, JuMP.VariableRef}}()
    shared_vars_dict["z"] = Dict{Tuple{Int64, Int64, Int64}, JuMP.VariableRef}()
    shared_vars_dict["p"] = Dict{Tuple{Int64, Int64, Int64}, JuMP.VariableRef}()
    shared_vars_dict["q"] = Dict{Tuple{Int64, Int64, Int64}, JuMP.VariableRef}()
    cut_arcs_from = PM.ref(pm, :cut_arcs_from)
    for (l,i,j) in cut_arcs_from
        shared_vars_dict["z"][(l,i,j)] = z[l]
        shared_vars_dict["p"][(l,i,j)] = p[(l,i,j)]
        shared_vars_dict["p"][(l,j,i)] = p[(l,j,i)]
        shared_vars_dict["q"][(l,i,j)] = q[(l,i,j)]
        shared_vars_dict["q"][(l,j,i)] = q[(l,j,i)]
    end
    return shared_vars_dict
end

function collect_split_vars(pm::PM.ACRPowerModel)
    wrr = PM.var(pm, :wrr)
    wri = PM.var(pm, :wri)
    wii = PM.var(pm, :wii)
    p  = PM.var(pm,  :p)
    q  = PM.var(pm,  :q)

    shared_vars_dict = Dict{String, Dict{Tuple, JuMP.VariableRef}}()
    shared_vars_dict["wrr"] = Dict{Tuple{Int64, Int64}, JuMP.VariableRef}()
    shared_vars_dict["wri"] = Dict{Tuple{Int64, Int64}, JuMP.VariableRef}()
    shared_vars_dict["wir"] = Dict{Tuple{Int64, Int64}, JuMP.VariableRef}()
    shared_vars_dict["wii"] = Dict{Tuple{Int64, Int64}, JuMP.VariableRef}()
    shared_vars_dict["p"] = Dict{Tuple{Int64, Int64, Int64}, JuMP.VariableRef}()
    shared_vars_dict["q"] = Dict{Tuple{Int64, Int64, Int64}, JuMP.VariableRef}()

    cut_arcs_from = PM.ref(pm, :cut_arcs_from)
    for (l,i,j) in cut_arcs_from
        if !((i,j) in keys(shared_vars_dict["wrr"]))
            shared_vars_dict["wrr"][(i,j)] = wrr[i,j]
            shared_vars_dict["wri"][(i,j)] = wri[i,j]
            shared_vars_dict["wir"][(i,j)] = wri[j,i]
            shared_vars_dict["wii"][(i,j)] = wii[i,j]
        end
        shared_vars_dict["p"][(l,i,j)] = p[(l,i,j)]
        shared_vars_dict["p"][(l,j,i)] = p[(l,j,i)]
        shared_vars_dict["q"][(l,i,j)] = q[(l,i,j)]
        shared_vars_dict["q"][(l,j,i)] = q[(l,j,i)]
    end
    return shared_vars_dict
end

function collect_split_vars(pm::PM.AbstractWRModel)
    wr = PM.var(pm, :wr)
    wi = PM.var(pm, :wi)
    p  = PM.var(pm,  :p)
    q  = PM.var(pm,  :q)

    shared_vars_dict = Dict{String, Dict{Tuple, JuMP.VariableRef}}()
    # shared_vars_dict["wr"] = Dict{Tuple{Int64, Int64}, JuMP.VariableRef}()
    # shared_vars_dict["wi"] = Dict{Tuple{Int64, Int64}, JuMP.VariableRef}()
    shared_vars_dict["p"] = Dict{Tuple{Int64, Int64, Int64}, JuMP.VariableRef}()
    shared_vars_dict["q"] = Dict{Tuple{Int64, Int64, Int64}, JuMP.VariableRef}()

    cut_arcs_from = PM.ref(pm, :cut_arcs_from)
    for (l,i,j) in cut_arcs_from
        # if !((i,j) in keys(shared_vars_dict["wr"]))
        #     shared_vars_dict["wr"][(i,j)] = wr[(i,j)]
        #     shared_vars_dict["wi"][(i,j)] = wi[(i,j)]
        # end
        shared_vars_dict["p"][(l,i,j)] = p[(l,i,j)]
        shared_vars_dict["p"][(l,j,i)] = p[(l,j,i)]
        shared_vars_dict["q"][(l,i,j)] = q[(l,i,j)]
        shared_vars_dict["q"][(l,j,i)] = q[(l,j,i)]
    end
    return shared_vars_dict
end

function collect_split_vars(pm::PM.AbstractBFModel)
    p  = PM.var(pm,  :p)
    q  = PM.var(pm,  :q)
    ccm = PM.var(pm, :ccm)

    shared_vars_dict = Dict{String, Dict}()
    shared_vars_dict["ccm"] = Dict{Int64, JuMP.VariableRef}()
    shared_vars_dict["p"] = Dict{Tuple{Int64, Int64, Int64}, JuMP.VariableRef}()
    shared_vars_dict["q"] = Dict{Tuple{Int64, Int64, Int64}, JuMP.VariableRef}()

    cut_arcs_from = PM.ref(pm, :cut_arcs_from)
    for (l,i,j) in cut_arcs_from
        shared_vars_dict["ccm"][l] = ccm[l]
        shared_vars_dict["p"][(l,i,j)] = p[(l,i,j)]
        shared_vars_dict["p"][(l,j,i)] = p[(l,j,i)]
        shared_vars_dict["q"][(l,i,j)] = q[(l,i,j)]
        shared_vars_dict["q"][(l,j,i)] = q[(l,j,i)]
    end
    return shared_vars_dict
end

function collect_split_vars(pm::PM.AbstractWRMModel)
    p = PM.var(pm, :p)
    q = PM.var(pm, :q)
    WR = PM.var(pm, :WR)
    WI = PM.var(pm, :WI)
    bus_ids = PM.ids(pm, :all_bus) # all buses from the full network
    lookup_w_index = Dict((bi,i) for (i,bi) in enumerate(bus_ids))

    shared_vars_dict = Dict{String, Dict{Tuple, JuMP.VariableRef}}()
    shared_vars_dict["p"] = Dict{Tuple{Int64, Int64, Int64}, JuMP.VariableRef}()
    shared_vars_dict["q"] = Dict{Tuple{Int64, Int64, Int64}, JuMP.VariableRef}()
    shared_vars_dict["WR"] = Dict{Tuple{Int64, Int64}, JuMP.VariableRef}()
    shared_vars_dict["WI"] = Dict{Tuple{Int64, Int64}, JuMP.VariableRef}()
    cut_arcs_from = PM.ref(pm, :cut_arcs_from)
    for (l,i,j) in cut_arcs_from # penalizing cut lines
        shared_vars_dict["p"][(l,i,j)] = p[(l,i,j)]
        shared_vars_dict["p"][(l,j,i)] = p[(l,j,i)]
        shared_vars_dict["q"][(l,i,j)] = q[(l,i,j)]
        shared_vars_dict["q"][(l,j,i)] = q[(l,j,i)]
    end
    bus = PM.ids(pm, :bus)
    for i in bus_ids, j in bus_ids
        shared_vars_dict["WR"][(i,j)] = WR[lookup_w_index[i], lookup_w_index[j]]
        shared_vars_dict["WI"][(i,j)] = WI[lookup_w_index[i], lookup_w_index[j]]
        if i != j
            shared_vars_dict["WR"][(j,i)] = WR[lookup_w_index[j], lookup_w_index[i]]
            shared_vars_dict["WI"][(j,i)] = WI[lookup_w_index[j], lookup_w_index[i]]
        end
    end
    return shared_vars_dict
end

######################################################
###             EXPERIMENTAL FUNCTIONS             ###
######################################################

function constraint_model_voltage_mod(pm::PM.AbstractPowerModel)
    PM.constraint_model_voltage(pm)
end

# function constraint_model_voltage_mod(pm::PM.AbstractWRMModel)
#     PM._check_missing_keys(PM.var(pm), [:WR,:WI], typeof(pm))

#     WR = PM.var(pm)[:WR]
#     WI = PM.var(pm)[:WI]

#     bus_ids = PM.ids(pm, :bus)
#     lookup_w_index = Dict((bi,i) for (i,bi) in enumerate(bus_ids))

#     local_buses = setdiff(PM.ids(pm, :bus), PM.ids(pm, :cut_bus))
#     local_w_idx = [lookup_w_index[i] for i in local_buses]

#     WR_local = WR[local_w_idx, local_w_idx]
#     WI_local = WI[local_w_idx, local_w_idx]
#     JuMP.@constraint(pm.model, [WR_local WI_local; -WI_local WR_local] in JuMP.PSDCone())
# end


text:
abstract type AbstractWindShearModel end

struct PowerLawWindShear{TF} <: AbstractWindShearModel

    # model parameter
    shear_exponent::TF

end

#TODO add log shear

function adjust_for_wind_shear(locz, reference_velocity, reference_height, ground_height, model::PowerLawWindShear)

    # initialize adjusted wind speed to zero
    adjusted_wind_speed = 0.0
    shear_exp = model.shear_exponent

    # check that the point of interest is above ground level
    if locz >= ground_height
        # adjusted wind speed for wind shear if point is above ground
        adjusted_wind_speed = reference_velocity*((locz-ground_height)/(reference_height-ground_height))^shear_exp
    end

    return adjusted_wind_speed
end

# function adjust_for_wind_shear(loc, point_velocity_no_shear, reference_height, ground_height, shear_model::CubicWindShear)

#     # initialize adjusted wind speed to zero
#     adjusted_wind_speed = 0.0

#     # check that the point of interest is above ground level
#     if loc[3] >= ground_height
#         # adjusted wind speed for wind shear if point is above ground
#         adjusted_wind_speed = point_velocity_no_shear*((loc[3]-ground_height)/(reference_height-ground_height))^shear_exp
#     else
#         # if the point of interest is below ground, set the wind speed to 0.0
#         adjusted_wind_speed = 0.0
#     end

#     return adjusted_wind_speed
# end


text:
using Escher
using Markdown
using Gadfly
using Distributions
using Reactive
using SerialPorts

function main(window)
    push!(window.assets, "layout2")
    push!(window.assets, "icons")
    push!(window.assets, "widgets")

    serial = serialport("/dev/ttyACM0", 250000)

    try
        serial = open(serial)
    catch err
        println("caught error, serial might be already open")
        println(err)
    end

    desc = md"""
    # Hello, World!
    You can write any **Markdown** text and place it in the middle of other
    tiles.
    """ |> pad(1em)

    αᵗ = Input(1.0)
    βᵗ = Input(1.0)
    z_axis_inp = Input(1.0)

    plot_beta(α, β) =
        plot(x -> pdf(Beta(α, β), x), 0, 1)

    t = tabs([hbox(icon("home")),
                      hbox(icon("info-outline")),
                      hbox(icon("settings")),
                      hbox(iconbutton("search"))])

    lift(αᵗ, βᵗ, z_axis_inp) do α, β, z_axis
        ui = vbox(
            maxwidth(3, t),
            vbox(md"## Dynamic plot",
                hbox(vbox(
                    hbox("Alpha: " |>
                        width(4em), slider(1:100) >>> αᵗ) |>
                        packacross(center),
                    hbox("Beta: "  |>
                        width(4em), slider(1:100) >>> βᵗ) |>
                        packacross(center)),
                    plot_beta(α,β) |> drawing((4*2/3)inch, 2inch)),
                  md"## Axes",
                  hbox("Z: " |> width(4em), slider(0:0.1:100) >>> z_axis_inp)
                    |> packacross(center)
            ) |> pad(2em))
        move_cmd = string("G1 Z",z_axis," F10000\n")
        write(serial, move_cmd)
        @show move_cmd
        ui
    end
end


text:
# module Bukdu

const env = Dict{Symbol, Any}(
    :server => nothing,
    :check_websocket => false,
    :check_server_sent_events => false,
    :sse_streams => Dict{RawFD,Any}(),
    :websockets => Dict{RawFD,Any}(),
)

function sse_streams()
    streams = Dict{RawFD,Any}()
    dirty = false
    for (key, http) in pairs(env[:sse_streams])
         s = Base.uv_status_string(http.stream.c.io)
         if s == "paused"
             streams[key] = http
         else
             dirty = true
         end
    end
    if dirty
        env[:sse_streams] = streams
    end
    values(streams)
end

function websockets()
    sockets = Dict{RawFD,Any}()
    dirty = false
    for (key, ws) in pairs(env[:websockets])
         s = Base.uv_status_string(ws.io)
         if s == "active"
             sockets[key] = ws
         else
             dirty = true
         end
    end
    if dirty
        env[:websockets] = sockets
    end
    values(sockets)
end

function _base_routing_handle(http::Deps.HTTP.Stream) # result
    request::Deps.Request = http.message
    request.body = read(http)
    route = Routing.handle(request)
    request_handler(route, request) # result
end

function routing_handle(http::Deps.HTTP.Stream)::Bool # needs_to_close
    result = _base_routing_handle(http)
    startwrite(http)                   #
    write(http, result.resp.body) #
    true
end

function full_routing_handle(http::Deps.HTTP.Stream)::Bool # needs_to_close
    if Deps.HTTP.WebSockets.is_upgrade(http.message)
        Deps.HTTP.WebSockets.upgrade(http) do ws::Deps.HTTP.WebSockets.WebSocket
            fd = Base._fd(Deps.HTTP.IOExtras.tcpsocket(ws.io))
            env[:websockets][fd] = ws
            while !eof(ws)
                data = readavailable(ws)
                write(ws,data)
            end
        end
        true
    else
        result = _base_routing_handle(http)
        if result.got isa EventStream
            fd = Base._fd(Deps.HTTP.IOExtras.tcpsocket(http.stream.c.io))
            env[:sse_streams][fd] = http
            needs_to_close = false
        else
            startwrite(http)              #
            write(http, result.resp.body) #
            needs_to_close = true
        end
        needs_to_close
    end
end

import Sockets: @ip_str
function start(port; host=ip"127.0.0.1")
    server = Deps.HTTP.Servers.Server(routing_handle, stdout)
    env[:server] = server
    @async serve(
                server,
                host,
                port,
                false) # verbose
end

function stop()
    server = env[:server]
    env[:server] = nothing
    env[:check_websocket] = false
    env[:check_server_sent_events] = false
    env[:sse_streams] = Dict{RawFD,Any}()
    env[:websockets] = Dict{RawFD,Any}()
    server isa Deps.HTTP.Servers.Server && put!(server.in, Deps.HTTP.Servers.KILL)
    nothing
end

print_listening_on(host, port) = @info "Bukdu Listening on: $host:$port"

#########################################################################
#
#  customizing HTTP.jl
#
#  original code from https://github.com/JuliaWeb/HTTP.jl/tree/master/src

import .Deps: HTTP
import .HTTP.Servers: Stream, ConnectionPool, iswritable, hasheader, setheader, writeheaders
import .HTTP.Streams: messagetowrite
import .HTTP.IOExtras: startwrite
import .HTTP.Servers: Server, https, RateLimit, Transaction, Streams, Stream, KILL, nosslconfig, Connection, nolimit, startread, closeread, closewrite, hasheader, isioerror, update!
import Sockets # Sockets.TCPServer
import Sockets: accept, IPAddr

# HTTP.jl - Servers.jl
# doesn't need to be closed when EventStream has set
function _handle_stream(f::Function, http::Stream)::Nothing
    needs_to_close = true
    try
        needs_to_close = f(http) # routing_handle full_routing_handle
    catch e
        if isopen(http) && !iswritable(http)
            @error Symbol(:server_, :_handle_stream) e stacktrace(catch_backtrace())
            http.message.response.status = 500
            startwrite(http)
            write(http, sprint(showerror, e))
        else
            rethrow(e)
        end
    end
    if needs_to_close # SSE
        closeread(http)
        closewrite(http)
    end
    nothing
end

import Dates: now
# HTTP.jl - Servers.jl
function _check_rate_limit(tcp;
                          ratelimits=nothing,
                          ratelimit::Rational{Int}=Int(10)//Int(1), kw...)
    ip = Sockets.getsockname(tcp)[1]
    rate = Float64(ratelimit.num)
    rl = get!(ratelimits, ip, RateLimit(rate, now()))
    update!(rl, ratelimit)
    if rl.allowance > rate
        # @warn Symbol(:server_, :_check_rate_limit) "throttling $ip" ## bukdu
        rl.allowance = rate
    end
    if rl.allowance < 1.0
        @warn Symbol(:server_, :_check_rate_limit) "discarding connection from $ip due to rate limiting"
        return false
    else
        rl.allowance -= 1.0
    end
    return true
end

# HTTP.jl - Servers.jl
function serve(server::Server{T, H}, host, port, verbose) where {T, H}

    tcpserver = Ref{Sockets.TCPServer}()

    @async begin
        while !isassigned(tcpserver)
            sleep(1)
        end
        while true
            val = take!(server.in)
            val == KILL && close(tcpserver[])
        end
    end

    f = env[:check_websocket] || env[:check_server_sent_events] ? full_routing_handle : routing_handle ## bukdu
    listen(f, host, port;
           tcpref=tcpserver,
           ssl=(T == https),
           sslconfig=server.options.sslconfig,
           verbose=verbose,
           tcpisvalid=server.options.ratelimit > 0 ? _check_rate_limit :
                                                     (tcp; kw...) -> true,
           ratelimits=Dict{IPAddr, RateLimit}(),
           ratelimit=server.options.ratelimit)

    return
end

# HTTP.jl - Servers.jl
listen(f, host, port; kw...) = listen(f, string(host), Int(port); kw...)
function listen(f::Function,
                host::String="127.0.0.1", port::Int=8081;
                ssl::Bool=false,
                require_ssl_verification::Bool=true,
                sslconfig::Deps.MbedTLS.SSLConfig=nosslconfig,
                pipeline_limit::Int=ConnectionPool.default_pipeline_limit,
                tcpisvalid::Function=(tcp; kw...)->true,
                tcpref::Ref{Sockets.TCPServer}=Ref{Sockets.TCPServer}(),
                kw...)

    if sslconfig === nosslconfig
        sslconfig = Deps.MbedTLS.SSLConfig(require_ssl_verification)
    end

    print_listening_on(host, port) ##
    tcpserver = Sockets.listen(Sockets.getaddrinfo(host), port)

    tcpref[] = tcpserver

    try
        while isopen(tcpserver)
            try
                io = accept(tcpserver)
            catch e
                if e isa Base.UVError
                    @warn Symbol(:server_, :listen) "$e"
                    break
                else
                    rethrow(e)
                end
            end
            if !tcpisvalid(io; kw...)
                close(io)
                continue
            end
            io = ssl ? getsslcontext(io, sslconfig) : io
            let io = Connection(host, string(port), pipeline_limit, 0, io)
                # @info "Accept:  $io" ##
                @async try
                    _handle_connection(f, io; kw...)
                catch e
                    @error Symbol(:server_, :listen) "Error:   $io" e stacktrace(catch_backtrace()) ##
                finally
                    close(io)
                    # @info "Closed:  $io" ##
                end
            end
        end
    catch e
        if typeof(e) <: InterruptException
            @warn Symbol(:server_, :listen) "Interrupted: listen($host,$port)"
        else
            rethrow(e)
        end
    finally
        close(tcpserver)
    end

    return
end

# HTTP.jl - Servers.jl
function _handle_connection(f::Function, c::Connection;
                           reuse_limit::Int=nolimit,
                           readtimeout::Int=0, kw...)

    wait_for_timeout = Ref{Bool}(true)
    if readtimeout > 0
        @async while wait_for_timeout[]
            @show inactiveseconds(c)
            if inactiveseconds(c) > readtimeout
                @warn Symbol(:server_, :_handle_connection) "Timeout: $c"
                writeheaders(c.io, Response(408, ["Connection" => "close"]))
                close(c)
                break
            end
            sleep(8 + rand() * 4)
        end
    end

    try
        count = 0
        while isopen(c)
            io = Transaction(c)
            final = count == reuse_limit
            _handle_transaction(f, io; final_transaction=final,
                                      kw...)
            if count == reuse_limit
                @info "close " count reuse_limit
                close(c)
            end
            count += 1
        end
    finally
        wait_for_timeout[] = false
    end
    return
end

# HTTP.jl - Servers.jl
function _handle_transaction(f::Function, t::Transaction;
                            final_transaction::Bool=false,
                            verbose::Bool=false, kw...)
    request = HTTP.Request()
    http = Streams.Stream(request, t)

    try
        startread(http)
    catch e
        # @show typeof(e)
        # @show fieldnames(e)
        if e isa EOFError && isempty(request.method)
            return
# FIXME https://github.com/JuliaWeb/HTTP.jl/pull/178#pullrequestreview-92547066
#        elseif !isopen(http)
#            @warn "Connection closed"
#            return
        elseif e isa HTTP.ParseError
            # @error Symbol(:server_, :_handle_transaction) e ## bukdu WebSocketError
            status = e.code == :HEADER_SIZE_EXCEEDS_LIMIT  ? 413 : 400
            write(t, Response(status, body = string(e.code)))
            close(t)
            return
        else
            rethrow(e)
        end
    end

    if verbose
        @info http.message
    end

    response = request.response
    #response.status = 200
    if final_transaction || hasheader(request, "Connection", "close")
        setheader(response, "Connection" => "close")
    end

    @async try
        _handle_stream(f, http)
    catch e
        if isioerror(e)
            @warn Symbol(:server_, :_handle_transaction) e
            #   ArgumentError("stream is closed or unusable")
            #   write: broken pipe (EPIPE)
        else
            # @error Symbol(:server_, :_handle_transaction) e stacktrace(catch_backtrace()) ## bukdu WebSocketError
        end
        close(t)
    end
    return
end

# module Bukdu

max_stars_repo_name: wookay/TestBukdu

max_issues_repo_name: wookay/TestBukdu

max_forks_repo_name: wookay/TestBukdu


text:
#
# Object headers
#

const HM_NIL = 0x00
const HM_DATASPACE = 0x01
const HM_LINK_INFO = 0x02
const HM_DATATYPE = 0x03
const HM_FILL_VALUE_OLD = 0x04
const HM_FILL_VALUE = 0x05
const HM_LINK_MESSAGE = 0x06
const HM_EXTERNAL_FILE_LIST = 0x07
const HM_DATA_LAYOUT = 0x08
const HM_BOGUS = 0x09
const HM_GROUP_INFO = 0x0a
const HM_FILTER_PIPELINE = 0x0b
const HM_ATTRIBUTE = 0x0c
const HM_OBJECT_COMMENT = 0x0d
const HM_SHARED_MESSAGE_TABLE = 0x0f
const HM_OBJECT_HEADER_CONTINUATION = 0x10
const HM_SYMBOL_TABLE = 0x11
const HM_MODIFICATION_TIME = 0x12
const HM_BTREE_K_VALUES = 0x13
const HM_DRIVER_INFO = 0x14
const HM_ATTRIBUTE_INFO = 0x15
const HM_REFERENCE_COUNT = 0x16

MESSAGE_TYPES = Dict(
    0x00 => "HM_NIL",
    0x01 => "HM_DATASPACE",
    0x02 => "HM_LINK_INFO",
    0x03 => "HM_DATATYPE",
    0x04 => "HM_FILL_VALUE_OLD",
    0x05 => "HM_FILL_VALUE",
    0x06 => "HM_LINK_MESSAGE",
    0x07 => "HM_EXTERNAL_FILE_LIST",
    0x08 => "HM_DATA_LAYOUT",
    0x09 => "HM_BOGUS",
    0x0a => "HM_GROUP_INFO",
    0x0b => "HM_FILTER_PIPELINE",
    0x0c => "HM_ATTRIBUTE",
    0x0d => "HM_OBJECT_COMMENT",
    0x0f => "HM_SHARED_MESSAGE_TABLE",
    0x10 => "HM_OBJECT_HEADER_CONTINUATION",
    0x11 => "HM_SYMBOL_TABLE",
    0x12 => "HM_MODIFICATION_TIME",
    0x13 => "HM_BTREE_K_VALUES",
    0x14 => "HM_DRIVER_INFO",
    0x15 => "HM_ATTRIBUTE_INFO",
    0x16 => "HM_REFERENCE_COUNT",
    )

const OH_ATTRIBUTE_CREATION_ORDER_TRACKED = 2^2
const OH_ATTRIBUTE_CREATION_ORDER_INDEXED = 2^3
const OH_ATTRIBUTE_PHASE_CHANGE_VALUES_STORED = 2^4
const OH_TIMES_STORED = 2^5

const OBJECT_HEADER_CONTINUATION_SIGNATURE = htol(0x4b48434f) # "OCHK"

struct ObjectStart
    signature::UInt32
    version::UInt8
    flags::UInt8
end
ObjectStart(flags::UInt8) = ObjectStart(OBJECT_HEADER_SIGNATURE, 2, flags)
define_packed(ObjectStart)

# Reads the start of an object including the signature, version, flags,
# and (payload) size. Returns the size.
# function read_obj_start(io::IO)
#     os = jlread(io, ObjectStart)
#     os.signature == OBJECT_HEADER_SIGNATURE || throw(InvalidDataException())
#     os.version == 2 || throw(UnsupportedVersionException())

#     if (os.flags & OH_TIMES_STORED) != 0
#         # Skip access, modification, change and birth times
#         skip(io, 128)
#     end
#     if (os.flags & OH_ATTRIBUTE_PHASE_CHANGE_VALUES_STORED) != 0
#         # Skip maximum # of attributes fields
#         skip(io, 32)
#     end

#     read_size(io, os.flags)
# end

function read_obj_start(io::IO)
    curpos = position(io)
    os = jlread(io, ObjectStart)
    if os.version == 2 && os.signature == OBJECT_HEADER_SIGNATURE
        if (os.flags & OH_TIMES_STORED) != 0
            # Skip access, modification, change and birth times
            skip(io, 128)
        end
        if (os.flags & OH_ATTRIBUTE_PHASE_CHANGE_VALUES_STORED) != 0
            # Skip maximum # of attributes fields
            skip(io, 32)
        end

        return read_size(io, os.flags), 2
    else
        seek(io, curpos)
        version = jlread(io, UInt8)
        version == 1 || throw(error("This should not have happened"))

        jlread(io, UInt8)
        num_messages = jlread(io, UInt16)
        #println("Reading $num_messages Messages")
        obj_ref_count = jlread(io, UInt32)
        #println("obj_ref_count = $obj_ref_count")
        obj_header_size = jlread(io, UInt32)
        #println("obj_header_size = $obj_header_size")
        #@info "If this errors, may need to skip to next 8bit alignment" maxlog = 1
        return obj_header_size, 1
    end
end

struct HeaderMessage
    msg_type::UInt8
    size::UInt16
    flags::UInt8
end
define_packed(HeaderMessage)

function isgroup(f::JLDFile, roffset::RelOffset)
    io = f.io
    seek(io, fileoffset(f, roffset))

    sz, version = read_obj_start(io)
    pmax::Int64 = position(io) + sz
    if version == 2
        while position(io) <= pmax-4
            msg = jlread(io, HeaderMessage)
            endpos = position(io) + msg.size
            if msg.msg_type == HM_LINK_INFO || msg.msg_type == HM_GROUP_INFO || msg.msg_type == HM_LINK_MESSAGE
                return true
            elseif msg.msg_type == HM_DATASPACE || msg.msg_type == HM_DATATYPE || msg.msg_type == HM_FILL_VALUE || msg.msg_type == HM_DATA_LAYOUT
                return false
            end
            seek(io, endpos)
        end
    elseif version == 1
        continuation_offset = -1
        continuation_length = 0
        chunk_end = pmax
        curpos = position(io)
        seek(io, curpos + 8 - mod1(curpos, 8))
        while true
            if continuation_offset != -1
                seek(io, continuation_offset)
                chunk_end = continuation_offset + continuation_length
                continuation_offset = -1
            end
            while position(io) <= chunk_end
                curpos = position(io)
                msg_type = jlread(io, UInt16)
                msg_size = jlread(io, UInt16)
                msg_flags = jlread(io, UInt8)
                jlread(io, UInt8); jlread(io, UInt16)
                msg = (;msg_type, size=msg_size, flags=msg_flags)
                endpos = position(io) + msg.size
                endpos = endpos + 8 - mod1(endpos, 8)
                if msg.msg_type == HM_LINK_INFO || msg.msg_type == HM_GROUP_INFO || msg.msg_type == HM_LINK_MESSAGE
                    return true
                elseif msg.msg_type == HM_DATASPACE || msg.msg_type == HM_DATATYPE || msg.msg_type == HM_FILL_VALUE || msg.msg_type == HM_DATA_LAYOUT
                    return false
                elseif msg_type == HM_OBJECT_HEADER_CONTINUATION
                    continuation_offset = chunk_start_offset = fileoffset(f, jlread(io, RelOffset))
                    continuation_length = jlread(io, Length)
                end
                seek(io, endpos)
            end
            continuation_offset == -1 && break
        end
    end
    return false
end

function print_header_messages(f::JLDFile, roffset::RelOffset)
    OffsetType = uintofsize(f.superblock.size_of_offsets)
    io = f.io
    chunk_start_offset::Int64 = fileoffset(f, roffset)
    seek(io, chunk_start_offset)

    # Test for V1 Obj header

    header_version = jlread(io, UInt8)
    if header_version == 1
        println("Object Header Message Version 1")
        cio = io
        skip(cio, 1)
        num_messages = jlread(cio, UInt16)
        obj_ref_count = jlread(cio, UInt32)
        obj_header_size = jlread(cio, UInt32)
        println("""
        Reading $num_messages Messages
            obj_ref_count = $obj_ref_count
            obj_header_size = $obj_header_size""")

        chunk_end = position(cio) + obj_header_size
        # Skip to nearest 8byte aligned position
        skip(cio, 4)
        #skip_to_aligned!(cio)
    elseif header_version == 2
        println("Object Header Message Version 2")

        seek(io, chunk_start_offset)
        cio = begin_checksum_read(io)
        sz = read_obj_start(cio)
        chunk_end = position(cio) + sz
    else
        throw(UnsupportedVersionException("Unknown object header version $header_version"))
    end
    # Messages
    continuation_message_goes_here::Int64 = -1
    links = OrderedDict{String,RelOffset}()

    continuation_offset::Int64 = -1
    continuation_length::Length = 0
    next_link_offset::Int64 = -1
    link_phase_change_max_compact::Int64 = -1
    link_phase_change_min_dense::Int64 = -1
    est_num_entries::Int64 = 4
    est_link_name_len::Int64 = 8
    chunk_end::Int64
    attrs = EMPTY_READ_ATTRIBUTES
    while true
        if continuation_offset != -1
            seek(io, continuation_offset)
            chunk_end = continuation_offset + continuation_length
            continuation_offset = -1
            if header_version == 2
                chunk_end = chunk_end #- 4 # leave space for checksum
                cio = begin_checksum_read(io)
                jlread(cio, UInt32) == OBJECT_HEADER_CONTINUATION_SIGNATURE || throw(InvalidDataException())
            end
        end

        while (curpos = position(cio)) <= chunk_end-4
            if header_version == 1
                println("now at position $(position(cio))")
                #skip(cio, 2)
                #println("skipped two bytes for no reason")
                msg_type = jlread(cio, UInt16)
                msg_size = jlread(cio, UInt16)
                flags = jlread(cio, UInt8)
                skip(cio, 3)
                endpos = curpos + 8 + msg_size
                endpos = endpos + 8 - mod1(endpos-chunk_start_offset, 8)
                msg = (; msg_type, size=msg_size, flags)
            else # version == 2
                msg = jlread(cio, HeaderMessage)
                endpos = curpos + jlsizeof(HeaderMessage) + msg.size
            end
            println("""
            Header Message:  $(MESSAGE_TYPES[msg.msg_type]) ($(msg.msg_type)))
                size: $(msg.size)
                flags: $(msg.flags)""")
            if msg.msg_type == HM_NIL
                if continuation_message_goes_here == -1 &&
                    chunk_end - curpos == CONTINUATION_MSG_SIZE
                    continuation_message_goes_here = curpos
                elseif endpos + CONTINUATION_MSG_SIZE == chunk_end
                    # This is the remaining space at the end of a chunk
                    # Use only if a message can potentially fit inside
                    # Single Character Name Link Message has 13 bytes payload
                    if msg.size >= 13
                        next_link_offset = curpos
                    end
                end
            else
                continuation_message_goes_here = -1
                if msg.msg_type == HM_LINK_INFO
                    fullmsg = read_link_info(cio, msg, OffsetType)
                    for (k,v) in pairs(fullmsg.fields)
                        println("    $k: $v")
                    end
                elseif msg.msg_type == HM_GROUP_INFO
                    if msg.size > 2
                        # Version Flag
                        jlread(io, UInt8) == 0 || throw(UnsupportedFeatureException())
                        flag = jlread(io, UInt8)
                        if flag%2 == 1 # first bit set
                            link_phase_change_max_compact = jlread(io, UInt16)
                            link_phase_change_min_dense = jlread(io, UInt16)
                            println("    link_phase_change_max_compact = $link_phase_change_max_compact")
                            println("    link_phase_change_min_dense = $link_phase_change_min_dense")
                        end
                        if (flag >> 1)%2 == 1 # second bit set
                            # Verify that non-default group size is given
                            est_num_entries = jlread(io, UInt16)
                            est_link_name_len = jlread(io, UInt16)
                            println("    est_num_entries = $est_num_entries")
                            println("    est_link_name_len = $est_link_name_len")
                        end
                    end
                elseif msg.msg_type == HM_LINK_MESSAGE
                    name, loffset = read_link(cio)
                    links[name] = loffset
                    println("   name = \"$name\"")
                    println("   offset = $(Int(loffset.offset))")
                elseif msg.msg_type == HM_OBJECT_HEADER_CONTINUATION
                    continuation_offset = chunk_start_offset = fileoffset(f, jlread(cio, RelOffset))
                    continuation_length = jlread(cio, Length)
                    # For correct behaviour, empty space can only be filled in the
                    # very last chunk. Forget about previously found empty space
                    next_link_offset = -1
                    println("""    offset = $(continuation_offset)\n    length = $(continuation_length)""")
                elseif msg.msg_type == HM_DATASPACE
                    dataspace = read_dataspace_message(cio)
                    println("    $dataspace")
                elseif msg.msg_type == HM_DATATYPE
                    datatype_class, datatype_offset = read_datatype_message(cio, f, (msg.flags & 2) == 2)
                    println("""    class: $class\n    offset: $datatype_offset""")
                elseif msg.msg_type == HM_FILL_VALUE_OLD
                    #(jlread(cio, UInt8) == 3 && jlread(cio, UInt8) == 0x09) || throw(UnsupportedFeatureException())
                elseif msg.msg_type == HM_FILL_VALUE
                    #(jlread(cio, UInt8) == 3 && jlread(cio, UInt8) == 0x09) || throw(UnsupportedFeatureException())
                elseif msg.msg_type == HM_DATA_LAYOUT
                    version = jlread(cio, UInt8)
                    println("""    version: $version""")
                    if version == 4 || version == 3
                        storage_type = jlread(cio, UInt8)
                        if storage_type == LC_COMPACT_STORAGE
                            data_length = jlread(cio, UInt16)
                            data_offset = position(cio)
                            println("""    type: compact storage\n    length: $length\n    offset: $(data_offset)""")
                        elseif storage_type == LC_CONTIGUOUS_STORAGE
                            data_offset = fileoffset(f, jlread(cio, RelOffset))
                            data_length = jlread(cio, Length)
                            println("""    type: contiguous storage\n    length: $length\n    offset: $(data_offset)""")

                        elseif storage_type == LC_CHUNKED_STORAGE
                            # TODO: validate this
                            flags = jlread(cio, UInt8)
                            dimensionality = jlread(cio, UInt8)
                            dimensionality_size = jlread(cio, UInt8)
                            skip(cio, Int(dimensionality)*Int(dimensionality_size))

                            chunk_indexing_type = jlread(cio, UInt8)
                            chunk_indexing_type == 1 || throw(UnsupportedFeatureException("Unknown chunk indexing type"))
                            data_length = jlread(cio, Length)
                            jlread(cio, UInt32)
                            data_offset = fileoffset(f, jlread(cio, RelOffset))
                            chunked_storage = true
                            println("""    type: chunked storage
                                    length: $length
                                    offset: $(data_offset)
                                    dimensionality: $dimensionality
                                    dimensionality_size: $dimensionality_size
                                    chunk indexing type: $chunk_indexing_type""")

                        else
                            throw(UnsupportedFeatureException("Unknown data layout"))
                        end
                    end
                elseif msg.msg_type == HM_FILTER_PIPELINE
                    version = jlread(cio, UInt8)
                    version == 2 || throw(UnsupportedVersionException("Filter Pipeline Message version $version is not implemented"))
                    nfilters = jlread(cio, UInt8)
                    nfilters == 1 || throw(UnsupportedFeatureException())
                    filter_id = jlread(cio, UInt16)
                    issupported_filter(filter_id) || throw(UnsupportedFeatureException("Unknown Compression Filter $filter_id"))
                elseif msg.msg_type == HM_SYMBOL_TABLE
                    v1_btree_address = jlread(cio, OffsetType)
                    local_heap_address = jlread(cio, OffsetType)
                    println("""    required for \"old style" groups\n    \
                            v1 B-Tree Adress: $(Int(v1_btree_address))\n    \
                            Local Heap Adress: $(Int(local_heap_address))""")
                elseif msg.msg_type == HM_ATTRIBUTE
                    if attrs === EMPTY_READ_ATTRIBUTES
                        attrs = ReadAttribute[read_attribute(cio, f)]
                    else
                        push!(attrs, read_attribute(cio, f))
                    end
                    attr = attrs[end]
                    println("""    name: "$(attr.name)"\n    \
                                datatype: $(DATATYPES[attr.datatype_class])""")
                    try
                        data = read_attr_data(f, attr)
                        println("""    data: "$data" """)
                    catch e
                        println("""    loading data failed""")
                    end
                elseif (msg.flags & 2^3) != 0
                    throw(UnsupportedFeatureException())
                end
            end
            seek(cio, endpos)
        end

        # Checksum
        seek(cio, chunk_end)
        if header_version == 2
            end_checksum(cio) == jlread(io, UInt32) || throw(InvalidDataException())
        end
        continuation_offset == -1 && break
    end
    nothing
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct PacketCaptureResult <: SwaggerModel
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    id::Any # spec type: Union{ Nothing, String } # spec name: id
    etag::Any # spec type: Union{ Nothing, String } # spec name: etag
    properties::Any # spec type: Union{ Nothing, PacketCaptureResultProperties } # spec name: properties

    function PacketCaptureResult(;name=nothing, id=nothing, etag=nothing, properties=nothing)
        o = new()
        validate_property(PacketCaptureResult, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(PacketCaptureResult, Symbol("id"), id)
        setfield!(o, Symbol("id"), id)
        validate_property(PacketCaptureResult, Symbol("etag"), etag)
        setfield!(o, Symbol("etag"), etag)
        validate_property(PacketCaptureResult, Symbol("properties"), properties)
        setfield!(o, Symbol("properties"), properties)
        o
    end
end # type PacketCaptureResult

const _property_map_PacketCaptureResult = Dict{Symbol,Symbol}(Symbol("name")=>Symbol("name"), Symbol("id")=>Symbol("id"), Symbol("etag")=>Symbol("etag"), Symbol("properties")=>Symbol("properties"))
const _property_types_PacketCaptureResult = Dict{Symbol,String}(Symbol("name")=>"String", Symbol("id")=>"String", Symbol("etag")=>"String", Symbol("properties")=>"PacketCaptureResultProperties")
Base.propertynames(::Type{ PacketCaptureResult }) = collect(keys(_property_map_PacketCaptureResult))
Swagger.property_type(::Type{ PacketCaptureResult }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_PacketCaptureResult[name]))}
Swagger.field_name(::Type{ PacketCaptureResult }, property_name::Symbol) =  _property_map_PacketCaptureResult[property_name]

function check_required(o::PacketCaptureResult)
    true
end

function validate_property(::Type{ PacketCaptureResult }, name::Symbol, val)
end


text:
tomlpath = joinpath(@__DIR__, "sbm_config.toml")

@testset "BMI" begin

    @testset "BMI functions" begin

        model = BMI.initialize(Wflow.Model, tomlpath)

        @testset "initialization and time functions" begin
            @test BMI.get_time_units(Wflow.Model) == "seconds since 1970-01-01T00:00:00"
            @test BMI.get_time_step(model) == 86400.0
            @test BMI.get_start_time(model) == 9.467712e8
            @test BMI.get_current_time(model) == 9.467712e8
            @test BMI.get_end_time(model) == 9.493632e8
        end

        @testset "model information functions" begin
            @test BMI.get_component_name(model) == "sbm"
            @test BMI.get_input_item_count(model) == 186
            @test BMI.get_output_item_count(model) == 186
            @test BMI.get_input_var_names(model)[[1, 5, 120, 186]] == [
                "vertical.Δt",
                "vertical.n_unsatlayers",
                "lateral.land.q_av",
                "lateral.river.reservoir.evaporation",
            ]
        end

        @testset "variable information functions" begin
            @test BMI.get_var_grid(model, "vertical.θₛ") == 4
            @test BMI.get_var_grid(model, "lateral.river.h") == 2
            @test BMI.get_var_grid(model, "lateral.river.reservoir.inflow") == 0
            @test_throws ErrorException BMI.get_var_grid(model, "lateral.river.lake.volume")
            # Vector{Float64} printing on Julia 1.6+
            @test BMI.get_var_type(model, "lateral.river.reservoir.inflow") in
                  ["Array{$Float,1}", "Vector{$Float}"]
            @test BMI.get_var_type(model, "vertical.n") == "Int64"
            @test BMI.get_var_units(model, "vertical.θₛ") == "-"
            @test BMI.get_var_itemsize(model, "lateral.subsurface.ssf") == sizeof(Float)
            @test BMI.get_var_nbytes(model, "vertical.n") == 8
            @test BMI.get_var_nbytes(model, "lateral.river.q") ==
                  length(model.lateral.river.q) * sizeof(Float)
            @test BMI.get_var_location(model, "lateral.river.q") == "node"
        end

        model = BMI.update(model)

        @testset "update and get and set functions" begin
            @test BMI.get_current_time(model) == 9.468576e8
            @test mean(BMI.get_value(model, "vertical.zi")) ≈ 276.3767651555451
            @test BMI.get_value_at_indices(model, "lateral.river.q", [1, 100, 5617]) ≈
                  [0.6211503865184697, 5.219305686635002, 0.026163746306482282]
            BMI.set_value(model, "vertical.zi", fill(300.0, length(model.vertical.zi)))
            @test mean(BMI.get_value(model, "vertical.zi")) == 300.0
            BMI.set_value_at_indices(model, "vertical.zi", [1], [250.0])
            @test BMI.get_value_at_indices(model, "vertical.zi", [1, 2]) == [250.0, 300.0]
        end

        @testset "model grid functions" begin
            @test BMI.get_grid_type(model, 0) == "unstructured"
            @test BMI.get_grid_rank(model, 0) == 2
            @test BMI.get_grid_size(model, 0) == 2
            @test BMI.get_grid_size(model, 2) == 5809
            @test BMI.get_grid_size(model, 4) == 50063
            @test BMI.get_grid_shape(model, 4) == 50063
            @test minimum(BMI.get_grid_x(model, 4)) ≈ 5.426666666666667f0
            @test maximum(BMI.get_grid_x(model, 4)) ≈ 7.843333333333344f0
            @test BMI.get_grid_x(model, 0) ≈ [5.760000000000002f0, 5.918333333333336f0]
            @test BMI.get_grid_y(model, 0) ≈ [48.92583333333333f0, 49.909166666666664f0]
            @test BMI.get_grid_node_count(model, 0) == 2
        end

        @testset "update until and finalize" begin
            time = BMI.get_current_time(model) + 2 * BMI.get_time_step(model)
            @test_logs (:info, "update model until 9.470304e8")
            model = BMI.update_until(model, time)
            BMI.finalize(model)
        end

    end

    @testset "BMI run SBM in parts" begin
        tomlpath = joinpath(@__DIR__, "sbm_gw.toml")
        model = BMI.initialize(Wflow.Model, tomlpath)

        # update the recharge part of the SBM model
        model = BMI.update(model, run = "sbm_until_recharge")

        @testset "recharge part of SBM" begin
            sbm = model.vertical
            @test sbm.interception[1] ≈ 0.6299999952316284f0
            @test sbm.ustorelayerdepth[1][1] ≈ 0.0f0
            @test sbm.snow[1] ≈ 3.1912317735997524f0
            @test sbm.recharge[5] ≈ -0.0727941579914808f0
            @test sbm.zi[5] ≈ 300.0f0
        end

        # set zi and exfiltwater from external source (e.g. a groundwater model)
        BMI.set_value(
            model,
            "lateral.subsurface.zi",
            fill(0.25, BMI.get_grid_size(model, 4)),
        )
        BMI.set_value(
            model,
            "lateral.subsurface.exfiltwater",
            fill(1.0e-5, BMI.get_grid_size(model, 4)),
        )
        # update SBM after subsurface flow
        model = BMI.update(model, run = "sbm_after_subsurfaceflow")

        @testset "SBM after subsurface flow" begin
            sbm = model.vertical
            sub = model.lateral.subsurface
            @test sbm.interception[1] ≈ 0.6299999952316284
            @test sbm.ustorelayerdepth[1][1] ≈ 0.0f0
            @test sbm.snow[1] ≈ 3.1912317735997524f0
            @test sbm.recharge[5] ≈ -0.0727941579914808f0
            @test sbm.zi[5] ≈ 250.0f0
            @test sub.zi[5] ≈ 0.25f0
            @test sub.exfiltwater[1] ≈ 1.0f-5
            @test sub.ssf[1] ≈ 0.0f0
        end

        BMI.finalize(model)
    end
end


text:
using Oceananigans.Operators: Δzᵃᵃᶜ

"""
Update the hydrostatic pressure perturbation pHY′. This is done by integrating
the `buoyancy_perturbation` downwards:

    `pHY′ = ∫ buoyancy_perturbation dz` from `z=0` down to `z=-Lz`
"""
@kernel function update_hydrostatic_pressure!(pHY′, grid, buoyancy, C)
    i, j = @index(Global, NTuple)

    @inbounds pHY′[i, j, grid.Nz] = - ℑzᵃᵃᶠ(i, j, grid.Nz+1, grid, z_dot_g_b, buoyancy, C) * Δzᵃᵃᶜ(i, j, grid.Nz+1, grid)

    @unroll for k in grid.Nz-1 : -1 : 1
        @inbounds pHY′[i, j, k] = pHY′[i, j, k+1] - ℑzᵃᵃᶠ(i, j, k+1, grid, z_dot_g_b, buoyancy, C) * Δzᵃᵃᶜ(i, j, k+1, grid)
    end
end


text:
using Optim
using HTTP
using GLM
using LinearAlgebra
using Random
using Statistics
using DataFrames
using Distributions
using CSV
using FreqTables
using CategoricalArrays
using Tables
using StatsModels
cd("C:/Users/jakek/Documents/econ_phd/structural_modeling_course/forked_class/ProblemSets/PS4-mixture/")

# 1

df = CSV.read("nlsw88t.csv")
X = [df.age df.white df.collgrad]
Z = hcat(df.elnwage1, df.elnwage2, df.elnwage3, df.elnwage4,
df.elnwage5, df.elnwage6, df.elnwage7, df.elnwage8)
y = df.occ_code

function mlogit(alpha, X, d,Z)
	T = promote_type(eltype(X),eltype(alpha), eltype(Z), eltype(d))
	Zbar = zeros(T, size(Z,1),1)
	u = zeros(T, size(Z,1),7)
	help = zeros(T, size(Z,1),7)
	loglike = zeros(T, 1)
	Zbar = Z .- Z[:,8]
	u = reduce(hcat,[X*alpha[(3*i-2):(3*i)] .+ Zbar[:,i] .*alpha[22] for i in 1:7])
	help = reduce(hcat,[d.==i for i in 1:7])
	loglike = sum(help .*u) - sum(log.(1 .+sum(exp.(u), dims=2)))
	return -loglike
end

# starting values - pull from last homework
start_ps3 = CSV.read("../PS3-gev/q1.csv")
starting = vec(reduce(vcat,[start_ps3[:,i] for i in 1:8])[Not([23 24]),:])

td = TwiceDifferentiable(b -> mlogit(b, X, y,Z), vec(starting); autodiff = :forward)
beta_hat_mlogit = optimize(td,starting, Optim.Options(g_tol=1e-5, iterations=100_000, show_trace=true))
H  = Optim.hessian!(td, beta_hat_mlogit.minimizer)
beta_hat_mlogit_se = sqrt.(diag(inv(H)))

res = vcat(beta_hat_mlogit.minimizer', beta_hat_mlogit_se' )
CSV.write("q1.csv", convert(DataFrame,res))

#2.  now the coefficient on wage is positive, reflecting that increased wages results in higher utility. This makes more sense than before.

#########################################
### 3.
#########################################
### a. Qudrature learning
include("lgwt.jl")
d = Normal(0,1)
nodes, weights = lgwt(7,-4,4)
sum(weights.*pdf.(d,nodes)) # pdf integration should be near 1.
sum(weights.*nodes.*pdf.(d,nodes)) # expectation shoudl bear near 0
### b. quadrature practice
d2 = Normal(0,2)
nodes, weights = lgwt(7,-10,10)
sum(weights.*nodes.^2 .*pdf.(d2,nodes))
nodes, weights = lgwt(10,-10,10)
sum(weights.*nodes.^2 .*pdf.(d2,nodes))
# the second option is much closer to the true variance, 2.

### c.
a= -10
b =10
d=10000000
nds = rand(d,1) .*(b-a) .+a
print((b-a)/d * sum(pdf.(d2,nds))) # should be close to 1.
print(sum(nds .^2 .* ((b-a)/d .* pdf.(d2,nds))))
print(sum(nds .* ((b-a)/d .* pdf.(d2,nds))))

#### 4. Quadrature.
function mixlogit(alpha, X, d,Z)
	include("lgwt.jl")
	T = promote_type(eltype(X),eltype(alpha), eltype(Z), eltype(d))
	Zbar = zeros(T, size(Z,1),1)
	u = repeat(zeros(T, size(Z,1),7),1,1,7)
	help = zeros(T, size(Z,1),7)
	loglike = zeros(T, 1)
	Zbar = Z .- Z[:,8]
	sig = exp(alpha[23])
	dist = Normal(alpha[22],sig)
	nodes, weights = lgwt(7,-5*sig+alpha[22],5*sig+alpha[22])
	# make 3d array, with quad points.
	for j in 1:7
		u[:,:,j] = exp.(reduce(hcat,[X*alpha[(3*i-2):(3*i)] .+ Zbar[:,i] .*nodes[j] for i in 1:7]))
	end
	usum =  sum(u,dims=2)
	help = repeat(reduce(hcat,[d.==i for i in 1:7]),1,1,7)
	premix =  prod(u .^help, dims=2) ./usum
	mix = reduce(hcat,[weights[i]*pdf(dist,nodes[i]) .* premix[:,:,i] for i in 1:7])
	loglike = sum(log.(sum(mix,dims=2)))
	return -loglike
end

starting = append!(vec(starting), 0)
td = TwiceDifferentiable(b -> mixlogit(b, X, y,Z), starting; autodiff = :forward)
beta_mixlogit = optimize(td,starting, Optim.Options(g_tol=1e-5, iterations=100_000, show_trace=true))
H  = Optim.hessian!(td, beta_mixlogit.minimizer)
#beta_hat_mlogit_se = sqrt.(diag(inv(H)))
# res = vcat(beta_mixlogit.minimizer', beta_hat_mlogit_se' )
CSV.write("q4.csv", convert(DataFrame,beta_mixlogit.minimizer'))

#### 5. simulation
function mixlogit_mc(alpha, X, d,Z)
	include("lgwt.jl")
	T = promote_type(eltype(X),eltype(alpha), eltype(Z), eltype(d))
	Zbar = zeros(T, size(Z,1),1)
	u = repeat(zeros(T, size(Z,1),7),1,1,1000000)
	help = zeros(T, size(Z,1),7)
	loglike = zeros(T, 1)
	Zbar = Z .- Z[:,8]
	dist = Normal(alpha[22],sig)
	sig = exp(alpha[23])
	nodes = rand(1000000,1) .*(2*sig*5) .+ alpha[22]
	weights = (2*sig*5)/1000000
	# make 3d array, with quad points.
	for j in 1:1000000
		u[:,:,j] = exp.(reduce(hcat,[X*alpha[(3*i-2):(3*i)] .+ Zbar[:,i] .*nodes[j] for i in 1:7]))
	end
	usum =  sum(u,dims=2)
	help = repeat(reduce(hcat,[d.==i for i in 1:7]),1,1,10)
	premix =  prod(u .^help, dims=2) ./usum
	mix = reduce(hcat,[weights[i]*pdf(dist,nodes[i]) .* premix[:,:,i] for i in 1:1000000])
	loglike = sum(mix)
	return -loglike
end


text:
function canonicalize_and_sign(env::AWSEnv, service, method, params)
    if env.sig_ver == 2
      signature_version_2(env, service, method, copy(params))
    elseif env.sig_ver == 4
      signature_version_4(env, service, method, copy(params))
    else
      error("invalid signature version $(env.sig_ver)")
    end
end
export canonicalize_and_sign

#=
import URIParser
function urlencode_params(params::Vector{Tuple})
    buf = IOBuffer()
    for param in params
        print(buf, URIParser.escape(param[1]))
        print(buf, '=')
        print(buf, URIParser.escape(param[2]))
        print(buf, '&')
    end
    if buf.size > 0
        buf.size -= 1  # trim off last &
    end
    return takebuf_string(buf)
end

function urlencode_form(params::Vector{Tuple})
    buf = IOBuffer()
    for param in params
        print(buf, URIParser.escape_form(param[1]))
        print(buf, '=')
        print(buf, URIParser.escape_form(param[2]))
        print(buf, '&')
    end
    if buf.size > 0
        buf.size -= 1  # trim off last &
    end
    return takebuf_string(buf)
end
=#

function signature_version_2(env::AWSEnv, service, use_post, request_parameters)
    if env.aws_token != ""
        push!(request_parameters, ("SecurityToken", env.aws_token))
    end
    push!(request_parameters, ("SignatureVersion", "2"))
    push!(request_parameters, ("SignatureMethod", "HmacSHA256"))
    sort!(request_parameters)

    canonical_querystring = Requests.format_query_str(request_parameters)

	method = use_post ? "POST" : "GET"

    string_to_sign = method * "\n" * ep_host(env, service) * "\n" * env.ep_path * "\n" * canonical_querystring

    signature = String(copy(base64encode(sign(env.aws_seckey, string_to_sign))))

    push!(request_parameters, ("Signature", signature))

    if (env.dbg) || (env.dry_run)
        info("Parameters:")
        for (k, v) in request_parameters
            info("  $k => $v")
        end
        info("--------\nString to sign:\n" * string_to_sign * "\n--------")
    end

    return (Tuple[], Requests.format_query_str(request_parameters))
end

function signature_version_4(env::AWSEnv, service, use_post, request_parameters)
    # inputs to request
    sort!(request_parameters)
    amzdate = get_utc_timestamp(; basic=true)
    datestamp = amzdate[1:searchindex(amzdate, "T")-1]
    payload = "" ##### Assume empty payload for the moment.

	if use_post
        ## payload = urlencode_form(request_parameters)
        payload = Requests.format_query_str(request_parameters)
        request_parameters = Array{Tuple}(0)
    end

	method = use_post ? "POST" : "GET"

    # Task 1: canonical request
    canonical_uri = env.ep_path
    canonical_querystring = Requests.format_query_str(request_parameters)
    ## canonical_querystring = urlencode_params(request_parameters)
    canonical_headers = "host:" * ep_host(env, service) * "\n" * "x-amz-date:" * amzdate * "\n"
    signed_headers = "host;x-amz-date"
    payload_hash = bytes2hex(Crypto.sha256(payload))
    canonical_request = method * "\n" * canonical_uri * "\n" * canonical_querystring * "\n" *
        canonical_headers * "\n" * signed_headers * "\n" * payload_hash

    # Task 2: string to sign
    algorithm = "AWS4-HMAC-SHA256"
	credential_scope = datestamp * "/" * env.region * "/" * service * "/" * "aws4_request"

    string_to_sign = algorithm * "\n" * amzdate * "\n" * credential_scope * "\n" * bytes2hex(Crypto.sha256(canonical_request))

    # Task 3: calculate the signature
    signing_key = get_signature_key(env.aws_seckey, datestamp, env.region, service, method)
    signature = bytes2hex(sign(signing_key, string_to_sign))

    # Task 4: add signing information to request ##### how to return header to caller?
    authorization_header = algorithm * " " * "Credential=" * env.aws_id * "/" * credential_scope * ", " *
        "SignedHeaders=" * signed_headers * ", " * "Signature=" * signature
    headers = [("Authorization", authorization_header),
              ("X-Amz-Date", amzdate)]
    if env.aws_token != ""
        push!(headers, ("X-Amz-Security-Token", env.aws_token))
    end

    if use_post
        push!(headers, ("Content-Type", "application/x-www-form-urlencoded"))
    end

    if (env.dbg) || (env.dry_run)
        info("Parameters:")
        for (k, v) in request_parameters
            info("  $k => $v")
        end
        info("--------\nCanonical request:\n" * canonical_request)
        info("--------\nString to sign:\n" * string_to_sign * "\n--------")
        info("Headers:")
        info("  Host: " * ep_host(env, service))
        for (k, v) in headers
            info("  $k: $v")
        end
    end

    return (headers, canonical_querystring, Requests.format_query_str(request_parameters))
    ## return (headers, canonical_querystring, urlencode_params(request_parameters))
end

function get_signature_key(key, datestamp, region, service, method)
    kDate = sign("AWS4" * key, datestamp)
    kRegion = sign(kDate, region)
    kService = sign(kRegion, service)
    kSigning = sign(kService, "aws4_request")
    return kSigning
end

sign(key, msg) = Crypto.hmacsha256_digest(msg, key)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

struct LogAnalyticsApi <: SwaggerApi
    client::Swagger.Client
end

"""

Export logs that show Api requests made by this subscription in the given time window to show throttling activities.
Param: parameters::RequestRateByIntervalInput (required)
Param: location::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: LogAnalyticsOperationResult
"""
function _swaggerinternal_logAnalyticsExportRequestRateByInterval(_api::LogAnalyticsApi, parameters, location::String, api_version::String, subscriptionId::String; _mediaType=nothing)

    _ctx = Swagger.Ctx(_api.client, "POST", LogAnalyticsOperationResult, "/subscriptions/{subscriptionId}/providers/Microsoft.Compute/locations/{location}/logAnalytics/apiAccess/getRequestRateByInterval", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "location", location)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function logAnalyticsExportRequestRateByInterval(_api::LogAnalyticsApi, parameters, location::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_logAnalyticsExportRequestRateByInterval(_api, parameters, location, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function logAnalyticsExportRequestRateByInterval(_api::LogAnalyticsApi, response_stream::Channel, parameters, location::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_logAnalyticsExportRequestRateByInterval(_api, parameters, location, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Export logs that show total throttled Api requests for this subscription in the given time window.
Param: parameters::ThrottledRequestsInput (required)
Param: location::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: LogAnalyticsOperationResult
"""
function _swaggerinternal_logAnalyticsExportThrottledRequests(_api::LogAnalyticsApi, parameters, location::String, api_version::String, subscriptionId::String; _mediaType=nothing)

    _ctx = Swagger.Ctx(_api.client, "POST", LogAnalyticsOperationResult, "/subscriptions/{subscriptionId}/providers/Microsoft.Compute/locations/{location}/logAnalytics/apiAccess/getThrottledRequests", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "location", location)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function logAnalyticsExportThrottledRequests(_api::LogAnalyticsApi, parameters, location::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_logAnalyticsExportThrottledRequests(_api, parameters, location, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function logAnalyticsExportThrottledRequests(_api::LogAnalyticsApi, response_stream::Channel, parameters, location::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_logAnalyticsExportThrottledRequests(_api, parameters, location, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export logAnalyticsExportRequestRateByInterval, logAnalyticsExportThrottledRequests


text:
using Test
using LinearAlgebra
using jInv.Mesh;
using FactoredEikonalFastMarching;
using Statistics
include("sensitivityTest.jl");

@testset "Sensitivity test" begin
	runSensitivityTest2D();
	runSensitivityTest3D();
	include("../examples/runExperiments.jl")
end


text:
# data/z-test.jl

using MIRT: ir_load_brainweb_t1_256
using Test: @test, @testset

@testset "brainweb" begin
    x = ir_load_brainweb_t1_256()
    @test extrema(x) == (0.0f0, 241.0f0)
end


text:
function read_arrays_jl_basic(filename)
  try
    keys = Int[]
    indexes = Int[]
    elements = Int[]
    for line in readlines(filename)
      s = split(line, "#")
      k, v = parse(Int, s[1]), parse.(Int, split(s[2], ","))
      m = length(v)
      keys = [keys; fill(k, m)]
      indexes = [indexes; 0:m-1]
      elements = [elements; v]
    end

    return keys, indexes, elements
  catch
    return Int[], Int[], Int[]
  end
end


text:
module Languages

using Unicode
using Niancat.Gameface

struct Word
    w::String
end

isinalphabet(c::Char) :: Bool = occursin(c, "abcdefghijklmnopqrstuvwxyzåäö")
ismeaningfuldiacritic(c::Char) :: Bool = occursin(c, "åäö")

function normalizechar(c::AbstractChar) :: String
    stripmark = !ismeaningfuldiacritic(c)
    Unicode.normalize(string(c), stripmark=stripmark)
end

function normalize(s::String) :: String
    lower = Unicode.normalize(s, casefold=true)
    normalizedchars = [normalizechar(c) for c in lower]

    normalizedword = join(normalizedchars)
    filter(isinalphabet, normalizedword)
end

function sortword(s::String) :: String
    ns = normalize(s)
    String(sort([c for c in ns]))
end

isanagram(s1::String, s2::String) = sortword(s1) == sortword(s2)

struct SwedishDictionary <: LanguageDictionary
    words::Set{Word}
    anagrams::Dict{String, Vector{String}}

    function SwedishDictionary(words::AbstractVector{<:String})
        normalizedwords = Word[]
        anagrams = Dict{String, Vector{String}}()
        for w in words
            nw = normalize(w)
            push!(normalizedwords, Word(nw))

            push!(get!(anagrams, sortword(nw), String[]), nw)
        end
        new(Set(normalizedwords), anagrams)
    end
end

Base.in(word::String, s::SwedishDictionary) = in(Word(normalize(word)), s.words)

findanagrams(dictionary::SwedishDictionary, s::String) = get(dictionary.anagrams, sortword(s), [])

normalizedword(::SwedishDictionary, word::String) :: String = normalize(word)

export SwedishDictionary, Word, isanagram, findanagrams, normalizedword

end


text:
using SharedArrays
using JLD
include("../examples/lorenz63.jl")
function ergodic_average(rho)
	s = [10., rho, 8/3]
	n = 100000000000
	n_runup = 5000
	z = 0.
	u = rand(3)
	u .= lorenz63(reshape(u,3,1),
			s, n_runup)[end,:,1]
	for i = 1:n
		u .= next(u, s)
		z += u[3]/n
		if i==n
		  @time u .= next(u,s)
		end
	end
	return z-rho
end
function compute_erg_avg()
	n_p = 250
	rho = LinRange(28.0, 28.25, n_p)
	z_rho = SharedArray{Float64, 1}(n_p)
	t = @distributed for i = 1:n_p
		@show i
		z_rho[i] = ergodic_average(rho[i])
	end
	wait(t)
	@show z_rho
	save("../data/erg_avg/Jar1_res.jld", "z_rho",
	     z_rho, "rho", rho)

end

max_stars_repo_name: nishaChandramoorthy/Lipschitz

max_issues_repo_name: nishaChandramoorthy/Lipschitz

max_forks_repo_name: nishaChandramoorthy/Lipschitz


text:
include("data/parser.jl")

import AbstractAlgebra
using BenchmarkTools
using Logging

global_logger(ConsoleLogger(stderr, Logging.Error))

BenchmarkTools.DEFAULT_PARAMETERS.samples = 2

function benchmark_system(system)
    @btime Groebner.groebner($s, ordering=:degrevlex)
end

function run_for_gleb()
    fns = [
        "gbSIWR",
        "gbSEAIJRC"
    ]

    for fn in fns
        system = load_system(fn)
        benchmark_system(system)
    end
end


text:
using DifferentialEquations
using Plots
n_iter=50
barrier=4
HH=range(.125001,stop=.499,length=n_iter)
exit_time=zeros(length(HH))
dE=zeros(length(HH))
max_time=exp(8)
tspan = (0.0,max_time)
pert=1e-6

ODE1(z,w)=conj(  im * w.*( 1 ./(z.^2-w.^2)+1 ./(1+w.^2) ))
ODE2(z,w)=conj(  im * z.*( 1 ./(w.^2-z.^2)+1 ./(1+z.^2) ))
condition(u,t,integrator)=abs(u[1])>barrier || abs(u[2])>barrier || abs(u[3])>barrier||abs(u[4])>barrier
affect!(integrator) = terminate!(integrator)
cb = DiscreteCallback(condition,affect!)
function Eq_of_M(du,u,p,t)
    du[1]=real(ODE1(u[1]+im*u[2],u[3]+im*u[4]))
    du[2]=imag(ODE1(u[1]+im*u[2],u[3]+im*u[4]))
    du[3]=real(ODE2(u[1]+im*u[2],u[3]+im*u[4]))
    du[4]=imag(ODE2(u[1]+im*u[2],u[3]+im*u[4]))
    return
end
# refernce
# z=X+iP W=Q+iY . Plane Q=0, P=0
# z=u[1]+i u[2] ; w=u[3]+i u[4]

for j=1:length(HH)
    H=HH[j]
    h=1/(2*H)
    Y_0=1/sqrt(h+1)
    u0= [0;0;pert;Y_0;0;0]
    prob = ODEProblem(Eq_of_M,u0,tspan)
    t,A=solve(prob, RK4(),reltol=1e-8,abstol=1e-10,maxiters=1e15,callback=cb)
    exit_time[j]=t.t[end]
    Z=A[1,:]+im *A[2,:];
    W=A[3,:]+im *A[4,:];
    N=length(t.t)
    Energy=zeros(N)
    for k=1:N
        Energy[k]=abs(1/(1+Z[k]^2)-1/(1+W[k]^2))
    end
    dE[j]=Energy[end]-Energy[1]

end

p=plot(HH,log.(exit_time),seriestype=:scatter,marker=:circle,msw = 0,ms=2,seriescolor = :black,legend=false)
savefig(p,"ExitLogLog1000b.pdf")
 plot(-HH,-log.(abs.(dE)))


text:
using Test, ApproxFun, SingularIntegralEquations

@testset "Nonlocal operators" begin
    α = 2.5
    δ = 0.1
    for sp in (Fourier(), SinSpace(), CosSpace(), Laurent(), Hardy{true}(PeriodicSegment()), Hardy{false}(PeriodicSegment()))
        L = NonlocalLaplacian(sp, α, δ)
        @test domain(L) == domain(sp)
        @test domainspace(L) == rangespace(L) == sp
        @test bandwidths(L) == (0,0)
    end

    NonlocalLaplacian(PeriodicSegment(), α, δ) === NonlocalLaplacian(Space(PeriodicSegment()), α, δ)

    L = NonlocalLaplacian(SinSpace(), α, δ)
    @test norm(L[10,10] + 9.836381919022900199967236210084047918127677209354580364785826568010653044652657e+01) < 100*norm(L[10,10])*eps()
    @test norm(L[1000,1000] + 1.670424124232011223460118736778899327539655416751949620343795973043019808627733e+05) < 100*norm(L[1000,1000])*eps()
end