text:
# --- Calling Embedded Julia from MATLAB --- #

# entry point for MATLAB calling Julia
jl_mex(plhs::Vector{Ptr{Cvoid}}, prhs::Vector{Ptr{Cvoid}}) = jl_mex_inner(plhs, prhs)

# runs julia function with mex function inputs, catching errors if they occur
function jl_mex_inner(plhs::Vector{Ptr{Cvoid}}, prhs::Vector{Ptr{Cvoid}})

    # get number of outputs
    nlhs = length(plhs)

    for i = 1:nlhs
        # default output is boolean false
        none = MATLAB.mxarray(false);
        # transfer ownership to MATLAB
        none.own = false
        # give pointer to MATLAB
        plhs[i] = none.ptr
    end

    try

        # extract function and arguments (function in first slot, arguments in remaining slots)
        fun = Core.eval(Main, Meta.parse(MATLAB.jstring(MATLAB.MxArray(prhs[1], false))))
        args = MATLAB.MxArray.(prhs[2:end], false)

        # call Julia function
        out = fun(args)

        out = (out isa Tuple ? out : (out,))::Tuple # wrap in tuple to simplify logic below

        # convert output(s) to MxArray type
        for i = 1:length(out)
            # stop early if max number of outputs is reached
            i > nlhs-1 && break
            # create MATLAB array for output
            mx = MATLAB.mxarray(out[i])
            # transfer ownership to MATLAB
            mx.own = false
            # give pointer to MATLAB
            plhs[1+i] = mx.ptr
        end

    catch exn

        # get backtrace
        bt = catch_backtrace()

        # create MATLABException from exception and backtrace
        mexn = MatlabException(exn, bt)

        # return MATLABexception in first output slot
        plhs[1] = mexn.ptr

    end

end

# --- Functions Used by MATLAB when calling embedded Julia --- #

# evaluates a Julia expressions
jl_eval(exprs::Vector{MATLAB.MxArray}) = ntuple(i -> Core.eval(Main, Meta.parse(MATLAB.jstring(exprs[i]))), length(exprs))

# Call a julia function, possibly with keyword arguments.
#
# The values in the args array are interpreted as follows:
#      Index = Meaning
# -----------------------------------------------------------
#          1 = the function to call
#          2 = an integer, npos, the number of positional arguments
#   3:2+npos = positional arguments
# 3+npos:end = keyword arguments, in keyword/value pairs
#
# If npos < 0, all arguments are assumed positional.
function jl_call_kw(args::Vector{MATLAB.MxArray})

    # process arguments
    nargs = length(args)

    # first argument is the function
    func = Meta.parse(MATLAB.jstring(args[1]))

    # second argument is the number of positional arguments
    npos = Int(MATLAB.jscalar(args[2]))::Int

    # construct the call expression
    expr = Expr(:call, func)

    # add positional arguments
    npos = npos < 0 ? nargs - 2 : npos # if npos is negative, all arguments are positional
    for i = 1:npos
        push!(expr.args, MATLAB.jvalue(args[2+i]))
    end

    # add keyword arguments
    nkw = div(nargs - 2 - npos, 2)
    for i = 1:nkw
        # assemble the key-value pair
        kw = Symbol(MATLAB.jstring(args[2+npos+(2*i-1)]))
        val = MATLAB.jvalue(args[2+npos+(2*i)])
        push!(expr.args, Expr(:kw, kw, val))
    end

    # return the evaluated expression
    return Core.eval(Main, expr)
end

# used for mimicing a basic Julia repl from the MATLAB console
input(prompt::String="julia> ") = call_matlab(1, "input", prompt, "s")[1]


text:
@testset "Periods" begin
    s = 1.0seconds
    m = 1.0minutes
    h = 1.0hours
    d = 1.0days
    y = 1.0years
    c = 1.0centuries
    @test s == AstroPeriod(seconds, 1.0)
    @test m == AstroPeriod(minutes, 1.0)
    @test h == AstroPeriod(hours, 1.0)
    @test d == AstroPeriod(days, 1.0)
    @test y == AstroPeriod(years, 1.0)
    @test c == AstroPeriod(centuries, 1.0)

    @test seconds(s) == 1.0seconds
    @test seconds(m) == 60.0seconds
    @test seconds(h) == 3600.0seconds
    @test seconds(d) == 86400.0seconds
    @test seconds(y) == 3.15576e7seconds
    @test seconds(c) == 3.15576e9seconds

    @test minutes(s) == (1.0 / 60.0)minutes
    @test minutes(m) == 1.0minutes
    @test minutes(h) == 60.0minutes
    @test minutes(d) == 1440.0minutes
    @test minutes(y) == 525960.0minutes
    @test minutes(c) == 5.2596e7minutes

    @test hours(s) == (1.0 / 3600.0)hours
    @test hours(m) == (1.0 / 60.0)hours
    @test hours(h) == 1.0hours
    @test hours(d) == 24.0hours
    @test hours(y) == 8766.0hours
    @test hours(c) == 876600.0hours

    @test days(s) == (1.0 / 86400.0)days
    @test days(m) == (1.0 / 1440.0)days
    @test days(h) == (1.0 / 24.0)days
    @test days(d) == 1.0days
    @test days(y) == 365.25days
    @test days(c) == 36525.0days

    @test years(s) == (1.0 / 3.15576e7)years
    @test years(m) == (1.0 / 525960.0)years
    @test years(h) == (1.0 / 8766.0)years
    @test years(d) == (1.0 / 365.25)years
    @test years(y) == 1.0years
    @test years(c) == 100.0years

    @test centuries(s) == (1.0 / 3.15576e9)centuries
    @test centuries(m) == (1.0 / 5.2596e7)centuries
    @test centuries(h) == (1.0 / 876600.0)centuries
    @test centuries(d) == (1.0 / 36525.0)centuries
    @test centuries(y) == (1.0 / 100.0)centuries
    @test centuries(c) == 1.0centuries

    @test zero(AstroPeriod{AstroTime.Periods.Year}) == 0.0years
    @test zero(1years) == 0years
    @test zero(1.0years) == 0.0years

    @test eltype(1.0years) == Float64
    @test eltype(typeof(1.0years)) == Float64

    @test 1.0seconds ≈ 1.0seconds
    @test 1.0 * seconds == seconds * 1.0
    @test 3 * 3.0seconds == 3.0seconds * 3
    @test 4.0seconds / 2 == 2.0seconds
    p = 1.0seconds
    @test -p == -1.0seconds
    @test 3.0seconds - 2.0seconds == 1.0seconds

    a = [1.0, 2.0, 3.0]
    @test a * seconds == [1.0seconds, 2.0seconds, 3.0seconds]
    @test a .* seconds == [1.0seconds, 2.0seconds, 3.0seconds]
    @test seconds * a == [1.0seconds, 2.0seconds, 3.0seconds]
    @test seconds .* a == [1.0seconds, 2.0seconds, 3.0seconds]
    @test Base.broadcastable(seconds) isa typeof(Ref(seconds))

    @test unit(1years) == years
    int_rng = 1seconds:3seconds
    @test step(int_rng) == 1seconds
    @test collect(int_rng) == [1seconds, 2seconds, 3seconds]
    float_rng = 1.0seconds:3.0seconds
    @test step(float_rng) == 1.0seconds
    @test collect(float_rng) == [1.0seconds, 2.0seconds, 3.0seconds]
    @test AstroPeriod{AstroTime.Periods.Second,Float64}(1.0seconds) == 1.0seconds

    @test Periods.name(seconds) == "seconds"
    @test Periods.name(minutes) == "minutes"
    @test Periods.name(hours) == "hours"
    @test Periods.name(days) == "days"
    @test Periods.name(years) == "years"
    @test Periods.name(centuries) == "centuries"

    @test string(1seconds) == "1.0 seconds"
    @test string(1minutes) == "1.0 minutes"
    @test string(1hours) == "1.0 hours"
    @test string(1days) == "1.0 days"
    @test string(1years) == "1.0 years"
    @test string(1centuries) == "1.0 centuries"
end

ext: jl

max_stars_repo_name: JuliaAstrodynamics/AstronomicalTime

max_issues_repo_name: JuliaAstrodynamics/AstronomicalTime

max_forks_repo_name: JuliaAstrodynamics/AstronomicalTime


text:
module DistributedFFT

using Distributed
using DistributedArrays
using DistributedArrays.SPMD
using FFTW

export brfft

function transpose3D(x)
    permutedims(x,(2,1,3))
end

