text:
"""
```julia
vpercentile!(A, p; dims)
```
Compute the `p`th percentile (where `p ∈ [0,100]`) of all elements in `A`,
optionally over dimensions specified by `dims`.

As `StatsBase.percentile`, but in-place, slightly vectorized, and supporting the `dims` keyword.

Be aware that, like `Statistics.median!`, this function modifies `A`, sorting or
partially sorting the contents thereof (specifically, along the dimensions specified
by `dims`, using either `quicksort!` or `quickselect!` depending on the size
of the array). Do not use this function if you do not want the contents of `A`
to be rearranged.

If reducing over multiple `dims`, these dimensions must be contiguous (i.e.
`dims=(2,3)` but not `dims=(1,3)`). Note also that specifying `dims` other than
`:` creates `view`s, with some nonzero performance cost.

## Examples
```julia
julia> using VectorizedStatistics

julia> A = [1 2 3; 4 5 6; 7 8 9]
3×3 Matrix{Int64}:
 1  2  3
 4  5  6
 7  8  9

 julia> vpercentile!(A, 50, dims=1)
 1×3 Matrix{Float64}:
  4.0  5.0  6.0

 julia> vpercentile!(A, 50, dims=2)
 3×1 Matrix{Float64}:
  2.0
  5.0
  8.0

 julia> vpercentile!(A, 50)
 5.0

 julia> A # Note that the array has been sorted
3×3 Matrix{Int64}:
 1  4  7
 2  5  8
 3  6  9
```
"""
vpercentile!(A, p::Number; dims=:) = _vquantile!(A, p/100, dims)
export vpercentile!

"""
```julia
vquantile!(A, q; dims)
```
Compute the `q`th quantile (where `q ∈ [0,1]`) of all elements in `A`,
optionally over dimensions specified by `dims`.

Similar to `StatsBase.quantile!`, but slightly vectorized, and supporting the `dims` keyword.

Be aware that, like `StatsBase.quantile!`, this function modifies `A`, sorting or
partially sorting the contents thereof (specifically, along the dimensions specified
by `dims`, using either `quicksort!` or `quickselect!` depending on the size
of the array). Do not use this function if you do not want the contents of `A`
to be rearranged.

If reducing over multiple `dims`, these dimensions must be contiguous (i.e.
`dims=(2,3)` but not `dims=(1,3)`). Note also that specifying `dims` other than
`:` creates `view`s, with some nonzero performance cost.

## Examples
```julia
julia> using VectorizedStatistics

julia> A = [1 2 3; 4 5 6; 7 8 9]
3×3 Matrix{Int64}:
 1  2  3
 4  5  6
 7  8  9

 julia> vquantile!(A, 0.5, dims=1)
 1×3 Matrix{Float64}:
  4.0  5.0  6.0

 julia> vquantile!(A, 0.5, dims=2)
 3×1 Matrix{Float64}:
  2.0
  5.0
  8.0

 julia> vquantile!(A, 0.5)
 5.0

 julia> A # Note that the array has been sorted
3×3 Matrix{Int64}:
 1  4  7
 2  5  8
 3  6  9
```
"""
vquantile!(A, q::Number; dims=:) = _vquantile!(A, q, dims)
export vquantile!

# Reduce one dim
_vquantile!(A, q::Real, dims::Int) = _vquantile!(A, q, (dims,))

# Reduce some dims
function _vquantile!(A::AbstractArray{T,N}, q::Real, dims::Tuple) where {T,N}
    iscontiguous(dims) || error("Only continuous `dims` are currently supported")
    sᵢ = size(A)
    sₒ = ntuple(Val(N)) do d
        ifelse(d ∈ dims, 1, sᵢ[d])
    end
    Tₒ = Base.promote_op(/, T, Int)
    B = similar(A, Tₒ, sₒ)
    _vquantile!(B, A, q, dims)
end

# Reduce all the dims!
_vquantile!(A, q::Real, ::Tuple{Colon}) = _vquantile!(A, q, :)
function _vquantile!(A, q::Real, ::Colon)
    iₗ, iᵤ = firstindex(A), lastindex(A)
    A, iₗ, iᵤ₋ = sortnans!(A, iₗ, iᵤ)
    if iᵤ₋ < iᵤ
        # Remove this if we'd rather ignore NaNs
        return A[iᵤ]
    end
    N₋ = iᵤ - iₗ
    iₚ = q*N₋ + iₗ
    iₚ₋ = floor(Int, iₚ)
    iₚ₊ = ceil(Int, iₚ)
    if N₋ < 384
        quicksort!(A, iₗ, iᵤ)
    else
        quickselect!(A, iₗ, iᵤ, iₚ₋)
        quickselect!(A, iₚ₊, iᵤ, iₚ₊)
    end
    f = iₚ - iₚ₋
    return f*A[iₚ₊] + (1-f)*A[iₚ₋]
end

# Generate customized set of loops for a given ndims and a vector
# `static_dims` of dimensions to reduce over
function staticdim_quantile_quote(static_dims::Vector{Int}, N::Int)
  M = length(static_dims)
  # `static_dims` now contains every dim we're taking the quantile over.
  Bᵥ = Expr(:call, :view, :B)
  Aᵥ = Expr(:call, :view, :A)
  reduct_inds = Int[]
  nonreduct_inds = Int[]
  # Firstly, build our expressions for indexing each array
  Aind = :(Aᵥ[])
  Bind = :(Bᵥ[])
  inds = Vector{Symbol}(undef, N)
  for n ∈ 1:N
    ind = Symbol(:i_,n)
    inds[n] = ind
    if n ∈ static_dims
      push!(reduct_inds, n)
      push!(Aᵥ.args, :)
      push!(Bᵥ.args, :(firstindex(B,$n)))
    else
      push!(nonreduct_inds, n)
      push!(Aᵥ.args, ind)
      push!(Bᵥ.args, :)
      push!(Bind.args, ind)
    end
  end
  # Secondly, build up our set of loops
  if !isempty(nonreduct_inds)
    firstn = first(nonreduct_inds)
    block = Expr(:block)
    loops = Expr(:for, :($(inds[firstn]) = indices((A,B),$firstn)), block)
    if length(nonreduct_inds) > 1
      for n ∈ @view(nonreduct_inds[2:end])
        newblock = Expr(:block)
        push!(block.args, Expr(:for, :($(inds[n]) = indices((A,B),$n)), newblock))
        block = newblock
      end
    end
    rblock = block
    # Push more things here if you want them at the beginning of the reduction loop
    push!(rblock.args, :(Aᵥ = $Aᵥ))
    push!(rblock.args, :($Bind = _vquantile!(Aᵥ, q, :)))
    # Put it all together
    return quote
      Bᵥ = $Bᵥ
      @inbounds $loops
      return B
    end
  else
    return quote
      Bᵥ = $Bᵥ
      Bᵥ[] = _vquantile!(A, q, :)
      return B
    end
  end
end

# Chris Elrod metaprogramming magic:
# Turn non-static integers in `dims` tuple into `StaticInt`s
# so we can construct `static_dims` vector within @generated code
function branches_quantile_quote(N::Int, M::Int, D)
  static_dims = Int[]
  for m ∈ 1:M
    param = D.parameters[m]
    if param <: StaticInt
      new_dim = _dim(param)::Int
      @assert new_dim ∉ static_dims
      push!(static_dims, new_dim)
    else
      t = Expr(:tuple)
      for n ∈ static_dims
        push!(t.args, :(StaticInt{$n}()))
      end
      q = Expr(:block, :(dimm = dims[$m]))
      qold = q
      ifsym = :if
      for n ∈ 1:N
        n ∈ static_dims && continue
        tc = copy(t)
        push!(tc.args, :(StaticInt{$n}()))
        qnew = Expr(ifsym, :(dimm == $n), :(return _vquantile!(B, A, q, $tc)))
        for r ∈ m+1:M
          push!(tc.args, :(dims[$r]))
        end
        push!(qold.args, qnew)
        qold = qnew
        ifsym = :elseif
      end
      # Else, if dimm ∉ 1:N, drop it from list and continue
      tc = copy(t)
      for r ∈ m+1:M
        push!(tc.args, :(dims[$r]))
      end
      push!(qold.args, Expr(:block, :(return _vquantile!(B, A, q, $tc))))
      return q
    end
  end
  return staticdim_quantile_quote(static_dims, N)
end

# Efficient @generated in-place quantile
@generated function _vquantile!(B::AbstractArray{Tₒ,N}, A::AbstractArray{T,N}, q::Real, dims::D) where {Tₒ,T,N,M,D<:Tuple{Vararg{Integer,M}}}
  branches_quantile_quote(N, M, D)
end
@generated function _vquantile!(B::AbstractArray{Tₒ,N}, A::AbstractArray{T,N}, q::Real, dims::Tuple{}) where {Tₒ,T,N}
  :(copyto!(B, A); return B)
end


text:
export hessian_operator

function hessian_operator(nlp, x::Vector{T}) where T
    n = nlp.meta.nvar
    temp = similar(x)
    return hess_op!(nlp, x, temp)
end

ext: jl

max_stars_repo_name: vepiteski/LSDescent

max_issues_repo_name: vepiteski/LSDescent

max_forks_repo_name: vepiteski/LSDescent


text:
using Weave, Compat
using Base.Test

cleanup = true

#Test chunk options and output formats
weave("documents/chunk_options.noweb", plotlib=nothing)
result = @compat readstring(open("documents/chunk_options.md"))
ref = @compat readstring(open("documents/chunk_options_ref.md"))
@test result == ref
cleanup && rm("documents/chunk_options.md")

weave("documents/chunk_options.noweb", doctype="tex", plotlib=nothing)
result = @compat readstring(open("documents/chunk_options.tex"))
ref = @compat readstring(open("documents/chunk_options_ref.tex"))
@test result == ref
cleanup && rm("documents/chunk_options.tex")

weave("documents/chunk_options.noweb", doctype="texminted", plotlib=nothing)
result = @compat readstring(open("documents/chunk_options.tex"))
ref = @compat readstring(open("documents/chunk_options_ref.texminted"))
@test result == ref
cleanup && rm("documents/chunk_options.tex")

weave("documents/chunk_options.noweb", doctype="rst", plotlib=nothing)
result = @compat readstring(open("documents/chunk_options.rst"))
ref = @compat readstring(open("documents/chunk_options_ref.rst"))
@test result == ref
cleanup && rm("documents/chunk_options.rst")

#Test out_path
weave("documents/chunk_options.noweb", doctype="rst",
      out_path="documents/outpath_options.rst" , plotlib=nothing)
result = @compat readstring(open("documents/outpath_options.rst"))
ref = @compat readstring(open("documents/chunk_options_ref.rst"))
@test result == ref
cleanup && rm("documents/outpath_options.rst")

#Test tangle
tangle("documents/chunk_options.noweb")
result = @compat readstring(open("documents/chunk_options.jl"))
ref = @compat readstring(open("documents/chunk_options_ref.jl"))
@test result == ref
cleanup && rm("documents/chunk_options.jl")

tangle("documents/chunk_options.noweb", out_path = "documents/outoptions.jl")
result = @compat readstring(open("documents/outoptions.jl"))
ref = @compat readstring(open("documents/chunk_options_ref.jl"))
@test result == ref
cleanup && rm("documents/outoptions.jl")

#Test functions and sandbox clearing
weave("documents/chunk_func.noweb", plotlib=nothing)
result = @compat readstring(open("documents/chunk_func.md"))
ref = @compat readstring(open("documents/chunk_func_ref.md"))
@test result == ref
cleanup && rm("documents/chunk_func.md")


