text:
import Test: @test, @testset, @test_throws
using Luna
import HDF5
using Distributed

@testset "Chunking" begin
    function contains_all_unique(chunks, x)
        contains_unique = []
        for xi in x
            c = count([xi in chi for chi in chunks])
            push!(contains_unique, (c==1))
        end
        return all(contains_unique)
    end
    pass = []
    for L = 5:11:896
        for n = 1:24
            x = collect(1:L)
            ch = Scans.chunks(x, n)
            push!(pass, contains_all_unique(ch, x))
        end
    end
    @test all(pass)
end

##
args_execs = Dict(["-l"] => Scans.LocalExec,
                  ["-r", "1:6"] => Scans.RangeExec,
                  ["-b", "2,1"] => Scans.BatchExec,
                  ["-q"] =>  Scans.QueueExec)
                #   ["-q", "-p", "4"] =>  Scans.QueueExec)
@testset "Scanning $arg" for (arg, exec) in pairs(args_execs)
for _ in eachindex(ARGS)
    pop!(ARGS)
end
push!(ARGS, arg...)

v = collect(1:10)
scan = Scan("scantest_cmdlineargs"; var=v)
@test scan.exec isa exec
runscan(scan) do scanidx, vi
    @test vi == v[scanidx]
    sleep(rand()) # avoid two processes finishing at precisely the same time
end
end

##
@testset "command-line args overwrite" begin
    for _ in eachindex(ARGS)
        pop!(ARGS)
    end
    push!(ARGS, "--range", "1:4")
    v = collect(1:10)
    scan = Scan("scantest_cmdlineargs_overwrite", Scans.LocalExec(); var=v)
    runscan(scan) do scanidx, vi
        # command line ARGS above should overwrite Scans.LocalExec passed in above
        @test scan.exec isa Scans.RangeExec
    end
    for _ in eachindex(ARGS)
        pop!(ARGS)
    end
end

##
try
@testset "scansave" begin
scan = Scan("scantest_scansave", Scans.LocalExec())
x = collect(1:8)
y = collect(1:6)
addvariable!(scan, :x, x)
addvariable!(scan, :y, y)

runscan(scan) do scanidx, xi, yi
    out = [xi * yi, (xi)^2, (yi)^2]
    slength = xi * yi
    e = fill(1.0*slength, slength)
    em = fill(1.0*slength, (2, slength))
    stats = Dict("energy" => e, "energym" => em)
    xx, yy = xi, yi
    Output.@scansave(scan, scanidx, Eω=out, stats=stats, keyword=[xx, yy])
end
HDF5.h5open("scantest_scansave_collected.h5", "r") do file
    @test read(file["scanvariables"]["x"]) == x
    @test read(file["scanvariables"]["y"]) == y
    @test read(file["scanorder"]) == ["x", "y"]
    out = read(file["Eω"])
    stats = read(file["stats"])
    pass = true
    for (ix, xi) in enumerate(x)
        for (iy, yi) in enumerate(y)
            pass = pass && (out[:, ix, iy] == [xi * yi, (xi)^2, (yi)^2])
            slength = xi * yi
            pass = pass && (stats["valid_length"][ix, iy] == slength)
            e = fill(1.0*slength, slength)
            em = fill(1.0*slength, (2, slength))
            pass = pass && (stats["energy"][1:slength, ix, iy] == e)
            pass = pass && all(isnan.(stats["energy"][slength+1:end, ix, iy]))
            pass = pass && (stats["energym"][:, 1:slength, ix, iy] == em)
            pass = pass && all(isnan.(stats["energym"][:, slength+1:end, ix, iy]))
            pass = pass && (file["keyword"][:, ix, iy] == [xi, yi])
        end
    end
    @test pass
    this = @__FILE__
    code = open(this, "r") do file
                read(file, String)
            end
    @test read(file["script"]) == this * "\n" * code
end
rm("scantest_scansave_collected.h5")
end
catch
rm("scantest_scansave_collected.h5")
end

##
@testset "ScanHDF5Output" begin
var1 = collect(range(1, length=5))
var2 = collect(1:3)
scan = Scan("scantest_hdf5output", Scans.LocalExec(); var1=var1, var2=var2)
files = String[]
runscan(scan) do scanidx, vi1, vi2
    out = Output.@ScanHDF5Output(scan, scanidx, 0, 1, 10)
    @test out["meta"]["scanarrays"]["var1"] == var1
    @test out["meta"]["scanarrays"]["var2"] == var2
    @test out["meta"]["scanvars"]["var1"] == vi1
    @test out["meta"]["scanvars"]["var2"] == vi2
    @test out["meta"]["scanshape"] == [length(var1), length(var2)]
    @test out["meta"]["scanorder"] == ["var1", "var2"]
    for z = range(0, 1; length=10)
        out([1.0, 1.0], z, 0.1, t -> [0.0, 0.0])
    end
    push!(files, out.fpath)
end
z, vvar1, vvar2 = Processing.scanproc() do output
    output["z"], output["meta"]["scanarrays"]["var1"], output["meta"]["scanarrays"]["var2"]
end
@test all(vvar1 .== var1)
@test all(vvar2 .== var2)
@test all(z .== collect(range(0, 1; length=10)))
rm.(files)
end

##
@testset "multi-process queue scan" begin
    ps = addprocs(2)
    @everywhere using Luna
    function worker()
        energies = collect(range(5e-6, 20e-6; length=16))
        scan = Scan("scantest_queue_multiproc", Scans.QueueExec(); energy=energies)
        idcs_run = Int[]
        runscan(scan) do scanidx, energy
            prop_capillary(125e-6, 3, :HeJ, 0.8; λ0=800e-9, τfwhm=10e-15, energy=energy,
                           trange=400e-15, λlims=(200e-9, 4e-6))
            push!(idcs_run, scanidx)
        end
        idcs_run
    end
    r2 = @spawnat ps[1] worker()
    r3 = @spawnat ps[2] worker()
    i2 = fetch(r2)
    i3 = fetch(r3)
    # check that both processes ran something
    @test (length(i2) > 0)
    @test (length(i3) > 0)
    push!(i2, i3...)
    for scanidx in 1:16
        @test count(i2 .== scanidx) == 1 # check that all indices have been run exactly once
    end
    rmprocs(ps)
end

@testset "multi-process queue scan with error" begin
    ps = addprocs(2)
    @everywhere using Luna
    # do it again but with one process giving an error
    scanname = "scantest_queue_multiproc_err"
    function worker_err()
        energies = collect(range(5e-6, 20e-6; length=16))
        scan = Scan(scanname, Scans.QueueExec(); energy=energies)
        idcs_run = Int[]
        runscan(scan) do scanidx, energy
            prop_capillary(125e-6, 3, :HeJ, 0.8; λ0=800e-9, τfwhm=10e-15, energy=energy,
                           trange=400e-15, λlims=(200e-9, 4e-6))
            if scanidx == 16
                error("This exception is expected as part of the test suite")
            end
            push!(idcs_run, scanidx)
        end
        idcs_run
    end
    r2 = @spawnat ps[1] worker_err()
    r3 = @spawnat ps[2] worker_err()
    i2 = fetch(r2)
    i3 = fetch(r3)
    # check that both processes ran something
    @test (length(i2) > 0)
    @test (length(i3) > 0)
    push!(i2, i3...)
    for scanidx in 1:15
        # check that all indices have been run once, except for the one with an error
        @test count(i2 .== scanidx) == 1
    end
    h = string(hash(scanname); base=16)
    qfile = joinpath(Utils.cachedir(), "qfile_$h.h5")
    @test !isfile(qfile) # check that scan completed fully and removed the queue file
    rmprocs(ps)
end

##
@testset "multi-process queue scan via exec" begin
    energies = collect(range(5e-6, 20e-6; length=16))
    scan = Scan("scantest_queue_multiproc_exec", Scans.QueueExec(4); energy=energies)
    td = joinpath(tempdir(), tempname())
    runscan(scan) do scanidx, energy
        println("running on $(myid())")
        prop_capillary(125e-6, 3, :HeJ, 0.8; λ0=800e-9, τfwhm=10e-15, energy=energy,
                       trange=400e-15, λlims=(200e-9, 4e-6),
                       filepath=joinpath(td, makefilename(scan, scanidx)))
        open(joinpath(td, "$(scanidx)_on_$(myid())"), "w") do io
            write(io, "$scanidx ran on $(myid())")
        end
    end
    # should be exactly 2 files per scanidx: output .h5 and "scanidx_on_procid"
    @test length(readdir(td)) == 2length(energies)
    rm(td; recursive=true)
end

# do it again to make sure we can run multiple multi-process scans in one session
@testset "multi-process queue scan via exec--again" begin
    energies = collect(range(5e-6, 20e-6; length=16))
    scan = Scan("scantest_queue_multiproc_exec_again", Scans.QueueExec(4); energy=energies)
    td = joinpath(tempdir(), tempname())
    runscan(scan) do scanidx, energy
        println("running on $(myid())")
        prop_capillary(125e-6, 3, :HeJ, 0.8; λ0=800e-9, τfwhm=10e-15, energy=energy,
                       trange=400e-15, λlims=(200e-9, 4e-6),
                       filepath=joinpath(td, makefilename(scan, scanidx)))
        open(joinpath(td, "$(scanidx)_on_$(myid())"), "w") do io
            write(io, "$scanidx ran on $(myid())")
        end
    end
    # should be exactly 2 files per scanidx: output .h5 and "scanidx_on_procid"
    @test length(readdir(td)) == 2length(energies)
    rm(td; recursive=true)
end

##
@testset "automatic ScanHDF5Output in prop_capillary scan" begin
    energies = collect(range(5e-6, 10e-6; length=4))
    scan = Scan("scantest_autofilename", Scans.LocalExec(); energy=energies)
    mktempdir() do td
        runscan(scan) do scanidx, energy
            prop_capillary(125e-6, 3, :HeJ, 0.8; λ0=800e-9, τfwhm=10e-15, energy=energy,
                        trange=400e-15, λlims=(200e-9, 4e-6), filepath=td, scan, scanidx)
        end
        @test length(readdir(td)) == length(energies)
    end
end

##
@testset "manual filename in ScanHDF5Output" begin
    energies = collect(range(5e-6, 10e-6; length=4))
    scan = Scan("scantest_manualfilename", Scans.LocalExec(); energy=energies)
    mktempdir() do td
        runscan(scan) do scanidx, energy
            prop_capillary(125e-6, 3, :HeJ, 0.8; λ0=800e-9, τfwhm=10e-15, energy=energy,
                        trange=400e-15, λlims=(200e-9, 4e-6), filepath=td, scan, scanidx,
                        filename="newname")
        end
        @test length(readdir(td)) == length(energies)
    @test all(startswith.(readdir(td), "newname"))
    end
end


text:
# ---
# title: 1145. Binary Tree Coloring Game
# id: problem1145
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Tree, Depth-first Search
# link: <https://leetcode.com/problems/binary-tree-coloring-game/description/>
# hidden: true
# ---
#
# Two players play a turn based game on a binary tree.  We are given the `root`
# of this binary tree, and the number of nodes `n` in the tree.  `n` is odd, and
# each node has a distinct value from `1` to `n`.
#
# Initially, the first player names a value `x` with `1 <= x <= n`, and the
# second player names a value `y` with `1 <= y <= n` and `y != x`.  The first
# player colors the node with value `x` red, and the second player colors the
# node with value `y` blue.
#
# Then, the players take turns starting with the first player.  In each turn,
# that player chooses a node of their color (red if player 1, blue if player 2)
# and colors an **uncolored** neighbor of the chosen node (either the left
# child, right child, or parent of the chosen node.)
#
# If (and only if) a player cannot choose such a node in this way, they must
# pass their turn.  If both players pass their turn, the game ends, and the
# winner is the player that colored more nodes.
#
# You are the second player.  If it is possible to choose such a `y` to ensure
# you win the game, return `true`.  If it is not possible, return `false`.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2019/08/01/1480-binary-tree-coloring-
# game.png)
#
#
#
#     Input: root = [1,2,3,4,5,6,7,8,9,10,11], n = 11, x = 3
#     Output: true
#     Explanation: The second player can choose the node with value 2.
#
#
#
#
# **Constraints:**
#
#   * `root` is the root of a binary tree with `n` nodes and distinct node values from `1` to `n`.
#   * `n` is odd.
#   * `1 <= x <= n <= 100`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end

ext: jl


text:
import Base: promote_rule, min, max

"""
Define universal infimum- and supremum values for `Integer` types.
"""
abstract type IntLimit <: Real end

struct Infimum <: IntLimit end
struct Supremum <: IntLimit end

const IntInf = Infimum()
const IntSup = Supremum()

name(::Infimum) = "-Inf"
name(::Supremum) = "Inf"

Base.show(io::IO, x::IntLimit) = print(io, "Int(" * name(x) * ")")

typeinf(::Type{T}) where T<:Union{Integer,Rational} = IntInf
typesup(::Type{T}) where T<:Union{Integer,Rational} = IntSup
typeinf(::Type{T}) where T<:Unsigned = zero(T)
typeinf(::Type{T}) where T<:Real = typemin(T)
typesup(::Type{T}) where T<:Real = typemax(T)

typeinf(::T) where T<:Real = typeinf(T)
typesup(::T) where T<:Real = typesup(T)