@views function alltoalltranspose3D(x, nx, recievearray; pids=workers())
    pstart = minimum(pids)
    np = length(pids)
    rank = myid() - pstart + 1
    ny = size(x,2)
    extrax = nx % np
    extray = ny % np
    if rank <= extrax
        haveextra = true
    else
        haveextra = false
    end
    if rank <= extray
        haveextray = true
    else
        haveextray = false
    end

    if haveextra
        dx = size(x,1) - 1
    else
        dx = size(x,1)
    end
    nz = size(x,3)
    dy = convert(Int64,floor(ny/np))

    xt = transpose3D(x)

    i=1
    while i < np
        i *= 2
    end
    if i == np
        pow2 = true
    else
        pow2 = false
    end

    if haveextra
        if haveextray
            offset = rank - 1
            recievearray[:,(1+dx*(rank-1)+offset):dx*rank+(offset+1),:] .=
                        xt[(1+dy*(rank-1)+offset):dy*rank+(offset+1),:,:]
        else
            offsetx = rank - 1
            offsety = extray
            recievearray[:,(1+dx*(rank-1)+offsetx):dx*rank+(offsetx+1),:] .=
                        xt[(1+dy*(rank-1)+offsety):dy*rank+(offsety),:,:]
        end
    else
        if haveextray
            offsetx = extrax
            offsety = rank - 1
            recievearray[:,(1+dx*(rank-1)+offsetx):dx*rank+(offsetx),:] .=
                        xt[(1+dy*(rank-1)+offsety):dy*rank+(offsety+1),:,:]
        else
            offsetx = extrax
            offsety = extray
            recievearray[:,(1+dx*(rank-1)+offsetx):dx*rank+offsetx,:] .=
                        xt[(1+dy*(rank-1)+offsety):dy*rank+offsety,:,:]
        end
    end

    for i = 1:(np-1)
        if pow2
            source = destination = ((rank-1) ⊻ i) + 1
        else
            source = (rank - i - 1 + np ) % np + 1
            destination = (rank + i - 1 + np) % np + 1
        end
        if destination <= extray
            offset = destination - 1
            SPMD.sendto(destination + pstart - 1, xt[(1+dy*(destination-1)+offset):dy*destination+(offset+1),:,:])
        else
            offset = extray
            SPMD.sendto(destination + pstart - 1, xt[(1+dy*(destination-1)+offset):dy*destination+offset,:,:])
        end

        if source <= extrax
            offset = source - 1
            recievearray[:,(1+dx*(source-1)+offset):dx*source+(offset+1),:] .= SPMD.recvfrom(source + pstart - 1)
        else
            offset = extrax
            recievearray[:,(1+dx*(source-1)+offset):dx*source+offset,:] .= SPMD.recvfrom(source + pstart - 1)
        end
        barrier(; pids=pids)
    end
end

function FFTW.brfft(D::DArray, d::Int, dims=1:ndims(D))
    @assert div(d,2) + 1 == size(D,dims[1])
    @assert length(dims) == 3 #Temporary
    pids = procs(D)
    np = length(pids)
    Dout = dzeros((d, size(D,dims[2]), size(D,dims[3])), pids, [np,1,1])
    spmd(_mybrfft, D, d, Dout, dims; pids=procs(D))
    return Dout
end

function _mybrfft(D, d, Dout, dims)
    function mygenerate_plans(D, z, dims)
        plan1 = plan_bfft(D[:L], Tuple(collect(dims[2:end])))
        plan2 = plan_brfft(z, size(z,3), (2,))
        return [plan1, plan2]
    end

    function myapply_plans(D, Dout, dims, plans, z)
        pids = procs(D)[:,1]
        x = plans[1]*D[:L]
        alltoalltranspose3D(x, size(D,dims[1]), z; pids=pids)
        #mul!(Dout[:L], plans[2], z)
        w = plans[2]*z
        alltoalltranspose3D(w, size(D,2), Dout[:L]; pids=pids)
    end

    pids = procs(D)
    pstart = minimum(pids)
    np = length(pids)
    rank = myid() - pstart + 1

    bcnmesh = size(D,dims[1])
    bymesh = size(D,dims[2])
    bzmesh = size(D,dims[3])

    transposehaveextra = (rank <= rem(size(D,dims[2]), np))
    if transposehaveextra
        transpose_size = (div(bymesh, np) + 1, bcnmesh, bzmesh)
    else
        transpose_size = (div(bymesh, np), bcnmesh, bzmesh)
    end

    z = zeros(Complex{Float64}, transpose_size)

    plans = mygenerate_plans(D, z, dims)

    myapply_plans(D, Dout, dims, plans, z)
end
end


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

Oid(id::Oid) = id
Oid(ptr::Ptr{Oid}) = unsafe_load(ptr)::Oid

function Oid(ptr::Ptr{UInt8})
    if ptr == C_NULL
        throw(ArgumentError("NULL pointer passed to Oid() constructor"))
    end
    oid_ptr = Ref(Oid())
    ccall((:git_oid_fromraw, :libgit2), Void, (Ptr{Oid}, Ptr{UInt8}), oid_ptr, ptr)
    return oid_ptr[]
end

function Oid(id::Array{UInt8,1})
    if length(id) != OID_RAWSZ
        throw(ArgumentError("invalid raw buffer size"))
    end
    return Oid(pointer(id))
end

function Oid(id::AbstractString)
    bstr = String(id)
    len = sizeof(bstr)
    oid_ptr  = Ref(Oid())
    err = if len < OID_HEXSZ
        ccall((:git_oid_fromstrn, :libgit2), Cint,
              (Ptr{Oid}, Ptr{UInt8}, Csize_t), oid_ptr, bstr, len)
    else
        ccall((:git_oid_fromstrp, :libgit2), Cint,
              (Ptr{Oid}, Cstring), oid_ptr, bstr)
    end
    err != 0 && return Oid()
    return oid_ptr[]
end

function Oid(ref::GitReference)
    isempty(ref) && return Oid()
    reftype(ref) != Consts.REF_OID && return Oid()
    oid_ptr = ccall((:git_reference_target, :libgit2), Ptr{UInt8}, (Ptr{Void},), ref.ptr)
    oid_ptr == C_NULL && return Oid()
    return Oid(oid_ptr)
end

function Oid(repo::GitRepo, ref_name::AbstractString)
    isempty(repo) && return Oid()
    oid_ptr  = Ref(Oid())
    @check ccall((:git_reference_name_to_id, :libgit2), Cint,
                    (Ptr{Oid}, Ptr{Void}, Cstring),
                     oid_ptr, repo.ptr, ref_name)
    return oid_ptr[]
end

function Oid(obj::Ptr{Void})
    oid_ptr = ccall((:git_object_id, :libgit2), Ptr{UInt8}, (Ptr{Void},), obj)
    oid_ptr == C_NULL && return Oid()
    return Oid(oid_ptr)
end

function Oid{T<:GitObject}(obj::T)
    obj === nothing && return Oid()
    return Oid(obj.ptr)
end

Base.hex(id::Oid) = join([hex(i,2) for i in id.val])

raw(id::Oid) = collect(id.val)

Base.string(id::Oid) = hex(id)

Base.show(io::IO, id::Oid) = print(io, "Oid($(string(id)))")

Base.hash(id::Oid, h::UInt) = hash(id.val, h)

cmp(id1::Oid, id2::Oid) = Int(ccall((:git_oid_cmp, :libgit2), Cint,
                                    (Ptr{Oid}, Ptr{Oid}), Ref(id1), Ref(id2)))

==(id1::Oid, id2::Oid) = cmp(id1, id2) == 0
Base.isless(id1::Oid, id2::Oid)  = cmp(id1, id2) < 0

function iszero(id::Oid)
    for i in 1:OID_RAWSZ
        id.val[i] != zero(UInt8) && return false
    end
    return true
end

Base.zero(::Type{Oid}) = Oid()

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia


text:
using SFML
#=------------------------------------------------------------------------------
                            Constants
------------------------------------------------------------------------------=#
t = 8*pi/4
a = 0*pi/8
td = 0.1 * rand()
ad = 0.2 * rand()
g = 9.8
δ = 1/90 # delta
#=------------------------------------------------------------------------------
------------------------------------------------------------------------------=#
f₁ = 100
f₂ = 200
Δt = 0.01
sim_time = 10
tₖ = 0:Δt:sim_time
#=------------------------------------------------------------------------------
                            window setup
------------------------------------------------------------------------------=#
# Setup the window with some anti aliasing
settings = ContextSettings()
settings.antialiasing_level = 3
window = RenderWindow(VideoMode(800, 600), "Rosette test", settings, window_defaultstyle)

# Framerate limit
set_framerate_limit(window, 60)
event = Event()

# Create one circle
#=------------------------------------------------------------------------------
                        create the shapes
------------------------------------------------------------------------------=#
circle = CircleShape()
set_radius(circle, 10)
set_fillcolor(circle, SFML.red)
set_origin(circle, Vector2f(10, 10))

#=------------------------------------------------------------------------------
                        principal loop
