text:
"""
    Life <: NeighborhoodRule

    Life(neighborhood, born=3, survive=(2, 3))

Rule for game-of-life style cellular automata. This is a demonstration of
Cellular Automata more than a seriously optimised game of life rule.

Cells becomes active if it is empty and the number of neightbors is a number in
the `born`, and remains active the cell is active and the number of neightbors is
in `survive`.

## Examples (gleaned from CellularAutomata.jl)

```julia
using DynamicGrids, Distributions

# Use `Binomial` to tweak the density random true values
init = Bool.(rand(Binomial(1, 0.5), 70, 70))
output = REPLOutput(init; tspan=1:100, fps=25, color=:red)

# Morley
sim!(output, Ruleset(Life(born=[3, 6, 8], survive=[2, 4, 5])))

# 2x2
sim!(output, Ruleset(Life(born=[3, 6], survive=[1, 2, 5])))

# Dimoeba
init = rand(Bool, 400, 300)
init[:, 100:200] .= 0
output = REPLOutput(init; tspan=1:100, fps=25, color=:blue, style=Braile())
sim!(output, Life(born=(3, 5, 6, 7, 8),  survive=(5, 6, 7, 8)))

## No death
sim!(output, Life(born=(3,),  survive=(0, 1, 2, 3, 4, 5, 6, 7, 8)))

## 34 life
sim!(output, Life(born=(3, 4), survive=(3, 4)))

# Replicator
init = fill(true, 300,300)
init[:, 100:200] .= false
init[10, :] .= 0
output = REPLOutput(init; tspan=1:100, fps=25, color=:yellow)
sim!(output, Life(born=(1, 3, 5, 7),  survive=(1, 3, 5, 7)))
nothing

# output

```

![REPL Life](https://raw.githubusercontent.com/cesaraustralia/DynamicGrids.jl/media/life.gif)
"""
struct Life{R,W,N,B,S,L} <: NeighborhoodRule{R,W}
    "A Neighborhood, usually Moore(1)"
    neighborhood::N
    "Int, Array, Tuple or Iterable of values that match sum(neighbors) when cell is empty"
    born::B
    "Int, Array, Tuple or Iterable of values that match sum(neighbors) when cell is full"
    survive::S
    lookup::L
end
function Life{R,W}(
    neighborhood::N, born::B, survive::S, lookup_
) where {R,W,N<:Neighborhood{<:Any,<:Any,L},B,S} where L
    lookup = ntuple(i -> (i - 1) in born, L+1), ntuple(i -> (i - 1) in survive, L+1)
    Life{R,W,N,B,S,typeof(lookup)}(neighborhood, born, survive, lookup)
end
function Life{R,W}(neighborhood, born, survive) where {R,W}
    Life{R,W}(neighborhood, born, survive, nothing)
end
function Life{R,W}(;
    neighborhood=Moore(1),
    born=Param(3, bounds=(0, 8)),
    survive=(Param(2, bounds=(0, 8)), Param(3, bounds=(0, 8))),
) where {R,W}
    Life{R,W}(neighborhood, born, survive, nothing)
end

function setwindow(r::Life{R,W,N,B,S,LU}, buffer) where {R,W,N,B,S,LU}
    hood = setwindow(r.neighborhood, buffer)
    Life{R,W,typeof(hood),B,S,LU}(hood, r.born, r.survive, r.lookup)
end

function applyrule(data, rule::Life, state, I)
    rule.lookup[state + 1][sum(neighbors(rule)) + 1]
end


text:
module PostProcess

using DataFrames
using JLD2
using PlotlyJS
using HypothesisTests
using Clustering
using Statistics

function ma(X, k)
    pad = zeros(eltype(X), k ÷ 2)
    return [pad; [mean(X[i:(i + k)]) for i in 1:(length(X) - k)]; pad]
end

function peak_population(X::AbstractVector; smooth=48 * 7 ÷ 2, threshold=0.0)
    # Smooth signal
    Y = X
    Y = ma(X, smooth)
    # Normlize
    MX = maximum(Y)
    mX = minimum(Y)
    Y = (Y .- mean(Y)) ./ std(Y)
    # Find the peaks
    return peaks = let r = Y .≥ threshold
        ranges = findall(isone, abs.(diff(r)))
        if isodd(length(ranges))
            push!(range, length(X))
        end
        map(Iterators.partition(ranges, 2)) do (a, b)
            _, offset = findmax(@view X[a:b])
            a + offset
        end
    end
end

function peak_population(file::AbstractString; kwargs...)
    peaks = jldopen(file) do f
        map(1:1000) do seed
            df = f["$seed"]
            peak_population(df.count_bph; kwargs...)
        end
    end
    peaks = let lens = length.(peaks)
        s = std(lens)
        m = mean(lens)
        a, b = m - 3s, m + 3s
        filter(x -> a ≤ length(x) ≤ b, peaks)
    end
    npeaks = maximum(length.(peaks))
    peaks = reduce(vcat, peaks)

    cl = kmeans(transpose(peaks), npeaks)
    peaks = map(1:npeaks) do k
        peaks[k .== cl.assignments]
    end
    mpeaks = mean.(peaks)
    perm = sortperm(mpeaks)
    return peaks[perm]
end

function bname(fpath)
    return split(fpath, r"[\\/]")[end]
end

"""
    metadata(file::AbstractString)

Get metadata from file
"""
function metadata(file::AbstractString)
    return jldopen(f -> f["metadata"], file, "r")
end

function merge_result(file::AbstractString)
    io = jldopen(file, "r")
    df = DataFrame(; step=0:2880)
    foreach(keys(io)) do k
        if all(isdigit.(collect(k)))
            seed = parse(Int, k)
            df_seed = io[k]
            for name in names(df_seed)
                if name === "step"
                    continue
                end
                df[!, "$(name)_$(string(seed; pad=4))"] = df_seed[!, name]
            end
        end
    end
    columns = ["step"; sort(setdiff(names(df), ["step"]))]
    close(io)
    return select!(df, columns)
end

function PlotlyJS.plot(filepath::AbstractString, attr)
    f = jldopen(filepath, "r")
    traces = map(1:1000) do seed
        key = string(seed)
        df = f[key]
        scatter(; x=df.step, y=df[:, attr], name="Seed $key")
    end
    filebasename = splitpath(filepath)[end]
    title = replace(filebasename, "_" => "\n")
    layout = Layout(; title="$title", xaxis_title="Step", yaxis_title="$attr")
    close(f)
    return plot(traces, layout)
end

function plot_bph(filepath)
    return plot(filepath, :count_bph)
end
function plot_rice(filepath)
    return plot(filepath, :food)
end

function test_rice(filepath::AbstractString, p0)
    f = jldopen(filepath, "r")
    passed = map(1:1000) do seed
        key = string(seed)
        df = f[key]
        df.food[end] < df.food[begin] ÷ 2 - 100
    end
    close(f)
    return t = BinomialTest(count(passed), length(passed), p0)
end

function is_flower_effective(t; alpha=0.05)
    # H0: flower is not effective -> p = p0
    # Ha: flower is effective -> p < p0 -> tail = left
    # good -> reject H0 -> alpha > pvalue
    #a, b = confint(t; tail=:right, level=1-alpha)
    return alpha ≥ pvalue(t; tail=:left)
end

function test_rice(files::AbstractVector, p0; alpha=0.05)
    metadata = map(files) do file
        jldopen(f -> f["metadata"], file)
    end
    df = DataFrame(metadata)
    df.test = map(files) do file
        test_rice(file, p0)
    end
    df = transform(
        df,
        :envmap => ByRow(bname) => :envmap,
        :test => ByRow(t -> t.x) => :k,
        :test => ByRow(t -> t.n) => :n,
        :test => ByRow(t -> is_flower_effective(t; alpha=alpha)) => :accept,
    )
    return df
end

function batch_test_rice(dir::AbstractString, p0; alpha, clean=false)
    files = joinpath.(dir, readdir(dir))
    files = filter(files) do f
        isfile(f) && endswith(f, ".jld2")
    end
    df = test_rice(files, p0; alpha=alpha)
    if clean
        transform!(df, :test => ByRow(is_flower_effective) => :h0)
        select(df, Not("test"))
    else
        df
    end
end

function batch_peak_population(dir::AbstractString; kwargs...)
    files = joinpath.(dir, readdir(dir))
    map(files) do file
        file => peak_population(file; kwargs...)
    end
end

export JLD2,
    PlotlyJS,
    ProgressMeter,
    HypothesisTest,
    plot,
    savefig,
    test_rice,
    peak_population,
    batch_test_rice,
    batch_peak_population

# end module
end

ext: jl

max_stars_repo_name: ndgnuh/RiceBPH

max_issues_repo_name: ndgnuh/RiceBPH

max_forks_repo_name: ndgnuh/RiceBPH


text:
using VegaDatasets
import IteratorInterfaceExtensions, TableTraits, TableShowUtils
using Test

@testset "VegaDatasets" begin

@testset "Core" begin

for filename in readdir(joinpath(@__DIR__,"..","data","data"))
    if splitext(filename)[2]==".json" || splitext(filename)[2]==".csv"
        if !in(splitdir(filename)[2], ["earthquakes.json", "graticule.json", "londonBoroughs.json", "londonTubeLines.json", "miserables.json", "sf-temps.csv", "us-10m.json", "world-110m.json"])
            vd = dataset(filename)

            @test IteratorInterfaceExtensions.isiterable(vd) == true
            @test TableTraits.isiterabletable(vd) == true

            cvd = collect(vd)

            @test isa(cvd, Array)
        end
    end
end

df = dataset("iris")

@test isa(df, VegaDatasets.VegaDataset)

df2 = dataset("airports")

@test isa(df2, VegaDatasets.VegaDataset)

@test_throws ErrorException dataset("doesnotexist")

@test_throws ErrorException dataset("7zip.png")

df3 = dataset("weather.csv")

@test isa(df3, VegaDatasets.VegaDataset)

df4 = dataset("weather.json")

@test isa(df4, VegaDatasets.VegaDataset)

@test sprint(show, dataset("lookup_groups")) == "9x2 Vega dataset\ngroup │ person\n──────┼───────\n1     │ Alan  \n1     │ George\n1     │ Fred  \n2     │ Steve \n2     │ Nick  \n2     │ Will  \n3     │ Cole  \n3     │ Rick  \n3     │ Tom   "

end

@testset "Show" begin

cars = dataset("cars")
@test sprint((stream,data)->show(stream, "text/html", data), cars) == "

Name	Miles_per_Gallon	Cylinders	Displacement	Horsepower	Weight_in_lbs	Acceleration	Year	Origin

"chevrolet chevelle malibu"	18.0	8	307.0	130	3504	12.0	"1970-01-01"	"USA"

"buick skylark 320"	15.0	8	350.0	165	3693	11.5	"1970-01-01"	"USA"

"plymouth satellite"	18.0	8	318.0	150	3436	11.0	"1970-01-01"	"USA"

"amc rebel sst"	16.0	8	304.0	150	3433	12.0	"1970-01-01"	"USA"

"ford torino"	17.0	8	302.0	140	3449	10.5	"1970-01-01"	"USA"

"ford galaxie 500"	15.0	8	429.0	198	4341	10.0	"1970-01-01"	"USA"

"chevrolet impala"	14.0	8	454.0	220	4354	9.0	"1970-01-01"	"USA"

"plymouth fury iii"	14.0	8	440.0	215	4312	8.5	"1970-01-01"	"USA"

"pontiac catalina"	14.0	8	455.0	225	4425	10.0	"1970-01-01"	"USA"

"amc ambassador dpl"	15.0	8	390.0	190	3850	8.5	"1970-01-01"	"USA"

⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮

... with 396 more rows.
"
@test sprint((stream,data)->show(stream, "application/vnd.dataresource+json", data), cars)[1:300] == "{\"schema\":{\"fields\":[{\"name\":\"Name\",\"type\":\"string\"},{\"name\":\"Miles_per_Gallon\",\"type\":\"number\"},{\"name\":\"Cylinders\",\"type\":\"integer\"},{\"name\":\"Displacement\",\"type\":\"number\"},{\"name\":\"Horsepower\",\"type\":\"integer\"},{\"name\":\"Weight_in_lbs\",\"type\":\"integer\"},{\"name\":\"Acceleration\",\"type\":\"number\"},{\"na"
@test showable("text/html", cars) == true
@test showable("application/vnd.dataresource+json", cars) == true

earthquakes = dataset("earthquakes")
graticule = dataset("graticule")
londonBoroughs = dataset("londonBoroughs")
londonTubeLines = dataset("londonTubeLines")
miserables = dataset("miserables")
sf_temps = dataset("sf-temps")
us_10m = dataset("us-10m")
world_110m = dataset("world-110m")
@test sprint(show, earthquakes) == ""
@test sprint(show, graticule) == ""
@test sprint(show, londonBoroughs) == ""
@test sprint(show, londonTubeLines) == ""
@test sprint(show, miserables) == ""
@test sprint(show, sf_temps) == ""
@test sprint(show, us_10m) == ""
@test sprint(show, world_110m) == ""

end

end


text:
using PortAudio

