text:
using Gtk
using Cairo
using BoardGames

mutable struct GuiData{G <: Game, M, B, STuple}
    movelist::Vector{M}
    boardlist::Vector{B}
    playable::Bool
    curr_visible_player::Int
    curr_visible_strategy::Vector{Int}
    playerboxes::Vector{GtkBox}
    strategyboxes::Vector{Vector{GtkBox}}
    elemententries::Vector{Vector{Vector{GtkEntry}}}
    playerstrategies::Vector{Int}
    allstrategies::Vector{STuple}
end

GuiData(::Type{G}, playerboxes, strategyboxes, elemententries, strategies) where G<:Game =
    GuiData{G, movetype(G), boardtype(G), typeof(strategies)}(
        movetype(G)[],
        boardtype(G)[initialboard(G())],
        false,
        0,
        ones(Int, nplayers(G)),
        playerboxes,
        strategyboxes,
        elemententries,
        ones(Int, nplayers(G)),
        [copy.(strategies) for i in 1:nplayers(G)]
    )

guidatagame(::GuiData{G,M,B}) where {G,M,B} = G

@guarded function mouseclick_cb(widgetptr::Ptr, eventptr::Ptr, user_data)
    widget = convert(GtkCanvas, widgetptr)
    event = unsafe_load(eventptr)

    (data, strategies) = user_data

    if data.playable
        ctx = Gtk.getgc(widget)
        w = Gtk.width(widget)
        h = Gtk.height(widget)

        x = event.x / w
        y = event.y / h
        move = correspondingmove(data.boardlist[end], x,y)

        if move in getmoves(data.boardlist[end])
            push!(data.movelist, move)
            push!(data.boardlist, play(data.boardlist[end], move))
        end

        println("x: $x, y: $y")

        playgame(data)

        reveal(widget)
    end
    nothing
end

@guarded function draw_cb(widgetptr::Ptr, cairoptr::Ptr, user_data)
    widget = convert(GtkCanvas, widgetptr)

    w, h = width(widget), height(widget)
    ctx = CairoContext(cairoptr)

    draw(ctx, user_data.boardlist[end], w, h)

    nothing
end

@guarded function play_cb(widgetptr::Ptr, user_data)
    widget = convert(GtkToggleButton, widgetptr)

    (data, strategies, canvas) = user_data

    data.playable = get_gtk_property(widget, "active", Bool)

    if data.playable
        playgame(data)
    end

    reveal(canvas)
    nothing
end

@guarded function restart_cb(widgetptr::Ptr, user_data)
    widget = convert(GtkButton, widgetptr)

    data, canvas, play_button = user_data

    data.movelist = []
    data.boardlist = [initialboard(guidatagame(data)())]

    ccall((:gtk_toggle_button_set_active , Gtk.libgtk), Ptr{GObject}, (Ptr{GObject}, Int),  play_button.handle, 0)

    reveal(canvas)
    nothing
end

@guarded function paramter_cb(widgetptr::Ptr, user_data)
    entry = convert(GtkEntry, widgetptr)

    strategy, name = user_data
end

@guarded function players_combobox_cb(widgetptr::Ptr, data)
    combobox = convert(GtkComboBoxText, widgetptr)

    if data.curr_visible_player != 0
        visible(data.playerboxes[data.curr_visible_player], false)
    end

    data.curr_visible_player = get_gtk_property(combobox, "active", Int)+1

    visible(data.playerboxes[data.curr_visible_player], true)

    nothing
end

@guarded function parameters_combobox_cb(widgetptr::Ptr, user_data)
    combobox = convert(GtkComboBoxText, widgetptr)

    (data, i) = user_data

    if data.curr_visible_strategy[i] != 0
        visible(data.strategyboxes[i][data.curr_visible_strategy[i]], false)
    end

    data.curr_visible_strategy[i] = get_gtk_property(combobox, "active", Int)+1

    visible(data.strategyboxes[i][data.curr_visible_strategy[i]], true)

    nothing
end

@guarded function set_parameters_cb(widgetptr::Ptr, user_data)
    data = user_data

    data.playerstrategies .= data.curr_visible_strategy

    #TODO:parameters

    for user in 1:length(data.playerboxes)
        for (sidx, strategy) in enumerate(data.allstrategies[user])
            for (pidx, parameter) in enumerate(getvarsnames(strategy))
                s = get_gtk_property(data.elemententries[user][sidx][pidx], :text, String)
                setvalue!(strategy, parameter, s)
            end
        end
    end

    nothing
end

function mainwindow(game::Type{<:Game}, strategies::Strategy...)
    n = nplayers(game)
    strategies = (UserStrategy(), strategies...)
    win = GtkWindow("A new window", 600, 600)

    mainhbox = GtkBox(:h)
    push!(win, mainhbox)

    w,h = defaultsize(game)
    canvas = GtkCanvas(w, h)
    push!(mainhbox, canvas)

    vbox = GtkBox(:v)
    push!(mainhbox, vbox)

    #side buttons
    restart = GtkButton("restart")
    play = GtkToggleButton("play")

    (
        strategy_frame,
        players_combobox,
        strategy_combobox,
        user_boxes,
        strategy_boxes,
        parameter_boxes,
        parameter_inputs
    ) = createframebox(n, strategies)

    set_parameters = GtkButton("set parameters")

    push!(vbox, strategy_frame)
    push!(vbox, set_parameters)
    push!(vbox, play)
    push!(vbox, restart)

    #visibility
    showall(win)

    for i in 1:length(user_boxes)
        visible(user_boxes[i], false)
        set_gtk_property!(strategy_combobox[i], :active, 0)
        for j in 2:length(strategy_boxes[i]) #show user
            visible(strategy_boxes[i][j], false)
        end
    end

    #callbacks
    data = GuiData(game, user_boxes, strategy_boxes, parameter_inputs, strategies)
    strategydict = Dict{String, Int}(zip(name.(strategies),1:length(strategies)))

    signal_connect(mouseclick_cb, canvas, "button_press_event", Nothing, (Ptr{Gtk.GdkEventButton},), false, (data, strategies))
    signal_connect(draw_cb, canvas, "draw", Nothing, (Ptr{Nothing},), false, data)
    signal_connect(play_cb, play, "toggled", Nothing, (), false, (data, strategies, canvas))
    signal_connect(restart_cb, restart, "clicked", Nothing, (), false, (data, canvas, play))
    signal_connect(players_combobox_cb, players_combobox, "changed", Nothing, (), false, data)
    for i in 1:n
        signal_connect(parameters_combobox_cb, strategy_combobox[i], "changed", Nothing, (), false, (data, i))
    end
    signal_connect(set_parameters_cb, set_parameters, "clicked", Nothing, (), false, data)
end

function createframebox(n, strategies)
    strategy_frame = GtkFrame("Strategies")
    strategy_frame_box = GtkBox(:v)
    push!(strategy_frame, strategy_frame_box)

    players_combobox = GtkComboBoxText()
    append!(players_combobox, "player " .* string.(1:n))
    push!(strategy_frame_box, players_combobox)

    strategy_combobox = [GtkComboBoxText() for i in 1:n]
    choices = map(name, strategies)
    for i in 1:n
        append!(strategy_combobox[i], choices)
    end

    parameter_inputs = Vector{Vector{GtkEntry}}[]
    parameter_names = Vector{Vector{GtkLabel}}[]

    user_boxes = GtkBox[]
    strategy_boxes = Vector{GtkBox}[]
    parameter_boxes = Vector{Vector{GtkBox}}[]

    for pidx in 1:n
        push!(parameter_inputs, [])
        push!(parameter_names, [])
        push!(strategy_boxes, [])
        push!(parameter_boxes, [])
        push!(user_boxes, GtkBox(:v))

        push!(user_boxes[pidx], strategy_combobox[pidx])

        for (sidx, strategy) in enumerate(strategies)
            push!(parameter_inputs[pidx], [])
            push!(parameter_names[pidx],  [])
            push!(parameter_boxes[pidx], [])
            push!(strategy_boxes[pidx], GtkBox(:v))

            parameters = getvalues(strategy)
            for (nidx, name) in enumerate(getvarsnames(strategy))
                push!(parameter_names[pidx][sidx], GtkLabel(name))

                push!(parameter_inputs[pidx][sidx], GtkEntry())
                set_gtk_property!(parameter_inputs[pidx][sidx][nidx], :text, string(parameters[nidx]))

                push!(parameter_boxes[pidx][sidx], GtkBox(:h))
                #put widgets in boxes
                push!(parameter_boxes[pidx][sidx][nidx], parameter_names[pidx][sidx][nidx])
                push!(parameter_boxes[pidx][sidx][nidx], parameter_inputs[pidx][sidx][nidx])
                push!(strategy_boxes[pidx][sidx], parameter_boxes[pidx][sidx][nidx])
            end
            push!(user_boxes[pidx], strategy_boxes[pidx][sidx])
        end
        push!(strategy_frame_box, user_boxes[pidx])
    end

    (
        strategy_frame,
        players_combobox,
        strategy_combobox,
        user_boxes,
        strategy_boxes,
        parameter_boxes,
        parameter_inputs
    )