text:
"""
    Adapt.@adapt_structure T

Define a method `adapt_structure(to, obj::T)` which calls `adapt_structure` on each field
of `obj` and constructs a new instance of `T` using the default constuctor `T(...)`.
"""
macro adapt_structure(T)
    names = fieldnames(Core.eval(__module__, T))
    quote
        function Adapt.adapt_structure(to, obj::$(esc(T)))
            $(esc(T))($([:(Adapt.adapt_structure(to, obj.$name)) for name in names]...))
        end
    end
end


text:
function pythagorean_triplets(N)
    results = []
    # a < b < c
    for a in 1:(N ÷ 3 - 1)
        for b in (a + 1):(N - 2a - 2)
            c = N - a - b
            # if a + b + c > N
            #     break
            # end
            if a^2 + b^2 == c^2
                push!(results, (a, b, c))
            end
        end
    end
    return results
end

max_stars_repo_name: twolodzko/Learning

max_issues_repo_name: twolodzko/Learning

max_forks_repo_name: twolodzko/Learning


text:
module PointGroup

using Molecules

tol = 1E-5

include("sea.jl")
include("moit.jl")
include("symmetry_finding.jl")
include("flowchart.jl")

end


text:
PROJECT_DIR = pwd()*"\\.."

using DataFrames
using Dates
#using HDF5, JLD

include("Params.jl")
include("Solution.jl")
include("LocalSearch.jl")
include("GeneticAlgorithm.jl")

#allinstances = readdir("../../instances/");
#allsolutions = readdir("../../solutions/");

#println("Number of instance files: ",length(allinstances))
#println("Number of solution files: ",length(allsolutions))

df = DataFrame(
    instance = String[],
    size = Int64[],
    seed = Int64[],
    finalCost = Int64[],
    executionTime = Float64[],
    numTrials = Int64[],
    bestKnownCost = Any[],
    gap = Any[]
)

timenow = Dates.now()
start_datetime = Dates.format(timenow, "yyyy-mm-dd_HHMMSS")
experiment_dir = joinpath(PROJECT_DIR,"experiments",start_datetime)
mkdir(experiment_dir)

allinstances = [
    "wil100.dat","tho40.dat","ste36c.dat","tai60b.dat","tai80b.dat",
    "sko49.dat","sko100d.dat","scr20.dat","nug24.dat","nug30.dat",
    "bur26a.dat","bur26e.dat","esc32e.dat","tai150b.dat","tho150.dat",
    "esc16h.dat","esc128.dat","chr15a.dat","chr20c.dat","wil50.dat",
    "lipa40b.dat","lipa50a.dat","lipa90b.dat","wil100.dat","els19.dat"
]

results_file = "aggregate_results.txt"

cpu_time = 120

cols_tuple = Tuple(string(sym) for sym in describe(df).variable)
open(joinpath(experiment_dir, results_file),"a+") do txtfile
   write(txtfile,join(cols_tuple,"\t"),"\n")
end

n_seed = 2:5
num_trials = 10^20

for seed = n_seed
    for instance_name in allinstances
        println("$instance_name; seed=$seed")
        instance_path = joinpath(PROJECT_DIR, "instances",instance_name)
        output_path = joinpath(experiment_dir,split(instance_name,'.')[1]*".jld")
        params = Params(instance_path, output_path, seed, cpu_time)

        Random.seed!(seed)
        sol, results = multistartLS(params, num_trials)

        println("LS solution:")
        printSolution(sol)
        println("Gap=",results["gap"])
        println()

        if occursin(" ", params.pathToOutput)
            throw("Your project path must not contain any spaces, because of the behaviour of the JLD package.")
        end
        #save(params.pathToOutput, "data", results) # for some unknown reason this does not work
        #save(joinpath("..\\..","experiments",start_datetime,split(instance_name,'.')[1]*".jld"), "data", results)

        new_tuple = (instance_name,
            params.datasetSize,
            seed,
            results["finalCost"],
            results["time"],
            results["numTrials"],
            params.solutionCost,
            results["gap"])
        push!(df, new_tuple)

        open(joinpath(experiment_dir, results_file),"a+") do txtfile
           write(txtfile,join(new_tuple,"\t"),"\n")
        end
    end
end

max_stars_repo_name: BSAraujo/QAP-Metaheuristics

max_issues_repo_name: BSAraujo/QAP-Metaheuristics

max_forks_repo_name: BSAraujo/QAP-Metaheuristics


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

mutable struct ResourceSkuCosts <: SwaggerModel
    meterID::Any # spec type: Union{ Nothing, String } # spec name: meterID
    quantity::Any # spec type: Union{ Nothing, Int64 } # spec name: quantity
    extendedUnit::Any # spec type: Union{ Nothing, String } # spec name: extendedUnit

    function ResourceSkuCosts(;meterID=nothing, quantity=nothing, extendedUnit=nothing)
        o = new()
        validate_property(ResourceSkuCosts, Symbol("meterID"), meterID)
        setfield!(o, Symbol("meterID"), meterID)
        validate_property(ResourceSkuCosts, Symbol("quantity"), quantity)
        setfield!(o, Symbol("quantity"), quantity)
        validate_property(ResourceSkuCosts, Symbol("extendedUnit"), extendedUnit)
        setfield!(o, Symbol("extendedUnit"), extendedUnit)
        o
    end
end # type ResourceSkuCosts

const _property_map_ResourceSkuCosts = Dict{Symbol,Symbol}(Symbol("meterID")=>Symbol("meterID"), Symbol("quantity")=>Symbol("quantity"), Symbol("extendedUnit")=>Symbol("extendedUnit"))
const _property_types_ResourceSkuCosts = Dict{Symbol,String}(Symbol("meterID")=>"String", Symbol("quantity")=>"Int64", Symbol("extendedUnit")=>"String")
Base.propertynames(::Type{ ResourceSkuCosts }) = collect(keys(_property_map_ResourceSkuCosts))
Swagger.property_type(::Type{ ResourceSkuCosts }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_ResourceSkuCosts[name]))}
Swagger.field_name(::Type{ ResourceSkuCosts }, property_name::Symbol) =  _property_map_ResourceSkuCosts[property_name]

function check_required(o::ResourceSkuCosts)
    true
end

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


text:
using Romeo, GLFW, GLAbstraction, Reactive, ModernGL, GLWindow, Color
function clear!(x::Vector{RenderObject})
    while !isempty(x)
        value = pop!(x)
        delete!(value)
    end
end
global const gaggaaa = Texture("pic.jpg")

function dropequal(a::Signal)
    is_equal = foldl((false, a.value), a) do v0, v1
        (v0[2] == v1, v1)
    end
    dropwhen(lift(first, is_equal), a.value, a)
end
global MemRefs = Any[]

function init_romeo()

    sourcecode_area = lift(Romeo.ROOT_SCREEN.area) do x
    	Rectangle(0, 0, div(x.w, 7)*3, x.h)
    end
    visualize_area = lift(Romeo.ROOT_SCREEN.area) do x
        Rectangle(div(x.w,7)*3, 0, div(x.w, 7)*3, x.h)
    end
    search_area = lift(visualize_area) do x
        Rectangle(x.x, x.y, x.w, div(x.h,10))
    end
    edit_area = lift(Romeo.ROOT_SCREEN.area) do x
    	Rectangle(div(x.w, 7)*6, 0, div(x.w, 7), x.h)
    end

    sourcecode_screen   = Screen(Romeo.ROOT_SCREEN, area=sourcecode_area)
    visualize_screen    = Screen(Romeo.ROOT_SCREEN, area=visualize_area)
    search_screen       = Screen(visualize_screen,  area=search_area)
    edit_screen         = Screen(Romeo.ROOT_SCREEN, area=edit_area)

    w_height = lift(Romeo.ROOT_SCREEN.area) do x
    	x.h
    end
    source_offset = lift(w_height) do x
        translationmatrix(Vec3(30,x-30,0))
    end
    w_height_search = lift(search_screen.area) do x
        x.h
    end
    search_offset = lift(w_height_search) do x
        translationmatrix(Vec3(30,x-30,0))
    end

    searchinput = Input("barplot")

    const sourcecode  = visualize("barplot = Float32[(sin(i/10f0) + cos(j/2f0))/4f0 \n for i=1:10, j=1:10]\n", model=source_offset, screen=sourcecode_screen)
    global const barplot = Float32[(sin(i/10f0) + cos(j/2f0))/4f0 for i=1:10, j=1:10]
    search            = visualize("barplot\n", model=search_offset, color=rgba(0.9,0,0.2,1), screen=search_screen)

    viz, source_text = edit(sourcecode[:text], sourcecode)
    viz, search_text = edit(search[:text], search)

    should_eval = dropequal(lift(Romeo.ROOT_SCREEN.inputs[:buttonspressed]) do keyset
        keyset == IntSet(GLFW.KEY_ENTER, GLFW.KEY_LEFT_CONTROL)
    end)

    soursupdate = lift(source_text, should_eval) do source, seval
        expr = parse(strip(source), raise=false)
        if seval
            if expr.head != :error
                try
                    eval(Main, expr)
                catch e
                    println(e)
                end
            else
                #println(expr)
            end
        end
        nothing
    end
    a = keepwhen(should_eval, nothing, soursupdate)
    lift(search_text, a) do x, _
        s = symbol(strip(x))
        if isdefined(s)
            value = eval(Main, s)
            if applicable(visualize, value)
                clear!(visualize_screen.renderlist)
                clear!(edit_screen.renderlist)
                push!(MemRefs, value)
                obj     = visualize(value, screen=visualize_screen)
                objedit = edit(obj,        screen=edit_screen)

                push!(visualize_screen.renderlist, obj)
                append!(edit_screen.renderlist, objedit)
            end
        end
        nothing
    end

    push!(sourcecode_screen.renderlist, sourcecode)
    push!(search_screen.renderlist, search)
    glClearColor(0,0,0,0)
end

init_romeo()
while Romeo.ROOT_SCREEN.inputs[:open].value
    Romeo.renderloop(Romeo.ROOT_SCREEN)
    sleep(0.0001)
end
GLFW.Terminate()


text:
module Day2

import DelimitedFiles

function read_input(filepath)
    return DelimitedFiles.readdlm(filepath, ' ', '\n')
end

function calculate_position_and_depth(course)
    position = 0
    depth = 0
    for row = 1:(size(course)[1])
        direction = course[row, 1]
        magnitude = course[row, 2]
        if direction == "forward"
            position += magnitude
        elseif direction == "up"
            depth -= magnitude
        elseif direction == "down"
            depth += magnitude
        end
    end
    return (position, depth)
end

function part1_answer()
    input_filepath = joinpath(@__DIR__, "input.txt")
    position, depth = calculate_position_and_depth(read_input(input_filepath))
    return position * depth
end

function calculate_position_and_depth_with_aim(course)
    position = 0
    depth = 0
    aim = 0
    for row = 1:(size(course)[1])
        direction = course[row, 1]
        magnitude = course[row, 2]
        if direction == "forward"
            position += magnitude
            depth += aim*magnitude
        elseif direction == "up"
            aim -= magnitude
        elseif direction == "down"
            aim += magnitude
        end
    end
    return (position, depth)
end

function part2_answer()
    input_filepath = joinpath(@__DIR__, "input.txt")
    position, depth = calculate_position_and_depth_with_aim(read_input(input_filepath))
    return position * depth
end

end


text:
using Test, Distributed