"""
Continuously read from the default audio input and plot an
ASCII level/peak meter
"""
function micmeter(metersize)
    mic = PortAudioStream(1, 0; latency = 0.1)

    signalmax = zero(eltype(mic))
    println("Press Ctrl-C to quit")
    while true
        block = read(mic, 512)
        blockmax = maximum(abs.(block)) # find the maximum value in the block
        signalmax = max(signalmax, blockmax) # keep the maximum value ever
        print("\r") # reset the cursor to the beginning of the line
        printmeter(metersize, blockmax, signalmax)
    end
end

"""
Print an ASCII level meter of the given size. Signal and peak
levels are assumed to be scaled from 0.0-1.0, with peak >= signal
"""
function printmeter(metersize, signal, peak)
    # calculate the positions in terms of characters
    peakpos = clamp(round(Int, peak * metersize), 0, metersize)
    meterchars = clamp(round(Int, signal * metersize), 0, peakpos - 1)
    blankchars = max(0, peakpos - meterchars - 1)

    for position in 1:meterchars
        printstyled(">", color = barcolor(metersize, position))
    end

    print(" "^blankchars)
    printstyled("|", color = barcolor(metersize, peakpos))
    print(" "^(metersize - peakpos))
end

"""
Compute the proper color for a given position in the bar graph. The first
half of the bar should be green, then the remainder is yellow except the final
character, which is red.
"""
function barcolor(metersize, position)
    if position / metersize <= 0.5
        :green
    elseif position == metersize
        :red
    else
        :yellow
    end
end

micmeter(80)

# function micmeter(metersize)
#     #meter size = 80
#     #max size = 1.0
#     #blockmax varaible max
#     mic = PortAudioStream(1, 0; latency = 0.1)
#     signalmax = zero(eltype(mic))
#     data = []
#     while true
#         block = read(mic, 512)
#         blockmax = maximum(abs.(block)) # find the maximum value in the block
#         signalmax = max(signalmax, blockmax) # keep the maximum value ever
#         println(size(data))
#         if size(data) == (16000,)
#             println("Emptied Array data")
#             senddata()
#             empty!(data)
#         end
#         append!(data, blockmax)
#         # println(metersize, " ",blockmax," ", signalmax)
#     end
# end

# micmeter(80)

max_stars_repo_name: ashwani-rathee/Sardor

max_issues_repo_name: ashwani-rathee/Sardor

max_forks_repo_name: ashwani-rathee/Sardor


text:
function Base.collect(enumerable::Enumerable, ::Type{T}) where {T}
    return T(enumerable)
end

function Base.collect(enumerable::Enumerable, f::Function)
    return f(enumerable)
end


text:
export nImgIn, nImgOut, nFeatIn, nFeatOut, nTheta, getOp, initTheta, AbstractConvKernel

abstract type AbstractConvKernel{T} <: AbstractMeganetElement{T} end

## All convKernel types are assumed to have fields nImage (size of the image) and sK (size of the Convolution Kernel)

function nImgIn(this::AbstractConvKernel)
    return [this.nImg[1]; this.nImg[2]; this.sK[3]]
end

function nImgOut(this::AbstractConvKernel)
   return [this.nImg[1]; this.nImg[2]; this.sK[4]]
end

function nFeatIn(this::AbstractConvKernel)
    return prod(nImgIn(this));
end
function nFeatOut(this::AbstractConvKernel)
   return prod(nImgOut(this));
end

function nTheta(this::AbstractConvKernel)
   return prod(this.sK);
end

function getOp(this::AbstractConvKernel{T},theta::Array{T}) where {T <: Number}

    m = prod(nImgOut(this))
    n = prod(nImgIn(this))

    A = LinearOperator{T}(m,n,false,false,
                        v -> Amv(this,theta,v),
                        Nullable{Function}(),
                        w -> ATmv(this,theta,w))
    return A
end

function initTheta(this::AbstractConvKernel{T}) where {T <: Number}

    sd    = T(0.01);
    theta = sd*randn(T,prod(this.sK));
    #id1 = find(theta>2*sd);
    #theta(id1[:]) = randn(numel(id1),1);

    #id2 = find(theta< -2*sd);
    #theta(id2(:)) = randn(numel(id2),1);
    #theta = max(min(2*sd, theta),-2*sd);
    return theta
end


text:
"""
Setup macro-style documentation datastructures.
"""
macro init()
    STRINGS  = esc(:__DOCILE_STRINGS__)
    METADATA = esc(:__DOCILE_METADATA__)
    quote
        if !isdefined(:__DOCILE_STRINGS__)
            const $STRINGS  = ObjectIdDict()
            const $METADATA = ObjectIdDict()
            nothing
        end
    end
end

docs()     = current_module().__DOCILE_STRINGS__
metadata() = current_module().__DOCILE_METADATA__

meta(docstring = ""; kwargs...) =
    (docstring, @compat(Dict{Symbol, Any}(kwargs)))

"""
Assign to docs to `n` and return a dictionary of keyword arguments.
"""
function data(n, docstr; kwargs...)
    dict = @compat(Dict{Symbol, Any}(kwargs))
    if isa(docstr, AbstractString)
        docs()[n] = docstr
    else
        docs()[n] = docstr[1]
        merge!(dict, docstr[2])
    end
    dict
end

"""
Get the symbolic name of an expression.
"""
function nameof(x::Expr)
    isa(x.args[1], Bool)    && return nameof(x.args[2])
    isexpr(x, :bitstype)    && return x.args[2]
    isexpr(x.args[1], :(.)) && return x.args[1]
    nameof(x.args[1])
end

nameof(s::Symbol) = s

"""
Extract the line number and documented object expression from a block.
"""
unblock(x::Expr) = isexpr(x, :block) ?
    (x.args[1].args[1], x.args[2]) :
    error("Invalid '@doc' syntax.")

"""
Extract the docstring and expression to be documented from an `->` expression.
"""
unarrow(x::Expr) = isexpr(x, :(->)) ?
    (x.args[1], x.args[2]) :
    error("Missing '->' syntax.")

"""
Return the category of an object.
"""
lateguess(::Function) = :function
lateguess(::Module)   = :module
lateguess(::Any)      = error("Unknown category for '@doc'ed object.")

"""
Main documentation generation routine.
"""
function doc(expr::Expr, generic = false)

    # Extract docstring, object, and line number from the expression.
    docs, object = unarrow(expr)
    line, object = unblock(object)

    codesource = :(($(line), @__FILE__))

    name     = nameof(object)
    category = Collector.getcategory(object)

    # Check for correct '@doc+' usage.
    (generic && category != :method) && error("'@doc+' can only be applied to methods.")

    # Easier passing of arguments to subsequent function calls.
    packed = (
        docs,
        object,
        category,
        name,
        codesource
        )

    # Branch to category-specific expression generation.
    generic                                   && return  genericdocs(packed)
    category == :symbol                       && return symbolicdocs(packed)
    category == :method                       && return   methoddocs(packed)
    category == :macro                        && return    macrodocs(packed)
    category in (:global, :typealias)         && return   globaldocs(packed)
    category in (:type, :abstract, :bitstype) && return     typedocs(packed)

    error("Cannot use '@doc' on object.")
end

function genericdocs(packed)
    docs, object, category, name, codesource = packed
    d, o, n = map(esc, (docs, object, name))
    quote
        @init
        $(o)
        metadata()[$(n)] = data($(n), $(d);
            category   = $(Expr(:quote, :function)),
            codesource = $(codesource),
            code       = $(Expr(:quote, object))
            )
        $(n)
    end
end
function symbolicdocs(packed)
    docs, object, category, name, codesource = packed
    d, n = map(esc, (docs, name))
    quote
        @init
        metadata()[$(n)] = data($(n), $(d);
            category   = lateguess($(n)),
            codesource = $(codesource)
            )
        $(n)
    end
end

defined(ex::Expr)  = :(isdefined($(ex.args[1]), $(Expr(:quote, ex.args[2].args[1]))))
defined(s::Symbol) = :(isdefined($(Expr(:quote, s))))

function methoddocs(packed)
    docs, object, category, name, codesource = packed
    d, o, n = map(esc, (docs, object, name))
    quote
        @init
        old = $(esc(defined(name))) ? Set{Method}(methods($(n))) : Set{Method}()
        $(o)
        for m in setdiff(Set{Method}(methods($(n))), old)
            metadata()[m] = data(m, $(d);
                category   = $(Expr(:quote, category)),
                codesource = $(codesource),
                code       = $(Expr(:quote, object))
                )
        end
        $(n)
    end
end
function macrodocs(packed)
    docs, object, category, name, codesource = packed
    name    = symbol(string("@", name))
    d, o, n = map(esc, (docs, object, name))
    quote
        @init
        $(o)
        metadata()[$(n)] = data($(n), $(d);
            category   = $(Expr(:quote, category)),
            codesource = $(codesource),
            code       = $(Expr(:quote, object)),
            signature  = $(Expr(:quote, object.args[1]))
            )
        $(n)
    end
end
function globaldocs(packed)
    docs, object, category, name, codesource = packed
    d, o, n = map(esc, (docs, object, name))
    quote
        @init
        $(o)
        metadata()[$(Expr(:quote, name))] = data(
            $(Expr(:quote, name)), $(d);
            category   = $(Expr(:quote, category)),
            codesource = $(codesource)
            )
        $(n)
    end
end
function typedocs(packed)
    docs, object, category, name, codesource = packed
    d, o, n = map(esc, (docs, object, name))
    quote
        @init
        $(o)
        metadata()[$(n)] = data($(n), $(d);
            category   = $(Expr(:quote, category)),
            codesource = $(codesource)
            )
        $(n)
    end
end

"""
Detect '@doc+' syntax in macro call.
"""
doc(s::Symbol, ex::Expr) = s ≡ :(+) ?
    doc(ex, true) :
    error("Invalid '@doc+' syntax.")

doc(other...) = throw(ArgumentError("Invalid '@doc' syntax."))

if VERSION < v"0.4-dev+6619"
"""
Document an object.

    @doc " ... " ->
    f(x) = x

    @doc+ " ... " ->
    g(x) = x

"""
macro doc(ex) isexpr(ex, :call) ? doc(ex.args...) : doc(ex) end
end


text:
## Restoration-Aware TE with multiple restoration options (lottery tickets) per failure scenario
function ARROW(GRB_ENV, edges, capacity, restored, flows, demand, scenarios, T, Tf, disturbance, verbose)
    nedges = size(edges,1)
    nflows = size(flows,1)
    ntunnels = size(T,1)  # all tunnels routing info on IP edge index
    noptions = size(restored[1], 1)
    nscenarios = size(scenarios,1)
    printstyled("\n** solving ARROW LP with $(noptions) actual lottery tickets...\n", color=:yellow)

    if verbose
        println("Arrow scenario number: ", size(restored,1))
        println("Arrow per-scenario options number: ", size(restored[1],1))
    end
    tunnel_num = 0 # tunnel number for each flow
    for x in 1:length(Tf)
        if size(Tf[x],1) > tunnel_num
            tunnel_num = size(Tf[x],1)
        end
    end

    # residual and restorable tunnels
    Tsf = []  # residual tunnels for all failure scenarios, will just repeat for different restoration options
    Taf = []  # restorable tunnels for all failure scenario and each restoration option
    for s in 1:nscenarios
        sft = []
        aft = []
        for r in 1:size(restored[s], 1)  # for every restore options in a scenario
            ssft = []
            aaft = []
            for f in 1:size(Tf,1)
                ft = []
                at = []
                for t in 1:size(Tf[f],1)
                    up = true
                    partial = false
                    if (length(T[Tf[f][t]]) == 0)
                        up = false
                    end
                    for e in T[Tf[f][t]]
                        if scenarios[s][e] == 0
                            if restored[s][r][e] > 0 # beyond binary failures, less than 1 are all failures
                                up = false
                                partial = true
                            else
                                up = false
                                partial = false
                            end
                        end
                    end
                    if up
                        push!(ft, t)
                    elseif partial
                        push!(at, t)
                    end
                end
                push!(ssft,ft)
                push!(aaft,at)
            end
            push!(sft,ssft)
            push!(aft,aaft)
        end
        push!(Tsf,sft)
        push!(Taf,aft)
    end

    # tunnel routing on IP links
    L = zeros(ntunnels, nedges)
    for t in 1:ntunnels
        for e in 1:nedges
            if in(e, T[t])
                L[t,e] = 1
            end
        end
    end

    model = Model(() -> Gurobi.Optimizer(GRB_ENV))
    set_optimizer_attribute(model, "OutputFlag", 0)
    set_optimizer_attribute(model, "Threads", 32)

    model_def_time = @elapsed begin
    ## The LP optimization model
    @variable(model, b[1:nflows] >= 0)
    @variable(model, a[1:nflows,1:tunnel_num] >= 0)  #tunnels bw
    @variable(model, v_restored[1:nscenarios,1:noptions,1:nedges])

    # Equation 1, flow bandwidth constraint
    for f in 1:nflows
        @constraint(model, sum(a[f,t] for t in 1:size(Tf[f],1)) >= b[f])
    end

    # Equation 2, link capacity constraint
    for e in 1:nedges
        @constraint(model, sum(a[f,t] * L[Tf[f][t],e] for f in 1:nflows for t in 1:size(Tf[f],1)) <= capacity[e])
    end

    # Equation 3, flow bandwidth no larger than demand
    for f in 1:nflows
        @constraint(model, b[f] <= demand[f])   #all allocated bandwidths must be less than the demand for that flow
    end

    # Equation 4, residual and restorable tunnels must be able to carry bandwidth for all scenarios
    for f in 1:nflows
        for s in 1:size(scenarios,1)
            for z in 1:size(restored[s], 1)
                @constraint(model, sum(a[f,t] for t in Tsf[s][z][f]) + sum(a[f,t] for t in Taf[s][z][f]) >= b[f])
            end
        end
    end

    # Equation 5, restorable tunnels and restorable bandwidth capacity
    for s in 1:nscenarios
        for e in 1:nedges
            for z in 1:noptions
                @constraint(model, sum(a[f,t] * L[Tf[f][t],e] for f in 1:nflows for t in Taf[s][z][f]) <= restored[s][z][e] + v_restored[s,z,e])
            end
        end
    end

    # Equation 6, bound for slack restorable bandwidth capacity
    for s in 1:nscenarios
        for z in 1:noptions
            @constraint(model, sum(v_restored[s,z,e] for e in 1:nedges) <= disturbance[s])
        end
    end

    @objective(model, Max, sum(b[i] for i in 1:size(b,1)))
    end  # model_def_time end

    optimize!(model)
    solve_runtime = solve_time(model)
    opt_runtime = solve_runtime + model_def_time

    return value.(a), value.(b), objective_value(model), value.(v_restored), solve_runtime, opt_runtime