------------------------------------------------------------------------------=#
while isopen(window)
    # Check for events
    while pollevent(window, event)
        if(get_type(event) == EventType.CLOSED)
            close(window)
        end
        # Allow the user to control the δ with the left and right arrow keys
        if(get_type(event) == EventType.KEY_PRESSED)
            if(get_key(event).key_code == KeyCode.LEFT)
                δ = -.1
            elseif(get_key(event).key_code == KeyCode.RIGHT)
                δ = .1
            end
        end
    end
    cycles = 5000
#=------------------------------------------------------------------------------
                        Simulate the Rosette
                    x(t) = (cos(f₁*t) + cos(f₂*t))
                    y(t) = (sin(f₁*t) - sin(f₂*t))
------------------------------------------------------------------------------=#
    tdd = 1
    add = 1

    for i = 1:cycles
        td += δ * tdd / cycles
        ad += δ * add / cycles
        t  += δ * td / cycles
        a  += δ * ad / cycles
    end
    x = 400 + 100 * (cos(t) + cos(a + t))
    y = 300 + 100 * (sin(t) - sin(a + t))
    # draw the circle
    set_position(circle, Vector2f(x, y))

    # Draw everything
    clear(window, SFML.black)
    draw(window, circle)
    display(window)
end


text:
using LinearAlgebra
using NearestNeighbors
n=200
K=10
μs = randn(n);
# generate true means
zs = sqrt(K) .* randn(n, K) .+ μs
# Observe 10 noisy observations for each mean
Zs = ReplicatedSample.(zs)

knn_cv_fit = StatsBase.fit(AuroraKNN(kKNN=100, loocv=true), Zs)

for j in Base.OneTo(K)
    knn_tmp_fit = StatsBase.fit(AuroraKNN(kKNN=knn_cv_fit.kskNN[j], loocv=false), Zs)
    @test knn_tmp_fit.μs_mat[:,j] ≈ knn_cv_fit.μs_mat[:,j]
end

# test Tree options

knn_fit_auto_loocv = StatsBase.fit(AuroraKNN(kKNN=22, loocv=true), Zs)
knn_fit_kd_loocv = StatsBase.fit(AuroraKNN(kKNN=22, loocv=true, tree=:kdtree), Zs)
knn_fit_ball_loocv = StatsBase.fit(AuroraKNN(kKNN=22, loocv=true, tree=:balltree), Zs)
knn_fit_brute_loocv = StatsBase.fit(AuroraKNN(kKNN=22, loocv=true, tree=:brutetree), Zs)

@test predict(knn_fit_auto_loocv) == predict(knn_fit_kd_loocv)
@test predict(knn_fit_kd_loocv) == predict(knn_fit_ball_loocv)
@test predict(knn_fit_kd_loocv) == predict(knn_fit_brute_loocv)

knn_fit_kd = StatsBase.fit(AuroraKNN(kKNN=20, loocv=false, tree=:kdtree), Zs)
knn_fit_ball = StatsBase.fit(AuroraKNN(kKNN=20, loocv=false, tree=:balltree), Zs)
knn_fit_brute = StatsBase.fit(AuroraKNN(kKNN=20, loocv=false, tree=:brutetree), Zs)

@test knn_fit_kd.kskNN == fill(20, 10)
@test predict(knn_fit_kd) ≈ predict(knn_fit_ball)
@test predict(knn_fit_kd) ≈ predict(knn_fit_brute)

for _ in Base.OneTo(10)
    i = sample(1:n)
    μs_tmp = zeros(K)
    for j in Base.OneTo(K)
        Zs_mat = Aurora.design_matrix(Zs, j)
        Y_vec = StatsBase.response.(Zs, j)
        dists_to_i = [norm(Zs_mat[:,i] .- zrow) for zrow in eachcol(Zs_mat)]
        sorted_idxs = sortperm(dists_to_i)
        @test sorted_idxs[1] == i
        μs_tmp[j] = mean(Y_vec[sorted_idxs[1:20]])
        @test  μs_tmp[j] ≈ knn_fit_kd.μs_mat[i, j]
    end
    @test mean(μs_tmp) ≈ knn_fit_kd.μs[i]
end

# test whether k=1 agrees with just sample mean

n=501
K=13
μs = sample(-1:1, n, replace=true)
# generate true means
zs =   sqrt(K) .* randn(n, K) .+ μs
# Observe 10 noisy observations for each mean
Zs = ReplicatedSample.(zs)

knn_cv_fit = StatsBase.fit(AuroraKNN(kKNN=1, loocv=true), Zs)
@test predict(knn_cv_fit) ≈ mean.(Zs)

knn_fit = StatsBase.fit(AuroraKNN(kKNN=1, loocv=false), Zs)
@test predict(knn_fit) ≈ mean.(Zs)

# Check the LOO NN code

kmax = 20
knn_cv_fit = StatsBase.fit(AuroraKNN(kKNN=kmax, loocv=true), Zs)
ks_opt = zeros(Int64, K)
mus_mat_opt = zeros(n, K)

for j in Base.OneTo(K)
    Zs_mat = Aurora.design_matrix(Zs, j)
    Y_vec = StatsBase.response.(Zs, j)

    kdd = KDTree(Zs_mat)

    loos = zeros(kmax)
    loos[1] = mean(abs2, Y_vec)

    for k in 1:(kmax-1)
        _idxs,_dists = knn(kdd, Zs_mat, k+1, true);
        mus_guess = [mean(Y_vec[_idx[2:end]]) for _idx in _idxs]
        loos[k+1] = mean(abs2, mus_guess .- Y_vec)
    end

    ks_opt[j] = argmin(loos)

    _idxs,_dists = knn(kdd, Zs_mat, ks_opt[j], true);
    mus_mat_opt[:,j] = [mean(Y_vec[_idx]) for _idx in _idxs]
end

@test ks_opt == knn_cv_fit.kskNN
@test mus_mat_opt ≈ knn_cv_fit.μs_mat


text:
## Exercise 8-1
## Write a function that takes a string as an argument and displays the letters backward, one per line.
println("Ans: ")

function reverse(s::AbstractString)
    s = s[end: -1: 1]

    for i in s
        println(i)
    end
end

reverse("Hello world")

println("End.")


text:
## Exercise 7-1
## Rewrite the function printn from Recursion using iteration instead of recursion.
println("Ans: ")

function printn(s, n::Int)
    while n > 0
        println(s)
        n -= 1
    end
end

printn("Hello iteration!", 4)

println("End.")


text:
# functions for running performance tests

using ODLCommonTools
using PDESolver
using EulerEquationMod
using AdvectionEquationMod
using Utils

"""
  Function that runs a test case and updates the performance history file.
  The case gets run twice, once for warmup, once for the real timing.

  The tests actually get run in a newly created directory called `tmpdir`.
  This directory is deleted if it already exists.  This directory should
  *not* be added to the git repo!

  **Inputs**

   * fname: the input file name
   * outname: the name of the performance history file
"""
function runtest(fname, outname)

  if isdir("./tmpdir")
    rm("./tmpdir", recursive=true)
  end

  mkdir("./tmpdir")
  cd("./tmpdir")
  fname = joinpath("..", fname)
  outname = joinpath("..", outname)

  # run the test in a try block so that even if something goes wrong, we
  # still cd back to the original directory (so we can run tests in a loop)
  try
    println("----- warming up -----")
    val, telapsed, tgc, gc_bytes = @time_all run_solver(fname)

    println("----- Final run -----")
    val, telapsed, tgc, gc_bytes = @time_all run_solver(fname)

    f = open(outname, "a+")
    tstr = getTimeString()
    branchname = getBranchName()
    println(f, tstr, " ", branchname, " ", telapsed, " ", gc_bytes)
    close(f)
  finally
    cd("..")
  end

  return nothing
end

"""
  Function to verify that none of the tests have the same input or output file
  name

  **Inputs**

   * vals: a n x 2 array of Strings, where each line contains the input file
           name followed by the output file name
"""
function checkInput(vals)

  println("size(vals) = ", size(vals))
  for i=1:size(vals, 1)
    println("i = ", i)
    for j=(i+1):size(vals, 1)
      println("j = ", j)
      @assert vals[j, 1] != vals[i, 1]
      @assert vals[j, 2] != vals[i, 2]

    end

    # check that files exist
    for d=1:2
      if !isfile(vals[i, d])
        error( "file $(vals[i, d]) does not exist")
      end
    end
  end

  return nothing
end

#runtest()


text:
using Test
using JuMP
using SemialgebraicSets
using PolyJuMP
using DynamicPolynomials