end

function playgame(data::GuiData)
    b = data.boardlist[end]
    p = playerturn(b)
    sidx = data.playerstrategies[p]

    if sidx == 1 || isempty(getmoves(b))
        return
    else
        s = data.allstrategies[p][sidx]
        m = getmove(b, s)
        b = play(b, m)

        push!(data.movelist, m)
        push!(data.boardlist, b)

        #TODO: Reveal?
        playgame(data)
    end
end


text:
a, b, c = 1, 2, 3       # three globals
(function ()
  a = 10                # introduces local and shadows
  global b = 20         # overwrites global
  d = 10c + 10          # new local, reading global
  local e = 5
  while true
    e = 50              # outer e because marked local
    f = 60              # local to while loop!
    break
  end
  @assert (a,b,c,d,e) == (10,20,3,40,50)
  @assert (try f catch e sqrt(100) end) == 10
end)()
@assert (a,b,c) == (1,20,3)

ext: jl

max_stars_repo_name: idrougge/ple

max_issues_repo_name: idrougge/ple

max_forks_repo_name: idrougge/ple


text:
#=
Functions for the time counter in the GUI
Reads time from the FPGA and displays it in real time in seconds
=#
function _make_clock()

    s_label=Label("0")
    m_label=Label("0")
    h_label=Label("0")
    sm_label=Label(":")
    mh_label=Label(":")

    frame_time=Frame("Time Elapsed")
    hbox_time=ButtonBox(:h)
    push!(frame_time,hbox_time)

    push!(hbox_time,h_label)
    push!(hbox_time,mh_label)
    push!(hbox_time,m_label)
    push!(hbox_time,sm_label)
    push!(hbox_time,s_label)

    mytime_widgets=mytime(0,h_label,0,m_label,0,s_label)

    (mytime_widgets,frame_time)
end

function update_time(rhd::RHD2000,han::Gui_Handles)

    total_seconds=convert(Int64,div(rhd.time[1,1],rhd.sr))

    this_h=div(total_seconds,3600)

    total_seconds=total_seconds - this_h*3600

    this_m=div(total_seconds,60)

    this_s=total_seconds - this_m*60

    if this_s != han.time.s
        setproperty!(han.time.s_l,:label,string(this_s))
        han.time.s=this_s
    end
    if this_m != han.time.m
        setproperty!(han.time.m_l,:label,string(this_m))
        han.time.m=this_m
    end
    if this_h != han.time.h
        setproperty!(han.time.h_l,:label,string(this_h))
        han.time.h=this_h
    end

    nothing
end


text:
using ACEds
using Test

@testset "ACEds.jl" begin
    # Write your tests here.
end


text:
module CSVLog
export beginLog, logData
using CSV
using DataFrames

function beginLog(pathName::String, firstRow::Dict)
    CSV.write(pathName, DataFrame(firstRow))
end

function logData(pathName::String, row::Dict)
    log = CSV.write(
    pathName,
    DataFrame(row),
    append=true
    )
end
end


text:
Random.seed!(0)

################################################################################
### Regularized least squares, with two variable blocks to make things weird
################################################################################

println("Testing: regularized least squares, with two variable blocks to make things weird")

m, n1, n2 = 30, 50, 100

A1 = randn(m, n1)
A2 = randn(m, n2)
b = randn(m)

lam1 = 0.2
lam2 = 1.0

# Solve with FPG

x1_fpg = Variable(n1)
x2_fpg = Variable(n2)
expr = ls(A1*x1_fpg + A2*x2_fpg - b) + lam1*norm(x1_fpg, 1) + lam2*norm(x2_fpg, 2)
prob = problem(expr)
@time sol = solve(prob, StructuredOptimization.FPG(tol=1e-10,verbose=0,maxit=20000))
println(sol)

# Solve with ZeroFPR

x1_zerofpr = Variable(n1)
x2_zerofpr = Variable(n2)
expr = ls(A1*x1_zerofpr + A2*x2_zerofpr - b) + lam1*norm(x1_zerofpr, 1) + lam2*norm(x2_zerofpr, 2)
prob = problem(expr)
@time sol = solve(prob, StructuredOptimization.ZeroFPR(tol=1e-10,verbose=0))
println(sol)

# Solve with PANOC

x1_panoc = Variable(n1)
x2_panoc = Variable(n2)
expr = ls(A1*x1_panoc + A2*x2_panoc - b) + lam1*norm(x1_panoc, 1) + lam2*norm(x2_panoc, 2)
prob = problem(expr)
@time sol = solve(prob, StructuredOptimization.PANOC(tol=1e-10,verbose=0))
println(sol)

# Solve with minimize, use default solver/options

x1 = Variable(n1)
x2 = Variable(n2)
@time sol = @minimize ls(A1*x1 + A2*x2 - b) + lam1*norm(x1, 1) + lam2*norm(x2, 2)
println(sol)

@test norm(~x1_fpg - ~x1_zerofpr, Inf)/(1+norm(~x1_zerofpr, Inf)) <= 1e-6
@test norm(~x2_fpg - ~x2_zerofpr, Inf)/(1+norm(~x2_zerofpr, Inf)) <= 1e-6
@test norm(~x1_fpg - ~x1_panoc, Inf)/(1+norm(~x1_panoc, Inf)) <= 1e-6
@test norm(~x2_fpg - ~x2_panoc, Inf)/(1+norm(~x2_panoc, Inf)) <= 1e-6
@test norm(~x1 - ~x1_zerofpr, Inf)/(1+norm(~x1_zerofpr, Inf)) <= 1e-3
@test norm(~x2 - ~x2_zerofpr, Inf)/(1+norm(~x2_zerofpr, Inf)) <= 1e-3

res = A1*~x1_fpg + A2*~x2_fpg - b
grad1 = A1'*res
grad2 = A2'*res
ind1_zero = (~x1_fpg .== 0)
subgr1 = lam1*sign.(~x1_fpg)
subdiff1_low, subdiff1_upp = copy(subgr1), copy(subgr1)
subdiff1_low[ind1_zero] .= -lam1
subdiff1_upp[ind1_zero] .= +lam1
subgr2 = lam2*(~x2_fpg/norm(~x2_fpg, 2))

@test maximum(subdiff1_low + grad1) <= 1e-6
@test maximum(-subdiff1_upp - grad1) <= 1e-6
@test norm(grad2 + subgr2) <= 1e-6

###############################################################################
## Lasso problem with known solution
###############################################################################

println("Testing: lasso problem with known solution")