end

## Restoration-Aware TE with one restoration option per failure scenario
function ARROW_ONE(GRB_ENV, edges, capacity, restored, flows, demand, scenarios, T, Tf, verbose)
    printstyled("\n** solving ARROW LP using one option...\n", color=:yellow)
    # println(restored)

    nedges = size(edges,1)
    nflows = size(flows,1)
    ntunnels = size(T,1)  # all tunnels routing info on IP edge index
    nscenarios = size(scenarios,1)

    if verbose
        println("Arrow scenario number: ", size(restored,1))
        # println("Arrow per-scenario options number: ", size(restored[1],1))
        # println("restored: ", restored)
    end

    tunnel_num = 0 # tunnel number for each flow
    for x in 1:length(Tf)
        if size(Tf[x],1) > tunnel_num
            tunnel_num = size(Tf[x],1)
        end
    end

    # residual and restorable tunnels
    Tsf = []  # residual tunnels for all failure scenarios, will just repeat for different restoration options
    Taf = []  # restorable tunnels for all failure scenario and each restoration option
    for s in 1:nscenarios
        ssft = []
        aaft = []
        for f in 1:size(Tf,1)
            ft = []
            at = []
            for t in 1:size(Tf[f],1)
                up = true
                partial = false
                if (length(T[Tf[f][t]]) == 0)
                    up = false
                end
                for e in T[Tf[f][t]]
                    if scenarios[s][e] == 0
                        if restored[s][e] > 0 # beyond binary failures, less than 1 are all failures
                            up = false
                            partial = true
                        else
                            up = false
                            partial = false
                        end
                    end
                end
                if up
                    push!(ft, t)
                elseif partial
                    push!(at, t)
                end
            end
            push!(ssft,ft)
            push!(aaft,at)
        end
        push!(Tsf,ssft)
        push!(Taf,aaft)
    end

    # tunnel routing on IP links
    L = zeros(ntunnels, nedges)
    for t in 1:ntunnels
        for e in 1:nedges
            if in(e, T[t])
                L[t,e] = 1
            end
        end
    end

    model = Model(() -> Gurobi.Optimizer(GRB_ENV))
    set_optimizer_attribute(model, "OutputFlag", 0)
    set_optimizer_attribute(model, "Threads", 32)

    model_def_time = @elapsed begin
    ## The LP optimization model
    @variable(model, b[1:nflows] >= 0)
    @variable(model, a[1:nflows,1:tunnel_num] >= 0)  #tunnels bw

    # Equation 7, flow bandwidth constraint
    for f in 1:nflows
        @constraint(model, sum(a[f,t] for t in 1:size(Tf[f],1)) >= b[f])
    end

    # Equation 8, link capacity constraint
    for e in 1:nedges
        @constraint(model, sum(a[f,t] * L[Tf[f][t],e] for f in 1:nflows for t in 1:size(Tf[f],1)) <= capacity[e])
    end

    # Equation 9, flow bandwidth no larger than demand
    for f in 1:nflows
        @constraint(model, b[f] <= demand[f])   #all allocated bandwidths must be less than the demand for that flow
    end

    # Equation 10, residual and restorable tunnels must be able to carry bandwidth for all scenarios
    for f in 1:nflows
        for s in 1:size(scenarios,1)
            # for z in 1:size(restored[s], 1)
            @constraint(model, sum(a[f,t] for t in Tsf[s][f]) + sum(a[f,t] for t in Taf[s][f]) >= b[f])
            # end
        end
    end

    # Equation 11, restorable tunnels and restorable bandwidth capacity
    for s in 1:nscenarios
        for e in 1:nedges
            @constraint(model, sum(a[f,t] * L[Tf[f][t],e] for f in 1:nflows for t in Taf[s][f]) <= restored[s][e])
        end
    end

    @objective(model, Max, sum(b[i] for i in 1:size(b,1)))
    end  # model_def_time end

    optimize!(model)
    solve_runtime = solve_time(model)
    opt_runtime = solve_runtime + model_def_time

    return value.(a), value.(b), objective_value(model), solve_runtime, opt_runtime
end

## Restoration-Aware TE with multiple restoration options (lottery tickets) and lottery ticket selection per failure scenario
function ARROW_BINARY(GRB_ENV, edges, capacity, restored, flows, demand, scenarios, T, Tf)
    println("solving ARROW Binary ILP")

    nedges = size(edges,1)
    nflows = size(flows,1)
    ntunnels = size(T,1)  # all tunnels routing info on IP edge index
    noptions = size(restored[1], 1)
    nscenarios = size(scenarios,1)

    println("Arrow scenario number: ", size(restored,1))
    println("Arrow per-scenario options number: ", size(restored[1],1))

    tunnel_num = 0 # tunnel number for each flow
    for x in 1:length(Tf)
        if size(Tf[x],1) > tunnel_num
            tunnel_num = size(Tf[x],1)
        end
    end

    # residual and restorable tunnels
    Tsf = []  # residual tunnels for all failure scenarios, will just repeat for different restoration options
    Taf = []  # restorable tunnels for all failure scenario and each restoration option
    for s in 1:nscenarios
        sft = []
        aft = []
        for r in 1:size(restored[s], 1)  # for every restore options in a scenario
            ssft = []
            aaft = []
            for f in 1:size(Tf,1)
                ft = []
                at = []
                for t in 1:size(Tf[f],1)
                    up = true
                    partial = false
                    if (length(T[Tf[f][t]]) == 0)
                        up = false
                    end
                    for e in T[Tf[f][t]]
                        if scenarios[s][e] == 0
                            if restored[s][r][e] > 0 # beyond binary failures, less than 1 are all failures
                                up = false
                                partial = true
                            else
                                up = false
                                partial = false
                            end
                        end
                    end
                    if up
                        push!(ft, t)
                    elseif partial
                        push!(at, t)
                    end
                end
                push!(ssft,ft)
                push!(aaft,at)
            end
            push!(sft,ssft)
            push!(aft,aaft)
        end
        push!(Tsf,sft)
        push!(Taf,aft)
    end

    # tunnel routing on IP links
    L = zeros(ntunnels, nedges)
    for t in 1:ntunnels
        for e in 1:nedges
            if in(e, T[t])
                L[t,e] = 1
            end
        end
    end

    model = Model(() -> Gurobi.Optimizer(GRB_ENV))
    set_optimizer_attribute(model, "OutputFlag", 0)
    set_optimizer_attribute(model, "Threads", 32)

    model_def_time = @elapsed begin
    @variable(model, b[1:nflows] >= 0)
    @variable(model, a[1:nflows,1:tunnel_num] >= 0)  #tunnels bw
    @variable(model, m[1:nscenarios,1:noptions] >= 0, Bin)

    ## Equation 28, flow bandwidth constraint
    for f in 1:nflows
        @constraint(model, sum(a[f,t] for t in 1:size(Tf[f],1)) >= b[f])
    end

    # Equation 29, link capacity constraint
    for e in 1:nedges
        @constraint(model, sum(a[f,t] * L[Tf[f][t],e] for f in 1:nflows for t in 1:size(Tf[f],1)) <= capacity[e])
    end

    # Equation 30, flow bandwidth no larger than demand
    for f in 1:nflows
        @constraint(model, b[f] <= demand[f])   #all allocated bandwidths must be less than the demand for that flow
    end

    # Equation 31, residual and restorable tunnels must be able to carry bandwidth for all scenarios
    for f in 1:nflows
        for s in 1:size(scenarios,1)
            for z in 1:size(restored[s], 1)
                @constraint(model, sum(a[f,t] for t in Tsf[s][z][f]) + sum(a[f,t] for t in Taf[s][z][f]) >= b[f] - 100000*(1-m[s,z]))
            end
        end
    end

    # Equation 32, restorable tunnels and the selected restorable bandwidth capacity
    for s in 1:nscenarios
        for e in 1:nedges
            for z in 1:noptions
                @constraint(model, sum(a[f,t] * L[Tf[f][t],e] for f in 1:nflows for t in Taf[s][z][f]) <= restored[s][z][e] + 100000*(1-m[s,z]))
            end
        end
    end

    # Equation 33, lottery ticket selection
    for s in 1:nscenarios
        @constraint(model, sum(m[s,z] for z in 1:noptions) == 1)
    end

    @objective(model, Max, sum(b[i] for i in 1:size(b,1)))
    end  # model_def_time end

    optimize!(model)
    solve_runtime = solve_time(model)
    opt_runtime = solve_runtime + model_def_time

    return value.(a), value.(b), objective_value(model), value.(m), solve_runtime, opt_runtime
end

max_stars_repo_name: hipersys-team/arrow

max_issues_repo_name: arrow-bow/arrow

max_forks_repo_name: arrow-bow/arrow


text:
"""
    fixpoint(g,p0)

adalah fungsi untuk mencari solusi dari `x=g(x)` dengan nilai tebakan awal `p0`

# Examples
```jl
julia> g(x) = exp(-x);

julia> pn, flag, M = fixpoint(g,0.5);

julia> pn
0.5671432633594872

julia> flag
0

julia> M
27×3 Array{Float64,2}:
  0.0  0.5       NaN
  1.0  0.606531    0.106531
  2.0  0.545239    0.0612914
  3.0  0.579703    0.0344639
  ⋮
 23.0  0.567143    4.09741e-7
 24.0  0.567143    2.32382e-7
 25.0  0.567143    1.31794e-7
 26.0  0.567143    7.4746e-8
```
return solusi `pn` dengan `flag` bernilai 0 jika metode bisection berhasil menemukan
solusi dan gagal jika tidak 0. Serta, matriks `M` yang berisi catatan proses tiap
iterasi `[k, pk, f(pk)]`

"""
function fixpoint(g,p0)
    # Definisikan nilai toleransi, maksimum iterasi dan tebakan awal
    delta = 10^-7;
    maxi = 100;
    flag = 1;
    pn = p0
    M = [0 pn NaN];
    for n = 2:maxi
        pn1 = pn;
        pn = g(pn1);
        err = abs(pn-pn1);
        relerr = err/(abs(pn)+eps());
        M = [M; [n-1 pn err]]
        if (err<delta) || (relerr<delta)
            flag = 0; break
        end
    end
    return pn, flag, M
end


text:
using Test
using CLIMAParameters
using CLIMAParameters.Planet

using CLIMAParameters.Atmos.Microphysics_0M
using CLIMAParameters.Atmos.Microphysics

@testset "Microphysics_0M" begin

  struct EarthParameterSet <: AbstractEarthParameterSet end
  ps = EarthParameterSet()

  @test !isnan(Microphysics_0M.τ_precip(ps))
  @test !isnan(Microphysics_0M.qc_0(ps))
  @test !isnan(Microphysics_0M.S_0(ps))

end