function plus_minus_test(optimizer,
                         config::MOI.Test.Config;
                         polymodule = nothing)
    atol = config.atol
    rtol = config.rtol

    model = _model(optimizer)

    if polymodule !== nothing
        setpolymodule!(model, polymodule)
    end

    @variable(model, α)
    @variable(model, β ≤ 1)

    @polyvar x y
    # x^2 ≤ 0 implies that x = 0
    #  VectorAffineFunction
    c1 = @constraint(model, α * (x + y) - β * y == 0, domain = @set x^2 ≤ 0)
    #  VectorOfVariables
    c2 = @constraint(model, α * x in PolyJuMP.ZeroPoly(),
                     domain = @set x^2 ≤ 0)

    @objective(model, Max, α)
    optimize!(model)

    @test termination_status(model) == MOI.OPTIMAL
    @test objective_value(model) ≈ 1.0 atol=atol rtol=rtol

    @test primal_status(model) == MOI.FEASIBLE_POINT
    @test value(α) ≈ 1.0 atol=atol rtol=rtol
    @test value(β) ≈ 1.0 atol=atol rtol=rtol
    @test value(UpperBoundRef(β)) ≈ 1.0 atol=atol rtol=rtol

    @test dual_status(model) == MOI.FEASIBLE_POINT
    @test dual(UpperBoundRef(β)) ≈ -1.0 atol=atol rtol=rtol

    _NonNegType(::ConstraintRef{<:JuMP.AbstractModel, MOI.ConstraintIndex{MOI.VectorOfVariables, PolyJuMP.PlusMinusSet{S}}}) where S = S

    @testset "Delete" begin
        F = MOI.VectorOfVariables
        G = MOI.VectorAffineFunction{Float64}
        NonNeg = _NonNegType(c2)
        S = PolyJuMP.PlusMinusSet{NonNeg}
        @test Set(MOI.get(model, MOI.ListOfConstraintTypesPresent())) == Set([
            (MOI.VariableIndex, MOI.LessThan{Float64}), (G, S), (F, S)])
        test_delete_bridge(model, c2, 2, ((F, NonNeg, 0), (G, NonNeg, 0)))
        test_delete_bridge(model, c1, 2, ((F, NonNeg, 0), (G, NonNeg, 0)))
    end
end

linear_tests["plus_minus"] = plus_minus_test


text:
using AMDB
using Base.Test
using IntervalSets

≅(a,b) = isequal(a,b)

# parsing

@test parsedtype(AMDB_Date) ≅ Date
@test tryparse(AMDB_Date, "00000000") ≅ Nullable{Date}()
@test tryparse(AMDB_Date, "19800101") ≅ Nullable(Date(1980, 1, 1))
@test tryparse(AMDB_Date, "19800100") ≅ Nullable(Date(1980, 1, 1))
@test tryparse(AMDB_Date, "") ≅ Nullable{Date}()
@test tryparse(AMDB_Date, "nonsense") ≅ Nullable{Date}()

@test parsedtype(Gender) ≅ Gender
@test tryparse(Gender, "M") ≅ Nullable(male)
@test tryparse(Gender, "F") ≅ Nullable(female)
@test tryparse(Gender, "") ≅ Nullable{Gender}()
@test tryparse(Gender, "nonsense") ≅ Nullable{Gender}()

@test parsedtype(MaybeEmpty{Int}) ≅ Nullable{Int}
@test tryparse(MaybeEmpty{Int}, "") ≅ Nullable(Nullable{Int}())
@test tryparse(MaybeEmpty{Int}, "42") ≅ Nullable(Nullable(42))
@test tryparse(MaybeEmpty{Int}, "42.2") ≅ Nullable{Nullable{Int}}()
@test tryparse(MaybeEmpty{Int}, "nonsense") ≅ Nullable{Nullable{Int}}()

@test parsedtype(MaybeEmpty{AMDB_Date}) ≅ Nullable{Date}
@test tryparse(MaybeEmpty{AMDB_Date}, "") ≅ Nullable(Nullable{Date}())
@test tryparse(MaybeEmpty{AMDB_Date}, "19700907") ≅
    Nullable{Nullable{Date}}(Date(1970,9,7))
@test tryparse(MaybeEmpty{AMDB_Date}, "42.2") ≅ Nullable{Nullable{Date}}()
@test tryparse(MaybeEmpty{AMDB_Date}, "nonsense") ≅ Nullable{Nullable{Date}}()

let ae = AutoEnum{Char}()
    v = map(x->get!(ae, x), ['a', 'b', 'a', 'a', 'c', 'b'])
    @test v == [1, 2, 1, 1, 3, 2]
    @test keys(ae) == ['a', 'b', 'c']
end

date_interval = ClosedInterval(Date(1980, 1, 1), Date(1980, 3, 1))
@test tryparse(PersonSpellParser, "1900;M;0;;0;;19800101;19800301;20010101;") ≅
    Nullable(PersonSpell(1900, male, 0, Nullable{Date}(), date_interval))
@test tryparse(AMPSpellParser, "1400;RE;19800101;19800301;T;") ≅
    Nullable(AMPSpell(1400, AMP.retired, date_interval))
@test tryparse(VMZSpellParser, "19800101;19800301;AL;A;20010101;") ≅
    Nullable(VMZSpell(VMZ.unemployed, VMZ.employment_Austria, date_interval))


text:
# This file is part of JuliaFMI.
# Licensed under MIT: https://github.com/AnHeuermann/JuliaFMI/blob/master/LICENSE.txt

# This file contains simple tests to validate JuliaFMI simulations.

"""
Define tests for unit testing
"""

thisDir = dirname(Base.source_path())
include("$(dirname(thisDir))/src/FMUSimulator.jl")

if Sys.iswindows()
    fmuTestDir = joinpath(thisDir, "winFMU")
elseif Sys.islinux()
    fmuTestDir = joinpath(thisDir, "linuxFMU")
elseif Sys.isapple()
    fmuTestDir = joinpath(thisDir, "darwinFMU")
else
    error("OS not supported for this tests.")
end

"""
    testHelloFMI20World()

Run most basic simulation test.

FMU generated for Modelica model:
```
    model HelloFMI20World
      Real x(start=1);
      parameter Real a=2;
    equation
      der(x) = a * x;
    end HelloFMI20World;
```
"""
function testHelloFMI20World()
    helloFMI20World = joinpath(fmuTestDir,"HelloFMI20World.fmu")
    main(helloFMI20World)
    return true
end

"""
    testBouncingBall()

Run simple test simulating a bouncing ball.

FMU generated from Modelica model:
```
    model BouncingBallFMI20
      parameter Real e=0.7 \"coefficient of restitution\";
      parameter Real g=9.81 \"gravity acceleration\";
      Real h(start=1) \"height of ball\";
      Real v \"velocity of ball\";
      Real v_new;
    equation
      der(v) = -g;
      der(h) = v;

      when h <= 0.0 then
        v_new = -e*pre(v);
        reinit(v, v_new);
      end when;
    end BouncingBallFMI20;
```
"""
function testBouncingBall()
    bouncingBall = joinpath(fmuTestDir,"BouncingBallFMI20.fmu")
    main(bouncingBall)
    return true
end

function testCauerLowPassAnalog()
    cauerLowPass = joinpath(fmuTestDir,"Modelica_Electrical_Analog_Examples_CauerLowPassAnalog.fmu")
    main(cauerLowPass)
    return true
end

max_stars_repo_name: JKRT/JuliaFMI

max_issues_repo_name: JKRT/JuliaFMI

max_forks_repo_name: JKRT/JuliaFMI


text:
# Stella, Themelis, Sopasakis, Patrinos, "A simple and efficient algorithm
# for nonlinear model predictive control", 56th IEEE Conference on Decision
# and Control (2017).

using Base.Iterators
using ProximalAlgorithms: LBFGS
using ProximalAlgorithms.IterationTools
using ProximalOperators: Zero
using LinearAlgebra
using Printf

struct PANOC_iterable{R <: Real, C <: Union{R, Complex{R}}, Tx <: AbstractArray{C}, Tf, TA, Tg}
    f::Tf             # smooth term
    A::TA             # matrix/linear operator
    g::Tg             # (possibly) nonsmooth, proximable term
    x0::Tx            # initial point
    alpha::R          # in (0, 1), e.g.: 0.95
    beta::R           # in (0, 1), e.g.: 0.5
    gamma::Maybe{R}   # stepsize parameter of forward and backward steps
    adaptive::Bool    # enforce adaptive stepsize even if L is provided
    memory::Int       # memory parameter for L-BFGS
end

mutable struct PANOC_state{R <: Real, Tx, TAx}
    x::Tx             # iterate
    Ax::TAx           # A times x
    f_Ax::R           # value of smooth term
    grad_f_Ax::TAx    # gradient of f at Ax
    At_grad_f_Ax::Tx  # gradient of smooth term
    gamma::R          # stepsize parameter of forward and backward steps
    y::Tx             # forward point
    z::Tx             # forward-backward point
    g_z::R            # value of nonsmooth term (at z)
    res::Tx           # fixed-point residual at iterate (= x - z)
    H::LBFGS.LBFGS_buffer{R} # variable metric
    tau::Maybe{R}     # stepsize (can be nothing since the initial state doesn't have it)
    # some additional storage:
    d::Tx
    Ad::TAx
    x_d::Tx
    Ax_d::TAx
    f_Ax_d::R
    grad_f_Ax_d::TAx
    At_grad_f_Ax_d::Tx
    z_curr::Tx