# benchmark uses several workers if available
@everywhere begin

    using BlackBoxOptimizationBenchmarking, Optim
    import BlackBoxOptimizationBenchmarking: minimizer, minimum, optimize
    const BBOB = BlackBoxOptimizationBenchmarking

    #interface for Optim
    pinit(D) = 10*rand(D).-5
    optimize(opt::Optim.AbstractOptimizer,f,D,run_length) =
        Optim.optimize(f, pinit(D), opt, Optim.Options(f_calls_limit=run_length,g_tol=1e-12))

    string(opt::Optim.AbstractOptimizer) = string(typeof(opt).name)
end

b = BBOB.benchmark(
    NelderMead(), BBOB.F1, [100,500,1000], 20, 2, 1e-6,
)

b = BBOB.benchmark(
    NelderMead(), 1, [100,500,1000], 20, 2, 1e-6,
)

b = BBOB.benchmark(
    [NelderMead(), GradientDescent()], [1, 2], [100,500,1000], 20, 2, 1e-6,
)


text:
function ComItAlg(A,B,Q,R,U,bet,n1,n2,cutoff=1.0)
#ComItAlg   Solves for optimization problem in commitment case.
#
#
#  Usage:     (M,C) = ComItAlg(A,B,Q,R,U,bet,n1,n2,cutoff)
#
#  Input:     A        nxn matrix, (n=n1+n2)
#             B        nxk matrix
#             Q        nxn matrix, symmetric
#             R        kxk matrix, symmetric
#             U        nxk matrix
#             bet      scalar, discount factor (eg 0.99)
#             n1       scalar, # of predetermined variables
#             n2       scalar, # of forward looking variables
#             cutoff   scalar, max modulus of backward looking variables
#
#  Output:    M        nxn matrix, [x1(t+1),p2(t+1)] = M*[x1(t),p2(t)] + e(t+1)
#             C        (n2+k+n1)xn matrix, [x2(t),u(t),p1(t)] = C*[x1(t),p2(t)]
#
#
#
#  Details:   Solve for optimization problem in commitment case. The economy
#             evolves as
#
#             x1(t+1)       =   A * x1(t)    + Bu(t) + e(t+1)
#             E(t)x2(t+1)           x2(t)              0
#
#             where x1(t) ("backward looking") has n1 elements,
#             and x2(t) ("forward looking") has n2 elements. The initial
#             value of x1, x1(0), is given by history. Let n=n1+n2.
#             u(t) is an kx1 vector of decision variables, which follows
#             the decision rule.
#
#             Define x(t)=[x1(t),x2(t)]. The loss function of the policy
#             maker is
#
#             Sum{ (bet^t)*[x(t)'Q*x(t)+2x(t)'U*u(t)+u(t)'R*u(t)],t=0,1,... }
#
#
#
#  Paul Söderlind, Paul.Soderlind@unisg.ch, Aug 2000, to Julia Jan 2016
#----------------------------------------------------------------------------

  Q = (Q + Q')/2                #to make symmetric
  R = (R + R')/2

  n = n1 + n2
  k = size(R,1)                #No. of control variables

  G =  [ [I             zeros(n,k)    zeros(n,n)   ];
         [zeros(n,n)    zeros(n,k)    (bet*A')     ];
         [zeros(k,n)    zeros(k,k)    (-B')        ]]

  D =  [ [A             B             zeros(n,n)   ];
         [(-bet*Q)      (-bet*U)      I            ];
         [U'            R             zeros(k,n)   ]]

  v1 = 1:n1                                 #reordering
  v2 = n1+1:n1+n2
  v3 = n1+n2+1:n1+n2+k+n1
  v4 = n1+n2+k+n1+1:n1+n2+k+n1+n2
  v  = [v1;v4;v2;v3]

  G = G[:,v]          #x1,x2,(u,p1),p2 -> x1,p2,x2,(u,p1)
  D = D[:,v]

  GS = schur(G,D)
  logconA = abs.(GS.beta) .<= (abs.(GS.alpha)*cutoff) #selecting stable eigenvalues
  ordschur!(GS,logconA)                                   #reordering, stable first

  if sum(logconA) < n
    @warn("Too few stable roots: no stable solution")
    return NaN,NaN
  elseif sum(logconA) > n
    @warn("Too many stable roots: inifite number of stable solutions")
    return NaN,NaN
  end

  Stt = GS.S[1:n,1:n]
  Ttt = GS.T[1:n,1:n]
  Zkt = GS.Z[1:n,1:n]
  Zlt = GS.Z[n+1:n+k+n,1:n]

  if cond(Zkt) > 1e+14
    @warn("Zkt is singular: rank condition for solution not satisfied")
    return NaN, NaN
  end

  Zkt_1 = inv(Zkt)         #inverting
  Stt_1 = inv(Stt)

  M = real(Zkt*Stt_1*Ttt*Zkt_1)      #[x1(t+1),p2(t+1)] = M*[x1(t),p2(t)]+e(t+1)
  C = real(Zlt*Zkt_1)                #[x2(t),u(t),p1(t)] =C*[x1(t),p2(t)]

  return M,C

end
#-----------------------------------------------------------------------

max_stars_repo_name: PaulSoderlind/PaulSoderlindCode

max_issues_repo_name: PaulSoderlind/PaulSoderlindCode

max_forks_repo_name: PaulSoderlind/PaulSoderlindCode


text:
struct SL2NormHandler end

gethandler(broadcasted, ::Type{typeof(LinearAlgebra.norm)}, ::Type{<:AbstractVector{<:SBitstream}}) =
    !broadcasted ? SL2NormHandler() :
                   error("Cannot generate hardware for broadcasted L2 norm.")
init_state(::SL2NormHandler) = (id = 0,)

function (handler::SL2NormHandler)(buffer, netlist, state, inputs, outputs)
    # update netlist with inputs
    setsigned!(netlist, inputs[1], true)

    # add output net to netlist
    setsigned!(netlist, outputs[1], true)

    write(buffer, """
        $stdcomment
        // BEGIN l2norm$(state.id)
        stoch_l2_norm #(
                .VEC_LEN($(netsize(inputs[1])[1]))
            ) l2norm$(state.id) (
                .CLK  (CLK),
                .nRST (nRST),
                .up   ($(name(inputs[1]))_p),
                .un   ($(name(inputs[1]))_m),
                .yp   ($(name(outputs[1]))_p),
                .yn   ($(name(outputs[1]))_m)
            );
        // END l2norm$(state.id)
        \n""")

    return buffer, (id = state.id + 1,)
end


text:
using Test
using SumOfSquares
using DynamicPolynomials

# Inspired from https://github.com/jump-dev/SumOfSquares.jl/issues/79
function concave_then_convex_cubic_test(optimizer, config::MOIT.TestConfig,
                                        MCT::Type{<:MOI.AbstractVectorSet})
    atol = config.atol
    rtol = config.rtol

    model = _model(optimizer)

    @polyvar x
    @variable(model, p, Poly(monomials(x, 0:3)))
    cone = ConvexPolyInnerCone{MCT}()
    cref_convex  = @constraint(model,  p in cone, domain = (@set x >= 0))
    cref_concave = @constraint(model, -p in cone, domain = (@set x <= 0))
    @constraint(model, p(x =>  2) ==  8)
    @constraint(model, p(x => -2) == -8)
    @constraint(model, p(x =>  1) ==  1)
    @constraint(model, p(x => -1) == -1)

    optimize!(model)

    @test termination_status(model) == MOI.OPTIMAL

    @test primal_status(model) == MOI.FEASIBLE_POINT
    @test value(p) ≈ x^3 atol=atol rtol=rtol
    @test_throws SumOfSquares.ValueNotSupported value(cref_convex)
    @test_throws SumOfSquares.ValueNotSupported value(cref_concave)

    # The monomials contain the variables created for the Hessian so we cannot
    # check them easily
    for μ in [dual(cref_convex), dual(cref_concave)]
        @test μ isa AbstractMeasure{Float64}
        @test length(moments(μ)) == 2
    end
    for μ in [moments(cref_convex), moments(cref_concave)]
        @test μ isa AbstractMeasure{Float64}
        @test length(moments(μ)) == 7
    end
end

function sos_concave_then_convex_cubic_test(optimizer, config)
    concave_then_convex_cubic_test(optimizer, config,
                                   MOI.PositiveSemidefiniteConeTriangle)
end
sd_tests["sos_concave_then_convex_cubic"] = sos_concave_then_convex_cubic_test
function sdsos_concave_then_convex_cubic_test(optimizer, config)
    concave_then_convex_cubic_test(optimizer, config,
                                   SumOfSquares.ScaledDiagonallyDominantConeTriangle)
end
soc_tests["sdsos_concave_then_convex_cubic"] = sdsos_concave_then_convex_cubic_test
function dsos_concave_then_convex_cubic_test(optimizer, config)
    concave_then_convex_cubic_test(optimizer, config,
                                   SumOfSquares.DiagonallyDominantConeTriangle)
end
linear_tests["dsos_concave_then_convex_cubic"] = dsos_concave_then_convex_cubic_test


text:
using Knet, Statistics, Random

# Defining the convolutional layer:
struct Conv; w; b; f; p; end
(c::Conv)(x) = c.f.(pool(conv4(c.w, dropout(x,c.p)) .+ c.b))
Conv(w1::Int,w2::Int,cx::Int,cy::Int,f=relu;pdrop=0) = Conv(param(w1,w2,cx,cy), param0(1,1,cy,1), f, pdrop)

# Defining dense layer
struct Dense; w; b; f; p; end
(d::Dense)(x) = d.f.(d.w * mat(dropout(x,d.p)) .+ d.b) # mat reshapes 4-D tensor to 2-D matrix so we can use matmul
Dense(i::Int,o::Int,f=relu;pdrop=0) = Dense(param(o,i), param0(o), f, pdrop)

struct Chain
    layers
    Chain(layers...) = new(layers)
end
(c::Chain)(x) = (for l in c.layers; x = l(x); end; x)
(c::Chain)(x,y) = nll(c(x),y)

struct BaselineRandomAgent
    chain::Chain
end
(b::BaselineRandomAgent)(img) = b.chain(img)
(b::BaselineRandomAgent)(img, actions) = sum(abs.(b(img) .- actions))


text:
d = SparseCat([1,2], [0.5, 0.5])
@test occursin("SparseCat", sprint(showdistribution, d))

d = SparseCat(1:50, fill(1/50, 50))
iob = IOBuffer()
io = IOContext(iob, :limit=>true, :displaysize=>(10, 7))
showdistribution(io, d)
str = String(take!(iob))
@test occursin("SparseCat", str)
@test occursin("<everything else>", str)

# test that it doesn't print <everything else> when there are enough lines
d = SparseCat([:a], 1.0)
iob = IOBuffer()
io = IOContext(iob, :limit=>true, :displaysize=>(10, 7))
showdistribution(io, d)
str = String(take!(iob))
@test occursin("SparseCat", str)
@test !occursin("<everything else>", str)


text:
using PyPlot
using Colors

export Circle, CircleAgentSpecification,
  get_block,
  mean_coverage,
  mean_area_coverage

# evaluation of the area coverage objective
struct Circle
  center::Array{Float64, 1}
  radius::Float64
end

struct CircleAgentSpecification
  sensor_radius
  station_radius
  num_sensors
end

# returns a vector of vectors of element indices abstracting the partition
# matroid

point_in_circle(center, radius, point) = norm(point - center) < radius
point_in_circle(circle::Circle, point) =
  point_in_circle(circle.center, circle.radius, point)

function point_covered(circles, x, y)
  #optimized point_in_circle code runs about ten times faster
  @inbounds for circle in circles
    if (x - circle.center[1])^2 + (y - circle.center[2])^2 < circle.radius^2
      return true
    end
  end

  false
end

mean_coverage(circles, points) =
  mean(map(point->point_covered(circles, point), points))

function mean_area_coverage(circles, discretization)
  num_covered = 0
  for x in range(0, stop=1, length=discretization),
    y in range(0, stop=1, length=discretization)

    num_covered += point_covered(circles, x, y)
  end

  num_covered / (discretization * discretization)
end

# agent generation
function rand_in_circle()
  r = sqrt(rand())
  theta = 2*pi*rand()
  r * [cos(theta);sin(theta)]
end

function make_agent(agent_specification::CircleAgentSpecification)
  agent_center = rand(2)
  make_agent(agent_specification, agent_center)
end
function make_agent(agent_specification::CircleAgentSpecification,
                    agent_center::Vector{Float64})
  make_sensor() = Circle(agent_center + agent_specification.station_radius *
                        rand_in_circle(),
                        agent_specification.sensor_radius)

  sensors = [make_sensor() for agent in 1:agent_specification.num_sensors]

  Agent(agent_center, agent_specification.station_radius, sensors)
end

function plot_circle(circle::Circle; x...)
  plot_circle(circle.center; radius = circle.radius, x...)
end

function plot_element(circle::Circle; color = "k", alpha = 0.3)
  center = circle.center
  scatter([center[1]], [center[2]], color = color, s = 4*agent_scale,
          marker = sensor_center, alpha = alpha)
end

function plot_agent(agent; color = "k")
  center = agent.center
  scatter([center[1]], [center[2]], color = color, s = 9*agent_scale,
          marker = "o", alpha = 0.3)

  plot_circle(Circle(agent.center, agent.radius);
              linestyle=station_style, color = color,
              linewidth = 3.0)
end

function plot_filled_circle(circle::Circle; color = "k", alpha = 0.3)
  ps = circle_points(circle.center, scale = 1, radius = circle.radius)
  full_circle = [ps ps[:,1]]

  ret = []

  append!(ret, fill(ps[1,:][:], ps[2,:][:], color = color, alpha=alpha, linewidth=0.0))

  append!(ret, plot_circle(circle, color="k", linewidth = 1.0))

  ret
end

function visualize_solution(circles::Array{Circle}, colors)
  map(circles, colors) do circle, color
    plot_filled_circle(circle, color = rgb_tuple(color), alpha=0.2)

    center = circle.center
    scatter([center[1]], [center[2]], color = rgb_tuple(color), s = 4*agent_scale,
            marker = selected_sensor, edgecolors="k")
  end
  nothing
end

max_stars_repo_name: mcorah/MultiAgentSensing

max_issues_repo_name: mcorah/MultiAgentSensing

max_forks_repo_name: mcorah/MultiAgentSensing


text:
using ModelingToolkit, StaticArrays, LinearAlgebra, LabelledArrays
using DiffEqBase, ForwardDiff
using Test

# Define some variables
@parameters t σ ρ β
@variables x(t) y(t) z(t)
D = Differential(t)

# Define a differential equation
eqs = [D(x) ~ σ * (y - x),
    D(y) ~ t * x * (ρ - z) - y,
    D(z) ~ x * y - β * z]

@named de = ODESystem(eqs)
ff = ODEFunction(de, [x, y, z], [σ, ρ, β], jac = true)

a = @SVector [1.0, 2.0, 3.0]
b = SLVector(x = 1.0, y = 2.0, z = 3.0)
c = [1.0, 2.0, 3.0]
p = SLVector(σ = 10.0, ρ = 26.0, β = 8 / 3)
@test ff(a, p, 0.0) isa SVector
@test typeof(ff(b, p, 0.0)) <: SLArray
@test ff(c, p, 0.0) isa Vector
@test ff(a, p, 0.0) == ff(b, p, 0.0)
@test ff(a, p, 0.0) == ff(c, p, 0.0)

@test ff.jac(a, p, 0.0) isa SMatrix
@test typeof(ff.jac(b, p, 0.0)) <: SMatrix
@test ff.jac(c, p, 0.0) isa Matrix
@test ff.jac(a, p, 0.0) == ff.jac(b, p, 0.0)
@test ff.jac(a, p, 0.0) == ff.jac(c, p, 0.0)

# Test similar_type
@test ff(b, p, ForwardDiff.Dual(0.0, 1.0)) isa SLArray
d = LVector(x = 1.0, y = 2.0, z = 3.0)
@test ff(d, p, ForwardDiff.Dual(0.0, 1.0)) isa LArray
@test ff.jac(b, p, ForwardDiff.Dual(0.0, 1.0)) isa SArray
@test eltype(ff.jac(b, p, ForwardDiff.Dual(0.0, 1.0))) <: ForwardDiff.Dual
@test ff.jac(d, p, ForwardDiff.Dual(0.0, 1.0)) isa Array
@inferred ff.jac(d, p, ForwardDiff.Dual(0.0, 1.0))
@test eltype(ff.jac(d, p, ForwardDiff.Dual(0.0, 1.0))) <: ForwardDiff.Dual


text:
using SimpleDrawing, LinearFractionalTransformations

export Circle

"""
    Circle

Represents a circle in the plane. Constructors:
* `circle(x,y,r)` based on center coordinates and radius
* `circle(z,r)` based on center given as a complex number and a radius
* `circle(a,b,c)` based on three complex numbers that represent points on the circle
"""
struct Circle <: Cline
    x::Float64
    y::Float64
    r::Float64
    function Circle(a::Real, b::Real, rad::Real)
        rad >= 0 || error("Radius must be nonnegative")
        new(a, b, rad)
    end
end

Circle(z::Complex, r::Real) = Circle(real(z), imag(z), r)

function Circle(a::Complex, b::Complex, c::Complex)
    z = find_center(a, b, c)
    !isinf(z) || error("Cannot construct a circle through three collinear points")
    r = abs(a - z)
    return Circle(z, r)
end

# copy constructor
Circle(C::Circle) = Circle(C.x, C.y, C.r)

# standard stuff

export radius, area, circumference, center

"""
    radius(C)
Return the radius of the circle.
"""
radius(C::Circle) = C.r

"""
    area(C)
Return the area of the circle.
"""
area(C::Circle) = π * C.r * C.r

"""
    circumference(C)
Return the circumference of the circle.
"""
circumference(C::Circle) = 2 * π * C.r

import Plots: center
"""
    center(C)
Return the center of the circle as a complex number.
"""
center(C::Circle) = Complex(C.x, C.y)

"""
    three_points(C)
Return a list of three points on the circle `C`.
"""
function three_points(C::Circle)
    z = center(C)
    r = radius(C)
    return [z + r * exp(k * (2π * im / 3)) for k = 0:2]
end
export three_points

import SimpleDrawing: draw
const _FILL_COLOR = :lightgray
const _LINE_COLOR = :black

"""
    draw(C::Circle, fill::Bool=false; args...)
Draw the circle `C`. With `fill=true`, fill in the interior in color.

Example:
`draw(C,true,color=:yellow, linecolor=:red, fillalpha=0.3)`
"""
function draw(C::Circle, fill::Bool = false; args...)
    if fill
        draw_disc(
            center(C),
            radius(C);
            linecolor = _LINE_COLOR,
            color = _FILL_COLOR,
            args...,
        )
    else
        draw_circle(center(C), radius(C); linecolor = _LINE_COLOR, args...)
    end
end

include("operations.jl")
include("kiss.jl")
include("inversion.jl")


text:
const default_index = """<!DOCTYPE html>
<html>
    <head>
        {%metas%}
        <title>{%title%}</title>
        {%favicon%}
        {%css%}
    </head>
    <body>
        {%app_entry%}
        <footer>
            {%config%}
            {%scripts%}
            {%renderer%}
        </footer>
    </body>
</html>"""

"""
    struct DashApp <: Any

Dash.jl's internal representation of a Dash application.

This `struct` is not intended to be called directly; developers should create their Dash application using the `dash` function instead.
"""
mutable struct DashApp
    root_path ::String
    is_interactive ::Bool
    config ::DashConfig
    index_string ::Union{String, Nothing}
    title ::String
    layout ::Union{Nothing, Component, Function}
    devtools ::DevTools
    callbacks ::Dict{Symbol, Callback}
    inline_scripts ::Vector{String}

    DashApp(root_path, is_interactive, config, index_string, title = "Dash") =
    new(root_path, is_interactive, config, index_string, title, nothing, DevTools(dash_env(Bool, "debug", false)), Dict{Symbol, Callback}(), String[])

end

#only name, index_string and layout are available to set
function Base.setproperty!(app::DashApp, property::Symbol, value)
    property == :name && return set_name!(app, value)
    property == :index_string && return set_index_string!(app, value)
    property == :layout && return set_layout!(app::DashApp, value)
    property == :title && return set_title!(app::DashApp, value)

    property in fieldnames(DashApp) && error("The property `$(property)` of `DashApp` is read-only")

    error("The property `$(property)` of `DashApp` does not exist.")
end

function set_name!(app::DashApp, name)
    setfield!(app, :name, name)
end

function set_title!(app::DashApp, title)
    setfield!(app, :title, title)
end

get_name(app::DashApp) = app.name

function set_layout!(app::DashApp, component::Union{Component,Function})
    setfield!(app, :layout, component)
end

get_layout(app::DashApp) = app.layout

function check_index_string(index_string::AbstractString)
    validate_index(
        "index_string", index_string,
        [
            "{%app_entry%}"=>r"{%app_entry%}",
            "{%config%}"=>r"{%config%}",
            "{%scripts%}"=>r"{%scripts%}"
        ]
    )
end
function set_index_string!(app::DashApp, index_string::AbstractString)
    check_index_string(index_string)
    setfield!(app, :index_string, index_string)
end

get_index_string(app::DashApp) = app.component.index_string

"""
Activate the dev tools, called by `run_server`.

If a parameter can be set by an environment variable, that is listed as:
  env: `DASH_****`
  Values provided here take precedence over environment variables.

Available dev_tools environment variables:
- `DASH_DEBUG`
- `DASH_UI`
- `DASH_PROPS_CHECK`
- `DASH_SERVE_DEV_BUNDLES`
- `DASH_HOT_RELOAD`
- `DASH_HOT_RELOAD_INTERVAL`
- `DASH_HOT_RELOAD_WATCH_INTERVAL`
- `DASH_HOT_RELOAD_MAX_RETRY`
- `DASH_SILENCE_ROUTES_LOGGING`
- `DASH_PRUNE_ERRORS`

# Arguments

- `debug::Bool` - Enable/disable all the dev tools unless overridden by the
       arguments or environment variables. Default is ``true`` when
       ``enable_dev_tools`` is called directly, and ``false`` when called
       via ``run_server``. env: ``DASH_DEBUG``.

- `dev_tools_ui::Bool` - Show the dev tools UI. env: ``DASH_UI``

- `dev_tools_props_check::Bool` - Validate the types and values of Dash
       component props. env: ``DASH_PROPS_CHECK``

- `dev_tools_serve_dev_bundles::Bool` - Serve the dev bundles. Production
    bundles do not necessarily include all the dev tools code.
    env: ``DASH_SERVE_DEV_BUNDLES``

- `dev_tools_hot_reload::Bool` - Activate hot reloading when app, assets,
    and component files change. env: ``DASH_HOT_RELOAD``

- `dev_tools_hot_reload_interval::Float64` - Interval in seconds for the
    client to request the reload hash. Default 3.
    env: ``DASH_HOT_RELOAD_INTERVAL``

- `dev_tools_hot_reload_watch_interval::Float64` - Interval in seconds for the
    server to check asset and component folders for changes.
    Default 0.5. env: ``DASH_HOT_RELOAD_WATCH_INTERVAL``

- `dev_tools_hot_reload_max_retry::Int` - Maximum number of failed reload
    hash requests before failing and displaying a pop up. Default 8.
    env: ``DASH_HOT_RELOAD_MAX_RETRY``
"""
function enable_dev_tools!(app::DashApp; debug = nothing,
            dev_tools_ui = nothing,
            dev_tools_props_check = nothing,
            dev_tools_serve_dev_bundles = nothing,
            dev_tools_hot_reload = nothing,
            dev_tools_hot_reload_interval = nothing,
            dev_tools_hot_reload_watch_interval = nothing,
            dev_tools_hot_reload_max_retry = nothing,
            dev_tools_silence_routes_logging = nothing,
            dev_tools_prune_errors = nothing)
    @env_default!(debug, Bool, true)
    setfield!(app, :devtools, DevTools(
        debug;
        ui = dev_tools_ui,
        props_check = dev_tools_props_check,
        serve_dev_bundles = dev_tools_serve_dev_bundles,
        hot_reload = dev_tools_hot_reload,
        hot_reload_interval = dev_tools_hot_reload_interval,
        hot_reload_watch_interval = dev_tools_hot_reload_watch_interval,
        hot_reload_max_retry = dev_tools_hot_reload_max_retry,
        silence_routes_logging = dev_tools_silence_routes_logging,
        prune_errors = dev_tools_prune_errors
    ))
end

get_devsetting(app::DashApp, name::Symbol) = getproperty(app.devtools, name)

get_setting(app::DashApp, name::Symbol) = getproperty(app.config, name)

get_assets_path(app::DashApp) = joinpath(app.root_path, get_setting(app, :assets_folder))

"""
    dash(;
            external_stylesheets,
            external_scripts,
            url_base_pathname,
            requests_pathname_prefix,
            routes_pathname_prefix,
            assets_folder,
            assets_url_path,
            assets_ignore,
            serve_locally,
            suppress_callback_exceptions,
            eager_loading ,
            meta_tags,
            index_string,
            assets_external_path,
            include_assets_files,
            show_undo_redo,
            compress,
            update_title
        )

Dash is a framework for building analytical web applications. No JavaScript required.

If a parameter can be set by an environment variable, that is listed as:
  env: `DASH_****`
  Values provided here take precedence over environment variables.

# Arguments
- `assets_folder::String` - a path, relative to the current working directory,
        for extra files to be used in the browser. Default `'assets'`. All .js and .css files will be loaded immediately unless excluded by `assets_ignore`, and other files such as images will be served if requested.

- `assets_url_path::String` - The local urls for assets will be:
        ``requests_pathname_prefix * assets_url_path * "/" * asset_path``
        where ``asset_path`` is the path to a file inside ``assets_folder``.
        Default ``'assets'`.

- `assets_ignore::String` - [EXPERIMENTAL] A regex, as a string to pass to ``Regex``, for
        assets to omit from immediate loading. Ignored files will still be
        served if specifically requested. You cannot use this to prevent access
        to sensitive files.
    :type assets_ignore: string

-  `assets_external_path::String` - [EXPERIMENTAL] an absolute URL from which to load assets.
        Use with ``serve_locally=false``. Dash can still find js and css to
        automatically load if you also keep local copies in your assets
        folder that Dash can index, but external serving can improve
        performance and reduce load on the Dash server.
        env: `DASH_ASSETS_EXTERNAL_PATH`

- `include_assets_files::Bool` - [EXPERIMENTAL] Default ``true``, set to ``false`` to prevent
        immediate loading of any assets. Assets will still be served if
        specifically requested. You cannot use this to prevent access
        to sensitive files.
        env: `DASH_INCLUDE_ASSETS_FILES`

- `url_base_pathname::String`: A local URL prefix to use app-wide.
        Default ``nothing``. Both `requests_pathname_prefix` and
        `routes_pathname_prefix` default to `url_base_pathname`.
        env: `DASH_URL_BASE_PATHNAME`

- `requests_pathname_prefix::String`: A local URL prefix for file requests.
        Defaults to `url_base_pathname`, and must end with
        `routes_pathname_prefix`
        env: `DASH_REQUESTS_PATHNAME_PREFIX`

- `routes_pathname_prefix::String`: A local URL prefix for JSON requests.
        Defaults to ``url_base_pathname``, and must start and end
        with ``'/'``.
        env: `DASH_ROUTES_PATHNAME_PREFIX`

- `serve_locally`: [EXPERIMENTAL] If `true` (default), assets and dependencies (Dash and Component js and css) will be served from local URLs. If `false` Dash will use CDN links where available.
        (Dash and Component js and css) will be served from local URLs.
        If ``false`` we will use CDN links where available.

- `meta_tags::Vector{Dict{String, String}}`: html <meta> tags to be added to the index page.
        Each dict should have the attributes and values for one tag, eg:
        ``Dict("name"=>"description", "content" => "My App")``

- `index_string::String`: Override the standard Dash index page.
        Must contain the correct insertion markers to interpolate various
        content into it depending on the app config and components used.
        See https://dash.plotly.com/external-resources for details.

- `external_scripts::Vector`: Additional JS files to load with the page.
        Each entry can be a String (the URL) or a Dict{String, String} with ``src`` (the URL)
        and optionally other ``<script>`` tag attributes such as ``integrity``
        and ``crossorigin``.

- `external_stylesheets::Vector`: Additional CSS files to load with the page.
        Each entry can be a String (the URL) or a Dict{String, String} with ``href`` (the URL)
        and optionally other ``<link>`` tag attributes such as ``rel``,
        ``integrity`` and ``crossorigin``.

- `suppress_callback_exceptions::Bool`: Default ``false``: check callbacks to
        ensure referenced IDs exist and props are valid. Set to ``true``
        if your layout is dynamic, to bypass these checks.
        env: `DASH_SUPPRESS_CALLBACK_EXCEPTIONS`

- `prevent_initial_callbacks::Bool`: Default ``false``: Sets the default value
        of ``prevent_initial_call`` for all callbacks added to the app.
        Normally all callbacks are fired when the associated outputs are first
        added to the page. You can disable this for individual callbacks by
        setting ``prevent_initial_call`` in their definitions, or set it
        ``true`` here in which case you must explicitly set it ``false`` for
        those callbacks you wish to have an initial call. This setting has no
        effect on triggering callbacks when their inputs change later on.

- `show_undo_redo::Bool`: Default ``false``, set to ``true`` to enable undo
        and redo buttons for stepping through the history of the app state.

- `compress::Bool`: Default ``true``, controls whether gzip is used to compress
        files and data served by HTTP.jl when supported by the client. Set to
        ``false`` to disable compression completely.

- `update_title::String`: Default ``Updating...``. Configures the document.title
        (the text that appears in a browser tab) text when a callback is being run.
        Set to '' if you don't want the document.title to change or if you
        want to control the document.title through a separate component or
        clientside callback.
"""
function dash(;
        external_stylesheets = ExternalSrcType[],
        external_scripts  = ExternalSrcType[],
        url_base_pathname = dash_env("url_base_pathname"),
        requests_pathname_prefix = dash_env("requests_pathname_prefix"),
        routes_pathname_prefix = dash_env("routes_pathname_prefix"),
        assets_folder = "assets",
        assets_url_path = "assets",
        assets_ignore = "",
        serve_locally = true,
        suppress_callback_exceptions = dash_env(Bool, "suppress_callback_exceptions", false),
        prevent_initial_callbacks = false,
        eager_loading = false,
        meta_tags = Dict{Symbol, String}[],
        index_string = default_index,
        assets_external_path = dash_env("assets_external_path"),
        include_assets_files = dash_env(Bool, "include_assets_files", true),
        show_undo_redo = false,
        compress = true,
        update_title = "Updating..."

    )

        check_index_string(index_string)
        config = DashConfig(
            external_stylesheets,
            external_scripts,
            pathname_configs(
                url_base_pathname,
                requests_pathname_prefix,
                routes_pathname_prefix
                )...,
            assets_folder,
            lstrip(assets_url_path, '/'),
            assets_ignore,
            serve_locally,
            suppress_callback_exceptions,
            prevent_initial_callbacks,
            eager_loading,
            meta_tags,
            assets_external_path,
            include_assets_files,
            show_undo_redo,
            compress,
            update_title
        )
        result = DashApp(app_root_path(), isinteractive(), config, index_string)
    return result
end


text:
summary_callback(integrator) = false # when used as condition; never call the summary callback during the simulation
summary_callback(u, t, integrator) = u_modified!(integrator, false) # the summary callback does nothing when called accidentally

"""
    SummaryCallback()

Create and return a callback that prints a human-readable summary of the simulation setup at the
beginning of a simulation and then resets the timer. When the returned callback is executed
directly, the current timer values are shown.
"""
function SummaryCallback()
  DiscreteCallback(summary_callback, summary_callback,
                   save_positions=(false,false),
                   initialize=initialize_summary_callback)
end

function Base.show(io::IO, cb::DiscreteCallback{<:Any, <:typeof(summary_callback)})
  @nospecialize cb # reduce precompilation time

  print(io, "SummaryCallback")
end

# Format a key/value pair for output from the SummaryCallback
function format_key_value_line(key::AbstractString, value::AbstractString, key_width, total_width;
                               indentation_level=0, guide='…', filler='…', prefix="│ ", suffix=" │")
  @assert key_width < total_width
  line  = prefix
  # Indent the key as requested (or not at all if `indentation_level == 0`)
  indentation = prefix^indentation_level
  reduced_key_width = key_width - length(indentation)
  squeezed_key = indentation * squeeze(key, reduced_key_width, filler=filler)
  line *= squeezed_key
  line *= ": "
  short = key_width - length(squeezed_key)
  if short <= 1
    line *= " "
  else
    line *= guide^(short-1) * " "
  end
  value_width = total_width - length(prefix) - length(suffix) - key_width - 2
  squeezed_value = squeeze(value, value_width, filler=filler)
  line *= squeezed_value
  short = value_width - length(squeezed_value)
  line *= " "^short
  line *= suffix

  @assert length(line) == total_width "should not happen: algorithm error!"

  return line
end
format_key_value_line(key, value, args...; kwargs...) = format_key_value_line(string(key), string(value), args...; kwargs...)

# Squeeze a string to fit into a maximum width by deleting characters from the center
function squeeze(message, max_width; filler::Char='…')
  @assert max_width >= 3 "squeezing works only for a minimum `max_width` of 3"

  length(message) <= max_width && return message

  keep_front = div(max_width, 2)
  keep_back  = div(max_width, 2) - (isodd(max_width) ? 0 : 1)
  remove_back  = length(message) - keep_front
  remove_front = length(message) - keep_back
  squeezed = (chop(message, head=0, tail=remove_back)
              * filler *
              chop(message, head=remove_front, tail=0))

  @assert length(squeezed) == max_width "`$(length(squeezed)) != $max_width` should not happen: algorithm error!"

  return squeezed
end

# Print a summary with a box around it with a given heading and a setup of key=>value pairs
function summary_box(io::IO, heading, setup=[])
  summary_header(io, heading)
  for (key, value) in setup
    summary_line(io, key, value)
  end
  summary_footer(io)
end

function summary_header(io, heading; total_width=100, indentation_level=0)
  total_width = get(io, :total_width, total_width)
  indentation_level = get(io, :indentation_level, indentation_level)

  @assert indentation_level >= 0 "indentation level may not be negative"

  # If indentation level is greater than zero, we assume the header has already been printed
  indentation_level > 0 && return

  # Print header
  println(io, "┌" * "─"^(total_width-2) * "┐")
  println(io, "│ " * heading * " "^(total_width - length(heading) - 4) * " │")
  println(io, "│ " * "═"^length(heading) * " "^(total_width - length(heading) - 4) * " │")
end

function summary_line(io, key, value; key_width=30, total_width=100, indentation_level=0)
  # Printing is not performance-critical, so we can use `@nospecialize` to reduce latency
  @nospecialize value # reduce precompilation time

  key_width = get(io, :key_width, key_width)
  total_width = get(io, :total_width, total_width)
  indentation_level = get(io, :indentation_level, indentation_level)

  s = format_key_value_line(key, value, key_width, total_width,
                            indentation_level=indentation_level)

  println(io, s)
end

function summary_footer(io; total_width=100, indentation_level=0)
  total_width = get(io, :total_width, 100)
  indentation_level = get(io, :indentation_level, 0)

  if indentation_level == 0
    s = "└" * "─"^(total_width-2) * "┘"
  else
    s = ""
  end

  print(io, s)
end

@inline increment_indent(io) = IOContext(io, :indentation_level => get(io, :indentation_level, 0) + 1)

# Print information about the current simulation setup
# Note: This is called *after* all initialization is done, but *before* the first time step
function initialize_summary_callback(cb::DiscreteCallback, u, t, integrator)

  mpi_isroot() || return nothing

  print_startup_message()

  io = stdout
  io_context = IOContext(io,
                         :compact => false,
                         :key_width => 30,
                         :total_width => 100,
                         :indentation_level => 0)

  semi = integrator.p
  show(io_context, MIME"text/plain"(), semi)
  println(io, "\n")
  mesh, equations, solver, _ = mesh_equations_solver_cache(semi)
  show(io_context, MIME"text/plain"(), mesh)
  println(io, "\n")
  show(io_context, MIME"text/plain"(), equations)
  println(io, "\n")
  show(io_context, MIME"text/plain"(), solver)
  println(io, "\n")

  callbacks = integrator.opts.callback
  if callbacks isa CallbackSet
    for cb in callbacks.continuous_callbacks
      show(io_context, MIME"text/plain"(), cb)
      println(io, "\n")
    end
    for cb in callbacks.discrete_callbacks
      # Do not show ourselves
      cb.affect! === summary_callback && continue

      show(io_context, MIME"text/plain"(), cb)
      println(io, "\n")
    end
  else
    show(io_context, MIME"text/plain"(), callbacks)
    println(io, "\n")
  end

  setup = Pair{String,Any}[
           "Start time" => first(integrator.sol.prob.tspan),
           "Final time" => last(integrator.sol.prob.tspan),
           "time integrator" => integrator.alg |> typeof |> nameof,
          ]
  summary_box(io, "Time integration", setup)
  println()

  reset_timer!(timer())

  return nothing
end

function (cb::DiscreteCallback{Condition,Affect!})(io::IO=stdout) where {Condition, Affect!<:typeof(summary_callback)}

  mpi_isroot() || return nothing

  TimerOutputs.complement!(timer())
  print_timer(io, timer(), title="Trixi.jl",
              allocations=true, linechars=:unicode, compact=false)
  println(io)
  return nothing
end


text:
# bias_variance-arwmh-logistic
#
# This script compares the Stein discrepancy of sample sequences drawn
# from asymptotically biased approximate random walk Metropolis-Hastings
# (ARWMH) chains (http://arxiv.org/pdf/1304.5299v4.pdf) and asymptotically
# exact RWMH chains.  By plotting the Stein discrepancy results against
# a measure of computation required per sample (e.g., the number of
# datapoint likelihood evaluations needed), we can quantify a bias-variance
# trade-off underlying the approximate MCMC scheme.
#
# The target distribution is a Bayesian logistic regression posterior
# under independent Gaussian priors based on a prostate cancer dataset with
# binary outcomes indicating whether cancer has spread to surrounding lymph
# nodes.

include("experiment_utils.jl")

# Set random number generator seed
@parseintcli seed "s" "seed" 37

using RDatasets
using StatsFuns: logistic
using SteinDiscrepancy: gsd
using SteinDistributions: SteinLogisticRegressionGaussianPrior, gradlogdensity
using SteinSamplers: runapproxmh

# Select a solver for Stein discrepancy optimization problem
#solver = "clp"
solver = "gurobi"

# Setup the data
data_set = "nodal"
nodal = dataset("boot", "nodal")
y = 2 * convert(Array{Float64},nodal[:R]) - 1
# X includes an intercept term
X = convert(Array{Float64, 2}, nodal[[:M, :Aged, :Stage, :Grade, :Xray, :Acid]])

# Parameter dimension
d = size(X,2)
# the max number of samples to run in the solver
n = 1000
# Associate sample weights
q = ones(n, 1) ./ n
# name of prior + posterior
distname = "logisticgaussianprior"
# setup distribution
target = SteinLogisticRegressionGaussianPrior(X, y)
# define gradlogp
function gradlogp(x::Array{Float64,1})
    gradlogdensity(target, x)
end
# intiial guess for beta
beta0 = zeros(d)
# burnin likelihoods
burninlikelihoods = 1e3
# max number of likelihoods
maxlikelihoods = 1e5
# batchsize for different algorithms [nodal set is ~50 rows]
batchsize = 2
# the epsilons to run the experiemnt
epsilons = [0.0, 0.1, 0.2]
# Sample sizes at which Stein discrepancy will be computed
ns = vcat(10:10:min(400,n), 420:20:min(600,n), 650:50:min(1000,n))

for epsilon = epsilons
    @printf("Generating samples for epsilon=%f\n", epsilon)

    (betas, numlikelihoods) =
        @setseed seed runapproxmh(target, beta0;
                                  epsilon=epsilon,
                                  numlikelihood=maxlikelihoods,
                                  batchsize=batchsize,
                                  proposalvariance=0.01)

    # Thin to n samples
    firstn = find([(l >= burninlikelihoods) for l in numlikelihoods])[1]
    lastn = length(numlikelihoods)
    idx = round(Int, linspace(firstn, lastn, min(lastn - firstn + 1, n)))

    # Keep track of predictive probabilities (to use as test functions for
    # external validation of bias-variance trade-off)
    predprobs = logistic(X * betas[idx,:]')
    scaledpredprobs = predprobs ./ vec(maximum(X, 2))

    for i in ns
        println("Beginning optimization for n=$(i), d=$(d), epsilon=$(epsilon)")
        # Compute Stein discrepancy for first i points
        betai = betas[idx[1:i],:]
        res = gsd(points=betai, gradlogdensity=gradlogp, solver=solver)
        println("\tn = $(i), objective = $(res.objectivevalue)")

        # Package and save results
        edges = res.edges
        edge_pairs = [vec(edges[i, :]) for i in 1:size(edges, 1)]

        instance_data = Dict{Any, Any}(
            "n" => i,
            "distname" => distname,
            "d" => d,
            "seed" => seed,
            "epsilon" => epsilon,
            "batchsize" => batchsize,
            "X" => betai,
            "xnorm" => res.xnorm,
            "edges" => edge_pairs,
            "g" => res.g,
            "gradg" => res.gradg,
            "objectivevalue" => vec(res.objectivevalue),
            "scaled.pred.probs" => vec(mean(scaledpredprobs[:,1:i], 2)),
            "numlikelihood" => numlikelihoods[idx[i]],
            "edgetime" => res.edgetime,
            "solvetime" => res.solvetime,
        )

        save_json(
            instance_data;
            dir="bias_variance-arwmh-logistic",
            distname=distname,
            dataset=data_set,
            n=i,
            burninlikelihoods=int(burninlikelihoods),
            epsilon=epsilon,
            batchsize=batchsize,
            d=d,
            seed=seed
        )
    end
end


text:
export fsort, fsort!

using SortingLab: fsortperm
using Tables

if VERSION >= v"1.3.0"
    import Base.Threads: @spawn
else
    macro spawn(_)
        println("DataConvenience: multithreading do not work in < Julia 1.3")
    end
end

using Base.Threads: @spawn

fsort(tbl, col::Symbol) = fsort(tbl, [col])
fsort!(tbl, col::Symbol) = fsort!(tbl, [col])

fsort(tbl, cols; threaded=true) = fsort!(copy(tbl), cols; threaded=threaded)

function fsort!(tbl, cols; threaded=true)
    @assert Tables.columnaccess(tbl)

    if threaded && (VERSION >= v"1.3")
        return _fsort_parallel!(tbl, cols)
    else
        return _fsort_single!(tbl, cols)
    end
end

function _fsort_parallel!(tbl, cols)
    @assert VERSION >= v"1.3"
    for col in reverse(cols)
        ordering = fsortperm(tbl[!, col])
        channel_lock = Channel{Bool}(length(names(tbl)))
        for c in names(tbl)
            @spawn begin
                v = tbl[!, c]
                @inbounds v .= v[ordering]
                put!(channel_lock, true)
            end
        end
        for _ in names(tbl)
            take!(channel_lock)
        end
    end
    tbl
end

function _fsort_single!(tbl, cols)
    for col in reverse(cols)
        ordering = fsortperm(tbl[!, col])
        for c in names(tbl)
            v = tbl[!, c]
            @inbounds v .= v[ordering]
        end
    end
    tbl
end

max_stars_repo_name: xiaodaigh/DataConvenience

max_issues_repo_name: xiaodaigh/DataConvenience

max_forks_repo_name: xiaodaigh/DataConvenience


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

@recipe function f(domain::PointSet{T,N}, data::AbstractVector) where {N,T}
  X = coordinates(domain)

  seriestype --> :scatter

  if N == 1
    X[1,:], data
  elseif N == 2
    aspect_ratio --> :equal
    marker_z --> data
    colorbar --> true
    X[1,:], X[2,:]
  elseif N == 3
    aspect_ratio --> :equal
    marker_z --> data
    colorbar --> true
    X[1,:], X[2,:], X[3,:]
  else
    @error "cannot plot in more than 3 dimensions"
  end
end

@recipe function f(domain::PointSet{T,N}) where {N,T}
  X = coordinates(domain)

  seriestype --> :scatter
  seriescolor --> :black
  legend --> false

  if N == 1
    @series begin
      X[1,:], fill(0, nelms(domain))
    end
  elseif N == 2
    aspect_ratio --> :equal
    @series begin
      X[1,:], X[2,:]
    end
  elseif N == 3
    aspect_ratio --> :equal
    @series begin
      X[1,:], X[2,:], X[3,:]
    end
  else
    @error "cannot plot in more than 3 dimensions"
  end
end


text:
# EXCLUDE FROM TESTING
using VectorEngine

mutable struct xm
      x::Int32
      m::Int64
end

@noinline function _pass_struct!(r::xm)
    @veprintf("r.x=%d, r.m=%ld\n", r.x, r.m)
    r.m = r.m + 1
    @veprintf("r.m=%ld\n", r.m)
end

function pass_struct!(r)
    _pass_struct!(r)
    return
end

a = xm(1, 100)

@veda pass_struct!(a)
synchronize()
@show a


text:
using Test, FinancialMonteCarlo, DualNumbers;
@show "Black Scholes Model"
S0 = 100.0;
K = 100.0;
T = 1.0;
d = 0.01;
D = 90.0;

Nsim = 10000;
Nstep = 30;
sigma = 0.2;
mc = MonteCarloConfiguration(Nsim, Nstep);
dr__ = 1e-7;
mc1 = MonteCarloConfiguration(Nsim, Nstep, FinancialMonteCarlo.AntitheticMC());
toll = 0.8
der_toll = 1e-2;

rfCurve = ZeroRate(dual(0.02, 1.0));
rfCurve1 = ZeroRate(0.02);

FwdData = Forward(T)
EUData = EuropeanOption(T, K)
AMData = AmericanOption(T, K)
BarrierData = BarrierOptionDownOut(T, K, D)
AsianFloatingStrikeData = AsianFloatingStrikeOption(T)
AsianFixedStrikeData = AsianFixedStrikeOption(T, K)
Model = BlackScholesProcess(sigma, Underlying(S0, d));

@show EuPrice = pricer(Model, rfCurve, mc, EUData);
@show EuPrice_der = rho(Model, rfCurve1, mc, EUData, dr__);

@show AmPrice = pricer(Model, rfCurve, mc, AMData);
@show AmPrice_der = rho(Model, rfCurve1, mc, AMData, dr__);

@test abs(EuPrice - 8.43005524824866) < toll
@test abs(EuPrice.epsilon - EuPrice_der) < der_toll

@test abs(AmPrice - 8.450489415187354) < toll
@test abs(AmPrice.epsilon - AmPrice_der) < der_toll

@show EuPrice = pricer(Model, rfCurve, mc1, EUData);
@show EuPrice_der = rho(Model, rfCurve1, mc1, EUData, dr__);
@show EuPrice_der_vec = rho(Model, rfCurve1, mc1, [EUData], dr__)[1];
tollanti = 0.6;

@test abs(EuPrice - 8.43005524824866) < tollanti
@test abs(EuPrice.epsilon - EuPrice_der) < der_toll
@test abs(EuPrice.epsilon - EuPrice_der_vec) < der_toll


text:
using BinaryBuilder
using Downloads
using LibGit2
import SHA: sha256

#
# Utilities
#

function purge()
    # remove old builds
    vers = TOML.parse(read(extra_versions_file(), String))
    rm(extra_versions_file())
    open(extra_versions_file(); append=true) do io
        for (ver, data) in vers
            @assert haskey(data, "file")
            @assert !haskey(data, "url")

            path = download_dir(data["file"])
            if round(now() - unix2datetime(mtime(path)), Dates.Week) < Dates.Week(1)
                println(io, "[\"$ver\"]")
                for (key, value) in data
                    println(io, "$key = $(repr(value))")
                end
                println(io)
            else
                rm(path)
                rm(path * ".sha256"; force=true)
            end
        end
    end
end

function hash_file(path)
    open(path, "r") do f
        bytes2hex(sha256(f))
    end
end

function installed_julia_dir(jp)
     jp_contents = readdir(jp)
     # Allow the unpacked directory to either be insider another directory (as produced by
     # the buildbots) or directly inside the mapped directory (as produced by the BB script)
     if length(jp_contents) == 1
         jp = joinpath(jp, first(jp_contents))
     end
     jp
end

#
# Versions
#

# a version is identified by a specific unique version, includes a hash to verify, and
# points to a local file or a remove resource.

read_versions() = merge(TOML.parse(read(versions_file(), String)),
                        TOML.parse(read(extra_versions_file(), String)))

"""
    prepare_julia(the_ver, dir=mktempdir())

Download and extract the specified version of Julia using the information provided in
`Versions.toml` to the directory `dir`.
"""
function prepare_julia(the_ver::VersionNumber, dir::String=mktempdir())
    vers = read_versions()
    for (ver, data) in vers
        ver == string(the_ver) || continue
        if haskey(data, "url")
            url = data["url"]

            file = get(data, "file", "julia-$ver.tar.gz")
            @assert !isabspath(file)
            file = download_dir(file)
            mkpath(dirname(file))

            Pkg.PlatformEngines.download_verify_unpack(url, data["sha"], dir;
                                                       tarball_path=file, force=true,
                                                       ignore_existence=true)
        else
            file = data["file"]
            !isabspath(file) && (file = download_dir(file))
            Pkg.PlatformEngines.verify(file, data["sha"])
            Pkg.PlatformEngines.unpack(file, dir)
        end
        return dir
    end
    error("Requested Julia version $the_ver not found")
end

function obtain_julia(spec::String)::VersionNumber
    try
        # maybe it's a named release
        # NOTE: check this first, because "1.3" could otherwise be interpreted as v"1.3.0"
        obtain_julia_release(spec)
    catch
        try
            # maybe it already refers to a version in Versions.toml
            version = VersionNumber(spec)
            prepare_julia(version)
            version
        catch
            # assume it points to something in our Git repository
            obtain_julia_build(spec)
        end
    end
end

#
# Releases
#

# a release is identified by name, points to an online resource, but does not have a version
# number. it will be downloaded, probed, and added to Versions.toml.

read_releases() = TOML.parsefile(releases_file())

# get the Julia version and short hash from a Julia installation
function get_julia_version(base)
    version = VersionNumber(read(`$(installed_julia_dir(base))/bin/julia -e 'print(Base.VERSION_STRING)'`, String))
    if version.prerelease != ()
        shorthash = read(`$(installed_julia_dir(base))/bin/julia -e 'print(Base.GIT_VERSION_INFO.commit_short)'`, String)
        if endswith(shorthash, '*')
            # strip the dirty marker
            shorthash = shorthash[1:end-1]
        end
        version = VersionNumber(string(version) * string("-", shorthash))
    end
    return version
end

"""
    version = obtain_julia_release(name::String)

Download Julia from an on-line source listed in Releases.toml as identified by `name`.
Returns the `version` (what other functions use to identify this build).
This version will be added to Versions.toml.
"""
function obtain_julia_release(name::String)
    releases = read_releases()
    @assert haskey(releases, name) "Julia release $name is not registered in Releases.toml"
    data = releases[name]

    # get the filename and extension from the url
    url = data["url"]
    filename = basename(url)
    if endswith(filename, ".tar.gz")
        ext = ".tar.gz"
        base = filename[1:end-7]
    else
        base, ext = splitext(filename)
    end

    # download
    filepath = download_dir(filename)
    duplicates = 0
    while ispath(filepath)
        duplicates += 1
        filepath = download_dir("$(base).$(duplicates)$(ext)")
    end
    mkpath(dirname(filepath))
    Downloads.download(url, filepath)

    # get version
    version = mktempdir() do install
        Pkg.PlatformEngines.unpack(filepath, install)
        get_julia_version(install)
    end

    versions = read_versions()
    if haskey(versions, string(version))
        @info "Julia $name (version $version) already available"
        rm(filepath)
    else
        # always use the hash of the downloaded file to force a check during `prepare_julia`
        filehash = hash_file(filepath)

        # rename to include the version
        filename = "julia-$version$ext"
        if ispath(download_dir(filename))
            @warn "Destination file $filename already exists, assuming it matches"
            rm(filepath)
        else
            mv(filepath, download_dir(filename))
        end

        # Update Versions.toml
        version_stanza = """
            ["$version"]
            file = "$filename"
            sha = "$filehash"
            """
        open(extra_versions_file(); append=true) do f
            println(f, version_stanza)
        end
    end

    return version
end

#
# Builds
#

# get a repository handle, cloning or updating the repository along the way
function get_repo(name)
    repo_path = download_dir(name)
    if !isdir(repo_path)
        @debug "Cloning $name to $repo_path..."
        repo = LibGit2.clone("https://github.com/$name", repo_path)
    else
        repo = LibGit2.GitRepo(repo_path)
        LibGit2.fetch(repo)

        # prune to get rid of nonexisting branches (like the pull/PR/merge ones)
        remote = LibGit2.get(LibGit2.GitRemote, repo, "origin")
        fo = LibGit2.FetchOptions(prune=true)
        LibGit2.fetch(remote, String[]; options=fo)
    end

    return repo
end

# look up a refspec in a git repository
function get_repo_commit(repo::LibGit2.GitRepo, spec)
    try
        # maybe it's a remote branch
        remote = LibGit2.get(LibGit2.GitRemote, repo, "origin")
        LibGit2.GitCommit(repo, "refs/remotes/origin/$spec")
    catch err
        isa(err, LibGit2.GitError) || rethrow()
        try
            # maybe it's a version tag
            LibGit2.GitCommit(repo, "refs/tags/v$spec")
        catch err
            isa(err, LibGit2.GitError) || rethrow()
            try
                # maybe it's a remote ref, and we need to fetch it first
                remote = LibGit2.get(LibGit2.GitRemote, repo, "origin")
                LibGit2.fetch(remote, ["+refs/$spec:refs/remotes/origin/$spec"])
                LibGit2.GitCommit(repo, "refs/remotes/origin/$spec")
            catch err
                isa(err, LibGit2.GitError) || rethrow()

                # give up and assume it's a commit or something
                LibGit2.GitCommit(repo, spec)
            end
        end
    end
end

# get the Julia version and short hash corresponding with a refspec in a Julia repo
function get_julia_repoversion(spec, repo_name)
    # get the Julia repo
    repo = get_repo(repo_name)
    commit = get_repo_commit(repo, spec)

    # lookup the version number and commit hash
    tree = LibGit2.peel(LibGit2.GitTree, commit)
    version = VersionNumber(chomp(LibGit2.content(tree["VERSION"])))
    hash = LibGit2.GitHash(commit)
    # FIXME: no way to get the short hash with LibGit2? It just uses the length argument.
    #shorthash = LibGit2.GitShortHash(hash, 7)
    shorthash = LibGit2.GitShortHash(chomp(read(`git -C $(download_dir(repo_name)) rev-parse --short $hash`, String)))
    version = VersionNumber(string(version) * "-" * string(shorthash))
    # NOTE: we append the hash to differentiate commits with identical VERSION files
    #       (we can't only do this for -DEV versions because of backport branches)

    return version, string(hash), string(shorthash)
end

function obtain_julia_build(spec::String="master", repo_name::String="JuliaLang/julia")
    version, hash, shorthash = get_julia_repoversion(spec, repo_name)
    versions = read_versions()
    if haskey(versions, string(version))
        return version
    end

    # try downloading it from the build bots
    url = "https://julialangnightlies.s3.amazonaws.com/bin/linux/x64/$(version.major).$(version.minor)/julia-$(shorthash)-linux64.tar.gz"
    try
        filename = basename(url)
        filepath = download_dir(filename)
        if ispath(filepath)
            @warn "Destination file $filename already exists, assuming it matches"
        else
            Downloads.download(url, filepath)
        end
        filehash = hash_file(filepath)

        # Update Versions.toml
        version_stanza = """
            ["$version"]
            file = "$filename"
            sha = "$filehash"
            """
        open(extra_versions_file(); append=true) do f
            println(f, version_stanza)
        end

        return version
    catch ex
        # if this was a download failure, proceed to build Julia ourselves
        isa(ex, Downloads.RequestError) || rethrow()
        bt = catch_backtrace()
        @error "Could not download Julia $spec (version $version), performing a build" exception=(ex,bt)
        perform_julia_build(spec, repo_name)
    end
end

"""
    version = perform_julia_build(spec::String="master";
                                  binarybuilder_args::Vector{String}=String[]
                                  buildflags::Vector{String}=String[])

Check-out and build Julia at git reference `spec` using BinaryBuilder.
Returns the `version` (what other functions use to identify this build).
This version will be added to Versions.toml.
"""
function perform_julia_build(spec::String="master", repo_name::String="JuliaLang/julia";
                             binarybuilder_args::Vector{String}=String[],
                             buildflags::Vector{String}=String[])
    version, hash, shorthash = get_julia_repoversion(spec, repo_name)
    if !isempty(buildflags)
        version = VersionNumber(version.major, version.minor, version.patch,
                                (version.prerelease...,
                                 "build-$(string(Base.hash(buildflags), base=16))"))
    end
    versions = read_versions()
    if haskey(versions, string(version))
        return version
    end

    # Define a Make.user
    bundled = mktempdir()
    open("$bundled/Make.user", "w") do io
        println(io, "JULIA_CPU_TARGET=generic;sandybridge,-xsaveopt,clone_all;haswell,-rdrnd,base(1)")
        for buildflag in buildflags
            println(io, buildflag)
        end
    end

    # Collection of sources required to build julia
    # NOTE: we need to check out Julia ourselves, because BinaryBuilder does not know how to
    #       fetch and checkout remote refs.
    repo = get_repo(repo_name)
    LibGit2.checkout!(repo, hash)
    repo_path = download_dir(repo_name)
    sources = [
        DirectorySource(repo_path),
        DirectorySource(bundled),
    ]

    # Pre-populate the srccache
    cd(repo_path) do
        Base.run(ignorestatus(`make -C deps getall NO_GIT=1`))
    end

    # Bash recipe for building across all platforms
    script = raw"""
        cd $WORKSPACE/srcdir
        mount -t devpts -o newinstance jrunpts /dev/pts
        mount -o bind /dev/pts/ptmx /dev/ptmx

        make -j${nproc}

        # prevent building documentation
        mkdir -p doc/_build/html/en
        touch doc/_build/html/en/index.html

        make install -j${nproc}

        cp LICENSE.md ${prefix}
        contrib/fixup-libgfortran.sh ${prefix}/lib/julia
        contrib/fixup-libstdc++.sh ${prefix}/lib ${prefix}/lib/julia
    """

    # These are the platforms we will build for by default, unless further
    # platforms are passed in on the command line
    platforms = [
        Linux(:x86_64, libc=:glibc)
    ]

    # The products that we will ensure are always built
    products = Product[
        ExecutableProduct("julia", :julia)
    ]

    # Dependencies that must be installed before this package can be built
    dependencies = []

    # Build the tarballs
    product_hashes = cd(joinpath(@__DIR__, "..", "deps")) do
        build_tarballs(binarybuilder_args, "julia", version, sources, script, platforms,
                       products, dependencies, preferred_gcc_version=v"7", skip_audit=true,
                       verbose=isdebug(:binarybuilder))
    end
    filepath, filehash = product_hashes[platforms[1]]
    filename = basename(filepath)
    if endswith(filename, ".tar.gz")
        ext = ".tar.gz"
        base = filename[1:end-7]
    else
        base, ext = splitext(filename)
    end

    # Copy the generated tarball to the downloads folder
    new_filepath = download_dir(filename)
    duplicates = 0
    while ispath(new_filepath)
        duplicates += 1
        new_filepath = download_dir("$(base).$(duplicates)$(ext)")
    end
    mv(filepath, download_dir(filename))

    # Update Versions.toml
    version_stanza = """
        ["$version"]
        file = "$filename"
        sha = "$filehash"
        """
    open(extra_versions_file(); append=true) do f
        println(f, version_stanza)
    end

    # clean-up
    rm(joinpath(dirname(@__DIR__), "deps", "build"); recursive=true, force=true)
    rm(joinpath(dirname(@__DIR__), "deps", "products"); recursive=true, force=true)

    return version
end


text:
"""
objective.jl
"""

############################################################################################
#                              OBJECTIVES                                                  #
############################################################################################

abstract type AbstractObjective end
state_dim(obj::AbstractObjective) = throw(ErrorException("state_dim not implemented"))
control_dim(obj::AbstractObjective) = throw(ErrorException("control_dim not implemented"))

"""$(TYPEDEF)
Objective: stores stage cost(s) and terminal cost functions

Constructors:
```julia
Objective(cost, cost_term, N)
```
"""
struct Objective{Tc} <: AbstractObjective
    cost::Vector{Tc}
    N::Int
end

# constructors
function Objective(cost::Vector{Tc}, N::Int, checks=true) where {
    Tc<:CostFunction}
    if checks
        @assert length(cost) == N
        for k = 1:N-1
            @assert cost[k].terminal == false
        end
        @assert cost[N].terminal == true
    end
    return Objective{Tc}(cost, N)
end

# methods
Base.copy(obj::Objective) = Objective(copy(obj.cost), obj.N)
Base.show(io::IO, obj::Objective) = print(io,"Objective")

@inline control_dim(obj::Objective) = control_dim(obj.cost[1])
@inline state_dim(obj::Objective) = state_dim(obj.cost[1])

@inline cost(obj::Objective, X::AbstractVector, U::AbstractVector, k::Int) = (
    cost(obj.cost[k], X, U, k)
)

@inline cost_derivatives!(E::QuadraticCost, obj::Objective, X::AbstractVector,
                          U::AbstractVector, k::Int) = (
                              cost_derivatives!(E, obj.cost[k], X, U, k)
)

# LQR objective
function LQRObjective(Q::AbstractMatrix, Qf::AbstractMatrix, R::AbstractMatrix,
                      xf::AbstractVector, n::Int, m::Int, N::Int, M, V)
    stage = LQRCost(Q, xf, R, M, V)
    terminal = LQRCost(Qf, xf, R, M, V; use_R=false, terminal=true)
    Tc = typeof(stage)
    cost = Vector{Tc}(undef, N)
    for k = 1:N-1
        cost[k] = stage
    end
    cost[N] = terminal
    return Objective(cost, N)
end


text:
module StokesEquationTests

using Gridap
using ForwardDiff
using LinearAlgebra
using Test
using GridapODEs.ODETools
using GridapODEs.TransientFETools
using Gridap.FESpaces: get_algebraic_operator

# using GridapODEs.ODETools: ThetaMethodLinear
import Gridap: ∇
import GridapODEs.TransientFETools: ∂t

θ = 0.5

u(x,t) = VectorValue(x[1],x[2])*t
u(t::Real) = x -> u(x,t)

p(x,t) = (x[1]-x[2])*t
p(t::Real) = x -> p(x,t)
q(x) = t -> p(x,t)

f(t) = x -> ∂t(u)(t)(x)-Δ(u(t))(x)+ ∇(p(t))(x)
g(t) = x -> (∇⋅u(t))(x)

domain = (0,1,0,1)
partition = (2,2)
model = CartesianDiscreteModel(domain,partition)

order = 2

reffeᵤ = ReferenceFE(lagrangian,VectorValue{2,Float64},order)
V0 = FESpace(
  model,
  reffeᵤ,
  conformity=:H1,
  dirichlet_tags="boundary"
)

reffeₚ = ReferenceFE(lagrangian,Float64,order-1)
Q = TestFESpace(
  model,
  reffeₚ,
  conformity=:H1,
  constraint=:zeromean
)

U = TransientTrialFESpace(V0,u)

P = TrialFESpace(Q)

Ω = Triangulation(model)
degree = 2*order
dΩ = Measure(Ω,degree)

#
a(u,v) = ∫(∇(u)⊙∇(v))dΩ
b((v,q),t) = ∫(v⋅f(t))dΩ + ∫(q*g(t))dΩ
m(ut,v) = ∫(ut⋅v)dΩ

X = TransientMultiFieldFESpace([U,P])
Y = MultiFieldFESpace([V0,Q])

res(t,(u,p),(ut,pt),(v,q)) = a(u,v) + m(ut,v) - ∫((∇⋅v)*p)dΩ + ∫(q*(∇⋅u))dΩ - b((v,q),t)
jac(t,(u,p),(ut,pt),(du,dp),(v,q)) = a(du,v) - ∫((∇⋅v)*dp)dΩ + ∫(q*(∇⋅du))dΩ
jac_t(t,(u,p),(ut,pt),(dut,dpt),(v,q)) = m(dut,v)

b((v,q)) = b((v,q),0.0)

mat((du1,du2),(v1,v2)) = a(du1,v1)+a(du2,v2)

U0 = U(0.0)
P0 = P(0.0)
X0 = X(0.0)
uh0 = interpolate_everywhere(u(0.0),U0)
ph0 = interpolate_everywhere(p(0.0),P0)
xh0 = interpolate_everywhere([uh0,ph0],X0)

op = TransientFEOperator(res,jac,jac_t,X,Y)

t0 = 0.0
tF = 1.0
dt = 0.1

ls = LUSolver()
odes = ThetaMethod(ls,dt,θ)
solver = TransientFESolver(odes)

sol_t = solve(solver,op,xh0,t0,tF)

l2(w) = w⋅w

tol = 1.0e-6
_t_n = t0

result = Base.iterate(sol_t)

for (xh_tn, tn) in sol_t
  global _t_n
  _t_n += dt
  uh_tn = xh_tn[1]
  ph_tn = xh_tn[2]
  e = u(tn) - uh_tn
  el2 = sqrt(sum( ∫(l2(e))dΩ ))
  e = p(tn) - ph_tn
  el2 = sqrt(sum( ∫(l2(e))dΩ ))
  @test el2 < tol
end

end #module


text:
module ODEsTests

using Test

@time @testset "ODETools" begin include("ODEsTests/runtests.jl") end

@time @testset "TransientFETools" begin include("TransientFEsTests/runtests.jl") end

# @time @testset "DiffEqsWrappers" begin include("DiffEqsWrappersTests/runtests.jl") end

# include("../bench/runbenchs.jl")

end #module


text:
include("plots.jl")

division(a, b) = a ./ b

function build_assets()
    plot_fibo(100)
    plot_factor(100, 2)
    plot_factor(1000, 2)
    plot_factor(100, 3)
    plot_factor(100, 5)
    plot_factor(10000, 23)
    plot_factor(100, 4)
    plot_factor_except(100, 2, 3)
    plot_factor_only(100, 2, 3)
    plot_factor(100, [2, 4])
    plot_factor(100, [2, 3])
    plot_distance(1000, 2, 4, 1, -)
    plot_distance(1000, 3, 9, 1, -)
    plot_distance(500, 2, 3, 1, -)
    plot_distance(1000, 2, 4, 1, division)
    plot_distance(1000, 3, 9, 1, division)
    plot_distance(500, 3, 9, 2, division)
    plot_distance(500, 2, 3, 1, division)
end

max_stars_repo_name: manuelelucchi/Fibonacci

max_issues_repo_name: manuelelucchi/Fibonacci

max_forks_repo_name: manuelelucchi/Fibonacci


text:
#compute the output of the neural net for input x0
#requires that the functions nnAct() and nnDAct() be defined,
#representing the activation function and its derivative, respectively
function nnForward(nn::neuralNet,x0::Array{Float64},nnAct::Function=nnAct,nnDAct::Function=nnDAct)
    #define inputs
    nn.x[nn.xIdx[1]] = x0;
    #loop over levels
    for l=1:nn.nLevels-1
        Wt = reshape(nn.w[nn.wIdx[l]],nn.nNeurons[l+1],nn.nNeurons[l]); #assumes that we've stored W' (row-major)
        z  = Wt*nn.x[nn.xIdx[l]] + nn.b[nn.bIdx[l]];
        nn.x[nn.xIdx[l+1]] = nnAct(z);
        nn.da[nn.bIdx[l]] = nnDAct(z); #save derivatives for later
    end

    return nn.x[nn.xIdx[nn.nLevels]]; #return NN output in case that's helpful
end

max_stars_repo_name: krometis/SimpleNeuralNet

max_issues_repo_name: krometis/SimpleNeuralNet

max_forks_repo_name: krometis/SimpleNeuralNet