promote_rule(::Type{<:IntLimit}, ::Type{BigInt}) = Union{IntLimit,BigInt}
promote_rule(::Type{<:IntLimit}, T::Type{<:Real}) = T

for t in [Int64, Int32, Int16, Int128, Int8, UInt64, UInt32, UInt16, UInt128, UInt8, Float64, Float32, BigFloat]
@eval begin
    (::Type{$t})(::Infimum) = typemin($t)
    (::Type{$t})(::Supremum) = typemax($t)
end
end

(::Type{BigInt})(x::Infimum) = x
(::Type{BigInt})(x::Supremum) = x

# integer arithmetic
import Base:    ==, isless, +, -, *, /, ÷, %

min(i::Infimum, ::BigInt) = i
min(::BigInt, i::Infimum) = i
min(::Supremum, r::BigInt) = r
min(r::BigInt, ::Supremum) = r
max(::Infimum, r::BigInt) = r
max(r::BigInt, ::Infimum) = r
max(s::Supremum, r::BigInt) = s
max(r::BigInt, s::Supremum) = s

==(a::IntLimit, b::IntLimit) = a === b
==(a::IntLimit, b::Real) = isinf(b) && (b < 0) == (a === IntInf)
==(b::Real, a::IntLimit) = isinf(b) && (b < 0) == (a === IntInf)

isless(::Infimum, b::Infimum) = false
isless(::Supremum, b::Supremum) = false
isless(::Supremum, b::Infimum) = false
isless(::Infimum, b::Supremum) = true
isless(a::Infimum, b::T) where T<:AbstractFloat = b != infimum(T)
isless(b::AbstractFloat, a::Infimum) = false
isless(a::Supremum, b::AbstractFloat) = false
isless(b::T, a::Supremum) where T<:AbstractFloat = b != supremum(T)

-(::Infimum) = IntSup
-(::Supremum) = IntInf

+(a::T, ::T) where T<:IntLimit = a
+(a::IntLimit, b::IntLimit) = throw(ArgumentError("inf - inf"))
+(a::IntLimit, ::Real) = a
+(::Real, a::IntLimit) = a
-(a::T, b::T) where T<:IntLimit = a + (-b)
-(a::Infimum, ::Supremum) = IntInf
-(a::Supremum, ::Infimum) = IntSup
-(a::IntLimit, ::Real) = a
-(::Real, a::IntLimit) = -a


text:
function readFamaFrenchRaw(url::ASCIIString)
    # download file, read it in as individual timematrices, delete all
    # downloaded files

    lines = downloadAndRemove(url)
    (startInd, endInd) = findDataBlocks(lines)
    (dataNames, varnames) = getDescriptionAndVariableName(lines, startInd)
    (vals, datesString) = splitOffDates(lines, startInd, endInd)
    valsArr = processValues(vals)
    dates = processDates(datesString)

    nDatasets = length(dates)
    data = Array(Any, nDatasets)
    for ii=1:nDatasets
        tm = Timematr(convert(DataFrame, valsArr[ii]), dates[ii])
        data[ii] = tm
    end

    return (data, dataNames, varnames)
end

function processDates(datesString)
    # transform from string to dates
    nDatasets = length(datesString)

    dates = Array(Any, nDatasets)
    for ii=1:nDatasets
        exampleDate = datesString[ii][1]

        if length(exampleDate) == 4
            dates[ii] = processYearlyDate(datesString[ii])
        elseif length(exampleDate) == 6
            dates[ii] = processMonthlyDate(datesString[ii])
        elseif length(exampleDate) == 8
            dates[ii] = processDailyDate(datesString[ii])
        end
    end
    return dates
end

function processDailyDate(datesArr)
    # 20030308 as string to date
    nObs = length(datesArr)
    dates = Array(Date, nObs)
    for ii=1:nObs
        yyyy = datesArr[ii][1:4]
        mm = datesArr[ii][5:6]
        dd = datesArr[ii][7:8]
        dates[ii] = Date(int(yyyy), int(mm), int(dd))
    end
    return dates
end

function processMonthlyDate(datesArr)
    # 200301 as string to date
    nObs = length(datesArr)
    dates = Array(Date, nObs)
    for ii=1:nObs
        yyyy = datesArr[ii][1:4]
        mm = datesArr[ii][5:6]
        datesBeginMonth = Date(int(yyyy), int(mm), 01)
        ## dd = lastdayofmonth(datesBeginMonth)
        ## dates[ii] = Date(int(yyyy), int(mm), dd)
        dates[ii] = lastdayofmonth(datesBeginMonth)
    end
    return dates
end

function processYearlyDate(datesArr)
    # 2003 as string to date
    nObs = length(datesArr)
    dates = Array(Date, nObs)
    for ii=1:nObs
        dates[ii] = Date(int(datesArr[ii]), 12, 31)
    end
    return dates
end

function processValues(vals)
    # transform from Array(Any, nData) to matrix
    nDatasets = length(vals)
    nVariables = length(vals[1][1])

    valsAsArr = Array(Any, nDatasets)
    for ii=1:nDatasets
        nObs = length(vals[ii])
        valsArr = Array(Float64, nObs, nVariables)

        for jj=1:nObs
            valsArr[jj, :] = vals[ii][jj]
        end
        valsAsArr[ii] = valsArr
    end
    return valsAsArr
end

function splitOffDates(lines, startInd, endInd)
    # for individual dataset, split off dates and values
    nDatasets = length(startInd)

    onlyDataLines = Array(Any, nDatasets)
    for ii=1:nDatasets
        onlyDataLines[ii] = lines[startInd[ii]:endInd[ii]]
    end

    vals = Array(Any, nDatasets)
    datesString = Array(Any, nDatasets)

    for jj=1:nDatasets

        data = onlyDataLines[jj]

        nObs = length(data)

        dats = Array(Union(ASCIIString,UTF8String), nObs)
        remainingPart = Array(Any, nObs)
        for ii=1:nObs
            ## get first digits as dates
            m = match(r"^ *([0-9]{4,})[^0-9]", data[ii])

            remainingPart[ii] = data[ii][(length(m.match)+1):end] |>
                                  strip |>
                                  x -> replace(x, r" +", ",") |> # replace
                                  # multiple whitespaces through
                                  # ","
                                  x -> split(x, ",") |>
                                  float64

            ## get matching dates
            dats[ii] = strip(m.match)
        end

        vals[jj] = remainingPart
        datesString[jj] = dats
    end

    return (vals, datesString)
end

function downloadAndRemove(url::ASCIIString)

    # get filename
    filepath = download(url)
    dirName = dirname(filepath)

    run(`unzip $filepath -d $dirName`)

    # get filename of unzipped file; cut off _TXT from filename
    extInd = basename(url) |>
             x -> searchindex(x, ".") |>
  				 x -> x - 4

    # read in file
    fnameLowCase = string(basename(url)[1:(extInd-1)], ".txt")
    completeFileName = joinpath(dirName, fnameLowCase)
    lines = []
    try
        fnameLowCase = string(basename(url)[1:(extInd-1)], ".txt")
        completeFileName = joinpath(dirName, fnameLowCase)
        lines = open(completeFileName) |>
    	         readlines
    catch
        fnameUpCase = string(basename(url)[1:(extInd-1)], ".TXT")
        completeFileName = joinpath(dirName, fnameUpCase)
        lines = open(completeFileName) |>
    	         readlines
    end

    # remove files from system again
    run(`rm $filepath`)
    run(`rm $completeFileName`)

    return lines
end

function findDataBlocks(lines)
    # split read in file into individual datasets
    nLines = length(lines)

    linesNoCharacters = Array(Bool, nLines)
    for ii=1:nLines
        linesNoCharacters[ii] = ismatch(r"^[ -.0-9]+\r\n$",  lines[ii])
    end

    # find first line of block of data values, and get respective name
    startInd = Array(Any, 0)
    endInd = Array(Any, 0)
    for ii=2:nLines
        # if false followed by true: data block start
        if linesNoCharacters[ii] & !(linesNoCharacters[ii-1])
            push!(startInd, ii)
        end

        # if true followed by false: data block end
        if !(linesNoCharacters[ii]) & linesNoCharacters[ii-1]
            push!(endInd, ii-1)
        end
    end

    # if no copyright at end of file and values right until the bottom
    if linesNoCharacters[end] == true
        push!(endInd, nLines)
    end

    return (startInd, endInd)
end

function getDescriptionAndVariableName(lines, startInd)
    # is two or three lines before start an empty line?
    descriptionLinesAbove = 3

    if ismatch(r"^\r\n$", lines[startInd[1]-3])
        descriptionLinesAbove = 2
    end

    # there are files where only one line of variable names exists
    if ismatch(r"^\r\n$", lines[startInd[1]-2])
        descriptionLinesAbove = 2
    end

    # get dataset names
    nDatasets = length(startInd)
    dataNames = Array(Symbol, nDatasets)
    for ii=1:nDatasets
        dataName = startInd[ii] - descriptionLinesAbove |>
                   x -> lines[x] |>
                   symbol
        dataNames[ii] = dataName
    end

    # get variable names
    if descriptionLinesAbove == 3 # two lines of variable names
        varnames = lines[(startInd[1]-2):(startInd[1]-1)]
    else
        varnames = lines[(startInd[1]-1)]
    end

    return (dataNames, varnames)
end


text:
using Metaheuristics
using Test
import Optim

@testset "Box-constrained" begin

    function run_methods(problem)

        desired_accuracy = 1e-4

        ff, bounds, optimums = Metaheuristics.TestProblems.get_problem(problem)

        # number of function evaluations
        f_calls = 0

        f(x) = begin
            f_calls += 1
            ff(x)
        end

        information = Information(f_optimum = 0.0)
        options = Options(f_tol = desired_accuracy, seed = 1)
        options_2 = Options(f_tol = desired_accuracy, seed = 1, store_convergence=true)
        options_3 = Options(f_tol = desired_accuracy, seed = 1)

        methods = [
                   ABC(options = options, information = information),

                   # differential evolution
                   DE(CR = 0.5, options = options_2, information = information),
                   DE(CR = 0.4, strategy=:best1, options = options_2, information = information),
                   DE(CR = 0.3, strategy=:rand2, options = options_2, information = information),
                   DE(CR = 0.5, strategy=:randToBest1, options = options_2, information = information),
                   DE(CR = 0.3, strategy=:best2, options = options_2, information = information),

                   ECA(options = options, information = information),
                   ECA(adaptive=true, resize_population=true, options = options, information = information),

                   PSO(options = options, information = information),
                   SA(options = options, information = information),
                   WOA(options = options, information = information),
                   MCCGA(options = options_3, information = information),
                  ]

        # CGSA
        for i = 1:10
            push!(methods, CGSA(chaosIndex=i,options = options, information = information))
        end

        for method in methods
            f_calls = 0
            res = optimize(f, bounds, method)
            show(IOBuffer(), res)
            show(IOBuffer(), "text/plain", res.population)
            show(IOBuffer(), "text/html", res.population)
            show(IOBuffer(), res.population[1])
            fitness = minimum( res )
            test_result(fitness, desired_accuracy)

            # also count the mumber of function evaluations
            @test f_calls == nfes(res)
        end
    end

    function run_methods_parallel_eval(problem)

        desired_accuracy = 1e-4

        f_calls = 0
        f(X) = begin
            f_calls += size(X,1)
            sum(X.^2,dims=2)
        end
        bounds = [-ones(5)'; ones(5)']
        options = Options(f_tol = desired_accuracy,
                          f_calls_limit=1000,
                          seed = 1,
                          parallel_evaluation=true
                         )

        methods = [
                   CGSA(options = options),
                   DE(options = options),
                   ECA(options = options),
                   PSO(options = options),
                   WOA(options = options),
                  ]

        for method in methods
            f_calls = 0
            res = optimize(f, bounds, method)
            @test f_calls == nfes(res)
        end
    end

    function test_result(fitness, tol)
        @test ≈(fitness, 0.0, atol=tol)
    end

    for problem in [:sphere]
        run_methods(problem)
        run_methods_parallel_eval(problem)
    end

end


text:
using BioFeatures
using DataFrames
# Make some coordinates
region_lengths = rand(100)*1000
starts = int( 1:1000:100000 )
stops = int(  starts + region_lengths )

# Read these into a BioFeature
regions = DataFrame(start = starts, stop = stops)
F = Features( regions )

seq_ids = keys( F.features )
s = nothing
for seq_id in keys(F.features)
   s = F.features[seq_id]
end

for v = values(s)
  println(v[:start])
end


text:
Qx,x = FlintQQ["x"]
K,a = NumberField(x^6+108)
OK = ring_of_integers(K)
lp = prime_decomposition(OK, 5)
P = lp[1][1]
L = Localization(OK, P)

@testset "OrderLocalization" begin

  @testset "Constructor" begin

    @test parent_type(OrdLocElem{nf_elem}) == OrdLoc{nf_elem}

    lp = prime_decomposition(OK, 3)
    P = lp[1][1]
    L = Localization(OK, P)

    @test elem_type(L) == OrdLocElem{nf_elem}
    @test nf(L) == K
    @test nf(L()) == K

end

#=
  @testset "Canonicalisation" begin
  for i in 1:300
    a = rand(L)
    c = canonical_unit(a)
    @test is_unit(c)
  end