end

PANOC_state(
    x::Tx, Ax::TAx, f_Ax::R, grad_f_Ax, At_grad_f_Ax, gamma::R, y, z, g_z, res, H, tau
) where {R, Tx, TAx} =
    PANOC_state{R, Tx, TAx}(
        x, Ax, f_Ax, grad_f_Ax, At_grad_f_Ax, gamma, y, z, g_z, res, H, tau,
        zero(x), zero(Ax), zero(x), zero(Ax), zero(R), zero(Ax), zero(x), zero(x)
    )

f_model(state::PANOC_state) = f_model(state.f_Ax, state.At_grad_f_Ax, state.res, state.gamma)

function Base.iterate(iter::PANOC_iterable{R}) where R
    x = iter.x0
    Ax = iter.A * x
    grad_f_Ax, f_Ax = gradient(iter.f, Ax)

    gamma = iter.gamma

    if gamma === nothing
        # compute lower bound to Lipschitz constant of the gradient of x ↦ f(Ax)
        xeps = x .+ R(1)
        grad_f_Axeps, f_Axeps = gradient(iter.f, iter.A*xeps)
        L = norm(iter.A' * (grad_f_Axeps - grad_f_Ax)) / R(sqrt(length(x)))
        gamma = iter.alpha/L
    end

    # compute initial forward-backward step
    At_grad_f_Ax = iter.A' * grad_f_Ax
    y = x - gamma .* At_grad_f_Ax
    z, g_z = prox(iter.g, y, gamma)

    # compute initial fixed-point residual
    res = x - z

    # initialize variable metric
    H = LBFGS.create(x, iter.memory)

    state = PANOC_state(x, Ax, f_Ax, grad_f_Ax, At_grad_f_Ax, gamma, y, z, g_z, res, H, nothing)

    return state, state
end

function Base.iterate(iter::PANOC_iterable{R}, state::PANOC_state{R, Tx, TAx}) where {R, Tx, TAx}
    Az, f_Az, grad_f_Az, At_grad_f_Az = nothing, nothing, nothing, nothing
    a, b, c = nothing, nothing, nothing

    f_Az_upp = f_model(state)

    # backtrack gamma (warn and halt if gamma gets too small)
    while iter.gamma === nothing || iter.adaptive == true
        if state.gamma < 1e-7 # TODO: make this a parameter, or dependent on R?
            @warn "parameter `gamma` became too small ($(state.gamma)), stopping the iterations"
            return nothing
        end
        Az = iter.A*state.z
        grad_f_Az, f_Az = gradient(iter.f, Az)
        tol = 10*eps(R)*(1 + abs(f_Az))
        if f_Az <= f_Az_upp + tol break end
        state.gamma *= 0.5
        state.y .= state.x .- state.gamma .* state.At_grad_f_Ax
        state.g_z = prox!(state.z, iter.g, state.y, state.gamma)
        state.res .= state.x .- state.z
        LBFGS.reset!(state.H)
        f_Az_upp = f_model(state)
    end

    # compute FBE
    FBE_x = f_Az_upp + state.g_z

    # update metric
    LBFGS.update!(state.H, state.x, state.res)

    # compute direction
    mul!(state.d, state.H, -state.res)

    # backtrack tau 1 → 0
    tau = R(1)
    mul!(state.Ad, iter.A, state.d)

    state.x_d .= state.x .+ state.d
    state.Ax_d .= state.Ax .+ state.Ad
    state.f_Ax_d = gradient!(state.grad_f_Ax_d, iter.f, state.Ax_d)
    mul!(state.At_grad_f_Ax_d, adjoint(iter.A), state.grad_f_Ax_d)

    state.x .= state.x_d
    state.Ax .= state.Ax_d
    state.f_Ax = state.f_Ax_d
    state.grad_f_Ax .= state.grad_f_Ax_d
    state.At_grad_f_Ax .= state.At_grad_f_Ax_d

    copyto!(state.z_curr, state.z)

    sigma = iter.beta * (0.5/state.gamma) * (1 - iter.alpha)
    tol = 10*eps(R)*(1 + abs(FBE_x))
    threshold = FBE_x - sigma * norm(state.res)^2 + tol

    for i = 1:20
        state.y .= state.x .- state.gamma .* state.At_grad_f_Ax
        state.g_z = prox!(state.z, iter.g, state.y, state.gamma)
        state.res .= state.x .- state.z
        FBE_x_new = f_model(state) + state.g_z

        if FBE_x_new <= threshold
            state.tau = tau
            return state, state
        end

        if Az === nothing Az = iter.A * state.z_curr end

        tau *= 0.5
        state.x .= tau .* state.x_d .+ (1-tau) .* state.z_curr
        state.Ax .= tau .* state.Ax_d .+ (1-tau) .* Az

        if ProximalOperators.is_quadratic(iter.f)
            # in case f is quadratic, we can compute its value and gradient
            # along a line using interpolation and linear combinations
            # this allows saving operations
            if grad_f_Az === nothing grad_f_Az, f_Az = gradient(iter.f, Az) end
            if At_grad_f_Az === nothing
                At_grad_f_Az = iter.A' * grad_f_Az
                c = f_Az
                b = real(dot(state.Ax_d .- Az, grad_f_Az))
                a = state.f_Ax_d - b - c
            end
            state.f_Ax = a * tau^2 + b * tau + c
            state.grad_f_Ax .= tau .* state.grad_f_Ax_d .+ (1-tau) .* grad_f_Az
            state.At_grad_f_Ax .= tau .* state.At_grad_f_Ax_d .+ (1-tau) .* At_grad_f_Az
        else
            # otherwise, in the general case where f is only smooth, we compute
            # one gradient and matvec per backtracking step
            state.f_Ax = gradient!(state.grad_f_Ax, iter.f, state.Ax)
            mul!(state.At_grad_f_Ax, adjoint(iter.A), state.grad_f_Ax)
        end
    end

    @warn "stepsize `tau` became too small ($(tau)), stopping the iterations"
    return nothing
end

# Solver

struct PANOC{R <: Real}
    alpha::R
    beta::R
    gamma::Maybe{R}
    adaptive::Bool
    memory::Int
    maxit::Int
    tol::R
    verbose::Bool
    freq::Int

    function PANOC{R}(; alpha::R=R(0.95), beta::R=R(0.5),
        gamma::Maybe{R}=nothing, adaptive::Bool=false, memory::Int=5,
        maxit::Int=1000, tol::R=R(1e-8), verbose::Bool=false, freq::Int=10
    ) where R
        @assert 0 < alpha < 1
        @assert 0 < beta < 1
        @assert gamma === nothing || gamma > 0
        @assert memory >= 0
        @assert maxit > 0
        @assert tol > 0
        @assert freq > 0
        new(alpha, beta, gamma, adaptive, memory, maxit, tol, verbose, freq)
    end
end

function (solver::PANOC{R})(
    x0::AbstractArray{C}; f=Zero(), A=I, g=Zero(), L::Maybe{R}=nothing
) where {R, C <: Union{R, Complex{R}}}

    stop(state::PANOC_state) = norm(state.res, Inf)/state.gamma <= solver.tol
    disp((it, state)) = @printf(
        "%5d | %.3e | %.3e | %.3e\n",
        it, state.gamma, norm(state.res, Inf)/state.gamma,
        (state.tau === nothing ? 0.0 : state.tau)
    )

    if solver.gamma === nothing && L !== nothing
        gamma = solver.alpha/L
    elseif solver.gamma !== nothing
        gamma = solver.gamma
    end

    iter = PANOC_iterable(
        f, A, g, x0,
        solver.alpha, solver.beta, solver.gamma, solver.adaptive, solver.memory
    )
    iter = take(halt(iter, stop), solver.maxit)
    iter = enumerate(iter)
    if solver.verbose iter = tee(sample(iter, solver.freq), disp) end

    num_iters, state_final = loop(iter)

    return state_final.z, num_iters

end

# Outer constructors

"""
    PANOC([gamma, adaptive, memory, maxit, tol, verbose, freq, alpha, beta])

Instantiate the PANOC algorithm (see [1]) for solving optimization problems
of the form

    minimize f(Ax) + g(x),

where `f` is smooth and `A` is a linear mapping (for example, a matrix).
If `solver = PANOC(args...)`, then the above problem is solved with

    solver(x0, [f, A, g, L])

Optional keyword arguments:

* `gamma::Real` (default: `nothing`), the stepsize to use; defaults to `alpha/L` if not set (but `L` is).
* `adaptive::Bool` (default: `false`), if true, forces the method stepsize to be adaptively adjusted even if `L` is provided (this behaviour is always enforced if `L` is not provided).
* `memory::Integer` (default: `5`), memory parameter for L-BFGS.
* `maxit::Integer` (default: `1000`), maximum number of iterations to perform.
* `tol::Real` (default: `1e-8`), absolute tolerance on the fixed-point residual.
* `verbose::Bool` (default: `true`), whether or not to print information during the iterations.
* `freq::Integer` (default: `10`), frequency of verbosity.
* `alpha::Real` (default: `0.95`), stepsize to inverse-Lipschitz-constant ratio; should be in (0, 1).
* `beta::Real` (default: `0.5`), sufficient decrease parameter; should be in (0, 1).

If `gamma` is not specified at construction time, the following keyword
argument can be used to set the stepsize parameter:

* `L::Real` (default: `nothing`), the Lipschitz constant of the gradient of x ↦ f(Ax).

References:

[1] Stella, Themelis, Sopasakis, Patrinos, "A simple and efficient algorithm
for nonlinear model predictive control", 56th IEEE Conference on Decision
and Control (2017).
"""
PANOC(::Type{R}; kwargs...) where R = PANOC{R}(; kwargs...)
PANOC(; kwargs...) = PANOC(Float64; kwargs...)


text:
"""
```
hair_plot(var, realized, histories, forecasts; kwargs...)

hair_plot(var, realized, initial_values, forecasts; plotroot = "",
    verbose = :low, kwargs...)
```

### Inputs

- `var::Symbol`: e.g. `:obs_gdp`
- `realized::DataFrame`: must contain realized values of `var`
- `histories::Vector{MeansBands}` (method 1) or
  `initial_values::Vector{Float64}` (method 2): vector of either historical
  `MeansBands` or initial forecast values (i.e. s_{T|T} or y_T). Needed to
  connect the forecast hairs to the realized data line
- `forecasts::Vector{MeansBands}`

### Keyword Arguments

- `plotroot::String`: if nonempty, plots will be saved in that directory
- `verbose::Symbol`

See `?hair` for additional keyword arguments, all of which can be passed into
`hair_plot`.

### Output

- `p::Plot`
"""
function hair_plot(var::Symbol, realized::DataFrame,
                   histories::Vector{MeansBands}, forecasts::Vector{MeansBands};
                   kwargs...)

    initial_values = map(history -> history.means[end, var], histories)
    hair_plot(var, realized, initial_values, forecasts; kwargs...)
end

function hair_plot(var::Symbol, realized::DataFrame,
                   initial_values::Vector{Float64}, forecasts::Vector{MeansBands};
                   plotroot::String = "",
                   verbose::Symbol = :low,
                   kwargs...)
    # Call recipe
    p = hair(var, realized, initial_values, forecasts; kwargs...)

    # Save plot
    if !isempty(plotroot)
        output_file = joinpath(plotroot, "hairplot_" * detexify(string(var)) * "." *
                               string(plot_extension()))
        save_plot(p, output_file, verbose = verbose)
    end

    return p
end

@userplot Hair

"""
```
hair(var, realized, initial_values, forecasts;
    hist_label = \"Realized\", forecast_label = \"Forecasts\",
    hist_color = :black, forecast_color = :red, forecast_palette = :none,
    tick_size = 2)
```

User recipe called by `hair_plot`.

### Inputs

- `var::Symbol`: e.g. `:obs_gdp`
- `initial_values::Vector{Float64}`: vector of initial forecast values (i.e. s_{T|T} or y_T). Needed to
  connect the forecast hairs to the realized data line
- `forecasts::Vector{MeansBands}`

### Keyword Arguments

- `hist_label::String`
- `forecast_label::String`
- `forecast_color`
- `forecast_palette`: if not `:none`, the hair colors will be chosen according
  to this palette; otherwise they will all be `forecast_color`. Values
  correspond to values of the Plots attribute `color_palette` (see
  docs.juliaplots.org/latest/attributes)
- `tick_size::Int`: x-axis (time) tick size in units of years

Additionally, all Plots attributes (see docs.juliaplots.org/latest/attributes)
are supported as keyword arguments.
"""
hair

@recipe function f(hp::Hair;
                   hist_label = "Realized",
                   forecast_label = "Forecasts",
                   hist_color = :black,
                   forecast_color = :red,
                   forecast_palette = :none,
                   tick_size = 2)
    # Error checking
    if length(hp.args) != 4 || typeof(hp.args[1]) != Symbol || typeof(hp.args[2]) != DataFrame ||
        !(typeof(hp.args[3]) <: AbstractVector) ||
        !(typeof(hp.args[4]) <: AbstractVector{MeansBands})

        error("hair must be given Tuple{Symbol, DataFrame, AbstractVector, AbstractVector{MeansBands}}. Got $(typeof(hf.args))")
    end

    var, realized, initial_values, forecasts = hp.args
    if length(initial_values) != length(forecasts)
        error("Lengths of initial_values ($length(initial_values)) and forecasts ($length(forecasts)) do not match")
    end

    # Assign date ticks
    date_ticks = Base.filter(x -> Dates.month(x) == 3,            realized[!, :date])
    date_ticks = Base.filter(x -> Dates.year(x) % tick_size == 0, date_ticks)
    xticks --> (date_ticks, map(Dates.year, date_ticks))

    # Realized series
    @series begin
        seriestype := :line
        linewidth := 2
        linecolor := hist_color
        label     := hist_label

        realized[!, :date], realized[!, var]
    end

    # Forecasts
    for (initial_value, forecast) in zip(initial_values, forecasts)
        @series begin
            seriestype := :line
            linewidth  := 1
            label      := forecast == forecasts[1] ? forecast_label : ""
            if forecast_palette == :none
                linecolor := forecast_color
            else
                palette   := forecast_palette
            end

            initial_date = iterate_quarters(forecast.means[1, :date], -1)
            x = vcat(initial_date,  forecast.means[!, :date])
            y = vcat(initial_value, forecast.means[!, var])
            x, y
        end
    end
end


text:
# Domain
x_max = 0.005

# Models
sigma = 0.0004
filter = (x, y) -> exp(-(x^2 + y^2)/2/sigma^2)
filter_dx = (x, y) -> -x/sigma^2*filter(x, y)
filter_dy = (x, y) -> -y/sigma^2*filter(x, y)
K = 2
v_max = 0.05
tau = 1/200

# Noise
sigma_noise = 0.001

# Medium
n_im = 500
dx = 0.0002

# Inverse problem
lp_norm = 1
jump_threshold = 0.2

# Sequence
p = 0.05
initial_position_generator =
    () -> x_max/2 + (1-2*rand()) * x_max/4 # so that it's not too close to the border


text:
@testset "Bus Constructors" begin
    tBus = Bus(nothing)
    tLoadZone = LoadZone(nothing)

    bus = Bus(
        1,
        "test",
        BusTypes.SLACK,
        0.0,
        0.0,
        (min = 0.0, max = 0.0),
        nothing,
        nothing,
        nothing,
    )
    @test PowerSystems.get_bustype(bus) == BusTypes.REF

    @test_throws(
        PowerSystems.DataFormatError,
        Bus(
            1,
            "test",
            BusTypes.ISOLATED,
            0.0,
            0.0,
            (min = 0.0, max = 0.0),
            nothing,
            nothing,
            nothing,
        )
    )
end

@testset "Generation Constructors" begin
    tThreePartCost = ThreePartCost(nothing)
    @test tThreePartCost isa PowerSystemType
    tTwoPartCost = TwoPartCost(nothing)
    @test tTwoPartCost isa PowerSystemType
    tThermalGen = ThermalStandard(nothing)
    @test tThermalGen isa PowerSystems.Component
    tHydroDispatch = HydroDispatch(nothing)
    @test tHydroDispatch isa PowerSystems.Component
    tHydroEnergyReservoir = HydroEnergyReservoir(nothing)
    @test tHydroEnergyReservoir isa PowerSystems.Component
    tRenewableFix = RenewableFix(nothing)
    @test tRenewableFix isa PowerSystems.Component
    tRenewableDispatch = RenewableDispatch(nothing)
    @test tRenewableDispatch isa PowerSystems.Component
    tRenewableDispatch = RenewableDispatch(nothing)
    @test tRenewableDispatch isa PowerSystems.Component
end

@testset "Source Constructors" begin
    tSource = Source(nothing)
    @test tSource isa PowerSystems.Component
end

@testset "Storage Constructors" begin
    tStorage = GenericBattery(nothing)
    @test tStorage isa PowerSystems.Component
end

@testset "Load Constructors" begin
    tPowerLoad = PowerLoad(nothing)
    @test tPowerLoad isa PowerSystems.Component
    tPowerLoadPF = PowerLoadPF(nothing)
    @test tPowerLoadPF isa PowerSystems.Component
    tPowerLoad = PowerLoad("init", true, Bus(nothing), nothing, 0.0, 0.0, 0.0, 0.0)
    @test tPowerLoad isa PowerSystems.Component
    tPowerLoadPF = PowerLoadPF("init", true, Bus(nothing), nothing, 0.0, 0.0, 1.0)
    @test tPowerLoadPF isa PowerSystems.Component
    tLoad = InterruptibleLoad(nothing)
    @test tLoad isa PowerSystems.Component
end

@testset "Branch Constructors" begin
    tLine = Line(nothing)
    @test tLine isa PowerSystems.Component
    tMonitoredLine = MonitoredLine(nothing)
    @test tMonitoredLine isa PowerSystems.Component
    tHVDCLine = HVDCLine(nothing)
    @test tHVDCLine isa PowerSystems.Component
    tVSCDCLine = VSCDCLine(nothing)
    @test tVSCDCLine isa PowerSystems.Component
    tTransformer2W = Transformer2W(nothing)
    @test tTransformer2W isa PowerSystems.Component
    tTapTransformer = TapTransformer(nothing)
    @test tTapTransformer isa PowerSystems.Component
    tPhaseShiftingTransformer = PhaseShiftingTransformer(nothing)
    @test tPhaseShiftingTransformer isa PowerSystems.Component
end

@testset "Service Constructors" begin
    tStaticReserve = StaticReserve{ReserveUp}(nothing)
    @test tStaticReserve isa PowerSystems.Service
    tVariableReserve = VariableReserve{ReserveDown}(nothing)
    @test tVariableReserve isa PowerSystems.Service
end

@testset "Forecast Constructors" begin
    tg = RenewableFix(nothing)
    forecast_data = PowerSystems.TimeSeries.TimeArray(
        [DateTime("01-01-01"), DateTime("01-01-01") + Hour(1)],
        [1.0, 1.0],
    )
    #Deterministic Tests
    tDeterministicForecast =
        PSY.Deterministic("scalingfactor", Hour(1), DateTime("01-01-01"), 24)
    @test tDeterministicForecast isa PowerSystems.Forecast
    tDeterministicForecast = PSY.Deterministic("scalingfactor", forecast_data)
    @test tDeterministicForecast isa PowerSystems.Forecast
    #Probabilistic Tests
    tProbabilisticForecast =
        PSY.Probabilistic("scalingfactor", Hour(1), DateTime("01-01-01"), [0.5, 0.5], 24)
    @test tProbabilisticForecast isa PowerSystems.Forecast
    tProbabilisticForecast = PSY.Probabilistic("scalingfactor", [1.0], forecast_data)
    @test tProbabilisticForecast isa PowerSystems.Forecast
    #Scenario Tests
    tScenarioForecast =
        PSY.ScenarioBased("scalingfactor", Hour(1), DateTime("01-01-01"), 2, 24)
    @test tScenarioForecast isa PowerSystems.Forecast
    tScenarioForecast = PSY.ScenarioBased("scalingfactor", forecast_data)
    @test tScenarioForecast isa PowerSystems.Forecast
end

@testset "Regulation Device" begin
    original_device = ThermalStandard(nothing)
    regulation = RegulationDevice(original_device)
    @test get_rating(regulation) == 0.0
    set_rating!(regulation, 10.0)
    @test get_rating(regulation) == 10.0
    regulation = RegulationDevice(original_device, droop = 0.5)
    @test get_droop(regulation) == 0.5
    @test get_participation_factor(regulation) == (up = 0.0, dn = 0.0)
    @test get_reserve_limit_up(regulation) == 0.0
    @test get_reserve_limit_dn(regulation) == 0.0
    @test get_inertia(regulation) == 0.0
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    Quadrangle(p1, p2, p3, p4)

A quadrangle with points `p1`, `p2`, `p3`, `p4`.
"""
struct Quadrangle{Dim,T,V<:AbstractVector{Point{Dim,T}}} <: Polygon{Dim,T}
  vertices::V
end

function measure(q::Quadrangle)
  vs = q.vertices
  Δ₁ = Triangle(view(vs, [1,2,3]))
  Δ₂ = Triangle(view(vs, [3,4,1]))
  measure(Δ₁) + measure(Δ₂)
end

function Base.in(p::Point, q::Quadrangle)
  vs = q.vertices
  Δ₁ = Triangle(view(vs, [1,2,3]))
  Δ₂ = Triangle(view(vs, [3,4,1]))
  p ∈ Δ₁ || p ∈ Δ₂
end


text:
# Autogenerated wrapper script for Unicorn_jll for x86_64-linux-musl
export libunicorn

JLLWrappers.@generate_wrapper_header("Unicorn")
JLLWrappers.@declare_library_product(libunicorn, "libunicorn.so.1")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libunicorn,
        "lib/libunicorn.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
function fill_with_nothing!(a)  # TODO: find more elegant way to do this
    for k = keys(a)
        a[k] = nothing
    end
end

function Radau_for_MechanismScenario(m::MechanismScenario{Dual{Type_Tag,Float64,NC}}) where {Type_Tag,NC}
    NX = num_x(m)
    return makeRadauIntegrator(m, NX, 1.0e-16, 2, NC)
end

as_static_vector(f::Wrench{T}) where {T} = vcat(angular(f), linear(f))
as_static_vector(f::Twist{T})  where {T} = vcat(angular(f), linear(f))

function get_bristle_d0(tm::TypedMechanismScenario{T}, bristle_id::BristleID) where {T}
    s = segments(tm.s)[bristle_id]
    return SVector{6,T}(s[1], s[2], s[3], s[4], s[5], s[6])
end
@inline get_bristle_d1(tm::TypedMechanismScenario{T}, bristle_id::BristleID) where {T} = segments(tm.ṡ)[bristle_id]

function find_mesh_id(ts::MeshCacheDict{MeshCache}, name::String)  # TODO: make this function more elegant
    id = MeshID(-9999)
    for k = keys(ts)
        if ts[k].name == name
            (id == MeshID(-9999)) || error("multiple")
            id = k
        end
    end
    (id == MeshID(-9999)) && error("no mesh found by name: $name")
    return id
end

function find_mesh_id(ts::MeshCacheDict{MeshCache}, mc::MeshCache)  # TODO: make this function more elegant
    id = MeshID(-9999)
    for k = keys(ts)
        if ts[k] == mc
            (id == MeshID(-9999)) || error("multiple")
            id = k
        end
    end
    (id == MeshID(-9999)) && error("no mesh found by name: $(mc.name)")
    return id
end

find_mesh_id(ts::MechanismScenario, input_2) = find_mesh_id(ts.MeshCache, input_2)

find_mesh(ts::MeshCacheDict{MeshCache}, name::String) = ts[find_mesh_id(ts, name)]
find_mesh(m::MechanismScenario, name::String) = find_mesh(m.MeshCache, name)


text:
@testset "TabularQApproximator" begin
    A = TabularQApproximator([1 2; 3 4])

    @test A(1) == [1.0, 2.0]
    @test A(2) == [3.0, 4.0]
    @test A(2, 2) == 4

    update!(A, (2, 2) => -1.0)
    @test A(2) == [3.0, 3.0]
    @test A(2, 2) == 3

    update!(A, 1 => [-1.0, -1.0])
    @test A(1) == [0.0, 1.0]
    @test A(1, 1) == 0
    @test A(1, 2) == 1
end


text:
# Copyright (c) 2015: Christof Stocker

"""
`Pegasos <: SvmSolver`

Description
============

Implementation of the primal sub-gradient solver for SVM
as it was proposed in (Shalev-Shwartz et al., 2007).

Usage
======

    model = svm(X, y, solver = Pegasos())

References
===========

- Shalev-shwartz, Shai, et al. "Pegasos: Primal Estimated sub-GrAdient SOlver for SVM."
Proceedings of the 24th International Conference on Machine Learning (ICML-07). 2007.
DOI=10.1145/1273496.1273598, http://dx.doi.org/10.1145/1273496.1273598,
"""
immutable Pegasos <: SvmSolver
    penalize_bias::Bool
end

Pegasos(; penalize_bias::Bool = false) = Pegasos(penalize_bias)

# ==========================================================================
# Implementation of primal sub-gradient solver for SVMs
# This particular method is specialized for
# - linear kernel
# - primal solution
# - dense arrays
# - univariate prediction

function fit{TKernel<:ScalarProductKernel, TLoss<:Union{MarginBasedLoss, DistanceBasedLoss}, TReg<:L2Penalty, INTERCEPT}(
        spec::CSVM{TKernel, TLoss, TReg},
        X::StridedMatrix, y⃗::StridedVector,
        solver::Pegasos,
        predmodel::UvPredicton{INTERCEPT},
        problem::PrimalProblem;
        ftol::Real = 1.0e-4,
        iterations::Int = 10000,
        show_trace::Bool = false,
        callback::Union{Function,Void} = nothing,
        nargs...)
    islipschitzcont_deriv(spec.loss) || throw(ArgumentError("The loss has to have a lipschitz continuous derivative"))

    # The risk functional needs to be strongly convex,
    # but since we limit the regularizer to L2Penalty,
    # it suffices to test the loss for convexity.
    isconvex(spec.loss) || throw(ArgumentError("The loss has to be convex"))

    # Get the size of the design matrix
    #   k ... number of features
    #   m ... number of observations
    k, m = size(X)

    # Do some housekeeping: use shorter variable names for the options
    lambda = Float64(1 / (spec.C * m))
    has_callback = typeof(callback) <: Function
    bias = INTERCEPT ? Float64(predmodel.bias) : zero(Float64)
    loss = spec.loss
    reg = L2Penalty(lambda)
    predictor = LinearPredictor(bias)
    risk = EmpiricalRisk(predictor, loss, reg, solver.penalize_bias)

    # Print log header if show_trace is set
    show_trace && print_iter_head()

    # Initialize w to all zero.
    # Note: w has length k+1 to have space for the potential intercept.
    #       If no intercept is fit, the bias term will be ignored
    w  = zeros(Float64, INTERCEPT ? k+1 : k)

    # Storage for the current gradient
    ▽  = zeros(Float64, INTERCEPT ? k+1 : k, 1)
    ▽⃗  = view(▽, 1:length(▽), 1)

    # Indicies into the observations of X (columns) for each iteration
    # They are random (uniform) for optimal average convergence.
    S = rand(1:m, iterations)

    # Precompute constants
    #   1/√(λ)
    one_sqrtlambda = 1 / sqrt(lambda)
    #   -1/λ
    minus_one_lambda = -1 / lambda

    iteration = 0; stopped = false
    @inbounds for (iteration, i) in enumerate(S)
        if stopped; break; end

        # Vectorview into the current observation i
        # We buffer it because we need it twice
        x⃗ᵢ = slice(X, :, i)

        # Compute current step size η
        #   -ηₜ = -1/(λ⋅t)
        minus_eta = minus_one_lambda / iteration

        # Compute prediction for the current observation
        #   ŷ = h(x⃗ᵢ,w⃗ₜ)
        ŷ = value(predictor, x⃗ᵢ, w) # Real number !

        # Compute the subgradient for the current observation
        #   ▽ₜ = λ⋅w⃗ₜ + L'(yᵢ, h(x⃗ᵢ,w⃗ₜ))⋅x⃗ᵢ
        grad!(▽, risk, x⃗ᵢ, w, y⃗[i], ŷ)

        # Update the coefficients using the the current subgradient
        #   w⃗ₜ₊₁ = w⃗ₜ - ηₜ⋅▽⃗ₜ
        axpy!(minus_eta, ▽⃗, w) # axpy!(a,X,Y) <=> Y = Y + a⋅X

        # Project the coefficients w on a ball
        #   nrm = (1/√(λ)) / ||w⃗ₜ₊₁||
        nrm = one_sqrtlambda / vecnorm(w, 2)
        if nrm < 1
            # w⃗ₜ₊₁ = min(1, nrm) ⋅ w⃗ₜ₊₁
            broadcast!(*, w, w, nrm)
        end

        # In case the user requested to print the learning process and/or provided
        # callback function this code will provide just that.
        # Note: If no callback function is provided then the compiler will be able
        #       to optimize the empty callback call away.
        if has_callback || show_trace
            # compute the risk of the current observation.
            f = value(risk, X, w, y⃗[i], ŷ)
            show_trace && print_iter(iteration, f, vecnorm(w, 2))
            stopped = _docallback(predmodel, callback, iteration, w, f, ▽)
        end
    end

    # Compute final objective value (risk) for the whole trainingset
    f = value(risk, X, w, y⃗)
    PrimalSolution(w, f, iteration, false)
end


text:
using AlgebraicPetri
using AlgebraicPetri.Epidemiology

using OrdinaryDiffEq
using LabelledArrays
using Plots

using Catlab.Theories
using Catlab.Graphics
using Catlab.CategoricalAlgebra
using Catlab.Programs.RelationalPrograms

display_uwd(ex) = to_graphviz(ex, box_labels=:name, junction_labels=:variable, edge_attrs=Dict(:len=>".75"));

sir = @relation (s, i, r) where (s, i, r) begin
    infection(s, i)
    recovery(i, r)
end
display_uwd(sir)

#-
p_sir = apex(oapply_epi(sir));
Graph(p_sir)

u0 = LVector(S=990, I=10, R=0);
t_span = (17.0,120.0)

γ = 1/14
β = t->begin
    policy_days = [20,60,120] .+ 17
    contact_rate = 0.05
    pol = findfirst(x->t<=x, policy_days) # array of days when policy changes
    growth_rate = pol == 1 ? 0.0 : (2^(1/((pol-1)*5)) - 1) # growth rate depending on policy
    return (growth_rate + γ) / 990 * (1-contact_rate) # calculate rate of infection
end
p = LVector(inf=β, rec=γ);

prob = ODEProblem(vectorfield(p_sir),u0,t_span,p)
sol = OrdinaryDiffEq.solve(prob,Tsit5())
plot(sol)
png("ode-chime.png")

max_stars_repo_name: AlgebraicJulia/AlgebraicPetri

max_issues_repo_name: AlgebraicJulia/AlgebraicPetri

max_forks_repo_name: AlgebraicJulia/AlgebraicPetri


text:
function beep(sound_type::String="alert")
    basename = joinpath(@__DIR__, "..", "sounds")

    path = joinpath(basename, string(sound_type, ".wav"))
    return(wavplay(path))
end

macro beep(st1, st2, expr)
    beep(st1)
    e = eval(expr)
    beep(st2)
    return e
end

macro beep(sound_type, expr)
    e = eval(expr)
    beep(sound_type)
    return e
end

macro beep(expr)
    e = eval(expr)
    beep()
    return e
end


text:
@testset "NoiseGrid" begin

  using DiffEqNoiseProcess, DiffEqBase, Test

  t = 0:0.001:1
  grid = exp.(t)
  W = NoiseGrid(t,grid)

  dt = 0.1
  calculate_step!(W,dt,nothing,nothing)

  for i in 1:10
    accept_step!(W,dt,nothing,nothing)
  end

  W = NoiseGrid(t,grid)
  prob = NoiseProblem(W,(0.0,1.0))
  sol = solve(prob;dt=0.1)

  @test !sol.step_setup
  @test_throws ErrorException accept_step!(sol,dt,nothing,nothing)

  t = 0:0.001:1
  grid = [[exp.(ti) for i in 1:8] for ti in t]
  W = NoiseGrid(t,grid)
  prob = NoiseProblem(W,(0.0,1.0))
  sol = solve(prob;dt=0.1)

  dt = 0.001
  t = 0:dt:1
  brownian_values = cumsum([0;[sqrt(dt)*randn() for i in 1:length(t)-1]])
  W = NoiseGrid(t,brownian_values)

  dt = 0.001
  t = 0:dt:1
  brownian_values2 = cumsum([[zeros(8)];[sqrt(dt)*randn(8) for i in 1:length(t)-1]])
  W = NoiseGrid(t,brownian_values2)
  prob = NoiseProblem(W,(0.0,1.0))
  sol = solve(prob;dt=0.1)

  dt = 1//1000
  t = 0:dt:1
  W = NoiseGrid(t,brownian_values)
  prob_rational = NoiseProblem(W,(0,1))
  sol = solve(prob_rational; dt=1//10)
end


text:
mutable struct PlaneWaveMWTD{T,F,P} <: TDFunctional{T}
  direction::P
  polarisation::P
  speedoflight::T
  amplitude::F
end

function planewave(polarisation,direction,amplitude,speedoflight)
    PlaneWaveMWTD(direction,polarisation,speedoflight,amplitude)
end

planewave(;signature, polarization, direction, speedoflight) =
    PlaneWaveMWTD(direction, polarization, speedoflight, signature)

scalartype(p::PlaneWaveMWTD) = eltype(p.polarisation)

*(a, pw::PlaneWaveMWTD) = PlaneWaveMWTD(
    pw.direction,
    a * pw.polarisation,
    pw.speedoflight,
    pw.amplitude
)

cross(k, pw::PlaneWaveMWTD) = PlaneWaveMWTD(
    pw.direction,
    k × pw.polarisation,
    pw.speedoflight,
    pw.amplitude
)

function (f::PlaneWaveMWTD)(r,t)
    t = cartesian(t)[1]
    dr = zero(typeof(t))
    for i in 1 : 3
        dr += r[i]*f.direction[i]
    end
    f.polarisation * f.amplitude(t - dr/f.speedoflight)
end

function integrate(f::BEAST.PlaneWaveMWTD)
    planewave(
        signature = integrate(f.amplitude),
        direction = f.direction,
        polarization = f.polarization,
        speedoflight = f.speedoflight)
end