m, n, nnz_x_star = 200, 100, 10
A = randn(m, n)
lam = 1.0
x_star = randn(n)
x_star[nnz_x_star+1:end] .= 0.0
y_star = lam*sign.(x_star)
b = A*x_star + A'\y_star
@test norm(A'*(A*x_star - b) + lam*sign.(x_star)) <= 1e-12

# Solve with PG

x_pg = Variable(n)
expr = ls(A*x_pg - b) + lam*norm(x_pg, 1)
prob = problem(expr)
@time sol = solve(prob, StructuredOptimization.PG(tol=1e-10,verbose=0))
println(sol)

@test norm(~x_pg - x_star, Inf) <= 1e-8
@test norm(A'*(A*~x_pg - b) + lam*sign.(~x_pg)) <= 1e-6

# Solve with FPG

x_fpg = Variable(n)
expr = ls(A*x_fpg - b) + lam*norm(x_fpg, 1)
prob = problem(expr)
@time sol = solve(prob, StructuredOptimization.FPG(tol=1e-10,verbose=0))
println(sol)

@test norm(~x_fpg - x_star, Inf) <= 1e-8
@test norm(A'*(A*~x_fpg - b) + lam*sign.(~x_fpg)) <= 1e-6

# Solve with ZeroFPR

x_zerofpr = Variable(n)
expr = ls(A*x_zerofpr - b) + lam*norm(x_zerofpr, 1)
prob = problem(expr)
@time sol = solve(prob, StructuredOptimization.ZeroFPR(tol=1e-10,verbose=0))
println(sol)

@test norm(~x_zerofpr - x_star, Inf) <= 1e-8
@test norm(A'*(A*~x_zerofpr - b) + lam*sign.(~x_zerofpr)) <= 1e-5

# Solve with PANOC

x_panoc = Variable(n)
expr = ls(A*x_panoc - b) + lam*norm(x_panoc, 1)
prob = problem(expr)
@time sol = solve(prob, StructuredOptimization.PANOC(tol=1e-10,verbose=0))
println(sol)

@test norm(~x_panoc - x_star, Inf) <= 1e-8
@test norm(A'*(A*~x_panoc - b) + lam*sign.(~x_panoc)) <= 1e-5

################################################################################
### Problem with smooth, non-quadratic term
################################################################################

println("Testing: problem with smooth, non-quadratic term")

m, n, nnz_x_orig = 200, 500, 10
A = randn(m, n)
lam = 1.0
x_orig = randn(n)
x_orig[nnz_x_orig+1:end] .= 0.0
b = A*x_orig + randn(m)

# Solve with PG

x_pg = Variable(n)
expr = smooth(norm(A*x_pg - b, 2)) + lam*norm(x_pg, 1)
prob = problem(expr)
@time sol = solve(prob, PG(tol=1e-8,verbose=0))
println(sol)

# Solve with FPG

x_fpg = Variable(n)
expr = smooth(norm(A*x_fpg - b, 2)) + lam*norm(x_fpg, 1)
prob = problem(expr)
@time sol = solve(prob, FPG(tol=1e-8,verbose=0))
println(sol)

# Solve with ZeroFPR

x_zerofpr = Variable(n)
expr = smooth(norm(A*x_zerofpr - b, 2)) + lam*norm(x_zerofpr, 1)
prob = problem(expr)
@time sol = solve(prob, StructuredOptimization.ZeroFPR(tol=1e-8,verbose=0))
println(sol)

# Solve with PANOC

x_panoc = Variable(n)
expr = smooth(norm(A*x_panoc - b, 2)) + lam*norm(x_panoc, 1)
prob = problem(expr)
@time sol = solve(prob, StructuredOptimization.PANOC(tol=1e-8,verbose=0))
println(sol)

# Solve with minimize, default solver/options

x = Variable(n)
@time sol = @minimize smooth(norm(A*x - b, 2)) + lam*norm(x, 1)
println(sol)

@test norm(~x_fpg - ~x_zerofpr, Inf)/(1+norm(~x_zerofpr, Inf)) <= 1e-6
@test norm(~x_fpg - ~x_panoc, Inf)/(1+norm(~x_panoc, Inf)) <= 1e-6
@test norm(~x - ~x_zerofpr, Inf)/(1+norm(~x_zerofpr, Inf)) <= 1e-3

################################################################################
### Box-constrained least-squares
################################################################################

println("Testing: box-constrained least-squares")

m, n = 500, 200
A = randn(m, n)
lb, ub = -1.0, 1.0
x_orig = 2.0*randn(n)
x_orig = max.(lb, min.(ub, x_orig))
b = A*x_orig + randn(m)

# Solve with PG

x_pg = Variable(n)
expr = ls(A*x_pg - b)
prob = problem(expr, x_pg in [lb, ub])
@time sol = solve(prob, StructuredOptimization.PG(tol=1e-8,verbose=0))
println(sol)

@test norm(~x_pg - max.(lb, min.(ub, ~x_pg)), Inf) <= 1e-12
@test norm(~x_pg - max.(lb, min.(ub, ~x_pg - A'*(A*~x_pg - b))), Inf)/(1+norm(~x_pg, Inf)) <= 1e-8

# Solve with FPG

x_fpg = Variable(n)
expr = ls(A*x_fpg - b)
prob = problem(expr, x_fpg in [lb, ub])
@time sol = solve(prob, StructuredOptimization.FPG(tol=1e-8,verbose=0))
println(sol)

@test norm(~x_fpg - max.(lb, min.(ub, ~x_fpg)), Inf) <= 1e-12
@test norm(~x_fpg - max.(lb, min.(ub, ~x_fpg - A'*(A*~x_fpg - b))), Inf)/(1+norm(~x_fpg, Inf)) <= 1e-8

# Solve with ZeroFPR

x_zerofpr = Variable(n)
expr = ls(A*x_zerofpr - b)
prob = problem(expr, x_zerofpr in [lb, ub])
@time sol = solve(prob, StructuredOptimization.ZeroFPR(tol=1e-8,verbose=0))
println(sol)

@test norm(~x_zerofpr - max.(lb, min.(ub, ~x_zerofpr)), Inf) <= 1e-12
@test norm(~x_zerofpr - max.(lb, min.(ub, ~x_zerofpr - A'*(A*~x_zerofpr - b))), Inf)/(1+norm(~x_zerofpr, Inf)) <= 1e-8

# Solve with PANOC

x_panoc = Variable(n)
expr = ls(A*x_panoc - b)
prob = problem(expr, x_panoc in [lb, ub])
@time sol = solve(prob, StructuredOptimization.PANOC(tol=1e-8,verbose=0))
println(sol)

@test norm(~x_panoc - max.(lb, min.(ub, ~x_panoc)), Inf) <= 1e-12
@test norm(~x_panoc - max.(lb, min.(ub, ~x_panoc - A'*(A*~x_panoc - b))), Inf)/(1+norm(~x_panoc, Inf)) <= 1e-8

# Solve with minimize, default solver/options

x = Variable(n)
@time sol = @minimize ls(A*x - b) st x in [lb, ub]
println(sol)

@test norm(~x - max.(lb, min.(ub, ~x)), Inf) <= 1e-12
@test norm(~x - max.(lb, min.(ub, ~x - A'*(A*~x - b))), Inf)/(1+norm(~x, Inf)) <= 1e-4

################################################################################
### Non-negative least-squares from a known solution
################################################################################

println("Testing: non-negative least-squares from a known solution")

# Lagrangian:
#
#   0.5||Ax-b||^2 + y'(x-z) + [z >= 0]
#
# Optimality conditions:
#
#   A'(Ax-b) + y = 0, or A'b = A'Ax + y
#   x = z
#   z >= 0
#   y <= 0
#   y'z = 0

m, n, nnz_x_star = 500, 200, 100
A = randn(m, n)
x_star = rand(n)
x_star[nnz_x_star+1:end] .= 0.0
y_star = -rand(n)
y_star[1:nnz_x_star] .= 0.0
b = A*x_star + A'\y_star

# Solve with PG

x_pg = Variable(n)
expr = ls(A*x_pg - b)
prob = problem(expr, x_pg >= 0.0)
@time sol = solve(prob, StructuredOptimization.PG(tol=1e-8,verbose=0))
println(sol)

@test all(~x_pg .>= 0.0)
@test norm(~x_pg - x_star, Inf)/(1+norm(x_star, Inf)) <= 1e-8

# Solve with FPG

x_fpg = Variable(n)
expr = ls(A*x_fpg - b)
prob = problem(expr, x_fpg >= 0.0)
@time sol = solve(prob, StructuredOptimization.FPG(tol=1e-8,verbose=0))
println(sol)

@test all(~x_fpg .>= 0.0)
@test norm(~x_fpg - x_star, Inf)/(1+norm(x_star, Inf)) <= 1e-8

# Solve with ZeroFPR

x_zerofpr = Variable(n)
expr = ls(A*x_zerofpr - b)
prob = problem(expr, x_zerofpr >= 0.0)
@time sol = solve(prob, StructuredOptimization.ZeroFPR(tol=1e-8,verbose=0))
println(sol)

@test all(~x_zerofpr .>= 0.0)
@test norm(~x_zerofpr - x_star, Inf)/(1+norm(x_star, Inf)) <= 1e-8

# Solve with PANOC

x_panoc = Variable(n)
expr = ls(A*x_panoc - b)
prob = problem(expr, x_panoc >= 0.0)
@time sol = solve(prob, StructuredOptimization.PANOC(tol=1e-8,verbose=0))
println(sol)

@test all(~x_zerofpr .>= 0.0)
@test norm(~x_zerofpr - x_star, Inf)/(1+norm(x_star, Inf)) <= 1e-8

# Solve with minimize, default solver/options

x = Variable(n)
@time sol = @minimize ls(A*x - b) st x >= 0.0
println(sol)

@test all(~x .>= 0.0)
@test norm(~x - x_star, Inf)/(1+norm(x_star, Inf)) <= 1e-6


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

JLLWrappers.@generate_main_file_header("lib4ti2")
JLLWrappers.@generate_main_file("lib4ti2", UUID("1493ae25-0f90-5c0e-a06c-8c5077d6d66f"))
end  # module lib4ti2_jll


text:
export gate_count

function gate_count(spec)
    ci, _ = code_yao(Semantic.main, spec; optimize = true)
    return gate_count(ci)
end

function gate_count(ci::CodeInfo)
    count = Dict{Symbol,Any}()
    for stmt in ci.code
        is_quantum_statement(stmt) || continue
        qt = quantum_stmt_type(stmt)

        if qt === :gate || qt === :ctrl
            gc = get!(count, qt, IdDict{Any,Int}())
            if stmt.head === :invoke
                gate, gt = obtain_gate_stmt(stmt.args[3], ci)
            elseif stmt.head === :call
                gate, gt = obtain_gate_stmt(stmt.args[2], ci)
            end

            if gate isa IntrinsicRoutine && isempty(gate.variables)
                gc[gate] = get(gc, gate, 0) + 1
            else
                gc[parent(gate)] = get(gc, parent(gate), 0) + 1
            end
        else
            count[qt] = get(count, qt, 0) + 1
        end
    end
    return count
end


text:
function element_singlevariable_test(
    model::MOI.ModelLike,
    config::MOIT.Config,
)
    MOI.empty!(model)

    @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.Integer)
    @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.EqualTo{Int})
    @test MOI.supports_constraint(model, MOI.VectorOfVariables, CP.Element{Int})

    x1, _ = MOI.add_constrained_variable(model, MOI.Integer())
    x2, _ = MOI.add_constrained_variable(model, MOI.Integer())

    c1 = MOI.add_constraint(model, x1, MOI.EqualTo(1))

    c2 = MOI.add_constraint(model, MOI.VectorOfVariables([x1, x2]), CP.Element([6, 5, 4]))

    @test MOI.is_valid(model, x1)
    @test MOI.is_valid(model, x2)
    @test MOI.is_valid(model, c1)
    @test MOI.is_valid(model, c2)

    if config.solve
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

        @test MOI.get(model, MOI.ResultCount()) >= 1
        @test MOI.get(model, MOI.VariablePrimal(), x1) == 1
        @test MOI.get(model, MOI.VariablePrimal(), x2) == 5
    end
end

function element_scalaraffinefunction_test(
    model::MOI.ModelLike,
    config::MOIT.Config,
)
    MOI.empty!(model)

    @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.Integer)
    @test MOI.supports_constraint(
        model,
        MOI.ScalarAffineFunction{Int},
        MOI.EqualTo{Int},
    )
    @test MOI.supports_constraint(
        model,
        MOI.VectorAffineFunction{Int},
        CP.Element{Int},
    )

    x1, _ = MOI.add_constrained_variable(model, MOI.Integer())
    x2, _ = MOI.add_constrained_variable(model, MOI.Integer())

    c1 = MOI.add_constraint(model, 1 * x1, MOI.EqualTo(1))

    c2 = MOI.add_constraint(model, MOIU.vectorize([x1, x2]), CP.Element([6, 5, 4]))

    @test MOI.is_valid(model, x1)
    @test MOI.is_valid(model, x2)
    @test MOI.is_valid(model, c1)
    @test MOI.is_valid(model, c2)

    if config.solve
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

        @test MOI.get(model, MOI.ResultCount()) >= 1
        @test MOI.get(model, MOI.VariablePrimal(), x1) == 1
        @test MOI.get(model, MOI.VariablePrimal(), x2) == 5
    end
end

const elementtests = Dict(
    "element_singlevariable" => element_singlevariable_test,
    "element_scalaraffinefunction" => element_scalaraffinefunction_test,
)

MOIT.@moitestset element


text:
using Test
import Forecast: acf

@testset "ccf" begin
    x = rand(100);
    @test_throws AssertionError acf(x; lag = 100)

    res = acf(x; lag = 20, type="cor");
    @test res isa CCF
    @test res.call == "ccf(x; type=\"cor\", lag=20, alpha=(0.95, 0.99))"
    @test res.type == "cor"
    @test res.auto
    @test res.lag == 20
    @test res.N == 100
    @test res.alpha == (0.95, 0.99)
    @test res.ci[1] ≈ 0.19900419155027554
    @test res.ci[2] ≈ 0.2615358405397182
end


text:
using JuLIP
using Base.Test
using JuLIP.Testing

verbose=true

julip_tests = [
   "testaux.jl";
   "testase.jl";
   "testanalyticpotential.jl";
   "testpotentials.jl";
   "testsolve.jl";
]

println("Starting JuLIP Tests")
println("=====================")

for test in julip_tests
   include(test)
end


text:
# all the plots for my presentation in IOR2017, Stavanger
using Plots
# include("rel_perms.jl")

# plot the fractional flow functions for a water-wet and an oil-wet system
include("rel_perms_real.jl")
# Water-wet Mixed-wet Oil-wet
# Swi 0.2 0.2 0.2
# Sor 0.1 0.2 0.3
# Ew=Krw(Sor) 0.6 0.8 0.8
# Eo=Kro(Swi) 0.8 0.8 0.5
# Nw 3 2 3
# No 2 2 3
krw0_ww = 0.33
krw0_ow = 0.28
kro0_ww = 0.79
kro0_ow = 0.72
nw_ww = 2.4
nw_ow= 2.4
no_ww = 2.0
no_ow= 1.6
sor_ww=0.145
sor_ow=0.37
sor_mw=0.08
swc_ww=0.08
swc_ow=0.08
teta_ow=deg2rad(180)
teta_ww=deg2rad(0)
teta_mw=deg2rad(90)
teta_prime=teta_mw-π/2 # see my report
a=(sor_mw-0.5*(sor_ww-sor_ow)*cos(teta_mw)-0.5*(sor_ww+sor_ow))/(cos(teta_mw)^2-1)
b=0.5*(sor_ww-sor_ow)
c=0.5*(sor_ww+sor_ow)-a

SF=0.0 # 1.0 is water-wet, 0.0 is oil-wet
sw_plot_ow=linspace(swc_ow,1-sor_ow,100)
sw_plot_ww=linspace(swc_ww,1-sor_ww,100)

plot(size=(500,400), xtickfont = font(12, "Courier"), ytickfont=font(12, "Courier"),
  ylabel="krw, kro", xlabel="Water saturation", legendfont=font(12, "Courier"),
  guidefont=font(12, "Courier"))
plot!(sw_plot_ow, krw.(sw_plot_ow, krw0_ow, sor_ow, swc_ow, nw_ow),
  xlims=(0.0,1.0), linewidth=3, linestyle=:dash, linecolor=:red,
  label="krw, Oil-wet")
plot!(sw_plot_ow, kro.(sw_plot_ow, kro0_ow, sor_ow, swc_ow, no_ow),
  xlims=(0.0,1.0), linewidth=3, linestyle=:dash, linecolor=:red,
  label="kro, Oil-wet")
plot!(sw_plot_ww, krw.(sw_plot_ww, krw0_ww, sor_ww, swc_ww, nw_ww),
  xlims=(0.0,1.0), linewidth=1, linestyle=:solid, linecolor=:blue,
  label="krw, Water-wet")
plot!(sw_plot_ww, kro.(sw_plot_ww, kro0_ww, sor_ww, swc_ww, no_ww),
  xlims=(0.0,1.0), linewidth=1, linestyle=:solid, linecolor=:blue,
  label="kro, Water-wet")
savefig("rel_perms_fathi.png")

max_stars_repo_name: simulkade/peteng

max_issues_repo_name: simulkade/peteng

max_forks_repo_name: simulkade/peteng


text:
using DataFrames, LowRankModels

println("cross validation example")
srand(5)

do_cv = true
do_cv_by_iter = true
do_reg_path = true
do_plot = false

m,n,k = 50,50,3
A = randn(m,k)*randn(k,n) + k*sprandn(m,n,.05)
losses = fill(HuberLoss(),n)
r = QuadReg(.1)
glrm = GLRM(A,losses,r,r,k+2)

if do_cv
    println("Computing cross validation error for each fold")
    params = Params(1.0, max_iter=100, abs_tol=0.0, min_stepsize=.001)
    train_error, test_error, train_glrms, test_glrms = cross_validate(glrm, nfolds=5, params=params)
    df = DataFrame(train_error = train_error, test_error = test_error)
end

if do_cv_by_iter
    println("Computing training and testing error as a function of iteration number")
    train_error, test_error = cv_by_iter(glrm)
    df = DataFrame(train_error = train_error, test_error = test_error)
end

if do_reg_path
    println("Computing regularization path")
    params = Params(1.0, max_iter=50, abs_tol=.00001, min_stepsize=.01)
    train_error, test_error, train_time, reg_params =
    regularization_path(glrm, params=params, reg_params=logspace(2,-2,15))
    df = DataFrame(train_error = train_error, test_error = test_error,
                   train_time = train_time, reg_param = reg_params)
	if do_plot
		p = plot(df, :reg_param, [:train_error, :test_error]; scale = :log, filename = None)
	end
end

println(df)


text:
# We reverse the order of comparisons here so that the result
# of x < y is equal to the result of Monomial(x) < Monomial(y)
Base.@pure Base.isless(v1::AbstractVariable, v2::AbstractVariable) = name(v1) > name(v2)
Base.isless(m1::AbstractTermLike, m2::AbstractTermLike) = isless(promote(m1, m2)...)

function Base.isless(t1::AbstractTerm, t2::AbstractTerm)
    if monomial(t1) < monomial(t2)
        true
    elseif monomial(t1) == monomial(t2)
        coefficient(t1) < coefficient(t2)
    else
        false
    end
end

# promoting multiplication is not a good idea
# For example a polynomial of Float64 * a polynomial of JuMP affine expression
# is a polynomial of JuMP affine expression but if we promote it would be a
# polynomial of quadratic expression
for op in [:+, :-, :(==)]
    @eval Base.$op(p1::APL, p2::APL) = $op(promote(p1, p2)...)
end
Base.isapprox(p1::APL, p2::APL; kwargs...) = isapprox(promote(p1, p2)...; kwargs...)

# @eval $op(p::APL, α) = $op(promote(p, α)...) would be less efficient
for (op, fun) in [(:+, :plusconstant), (:-, :minusconstant), (:*, :multconstant), (:(==), :eqconstant)]
    @eval Base.$op(p::APL, α) = $fun(p, α)
    @eval Base.$op(α, p::APL) = $fun(α, p)
end

# Special case AbstractArrays of APLs
# We add these instead of relying on the broadcasting API since the above method definitions are very wide.
# In particular, there is support for Matrices as coefficents. In order to avoid issues like #104 we therefore
# explicitly define this (instead of implictly getting unexpected results).
for op in [:+, :-, :*]
    @eval Base.$op(p::APL, A::AbstractArray{<:APL}) = map(f -> $op(p, f), A)
    @eval Base.$op(A::AbstractArray{<:APL}, p::APL) = map(f -> $op(f, p), A)
end
Base.:/(A::AbstractArray{<:APL}, p::APL) = map(f -> f / p, A)

Base.isapprox(p::APL, α; kwargs...) = isapprox(promote(p, α)...; kwargs...)
Base.isapprox(α, p::APL; kwargs...) = isapprox(promote(p, α)...; kwargs...)

Base.:-(m::AbstractMonomialLike) = (-1) * m
Base.:-(t::AbstractTermLike) = (-coefficient(t)) * monomial(t)
Base.:-(p::APL) = polynomial((-).(terms(p)))
Base.:+(p::Union{APL, RationalPoly}) = p
Base.:*(p::Union{APL, RationalPoly}) = p

# Avoid adding a zero constant that might artificially increase the Newton polytope
# Need to add polynomial conversion for type stability
plusconstant(p::APL{S}, α::T)  where {S, T} = iszero(α) ? polynomial( p, Base.promote_op(+, S, T)) : p + constantterm(α, p)
plusconstant(α::S, p::APL{T})  where {S, T} = iszero(α) ? polynomial( p, Base.promote_op(+, S, T)) : constantterm(α, p) + p
minusconstant(p::APL{S}, α::T) where {S, T} = iszero(α) ? polynomial( p, Base.promote_op(-, S, T)) : p - constantterm(α, p)
minusconstant(α::S, p::APL{T}) where {S, T} = iszero(α) ? polynomial(-p, Base.promote_op(-, S, T)) : constantterm(α, p) - p

# Coefficients and variables commute
multconstant(α, v::AbstractVariable) = multconstant(α, monomial(v)) # TODO linear term
multconstant(m::AbstractMonomialLike, α) = multconstant(α, m)

_multconstant(α, f, t::AbstractTermLike) = mapcoefficientsnz(f, t)
function _multconstant(α::T, f, p::AbstractPolynomial{S}) where {S, T}
    if iszero(α)
        zero(polynomialtype(p, Base.promote_op(*, T, S)))
    else
        mapcoefficientsnz(f, p)
    end
end
_multconstant(α, f, p::AbstractPolynomialLike) = _multconstant(α, f, polynomial(p))

multconstant(α, p::AbstractPolynomialLike) = _multconstant(α, β -> α*β, p)
multconstant(p::AbstractPolynomialLike, α) = _multconstant(α, β -> β*α, p)

Base.:*(m1::AbstractMonomialLike, m2::AbstractMonomialLike) = mapexponents(+, m1, m2)
#Base.:*(m1::AbstractMonomialLike, m2::AbstractMonomialLike) = *(monomial(m1), monomial(m2))

Base.:*(m::AbstractMonomialLike, t::AbstractTermLike) = coefficient(t) * (m * monomial(t))
Base.:*(t::AbstractTermLike, m::AbstractMonomialLike) = coefficient(t) * (monomial(t) * m)
Base.:*(t1::AbstractTermLike, t2::AbstractTermLike) = (coefficient(t1) * coefficient(t2)) * (monomial(t1) * monomial(t2))

Base.:*(t::AbstractTermLike, p::APL) = polynomial(map(te -> t * te, terms(p)))
Base.:*(p::APL, t::AbstractTermLike) = polynomial(map(te -> te * t, terms(p)))
Base.:*(p::APL, q::APL) = polynomial(p) * polynomial(q)

# guaranteed that monomial(t1) > monomial(t2)
function _polynomial_2terms(t1::TT, t2::TT, ::Type{T}) where {TT<:AbstractTerm, T}
    if iszero(t1)
        polynomial(t2, T)
    elseif iszero(t2)
        polynomial(t1, T)
    else
        polynomial(termtype(TT, T)[t1, t2], SortedUniqState())
    end
end
for op in [:+, :-]
    @eval begin
        Base.$op(t1::AbstractTermLike, t2::AbstractTermLike) = $op(term(t1), term(t2))
        Base.$op(t1::AbstractTerm, t2::AbstractTerm) = $op(promote(t1, t2)...)
        function Base.$op(t1::TT, t2::TT) where {T, TT <: AbstractTerm{T}}
            S = Base.promote_op($op, T, T)
            # t1 > t2 would compare the coefficient in case the monomials are equal
            # and it will throw a MethodError in case the coefficients are not comparable
            if monomial(t1) == monomial(t2)
                polynomial($op(coefficient(t1), coefficient(t2)) * monomial(t1), S)
            elseif monomial(t1) > monomial(t2)
                ts = _polynomial_2terms(t1, $op(t2), S)
            else
                ts = _polynomial_2terms($op(t2), t1, S)
            end
        end
    end
end

LinearAlgebra.adjoint(v::AbstractVariable) = v
LinearAlgebra.adjoint(m::AbstractMonomial) = m
LinearAlgebra.adjoint(t::AbstractTerm) = LinearAlgebra.adjoint(coefficient(t)) * monomial(t)
LinearAlgebra.adjoint(p::AbstractPolynomialLike) = polynomial(map(LinearAlgebra.adjoint, terms(p)))
LinearAlgebra.adjoint(r::RationalPoly) = adjoint(numerator(r)) / adjoint(denominator(r))

LinearAlgebra.transpose(v::AbstractVariable) = v
LinearAlgebra.transpose(m::AbstractMonomial) = m
LinearAlgebra.transpose(t::AbstractTerm) = LinearAlgebra.transpose(coefficient(t)) * monomial(t)
LinearAlgebra.transpose(p::AbstractPolynomialLike) = polynomial(map(LinearAlgebra.transpose, terms(p)))
LinearAlgebra.transpose(r::RationalPoly) = transpose(numerator(r)) / transpose(denominator(r))

LinearAlgebra.dot(p1::AbstractPolynomialLike, p2::AbstractPolynomialLike) = p1' * p2
LinearAlgebra.dot(x, p::AbstractPolynomialLike) = x' * p
LinearAlgebra.dot(p::AbstractPolynomialLike, x) = p' * x

# Amazingly, this works! Thanks, StaticArrays.jl!
"""
Convert a tuple of variables into a static vector to allow array-like usage.
The element type of the vector will be Monomial{vars, length(vars)}.
"""
Base.vec(vars::Tuple{Vararg{AbstractVariable}}) = [vars...]
# vec(vars::Tuple{Vararg{<:TypedVariable}}) = SVector(vars)

# https://github.com/JuliaLang/julia/pull/23332
Base.:^(x::AbstractPolynomialLike, p::Integer) = Base.power_by_squaring(x, p)


text:
export DecomposedDiscretePost

"""
    DecomposedDiscretePost <: DiscretePost

Textbook implementation of a discrete post operator, but with decomposed intersections.

### Fields

- `options` -- an `Options` structure that holds the algorithm-specific options

### Algorithm

The algorithm is based on [Reachability analysis of linear hybrid systems via block
decomposition](https://arxiv.org/pdf/1905.02458.pdf).
This discrete-post operator can only be used in combination with the continuous-post operator [`BFFPS19`](@ref).
"""
struct DecomposedDiscretePost <: DiscretePost
    options::Options

    function DecomposedDiscretePost(𝑂::Options)
        𝑂copy = copy(𝑂)
        # TODO: Check why it takes always default value for convex_hull
        check_aliases_and_add_default_value!(𝑂.dict, 𝑂copy.dict, [:overapproximation], Hyperrectangle)
        check_aliases_and_add_default_value!(𝑂.dict, 𝑂copy.dict, [:out_vars], Vector{Int}())

        return new(𝑂copy)
    end
end

# convenience constructor from pairs of symbols
DecomposedDiscretePost(𝑂::Pair{Symbol,<:Any}...) = DecomposedDiscretePost(Options(Dict{Symbol,Any}(𝑂)))

# default options for the DecomposedDiscretePost discrete post operator
DecomposedDiscretePost() = DecomposedDiscretePost(Options())

init(𝒫::DecomposedDiscretePost, 𝒮::AbstractSystem, 𝑂::Options) = init!(𝒫, 𝒮, copy(𝑂))

# TODO: use 𝑂 only?
function init!(𝒫::DecomposedDiscretePost, 𝒮::AbstractSystem, 𝑂::Options)
    𝑂[:n] = statedim(𝒮, 1)

    # solver-specific options (adds default values for unspecified options)
    𝑂out = validate_solver_options_and_add_default_values!(𝑂)

    return 𝑂out
end

function tube⋂inv(𝒫::DecomposedDiscretePost,
                  reach_tube::Vector{<:SparseReachSet{<:LazySet{N}}},
                  invariant,
                  start_interval
                 ) where {N}

    Rsets = Vector{AbstractReachSet{<:LazySet{N}}}()
    @inbounds for reach_set in reach_tube
        # intersection with invariant is computed inside BFFPS19 CPost operator
        push!(Rsets,
              substitute(reach_set,
                         time_start=time_start(reach_set) + start_interval[1],
                         time_end=time_end(reach_set) + start_interval[2]))
    end

    return Rsets
end

function post(𝒫::DecomposedDiscretePost,
              HS::HybridSystem,
              waiting_list::Vector{Tuple{Int, <:AbstractReachSet{LazySet{N}}, Int}},
              passed_list,
              source_loc_id,
              tube⋂inv,
              jumps,
              options
             ) where {N}
    jumps += 1
    n = options[:n]
    oa = 𝒫.options[:overapproximation]
    n_lowdim = length(𝒫.options[:temp_vars])
    source_invariant = HS.modes[source_loc_id].X
    inv_isa_Hrep, inv_isa_H_polytope = get_Hrep_info(source_invariant)
    for trans in out_transitions(HS, source_loc_id)
        info("Considering transition: $trans")
        target_loc_id = target(HS, trans)
        target_loc = HS.modes[target(HS, trans)]
        target_invariant = target_loc.X
        constrained_map = resetmap(HS, trans)
        guard = stateset(constrained_map)
        # perform jumps
        post_jump = Vector{ReachSet{LazySet{N}}}()
        sizehint!(post_jump, length(tube⋂inv))
        for reach_set in tube⋂inv
            if (dim(set(reach_set)) == n_lowdim && n_lowdim < n)
                continue
            end
            # check intersection with guard
            R⋂G = Intersection(set(reach_set), guard)
            if isempty(R⋂G)
                continue
            end
            R⋂G = overapproximate(R⋂G, CartesianProductArray, oa)

            # apply assignment
            A⌜R⋂G⌟ = apply_assignment(𝒫, constrained_map, R⋂G)
            A⌜R⋂G⌟ = overapproximate(A⌜R⋂G⌟, CartesianProductArray, oa)

            # intersect with target invariant
            A⌜R⋂G⌟⋂I = Intersection(target_invariant, A⌜R⋂G⌟)
            if isempty(A⌜R⋂G⌟⋂I)
                continue
            end

            A⌜R⋂G⌟⋂I = overapproximate(A⌜R⋂G⌟⋂I, CartesianProductArray, oa)

            # store result
            push!(post_jump, ReachSet{LazySet{N}}(A⌜R⋂G⌟⋂I,
                                                  time_start(reach_set),
                                                  time_end(reach_set)))
        end

        postprocess(𝒫, HS, post_jump, options, waiting_list, passed_list,
            target_loc_id, jumps)
    end
end


text:
using Distributed
using SharedArrays
using ArgParse

@everywhere using LatentOrderBookModel

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table s begin
        "replications"
            arg_type = Int
            default = 2
        "surface_points"
            arg_type = Int
            default = 2
        "param1_name"
            arg_type = String
            default = "D"
        "param1_value_lower"
            arg_type = Float64
            default = 1.0
        "param1_value_upper"
            arg_type = Float64
            default = 5.0
        "param2_name"
            arg_type = String
            default = "σ"
        "param2_value_lower"
            arg_type = Float64
            default = 0.1
        "param2_value_upper"
            arg_type = Float64
            default = 4.0
    end
    return parse_args(s)
end

function main()
    parsed_args = parse_commandline()
    cpus_allocated = nworkers()
    println("$cpus_allocated cores allocated for job")

    replications = parsed_args["replications"]
    surface_points = parsed_args["surface_points"]
    println("Number of price paths generated: $(surface_points^2*replications)")
    println("Surface points: $surface_points and Replications: $replications")
    param1_name = parsed_args["param1_name"]
    param1_range = [parsed_args["param1_value_lower"], parsed_args["param1_value_upper"]]
    param1_values = repeat(repeat(collect(range(param1_range..., length=surface_points)),inner=surface_points),replications)

    param2_name = parsed_args["param2_name"]
    param2_range = [parsed_args["param2_value_lower"], parsed_args["param2_value_upper"]]
    param2_values = repeat(repeat(collect(range(param2_range..., length=surface_points)),outer=surface_points),replications)
    println("Creating surface plot for parameter $param1_name [$(param1_range[1]), $(param1_range[2])] and parameter $param2_name [$(param2_range[1]), $(param2_range[2])]")

    fixed_params=Dict()
    fixed_params["num_paths"] = replications
    fixed_params["T"]=2299
    fixed_params["p₀"]=238.745
    fixed_params["M"]=250
    fixed_params["β"]=0.1
    fixed_params["L"]=150
    fixed_params["σ"]=2.0
    fixed_params["D"]=5.0
    fixed_params["nu"]=0.0
    fixed_params["λ"]=1.0
    fixed_params["μ"]=0.002

    objective_surface = ObjectiveSurface(
        param1_name, param1_values,
        param2_name, param2_values,
        surface_points, fixed_params)

    random_seed = Int(rand(UInt32))
    price_paths = @time objective_surface(random_seed, true)
    println("Random seed is $random_seed")
    println(size(price_paths))
    io = open("$param1_name" * "_vs_" * "$param2_name" * "_price_paths.csv", "w")
    for i in 1:size(price_paths,1)
        write(io, join(Array(price_paths[i,:]), ", "), "\n")
    end
    close(io)
end

if "" != PROGRAM_FILE && occursin(PROGRAM_FILE, @__FILE__)
    main()
end

max_stars_repo_name: GantZA/UCT-HPC-Resources

max_issues_repo_name: GantZA/UCT-HPC-Resources

max_forks_repo_name: GantZA/UCT-HPC-Resources


text:
#--------------------------------------------------------------------
# Spacetime Discretization methods in Julia
# Soham 01-2018
#--------------------------------------------------------------------

function testextractPatchCoeffs(Nx::Int, Ny::Int)::Float64
    coeffs  = randn(Nx+1,Ny+1)
    vndmx   = vandermonde(Nx)
    vndmy   = vandermonde(Ny)
    patch   = Patch([1,1], vndmx'*coeffs*vndmy)
    ncoeffs = extractPatchCoeffs(patch)
    return L2norm(ncoeffs, coeffs, chebweights(Nx), chebweights(Ny))
end

function testinterpolatePatch(Nx::Int, Ny::Int)::Float64
    coeffs  = randn(Nx+1,Ny+1)
    vndmx   = vandermonde(Nx)
    vndmy   = vandermonde(Ny)
    patch   = Patch([1,1], vndmx'*coeffs*vndmy)
    npatch  = interpolatePatch(patch, Nx, Ny)
    return L2norm(npatch.value, patch.value, chebweights(Nx), chebweights(Ny))
end

function testinterpolatePatch(Nx::Int, Ny::Int, fn::Function)::Float64
    fpatch = Float64[fn(x,y) for x in chebgrid(Nx), y in chebgrid(Ny)]
    patch  = Patch([1,1], fpatch)
    npatch = interpolatePatch(patch, 2*Nx, 2*Ny)
    fpatch2N = Float64[fn(x,y) for x in chebgrid(2*Nx), y in chebgrid(2*Ny)]
    return L2norm(npatch.value, fpatch2N, chebweights(2*Nx), chebweights(2*Ny))
end

function testcalcPatch(bnd1fn::Function, bnd2fn::Function, rhsfn::Function, Nx::Int, Ny::Int)::Float64
    bndx = getPatchIC(bnd1fn, 0, Ny, 1, 1)
    bndy = getPatchIC(bnd2fn, 1, Ny, 1, 1)
    rhs  = RHS(rhsfn, Nx, Ny, 1, [1,1])
    dop  = derivOP(Nx,Ny)
    bop  = boundaryOP(Nx,Ny)
    npatch = calcPatch(bndx, bndy, rhs, dop, bop, [1,1])
    fpatch = Float64[bnd1fn(x) + bnd1fn(y) for x in chebgrid(Nx), y in chebgrid(Ny)]
    return L2norm(npatch.value, fpatch, chebgrid(Nx), chebgrid(Ny))
end

fpatch = Float64[(x-1)^2 + (y-1)^2 for x in 1:10, y in 1:15]

@test getPatchIC(x->x, 0, 4, 2, 1).value ≈ (chebgrid(4) + 1)/2
@test getPatchIC(x->x, 1, 6, 2, 1).value ≈ (chebgrid(6) + 1)/2
@test getPatchIC(x->x, 1, 9, 2, 2).value ≈ (chebgrid(9) - 1)/2
@test getPatchBnd(Patch([1,1], fpatch), 0).value == Float64[(x-1)^2 + 14^2 for x in 1:10]
@test getPatchBnd(Patch([1,1], fpatch), 1).value == Float64[(y-1)^2 + 9^2 for y in 1:15]
@test testextractPatchCoeffs(7,3) < 1e-14
@test testinterpolatePatch(12,5) < 1e-14
@test testinterpolatePatch(12, 13, (x,y)->x^2 + y^3 + y^2*x^3) < 1e-14
@test testinterpolatePatch(29, 31, (x,y)-> sin(pi*x) + exp(-y^2)) < 1e-14
@test testcalcPatch(x->(x^9-1), y->(y^9-1), (x,y)->0, 9, 9) < 1e-13


text:
"Generate a Jupyter notebook for a given exercise."
function generate_notebook(slug)
    # readlines returns an array of strings.
    # Jupyter treats arrays of strings as multiline strings
    # but does not insert \n between elements of the array,
    # so it must be appended to each line manually.
    readme = readlines(joinpath("exercises", "$slug", "README.md")) .* "\n"
    stub = readlines(joinpath("exercises", "$slug", "$slug.jl")) .* "\n"
    tests = readlines(joinpath("exercises", "$slug", "runtests.jl")) .* "\n"

    # mark stub cell as to-be-submitted
    pushfirst!(stub, "# submit\n")

    # remove trailing newline at the end of each cell
    readme[end] = readme[end][1:end - 1]
    tests[end] = tests[end][1:end - 1]
    # not all exercises have stubs
    if !isempty(stub)
        stub[end] = stub[end][1:end - 1]
    end

    # remove include($stub.jl) from the tests
    for (i, line) in enumerate(tests)
        tests[i] = replace(line, "include(\"$slug.jl\")" => "# include(\"$slug.jl\")")
    end

    """
    {
        "cells": [
            {
                "cell_type": "markdown",
                "metadata": {},
                "source": $readme
            },
            {
                "cell_type": "markdown",
                "metadata": {},
                "source": ["## Your solution"]
            },
            {
                "cell_type": "code",
                "execution_count": null,
                "metadata": {},
                "outputs": [],
                "source": $stub
            },
            {
                "cell_type": "markdown",
                "metadata": {},
                "source": ["## Test suite"]
            },
            {
                "cell_type": "code",
                "execution_count": null,
                "metadata": {},
                "outputs": [],
                "source": $tests
            },
            {
                "cell_type": "markdown",
                "metadata": {},
                "source": [
                    "## Prepare submission\\n",
                    "To submit your exercise, you need to save your solution in a file called `$slug.jl` before using the CLI.\\n",
                    "You can either create it manually or use the following functions, which will automatically write every notebook cell that starts with `# submit` to the file `$slug.jl`.\\n"
                ]
            },
            {
                "cell_type": "code",
                "execution_count": null,
                "metadata": {},
                "outputs": [],
                "source": [
                    "# using Pkg; Pkg.add(\\\"Exercism\\\")\\n",
                    "# using Exercism\\n",
                    "# Exercism.create_submission(\\\"$slug\\\")"
                ]
            }
        ],
        "metadata": {
            "kernelspec": {
                "display_name": "Julia 1.3.0",
                "language": "julia",
                "name": "julia-1.3"
            },
            "language_info": {
                "file_extension": ".jl",
                "mimetype": "application/julia",
                "name": "julia",
                "version": "1.3.0"
            }
        },
        "nbformat": 4,
        "nbformat_minor": 2
    }
    """
end

function generate_notebooks()
    @info "Generating notebooks..."
    for slug in readdir("exercises")
        @debug "  $slug"
        outfile = joinpath("exercises", slug, "$slug.ipynb")
        write(outfile, generate_notebook(slug))
    end
    @info "Notebooks successfully generated"
end

generate_notebooks()


text:
function QuadraticCost(Q, R)
    return function (x, u)
        x'*Q*x + u'*R*u
    end
end


text:
using Oceananigans.Solvers
using Oceananigans.Operators
using Oceananigans.Architectures
using Oceananigans.Fields: Field, ZReducedField

import Oceananigans.Solvers: solve!
import Oceananigans.Architectures: architecture

struct PCGImplicitFreeSurfaceSolver{V, S, R}
    vertically_integrated_lateral_areas :: V
    preconditioned_conjugate_gradient_solver :: S
    right_hand_side :: R
end

architecture(solver::PCGImplicitFreeSurfaceSolver) =
    architecture(solver.preconditioned_conjugate_gradient_solver)

"""
    PCGImplicitFreeSurfaceSolver(grid, settings)

Return a solver based on a preconditioned conjugate gradient method for the elliptic equation

```math
[∇ ⋅ H ∇ - Az / (g Δt²)] ηⁿ⁺¹ = (∇ʰ ⋅ Q★ - Az ηⁿ / Δt) / (g Δt)
```

representing an implicit time discretization of the linear free surface evolution equation
for a fluid with variable depth `H`, horizontal areas `Az`, barotropic volume flux `Q★`, time
step `Δt`, gravitational acceleration `g`, and free surface at time-step `n` `ηⁿ`.
"""
function PCGImplicitFreeSurfaceSolver(grid::AbstractGrid, gravitational_acceleration::Number, settings)
    # Initialize vertically integrated lateral face areas
    ∫ᶻ_Axᶠᶜᶜ = Field{Face, Center, Nothing}(grid)
    ∫ᶻ_Ayᶜᶠᶜ = Field{Center, Face, Nothing}(grid)

    vertically_integrated_lateral_areas = (xᶠᶜᶜ = ∫ᶻ_Axᶠᶜᶜ, yᶜᶠᶜ = ∫ᶻ_Ayᶜᶠᶜ)

    compute_vertically_integrated_lateral_areas!(vertically_integrated_lateral_areas, grid)

    right_hand_side = Field{Center, Center, Nothing}(grid)

    # Set maximum iterations to Nx * Ny if not set
    settings = Dict{Symbol, Any}(settings)
    maximum_iterations = get(settings, :maximum_iterations, grid.Nx * grid.Ny)
    settings[:maximum_iterations] = maximum_iterations

    solver = PreconditionedConjugateGradientSolver(implicit_free_surface_linear_operation!;
                                                   template_field = right_hand_side,
                                                   settings...)

    return PCGImplicitFreeSurfaceSolver(vertically_integrated_lateral_areas, solver, right_hand_side)
end

build_implicit_step_solver(::Val{:PreconditionedConjugateGradient}, grid, gravitational_acceleration, settings) =
    PCGImplicitFreeSurfaceSolver(grid, gravitational_acceleration, settings)

#####
##### Solve...
#####

function solve!(η, implicit_free_surface_solver::PCGImplicitFreeSurfaceSolver, rhs, g, Δt)
    #=
    # Somehow take an explicit step first to produce an improved initial guess
    # for η for the iterative solver.
    event = explicit_ab2_step_free_surface!(free_surface, model, Δt, χ)
    wait(device(model.architecture), event)
    =#

    ∫ᶻA = implicit_free_surface_solver.vertically_integrated_lateral_areas
    solver = implicit_free_surface_solver.preconditioned_conjugate_gradient_solver

    # solve!(x, solver, b, args...) solves A*x = b for x.
    solve!(η, solver, rhs, ∫ᶻA.xᶠᶜᶜ, ∫ᶻA.yᶜᶠᶜ, g, Δt)

    return nothing
end

function compute_implicit_free_surface_right_hand_side!(rhs,
                                                        implicit_solver::PCGImplicitFreeSurfaceSolver,
                                                        g, Δt, ∫ᶻQ, η)

    solver = implicit_solver.preconditioned_conjugate_gradient_solver
    arch = architecture(solver)
    grid = solver.grid

    event = launch!(arch, grid, :xy,
                    implicit_free_surface_right_hand_side!,
                    rhs, grid, g, Δt, ∫ᶻQ, η,
		            dependencies = device_event(arch))

    return event
end

""" Compute the divergence of fluxes Qu and Qv. """
@inline flux_div_xyᶜᶜᵃ(i, j, k, grid, Qu, Qv) = δxᶜᵃᵃ(i, j, k, grid, Qu) + δyᵃᶜᵃ(i, j, k, grid, Qv)

@kernel function implicit_free_surface_right_hand_side!(rhs, grid, g, Δt, ∫ᶻQ, η)
    i, j = @index(Global, NTuple)
    Az = Azᶜᶜᵃ(i, j, 1, grid)
    δ_Q = flux_div_xyᶜᶜᵃ(i, j, 1, grid, ∫ᶻQ.u, ∫ᶻQ.v)
    @inbounds rhs[i, j, 1] = (δ_Q - Az * η[i, j, 1] / Δt) / (g * Δt)
end

"""
Returns `L(ηⁿ)`, where `ηⁿ` is the free surface displacement at time step `n`
and `L` is the linear operator that arises
in an implicit time step for the free surface displacement `η`.

(See the docs section on implicit time stepping.)
"""
function implicit_free_surface_linear_operation!(L_ηⁿ⁺¹, ηⁿ⁺¹, ∫ᶻ_Axᶠᶜᶜ, ∫ᶻ_Ayᶜᶠᶜ, g, Δt)
    grid = L_ηⁿ⁺¹.grid
    arch = architecture(L_ηⁿ⁺¹)

    fill_halo_regions!(ηⁿ⁺¹, arch)

    event = launch!(arch, grid, :xy, _implicit_free_surface_linear_operation!,
                    L_ηⁿ⁺¹, grid,  ηⁿ⁺¹, ∫ᶻ_Axᶠᶜᶜ, ∫ᶻ_Ayᶜᶠᶜ, g, Δt,
                    dependencies = device_event(arch))

    wait(device(arch), event)

    fill_halo_regions!(L_ηⁿ⁺¹, arch)

    return nothing
end

# Kernels that act on vertically integrated / surface quantities
@inline ∫ᶻ_Ax_∂x_ηᶠᶜᶜ(i, j, k, grid, ∫ᶻ_Axᶠᶜᶜ, η) = @inbounds ∫ᶻ_Axᶠᶜᶜ[i, j, k] * ∂xᶠᶜᵃ(i, j, k, grid, η)
@inline ∫ᶻ_Ay_∂y_ηᶜᶠᶜ(i, j, k, grid, ∫ᶻ_Ayᶜᶠᶜ, η) = @inbounds ∫ᶻ_Ayᶜᶠᶜ[i, j, k] * ∂yᶜᶠᵃ(i, j, k, grid, η)

"""
Compute the horizontal divergence of vertically-uniform quantity using
vertically-integrated face areas `∫ᶻ_Axᶠᶜᶜ` and `∫ᶻ_Ayᶜᶠᶜ`.
"""
@inline Az_∇h²ᶜᶜᵃ(i, j, k, grid, ∫ᶻ_Axᶠᶜᶜ, ∫ᶻ_Ayᶜᶠᶜ, η::ZReducedField) =
    (δxᶜᵃᵃ(i, j, k, grid, ∫ᶻ_Ax_∂x_ηᶠᶜᶜ, ∫ᶻ_Axᶠᶜᶜ, η) +
     δyᵃᶜᵃ(i, j, k, grid, ∫ᶻ_Ay_∂y_ηᶜᶠᶜ, ∫ᶻ_Ayᶜᶠᶜ, η))

"""
    _implicit_free_surface_linear_operation!(L_ηⁿ⁺¹, grid, ηⁿ⁺¹, ∫ᶻ_Axᶠᶜᶜ, ∫ᶻ_Ayᶜᶠᶜ, g, Δt)

Return the left side of the "implicit η equation"

```math
(∇ʰ⋅H∇ʰ - 1/gΔt²) ηⁿ⁺¹ = 1 / (gΔt) ∇ʰ ⋅ Q★ - 1 / (gΔt²) ηⁿ
----------------------
        ≡ L_ηⁿ⁺¹
```

which is derived from the discretely summed barotropic mass conservation equation,
and arranged in a symmetric form by multiplying by horizontal areas Az:

```
δⁱÂʷ∂ˣηⁿ⁺¹ + δʲÂˢ∂ʸηⁿ⁺¹ - 1/gΔt² Az ηⁿ⁺¹ = 1 / (gΔt) (δⁱÂʷu̅ˢᵗᵃʳ + δʲÂˢv̅ˢᵗᵃʳ) - 1 / gΔt² Az ηⁿ
```

where  ̂ indicates a vertical integral, and
       ̅ indicates a vertical average
"""
@kernel function _implicit_free_surface_linear_operation!(L_ηⁿ⁺¹, grid, ηⁿ⁺¹, ∫ᶻ_Axᶠᶜᶜ, ∫ᶻ_Ayᶜᶠᶜ, g, Δt)
    i, j = @index(Global, NTuple)
    Az = Azᶜᶜᵃ(i, j, 1, grid)
    @inbounds L_ηⁿ⁺¹[i, j, 1] = Az_∇h²ᶜᶜᵃ(i, j, 1, grid, ∫ᶻ_Axᶠᶜᶜ, ∫ᶻ_Ayᶜᶠᶜ, ηⁿ⁺¹) - Az * ηⁿ⁺¹[i, j, 1] / (g * Δt^2)
end