@testset "Microphysics" begin

  struct EarthParamSet <: AbstractEarthParameterSet end

  earth = EarthParamSet()

  # Test that all methods are callable, and that nothing returns NaNs
  @test !isnan(Microphysics.C_drag(earth))
  @test !isnan(K_therm(earth))
  @test !isnan(D_vapor(earth))
  @test !isnan(ν_air(earth))
  @test !isnan(N_Sc(earth))

  @test !isnan(τ_cond_evap(earth))

  @test !isnan(τ_sub_dep(earth))
  @test !isnan(r_ice_snow(earth))
  @test !isnan(n0_ice(earth))
  @test !isnan(r0_ice(earth))
  @test !isnan(me_ice(earth))
  @test !isnan(m0_ice(earth))

  @test !isnan(q_liq_threshold(earth))
  @test !isnan(τ_acnv_rai(earth))
  @test !isnan(a_vent_rai(earth))
  @test !isnan(b_vent_rai(earth))
  @test !isnan(n0_rai(earth))
  @test !isnan(r0_ice(earth))
  @test !isnan(me_rai(earth))
  @test !isnan(ae_rai(earth))
  @test !isnan(ve_rai(earth))
  @test !isnan(m0_rai(earth))
  @test !isnan(a0_rai(earth))

  @test !isnan(q_ice_threshold(earth))
  @test !isnan(τ_acnv_sno(earth))
  @test !isnan(a_vent_sno(earth))
  @test !isnan(b_vent_sno(earth))
  @test !isnan(μ_sno(earth))
  @test !isnan(ν_sno(earth))
  @test !isnan(r0_sno(earth))
  @test !isnan(me_sno(earth))
  @test !isnan(ae_sno(earth))
  @test !isnan(ve_sno(earth))
  @test !isnan(m0_sno(earth))
  @test !isnan(a0_sno(earth))
  @test !isnan(v0_sno(earth))

  @test !isnan(E_liq_rai(earth))
  @test !isnan(E_liq_sno(earth))
  @test !isnan(E_ice_rai(earth))
  @test !isnan(E_ice_sno(earth))
  @test !isnan(E_rai_sno(earth))

  # Correctness / relations

  # Atmos.Microphysics
  @test N_Sc(earth)           ≈ ν_air(earth)/D_vapor(earth)
  @test m0_rai(earth) ≈ 4/3. * π * ρ_cloud_liq(earth) * r0_rai(earth)^me_rai(earth)
  @test m0_ice(earth) ≈ 4/3. * π * ρ_cloud_ice(earth) * r0_ice(earth)^me_ice(earth)

  @test χm_ice(earth) ≈ 1
  @test Δm_ice(earth) ≈ 0

  @test χm_rai(earth) ≈ 1
  @test Δm_rai(earth) ≈ 0
  @test χa_rai(earth) ≈ 1
  @test Δa_rai(earth) ≈ 0
  @test χv_rai(earth) ≈ 1
  @test Δv_rai(earth) ≈ 0

  @test χm_sno(earth) ≈ 1
  @test Δm_sno(earth) ≈ 0
  @test χa_sno(earth) ≈ 1
  @test Δa_sno(earth) ≈ 0
  @test χv_sno(earth) ≈ 1
  @test Δv_sno(earth) ≈ 0

end

# Adding convenience methods:
abstract type Phase end
struct Liquid <: Phase end
struct Rain   <: Phase end
struct Snow   <: Phase end

E(earth::AbstractEarthParameterSet, ::Liquid, ::Rain) = E_liq_rai(earth)
E(earth::AbstractEarthParameterSet, ::Liquid, ::Snow) = E_liq_sno(earth)

@testset "Microphysics" begin

  struct EarthParamSet <: AbstractEarthParameterSet end

  earth = EarthParamSet()

  @test !isnan(E(earth, Liquid(), Rain()))
  @test !isnan(E(earth, Liquid(), Snow()))

end


text:
using DataDeps
using Plots
using OceanParameterizations
using FreeConvection
using FreeConvection: coarse_grain

ENV["DATADEPS_ALWAYS_ACCEPT"] = "true"
ENV["GKSwstype"] = "100"

Nz = 32

## Register data dependencies

@info "Registering data dependencies..."
for dd in FreeConvection.LESBRARY_DATA_DEPS
    DataDeps.register(dd)
end

## Load data

@info "Loading data..."
datasets = Dict{Int,Any}(
    1 => NCDstack(datadep"free_convection_8days_Qb1e-8/statistics.nc"),
    2 => NCDstack(datadep"free_convection_8days_Qb2e-8/statistics.nc"),
    3 => NCDstack(datadep"free_convection_8days_Qb3e-8/statistics.nc"),
    4 => NCDstack(datadep"free_convection_8days_Qb4e-8/statistics.nc"),
    5 => NCDstack(datadep"free_convection_8days_Qb5e-8/statistics.nc"),
    6 => NCDstack(datadep"free_convection_8days_Qb6e-8/statistics.nc")
)

## Add surface fluxes to data

@info "Inserting surface fluxes..."
datasets = Dict{Int,Any}(id => add_surface_fluxes(ds) for (id, ds) in datasets)

## Coarse grain training data

@info "Coarse graining data..."
coarse_datasets = Dict{Int,Any}(id => coarse_grain(ds, Nz) for (id, ds) in datasets)

## Split into training and testing data

@info "Partitioning data into training and testing datasets..."

ids_train = [1, 2, 4, 6]
ids_test = [3, 5]

training_datasets = Dict(id => datasets[id] for id in ids_train)
testing_datasets = Dict(id => datasets[id] for id in ids_test)

coarse_training_datasets = Dict(id => coarse_datasets[id] for id in ids_train)
coarse_testing_datasets = Dict(id => coarse_datasets[id] for id in ids_test)

## Pull out input (T) and output (wT) training data

@info "Wrangling training data..."
input_training_data = wrangle_input_training_data(coarse_training_datasets)
output_training_data = wrangle_output_training_data(coarse_training_datasets)

## Feature scaling

@info "Scaling features..."

T_training_data = reduce(hcat, input.temperature for input in input_training_data)
wT_training_data = output_training_data

@assert size(wT_training_data, 1) == size(T_training_data, 1) + 1
@assert size(wT_training_data, 2) == size(T_training_data, 2)

T_scaling = ZeroMeanUnitVarianceScaling(T_training_data)
wT_scaling = ZeroMeanUnitVarianceScaling(wT_training_data)

## Optimizing KPP parameters

@info "Optimizing KPP parameters..."
true_solutions = Dict(id => T_scaling.(ds[:T].data) for (id, ds) in coarse_datasets)

# eki_ensemble_size = 10
# eki_iterations = 10
# eki, eki_loss_history = optimize_kpp_parameters(coarse_training_datasets, true_solutions, T_scaling,
#                                                 ensemble_members=eki_ensemble_size, iterations=eki_iterations)

# plot(eki_loss_history, linewidth=3, linealpha=0.8, yaxis=:log,
#      label="", xlabel="EKI iteration", ylabel="mean squared error",
#      title="Optimizing KPP: EKI particle loss", grid=false, framestyle=:box, dpi=200)

# savefig("eki_kpp_loss_history.png")

# kwargs = (label="", grid=false, framestyle=:box)

# anim = @animate for n in 1:eki_iterations
#     h1 = histogram(eki.u[n][:, 1], bins=0:0.1:1, xlims=(0, 1), xlabel="CSL", title="Optimizing KPP parameters"; kwargs...)
#     h2 = histogram(eki.u[n][:, 2], bins=0:0.8:8, xlims=(0, 8), xlabel="CNL", title="EKI iteration $n"; kwargs...)
#     h3 = histogram(eki.u[n][:, 3], bins=0:0.6:6, xlims=(0, 6), xlabel="CbT"; kwargs...)
#     h4 = histogram(eki.u[n][:, 4], bins=0:0.5:5, xlims=(0, 5), xlabel="CKE"; kwargs...)
#     plot(h1, h2, h3, h4, layout=(2, 2), dpi=200)
# end

# gif(anim, "eki_kpp_parameters_histograms.gif", fps=1)

using Random, Distributions, LinearAlgebra
# numerical gradient
function gradient(L, params; δ = 1e-4 .* ones(length(params)))
    ∇L = zeros(length(params))
    e = I + zeros(length(params), length(params))
    Lcurrent = L(params)
    for i in eachindex(params)
        up = L(params + δ[i] * e[i,:])
        ∇L[i] = (up - Lcurrent)/(δ[i])
    end
    return ∇L
end

# Define Method structs
struct RandomPlugin{S,T,V,U}
    priors::S
    fcalls::T
    seed::V
    progress::U
end

function RandomPlugin(priors, fcalls::Int; seed = 1234, progress = true)
    return RandomPlugin(priors, fcalls, 1234, true)
end

struct RandomLineSearch{I, T, B}
    linesearches::I
    linelength::I
    linebounds::T
    progress::B
    seed::I
end

function RandomLineSearch(linesearches, linelength)
    return RandomLineSearch(linesearches, linelength, (-0.1,1), true, 1234)
end
function RandomLineSearch(; linesearches = 10, linelength = 10, linebounds = (-0.1,1), progress = true, seed = 1234)
    return RandomLineSearch(linesearches, linelength, linebounds, progress, seed)
end

# Define Helper functions

# RandomPlugin
function priorloss(ℒ, fcalls, priors; seed = 1234, progress = true)
    Random.seed!(seed)
    losses = []
    vals = []
    for i in 1:fcalls
        guessparams = rand.(priors)
        push!(vals, guessparams)
        currentloss = ℒ(guessparams)
        push!(losses, currentloss)
        if progress
            println("iteration " * string(i))
        end
    end
    return losses, vals
end

# LineSearch
function randomlinesearch(ℒ, ∇ℒ, bparams; linesearches = 10, linelength = 10, linebounds = (-0.1, 1.0), progress = true, seed = 1234)
    params = copy(bparams)
    Random.seed!(seed)
    for i in 1:linesearches
        αs = [rand(Uniform(linebounds...), linelength-1)..., 0.0]
        ∇L = ∇ℒ(params)
        LS = [ℒ(params - α * ∇L) for α in αs]
        b = argmin(LS)
        if progress
            println("best after line search ", i)
        end
        params .= params - αs[b] * ∇L
        if αs[b] == 0
            αs .= 0.5 .* αs
            if progress
                println("staying the same at ")
            end
        end
        if progress
            println((params, LS[b]))
        end
    end
    return params
end

# Define optimize functions
function optimize(ℒ, method::RandomPlugin; history = false, printresult = true)
    losses, vals = priorloss(ℒ, method.fcalls,
                       method.priors, seed = method.seed,
                       progress = method.progress)
    indmin = argmin(losses)
    if printresult
        println("The minimum loss is ", losses[indmin])
        println("The minimum argument is ", vals[indmin])
        println("This occured at iteration ", indmin)
    end
    if history == false
        return vals[indmin]
    else
        return vals, losses
    end
end
optimize(ℒ, initialguess, method::RandomPlugin; history = false, printresult = true) = optimize(ℒ, method::RandomPlugin; history = false, printresult = true)

function optimize(ℒ, ∇ℒ, params, method::RandomLineSearch; history = false, printresult = true)
    bestparams = randomlinesearch(ℒ, ∇ℒ, params; linesearches = method.linesearches,
                    linelength = method.linelength,
                    linebounds = method.linebounds,
                    progress = method.progress,
                    seed = method.seed)
    return bestparams
end

using Flux: mse
import OceanTurb

## Example
# loss function

# loss(x) = (x[1] - 1)^2 + (x[2] - 2)^2 + (x[3] - 3)^2 + (x[4]-4)^2

# loss(p) = mean(
#     mse(T_scaling.(free_convection_kpp(ds, parameters=OceanTurb.KPP.Parameters(CSL=p[1], CNL=p[2], Cb_T=p[3], CKE=p[4])).T),
#         true_solutions[id])
#     for (id, ds) in datasets
# )

loss(p) = mse(T_scaling.(free_convection_kpp(coarse_datasets[6], parameters=OceanTurb.KPP.Parameters(CSL=p[1], CNL=p[2], Cb_T=p[3], CKE=p[4])).T),
                         true_solutions[6])

# First construct global search
# Create Prior
lower = [0.0, 0.0, 0.0,  0.0]
upper = [1.0, 10.0, 10.0, 10.0]
priors = Uniform.(lower, upper)
# Determine number of function calls
functioncalls = 1000
# Define Method
method = RandomPlugin(priors, functioncalls)
# Optimize
minparam = optimize(loss, method)

# Next do gradient descent
# construct numerical gradient
∇loss(params) = gradient(loss, params)
# optimize choosing minimum from the global search for refinement
best_params = minparam
method  = RandomLineSearch(linebounds = (0, 1e-0/norm(∇loss(best_params))), linesearches = 20)
bestparam = optimize(loss, ∇loss, best_params, method)


text:
##################################################
# Helper functions
##################################################

_elementtuple(t::HomoNuclearTable) = (t.element_symbol, t.element_symbol)
_elementtuple(t::HeteroNuclearTable) = (t.first_element_symbol, t.second_element_symbol)

##################################################

"""
    loaddir([T], dirpath::String)

Parse a directory at `dirpath` containing SKF files and return a
`SlaterKosterTable` derived from its contents.

See also [`load`](@ref SlaterKoster.load).

"""
loaddir(dirpath) = loaddir(Float64, dirpath)
function loaddir(T::Type, dirpath)
    data = Dict{Tuple{Symbol,Symbol},PrimitiveTable{T}}()
    for file in filter(f -> occursin(r"(skf|SKF)$", f), Base.readdir(dirpath))
        table = load(T, joinpath(dirpath, file))
        data[_elementtuple(table)] = table
    end
    return SlaterKosterTable(data)
end