end

    @testset "Valuation" begin

      @test valuation(L(5)) == 1
      @test valuation(L(5//5)) == 0
      @test valuation(L(5//12)) == 1
      @test valuation(L(500)) == 3

      for i in 1:300
        a = rand(L)
        while a == L()
          a = rand(L)
        end
        @test valuation(a) >= 0
      end

    end
    =#

    @testset "Parent object call overloading" begin

      @test parent(L(3)) == L
      @test L(6//3) == L(2)
      @test L(K(5)) == L(5)
      @test L(OK(5)) == L(5)
      @test_throws ErrorException L(24//90)
    end

    @testset "GCDX" begin
      for i in 1:550
        a = rand(L)
        b = rand(L)
        (g,u,v) = gcdx(a,b)
        @test g == u*a + v*b
      end
    end

    @testset "GCD" begin
      for i in 1:550
        a = rand(L)
        b = rand(L)
        g = gcd(a,b)
        @test divides(a,g)[1] && divides(b,g)[1]
      end
    end

    @testset "Exact division" begin

      @test divides(L(15),L(3)) == (true, L(5))
      @test divides(L(15),L(30)) == (true, L(1//2))
      @test divides(L(15),L(25)) == (false, L(0))
      @test_throws ErrorException divexact(L(8),L(10))
      for i in 1:300
        a = rand(L)
        b = rand(L)
        d = divides(a,b)
        if d[1]
          @test a == d[2] * b
        end
      end
    end

    @testset "Inversion" begin

      @test inv(L(8)) == L(1//8)
      @test inv(L(23)) == L(1//23)
      @test_throws ErrorException inv(L(40))
      @test_throws DivideError inv(L())
    end

    @testset "Binary operators" begin

      @test L(19) + L(1) == L(20)
      @test L(19) - L(1) == L(18)
      @test L(3) * L(4//2) == L(6)

      @test L(18//2) + L(2//1) == L(11)
      @test L(18//3) - L(1) == L(5)
      @test L(32) * L(4) == L(128)
    end

    @testset "Basic manipulation" begin

      @test iszero(L()) == true
      @test iszero(L(0)) == true
      @test isone(L(2//3)) == false
      @test is_unit(L(24)) == true

      @test iszero(L(4//3)) == false
      @test isone(L(13//13)) == true
      @test is_unit(L(50//2)) == false
    end
end


text:
"""
`conv_mode!(m::KnetModule, mode::Int)` sets all convolutions
to `mode=mode` in a module.
"""
function conv_mode!(m::KnetModule, mode::Int)
    for m in modules(m)
        if isa(m, Conv)
            m.opt = filter(x->x[1]!==:mode, m.opt)
            push!(m.opt, (:mode, mode))
        end
    end
end

"""
`kaiming(et, h, w, i, o)` Default initialization for conv
"""
kaiming(et, dims...) =
    et(sqrt(2 / prod(dims[1:end-1]))) .* randn(et, dims...)

"""
`Conv`: Performs the operation conv4(w, x; opt...) [.+ b]

# Constructors
    `Conv(r::Int, c::Int, i::Int, o::Int; kwargs)`

# Keywords (Constructor)
     `bias=true` Whether or not to add bias
     `winit=kaiming`  Weight initialization `(eltype, dims...)->array`
     `binit=zeros`    Bias initialization `(eltype, dims...)->array`
     `dtype=Float32`  Element type of weight and bias
     `opt...`         See `conv4` kwargs

# Fields
    `w`: Weight
    `b`: Bias
    `opt`: kwargs of `conv4`, provided in initialization

# Forward execution
    `c(ctx, x)
"""
type Conv <: KnetModule
    w::Param
    b::Union{Param, Void}
    opt
end

function Conv(dims...;
              winit=kaiming,
              binit=zeros,
              bias=true,
              dtype=Float32,
              opt...)
    w = Param(winit(dtype, dims...))
    b = nothing
    if bias
        b = Param(binit(dtype,
                        [1 for i=1:length(dims)-2]...,
                        dims[end], 1))
    end
    return Conv(w, b, opt)
end

function (c::Conv)(ctx, x; o...) #forward(ctx, c::Conv, x)
    o = conv4(val(ctx, c.w), x;
              c.opt..., o...)
    if c.b !== nothing
        o = o .+ val(ctx, c.b)
    end
    return o
end


text:
println("\n\nmake_d_severity: start")

prepend!(LOAD_PATH, ["Project.toml"])
import LibPQ

DISASTER_PASSWORD = read("./etc/disaster-pass", String)

postgres = LibPQ.Connection("host=127.0.0.1 port=5432 dbname=disaster user=disaster password=$DISASTER_PASSWORD")

sql =
"""
drop table if exists public.d_severity;
create table if not exists public.d_severity (
    d_severity_id  integer primary key
);

do \$body\$
declare
    _value  integer;
begin
    for _value in select generate_series(1, 10, 1) loop
        -- raise notice '_value=%', _value;
        insert into public.d_severity (
            d_severity_id
        ) values (
            _value
        ) on conflict do nothing;
    end loop;

end; \$body\$ language plpgsql;
"""
LibPQ.execute(postgres, sql)
LibPQ.close(postgres)

println("make_d_severity: stop")


text:
#!/usr/bin/env julia
using Shell
using Glob
using Suppressor
using Dates

#=
[0-RUN03211300]job-001.sh
[0-ORZ03211301][0-RUN03211300]job-001.sh
=#

runners = filter(x -> strip(x) != "", split(read(ignorestatus(pipeline(
    `ps -U $(ENV["USER"]) -ux`, `grep '[0-9]*:[0-9]*\s*julia\s*.*runner.jl'`)), String), "\n"))

if length(runners) == 0
    error("cannot detect self, the programe must be wrong!")
elseif length(runners) > 1
    println(join(runners, "\n"))
    warn("runner already running, exiting...")
    exit()
end

const jobroot = expanduser("~/jobs")
const jobqueue = joinpath(jobroot, "queue")
const jobdone = joinpath(jobroot, "done")
const jobstop = joinpath(jobroot, "stop")
const jobtrash = joinpath(jobroot, "trash")
const jobresume = joinpath(jobroot, "resume")
const joblog = joinpath(jobroot, "log")
const logfile = joinpath(joblog, "runner.log")
mkpath.([jobroot, jobqueue, jobdone, jobstop, jobtrash, jobresume, joblog])

@suppress Base.println(xs...) = open(f -> (println(f, "[$(now())] ", xs...);
    println(stdout, "[$(now())] ", xs...)), logfile, "a")

println("===== start =====")

const ngpu = parse(Int, read(pipeline(`nvidia-smi -L`, `wc -l`), String))
println("$ngpu GPUs detected.")

RETRY = 0
if length(ARGS) == 1
    VISIBLE_GPU = parse.(Int, split(ARGS[1], ","))
elseif length(ARGS) == 2 && ARGS[1] == "--retry"
    RETRY = parse(Int, ARGS[2])
    VISIBLE_GPU = collect(0:ngpu-1)
elseif length(ARGS) == 4 && ARGS[1] == "--retry" && ARGS[3] == "--gpu"
    RETRY = parse(Int, ARGS[2])
    VISIBLE_GPU = parse.(Int, split(ARGS[4], ","))
else
    VISIBLE_GPU = collect(0:ngpu-1)
end
println("visible gpu: $(VISIBLE_GPU).")

RETRY = min(RETRY, 3)

const CLEAN_TICK = 100
const ticks = CLEAN_TICK * ones(ngpu)
const MAX_BRACKETS = 12

timestamp() = Dates.format(now(), "mmddHHMM")

"""
    an array of gpu status in which `true` means free
"""
function gpustatus()
    stats = ""
    for n in 1:100
        try
            stats = read(`nvidia-smi`, String)
            break
        catch
            print("nvidia-smi err, retry $n...")
            sleep(1 + rand())
        end
    end
    m = match(r"GPU\s*PID\s*Type\s*Process", stats)
    processes = stats[m.offset:end]
    gpus = trues(ngpu)
    for m in eachmatch(r"\|\s+(\d+)\s+\d+\s+C", processes)
        gpus[parse(Int, m.captures[1]) + 1] = false
    end
    return gpus
end

"""
    jobs are `.sh` or `.bk` files in job queue dir.
"""
function nextjob()
    jobs = filter(f->endswith(f, ".sh") || endswith(f, ".sh.bk"), readdir(jobqueue))
    return length(jobs) > 0 ? jobs[1] : nothing
end

"""
    next available gpu.

    gpu is available only if
    * gpu is free, and
    * no job claim that gpu
"""
function nextgpu()
    gpus = gpustatus()
    for i in findall(gpus)
        i - 1 ∉ VISIBLE_GPU && continue
        if(length(glob([Regex("^\\[$(i-1)-RUN\\d+].*\\.sh\$")], jobroot)) == 0)
            return i - 1
        end
    end
    return nothing
end

"""
    stop a job by throwing its .bk file to stop folder.
"""
function stop_job(job)
    println("try stop $job")
    jobstr = replace(replace(job, "[" => "\\["), "]" => "\\]")
    ps = filter(x -> strip(x) != "", split(read(ignorestatus(pipeline(
        `ps -aux`, `grep $jobstr`)), String), "\n"))
    if length(ps) == 1
        m = match(r"\S+\s+(\d+)\s", string(ps))
        try
            println("stopping job $job by cmd $(`kill -9 -$(m.captures[1])`)")
            run(`kill -9 -$(m.captures[1])`)
            println("stopped job $job")
            # move files
            newname = "[STOP$(timestamp())]$job"
            mv(joinpath(jobstop, "$job.bk"), joinpath(jobroot, "$newname.bk"))
            println("mv .bk to $newname.bk")
            if isfile(joinpath(jobroot, job))
                mv(joinpath(jobroot, job), joinpath(jobroot, newname))
                println("mv script to $newname")
            end
            if isfile(joinpath(jobroot, "$job.log"))
                mv(joinpath(jobroot, "$job.log"), joinpath(jobroot, "$newname.log"))
                println("mv log to $newname.log")
            end
        catch err
            println("error when stopping $job, message: $err")
            newname = "[STOPFAIL$(timestamp())]$job"
            mv(joinpath(jobstop, "$job.bk"), joinpath(jobstop, "$newname.bk"))
            println("mv .bk to $newname.bk")
        end
    else
        println("already stopped job $job")
        newname = "[STOPPED$(timestamp())]$job"
        # mv file
        mv(joinpath(jobstop, "$job.bk"), joinpath(jobstop, "$newname.bk"))
        println("mv .bk to $newname.bk")
    end
end

function check_stop()
    stoplist = glob([Regex("^\\[\\d+-RUN\\d+].*\\.sh.bk\$")], jobstop)
    for s in stoplist
        stop_job(basename(s)[1:end-3])
    end
end

function check_resume()
    resumelist = glob(glob"*.sh.bk", jobresume)
    for s in resumelist
        script = strip(read(s, String))
        if length(script) > 2 && script[1] == '"' && script[end] == '"'
            script = strip(script[2:end-1])
        elseif script == ""
            script = "exit 1\necho Empty script!"
        end
        if !occursin("--restore", script)
            script *= " --restore"
        end
        jobname = basename(s)[1:end-3]
        f = open(joinpath(jobqueue, "[RES$(timestamp())]$jobname"), "w")
        println(f, script)
        close(f)
        rm(s)
        println("generate restore script $jobname")
    end
end

"""
    process job.

    * move to job root
    * backup
    * generate running script
        * redirect stdout/stderr
        * if success, mark as done and move self to job done
        * if fail, mark as err
"""
function process_job(job, gpu)
    println("processing $job on gpu $gpu")
    # check job type: normal or bk
    if endswith(job, ".sh.bk")
        jobname = replace(replace(job, r"\[.*\]" => ""), ".bk" => "")
    else
        jobname = job
    end
    jobname = "[$gpu-RUN$(timestamp())]$jobname"
    brackets = collect(eachmatch(r"\[.*?\]", jobname))
    if length(brackets) > MAX_BRACKETS
        cutstart = brackets[MAX_BRACKETS + 1].offset
        cutend = brackets[end].offset + length(brackets[end].match) - 1
        jobname = jobname[1:cutstart-1] * jobname[cutend+1:end]
    end
    jobfile = joinpath(jobroot, jobname)

    # backup script
    mv(joinpath(jobqueue, job), "$jobfile.bk")
    println("backup to $jobname.bk")

    # build running script
    f = open(jobfile, "w")
    script = replace(strip(read("$jobfile.bk", String)), "\$GPU" => gpu)
    if length(script) > 2 && script[1] == '"' && script[end] == '"'
        script = strip(script[2:end-1])
    elseif script == ""
        script = "exit 1\necho Empty script!"
    end
    lines = split(script, "\n", keepempty=false)
    lines[end] = "stdbuf -oL " * lines[end]
    script = join(lines, "\n")
    ismove = length(collect(eachmatch(r"ERR", jobname))) < RETRY

    println(f, """
#!/usr/bin/sh
# redirect output to log file
$script >> '$jobfile.log' 2>&1
# post-process
if [ \$? -eq 0 ]; then
    DATE=\$(date +%m%d%H%M)
    mv '$jobfile' "$jobdone/[DONE\$DATE]$jobname"
    mv '$jobfile.bk' "$jobdone/[DONE\$DATE]$jobname.bk"
    mv '$jobfile.log' "$jobdone/[DONE\$DATE]$jobname.log"
    echo OK
else
    DATE=\$(date +%m%d%H%M)
    mv '$jobfile' "$jobroot/[$gpu-ERR\$DATE]$jobname"
    mv '$jobfile.bk' "$jobroot/[$gpu-ERR\$DATE]$jobname.bk"
    mv '$jobfile.log' "$jobroot/[$gpu-ERR\$DATE]$jobname.log"
    if [ $ismove = true ]; then
        cp "$jobroot/[$gpu-ERR\$DATE]$jobname.bk" "$jobresume"
    fi
    echo FAIL
fi""")
    close(f)
    run(`chmod +x $jobfile`)
    println("generate script $jobname")

    # execute script
    Shell.runfile(jobfile, background=true)
    ticks[gpu+1] = CLEAN_TICK
    println("execute script $jobname")
end

"""
    if a gpu is free for `CLEAN_TICK` seconds and a running job own that gpu,
    then mark this job as unknown.
"""
function check_unknown()
    stats = gpustatus()
    for (i, s) in enumerate(stats)
        if s
            suspects = glob([Regex("^\\[$(i-1)-RUN\\d+].*\\.sh\$")], jobroot)
            if length(suspects) > 0
                ticks[i] -= 1
                if ticks[i] == 0
                    for s in suspects
                        println("unknown detected $s")
                        newname = "[$(i-1)-ORZ$(timestamp())]$(basename(s))"
                        # rename to unknown
                        mv(s, joinpath(jobroot, newname))
                        println("rename script to $newname")
                        # rename log to unknown
                        logname = "$s.log"
                        if isfile(logname)
                            mv(logname, joinpath(jobroot, "$newname.log"))
                            println("rename log to $newname.log")
                        end
                        # rename bk to unknown
                        bkname = "$s.bk"
                        if isfile(bkname)
                            mv(bkname, joinpath(jobroot, "$newname.bk"))
                            println("rename backup to $newname.bk")
                        end
                    end
                else
                    continue
                end
            end
        end
        ticks[i] = CLEAN_TICK
    end
end

while true
    check_stop()
    check_unknown()
    check_resume()

    job = nextjob()

    if job != nothing
        gpu = nextgpu()
        if gpu != nothing
            process_job(job, gpu)
        end
    end

    sleep(1)
end

max_stars_repo_name: innerlee/runner

max_issues_repo_name: innerlee/runner

max_forks_repo_name: innerlee/runner


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

"""
    @decision(model, expr, args..., kw_args...)

Add a decision variable to `model` described by the expression `expr`. If used inside a [`@stage`](@ref) block, the created variable can be used in subsequent stage blocks. [`@recourse`](@ref) should be used to mark decisions in the final stage. See `@variable` for syntax details.

## Examples

```julia
@decision(model, x >= 40)
```

See also [`@recourse`](@ref), [`@parameters`](@ref), [`@uncertain`](@ref), [`@stage`](@ref)
"""
macro decision(def...) @warn "@decision should be used inside a @stage block." end
"""
    @recourse(model, expr, args..., kw_args...)

Add a recourse decision variable to `model` described by the expression `expr`. Replaces [`@decision`](@ref) in the [`@stage`](@ref) block of the final stage, and can only be used there. See `@variable` for syntax details.

## Examples

```julia
@recourse(model, 0 <= y <= 1)
```

See also [`@decision`](@ref), [`@parameters`](@ref), [`@uncertain`](@ref), [`@stage`](@ref)
"""
macro recourse(def...) @warn "@recourse should be used inside the @stage block of the final stage." end
macro _decision(args...)
    args = [args...]
    known = false
    recourse = false
    # Check context of decision definition
    if args[1] isa Symbol && args[1] == :known
        known = true
    elseif args[1] isa Symbol && args[1] == :recourse
        recourse = true
    else
        args[1] == :decision || error("Incorrect usage of @_decision.")
    end
    # Remove context from args
    deleteat!(args, 1)
    # Cache stage and remove from args
    stage = args[1]
    deleteat!(args, 1)
    # Move any set definitions
    set = NoSpecifiedConstraint()
    for (i, arg) in enumerate(args)
        set, arg, found = extract_set(arg)
        if found
            if arg == :()
                deleteat!(args, i)
            else
                args[i] = arg
            end
            break
        end
    end
    # Return expression based on context
    if known
        return esc(:(@variable $((args)...) set = StochasticPrograms.KnownSet($stage)))
    else
        return esc(:(@variable $((args)...) set = StochasticPrograms.DecisionSet($stage, constraint = $set, is_recourse = $recourse)))
    end
end
"""
    @known(def)

Annotate each decision taken in the previous stage. Any [`@decision`](@ref) included in a [`@stochastic_model`](@ref) definition will implicitly add `@known` annotations to subsequent stages.

## Examples

```julia
@known x₁, x₂
```

See also [`@decision`](@ref), [`@parameters`](@ref), [`@uncertain`](@ref), [`@stage`](@ref)
"""
macro known(model_name, args...)
    code = Expr(:block)
    for var in args
        varkey = Meta.quot(var)
        push!(code.args, :($(esc(var)) = $(esc(model_name))[$varkey]))
    end
    return code
end


text:
import RowEchelon

#=
Homework 6
Section 2.5
Problems 4
=#

# Create the matrix and the LU matrices
Ab = [1 -1 2 0 ; 1 -3 1 -5 ; 3 7 5 7]
L = [1 0 0 ; 1 1 0; 3 -5 1]
U = [1 -1 2 ; 0 -2 -1 ; 0 0 -6]

# Row Reduce A
Arref = RowEchelon.rref(Ab)
print("ROW REDUCED SOLUTION: \t")
println(Arref)

# Find Ly = b
# Find Ux = y
# Append b to L
Lb = [1 0 0 0; 1 1 0 -5; 3 -5 1 7]
Lbrref = RowEchelon.rref(Lb)

# Append y to U
Uy = [1 -1 2 Lbrref[1,4] ; 0 -2 -1 Lbrref[2,4]; 0 0 -6 Lbrref[3,4]]
Uyrref = RowEchelon.rref(Uy)
print("LU SOLUTION: \t \t")
println(Uyrref)


text:
module Logger

export OceananigansLogger

using Dates
using Logging
using Crayons

import Logging: shouldlog, min_enabled_level, catch_exceptions, handle_message

const RED    = Crayon(foreground=:red)
const YELLOW = Crayon(foreground=:light_yellow)
const CYAN   = Crayon(foreground=:cyan)
const BLUE   = Crayon(foreground=:blue)

const BOLD      = Crayon(bold=true)
const UNDERLINE = Crayon(underline=true)

struct OceananigansLogger <: Logging.AbstractLogger
              stream :: IO
           min_level :: Logging.LogLevel
      message_limits :: Dict{Any,Int}
    show_info_source :: Bool
end

"""
    OceananigansLogger(stream::IO=stdout, level=Logging.Info; show_info_source=false)

Based on Logging.SimpleLogger, it tries to log all messages in the following format:

    [yyyy/mm/dd HH:MM:SS.sss] log_level message [-@-> source_file:line_number]

where the source of the message between the square brackets is included only if
`show_info_source=true` or if the message is not an info level message.
"""
OceananigansLogger(stream::IO=stdout, level=Logging.Info; show_info_source=false) =
    OceananigansLogger(stream, level, Dict{Any,Int}(), show_info_source)

shouldlog(logger::OceananigansLogger, level, _module, group, id) =
    get(logger.message_limits, id, 1) > 0

min_enabled_level(logger::OceananigansLogger) = logger.min_level

catch_exceptions(logger::OceananigansLogger) = false

function level_to_string(level)
    level == Logging.Error && return "ERROR"
    level == Logging.Warn  && return "WARN "
    level == Logging.Info  && return "INFO "
    level == Logging.Debug && return "DEBUG"
    return string(level)
end

function level_to_crayon(level)
    level == Logging.Error && return RED
    level == Logging.Warn  && return YELLOW
    level == Logging.Info  && return CYAN
    level == Logging.Debug && return BLUE
    return identity
end

function handle_message(logger::OceananigansLogger, level, message, _module, group, id,
                                filepath, line; maxlog = nothing, kwargs...)

    if !isnothing(maxlog) && maxlog isa Int
        remaining = get!(logger.message_limits, id, maxlog)
        logger.message_limits[id] = remaining - 1
        remaining > 0 || return nothing
    end

    buf = IOBuffer()
    iob = IOContext(buf, logger.stream)

    level_name = level_to_string(level)
    crayon = level_to_crayon(level)

    module_name = something(_module, "nothing")
    file_name   = something(filepath, "nothing")
    line_number = something(line, "nothing")
    msg_timestamp = Dates.format(Dates.now(), "[yyyy/mm/dd HH:MM:SS.sss]")

    formatted_message = "$(crayon(msg_timestamp)) $(BOLD(crayon(level_name))) $message"
    if logger.show_info_source || level != Logging.Info
        formatted_message *= " $(BOLD(crayon("-@->"))) $(UNDERLINE("$file_name:$line_number"))"
    end

    println(iob, formatted_message)
    write(logger.stream, take!(buf))

    return nothing
end

end # module


text:
module StopBoostTrigger

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Stop Boost" => Ahorn.EntityPlacement(
        Maple.StopBoostTrigger,
        "rectangle"
    )
)

end


text:
# In this tutorial, we will learn
#
#  - How to approximate vector-valued problems
#  - How to solve problems with complex constitutive laws
#  - How to impose Dirichlet boundary conditions only in selected components
#  - How to impose Dirichlet boundary conditions described by more than one function
#
# ## Problem statement
#
# In this tutorial, we detail how to solve a linear elasticity problem defined on the 3D domain depicted in next figure.
#
# ![](../assets/elasticity/solid.png)
#
# We impose the following boundary conditions. All components of the displacement vector  are constrained to zero on the surface $\Gamma_{\rm G}$, which is marked in green in the figure. On the other hand, the first component of the displacement vector is prescribed to the value $\delta\doteq 5$mm on the surface $\Gamma_{\rm B}$, which is marked in blue. No body or surface forces are included in this example. Formally, the PDE to solve is
#
# ```math
# \left\lbrace
# \begin{aligned}
# -∇\cdot\sigma(u) = 0 \ &\text{in} \ \Omega,\\
# u = 0 \ &\text{on}\ \Gamma_{\rm G},\\
# u_1 = \delta \ &\text{on}\ \Gamma_{\rm B},\\
# \sigma(u)\cdot n = 0 \ &\text{on}\  \Gamma_{\rm N}.\\
# \end{aligned}
# \right.
# ```
#
# The variable $u$ stands for the unknown displacement vector, the vector $n$ is the unit outward normal to the Neumann boundary $\Gamma_{\rm N}\doteq\partial\Omega\setminus\left(\Gamma_{\rm B}\cup\Gamma_{\rm G}\right)$ and $\sigma(u)$ is the stress tensor defined as
# ```math
# \sigma(u) \doteq \lambda\ {\rm tr}(\varepsilon(u)) \ I +2 \mu \  \varepsilon(u),
# ```
# where $I$ is the 2nd order identity tensor, and $\lambda$ and $\mu$ are the *Lamé parameters* of the material. The operator $\varepsilon(u)\doteq\frac{1}{2}\left(\nabla u + (\nabla u)^t \right)$ is the symmetric gradient operator (i.e., the strain tensor). Here, we consider material parameters corresponding to aluminum with Young's modulus $E=70\cdot 10^9$ Pa and Poisson's ratio $\nu=0.33$. From these values, the Lamé parameters are obtained as $\lambda = (E\nu)/((1+\nu)(1-2\nu))$ and $\mu=E/(2(1+\nu))$.
#
#
# ## Numerical scheme
#
# As in previous tutorial, we use a conventional Galerkin FE method with conforming Lagrangian FE spaces. For this formulation, the weak form is: find $u\in U$ such that $ a(u,v) = 0 $ for all $v\in V_0$, where $U$ is the subset of functions in $V\doteq[H^1(\Omega)]^3$ that fulfill the Dirichlet boundary conditions of the problem, whereas $V_0$ are functions in $V$ fulfilling $v=0$ on $\Gamma_{\rm G}$ and $v_1=0$ on $\Gamma_{\rm B}$. The bilinear form of the problem is
# ```math
# a(u,v)\doteq \int_{\Omega} \varepsilon(v) : \sigma(u) \ {\rm d}\Omega.
# ```
#
# The main differences with respect to previous tutorial is that we need to deal with a vector-valued problem, we need to impose different prescribed values on the Dirichlet boundary, and the integrand of the bilinear form $a(\cdot,\cdot)$ is more complex as it involves the symmetric gradient operator and the stress tensor. However, the implementation of this numerical scheme is still done in a user-friendly way since all these features can be easily accounted for with the abstractions in the library.
#
# ## Discrete model
#
# We start by loading the discrete model from a file

using Gridap
model = DiscreteModelFromFile("../models/solid.json")

# In order to inspect it, write the model to vtk

writevtk(model,"model")

# and open the resulting files with Paraview. The boundaries $\Gamma_{\rm B}$ and $\Gamma_{\rm G}$ are identified  with the names `"surface_1"` and `"surface_2"` respectively.  For instance, if you visualize the faces of the model and color them by the field `"surface_2"` (see next figure), you will see that only the faces on $\Gamma_{\rm G}$ have a value different from zero.
#
# ![](../assets/elasticity/solid-surf2.png)
#
# ## Vector-valued FE space
#
# The next step is the construction of the FE space. Here, we need to build a vector-valued FE space, which is done as follows:

order = 1

reffe = ReferenceFE(lagrangian,VectorValue{3,Float64},order)
V0 = TestFESpace(model,reffe;
  conformity=:H1,
  dirichlet_tags=["surface_1","surface_2"],
  dirichlet_masks=[(true,false,false), (true,true,true)])

# As in previous tutorial, we construct a continuous Lagrangian interpolation of order 1. The vector-valued interpolation is selected via the option `valuetype=VectorValue{3,Float64}`, where we use the type `VectorValue{3,Float64}`, which is the way Gridap represents vectors of three `Float64` components. We mark as Dirichlet the objects identified with the tags `"surface_1"` and `"surface_2"` using the `dirichlet_tags` argument. Finally, we chose which components of the displacement are actually constrained on the Dirichlet boundary via the `dirichlet_masks` argument. Note that we constrain only the first component on the boundary $\Gamma_{\rm B}$ (identified as `"surface_1"`), whereas we constrain all components on $\Gamma_{\rm G}$ (identified as `"surface_2"`).
#
# The construction of the trial space is slightly different in this case. The Dirichlet boundary conditions are described with two different functions, one for boundary $\Gamma_{\rm B}$ and another one for $\Gamma_{\rm G}$. These functions can be defined as

g1(x) = VectorValue(0.005,0.0,0.0)
g2(x) = VectorValue(0.0,0.0,0.0)

# From functions `g1` and `g2`, we define the trial space as follows:

U = TrialFESpace(V0,[g1,g2])

# Note that the functions `g1` and `g2` are passed to the `TrialFESpace` constructor in the same order as the boundary identifiers are passed previously in the `dirichlet_tags` argument of the `TestFESpace` constructor.
#
# ## Constitutive law
#
# Once the FE spaces are defined, the next step is to define the weak form.  In this example, the construction of the weak form requires more work than in previous tutorial since we need to account for the constitutive law that relates strain and stress.  The symmetric gradient operator is represented by the function `ε` provided by Gridap (also available as `symmetric_gradient`). However, function `σ` representing the stress tensor is not predefined in the library and it has to be defined ad-hoc by the user, namely

const E = 70.0e9
const ν = 0.33
const λ = (E*ν)/((1+ν)*(1-2*ν))
const μ = E/(2*(1+ν))
σ(ε) = λ*tr(ε)*one(ε) + 2*μ*ε

# Function `σ` takes a strain tensor `ε`(one can interpret this strain as the strain at an arbitrary integration point) and computes the associated stress tensor using the Lamé operator.  Note that the implementation of function `σ` is very close to its mathematical definition.
#
#  ## Weak form
#
#  As seen in previous tutorials, in order to define the week form we need to build the integration mesh and the corresponding measure

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

#  From these objects and the constitutive law previously defined, we can write the weak form as follows

a(u,v) = ∫( ε(v) ⊙ (σ∘ε(u)) )*dΩ
l(v) = 0

# Note that we have composed function `σ` with the strain field `ε(u)` in order to compute the stress field associated with the trial function `u`. The linear form is simply `l(v) = 0` since there are not external forces in this example.
#
# ## Solution of the FE problem
#
# The remaining steps for solving the FE problem are essentially the same as in previous tutorial.

op = AffineFEOperator(a,l,U,V0)
uh = solve(op)

# Note that we do not have explicitly constructed a `LinearFESolver` in order to solve the FE problem. If a `LinearFESolver` is not passed to the `solve` function, a default solver (LU factorization) is created and used internally.
#
# Finally, we write the results to a file. Note that we also include the strain and stress tensors into the results file.

writevtk(Ω,"results",cellfields=["uh"=>uh,"epsi"=>ε(uh),"sigma"=>σ∘ε(uh)])

# It can be clearly observed (see next figure) that the surface  $\Gamma_{\rm B}$ is pulled in $x_1$-direction and that the solid deforms accordingly.
#
# ![](../assets/elasticity/disp_ux_40.png)

# ## Multi-material problems
#
# We end this tutorial by extending previous code to deal with multi-material problems. Let us assume that the piece simulated before is now made of 2 different materials (see next figure). In particular, we assume that the volume depicted in dark green is made of aluminum, whereas the volume marked in purple is made of steel.
#
# ![](../models/solid-mat.png)
#
# The two different material volumes are properly identified in the model we have previously loaded. To check this, inspect the model with Paraview (by writing it to vtk format as done before). Note that the volume made of aluminum is identified as `"material_1"`, whereas the volume made of steel is identified as `"material_2"`.
#
# In order to build the constitutive law for the bi-material problem, we need a vector that contains information about the material each cell in the model is composed. This is achieved by these lines

using Gridap.Geometry
labels = get_face_labeling(model)
dimension = 3
tags = get_face_tag(labels,dimension)

# Previous lines generate a vector, namely `tags`, whose length is the number of cells in the model and for each cell contains an integer that identifies the material of the cell.  This is almost what we need. We also need to know which is the integer value associated with each material. E.g., the integer value associated with `"material_1"` (i.e. aluminum) is retrieved as

const alu_tag = get_tag_from_name(labels,"material_1")

# Now, we know that cells whose corresponding value in the `tags` vector is `alu_tag` are made of aluminum, otherwise they are made of steel (since there are only two materials in this example).
#
# At this point, we are ready to define the multi-material constitutive law. First, we define the material parameters for aluminum and steel respectively:

function lame_parameters(E,ν)
  λ = (E*ν)/((1+ν)*(1-2*ν))
  μ = E/(2*(1+ν))
  (λ, μ)
end

const E_alu = 70.0e9
const ν_alu = 0.33
const (λ_alu,μ_alu) = lame_parameters(E_alu,ν_alu)

const E_steel = 200.0e9
const ν_steel = 0.33
const (λ_steel,μ_steel) = lame_parameters(E_steel,ν_steel)

# Then, we define the function containing the constitutive law:

function σ_bimat(ε,tag)
  if tag == alu_tag
    return λ_alu*tr(ε)*one(ε) + 2*μ_alu*ε
  else
    return λ_steel*tr(ε)*one(ε) + 2*μ_steel*ε
  end
end

# Note that in this new version of the constitutive law, we have included a third argument that represents the integer value associated with a certain material. If the value corresponds to the one for aluminum (i.e., `tag == alu_tag`), then, we use the constitutive law for this material, otherwise, we use the law for steel.
#
# Since we have constructed a new constitutive law, we need to re-define the bilinear form of the problem:

a(u,v) = ∫( ε(v) ⊙ (σ_bimat∘(ε(u),tags)) )*dΩ

# In previous line, pay attention in the usage of the new constitutive law `σ_bimat`. Note that we have passed the vector `tags` containing the material identifiers in the last argument of the function`.
#
# At this point, we can build the FE problem again and solve it

op = AffineFEOperator(a,l,U,V0)
uh = solve(op)

# Once the solution is computed, we can store the results in a file for visualization. Note that, we are including the stress tensor in the file (computed with the bi-material law).

writevtk(Ω,"results_bimat",cellfields=
  ["uh"=>uh,"epsi"=>ε(uh),"sigma"=>σ_bimat∘(ε(uh),tags)])

#  Tutorial done!

max_stars_repo_name: eschnett/Tutorials

max_issues_repo_name: eschnett/Tutorials

max_forks_repo_name: eschnett/Tutorials


text:
# Computes μ(mean) and σ(standard deviation) for each kind of component

function getStepData(IT,f)
    getStepData(IT)
end
function getStepData(IT)
    N = IT.obs
    Nel = IT.nelements[STEP]
    μ = Vector{Float64}(Nel)
    σ = Vector{Float64}(Nel)

    for i = 1:Nel
        μ[i] = (N-i)/N
        σ[i] = sqrt(( i*μ[i]^2+(N-i)*(1-μ[i])^2 )/N)
    end
    return σ,μ
end

function getSlopeData(IT,f)
    getSlopeData(IT)
end
function getSlopeData(IT)
    N = IT.obs
    Nel = IT.nelements[SLOPE]
    μ = Vector{Float64}(Nel)
    σ = Vector{Float64}(Nel)

    for i = 1:Nel
        μ[i] = (
        0.5*(N+1)^2 - 0.5*N - i*(N+1) - 0.5*(i+1)^2 + 0.5*i + i*(i+1)
        )/N
        aux =
        (i*μ[i]^2+(N+1)*i^2+(2*(N+1))*i*μ[i]+(N+1)*μ[i]^2
        -i*(N+1)^2+i*(N+1)-μ[i]*(N+1)^2+μ[i]*(N+1)+(1/3)*(N+1)^3
        -(1/2)*(N+1)^2+(1/6)*N-(i+1)*i^2-(2*(i+1))*i*μ[i]
        -(i+1)*μ[i]^2+i*(i+1)^2-i*(i+1)+μ[i]*(i+1)^2
        -μ[i]*(i+1)-(1/3)*(i+1)^3+(1/2)*(i+1)^2-(1/6)*i)/N
        σ[i] = sqrt(aux)
    end
    return σ,μ
end

function getSpikeData(IT,f)
    getSpikeData(IT)
end
function getSpikeData(IT)
    N = IT.obs
    Nel = IT.nelements[SPIKE]
    μ = Vector{Float64}(Nel)
    σ = Vector{Float64}(Nel)

    for i = 1:Nel
        μ[i] = 1/N
        σ[i] = sqrt( ((N-1)*μ[i]^2 + (1-μ[i])^2)/N )
    end
    return σ,μ
end

function getSineData(IT, ω)
    N = IT.obs
    N_frequencies = IT.nelements[4]
    μ = Vector{Float64}(N_frequencies)
    σ = Vector{Float64}(N_frequencies)

    for i = 1:N_frequencies
        sin_ω = sin(ω[i])
        cos_ω = cos(ω[i])
        sin_N1ω = sin((N+1)*ω[i])
        cos_N1ω = cos((N+1)*ω[i])

        if abs(sin_ω)<1e-8
            println("sin ωs")
            println("$(ω[i]) $i")
        end
        if abs(cos_ω)<1e-8
            println("cos ωs")
            println("$(ω[i]) $i")
        end

        μ[i] =  (
        (1/2)*sin_ω*cos_N1ω/(cos_ω-1)-(1/2)*sin_N1ω
        -(1/2)*sin_ω*cos_ω/(cos_ω-1)+(1/2)*sin_ω
        )/N

        σ[i] = sqrt( (

        (1/2+μ[i]^2)*(N+1)+(1/2)*cos_N1ω^2
        -(1/2)*cos_ω*sin_N1ω*cos_N1ω/sin_ω
        -μ[i]*sin_ω*cos_N1ω/(cos_ω-1)+sin_N1ω*μ[i]
        -1/2-μ[i]^2+μ[i]*sin_ω*cos_ω/(cos_ω-1)-μ[i]*sin_ω

      )/N )
    end
    return σ,μ
end

function getCosData(IT, ω)
    N = IT.obs
    N_frequencies = IT.nelements[5]
    μ = Vector{Float64}(N_frequencies)
    σ = Vector{Float64}(N_frequencies)

    for i = 1:N_frequencies
        sin_ω = sin(ω[i])
        cos_ω = cos(ω[i])
        sin_N1ω = sin((N+1)*ω[i])
        cos_N1ω = cos((N+1)*ω[i])

        if abs(sin_ω)<1e-8
            println("sin ωs")
            println("$(ω[i]) $i")
        end
        if abs(cos_ω)<1e-8
            println("cos ωs")
            println("$(ω[i]) $i")
        end

        μ[i] =  (
                 -(1/2)*cos_N1ω-(1/2)*sin_ω*sin_N1ω/(cos_ω-1)
                 +(1/2)*cos_ω+(1/2)*sin_ω^2/(cos_ω-1)
                )/N

                (1/2+μ[i]^2)*(N+1)-(1/2)*cos_N1ω^2
                +(1/2)*cos_ω*sin_N1ω*cos_N1ω/sin_ω
                +cos_N1ω*μ[i]+μ[i]*sin_ω*sin_N1ω/(cos_ω-1)
                -1/2-μ[i]^2-cos_ω*μ[i]-μ[i]*sin_ω^2/(cos_ω-1)

        σ[i] = sqrt( (
                (1/2+μ[i]^2)*(N+1)-(1/2)*cos_N1ω^2
                +(1/2)*cos_ω*sin_N1ω*cos_N1ω/sin_ω
                +cos_N1ω*μ[i]+μ[i]*sin_ω*sin_N1ω/(cos_ω-1)
                -1/2-μ[i]^2-cos_ω*μ[i]-μ[i]*sin_ω^2/(cos_ω-1)
                )/N )
    end
    return σ,μ
end

getData = Vector{Function}(5)
getData[STEP] = getStepData
getData[SPIKE] = getSpikeData
getData[SLOPE] = getSlopeData
getData[SIN] = getSineData
getData[COS] = getCosData

max_stars_repo_name: joaquimg/ComponentFilterCD

max_issues_repo_name: joaquimg/ComponentFilterCD

max_forks_repo_name: joaquimg/ComponentFilterCD


text:
using PyPlot
using JLD

function main(code::AbstractString)
    log_path_record = "$(homedir())/open_loop/output-record-$(code).jld"
    #log_path_record = "$(homedir())/simulations/output-record-$(code).jld"
    d_rec = load(log_path_record)
    L_b = 0.125

    imu_meas    = d_rec["imu_meas"]
    gps_meas    = d_rec["gps_meas"]
    cmd_log     = d_rec["cmd_log"]
    cmd_pwm_log = d_rec["cmd_pwm_log"]
    vel_est     = d_rec["vel_est"]
    pos_info    = d_rec["pos_info"]

    t0 = max(cmd_pwm_log.t[1],vel_est.t[1],imu_meas.t[1])
    t_end = min(cmd_pwm_log.t[end],vel_est.t[end],imu_meas.t[end])

    t = t0+0.1:.02:t_end-0.1
    v = zeros(length(t))
    #v2 = copy(v)
    psiDot = zeros(length(t))
    cmd = zeros(length(t))
    for i=1:length(t)
        #v[i] = vel_est.z[t[i].>vel_est.t,1][end]
        v[i] = pos_info.z[t[i].>pos_info.t,15][end]
        psiDot[i] = imu_meas.z[t[i].>imu_meas.t,3][end]
        cmd[i] = cmd_pwm_log.z[t[i].>cmd_pwm_log.t,2][end]
    end
    v_x = real(sqrt(complex(v.^2-psiDot.^2*L_b^2)))
    v_y = L_b*psiDot

    # idx = v_x.>1.1
    # psiDot = psiDot[idx]
    # v_x = v_x[idx]
    # cmd = cmd[idx]
    delta = atan2(psiDot*0.25,v_x)

    figure(1)
    ax2=subplot(211)
    plot(t-t0,delta,cmd_log.t-t0,cmd_log.z[:,2])
    grid("on")
    xlabel("t [s]")
    legend(["delta_true","delta_input"])
    subplot(212,sharex=ax2)
    plot(cmd_pwm_log.t-t0,cmd_pwm_log.z[:,2])
    grid("on")

    figure(3)
    plot(t-t0,delta,cmd_pwm_log.t-t0,(cmd_pwm_log.z[:,2]-90)/89.0)
    grid("on")
    xlabel("t [s]")
    ylabel("d_f [rad]")
    legend(["delta_true","delta_input"])

    figure(2)
    ax1=subplot(211)
    plot(cmd_pwm_log.t-t0,floor(cmd_pwm_log.z[:,2]))
    grid("on")
    subplot(212,sharex=ax1)
    plot(vel_est.t-t0,vel_est.z,vel_est.t-t0,mean(vel_est.z[:,2:3],2,),"--")
    grid("on")
    legend(["v","v_fl","v_fr","v_bl","v_br","v_mean"])
end

max_stars_repo_name: RoyaFiroozi/barc

max_issues_repo_name: RoyaFiroozi/barc

max_forks_repo_name: RoyaFiroozi/barc


text:
module JSONRPC

using JSON
import Base.parse

export Request, Response, parse_params

mutable struct Request{method,Tparams}
    id::Union{Nothing,Union{String,Int64}}
    params::Tparams

    function Request{method,Tparams}(id, params) where {method,Tparams}
        return new{method,Tparams}(id, params)
    end

    function Request{method,Tparams}(id::Int32, params) where {method,Tparams}
        return new{method,Tparams}(convert(Int64, id), params)
    end
end

mutable struct Error
end

mutable struct Response{Tresult}
    id::Union{String,Int64}
    result::Union{Nothing,Tresult}
    error::Union{Nothing,Error}
end

Response(id, result) = Response(id, result, nothing)

mutable struct Notification{method,Tparams}
    params::Union{Nothing,Tparams}
end

function parse_params end

function parse(::Type{Request}, message_dict::Dict)
    if message_dict["jsonrpc"] != "2.0"
        error("Invalid JSON-RPC version")
    end
    id = get(message_dict, "id", nothing)
    method = Val{Symbol(message_dict["method"])}
    params = message_dict["params"]

    params_instance = parse_params(method, params)

    ret = Request{method,typeof(params_instance)}(id, params_instance)

    return ret
end

function JSON.json(request::Request{method,Tparams}) where {method, Tparams}
    request_dict = Dict()
    request_dict["jsonrpc"] = "2.0"
    request_dict["method"] = string(method.parameters[1])
    if !(request.id isa Nothing)
        request_dict["id"] = request.id
    end
    request_dict["params"] = request.params
    return JSON.json(request_dict)
end

function JSON.json(response::Response{TResult}) where {TResult}
    response_dict = Dict()
    response_dict["jsonrpc"] = "2.0"
    response_dict["id"] = response.id
    if !(response.result isa Nothing)
        response_dict["result"] = response.result
    elseif !(response.error isa Nothing)
        error("Not yet implemented")
    else
        error("Invalid JSON-RPC response object.")
    end
    return JSON.json(response_dict)
end

function JSON.json(response::Notification{method,Tparams}) where {method, Tparams}
    notification_dict = Dict()
    notification_dict["jsonrpc"] = "2.0"
    notification_dict["method"] = string(method.parameters[1])
    notification_dict["params"] = response.params
    return JSON.json(notification_dict)
end

end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg, BinaryBuilderBase

name = "WhiteboxTools"
version = v"2.0.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/jblindsay/whitebox-tools/archive/refs/tags/v$(version).tar.gz", "18705fc948bdb2f96cd816e5a72d36b9cc460aa8c910383d23fdbd61641aab60")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/whitebox-tools-*/

cargo build --release
mkdir -p "${bindir}"
find target/${rust_target}/release/ -maxdepth 1 -type f -executable | xargs -I '{}' mv {} "${bindir}"
install_license LICENSE.txt
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(; experimental=true)
filter!(p -> !Sys.iswindows(p) || arch(p) != "i686", platforms)

# The products that we will ensure are always built
products = [
    ExecutableProduct("quinn_flow_accumulation", :quinn_flow_accumulation),
    ExecutableProduct("qin_flow_accumulation", :qin_flow_accumulation),
    ExecutableProduct("raster_calculator", :raster_calculator),
    ExecutableProduct("rho8_flow_accumulation", :rho8_flow_accumulation),
    ExecutableProduct("edge_contamination", :edge_contamination),
    ExecutableProduct("whitebox_tools", :whitebox_tools),
    ExecutableProduct("split_vector_lines", :split_vector_lines),
    ExecutableProduct("exposure_towards_wind_flux", :exposure_towards_wind_flux),
    ExecutableProduct("conditional_evaluation", :conditional_evaluation),
]

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

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6", compilers = [:rust, :c], lock_microarchitecture = false)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
# MGH problem 35 - Chebyquad function
#
#   Source:
#   J. J. Moré, B. S. Garbow and K. E. Hillstrom
#   Testing Unconstrained Optimization Software
#   ACM Transactions on Mathematical Software, 7(1):17-41, 1981
#
# A. Montoison, Montreal, 05/2018.

export mgh35

"Chebyquad function"
function mgh35(m :: Int = 10, n :: Int = 10)
  if m < n
    warn(": number of function must be ≥ number of variables. Adjusting to m = n")
    m = n
  end

  nls  = Model()
  x0   = collect(1:n)/(n+1)
  @variable(nls, x[i=1:n], start = x0[i])

  function Tsim(x, n)
    if n == 0
      return 1
    elseif n == 1
      return x
    else
      return 2*x*Tsim(x,n-1) - Tsim(x,n-2)
    end
  end

  Ts = Vector{Function}(undef, n)
  Tnames = Vector{Symbol}(undef, n)
  for i = 1:n
    Ts[i] = x -> Tsim(2*x-1, i)
    Tnames[i] = gensym()
    JuMP.register(nls, Tnames[i], 1, Ts[i], autodiff=true)
  end

  I = [i%2 == 0 ? -1/(i^2-1) : 0 for i = 1:n]
  @NLexpression(nls, F[i=1:n], sum($(Tnames[i])(x[j]) for j = 1:n)/n - I[i])

  return MathProgNLSModel(nls, F, name="mgh35")
end


text:
# # Logistic regression

using TransformVariables, LogDensityProblems, DynamicHMC, DynamicHMC.Diagnostics
using MCMCDiagnostics
using Parameters, Statistics, Random, Distributions, StatsFuns
import ForwardDiff              # use for AD

"""
Logistic regression.

For each draw, ``logit(Pr(yᵢ == 1)) ∼ Xᵢ β``. Uses a `β ∼ Normal(0, σ)` prior.

`X` is supposed to include the `1`s for the intercept.
"""
struct LogisticRegression{Ty, TX, Tσ}
    y::Ty
    X::TX
    σ::Tσ
end

function (problem::LogisticRegression)(θ)
    @unpack y, X, σ = problem
    @unpack β = θ
    loglik = sum(logpdf.(Bernoulli.(logistic.(X*β)), y))
    logpri = sum(logpdf.(Ref(Normal(0, σ)), β))
    loglik + logpri
end

# Make up parameters, generate data using random draws.

N = 1000
β = [1.0, 2.0]
X = hcat(ones(N), randn(N))
y = rand.(Bernoulli.(logistic.(X*β)));

# Create a problem, apply a transformation, then use automatic differentiation.

p = LogisticRegression(y, X, 10.0)   # data and (vague) priors
t = as((β = as(Array, length(β)), )) # identity transformation, just to get the dimension
P = TransformedLogDensity(t, p)      # transformed
∇P = ADgradient(:ForwardDiff, P)

# Sample using NUTS, random starting point.

results = mcmc_with_warmup(Random.GLOBAL_RNG, ∇P, 1000);

# Extract the posterior. (Here the transformation was not really necessary).

β_posterior = first.(transform.(t, results.chain));

# Check that we recover the parameters.

mean(β_posterior)

# Quantiles

qs = [0.05, 0.25, 0.5, 0.75, 0.95]
quantile(first.(β_posterior), qs)

quantile(last.(β_posterior), qs)

# Check that mixing is good.

ess = vec(mapslices(effective_sample_size, reduce(hcat, β_posterior); dims = 2))


text:
### Multi-Lane intersection with fixed obstacles

"""
Container type for the environment parameters
"""
@with_kw mutable struct TInterParams
    # geometry
    x_min::Float64 = -20.0
    x_max::Float64 = 20.0
    y_min::Float64 = -20.0
    inter_x::Float64 = 0.
    inter_y::Float64 = 0.
    inter_width::Float64 = 6.0
    inter_height::Float64 = 1.0
    lane_width::Float64 = 3.0
    nlanes_main::Int64 = 2
    nlanes::Int64 = 1

    stop_line::Float64 = 10. # in m along the ego car lane
    # cars
    speed_limit::Float64 = 10.0 # in m/s
    car_rate::Float64 = 0.3 # probability of a car appearing every time steps
end

"""
Generate a T intersection given a geometry specify by TInterParams
"""
function gen_T_roadway(params::TInterParams)
    #=
    ________________________________________________________
        <---- seg 2     main street <----- seg 1
        -----> seg 3                 -----> seg 4
    _________________            ___________________________
                    | seg . seg  |
                    |  5  .  6   |
                    |     .      |
                    |     .      |
                    |     .      |

    =#
    roadway = Roadway()
    # aliases
    x_min = params.x_min
    x_max = params.x_max
    y_min = params.y_min
    inter_x = params.inter_x
    inter_y = params.inter_y
    inter_width = params.inter_width
    inter_height = params.inter_height
    lane_width = params.lane_width
    nlanes_main = params.nlanes_main
    nlanes = params.nlanes

    seg1 = VecSE2(x_max, inter_y+nlanes_main*lane_width - lane_width/2, pi)
    length1 = x_max - inter_x - inter_width
    add_line!(seg1, nlanes_main, length1, roadway)

    seg2 = VecSE2(inter_x - inter_width, inter_y+nlanes_main*lane_width - lane_width/2, pi)
    length2 =  inter_x - inter_width - x_min
    add_line!(seg2, nlanes_main, length2, roadway)

    seg3 = VecSE2(x_min, inter_y-nlanes_main*lane_width + lane_width/2, 0.0)
    length3 = inter_x - x_min - inter_width
    add_line!(seg3, nlanes_main, length3, roadway)

    seg4 = VecSE2(inter_x + inter_width, inter_y-nlanes_main*lane_width + lane_width/2, 0.0)
    length4 = x_max - inter_x - inter_width
    add_line!(seg4, nlanes_main, length4, roadway)

    seg5 = VecSE2(inter_x-nlanes*lane_width + lane_width/2, inter_y - nlanes_main*lane_width - inter_height, -pi/2)
    length5 = inter_y - nlanes_main*lane_width - y_min - inter_height
    add_line!(seg5, nlanes, length5, roadway)

    seg6 = VecSE2(inter_x+nlanes*lane_width - lane_width/2, y_min, pi/2)
    length6 = inter_y - nlanes_main*lane_width - y_min - inter_height
    add_line!(seg6, nlanes, length6, roadway)

    connections = [Connection(1,2), Connection(3,4)]
    for i=1:nlanes_main
        push!(connections, Connection(1, 5, 0, [(i, 1)]))
        push!(connections, Connection(3, 5, 0, [(i, 1)]))
        push!(connections, Connection(6, 2, 0, [(1, i)]))
        push!(connections, Connection(6, 4, 0, [(1, i)]))
    end

    add_junction!(Junction(connections), roadway)
    return roadway
end

function build_obstacles(params::TInterParams)
    return
end

"""
T intersection environment
"""
mutable struct IntersectionEnv <: OccludedEnv
    roadway::Roadway
    obstacles::Vector{ConvexPolygon}
    params::TInterParams
    lane_map::Dict{String, Lane}
end

function IntersectionEnv(params::TInterParams = TInterParams())
    roadway = gen_T_roadway(params)
    lane_map = Dict{String, Lane}()
    obstalces = Vector{ConvexPolygon}[]
    return IntersectionEnv(roadway, Vector{ConvexPolygon}[], params, lane_map)
end

max_stars_repo_name: MaximeBouton/AutomotivePOMDPs

max_issues_repo_name: MaximeBouton/AutomotivePOMDPs

max_forks_repo_name: MaximeBouton/AutomotivePOMDPs


text:
import Base.max
function max(pot::ArrayPotential,variables;MaxOver=true,ReturnState=false,Ind2Sub=false)
    #MAX Maximise a multi-dimensional array over a set of dimenions
    # maxval=max(x,variables)
    # find the values and states that maximize the multidimensional array x
    # over the dimensions in maxover
    #
    if !MaxOver
        variables = setdiff(pot.variables,variables)
    end

    newvars=setdiff(pot.variables,variables)
    max_variables=memberinds(variables,pot.variables)

    maxval,maxind=findmax(pot.content,max_variables) ## PROBABLY SHOULD CHANGE TO THE MAX FUNCTION -- NEED TO ENSURE THIS CAN WORK FOR ANY POTENTIAL

    #maxval=squeeze(maxval,find(memberinds(size(maxval),1)))
    maxval=squeeze(maxval,tuple(find(memberinds(size(maxval),1))...))

    outpot=deepcopy(pot)
    outpot.variables=newvars; outpot.content=maxval

    if !ReturnState
        return outpot
    end

    if !Ind2Sub
        return outpot, maxind
    end

    s=ind2sub(size(pot.content),maxind[:]) # return in more readable form
    maxstate=zeros(Int64,length(s[1]),length(s))
    for i=1:length(pot.variables)
        maxstate[:,i]=s[i]
    end

    return outpot, maxstate

end

##################################################################
# PotArray type functions:

export numstates
function numstates{T<:ArrayPotential}(A::T)

    if isavector(A.content)
        n=prod(size(A.content));
    else
        n=size(A.content);
    end
    return n
end

export standardise
function standardise{T<:ArrayPotential}(A::T)
    # If the ArrayPotential is a vector, make this a column vector, otherwise leave unchanged
    if isavector(A.content)
        A.content=vec(A.content)
        end
    return A
end


text:
##################
# Games creation #
##################

"""
generate_game(payoffm...)

Returns given payoff matrices for given number of players with given number of actions

**Input parameters**
* `payoffm` - payoff matrices
"""
function generate_game(payoffm...)
    #!all(@. collect(payoffm) isa Array{<: Real}) && return @error "All arguments have to be arrays"
    !all(y->y==size.(payoffm)[1], size.(payoffm)) && return @error "All arrays should have same dimension"
    !all(size.(payoffm) .|> length .== length(payoffm)) && return @error "Dimensions of payoff matrices should equal to number of players (input arguments)"
    return Dict("player"*string(i)=>collect(payoffm)[i] for i in 1:length(payoffm))
end

"""
random_2players_game(dist::Distributions.Sampleable,p1_size::Int,p2_size::Int)

Returns random payoff matrices for 2 players with given number of actions

**Input parameters**
* `dist::Distributions.Sampleable` - distribution from which payoffs are sampled
* `p1_size::Int` - number of actions of first player
* `p2_size::Int` - number of actions of second player
"""
function random_2players_game(dist::Distributions.Sampleable,
                        p1_size::Int,
                        p2_size::Int)
    return Dict("player"*string(i)=>rand(dist, p1_size, p2_size) for i in 1:2)
end

"""
random_nplayers_game(dist::Distributions.Sampleable, size::Union{Array{Int,1},Tuple})

Returns random payoffs for n players with given number of actions

**Input parameters**
* `dist::Distributions.Sampleable` - distribution from which payoffs are sampled
* `size::Union{Array{Int,1},Tuple}` - vector or tuple of number of players' actions
"""
function random_nplayers_game(dist::Distributions.Sampleable,
                        size::Union{Array{Int,1},Tuple})
    return Dict("player"*string(i)=>rand(dist, size...) for i in 1:length(size))
end


text:
module MakieLayout

using Random
using PlotUtils
using AbstractPlotting
using Printf
import Showoff

include("geometry_integration.jl")
include("layout_engine.jl")
include("makie_integration.jl")

export LayoutedAxis
export linkxaxes!
export linkyaxes!
export GridLayout
export AxisLayout
export BBox
export outersolve
export shrinkbymargin
export applylayout

end # module


text:
using Wflow

tomlpath = joinpath(@__DIR__, "sediment_config.toml")
config = Wflow.Config(tomlpath)

model = Wflow.initialize_sediment_model(config)
@unpack network = model

model = Wflow.update(model)

@testset "first timestep sediment model (vertical)" begin
    eros = model.vertical

    @test eros.erosov[1] ≈ 0.9f0
    @test model.clock.iteration == 2
    @test mean(eros.leaf_area_index) ≈ 1.7120018886212223f0
    @test eros.dmsand[1] == 200.0f0
    @test eros.dmlagg[1] == 500.0f0
    @test mean(eros.interception) ≈ 0.4767846753916875f0
    @test mean(eros.soilloss) ≈ 0.008596682196335555f0
end

# run the second timestep
model = Wflow.update(model)

@testset "second timestep sediment model (vertical)" begin
    eros = model.vertical

    @test mean(eros.soilloss) ≈ 0.07601393657280235f0
    @test mean(eros.erosclay) ≈ 0.0022388720384961766f0
    @test mean(eros.erossand) ≈ 0.02572046244407882f0
    @test mean(eros.eroslagg) ≈ 0.022126541806118796f0
    @test mean(eros.TCsed) == 0.0
    @test mean(eros.TCsilt) ≈ 1.0988158364353527f6
    @test mean(eros.TCsand) ≈ 1.0987090622888755f6
    @test mean(eros.TCclay) ≈ 1.0992655197016734f6
    @test mean(eros.TCsilt) ≈ 1.0988158364353527f6
end

@testset "second timestep sediment model (lateral)" begin
    lat = model.lateral

    @test mean(lat.land.inlandsed) ≈ 0.07463801685030906f0
    @test mean(lat.land.inlandclay) ≈ 0.0022367786781657497f0
    @test mean(lat.land.inlandsand) ≈ 0.02519222037812127f0
    @test mean(lat.land.olclay) ≈ 0.006443036462118322f0

    @test mean(lat.river.SSconc) ≈ 0.8259993252994058f0
    @test mean(lat.river.inlandclay) ≈ 0.01980468760667709f0
    @test lat.river.h_riv[network.river.order[end]] ≈ 0.006103649735450745f0
    @test lat.river.outclay[5649] ≈ 2.359031898208781f-9
end

Wflow.close_files(model)

benchmark = @benchmark Wflow.run(tomlpath)
trialmin = BenchmarkTools.minimum(benchmark)
println("Sediment Model update (run)")
print_benchmark(trialmin)


text:
using Test
#using Distributions
#import LessOLS: Sample, ols, yhat, r2, desc

# Testing workflow discovered:
# https://discourse.julialang.org/t/mac-resolving-package-versions-takes-a-long-time/16905/13

include("test_sample.jl")
include("test_ols.jl")
include("test_norris.jl")


text:
#InspectDR: Event handling with Gtk
#-------------------------------------------------------------------------------

import Gtk: GConstants.GdkModifierType, Gtk.GConstants.GdkScrollDirection

#==Constants
===============================================================================#
const MODIFIER_ALT = GdkModifierType.GDK_MOD1_MASK #Is this bad? How do I query?
const MODIFIER_CTRL = GdkModifierType.GDK_CONTROL_MASK
const MODIFIER_SHIFT = GdkModifierType.GDK_SHIFT_MASK

#Other modifiers are ignored... ex: numlock (Mod2)
const MODIFIERS_SUPPORTED = MODIFIER_ALT|MODIFIER_CTRL|MODIFIER_SHIFT

#Symbols not currently supported by Gtk package:
const GdkKeySyms_Plus = 0xffab
const GdkKeySyms_Minus = 0xffad

#==Main types
===============================================================================#
struct ISNormal <: InputState; end #Default state

#==Helper functions
===============================================================================#
modifiers_pressed(eventstate, modmask) = (modmask == (MODIFIERS_SUPPORTED & eventstate))

#Returns 0 or strip index:
#TODO: refer to graphbblist only - instead of pwidget???
function hittest(pwidget::PlotWidget, x::Float64, y::Float64)
	infolist = pwidget.plotinfo.strips #WANTCONST
	for i = 1:length(infolist)
		if isinside(infolist[i].graphbb, x, y)
			return i
		end
	end
	return 0
end

#==Wrapper functions
===============================================================================#
raiseevent_plothover(pwidget::PlotWidget, event::Gtk.GdkEventMotion) =
	raiseevent(pwidget.eh_plothover, pwidget, event.x, event.y)

#==Low-level event handlers
===============================================================================#
function handleevent_plothover(pwidget::PlotWidget, event::Gtk.GdkEventMotion, istrip::Int)
	x, y = event.x, event.y

	if istrip > 0
		xf = Transform2D(pwidget.plotinfo, istrip)
		pt = map2axis(xf, Point2D(x, y))
		ixf = InputXfrm2D(pwidget.plotinfo, istrip)
		pwidget.mouseover.pos = axis2read(pt, ixf) #TODO: do we want "read"able coords?
	else
		pwidget.mouseover.pos = nothing
	end

	raiseevent_plothover(pwidget, event) #Inform owner
	return
end

#Update pwidget.mouseover & raise event to owner
function handleevent_plothover(pwidget::PlotWidget, event::Gtk.GdkEventMotion)
	x, y = event.x, event.y
	istrip = hittest(pwidget, x, y)
	pwidget.mouseover.istrip = istrip
	handleevent_plothover(pwidget, event, istrip)
	return
end

#Last line


text:
@everywhere begin
    @views avg(A) = 0.25 .* ( A[1:end-1,1:end-1] .+ A[2:end,1:end-1] .+ A[1:end-1,2:end] .+ A[2:end,2:end] )
    @views avg_x(A) = 0.5 .* ( A[1:end-1,:] .+ A[2:end,:] )
    @views avg_y(A) = 0.5 .* ( A[:,1:end-1] .+ A[:,2:end] )
    @views diff_x(A) = (A[begin + 1:end, :] .- A[1:end - 1, :])
    @views diff_y(A) = (A[:, begin + 1:end] .- A[:, 1:end - 1])
    @views inn(A) = A[2:end-1,2:end-1]
    end # @everywhere

    function glacier_evolution(glacier_list, ensemble=ensemble)

        # Initialize all matrices for the solver
        S, dSdx, dSdy = zeros(Float64,nx,ny),zeros(Float64,nx-1,ny),zeros(Float64,nx,ny-1)
        dSdx_edges, dSdy_edges, ∇S = zeros(Float64,nx-1,ny-2),zeros(Float64,nx-2,ny-1),zeros(Float64,nx-1,ny-1)
        H, D, dH, Fx, Fy = zeros(Float64,nx,ny), zeros(Float64,nx-1,ny-1),zeros(Float64,nx-2,ny-2),zeros(Float64,nx-1,ny-2),zeros(Float64,nx-2,ny-1)
        V, Vx, Vy = zeros(Float64,nx-1,ny-1),zeros(Float64,nx-1,ny-1),zeros(Float64,nx-1,ny-1)

        # Gather simulation parameters
        current_year = 0
        context = ArrayPartition([A], B, S, dSdx, dSdy, D, dSdx_edges, dSdy_edges, ∇S, Fx, Fy, Vx, Vy, V, [current_year])

        function prob_iceflow_func(prob, i, repeat, context, glacier_list) # closure

            H = glacier_gd[i].distributed_thickness.data # initial ice thickness conditions for forward model
            B = glacier_gd[i].topo.data - glacier_gd.distributed_thickness.data # bedrock

            println("Processing glacier #$i ≈ ", mean(temp_series[i]))
            context.x[2] .= B # Bedrock for current glacier

            return remake(prob, p=context)
        end

        prob_func(prob, i, repeat) = prob_iceflow_func(prob, i, repeat, context, temp_series) # closure

        # Perform reference simulation with forward model
        println("Running forward PDE ice flow model...\n")
        iceflow_prob = ODEProblem(iceflow!,H,(0.0,t₁),context)
        ensemble_prob = EnsembleProblem(iceflow_prob, prob_func = prob_func)
        iceflow_sol = solve(ensemble_prob, ROCK4(), ensemble, trajectories = length(temp_series),
                            pmap_batch_size=length(temp_series), reltol=1e-6,
                            progress=true, saveat=1.0, progress_steps = 50)

        # Save only matrices
        idx = 1
        for result in iceflow_sol
            if idx == 1
                H_refs = result.u[end]
            else
                @views H_refs = cat(H_refs, result.u[end], dims=3)
            end
            idx += 1
        end

        return H_refs
    end

    function generate_ref_dataset(temp_series, H₀, ensemble=ensemble)
        # Compute reference dataset in parallel
        H = deepcopy(H₀)

        # Initialize all matrices for the solver
        S, dSdx, dSdy = zeros(Float64,nx,ny),zeros(Float64,nx-1,ny),zeros(Float64,nx,ny-1)
        dSdx_edges, dSdy_edges, ∇S = zeros(Float64,nx-1,ny-2),zeros(Float64,nx-2,ny-1),zeros(Float64,nx-1,ny-1)
        D, dH, Fx, Fy = zeros(Float64,nx-1,ny-1),zeros(Float64,nx-2,ny-2),zeros(Float64,nx-1,ny-2),zeros(Float64,nx-2,ny-1)
        V, Vx, Vy = zeros(Float64,nx-1,ny-1),zeros(Float64,nx-1,ny-1),zeros(Float64,nx-1,ny-1)
        A = 2e-16
        α = 0                       # Weertman-type basal sliding (Weertman, 1964, 1972). 1 -> sliding / 0 -> no sliding
        C = 15e-14                  # Sliding factor, between (0 - 25) [m⁸ N⁻³ a⁻¹]

        # Gather simulation parameters
        current_year = 0
        context = ArrayPartition([A], B, S, dSdx, dSdy, D, copy(temp_series[5]), dSdx_edges, dSdy_edges, ∇S, Fx, Fy, Vx, Vy, V, C, α, [current_year])

        function prob_iceflow_func(prob, i, repeat, context, temp_series) # closure

            println("Processing temp series #$i ≈ ", mean(temp_series[i]))
            context.x[7] .= temp_series[i] # We set the temp_series for the ith trajectory

            return remake(prob, p=context)
        end

        prob_func(prob, i, repeat) = prob_iceflow_func(prob, i, repeat, context, temp_series) # closure

        # Perform reference simulation with forward model
        println("Running forward PDE ice flow model...\n")
        iceflow_prob = ODEProblem(iceflow!,H,(0.0,t₁),context)
        ensemble_prob = EnsembleProblem(iceflow_prob, prob_func = prob_func)

                    iceflow_sol = solve(ensemble_prob, ROCK4(), ensemble, trajectories = n_trajectories,
                            pmap_batch_size=length(temp_series), reltol=1e-6, save_everystep=false,
                            progress=true, saveat=1.0, progress_steps = 50)

        # Distributed progress bars for solving PDEs in paralell
        # @sync begin
        #     @async begin
        #         tasksdone = 0
        #         while tasksdone < n_trajectories
        #             tasksdone += take!(channel)
        #             update!(progress, tasksdone)
        #         end
        #     end
        #     @async begin
        #         iceflow_sol = solve(ensemble_prob, ROCK4(), ensemble, trajectories = n_trajectories,
        #                     pmap_batch_size=length(temp_series), reltol=1e-6, save_everystep=false,
        #                     progress=true, saveat=1.0, progress_steps = 50)
        #     end
        # end

        # Save only matrices
        idx = 1
        H_refs = []
        for result in iceflow_sol
            if idx == 1
                H_refs = result.u[end]
            else
                @views H_refs = cat(H_refs, result.u[end], dims=3)
            end
            idx += 1
        end

        return H_refs
    end

    function train_iceflow_UDE(H₀, UA, θ, train_settings, H_refs, temp_series)
        H = deepcopy(H₀)
        current_year = 0.0
        optimizer = train_settings[1]
        epochs = train_settings[2]
        # Tuple with all the temp series and H_refs
        context = (B, H, current_year, temp_series)
        loss(θ) = loss_iceflow(θ, context, UA, H_refs) # closure

        # Debugging
        # println("Gradients: ", gradient(loss, θ))
        # @infiltrate

        println("Training iceflow UDE...")
        # println("Using solver: ", solver)
        iceflow_trained = DiffEqFlux.sciml_train(loss, θ, optimizer, cb=callback, maxiters = epochs)

        return iceflow_trained
    end

    @everywhere begin

    callback = function (θ,l) # callback function to observe training
        println("Epoch #$current_epoch - Loss H: ", l)

        pred_A = predict_A̅(UA, θ, collect(-20.0:0.0)')
        pred_A = [pred_A...] # flatten
        true_A = A_fake(-20.0:0.0)

        plot(true_A, label="True A")
        plot_epoch = plot!(pred_A, label="Predicted A")
        savefig(plot_epoch,joinpath(root_plots,"training","epoch$current_epoch.png"))
        global current_epoch += 1

        false
    end

    function loss_iceflow(θ, context, UA, H_refs)
        H_preds = predict_iceflow(θ, UA, context)

        # Compute loss function for the full batch
        l_H = 0.0

        for i in 1:length(H_preds)
            H_ref = H_refs[:,:,i]
            H = H_preds[i].u[end]
            l_H += Flux.Losses.mse(H[H .!= 0.0], H_ref[H.!= 0.0]; agg=mean)
        end

        l_H_avg = l_H/length(H_preds)

        return l_H_avg
    end

    function output_func(sol, i)
        put!(channel, 1)
        sol, false
    end

    function predict_iceflow(θ, UA, context, ensemble=ensemble)

        function prob_iceflow_func(prob, i, repeat, context, UA) # closure
            # B, H, current_year, temp_series)
            temp_series = context[4]

            # println("Processing temp series #$i ≈ ", mean(temp_series[i]))
            # We add the ith temperature series
            iceflow_UDE_batch(H, θ, t) = iceflow_NN(H, θ, t, context, temp_series[i], UA) # closure

            return remake(prob, f=iceflow_UDE_batch)
        end

        prob_func(prob, i, repeat) = prob_iceflow_func(prob, i, repeat, context, UA)

        # (B, H, current_year, temp_series)
        H = context[2]
        tspan = (0.0,t₁)

        iceflow_UDE(H, θ, t) = iceflow_NN(H, θ, t, context, temp_series[5], UA) # closure
        iceflow_prob = ODEProblem(iceflow_UDE,H,tspan,θ)
        ensemble_prob = EnsembleProblem(iceflow_prob, prob_func = prob_func)
        # ensemble_prob = EnsembleProblem(iceflow_prob, prob_func = prob_func, output_func = output_func)

        H_pred = solve(ensemble_prob, solver, ensemble, trajectories = length(temp_series),
            pmap_batch_size=batch_size, u0=H, p=θ, reltol=1e-6, save_everystep=false,
            progress=true, progress_steps = 10)

        # Distributed progress bars for solving UDEs in paralell
        # @sync begin
        #     @async begin
        #         tasksdone = 0
        #         println("tasksdone: ", tasksdone)
        #         while tasksdone < n_trajectories
        #             println("here")
        #             tasksdone += take!(channel)
        #             println("there")
        #             update!(progress, tasksdone)
        #             println("out")
        #         end
        #     end
        #     @async begin
        #         println("solving")
        #         H_pred = solve(ensemble_prob, solver, ensemble, trajectories = length(temp_series),
        #                 pmap_batch_size=batch_size, u0=H, p=θ, reltol=1e-6,
        #                 sensealg = InterpolatingAdjoint(autojacvec=ZygoteVJP()), save_everystep=false,
        #                 progress=true, progress_steps = 10)
        #     end
        # end

        return H_pred
    end

    function iceflow!(dH, H, context,t)
        # Unpack parameters
        #A, B, S, dSdx, dSdy, D, temps, dSdx_edges, dSdy_edges, ∇S, Fx, Fy, Vx, Vy, V, C, α, current_year
        current_year = Ref(context.x[18])
        A = Ref(context.x[1])

        # Get current year for MB and ELA
        year = floor(Int, t) + 1
        if year != current_year[] && year <= t₁
            temp = Ref{Float64}(context.x[7][year])
            A[] .= A_fake(temp[])
            current_year[] .= year
        end

        # Compute the Shallow Ice Approximation in a staggered grid
        SIA!(dH, H, context)
    end

    function iceflow_NN(H, θ, t, context, temps, UA)

        year = floor(Int, t) + 1
        if year <= t₁
            temp = temps[year]
        else
            temp = temps[year-1]
        end

        A = predict_A̅(UA, θ, [temp]) # FastChain prediction requires explicit parameters

        # Compute the Shallow Ice Approximation in a staggered grid
        return SIA(H, A, context)
    end

    """
        SIA!(dH, H, context)

    Compute a step of the Shallow Ice Approximation PDE in a forward model
    """
    function SIA!(dH, H, context)
        # Retrieve parameters
        #A, B, S, dSdx, dSdy, D, norm_temps, dSdx_edges, dSdy_edges, ∇S, Fx, Fy, Vx, Vy, V, C, α, current_year, H_ref, H, UA, θ
        A = context.x[1]
        B = context.x[2]
        S = context.x[3]
        dSdx = context.x[4]
        dSdy = context.x[5]
        D = context.x[6]
        dSdx_edges = context.x[8]
        dSdy_edges = context.x[9]
        ∇S = context.x[10]
        Fx = context.x[11]
        Fy = context.x[12]

        # Update glacier surface altimetry
        S .= B .+ H

        # All grid variables computed in a staggered grid
        # Compute surface gradients on edges
        dSdx .= diff_x(S) / Δx
        dSdy .= diff_y(S) / Δy
        ∇S .= (avg_y(dSdx).^2 .+ avg_x(dSdy).^2).^((n - 1)/2)

        Γ = 2 * A * (ρ * g)^n / (n+2) # 1 / m^3 s
        D .= Γ .* avg(H).^(n + 2) .* ∇S

        # Compute flux components
        dSdx_edges .= diff_x(S[:,2:end - 1]) / Δx
        dSdy_edges .= diff_y(S[2:end - 1,:]) / Δy
        Fx .= .-avg_y(D) .* dSdx_edges
        Fy .= .-avg_x(D) .* dSdy_edges

        #  Flux divergence
        inn(dH) .= .-(diff_x(Fx) / Δx .+ diff_y(Fy) / Δy) # MB to be added here
    end

    # Function without mutation for Zygote, with context as an ArrayPartition
    function SIA(H, A, context)
        # Retrieve parameters
        B = context[1]

        # Update glacier surface altimetry
        S = B .+ H

        # All grid variables computed in a staggered grid
        # Compute surface gradients on edges
        dSdx = diff_x(S) / Δx
        dSdy = diff_y(S) / Δy
        ∇S = (avg_y(dSdx).^2 .+ avg_x(dSdy).^2).^((n - 1)/2)

        Γ = 2 * A * (ρ * g)^n / (n+2) # 1 / m^3 s
        D = Γ .* avg(H).^(n + 2) .* ∇S

        # Compute flux components
        dSdx_edges = diff_x(S[:,2:end - 1]) / Δx
        dSdy_edges = diff_y(S[2:end - 1,:]) / Δy
        Fx = .-avg_y(D) .* dSdx_edges
        Fy = .-avg_x(D) .* dSdy_edges

        #  Flux divergence
        @tullio dH[i,j] := -(diff_x(Fx)[pad(i-1,1,1),pad(j-1,1,1)] / Δx + diff_y(Fy)[pad(i-1,1,1),pad(j-1,1,1)] / Δy) # MB to be added here

        return dH
    end

    function A_fake(temp, noise=false)
        A = @. minA + (maxA - minA) * ((temp-minT)/(maxT-minT) )^2
        if noise
            A = A .+ randn(rng_seed(), length(temp)).*4e-17
        end

        return A
    end

    predict_A̅(UA, θ, temp) = UA(temp, θ) .* 1e-16

    function fake_temp_series(t, means=Array{Float64}([0,-2.0,-3.0,-5.0,-10.0,-12.0,-14.0,-15.0,-20.0]))
        temps, norm_temps, norm_temps_flat = [],[],[]
        for mean in means
           push!(temps, mean .+ rand(t).*1e-1) # static
           append!(norm_temps_flat, mean .+ rand(t).*1e-1) # static
        end

        # Normalise temperature series
        norm_temps_flat = Flux.normalise([norm_temps_flat...]) # requires splatting

        # Re-create array of arrays
        for i in 1:t₁:length(norm_temps_flat)
            push!(norm_temps, norm_temps_flat[i:i+(t₁-1)])
        end

        return temps, norm_temps
    end

    end # @everywhere


text:
using CodecZlib: GzipDecompressorStream

using CSV

using JuliaDB

using IndexedTables

function make_vcf_ndsparse(vcf_gz_file_path::String)

    start_time = now()

    vcf_file_path = string(Array(split(vcf_gz_file_path, ".gz"))[1])

    if isfile(vcf_file_path) == false

        io_gz = open(vcf_gz_file_path)

        io = GzipDecompressorStream(io_gz)

        file = CSV.File(io, comment="##", delim='\t', header=1)

        CSV.write(vcf_file_path, delim='\t', file)

        close(io)

    end

    vcf_ndsparse = loadndsparse(
        vcf_file_path,
        delim='\t',
        header_exists=false,
        colnames=[:CHROM, :POS, :ID, :REF, :ALT, :QUAL, :FILTER, :INFO, :FORMAT, :GERM],
        )

    n_rows = length(rows(vcf_ndsparse))

    vcf_ndsparse = vcf_ndsparse[2:n_rows]

    return vcf_ndsparse

    end_time = now()

    println("\nTook $(canonicalize(Dates.CompoundPeriod(end_time - start_time))).\n")

end

export make_vcf_ndsparse