"""
    load([T], filepath::String; [element_symbols::NTuple{2,String}])

Parse an SKF located at `filepath` and return a primitive table
derived from its contents.  The element symbols will be inferred from
the file name if it follows the standard `X-Y.skf` format, where `X`,
and `Y` are element symbols.  Alternatively, if a tuple of symbols is
provided to the keyword argument `element_symbols`, the element
symbols listed in the tuple will be used instead.  If `T` is supplied,
floating-point data will be parsed using the type `T` instead of the
default `Float64`.

See also [`loaddir`](@ref SlaterKoster.loaddir).

"""
load(filepath; element_symbols=(nothing, nothing)) = load(Float64, filepath, element_symbols=element_symbols)
function load(T::Type, filepath::String; element_symbols=(nothing, nothing))
    dirname, filename = splitdir(filepath)
    if element_symbols === (nothing, nothing)
        # Try to derive the element symbols from the file name by
        # matching file name patterns of known datasets
        m = match(r"^([A-Z][a-z]*)[-]*([A-Z][a-z]*)\.(skf|SKF)$", filename)
        m === nothing && error("could not determine element symbols from SKF file ($filename); please specify them manually using the keyword argument `element_symbols`")
        first_element_symbol, second_element_symbol, _ = m.captures
    else
        first_element_symbol, second_element_symbol = element_symbols
    end
    lines = Base.readlines(filepath)
    if first_element_symbol == second_element_symbol
        if startswith(lines[1], "@")
            return HomoNuclearTable(
                Symbol(first_element_symbol),
                _readextendedformat_homonuclear(T, lines)...)
        else
            return HomoNuclearTable(
                Symbol(first_element_symbol),
                _readsimpleformat_homonuclear(T, lines)...)
        end
    else
        if startswith(lines[1], "@")
            return HeteroNuclearTable(
                Symbol(first_element_symbol),
                Symbol(second_element_symbol),
                _readextendedformat_heteronuclear(T, lines)...)
        else
            return HeteroNuclearTable(
                Symbol(first_element_symbol),
                Symbol(second_element_symbol),
                _readsimpleformat_heteronuclear(T, lines)...)
        end
    end
end

function parseheader(T::Type, line)
    words = filter(word -> !all(isspace, word), split(line, x -> isspace(x) || x === ','))
    return parse(T, words[1]), parse(Int, words[2])
end

function parseline(T::Type, line)
    words = filter(word -> !all(isspace, word), split(line, x -> isspace(x) || x === ','))
    i = 1
    while i <= length(words)
        m = match(r"([0-9]+)\*([0-9\.]+)", words[i])
        if m !== nothing
            deleteat!(words, i)
            for _ in 1:parse(Int, m[1])
                insert!(words, i, m[2])
            end
        end
        i += 1
    end
    return parse.(T, words)
end

function _readsimpleformat_homonuclear(T, lines)
    grid_dist, n_grid_points = parseheader(T, lines[1])

    energy_table = DataFrame(
        Ed = T[], Ep = T[], Es = T[],
        SPE = T[],
        Ud = T[], Up = T[], Us = T[],
        fd = T[], fp = T[], fs = T[])
    push!(energy_table, parseline(T, lines[2]))

    mass_table = DataFrame(
        mass = T[],
        c2 = T[], c3 = T[], c4 = T[], c5 = T[],
        c6 = T[], c7 = T[], c8 = T[], c9 = T[],
        rcut = T[],
        d1 = T[], d2 = T[], d3 = T[], d4 = T[], d5 = T[],
        d6 = T[], d7 = T[], d8 = T[], d9 = T[], d10 = T[])
    push!(mass_table, parseline(T, lines[3]))

    integral_table = DataFrame(
        dist = T[],
        Hdd0 = T[], Hdd1 = T[], Hdd2 = T[], Hpd0 = T[], Hpd1 = T[], Hpp0 = T[],
        Hpp1 = T[], Hsd0 = T[], Hsp0 = T[], Hss0 = T[],
        Sdd0 = T[], Sdd1 = T[], Sdd2 = T[], Spd0 = T[], Spd1 = T[], Spp0 = T[],
        Spp1 = T[], Ssd0 = T[], Ssp0 = T[], Sss0 = T[])

    for (i, line) in enumerate(lines[4:4+n_grid_points-1])
        push!(integral_table, vcat([i*grid_dist], parseline(T, line)))
    end
    return grid_dist, n_grid_points, energy_table, mass_table, integral_table
end

function _readsimpleformat_heteronuclear(T, lines)
    grid_dist, n_grid_points = parseheader(T, lines[1])

    mass_table = DataFrame(
        mass = T[],
        c2 = T[], c3 = T[], c4 = T[], c5 = T[],
        c6 = T[], c7 = T[], c8 = T[], c9 = T[],
        rcut = T[],
        d1 = T[], d2 = T[], d3 = T[], d4 = T[], d5 = T[],
        d6 = T[], d7 = T[], d8 = T[], d9 = T[], d10 = T[])
    push!(mass_table, parseline(T, lines[2]))

    integral_table = DataFrame(
        dist = T[],
        Hdd0 = T[], Hdd1 = T[], Hdd2 = T[], Hpd0 = T[], Hpd1 = T[], Hpp0 = T[],
        Hpp1 = T[], Hsd0 = T[], Hsp0 = T[], Hss0 = T[],
        Sdd0 = T[], Sdd1 = T[], Sdd2 = T[], Spd0 = T[], Spd1 = T[], Spp0 = T[],
        Spp1 = T[], Ssd0 = T[], Ssp0 = T[], Sss0 = T[])

    for (i, line) in enumerate(lines[3:3+n_grid_points-1])
        push!(integral_table, vcat([i*grid_dist], parseline(T, line)))
    end
    return grid_dist, n_grid_points, mass_table, integral_table
end

function _readextendedformat_homonuclear(T, lines)
    grid_dist, n_grid_points = parseheader(T, lines[1])

    energy_table = DataFrame(
        Ef = T[], Ed = T[], Ep = T[], Es = T[],
        SPE = T[],
        Uf = T[], Ud = T[], Up = T[], Us = T[],
        ff = T[], fd = T[], fp = T[], fs = T[])
    push!(energy_table, parseline(T, lines[3]))

    mass_table = DataFrame(
        mass = T[],
        c2 = T[], c3 = T[], c4 = T[], c5 = T[],
        c6 = T[], c7 = T[], c8 = T[], c9 = T[],
        rcut = T[],
        d1 = T[], d2 = T[], d3 = T[], d4 = T[], d5 = T[],
        d6 = T[], d7 = T[], d8 = T[], d9 = T[], d10 = T[])
    push!(mass_table, parseline(T, lines[4]))

    integral_table = DataFrame(
        dist = T[],
        Hff0 = T[], Hff1 = T[], Hff2 = T[], Hff3 = T[],
        Hdf0 = T[], Hdf1 = T[], Hdf2 = T[], Hdd0 = T[],
        Hdd1 = T[], Hdd2 = T[], Hpf0 = T[], Hpf1 = T[],
        Hpd0 = T[], Hpd1 = T[], Hpp0 = T[], Hpp1 = T[],
        Hsf0 = T[], Hsd0 = T[], Hsp0 = T[], Hss0 = T[],
        Sdd1 = T[], Sdd2 = T[], Sff0 = T[], Sff1 = T[],
        Sff2 = T[], Sff3 = T[], Sdf0 = T[], Sdf1 = T[],
        Sdf2 = T[], Sdd0 = T[], Spf0 = T[], Spf1 = T[],
        Spd0 = T[], Spd1 = T[], Spp0 = T[], Spp1 = T[],
        Ssf0 = T[], Ssd0 = T[], Ssp0 = T[], Sss0 = T[])

    for (i, line) in enumerate(lines[5:5+n_grid_points-1])
        push!(integral_table, vcat([i*grid_dist], parseline(T, line)))
    end
    return grid_dist, n_grid_points, energy_table, mass_table, integral_table
end

function _readextendedformat_heteronuclear(T, lines)
    grid_dist, n_grid_points = parseheader(T, lines[1])

    mass_table = DataFrame(
        mass = T[],
        c2 = T[], c3 = T[], c4 = T[], c5 = T[],
        c6 = T[], c7 = T[], c8 = T[], c9 = T[],
        rcut = T[],
        d1 = T[], d2 = T[], d3 = T[], d4 = T[], d5 = T[],
        d6 = T[], d7 = T[], d8 = T[], d9 = T[], d10 = T[])
    push!(mass_table, parseline(T, lines[3]))

    integral_table = DataFrame(
        dist = T[],
        Hff0 = T[], Hff1 = T[], Hff2 = T[], Hff3 = T[],
        Hdf0 = T[], Hdf1 = T[], Hdf2 = T[], Hdd0 = T[],
        Hdd1 = T[], Hdd2 = T[], Hpf0 = T[], Hpf1 = T[],
        Hpd0 = T[], Hpd1 = T[], Hpp0 = T[], Hpp1 = T[],
        Hsf0 = T[], Hsd0 = T[], Hsp0 = T[], Hss0 = T[],
        Sdd1 = T[], Sdd2 = T[], Sff0 = T[], Sff1 = T[],
        Sff2 = T[], Sff3 = T[], Sdf0 = T[], Sdf1 = T[],
        Sdf2 = T[], Sdd0 = T[], Spf0 = T[], Spf1 = T[],
        Spd0 = T[], Spd1 = T[], Spp0 = T[], Spp1 = T[],
        Ssf0 = T[], Ssd0 = T[], Ssp0 = T[], Sss0 = T[])

    for (i, line) in enumerate(lines[4:4+n_grid_points-1])
        push!(integral_table, vcat([i*grid_dist], parseline(T, line)))
    end
    return grid_dist, n_grid_points, mass_table, integral_table
end


text:
#=
    Here we check that the result of computation is
    of the desired type in respect to input types
=#

using AbstractAlgebra

@testset "AbstractAlgebra io consistency" begin
    R, (x, y) = PolynomialRing(GF(2^31 - 1), ["x", "y"], ordering=:lex)
    fs = [x^2*y + 3, (2^31 - 5)*x - (2^31 - 4)*y]
    gb = Groebner.groebner(fs)
    @test parent(gb[1]) == R

    R, (x, y) = PolynomialRing(QQ, ["x", "y"], ordering=:lex)
    fs = [x^2*y + 3//2, (2^31 - 5)*x - (2^31 - 4)*y]
    gb = Groebner.groebner(fs)
    @test parent(gb[1]) == R

    R, (x, y) = PolynomialRing(ZZ, ["x", "y"], ordering=:lex)
    fs = [x^2*y + 3, (2^31 - 5)*x - (2^31 - 4)*y]
    gb = Groebner.groebner(fs)
    @test parent(gb[1]) == R

    R, (x, y) = PolynomialRing(GF(2^31 - 1), ["x", "y"], ordering=:degrevlex)
    fs = [x^2*y + 3, (2^31 - 5)*x - (2^31 - 4)*y]
    gb = Groebner.groebner(fs)
    @test parent(gb[1]) == R

    R, (x, y) = PolynomialRing(QQ, ["x", "y"], ordering=:degrevlex)
    fs = [x^2*y + 3//2, (2^31 - 5)*x - (2^31 - 4)*y]
    gb = Groebner.groebner(fs)
    @test parent(gb[1]) == R

    R, (x, y) = PolynomialRing(ZZ, ["x", "y"], ordering=:degrevlex)
    fs = [x^2*y + 3, (2^31 - 5)*x - (2^31 - 4)*y]
    gb = Groebner.groebner(fs)
    @test parent(gb[1]) == R
end


text:
using Flux
using Flux: train!, @epochs
using Flux: Data.DataLoader

actual(x) = 4x + 2

x_train, x_test = hcat(0:5...), hcat(6:10...)
y_train, y_test = actual.(x_train), actual.(x_test)
data = DataLoader((x_train, y_train), shuffle=true)

model = Dense(1, 1)
loss(x, y) = Flux.Losses.mse(model(x), y)

opt = Descent()

@epochs 200 begin
  train!(loss, params(model), data, opt)
  @show(loss(x_test, y_test))
  @show(params(model))
end


text:
module Powseeker

using StaticArrays
using ParticleFilters
using POMDPs
using Plots
using Parameters
using POMDPToolbox

import POMDPs: generate_sr, generate_s
import POMDPs: isterminal, initial_state_distribution, actions
import POMDPs: observation
import POMDPs: pdf
import POMDPs: action
import Base: rand, eltype

const Vec2 = SVector{2,Float64}

export
    PowseekerMDP,
    PowseekerPOMDP,
    SkierState,
    GPSOrAngle,
    SkierObs,
    SkierUnif,

    peaks,
    scaled_peaks,
    grad_peaks,
    grad_scaled_peaks,
    peaks2,
    scaled_peaks2,
    grad_peaks2,
    grad_scaled_peaks2,
    plot_peaks,
    plot_scaled_peaks,

    mdp,

    Downhill,
    Topout,
    SkiOver,
    RandomlyCheckGPS

include("peaks.jl")

# package code goes here
immutable SkierState
    time::Int
    pos::Vec2
    psi::Float64
end

immutable GPSOrAngle
    gps::Bool
    angle::Float64
end

immutable SkierObs
    time::Int
    pos::Nullable{Vec2}
    dir::Nullable{Float64}
    dist::Float64
    grad::Vec2

    SkierObs(t::Int, p::Vec2, dir::Float64) = new(t, Nullable{Vec2}(p), Nullable{Float64}(dir))
    SkierObs(t::Int, d::Float64, g::Vec2) = new(t, Nullable{Vec2}(), Nullable{Float64}(), d, g)
end

@with_kw immutable PowseekerMDP{G<:Function} <: MDP{SkierState, Float64}
    topology::Function      = scaled_peaks2
    gradient::G             = grad_scaled_peaks2
    duration::Int           = 60
    dt::Float64             = 2.0
    step_len::Float64       = 60.0
    force::Float64          = 4.0
    terminal_vel::Float64   = 50.0
    max_flat_speed::Float64 = 4.0
    psi_std::Float64        = deg2rad(0.5)
    vel_std::Float64        = 0.05
    xlim::NTuple{2,Float64} = (-4000.0, 4000.0)
    ylim::NTuple{2,Float64} = (-4000.0, 4000.0)
end

@with_kw immutable PowseekerPOMDP{G<:Function} <: POMDP{SkierState, GPSOrAngle, SkierObs}
    mdp::PowseekerMDP{G}    = PowseekerMDP()
    dist_std_frac::Float64  = 0.3
    grad_std::Float64       = 0.2
    gps_std::Float64        = 50.0
    compass_std::Float64    = deg2rad(10.0)
end

typealias PowseekerProblem Union{PowseekerMDP, PowseekerPOMDP}
mdp(p::PowseekerMDP) = p
mdp(p::PowseekerPOMDP) = p.mdp

function generate_sr(pp::PowseekerProblem, s::SkierState, a::Float64, rng::AbstractRNG)
    p = mdp(pp)
    dt = p.dt
    mfs = p.max_flat_speed
    vel = 0.0
    pos = s.pos
    psi = s.psi + a
    isteps = round(Int, p.step_len/dt)
    r = 0.0
    for i in 1:isteps
        dir = SVector(cos(psi), sin(psi))
        d = SkiLineDynamics(p, pos, dir)
        vec = rk4step(d, vel)
        vel = vec[1]+randn(rng)*p.vel_std
        pos += vec[2]*dir
        psi += randn(rng)*p.psi_std
        r += exp(vel)
        # r += vel^8
    end
    sp = SkierState(s.time+1, pos, psi)
    r /= isteps
    return (sp, r)
end

function generate_s(p::PowseekerProblem, s::SkierState, a::Float64, rng::AbstractRNG)
    sp, r = generate_sr(mdp(p), s, a, rng)
    return sp
end

function generate_sr(p::PowseekerPOMDP, s::SkierState, a::GPSOrAngle, rng::AbstractRNG)
    if a.gps
        return SkierState(s.time+1, s.pos, s.psi), 0.0
    else
        return generate_sr(mdp(p), s, a.angle, rng)
    end
end

function generate_s(p::PowseekerPOMDP, s::SkierState, a::GPSOrAngle, rng::AbstractRNG)
    if a.gps
        return SkierState(s.time+1, s.pos, s.psi)
    else
        return generate_s(mdp(p), s, a.angle, rng)
    end
end

function observation(p::PowseekerPOMDP, s::SkierState, a::GPSOrAngle, sp::SkierState)
    return SkierObsDist(p, a.gps, s, sp)
end

isterminal(p::PowseekerProblem, s::SkierState) = s.time > mdp(p).duration

include("runge-kutta.jl")
include("actions.jl")
include("observations.jl")
include("initial.jl")
include("heuristics.jl")
include("visualization.jl")

end # module


text:
# Julia wrapper for header: /home/astukalov/.julia/v0.5/Sundials/deps/usr/include/nvector/nvector_parallel.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

# Julia wrapper for header: /home/astukalov/.julia/v0.5/Sundials/deps/usr/include/nvector/nvector_serial.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function __N_VNew_Serial(vec_length::Clong)
    ccall((:N_VNew_Serial,libsundials_nvecserial),N_Vector,(Clong,),vec_length)
end

function N_VNew_Serial(vec_length)
    __N_VNew_Serial(convert(Clong,vec_length))
end

function __N_VNewEmpty_Serial(vec_length::Clong)
    ccall((:N_VNewEmpty_Serial,libsundials_nvecserial),N_Vector,(Clong,),vec_length)
end

function N_VNewEmpty_Serial(vec_length)
    __N_VNewEmpty_Serial(convert(Clong,vec_length))
end

function __N_VMake_Serial(vec_length::Clong,v_data::Ptr{realtype})
    ccall((:N_VMake_Serial,libsundials_nvecserial),N_Vector,(Clong,Ptr{realtype}),vec_length,v_data)
end

function N_VMake_Serial(vec_length,v_data)
    __N_VMake_Serial(convert(Clong,vec_length),pointer(v_data))
end

function __N_VCloneVectorArray_Serial(count::Cint,w::N_Vector)
    ccall((:N_VCloneVectorArray_Serial,libsundials_nvecserial),Ptr{N_Vector},(Cint,N_Vector),count,w)
end

function N_VCloneVectorArray_Serial(count,w)
    __w = convert(NVector,w)
    __N_VCloneVectorArray_Serial(convert(Cint,count),convert(N_Vector,__w))
end

function __N_VCloneVectorArrayEmpty_Serial(count::Cint,w::N_Vector)
    ccall((:N_VCloneVectorArrayEmpty_Serial,libsundials_nvecserial),Ptr{N_Vector},(Cint,N_Vector),count,w)
end

function N_VCloneVectorArrayEmpty_Serial(count,w)
    __w = convert(NVector,w)
    __N_VCloneVectorArrayEmpty_Serial(convert(Cint,count),convert(N_Vector,__w))
end

function __N_VDestroyVectorArray_Serial(vs::Ptr{N_Vector},count::Cint)
    ccall((:N_VDestroyVectorArray_Serial,libsundials_nvecserial),Void,(Ptr{N_Vector},Cint),vs,count)
end

function N_VDestroyVectorArray_Serial(vs,count)
    __N_VDestroyVectorArray_Serial(pointer(vs),convert(Cint,count))
end

function __N_VPrint_Serial(v::N_Vector)
    ccall((:N_VPrint_Serial,libsundials_nvecserial),Void,(N_Vector,),v)
end

function N_VPrint_Serial(v)
    __v = convert(NVector,v)
    __N_VPrint_Serial(convert(N_Vector,__v))
end

function __N_VCloneEmpty_Serial(w::N_Vector)
    ccall((:N_VCloneEmpty_Serial,libsundials_nvecserial),N_Vector,(N_Vector,),w)
end

function N_VCloneEmpty_Serial(w)
    __w = convert(NVector,w)
    __N_VCloneEmpty_Serial(convert(N_Vector,__w))
end

function __N_VClone_Serial(w::N_Vector)
    ccall((:N_VClone_Serial,libsundials_nvecserial),N_Vector,(N_Vector,),w)
end

function N_VClone_Serial(w)
    __w = convert(NVector,w)
    __N_VClone_Serial(convert(N_Vector,__w))
end

function __N_VDestroy_Serial(v::N_Vector)
    ccall((:N_VDestroy_Serial,libsundials_nvecserial),Void,(N_Vector,),v)
end

function N_VDestroy_Serial(v)
    __v = convert(NVector,v)
    __N_VDestroy_Serial(convert(N_Vector,__v))
end

function __N_VSpace_Serial(v::N_Vector,lrw::Ptr{Clong},liw::Ptr{Clong})
    ccall((:N_VSpace_Serial,libsundials_nvecserial),Void,(N_Vector,Ptr{Clong},Ptr{Clong}),v,lrw,liw)
end

function N_VSpace_Serial(v,lrw,liw)
    __v = convert(NVector,v)
    __N_VSpace_Serial(convert(N_Vector,__v),pointer(lrw),pointer(liw))
end

function __N_VGetArrayPointer_Serial(v::N_Vector)
    ccall((:N_VGetArrayPointer_Serial,libsundials_nvecserial),Ptr{realtype},(N_Vector,),v)
end

function N_VGetArrayPointer_Serial(v)
    __v = convert(NVector,v)
    __N_VGetArrayPointer_Serial(convert(N_Vector,__v))
end

function __N_VSetArrayPointer_Serial(v_data::Ptr{realtype},v::N_Vector)
    ccall((:N_VSetArrayPointer_Serial,libsundials_nvecserial),Void,(Ptr{realtype},N_Vector),v_data,v)
end

function N_VSetArrayPointer_Serial(v_data,v)
    __v = convert(NVector,v)
    __N_VSetArrayPointer_Serial(pointer(v_data),convert(N_Vector,__v))
end

function __N_VLinearSum_Serial(a::realtype,x::N_Vector,b::realtype,y::N_Vector,z::N_Vector)
    ccall((:N_VLinearSum_Serial,libsundials_nvecserial),Void,(realtype,N_Vector,realtype,N_Vector,N_Vector),a,x,b,y,z)
end

function N_VLinearSum_Serial(a,x,b,y,z)
    __x = convert(NVector,x)
    __y = convert(NVector,y)
    __z = convert(NVector,z)
    __N_VLinearSum_Serial(a,convert(N_Vector,__x),b,convert(N_Vector,__y),convert(N_Vector,__z))
end

function __N_VConst_Serial(c::realtype,z::N_Vector)
    ccall((:N_VConst_Serial,libsundials_nvecserial),Void,(realtype,N_Vector),c,z)
end

function N_VConst_Serial(c,z)
    __z = convert(NVector,z)
    __N_VConst_Serial(c,convert(N_Vector,__z))
end

function __N_VProd_Serial(x::N_Vector,y::N_Vector,z::N_Vector)
    ccall((:N_VProd_Serial,libsundials_nvecserial),Void,(N_Vector,N_Vector,N_Vector),x,y,z)
end

function N_VProd_Serial(x,y,z)
    __x = convert(NVector,x)
    __y = convert(NVector,y)
    __z = convert(NVector,z)
    __N_VProd_Serial(convert(N_Vector,__x),convert(N_Vector,__y),convert(N_Vector,__z))
end

function __N_VDiv_Serial(x::N_Vector,y::N_Vector,z::N_Vector)
    ccall((:N_VDiv_Serial,libsundials_nvecserial),Void,(N_Vector,N_Vector,N_Vector),x,y,z)
end

function N_VDiv_Serial(x,y,z)
    __x = convert(NVector,x)
    __y = convert(NVector,y)
    __z = convert(NVector,z)
    __N_VDiv_Serial(convert(N_Vector,__x),convert(N_Vector,__y),convert(N_Vector,__z))
end

function __N_VScale_Serial(c::realtype,x::N_Vector,z::N_Vector)
    ccall((:N_VScale_Serial,libsundials_nvecserial),Void,(realtype,N_Vector,N_Vector),c,x,z)
end

function N_VScale_Serial(c,x,z)
    __x = convert(NVector,x)
    __z = convert(NVector,z)
    __N_VScale_Serial(c,convert(N_Vector,__x),convert(N_Vector,__z))
end

function __N_VAbs_Serial(x::N_Vector,z::N_Vector)
    ccall((:N_VAbs_Serial,libsundials_nvecserial),Void,(N_Vector,N_Vector),x,z)
end

function N_VAbs_Serial(x,z)
    __x = convert(NVector,x)
    __z = convert(NVector,z)
    __N_VAbs_Serial(convert(N_Vector,__x),convert(N_Vector,__z))
end

function __N_VInv_Serial(x::N_Vector,z::N_Vector)
    ccall((:N_VInv_Serial,libsundials_nvecserial),Void,(N_Vector,N_Vector),x,z)
end

function N_VInv_Serial(x,z)
    __x = convert(NVector,x)
    __z = convert(NVector,z)
    __N_VInv_Serial(convert(N_Vector,__x),convert(N_Vector,__z))
end

function __N_VAddConst_Serial(x::N_Vector,b::realtype,z::N_Vector)
    ccall((:N_VAddConst_Serial,libsundials_nvecserial),Void,(N_Vector,realtype,N_Vector),x,b,z)
end

function N_VAddConst_Serial(x,b,z)
    __x = convert(NVector,x)
    __z = convert(NVector,z)
    __N_VAddConst_Serial(convert(N_Vector,__x),b,convert(N_Vector,__z))
end

function __N_VDotProd_Serial(x::N_Vector,y::N_Vector)
    ccall((:N_VDotProd_Serial,libsundials_nvecserial),realtype,(N_Vector,N_Vector),x,y)
end

function N_VDotProd_Serial(x,y)
    __x = convert(NVector,x)
    __y = convert(NVector,y)
    __N_VDotProd_Serial(convert(N_Vector,__x),convert(N_Vector,__y))
end

function __N_VMaxNorm_Serial(x::N_Vector)
    ccall((:N_VMaxNorm_Serial,libsundials_nvecserial),realtype,(N_Vector,),x)
end

function N_VMaxNorm_Serial(x)
    __x = convert(NVector,x)
    __N_VMaxNorm_Serial(convert(N_Vector,__x))
end

function __N_VWrmsNorm_Serial(x::N_Vector,w::N_Vector)
    ccall((:N_VWrmsNorm_Serial,libsundials_nvecserial),realtype,(N_Vector,N_Vector),x,w)
end

function N_VWrmsNorm_Serial(x,w)
    __x = convert(NVector,x)
    __w = convert(NVector,w)
    __N_VWrmsNorm_Serial(convert(N_Vector,__x),convert(N_Vector,__w))
end

function __N_VWrmsNormMask_Serial(x::N_Vector,w::N_Vector,id::N_Vector)
    ccall((:N_VWrmsNormMask_Serial,libsundials_nvecserial),realtype,(N_Vector,N_Vector,N_Vector),x,w,id)
end

function N_VWrmsNormMask_Serial(x,w,id)
    __x = convert(NVector,x)
    __w = convert(NVector,w)
    __id = convert(NVector,id)
    __N_VWrmsNormMask_Serial(convert(N_Vector,__x),convert(N_Vector,__w),convert(N_Vector,__id))
end

function __N_VMin_Serial(x::N_Vector)
    ccall((:N_VMin_Serial,libsundials_nvecserial),realtype,(N_Vector,),x)
end

function N_VMin_Serial(x)
    __x = convert(NVector,x)
    __N_VMin_Serial(convert(N_Vector,__x))
end

function __N_VWL2Norm_Serial(x::N_Vector,w::N_Vector)
    ccall((:N_VWL2Norm_Serial,libsundials_nvecserial),realtype,(N_Vector,N_Vector),x,w)
end

function N_VWL2Norm_Serial(x,w)
    __x = convert(NVector,x)
    __w = convert(NVector,w)
    __N_VWL2Norm_Serial(convert(N_Vector,__x),convert(N_Vector,__w))
end

function __N_VL1Norm_Serial(x::N_Vector)
    ccall((:N_VL1Norm_Serial,libsundials_nvecserial),realtype,(N_Vector,),x)
end

function N_VL1Norm_Serial(x)
    __x = convert(NVector,x)
    __N_VL1Norm_Serial(convert(N_Vector,__x))
end

function __N_VCompare_Serial(c::realtype,x::N_Vector,z::N_Vector)
    ccall((:N_VCompare_Serial,libsundials_nvecserial),Void,(realtype,N_Vector,N_Vector),c,x,z)
end

function N_VCompare_Serial(c,x,z)
    __x = convert(NVector,x)
    __z = convert(NVector,z)
    __N_VCompare_Serial(c,convert(N_Vector,__x),convert(N_Vector,__z))
end

function __N_VInvTest_Serial(x::N_Vector,z::N_Vector)
    ccall((:N_VInvTest_Serial,libsundials_nvecserial),Cint,(N_Vector,N_Vector),x,z)
end

function N_VInvTest_Serial(x,z)
    __x = convert(NVector,x)
    __z = convert(NVector,z)
    __N_VInvTest_Serial(convert(N_Vector,__x),convert(N_Vector,__z))
end

function __N_VConstrMask_Serial(c::N_Vector,x::N_Vector,m::N_Vector)
    ccall((:N_VConstrMask_Serial,libsundials_nvecserial),Cint,(N_Vector,N_Vector,N_Vector),c,x,m)
end

function N_VConstrMask_Serial(c,x,m)
    __c = convert(NVector,c)
    __x = convert(NVector,x)
    __m = convert(NVector,m)
    __N_VConstrMask_Serial(convert(N_Vector,__c),convert(N_Vector,__x),convert(N_Vector,__m))
end

function __N_VMinQuotient_Serial(num::N_Vector,denom::N_Vector)
    ccall((:N_VMinQuotient_Serial,libsundials_nvecserial),realtype,(N_Vector,N_Vector),num,denom)
end

function N_VMinQuotient_Serial(num,denom)
    __num = convert(NVector,num)
    __denom = convert(NVector,denom)
    __N_VMinQuotient_Serial(convert(N_Vector,__num),convert(N_Vector,__denom))
end


text:
function ConstructionHeuristic(order::Array{NamedTuple{(:j, :prob),Tuple{Int, Float64}}, 1}, crit::String, tasks_by_w::Dict{Int, LTask}, prec::Dict{Int, Array}, bj::Array{Int, 1}, LS::LoadingSequence, TIME::Timer, QC::Array{QuayCrane, 1}, CTS::Constants)
    for tuple in order
        target_bay = tuple.j
        #select a task within the target_bay
        for p in bay_to_pos(target_bay, bj)
            if check_prec(LS, tasks_by_w[p], prec) == true
                if check_loaded_filled(LS, tasks_by_w[p]) == true
                    target_task = tasks_by_w[p]
                    #check if any crane that can perform for this task
                    available_cranes=Array{Tuple{Int, Number}, 1}()
                    for q in TIME.available_cranes
                        #randomly select the cranes based on travel time
                        if target_task.b in QC[q].available_bays
                            if (target_task.b-QC[q].current_bay) == 0
                                available_cranes = [(q, 1)]
                                break
                            else
                                if crit == "dist"
                                    push!(available_cranes, (q, (abs(target_task.b-QC[q].current_bay)^2)/CTS.J*rand()))
                                else
                                    push!(available_cranes, (q, (abs(target_task.b-init_bay(q, CTS)^2)*rand())))
                                end
                                #push!(available_cranes, (q, (abs(target_task.b-QC[q].current_bay)^2)/CTS.J))
                                #push!(available_cranes, (q, (abs(target_task.b-init_bay(q, CTS)^2))))
                            end
                        end
                    end
                    if length(available_cranes) == 0
                        break
                    end

                    #try to add task if there are no clearance problems
                    for q in sort(available_cranes, by=last)
                        if check_clearance(QC, q[1], target_bay, CTS) == true && length(QC[q[1]].task_buffer) == 0
                            add_task(target_task, q[1], TIME, LS, QC, CTS)
                            return("New task successfuly added. \n---------------")
                        end
                    end

                    for q in available_cranes
                        (check, move_q, move_bay, move_status) = check_double_clearance(QC, q[1], target_bay, CTS)
                        if  check == true && length(QC[q[1]].task_buffer) == 0
                            add_task_move(target_task, q[1], move_q, move_bay, move_status, LS, TIME, QC, CTS)
                            return("New task + neighbour move successfuly added. \n---------------")
                        end
                    end

                    for q in available_cranes
                        (check, move_q, move_bay, move_q2, move_bay2) = check_triple_clearance(QC, q[1], target_bay, CTS)
                        if  check == true && length(QC[q[1]].task_buffer) == 0
                            add_move(move_q2, move_bay2, LS, TIME, QC, CTS)
                            add_task_move(target_task, q[1], move_q, move_bay, "idle", LS, TIME, QC, CTS)
                            return("New task + neighbour move successfuly added. \n---------------")
                        end
                    end
                end
            end
        end
    end
    for q in QC
        if q.status != "idle"
            return("Next time period")
        end
    end
    if CTS.Q >= 4
        stucked_bays(order, tasks_by_w, prec, bj, LS, TIME, QC, CTS)
    end
    return("Next time period")
end

function randomizedConstructionHeuristic(ind::String, crit::String, tasks_by_w::Dict{Int, LTask}, prec::Dict{Int, Array}, bj::Array{Int, 1}, LS::LoadingSequence, TIME::Timer, QC::Array{QuayCrane, 1}, CTS::Constants)
    #check if all cranes are busy
    if length(TIME.available_cranes) == 0
        return("All cranes are busy")
    end
    #check if LS is completed
    if LS.tasks_left == 0
        return("LS is completed")
    end

    #order according to indicator
    work_load = order_by_maximal_bay(false, tasks_by_w, bj, LS, CTS)
    if ind == "dist"
        order = order_by_dist(true, work_load, QC, CTS)
    elseif ind == "number"
        work_load_dist = order_by_dist(false, work_load, QC, CTS)
        order = order_by_number(true, work_load_dist, QC, CTS)
    end

    output = ConstructionHeuristic(order, crit, tasks_by_w, prec, bj, LS, TIME, QC, CTS)
    return(output)
end

function deterministicConstructionHeuristic(ind::String, crit::String, tasks_by_w::Dict{Int, LTask}, prec::Dict{Int, Array}, bj::Array{Int, 1}, LS::LoadingSequence, TIME::Timer, QC::Array{QuayCrane, 1}, CTS::Constants)
    #check if all cranes are busy
    if length(TIME.available_cranes) == 0
        return("All cranes are busy")
    end
    #check if LS is completed
    if LS.tasks_left == 0
        return("LS is completed")
    end

    #order according to indicator
    work_load = order_by_maximal_bay(false, tasks_by_w, bj, LS, CTS)
    if ind == "minimal"
        order = reverse!(work_load)
    elseif ind == "number"
        work_load_dist = order_by_dist(false, work_load, QC, CTS)
        order = order_by_number(false, work_load_dist, QC, CTS)
    end

    output = ConstructionHeuristic(order, crit, tasks_by_w, prec, bj, LS, TIME, QC, CTS)
    return(output)
end

function stucked_bays(order::Array{NamedTuple{(:j, :prob),Tuple{Int, Float64}}, 1}, tasks_by_w::Dict{Int, LTask}, prec::Dict{Int, Array}, bj::Array{Int, 1}, LS::LoadingSequence, TIME::Timer, QC::Array{QuayCrane, 1}, CTS::Constants)
    if LS.tasks_left != 0
        for tuple in order
            target_bay = tuple.j
            #select a task within the target_bay
            for p in bay_to_pos(target_bay, bj)
                if check_prec(LS, tasks_by_w[p], prec) == true
                    if check_loaded_filled(LS, tasks_by_w[p]) == true
                        target_task = tasks_by_w[p]
                        #check if any crane that can perform for this task
                        if target_bay in QC[1].available_bays && 1 in TIME.available_cranes
                            add_move(3, target_bay + CTS.delta + 3, LS, TIME, QC, CTS)
                            add_move(4, target_bay + CTS.delta + 5, LS, TIME, QC, CTS)
                            add_task_move(target_task, 1, 2, target_bay + CTS.delta + 1, "idle", LS, TIME, QC, CTS)
                            return("New task successfuly added. \n---------------")
                        elseif target_bay in QC[CTS.Q].available_bays && CTS.Q in TIME.available_cranes
                            add_move(2, target_bay - CTS.delta - 3, LS, TIME, QC, CTS)
                            add_move(1, target_bay - CTS.delta - 5, LS, TIME, QC, CTS)
                            add_task_move(target_task, 4, 3, target_bay - CTS.delta - 1, "idle", LS, TIME, QC, CTS)
                            return("New task successfuly added. \n---------------")
                        end
                    end
                end
            end
        end
        return("Next time period")
    else
        return("LS is completed")
    end
end

max_stars_repo_licenses: Zlib MIT

max_issues_repo_licenses: Zlib MIT

max_forks_repo_licenses: Zlib MIT


text:
mis_animales_favoritos = ("pingüino", "gato", "petauro_del_azúcar")

myfavoriteanimals[1]

mis_animales_favoritos[1] = "nutria"

miagenda = Dict("Jenny" => "867-5309", "Cazafantasmas" => "555-2368")

miagenda["Jenny"]

miagenda["Kramer"] = "555-FILK"

miagenda

pop!(miagenda, "Kramer")

miagenda

miagenda[1]

myfriends = ["Ted", "Robyn", "Barney", "Lily", "Marshall"]

fibonacci = [1, 1, 2, 3, 5, 8, 13]

mezcla = [1, 1, 2, 3, "Ted", "Robyn"]

myfriends[3]

myfriends[3] = "Baby Bop"

push!(fibonacci, 21)

pop!(fibonacci)

fibonacci

favoritos = [["koobideh", "chocolate", "eggs"],["penguins", "cats", "sugargliders"]]

numbers = [[1, 2, 3], [4, 5], [6, 7, 8, 9]]

rand(4, 3)

rand(4, 3, 2)

fibonacci

somenumbers = fibonacci

somenumbers[1] = 404

fibonacci

# Primero restauramos a fibonnaci
fibonacci[1] = 1
fibonacci

somemorenumbers = copy(fibonacci)

somemorenumbers[1] = 404

fibonacci

max_stars_repo_path: es-es/intro-to-julia-ES/plutos/3. Estructuras de datos.jl

max_stars_repo_name: miguelraz/JuliaTutorials

max_issues_repo_name: miguelraz/JuliaTutorials

max_forks_repo_name: miguelraz/JuliaTutorials


text:
# Test extracting way, building, and feature classes from OSM data
module TestClasses

using OpenStreetMap
using Base.Test

MAP_FILENAME = "tech_square.osm"

# Load and crop map to file bounds
nodes, hwys, builds, feats = getOSMData(MAP_FILENAME)
bounds = getBounds(parseMapXML(MAP_FILENAME))
cropMap!(nodes, bounds, highways=hwys, buildings=builds, features=feats, delete_nodes=true)

roads = roadways(hwys)
peds = walkways(hwys)
cycles = cycleways(hwys)
bldg_classes = classify(builds)
feat_classes = classify(feats)

@test length(roads) == 38
@test length(peds) == 55
@test length(cycles) == 36
@test length(bldg_classes) == 17
@test length(feat_classes) == 4

for key in keys(feat_classes)
    @test feat_classes[key] == 1
end

end # module TestClasses


text:
function logunif(min::Float64,max::Float64, n::Int64 = 1)
    exps = (max - min).*rand(n) .+ min
    r = 10.0 .^ exps
    return r
end

function Sample_MoBlueParameters(total_mass::Int64, total_time::Float64, n_sims::Int64)

    @everywhere include("Gillespie.jl")
    sleep_time = logunif(-2.0, 0.0, n_sims)
    k_f =1E-2 #logunif(-6.0, -1.0)
    stable_backward = 1.0#logunif(-6.0,0.0)
    dimerization_ratio = logunif(-5.0,2.0, n_sims)
    mo36_enhance = 1.0 #logunif(0.0, 5.0)
    ball_growth_multiplier =1000.0 #logunif(0.0, 6.0)
    wheel_growth_multiplier =1000.0

    @sync @distributed for i=1:n_sims
        sleep(sleep_time[i])
        Run_MoBlue(total_mass, total_time, k_f, stable_backward, dimerization_ratio[i], mo36_enhance, wheel_growth_multiplier, ball_growth_multiplier)
    end

end


text:
using Base: _maxlength, Bool
using FileIO

# run(`cd GameOfLife`)
files = readdir()
txt_files = filter(endswith(".cells"), files)
for (i, file) in enumerate(txt_files)
    # println(typeof(file))
    filename = split(file, ".")[1]
    cp(file, "$(filename).txt", force = true)
end

content = []
open("input.txt") do f
    line = 0
    global content = []
    while !eof(f)
        s = readline(f)
        if length(s) != 0 && s[1] == '!'
            continue
        else
            line += 1
            println("$s")
        end
        if length(s) >= maxlength
            global maxlength = length(s)
        end
        push!(content,s)
    end

end
# close(f)

finalresult = BitArray[]
for i in 1:length(content)
    global finalresult;
    local1 = []
    if length(content[i])<maxlength
        content[i]=rpad.(content[i],maxlength,"0")
    end
    for j in 1:maxlength
        if content[i][j] == 'O'
            push!(local1,1)
        else
            push!(local1,0)
        end
    end
    push!(finalresult, local1)
end
as = finalresult
println(as)
final = open("result.js", "w")
array = "$as"
# writing to a file using write() method
write(final, "var output= $(array[9:end])")
close(final)

max_stars_repo_name: ashwani-rathee/Sigil

max_issues_repo_name: ashwani-rathee/Sigil

max_forks_repo_name: ashwani-rathee/Sigil


text:
export evaluate, gradient, gradient!,
    evaluate_and_gradient, evaluate_and_gradient!,
    differentiate_parameters, differentiate_parameters!,
    hessian!, hessian, gradient_and_hessian, gradient_and_hessian!

import Base: @propagate_inbounds

@doc """
    evaluate(f::Polynomial, x)

Evaluate the polynomial `f` at `x`.
""" evaluate(f::Polynomial, x)

@doc """
    evaluate(f::Polynomial, x, p)

Evaluate `f` at `x` with parameters `p`.
"""  evaluate(f::Polynomial, x, p)

(f::Polynomial)(x::AbstractVector) = evaluate(f, x)
(f::Polynomial)(x::AbstractVector, p) = evaluate(f, x, p)

function evaluate_impl(f::Type{Polynomial{T, E, P}}) where {T, E, P}
    if P == Nothing
        access_input = x_
    else
        n = size(P, 1)
        access_input(i) = i ≤ n ? :(p[$i]) : :(x[$(i-n)])
    end

    boundschecks = [(:@boundscheck length(x) ≥ $(size(E)[1]))]
    if P != Nothing
        push!(boundschecks, (:@boundscheck length(p) ≥ $(size(P)[1])))
    end

    quote
        $(boundschecks...)
        c = coefficients(f)
        out = @inbounds begin
            $(generate_evaluate(exponents(P, E), T, access_input))
        end
        out
    end
end

@generated function evaluate(f::Polynomial{T, E, Nothing}, x::AbstractVector) where {T, E}
    evaluate_impl(f)
end
@generated function evaluate(f::Polynomial{T, E, P}, x::AbstractVector, p) where {T, E, P}
    evaluate_impl(f)
end

function _val_gradient_impl(f::Type{Polynomial{T, E, P}}) where {T, E, P}
    if P == Nothing
        access_input = x_
    else
        n = size(P, 1)
        access_input(i) = i ≤ n ? :(p[$i]) : :(x[$(i-n)])
    end
    boundschecks = [(:@boundscheck length(x) ≥ $(size(E)[1]))]
    if P != Nothing
        push!(boundschecks, :(@boundscheck length(p) ≥ $(size(P)[1])))
    end

    quote
        $(boundschecks...)
        c = coefficients(f)
        val, grad = @inbounds begin
            $(generate_gradient(exponents(E), exponents(P), T, access_input))
        end
        val, grad
    end
end

@generated function _val_gradient(f::Polynomial{T, E, Nothing}, x::AbstractVector) where {T, E}
    _val_gradient_impl(f)
end

@generated function _val_gradient(f::Polynomial, x::AbstractVector, p)
    _val_gradient_impl(f)
end

@propagate_inbounds @inline function _gradient(f, x)
    _, grad = _val_gradient(f, x)
    grad
end

@propagate_inbounds @inline function _gradient(f, x, p)
    _, grad = _val_gradient(f, x, p)
    grad
end

@doc """
    gradient(f::Polynomial, x)

Evaluate the gradient of the polynomial `f` at `x`.
"""
function gradient(f::Polynomial{T, E, Nothing}, x::AbstractVector) where {T, E}
    Vector(_gradient(f, x))
end
function gradient(f::Polynomial{T, E, Nothing}, x::SVector) where {T, E}
    _gradient(f, x)
end

"""
    gradient(f::Polynomial, x, p)

Evaluate the gradient of the polynomial `f` at `x` with parameters `p`.
"""
gradient(f::Polynomial, x::AbstractVector, p) = Vector(_gradient(f, x, p))
gradient(f::Polynomial, x::SVector, p) = _gradient(f, x, p)

"""
    gradient!(u, f::Polynomial, x)

Evaluate the gradient of the polynomial `f` at `x` and store the result in `u`.

    gradient!(u, f::Polynomial, x, p)

Evaluate the gradient of the polynomial `f` at `x` with parameters `p` and store the result in `u`.
"""
@propagate_inbounds function gradient!(u::AbstractVector, f::Polynomial, x::AbstractVector)
    u .= _gradient(f, x)
    u
end
@propagate_inbounds function gradient!(u::AbstractVector, f::Polynomial, x::AbstractVector, p)
    u .= _gradient(f, x, p)
    u
end

"""
    evaluate_and_gradient(f::Polynomial, x)

Evaluate the polynomial `f` and its gradient at `x`. Returns a tuple.

    evaluate_and_gradient(f::Polynomial, x, p)

Evaluate the polynomial `f` and its gradient at `x` with parameters `p`. Returns a tuple.
"""
function evaluate_and_gradient(f::Polynomial{T, E, Nothing}, x::SVector) where {T, E}
    _val_gradient(f, x)
end
evaluate_and_gradient(f::Polynomial, x::SVector, p) = _val_gradient(f, x, p)

function evaluate_and_gradient(f::Polynomial{T, E, Nothing}, x::AbstractVector) where {T, E}
    val, grad = _val_gradient(f, x)
    val, Vector(grad)
end
function evaluate_and_gradient(f::Polynomial, x::AbstractVector, p)
    val, grad = _val_gradient(f, x, p)
    val, Vector(grad)
end

"""
    evaluate_and_gradient!(u, f::Polynomial, x)

Evaluate the polynomial `f` and its gradient at `x`. Stores the gradient in `u` and
returns the `f(x)`.

    evaluate_and_gradient!(u, f::Polynomial, x, p)

Evaluate the polynomial `f` and its gradient at `x` with parameters `p`. Stores the gradient in `u` and
returns the `f(x)`.
"""
@propagate_inbounds function evaluate_and_gradient!(u::AbstractVector, f::Polynomial{T, E, Nothing}, x::AbstractVector) where {T, E}
    val, grad = _val_gradient(f, x)
    u .= grad
    val
end
@propagate_inbounds function evaluate_and_gradient!(u::AbstractVector, f::Polynomial, x::AbstractVector, p)
    val, grad = _val_gradient(f, x, p)
    u .= grad
    val
end

#############
## Hessian ##
#############

function _gradient_hessian_impl(f::Type{Polynomial{T, E, P}}; return_grad=true) where {T, E, P}
    if P == Nothing
        access_input = x_
    else
        n = size(P, 1)
        access_input(i) = i ≤ n ? :(p[$i]) : :(x[$(i-n)])
    end

    exps = exponents(E)
    params = exponents(P)

    n, m = size(E)

    code = []
    ∇ = [Symbol("∇", i) for i=1:n]
    vals = [Symbol("val", i) for i=1:n]

    # for each variable we simulate the derivative
    for i in 1:n
        coeffsᵢ = Union{Symbol, Expr}[]
        expsᵢ = Vector{eltype(exps)}()
        paramsᵢ = P === nothing ? nothing : Vector{eltype(exps)}()
        mᵢ = 0

        # For each term make the symbolic derivative
        for j in 1:m
            if exps[i,j] == 0
                continue
            end
            for k in 1:n
                if i == k
                    push!(expsᵢ, exps[k, j] - 1)
                    if exps[i,j] > 1
                        push!(coeffsᵢ, :($(exps[i,j]) * c[$j]))
                    else
                        push!(coeffsᵢ, :(c[$j]))
                    end
                else
                    push!(expsᵢ, exps[k, j])
                end
            end
            if params !== nothing
                for k in 1:size(params, 1)
                    push!(paramsᵢ, params[k, j])
                end
            end

            mᵢ += 1
        end
        if mᵢ == 0
            zero_tuple = Expr(:tuple, (:(zero($T)) for _=1:n)...)
            val_grad_codeᵢ = :((zero($T), SVector($zero_tuple)))
        else
            Eᵢ = reshape(expsᵢ, (n, mᵢ))
            Pᵢ = P === nothing ? nothing : reshape(paramsᵢ, (size(paramsᵢ, 1), mᵢ))
            val_grad_codeᵢ = generate_gradient(Eᵢ, Pᵢ, T, access_input, coeffsᵢ)
        end

        push!(code, :($(Expr(:tuple, vals[i], ∇[i])) = begin $(val_grad_codeᵢ) end))
    end

    push!(code, :(grad = SVector($(Expr(:tuple, vals...)))))
    push!(code, :(hess = assemble_matrix(SVector($(Expr(:tuple, (∇[i] for i=1:n)...))))))
    quote
        @boundscheck length(x) ≥ $(size(E)[1])
        c = coefficients(f)
        $(code...)
        grad, hess
    end
end

@generated function gradient_and_hessian(f::Polynomial{T,E,Nothing}, x::AbstractVector) where {T, E}
    _gradient_hessian_impl(f)
end
@generated function gradient_and_hessian(f::Polynomial, x::AbstractVector, p)
    _gradient_hessian_impl(f)
end

@doc """
     gradient_and_hessian(F::PolynomialSystem, x, [p])

Compute the gradient and Hessian of `f` at `x`.
""" gradient_and_hessian(F::Polynomial, x)

"""
    gradient_and_hessian!(u, U, f::Polynomial, x, [p])

Compute the gradient and Hessian of `f` at `x` and store them in `u` and `U`.
"""
function gradient_and_hessian!(u, U, f::Polynomial, xp...)
    grad, hess = gradient_and_hessian(f, xp...)
    u .= grad
    U .= hess
    nothing
end

"""
    hessian(f::Polynomial, x, [p])

Compute the hessian of `f` at `x`.
"""
function hessian(f::Polynomial, xp...)
    last(gradient_and_hessian(f, xp...))
end

"""
    hessian!(U, f::Polynomial, x, [p])

Compute the hessian of `f` at `x` and store it in `U`.
"""
function hessian!(U, f::Polynomial, xp...)
    U .= hessian(f, xp...)
    U
end

################
## Parameters ##
################
function _differentiate_parameters_impl(f::Type{Polynomial{T, E, P}}) where {T, E, P}
    @assert P != Nothing

    # The role of E and P is interchanged
    n = size(E, 1)
    access_input(i) = i ≤ n ? :(x[$i]) : :(p[$(i-n)])

    quote
        @boundscheck length(x) ≥ $(size(E, 1))
        @boundscheck length(p) ≥ $(size(P, 1))
        c = coefficients(f)
        @inbounds begin
            $(generate_differentiate_parameters(exponents(E), exponents(P), T, access_input))
        end
    end
end

@generated function _differentiate_parameters(f::Polynomial, x::AbstractVector, p)
    _differentiate_parameters_impl(f)
end

@doc """
    differentiate_parameters(f::Polynomial, x, p)

Evaluate the gradient of the polynomial `f` w.r.t. the parameters at `x` with parameters `p`.
"""
differentiate_parameters(f::Polynomial, x, p) = Vector(_differentiate_parameters(f, x, p))
differentiate_parameters(f::Polynomial, x::SVector, p) = _differentiate_parameters(f, x, p)

@doc """
    differentiate_parameters!(u, f::Polynomial, x, p)

Evaluate the gradient of the polynomial `f` w.r.t. the parameters at `x` with parameters `p` and store the result in `u`.
"""
function differentiate_parameters!(u, f::Polynomial, x, p)
    u .= _differentiate_parameters(f, x, p)
    u
end