text:
using Jenks
using Test

@testset "RandnLin" begin
    data = randn(10_000)
    n = 10
    JR = JenksClassification(n,data,feedback=false)
    @test minimum(data) < JR.centres[1]
    @test maximum(data) > JR.centres[end]
    @test minimum(data) == JR.bounds[1]
    @test maximum(data) == JR.bounds[end]
    @test JR.breaks[1] == 1
    @test JR.breaks[end] == length(data)+1
    @test JR.ARE > 0.0
    @test JR.GVF >= 0.0 && JR.GVF <= 1.0
    @test sum(JR.n_in_class) == length(data)
end

@testset "RandnSq" begin
    data = randn(10_000)
    n = 10
    JR = JenksClassification(n,data,errornorm=2,feedback=false)
    @test minimum(data) < JR.centres[1]
    @test maximum(data) > JR.centres[end]
    @test minimum(data) == JR.bounds[1]
    @test maximum(data) == JR.bounds[end]
    @test JR.breaks[1] == 1
    @test JR.breaks[end] == length(data)+1
    @test JR.ARE > 0.0
    @test JR.GVF >= 0.0 && JR.GVF <= 1.0
    @test sum(JR.n_in_class) == length(data)
end

@testset "RandnSq" begin
    data = randn(10_000)
    n = 10
    JR = JenksClassification(n,data,errornorm=1,feedback=false,inimode="rand")
    @test minimum(data) < JR.centres[1]
    @test maximum(data) > JR.centres[end]
    @test minimum(data) == JR.bounds[1]
    @test maximum(data) == JR.bounds[end]
    @test JR.breaks[1] == 1
    @test JR.breaks[end] == length(data)+1
    @test JR.ARE > 0.0
    @test JR.GVF >= 0.0 && JR.GVF <= 1.0
    @test sum(JR.n_in_class) == length(data)
end

ext: jl


text:
"""
`ExponentiatedKernel([ρ=1])`
The exponentiated kernel is a Mercer kernel given by:
```
    κ(x,y) = exp(ρ²xᵀy)
```
"""
struct ExponentiatedKernel{T,Tr} <: Kernel{T,Tr}
    transform::Tr
    metric::DotProduct
    function ExponentiatedKernel{T,Tr}(transform::Tr) where {T,Tr<:Transform}
        return new{T,Tr}(transform,DotProduct())
    end
end
@inline kappa(κ::ExponentiatedKernel, xᵀy::T) where {T<:Real} = exp(xᵀy)


text:
# ----------------------------
# Advection-Diffusion Example
# ----------------------------

using KitBase, Plots

set = Setup(
    "scalar", # matter
    "advection", # case
    "1d0f0v", # space
    "gks", # flux
    "", # collision: for scalar conservation laws there are none
    1, # species
    1, # interpolation order
    "vanleer", # limiter
    "period", # boundary
    0.5, # cfl
    1.0, # simulation time
)
pSpace = PSpace1D(0.0, 1.0, 100, 1)
vSpace = nothing
property = Scalar(1.0, 1e-6)
ib = IB(x -> sin(2π * x), property)

ks = SolverSet(set, pSpace, vSpace, property, ib)
ctr, face = init_fvm(ks, ks.ps)

t = 0.0
dt = KitBase.timestep(ks, ctr, t)
nt = ks.set.maxTime ÷ dt |> Int
anim = @animate for iter = 1:nt
    reconstruct!(ks, ctr)
    evolve!(ks, ctr, face, dt)
    update!(ks, ctr, face, dt, 0.0)

    plot(ks, ctr, ylims = [-1, 1])
end

gif(anim, "advection.gif", fps = 45)


text:
using UUIDs, Dates

include("filtering.jl")

# big numbers may fail to display properly when deserialized in JS
_is_javascript_safe(x::Real) = false

function _is_javascript_safe(x::Integer)
    min_safe_int = -(Int64(2)^53-1)
    max_safe_int = Int64(2)^53-1
    min_safe_int < x < max_safe_int
end

function _is_javascript_safe(x::AbstractFloat)
    min_safe_float = -(Float64(2)^53-1)
    max_safe_float = Float64(2)^53-1
    min_safe_float < x < max_safe_float
end

# ag-grid special cases `.` in fieldnames; we need to special case that here
col_name_fixer(name) = replace(string(name), '.' => '_')

# loading DataValues will add an overload for this
json_sprint(x) = sprint(print, x)

schema_type(::Type{Union{}}) = "string"
schema_type(::Type{T}) where {T} = "string"
schema_type(::Type{T}) where {T <: AbstractFloat} = "number"
schema_type(::Type{T}) where {T <: Integer} = "integer"
schema_type(::Type{T}) where {T <: Bool} = "boolean"
schema_type(::Type{T}) where {T <: Dates.Time} = "time"
schema_type(::Type{T}) where {T <: Dates.Date} = "date"
schema_type(::Type{T}) where {T <: Dates.DateTime} = "datetime"
schema_type(::Type{T}) where {T <: AbstractString} = "string"

ag_schema_type(::Type{T}) where {T} = nothing
ag_schema_type(::Type{T}) where {T <: Number} = "numericColumn"

ag_filter_type(::Type{T}) where {T}= true
ag_filter_type(::Type{T}) where {T <: Number} = "agNumberColumnFilter"
ag_filter_type(::Type{T}) where {T <: Union{Dates.Date, Dates.DateTime}} = "agDateColumnFilter"

# for small tables only
function print_table(io::IO, source, col_names, fixed_col_names, col_types, title = "")
    ctx = JSON.Writer.CompactContext(io)

    JSON.begin_object(ctx)

    JSON.show_pair(ctx, JSON.StandardSerialization(), "name", title)

    JSON.show_key(ctx, "schema")
    print_schema(ctx, col_names, fixed_col_names, col_types, filterable = true)

    JSON.show_key(ctx, "data")
    print_body(ctx, source, fixed_col_names)

    JSON.end_object(ctx)
end

# https://specs.frictionlessdata.io/table-schema
function print_schema(ctx, col_names, fixed_col_names, col_types; filterable = false, sortable = false)
    ser = JSON.StandardSerialization()

    JSON.begin_object(ctx)
    JSON.show_key(ctx, "fields")
    JSON.begin_array(ctx)
    for i in 1:length(col_names)
        JSON.delimit(ctx)
        JSON.begin_object(ctx)
        # standard fields
        JSON.show_pair(ctx, ser, "name", fixed_col_names[i])
        JSON.show_pair(ctx, ser, "title", col_names[i])
        JSON.show_pair(ctx, ser, "type", schema_type(col_types[i]))
        # custom fields
        JSON.show_pair(ctx, ser, "ag_type", ag_schema_type(col_types[i]))
        JSON.show_pair(ctx, ser, "ag_filter", filterable ? ag_filter_type(col_types[i]) : false)
        JSON.show_pair(ctx, ser, "ag_sortable", sortable)
        JSON.end_object(ctx)
    end
    JSON.end_array(ctx)
    JSON.end_object(ctx)
end

function print_body(ctx, source, fixed_col_names; first = 1, last = typemax(Int64))
    JSON.begin_array(ctx)
    i = 0
    for row in source
        i += 1
        i < first && continue
        i > last && break

        JSON.delimit(ctx)
        JSON.begin_object(ctx)
        for col = 1:length(fixed_col_names)
            print_el(ctx, fixed_col_names[col], row[col])
        end
        JSON.end_object(ctx)
    end
    JSON.end_array(ctx)
end

function print_el(ctx, name, val)
    ser = JSON.StandardSerialization()
    if val isa Real && isfinite(val) && _is_javascript_safe(val)
        JSON.show_pair(ctx, ser, name, val)
    elseif val === nothing || val === missing
        JSON.show_pair(ctx, ser, name, repr(val))
    else
        JSON.show_pair(ctx, ser, name, json_sprint(val))
    end
end

_isiterabletable = x -> false
_getiterator = x -> x

const tabletraits_uuid = UUIDs.UUID("3783bdb8-4a98-5b6b-af9a-565f29a5fe9c")
const datavalues_uuid = UUIDs.UUID("e7dc6d0d-1eca-5fa6-8ad6-5aecde8b7ea5")
function on_pkg_load(pkg)
    if  pkg.uuid == tabletraits_uuid
        TableTraits = Base.require(pkg)

        global _isiterabletable = TableTraits.isiterabletable
        global _getiterator = TableTraits.getiterator
    elseif pkg.uuid == datavalues_uuid
        DataValues = Base.require(pkg)

        eval(
            quote
                function json_sprint(val::$(DataValues.DataValue))
                    $(DataValues.isna)(val) ? "null" : json_sprint(val[])
                end
            end
        )
    end
end

# these assume that table elements are reasonably small
const MAX_SYNC_TABLE_ELEMENTS = 100_000
const MAX_CACHE_TABLE_ELEMENTS = 10_000_000

# make a copy of medium size tables (MAX_SYNC_TABLE_ELEMENTS < #el < MAX_CACHE_TABLE_ELEMENTS)
# store a reference for big tables
# (column_names, column_types, table_iterator, table_length, table_indexable)
const TABLES = Dict{UUID, Tuple{Any, Any, Any, Int, Bool}}()

# special-case vectors for a few known eltypes
showtable(table::AbstractVector{<:Union{Number, AbstractString, Date, DateTime, Time, AbstractVector}}, title = "") = showtable(reshape(table, :, 1), title)

function showtable(table::T, title = "") where T
    if showable("application/vnd.dataresource+json", table)
        return _display(InlineDisplay(), table)
    end

    iter =  _getiterator(table)

    if T <: AbstractMatrix
        col_names = [string(i) for i in 1:size(table, 2)]
        col_types = [eltype(table) for _ in 1:size(table, 2)]
        # transform matrix to iterator over its rows
        iter = (view(table, i, :) for i in axes(table, 1)) # eachrow
    elseif Base.IteratorEltype(iter) isa Base.EltypeUnknown
        first_el = first(iter)
        col_names = String.(propertynames(first_el))
        col_types = [fieldtype(typeof(first_el), i) for i = 1:length(col_names)]
    else
        col_names = String.(fieldnames(eltype(iter)))
        col_types = [fieldtype(eltype(iter), i) for i = 1:length(col_names)]
    end
    fixed_col_names = col_name_fixer.(col_names)

    if length(fixed_col_names) == 0
        throw(ErrorException("Input table does not seem to have columns."))
    end

    if Base.haslength(iter)
        tablelength = length(iter)
        els = tablelength * length(col_names)
        async = els > MAX_SYNC_TABLE_ELEMENTS
        should_copy = els < MAX_CACHE_TABLE_ELEMENTS
        indexable = els < MAX_CACHE_TABLE_ELEMENTS
    else
        tablelength = 0
        async = true
        should_copy = false
        indexable = false
    end

    if async
        id = uuid4()
        if should_copy
            TABLES[id] = (col_names, fixed_col_names, collect(iter), tablelength, indexable)
        else
            TABLES[id] = (col_names, fixed_col_names, iter, tablelength, indexable)
        end

        io = IOBuffer()
        ctx = JSON.Writer.CompactContext(io)
        print_schema(ctx, col_names, fixed_col_names, col_types, sortable = should_copy, filterable = should_copy)

        schema = JSON.JSONText(String(take!(io)))

        payload = (
            schema = schema,
            rowCount = tablelength,
            name = title,
            id = string(id)
        )
        sendDisplayMsg("application/vnd.dataresource+lazy", JSON.json(payload))
    else
        io = IOBuffer()
        print_table(io, iter, col_names, fixed_col_names, col_types, title)
        sendDisplayMsg("application/vnd.dataresource+json", String(take!(io)))
    end
end

function get_table_data(conn, params::GetTableDataRequest)
    id = UUID(params.id)
    if !haskey(TABLES, id)
        return JSONRPC.JSONRPCError(-32600, "Table not found.", nothing)
    end

    col_names, fixed_col_names, table, tablelength, indexable = TABLES[id]
    will_filter = !isempty(params.filterModel)
    will_sort = !isempty(params.sortModel)

    # this will only be called for medium-sized tables
    if will_filter
        filter_hash = hash(params.filterModel)
        # we have already filtered this table with the specified filterModel
        if haskey(FILTER_CACHE, id) && first(FILTER_CACHE[id]) == filter_hash
            table = last(FILTER_CACHE[id])
        else
            filt = generate_filterer(params.filterModel, fixed_col_names)
            table = filter(r -> all(f -> f(r), filt), table)
            FILTER_CACHE[id] = (filter_hash, table)
        end
        tablelength = length(table)
    end
    if will_sort
        sort_hash = hash(params.sortModel)
        if get(SORTED_CACHE, id, 0x0) != sort_hash
            sorter = generate_sorter(params.sortModel, fixed_col_names)
            Base.invokelatest(sort!, table, lt = sorter)
            SORTED_CACHE[id] = sort_hash
        end
    end

    io = IOBuffer()
    ctx = JSON.Writer.CompactContext(io)
    if indexable
        part = isempty(table) ? table : @view table[min(params.startRow + 1, end):min(params.endRow + 1, end)]
        Base.invokelatest(print_body, ctx, part, fixed_col_names)
    else
        Base.invokelatest(print_body, ctx, table, fixed_col_names; first = params.startRow + 1, last = params.endRow + 1)
    end

    return (
        rows = JSON.JSONText(String(take!(io))),
        lastRow = tablelength
    )
end

function clear_lazy_table(conn, params::NamedTuple{(:id,), Tuple{String}})
    id = UUID(params.id)
    delete!(TABLES, id)
    delete!(FILTER_CACHE, id)
    delete!(SORTED_CACHE, id)
end

max_stars_repo_name: JaapHaitsma/julia-vscode

max_issues_repo_name: JaapHaitsma/julia-vscode

max_forks_repo_name: JaapHaitsma/julia-vscode


text:
# +
include("epidemiology.jl")

using OrdinaryDiffEq
using Catlab.Graphics.Graphviz: run_graphviz

# +
save_wd(ex, fname::AbstractString) = begin
    g = display_wd(ex)
    open(fname, "w") do io
        run_graphviz(io, g, format="svg")
    end
end
save_graph(g, fname::AbstractString) = begin
    open(fname, "w") do io
        run_graphviz(io, g, format="svg")
    end
end

splot(soln, fname) = begin
    p = plot(soln, linewidth=5, legend=false)
    savefig(p, fname)
end

splot(soln, t, fname) = begin
    p = plot(t, soln, linewidth=5, legend=false)
    savefig(p, fname)
end

splotchannels(sol, dir) = begin
    mkpath(dir)
    jp(p) = joinpath(dir, p)
    splot(sol, jp("solution.svg"))
    pltdims = [3,8,13]
    splot(sol[pltdims,:]', sol.t, jp("infecteds.svg"))
    splot(sum(sol(sol.t, idxs=pltdims), dims=1)', sol.t, jp("infecteds_sum.svg"))
    pltdims = [3,8,13] .+ 1
    splot(sol[pltdims,:]', sol.t, jp( "recovered.svg" ))
    splot(sum(sol(sol.t, idxs=pltdims), dims=1)', sol.t, jp( "recovered_sum.svg" ))
end
# -

save_wd(sir, "sir.svg")
save_graph(Graph(p_sir), "sir_graph.svg")

# define initial states and transition rates
u0 = [10.0, 1, 0]
p = [0.4, 0.4]
# create and solve ODE problem
prob = ODEProblem(p_sir,u0,(0.0,7.5),p)
sol = OrdinaryDiffEq.solve(prob,Tsit5())
# visualize the solution
splot(sol, "sir_soln.svg")

save_wd(seir, "seir.svg")
save_graph(Graph(p_seir), "seir_graph.svg")

# define initial states and transition rates
u0 = [10.0, 1, 0, 0]
p = [0.9, 0.2, 0.5]
# create and solve ODE problem
prob = ODEProblem(p_seir,u0,(0.0,15.0),p)
sol = OrdinaryDiffEq.solve(prob,Tsit5())
# visualize the solution
splot(sol, "seir_soln.svg")

save_wd(seird, "seird.svg")
save_graph(Graph(p_seird), "seird_graph.svg")

# define initial states and transition rates
u0 = [10.0, 1, 0, 0, 0]
p = [0.9, 0.2, 0.5, 0.1]
# create and solve ODE problem
prob = ODEProblem(p_seird,u0,(0.0,15.0),p)
sol = OrdinaryDiffEq.solve(prob,Tsit5())
# visualize the solution
splot(sol, "seird_soln.svg")

include("covid.jl")

[save_graph(Graph(decoration(p)), pname) for (p,pname) in [(spontaneous_petri, "spont.svg"),
                                               (transmission_petri, "transmission.svg"),
                                               (exposure_petri, "exposure.svg"),
                                               (travel_petri, "travel.svg")]]

save_wd(seird_city, "seird_city.svg")
save_graph(Graph(decoration(F(seird_city))), "seird_graph.svg")

save_wd(seird_3, "seird_3.svg")
save_graph(Graph(p_seird_3), "seird_3_graph.svg")

pltdims = [4,9, 14]
# Define time frame, 3 months
tspan = (0.0,90.0)
# Define initial states
u0 = zeros(Float64, base(pc_seird_3).n)
u0[1]  = 10000
u0[6]  = 10000
u0[11] = 10000
u0[2]  = 1
# Define transition rates
params = seirdparams(3, 5)
# Generate and solve resulting ODE
prob = ODEProblem(p_seird_3,u0,tspan,params)
sol = OrdinaryDiffEq.solve(prob,Tsit5())
splotchannels(sol, "seird")

function drawsol(f, a, b, p, dir)
    tspan = (0, 240.0)
    println("INFO: Processing $dir")
    dynparams = waveparams(f, a/sum(u0), b/sum(u0),p)
    mkpath(dir)
    prob = ODEProblem(p_seird_3,u0,tspan, dynparams)
    sol = OrdinaryDiffEq.solve(prob,Tsit5(), saveat=1:1:tspan[2])
    splotchannels(sol, dir)
    savefig(plot(sol.t, f(a,b,p), ylim=(0,12), legend=false, linewidth=5), "$dir/transmissability.svg")
end

a,b = 10, 1
dir = "dynamic/mitigation/"
drawsol(asymptotic,  a,b,1/6, "$dir/asymptotic")
drawsol(triangleasm, a,b,25 , "$dir/trianglesm")
drawsol(coswave,     a,b,25 , "$dir/coswave")
drawsol(sincwave,    a,b,100, "$dir/sincwave")
drawsol(modsincwave, a,b,100, "$dir/modsincwave")

a,b = 10, 1/2
dir = "dynamic/containment"
drawsol(asymptotic,  a,b,1/6, "$dir/asymptotic")
drawsol(triangleasm, a,b, 25, "$dir/trianglesm")
drawsol(coswave,     a,b, 25, "$dir/coswave")
drawsol(sincwave,    a,b,100, "$dir/sincwave")
drawsol(modsincwave, a,b,100, "$dir/modsincwave")


text:
mutable struct SLOB
    num_paths::Int64
    T::Int64
    p₀::Float64
    M::Int64
    L::Float64
    D::Float64
    σ::Float64
    nu::Float64
    α::Float64
    source_term::SourceTerm
    x::Array{Float64, 1}
    Δx::Float64
    Δt::Float64
end
function SLOB(num_paths::Int64, T::Int64, p₀::Float64,
    M::Int64, L::Real, D::Float64, σ::Float64, nu::Float64,
    α::Float64, source_term::SourceTerm)

    x₀ = p₀ - 0.5*L
    x_m = p₀ + 0.5*L
    @assert x₀ >= 0
    x = collect(Float64, range(x₀, stop=x_m, length=M+1))
    Δx = L/M
    Δt = (Δx^2) / (2.0*D)
    return SLOB(num_paths, T, p₀, M, L, D, σ, nu, α, source_term, x, Δx, Δt)
end

SLOB(dict) = SLOB(
    dict["num_paths"], dict["T"], dict["p₀"], dict["M"],
    dict["L"], dict["D"], dict["σ"], dict["nu"], dict["α"],
    SourceTerm(dict["λ"], dict["μ"]))

SLOB(;num_paths::Int64=1, T::Int64=100, p₀::Real=100.0,
    M::Int64=100, L::Real=50.0, D::Real=4.0, σ::Real=1.0,
    nu::Real=0.1, α::Real=20.0, λ::Real=1.0, μ::Real=0.5) =
    SLOB(num_paths, T, p₀, M, L, D, σ, nu, α,
    SourceTerm(λ, μ))

function (slob::SLOB)(seed::Int=-1)
    if seed == -1
            seeds = Int.(rand(MersenneTwister(), UInt32, slob.num_paths))
        else
            seeds = Int.(rand(MersenneTwister(seed), UInt32, slob.num_paths))
    end

    time_steps = get_time_steps(slob.T, slob.Δt)

    mid_price_paths = ones(Float64, slob.T + 1, slob.num_paths)
    for path in 1:slob.num_paths
        Random.seed!(seeds[path])
        mid_price_paths[:, path] = dtrw_solver(slob)
    end

    return mid_price_paths
end

function (slob::SLOB)(debug::Bool, seed::Int=-1)
    logger = FileLogger(Dict(Logging.Info => "info.log", Logging.Error => "error.log"), append=false)
    global oldglobal = global_logger(logger)
    if seed == -1
            seeds = Int.(rand(MersenneTwister(), UInt32, slob.num_paths))
        else
            seeds = Int.(rand(MersenneTwister(seed), UInt32, slob.num_paths))
    end
    time_steps = get_time_steps(slob.T, slob.Δt)

    raw_price_paths = ones(Float64, time_steps + 1, slob.num_paths)

    sample_price_paths = ones(Float64, slob.T + 1, slob.num_paths)

    lob_densities = zeros(Float64, slob.M+1, time_steps + 1, slob.num_paths)
    P⁺s = ones(Float64, time_steps, slob.num_paths)
    P⁻s = ones(Float64, time_steps, slob.num_paths)
    Ps = ones(Float64, time_steps, slob.num_paths)
    for path in 1:slob.num_paths
        Random.seed!(seeds[path])
        @info "path $path with seed $(seeds[path])"
        lob_densities[:, :, path], raw_price_paths[:, path],
            sample_price_paths[:, path], P⁺s[:, path],
            P⁻s[:, path], Ps[:, path] = dtrw_solver(debug, slob)
    end

    return lob_densities, raw_price_paths, sample_price_paths, P⁺s, P⁻s, Ps
end

function get_time_steps(T::Int, Δt::Float64)
    time_steps =  T / Δt
    return floor(Int, time_steps + eps(time_steps))
end

function sample_mid_price_path(slob::SLOB, price_path)
    Δt = slob.Δt
    get_mid_price_inds(x) = get_time_steps(x, Δt)
    mid_price_inds = get_mid_price_inds.(0:slob.T) .+ 1
    mid_prices = price_path[mid_price_inds]
    return mid_prices
end


text:
# From: https://thenewphalls.wordpress.com/tag/julia/

(errRead, errWrite) = redirect_stderr()
(outRead, outWrite) = redirect_stdout()

function flushlog()
    global errRead
    global errWrite
    global outRead
    global outWrite
    dir = "/home/nlg-05/dy_052/kuparser/dev/test"

    close(errWrite)
    err = readavailable(errRead)
    close(errRead)
    errfile = open("$(dir)/$(myid()).err", "a")
    write(errfile, err)
    close(errfile)

    close(outWrite)
    out = readavailable(outRead)
    close(outRead)
    outfile = open("$(dir)/$(myid()).out", "a")
    write(outfile, out)
    close(outfile)
end

atexit(flushlog)


text:
#
# Nested Monte Carlo Search choice model
#
# constructor usage:
#
#		NCMSChoiceModel(choicemodel, fitnessfn, samplesize)
#
# where:
#	 choicemodel is an instance of another choice model (e.g. SamplerChoiceModel) used for simulating the outcome of choices
#  fitnessfn is a fitness function that takes a generated object as a parameter and returns a fitness where lower values are better
#  samplesize is the number of choices to sample when deciding on the value of each choice point
#
# example:
#		scm = SamplerChoiceModel(gn)
#		f = abs(size(x)-64)
#		ncm = NMCSChoiceModel(scm, f, 2)
#
# To implement higher level NMCS, use an NMCS instance as the policy choice model.  For example:
#		ncm2 = NMCSChoiceModel(NMCSChoiceModel(scm,f,2),f,2)
# specifies a 2-level NMCS
# While a little verbose, this can allow different sample sizes at different levels (and even different fitness functions).
#
# Note that NMCSChoiceModel does not require a fitness function that handles the termination exception (or, equivalently
# when using robustgen, nothing as the value returned by the generator).  Instead this is handled internally by the
# choice model.  However, the generator may still raise this exception (or return nothing when using robustgen) if
# all simulations at a particular choice points are terminated by the exception
#

type NMCSChoiceModel <: ChoiceModel
	policychoicemodel::ChoiceModel
	fitnessfunction::Function
	samplesize::Int 									# the (maximum) number of samples to take at each choice point
	minimumsamplesize::Int								# the minimum sample size (comes into effect when totalsamplelimit or ruledepthlimit is reached)
	totalsamplelimit::Int								# limit the total number of samples to take (does not include samples taken by nested samplers)
	ruledepthlimit::Int									# limit on generator rule depth
	totalsamplecount::Int								# number of samples taken in total so far
	bestfitness::Real 									# lower is better
	bestgodelsequence::Vector{Real} 					# the best godel sequence found so far
	besttracesequence::Vector{Dict}						# the traces for that sequence from the underlying choice model
	function NMCSChoiceModel(policychoicemodel::ChoiceModel, fitnessfunction::Function, samplesize::Int=1, minimumsamplesize::Int=0, totalsamplelimit::Int=typemax(Int), ruledepthlimit::Int=typemax(Int))
		(minimumsamplesize >= 0) || error("Minimum sample size must be 0 or more")
		(minimumsamplesize <= samplesize) || error("Minimum sample size must be less or equal to the sample size")
		(totalsamplelimit >= 0) || error("Total sample count must be 0 or more")
		(ruledepthlimit >= 0) || error("Rule depth limit must be 0 or more")
		new(deepcopy(policychoicemodel), fitnessfunction, samplesize, minimumsamplesize, totalsamplelimit, ruledepthlimit, 0, +Inf, Real[], Dict[])
	end
end

function NMCSChoiceModel(g::Generator, fitnessfunction::Function, samplesize::Int=1; minimumsamplesize::Int=0, totalsamplelimit::Int=typemax(Int), ruledepthlimit::Int=typemax(Int))
	NMCSChoiceModel(choicemodel(g), fitnessfunction, samplesize, minimumsamplesize, totalsamplelimit, ruledepthlimit)
end

function setnmcschoicemodel!(g::Generator, fitnessfunction::Function, samplesize::Int=1; minimumsamplesize::Int=0, totalsamplelimit::Int=typemax(Int), ruledepthlimit::Int=typemax(Int))
	Base.depwarn("Use setchoicemodel!(g, NMCSChoiceModel(g, ...)) in place of setnmcschoicemodel!(g, ...)", :setnmcschoicemodel!)
	setchoicemodel!(g, NMCSChoiceModel(g, fitnessfunction, samplesize; minimumsamplesize=minimumsamplesize, totalsamplelimit=totalsamplelimit, ruledepthlimit=ruledepthlimit))
end

# reset any state
function resetstate!(cm::NMCSChoiceModel)
	resetstate!(cm.policychoicemodel)
	cm.totalsamplecount = 0
	cm.bestfitness = +Inf
	cm.bestgodelsequence = (Real)[]
	cm.besttracesequence =  Dict[]
end

function godelnumber(cm::NMCSChoiceModel, cc::ChoiceContext)
	existinggodelsequence = cc.derivationstate.godelsequence
	existingtracesequence = map(t->t[2], cc.derivationstate.cmtrace) # cmtrace is tuples of cpid and trace for that cp; need here just the trace
	s = 0 # local sample count (for this choice point)
	# if total sample count is exceeded, no further samples are taken, and instead godel number
	# is taken from the best godel sequence found so far
	validruledepth = (cm.ruledepthlimit == typemax(Int)) ? true : getruledepth(cc.derivationstate) <= cm.ruledepthlimit
	while (s < cm.samplesize) && ((s < cm.minimumsamplesize) || ((validruledepth) && (cm.totalsamplecount < cm.totalsamplelimit)))
		s += 1
		cm.totalsamplecount += 1
		simulationgenerator = deepcopy(cc.derivationstate.generator)
		simulationcm = NMCSSimulationChoiceModel(deepcopy(cm.policychoicemodel), deepcopy(existinggodelsequence), deepcopy(existingtracesequence))
		setchoicemodel!(simulationgenerator, simulationcm)
		result, state = nothing, nothing
		try
			result, state = generate(simulationgenerator; maxchoices=cc.derivationstate.maxchoices)
		catch e
		  if isa(e,GenerationTerminatedException)
				continue # skip the remainder of this loop iteration
			else
				throw(e)
			end
		end
		fitness = cm.fitnessfunction(result)
		if fitness <= cm.bestfitness # <= rather than < in order to improve diversity: if new sequence is just as good, prefer it to the stored best
			cm.bestfitness = fitness
			cm.bestgodelsequence = deepcopy(state.godelsequence)
			cm.besttracesequence = map(t->t[2], state.cmtrace) # cmtrace is tuples of cpid and trace for that cp; need here just the trace
		end
	end
	if isempty(cm.bestgodelsequence)
		# if all of the loop iterations above caught a GenerationTerminatedException AND this is the first choice point, then no bestsequence
		# is set - throw an exception since we cannot return a choice for the choice point
		throw(GenerationTerminatedException("for all simulations run at the first choice point in NMCS, the number of the choices made exceeded $(cc.derivationstate.maxchoices): specify a larger value of maxchoices as a parameter to generate, or increase the NMCS sample size"))
	end
	gn = cm.bestgodelsequence[length(cc.derivationstate.godelsequence)+1]
	trace = cm.besttracesequence[length(cc.derivationstate.cmtrace)+1]
	gn, trace
end

setparams!(cm::NMCSChoiceModel, params) = setparams!(cm.policychoicemodel, params)
getparams(cm::NMCSChoiceModel) = getparams(cm.policychoicemodel)
paramranges(cm::NMCSChoiceModel) = paramranges(cm.policychoicemodel)

type NMCSSimulationChoiceModel <: ChoiceModel
	policychoicemodel::ChoiceModel
	presetgodelsequence::Vector{Real}
	presettracesequence::Vector{Dict}
end

function godelnumber(cm::NMCSSimulationChoiceModel, cc::ChoiceContext)
	if isempty(cm.presetgodelsequence)
		gn, trace = godelnumber(cm.policychoicemodel, cc)
	else
		gn, trace = shift!(cm.presetgodelsequence), shift!(cm.presettracesequence)
	end
	gn, trace
end

# reset any state
function resetstate!(cm::NMCSSimulationChoiceModel)
	resetstate!(cm.policychoicemodel)
end

show(io::IO, cm::NMCSChoiceModel) = print(io, "NMCS choice model (policy: ", cm.policychoicemodel, ")")


text:
# AUTO GENERATED FILE - DO NOT EDIT

export vtk_mesh

"""
    vtk_mesh(;kwargs...)

A Mesh component.
Mesh is exposing a vtkPolyData to a downstream filter
It takes the following set of properties:
  - state: { mesh: { ...polydata-props }, field: { ...dataArray } }
Keyword arguments:
- `id` (String; optional): The ID used to identify this component.
- `port` (Real; optional): downstream connection port
- `state` (Dict; optional): State of the mesh
"""
function vtk_mesh(; kwargs...)
        available_props = Symbol[:id, :port, :state]
        wild_props = Symbol[]
        return Component("vtk_mesh", "Mesh", "dash_vtk", available_props, wild_props; kwargs...)
end

max_stars_repo_name: rodrigomologni/dash-vtk

max_issues_repo_name: rodrigomologni/dash-vtk

max_forks_repo_name: rodrigomologni/dash-vtk


text:
module SubspaceInference

using Flux

using Flux: Data.DataLoader
using Flux: @epochs

using LinearAlgebra
using LowRankApprox
using Statistics
using PyPlot
using Distributions
using AdvancedMH, AdvancedHMC, AdvancedVI
using MCMCChains
using ForwardDiff, Zygote, ReverseDiff
using DiffEqFlux
using DifferentialEquations
using DiffResults
using StructArrays
using ManifoldLearning
using LazyArrays, Memoization, Turing, DistributionsAD

###########
# Exports #
###########
export  subspace_construction,
		auto_encoder_subspace,
		diffusion_subspace,
		subspace_inference,
		sub_inference,
		autoencoder_inference,
		auto_inference,
		turing_inference

#include functions
include("plotting.jl")
include("subspace_construction.jl")
include("libs.jl")
include("space_inference.jl")
include("uncertainty.jl")

end # module


text:
"""
random(a, b [, T = Float64])

Returns a random number between `a` and `b`. The result is converted to the type `T`.
"""
random(a::Real, b::Real, ::Type{T} = Float64) where {T <: Real} = convert(T, (b - a) * rand(Float64) + a)

"""
random(len [, T = Float64])

Returns a random number between `0` and `len`. The result is converted to the type `T`.
"""
random(len::Real, ::Type{T} = Float64) where {T} = random(0, len, T)

"""
        point(v::Vector)

Converts a [`Vector`](@ref) to a [`Point`](@ref) object.
"""
point(v::Vector{<:Real}) = Luxor.Point(v.x, v.y)

"""
        mapr(s::Real, intervals::Matrix)

Return the value `s` linearly mapped from the intervals in `intervals`.

# Examples

```julia
julia> mapr(1, [0 1; 0 10])
10.0

julia> mapr(0.5, [0 1; 0 10])
5.0

julia> mapr(-2, [-5 5; 0 5])
1.5
```
"""
function mapr(s::Real, intervals::Matrix{<:Real})
    a₁, a₂ = intervals[1, :]
    b₁, b₂ = intervals[2, :]
    return b₁ + (s - a₁) * (b₂ - b₁) / (a₂ - a₁)
end

"""
        s ⟶ [a₁ a₂; b₁ b₂]

Return the value `s` linearly mapped from the `a` interval to the `b` interval.

# Examples

```julia
julia> 1 ⟶ [0 1; 0 10]
10.0

julia> 0.5 ⟶ [0 1; 0 10]
5.0

julia> -2 ⟶ [-5 5; 0 5]
1.5
```
"""
⟶(s::Real, intervals::Matrix{<:Real}) = mapr(s, intervals)


text:
abstract type EnergyBC end

"""
    Insulating() :: EnergyBC

No energy flux across the boundary.
"""
struct Insulating <: EnergyBC end
function atmos_energy_boundary_state!(nf, bc_energy::Insulating, atmos, args...) end
function atmos_energy_normal_boundary_flux_diffusive!(
    nf,
    bc_energy::Insulating,
    atmos,
    args...,
) end

"""
    PrescribedTemperature(fn) :: EnergyBC

Prescribe the temperature at the boundary by `fn`, a function with signature
`fn(state, aux, t)` returning the temperature (in K).
"""
struct PrescribedTemperature{FN} <: EnergyBC
    fn::FN
end
function atmos_energy_boundary_state!(
    nf,
    bc_energy::PrescribedTemperature,
    atmos,
    state⁺,
    aux⁺,
    n,
    state⁻,
    aux⁻,
    bctype,
    t,
    args...,
)
    FT = eltype(aux⁻)
    _T_0::FT = T_0(atmos.param_set)
    _cv_d::FT = cv_d(atmos.param_set)

    T = bc_energy.fn(state⁻, aux⁻, t)
    E_int⁺ = state⁺.ρ * _cv_d * (T - _T_0)
    state⁺.ρe =
        E_int⁺ + state⁺.ρ * gravitational_potential(atmos.orientation, aux⁻)
end
function atmos_energy_normal_boundary_flux_diffusive!(
    nf,
    bc_energy::PrescribedTemperature,
    atmos,
    fluxᵀn,
    n⁻,
    state⁻,
    diff⁻,
    hyperdiff⁻,
    aux⁻,
    state⁺,
    diff⁺,
    hyperdiff⁺,
    aux⁺,
    bctype,
    t,
    args...,
)

    # TODO: figure out a better way...
    ν, D_t, _ = turbulence_tensors(atmos, state⁻, diff⁻, aux⁻, t)
    d_h_tot = -D_t .* diff⁻.∇h_tot
    nd_h_tot = dot(n⁻, d_h_tot)
    # both sides involve projections of normals, so signs are consistent
    fluxᵀn.ρe += nd_h_tot * state⁻.ρ
end

"""
    PrescribedEnergyFlux(fn) :: EnergyBC

Prescribe the net inward energy flux across the boundary by `fn`, a function
with signature `fn(state, aux, t)`, returning the flux (in W/m^2).
"""
struct PrescribedEnergyFlux{FN} <: EnergyBC
    fn::FN
end
function atmos_energy_boundary_state!(
    nf,
    bc_energy::PrescribedEnergyFlux,
    atmos,
    args...,
) end
function atmos_energy_normal_boundary_flux_diffusive!(
    nf,
    bc_energy::PrescribedEnergyFlux,
    atmos,
    fluxᵀn,
    n⁻,
    state⁻,
    diff⁻,
    hyperdiff⁻,
    aux⁻,
    state⁺,
    diff⁺,
    hyperdiff⁺,
    aux⁺,
    bctype,
    t,
    args...,
)

    # DG normal is defined in the outward direction
    # we want to prescribe the inward flux
    fluxᵀn.ρe -= bc_energy.fn(state⁻, aux⁻, t)
end

max_stars_repo_name: yairchn/CLIMA

max_issues_repo_name: yairchn/CLIMA

max_forks_repo_name: yairchn/CLIMA


text:
### A Pluto.jl notebook ###
# v0.14.2

using Markdown
using InteractiveUtils

# ╔═╡ 3fa5970e-8d82-11eb-302e-d53a453e984f
begin
	import Pkg;
	Pkg.add("SpecialFunctions");
	using SpecialFunctions
end

# ╔═╡ 2e8a0d14-8d6a-11eb-3a36-01676cf20447
begin
	Pkg.add(["VegaLite", "VegaDatasets"]);
	using VegaLite, VegaDatasets
end

# ╔═╡ 8d4852d8-88d4-11eb-10e6-51c750d36b54
begin
	Pkg.add("PlutoUI")
	using PlutoUI
end

# ╔═╡ f1261dab-6736-4b0a-803b-b2d2836ca24a
Pkg.add("ColorSchemes")

# ╔═╡ 1ccb3a84-88d4-11eb-2499-91af66e78e89
using Plots

# ╔═╡ a400134c-8dc7-11eb-3a29-f99c5910de8c
using DataFrames

# ╔═╡ 47452d72-88d6-11eb-27ef-bbc1d061060d
md"""
## Resources

* Plots Docs [Useful Tips](https://docs.juliaplots.org/latest/basics/#Useful-Tips)
* Plots Docs [Examples](https://docs.juliaplots.org/latest/generated/gr/)
* [Plot Attributes](http://docs.juliaplots.org/latest/generated/attributes_plot/)
* [Axis Attributes](http://docs.juliaplots.org/latest/generated/attributes_axis/#Axis)
* [Color Names](http://juliagraphics.github.io/Colors.jl/stable/namedcolors/)
"""

# ╔═╡ 2a2338e4-88d4-11eb-0e6b-03609b700baa
md"""
## Loading Time

Plots loading time sucks big time.
This is a known problem.  We know, we know, Julia is fast, just grab some coffee ....
☕☕☕☕☕☕☕☕☕
"""

# ╔═╡ 4316799c-88d4-11eb-0e83-a7df319711ad
md"""
## Hello World
"""

# ╔═╡ 6fcad280-88d4-11eb-310f-5f89a5d21f9b
plot(rand(4))

# ╔═╡ 63da4054-88d5-11eb-29c9-b3648548c367
md"""
## Adding Stuff
"""

# ╔═╡ 6b55e978-88d5-11eb-1fea-739e61bbb35a
begin
		plot!(rand(4))
		scatter!(rand(4))
end

# ╔═╡ 198eb00c-88d5-11eb-3d3c-3963d197f0e0
md"""
## Removing Stuff
"""

# ╔═╡ 2b4cc1d0-88d5-11eb-0afd-3988abd9a870
plot!(legend=false, axis=false, grid=false, ticks=false )

# ╔═╡ dd81f7a5-1acf-430d-9d25-ab2877bf34cf
md"""
## Lines and points
"""

# ╔═╡ cf89db1e-88d7-11eb-2a92-850c7d46a296
md"""
## Square Aspect Ratio
"""

# ╔═╡ da7b4c10-88d7-11eb-011e-dbb639e6fa2b
begin
	v = rand(4)
	plot(v, ratio=1, legend=false)
	scatter!(v)
end

# ╔═╡ cca58686-88d8-11eb-2fe4-d7d378c5408a
md"""
##  Matrices with color (heatmap)
yflip = true places the 2,1 entry where you want it, since you want to interchange xy with ij and i runs down and y runs up.
"""

# ╔═╡ d2c81aec-88d8-11eb-0e25-57a2068c3bf9
A = [ 1 1000 1 ; 1 1 1; 1 1 1]

# ╔═╡ e3f7d1ea-88d8-11eb-21bf-cbcdec1e8810
heatmap(A, ratio=1, yflip=true, legend=false, axis=false, grid=false,  ticks=false)

# ╔═╡ b059780e-88db-11eb-028c-6b07355fb1ab
heatmap( rand(10,10), clim=(0,1), ratio=1, legend=false, axis=false, ticks=false)

# ╔═╡ 48561d5a-8991-11eb-2a55-db793e4c9fea
begin
	MM = [ 0 1 0; 0 0 0; 1 0 0]

	whiteblack = [RGBA(1,1,1,0), RGB(0,0,0)]
	 heatmap(c=whiteblack, MM, aspect_ratio = 1, ticks=.5:3.5, lims=(.5,3.5), gridalpha=1, legend=false, axis=false, ylabel="i", xlabel="j")
end

# ╔═╡ 1bc15efe-8a44-11eb-11e4-5349f4279202
begin
	p=plot(1:4,guideposition=:top)
end

# ╔═╡ 63a807fe-8a44-11eb-3a3f-23fdf0804149
begin
	p.attr[:foreground_color] = RGB(1,0,0)
	p.attr[:foreground_color]
end

# ╔═╡ a0cd1d04-8a44-11eb-248c-1f5d699d60d2
p

# ╔═╡ b3aba314-8a44-11eb-0ea1-c9aafbecbd7e
begin
		q = (p.series_list[1]).plotattributes
	    q.explicit[:linecolor]=:green
end

# ╔═╡ 5bb7a350-8a45-11eb-15eb-d76bc70ab2bc
p

# ╔═╡ 15e03bd0-898c-11eb-3424-6fc74206ce9a
md"""
## Colors
"""

# ╔═╡ 855cae06-898e-11eb-36d5-0ddaf08022e9
Colors.color_names

# ╔═╡ 1d8dd040-898c-11eb-387b-3143b3997eee
mycolors = [colorant"lightslateblue",colorant"limegreen",colorant"red"]

# ╔═╡ 3d2c46d0-898b-11eb-3073-e1012b3ebb67
begin

	AA = [i for i=50:300, j=1:100]

	heatmap(AA, c=mycolors, clim=(1,300))

end

# ╔═╡ 0b8a01aa-8993-11eb-2e43-074f49edc175
md"""
## Area under curves
"""

# ╔═╡ e2046a88-8991-11eb-0906-2f0e89a943c3
begin
	y = rand(10)
	plot(y, fillrange= y.*0 .+ .5, label= "above/below 1/2", legend =:topleft)
end

# ╔═╡ a990150c-8992-11eb-22f5-bfdf98f17298
begin
	x = LinRange(0,2,100)
	y1 = exp.(x)
	y2 = exp.(1.3 .* x)

	plot(x, y1, fillrange = y2, fillalpha = 0.35, c = 1, label = "Confidence band", legend = :topleft)
end

# ╔═╡ 940bd6dc-8c12-11eb-09a2-89206356c143
let
	x = -3:.01:3
	areaplot(x, exp.(-x.^2/2)/√(2π),alpha=.25,legend=false)
end

# ╔═╡ 78727c02-8c13-11eb-16b7-bd7d4414577d
begin
	M = [1 2 3; 7 8 9; 4 5 6; 0 .5 1.5]
	areaplot(1:3, M, seriescolor = [:red :green :blue ], fillalpha = [0.2 0.3 0.4])
end

# ╔═╡ b8435896-8d81-11eb-3d44-7db7b987f992
let
	f = x->exp(-x^2/2)/√(2π)
	δ = .01
	plot()
	x = √2 .* erfinv.(2 .*(δ/2 : δ : 1) .- 1)
	areaplot( x, f.(x), seriescolor=[ :red,:blue], legend=false)
	plot!( f, x,c=:black)
end

# ╔═╡ d2c45794-8c12-11eb-17bb-a37bae10d084
md"""
## Shapes
"""

# ╔═╡ dcafeb72-8bce-11eb-23ae-37cad92e0f82
begin
	rectangle(w, h, x, y) = Shape(x .+ [0,w,w,0], y .+ [0,0,h,h])
	circle(r,x,y) = (θ = LinRange(0,2π,500); (x.+r.*cos.(θ), y.+r.*sin.(θ)))
end

# ╔═╡ d7c8534e-8c12-11eb-211f-6b8718e52dd7
begin
	plot(circle(5,0,0), ratio=1, c=:red, fill=true)
	plot!(rectangle(5*√2,5*√2,-2.5*√2,-2.5*√2),c=:white,fill=true,legend=false)
end

# ╔═╡ c798a01a-8b51-11eb-2517-c7f45234db58
md"""
## Editing Plots manually
"""

# ╔═╡ dbfaa24e-8b51-11eb-3c20-bf7919c3a167
pl = plot(1:4,[1, 4, 9, 16])

# ╔═╡ ebbb02a0-8b51-11eb-08ef-4bb2f55782c1
pl.attr

# ╔═╡ f9490ca8-8b49-11eb-3bab-a3569a181f63
pl.series_list[1]

# ╔═╡ 131aa286-8b4a-11eb-3b4b-fb34a39cc7fb
pl[:size]=(300,200)

# ╔═╡ 111bec42-8b4a-11eb-08ca-e194014175f7
pl

# ╔═╡ f5df0246-8d80-11eb-0fb8-e79a0d5d7a23
md"""
##  VegaLite
"""

# ╔═╡ 78f0c7c6-8d6a-11eb-29f4-af99ddf71960
@vlplot(:point, rand(10), rand(10) )

# ╔═╡ 129de546-8d81-11eb-2ea0-a3e60e88ae13
begin
	dataset("zipcodes") |>
	@vlplot(
	    :circle,
	    width=500, height=300,
	    transform=[{calculate="substring(datum.zip_code, 0, 1)", as=:digit}],
	    projection={type=:albersUsa},
	    longitude=:longitude,
	    latitude=:latitude,
	    size={value=1},
	    color="digit:n"
	)
end

# ╔═╡ 837a7086-8dc7-11eb-3397-e10d6f40481f
 d = DataFrame(dataset("zipcodes"))

# ╔═╡ c0494eba-8dc7-11eb-3c1a-1f19979a4f4a
scatter( d[!,:longitude], d[!,:latitude], m=:., ms=1, xlim=(-150,-50), ylim=(20,55))

# ╔═╡ 3395ee1c-8d81-11eb-3e94-b1373d69dc93

# ╔═╡ 6dbc29fc-8da3-11eb-3196-e1e235eb6e86
begin

	struct RankOneMatrix{T}
	  v::AbstractVector{T}
	   w::AbstractVector{T}
	end

end

# ╔═╡ 2e9d597e-8da5-11eb-114b-5f87731d4ced

# ╔═╡ 41e5e100-8da4-11eb-3bb8-d9d414af8e2a
RankOneMatrix( rand(3) , rand(3) )

# ╔═╡ e9b2b30e-8da4-11eb-04b3-8ba421898f0a
methods(RankOneMatrix)

# ╔═╡ 311ce496-8dad-11eb-330c-d3757a6fac8f
begin
	xx = .1:.1:1
	plot(xx.^2, xaxis=:log, yaxis=:log)
end

# ╔═╡ 5ae674d6-8dad-11eb-3abf-6fdbbb6cfe43
begin

	plot(exp.(x), yaxis=:log)
end

# ╔═╡ eb3e721e-88d4-11eb-1f09-cfba69f498d4
TableOfContents()

# ╔═╡ Cell order:
# ╠═1ccb3a84-88d4-11eb-2499-91af66e78e89
# ╟─47452d72-88d6-11eb-27ef-bbc1d061060d
# ╟─2a2338e4-88d4-11eb-0e6b-03609b700baa
# ╟─4316799c-88d4-11eb-0e83-a7df319711ad
# ╠═6fcad280-88d4-11eb-310f-5f89a5d21f9b
# ╟─63da4054-88d5-11eb-29c9-b3648548c367
# ╠═6b55e978-88d5-11eb-1fea-739e61bbb35a
# ╟─198eb00c-88d5-11eb-3d3c-3963d197f0e0
# ╠═2b4cc1d0-88d5-11eb-0afd-3988abd9a870
# ╟─dd81f7a5-1acf-430d-9d25-ab2877bf34cf
# ╟─cf89db1e-88d7-11eb-2a92-850c7d46a296
# ╠═da7b4c10-88d7-11eb-011e-dbb639e6fa2b
# ╟─cca58686-88d8-11eb-2fe4-d7d378c5408a
# ╠═d2c81aec-88d8-11eb-0e25-57a2068c3bf9
# ╠═e3f7d1ea-88d8-11eb-21bf-cbcdec1e8810
# ╠═b059780e-88db-11eb-028c-6b07355fb1ab
# ╠═48561d5a-8991-11eb-2a55-db793e4c9fea
# ╠═1bc15efe-8a44-11eb-11e4-5349f4279202
# ╠═63a807fe-8a44-11eb-3a3f-23fdf0804149
# ╠═a0cd1d04-8a44-11eb-248c-1f5d699d60d2
# ╠═b3aba314-8a44-11eb-0ea1-c9aafbecbd7e
# ╠═5bb7a350-8a45-11eb-15eb-d76bc70ab2bc
# ╟─15e03bd0-898c-11eb-3424-6fc74206ce9a
# ╠═855cae06-898e-11eb-36d5-0ddaf08022e9
# ╠═1d8dd040-898c-11eb-387b-3143b3997eee
# ╠═3d2c46d0-898b-11eb-3073-e1012b3ebb67
# ╟─0b8a01aa-8993-11eb-2e43-074f49edc175
# ╠═e2046a88-8991-11eb-0906-2f0e89a943c3
# ╠═a990150c-8992-11eb-22f5-bfdf98f17298
# ╠═940bd6dc-8c12-11eb-09a2-89206356c143
# ╠═78727c02-8c13-11eb-16b7-bd7d4414577d
# ╠═3fa5970e-8d82-11eb-302e-d53a453e984f
# ╠═b8435896-8d81-11eb-3d44-7db7b987f992
# ╟─d2c45794-8c12-11eb-17bb-a37bae10d084
# ╠═dcafeb72-8bce-11eb-23ae-37cad92e0f82
# ╠═d7c8534e-8c12-11eb-211f-6b8718e52dd7
# ╠═c798a01a-8b51-11eb-2517-c7f45234db58
# ╠═dbfaa24e-8b51-11eb-3c20-bf7919c3a167
# ╠═ebbb02a0-8b51-11eb-08ef-4bb2f55782c1
# ╠═f9490ca8-8b49-11eb-3bab-a3569a181f63
# ╠═131aa286-8b4a-11eb-3b4b-fb34a39cc7fb
# ╠═111bec42-8b4a-11eb-08ca-e194014175f7
# ╟─f5df0246-8d80-11eb-0fb8-e79a0d5d7a23
# ╠═2e8a0d14-8d6a-11eb-3a36-01676cf20447
# ╠═78f0c7c6-8d6a-11eb-29f4-af99ddf71960
# ╠═129de546-8d81-11eb-2ea0-a3e60e88ae13
# ╠═837a7086-8dc7-11eb-3397-e10d6f40481f
# ╠═c0494eba-8dc7-11eb-3c1a-1f19979a4f4a
# ╠═a400134c-8dc7-11eb-3a29-f99c5910de8c
# ╟─3395ee1c-8d81-11eb-3e94-b1373d69dc93
# ╠═6dbc29fc-8da3-11eb-3196-e1e235eb6e86
# ╠═2e9d597e-8da5-11eb-114b-5f87731d4ced
# ╠═41e5e100-8da4-11eb-3bb8-d9d414af8e2a
# ╠═e9b2b30e-8da4-11eb-04b3-8ba421898f0a
# ╠═311ce496-8dad-11eb-330c-d3757a6fac8f
# ╠═5ae674d6-8dad-11eb-3abf-6fdbbb6cfe43
# ╠═8d4852d8-88d4-11eb-10e6-51c750d36b54
# ╟─eb3e721e-88d4-11eb-1f09-cfba69f498d4
# ╠═f1261dab-6736-4b0a-803b-b2d2836ca24a


text:
module ReinforcementLearning

include("Environment.jl"); export AbsEnvironment, AbsAction, transfer, getInitialState, getAllStates, AbsState, getActions, isTerminal, getSuccessors, render_env
include("Policy.jl"); export Policy, actionsAndProbs, prettyprint
include("Agent.jl"); export AbsAgent, play, observe, RandomAgent, QLearner, SarsaLearner, SarsaLambdaLearner, PolicyAgent, MonteCarloAgent
include("Engine.jl"); export playEpisode
include("MDP.jl"); export MDP, State, Action, printmdp
include("DP.jl"); export Policy, actionsAndProbs, iterative_policy_evaluation, policy_iteration, synchronous_value_iteration, gauss_seidel_value_iteration
try
    include("GymEnvironment.jl"); export GymState, GymAction, GymEnv, sample, monitor_start, monitor_close
catch
    info("You need to install openai gym to use gym specific environments.")
    info("Follow the instructions: https://github.com/openai/gym")
end
end # module

max_stars_repo_name: ozanarkancan/ReinforcementLearning

max_issues_repo_name: ozanarkancan/ReinforcementLearning

max_forks_repo_name: ozanarkancan/ReinforcementLearning


text:
using BioServices.EUtils
using BioMedQuery.PubMed
using LightXML

"""
    export_citation(pmid::Int64, citation_type, output_file,verbose)

Export, to an output file, the citation for PubMed article identified by the given pmid

## Arguments
* `citation_type::String`: At the moment supported types include: "endnote", "bibtex"
"""
function export_citation(pmid::Int64, citation_type, output_file, overwrite = true, verbose=false)

    efetch_response = efetch(db = "pubmed", tool = "BioJulia", retmode = "xml", rettype = "null", id = [pmid])

    if verbose
        xdoc = parse_string(String(efetch_response.body))
        save_file(xdoc, "./efetch.xml")
    end

    #convert xml to ezxml doc
    efetch_doc = root(parse_string(String(efetch_response.body)))
    citation_output = PubMed.CitationOutput(citation_type, output_file, overwrite)
    save_efetch!(citation_output, efetch_doc, verbose)
end

"""
    export_citation(pmids::Vector{Int64}, citation_type, output_file,verbose)

Export, to an output file, the citation for collection of PubMed articles identified by the given pmids

## Arguments
* `citation_type::String`: At the moment supported types include: "endnote", "bibtex"
"""
function export_citation(pmids::Vector{Int64}, citation_type, output_file, overwrite=true, verbose=false)
    efetch_response = efetch(db = "pubmed", tool = "BioJulia", retmode = "xml", rettype = "null", id = pmids)

    if verbose
        xdoc = parse_string(efetch_response.body)
        save_file(xdoc, "./efetch.xml")
    end

    #convert xml to ezxml doc
    efetch_doc = root(parse_string(String(efetch_response.body)))

    citation_output = PubMed.CitationOutput(citation_type, output_file, overwrite)
    save_efetch!(citation_output, efetch_doc, verbose)
end


text:
arr = [1,2,3,4,5]

map(x -> x * 10, arr) |> println

max_stars_repo_name: mkanenobu/trashbox

max_issues_repo_name: mkanenobu/trashbox

max_forks_repo_name: mkanenobu/trashbox


text:
using Test, OrdinaryDiffEq, HyperbolicDiffEq, DiffEqCallbacks
using LinearAlgebra

function save_func(u, t, integrator)
    meshx = integrator.f.f.meshx
    basis = integrator.f.f.basis
    balance_law = integrator.f.f.balance_law
    integrate(u->IntegralQuantitiesBurgers(u,balance_law), u, meshx, basis)
end

function compute_solution(balance_law, p, N, basis_type, fvol_type, fnum_type, cfl,
                            mesh_type=UniformPeriodicMesh1D, parallel=Val{:serial}())
    xmin, xmax = -1., 1.
    meshx = mesh_type(xmin, xmax, N)
    if basis_type == :LobattoLegendre
        basis = LobattoLegendre(p)
    elseif basis_type == :GaussLegendre
        basis = GaussLegendre(p)
    else
        error("basis_type $basis_type unknown.")
    end

    if fvol_type == :Central
        fvol = CentralFlux()
    elseif fvol_type == :EC
        fvol = EnergyConservativeFlux()
    else
        error("fvol_type $fvol_type unknown.")
    end
    if fnum_type == :Central
        fnumint = CentralFlux()
    elseif fnum_type == :EC
        fnumint = EnergyConservativeFlux()
    elseif fnum_type == :Godunov
        fnumint = GodunovFlux()
    elseif fnum_type == :LLF
        fnumint = LocalLaxFriedrichsFlux()
    elseif fnum_type == :HLL
        fnumint = HartenLaxVanLeerFlux()
    else
        error("fnum_type $fnum_type unknown.")
    end
    if mesh_type <: UniformPeriodicMesh1D
        semidisc = UniformPeriodicFluxDiffDisc1D(balance_law, meshx, basis, fvol,
                                                    fnumint, parallel)
    else
        semidisc = UniformFluxDiffDisc1D(balance_law, meshx, basis, fvol, fnumint,
                                            GodunovFlux(), zero, zero, parallel)
    end

    if typeof(balance_law) <: ConstantLinearAdvection
        tspan = (0., 1.0)
    elseif typeof(balance_law) <: Burgers
        tspan = (0., 0.3)
    elseif typeof(balance_law) <: Cubic
        tspan = (0., 0.1)
    elseif typeof(balance_law) <: Quartic
        tspan = (0., 0.065)
    elseif typeof(balance_law) <: Quintic
        tspan = (0., 0.045)
    elseif typeof(balance_law) <: Sextic
        tspan = (0., 0.035)
    elseif typeof(balance_law) <: Septic
        tspan = (0., 0.025)
    elseif typeof(balance_law) <: Octic
        tspan = (0., 0.020)
    else
        error("Balance law $balance_law not supported.")
    end

    u₀ = sinpi

    ode = semidiscretise(semidisc, u₀, tspan)
    tstops = range(tspan[1], tspan[end], length=2)
    maxdt = (u,p,t) -> max_dt(t, u, semidisc, cfl)
    saved_values = SavedValues(Vector{Float64}(), Vector{IntegralQuantitiesBurgers{Float64}}())
    cb = CallbackSet(StepsizeLimiter(maxdt), SavingCallback(save_func, saved_values))
    sol = solve(ode, SSPRK104(), dt=maxdt(ode.u0,nothing,tspan[1]), saveat=tstops, tstops=tstops,
                dense=false, callback=cb)
    sol, saved_values
end

function mass_and_energy_differences(balance_law, basis_type, cfl)
    p = 4
    N = 20
    fvol_type = :EC
    fnum_type = :EC

    sol, saved_values = compute_solution(balance_law, p, N, basis_type, fvol_type, fnum_type, cfl)
    mass_diff = saved_values.saveval[1].mass-saved_values.saveval[end].mass
    ener_diff = saved_values.saveval[1].energy-saved_values.saveval[end].energy

    mass_diff, ener_diff
end

mass_diff, ener_diff = mass_and_energy_differences(ConstantLinearAdvection(), :LobattoLegendre, 0.9)
@test abs(mass_diff) < 10*eps()
@test ener_diff < 5.e-8
mass_diff, ener_diff = mass_and_energy_differences(ConstantLinearAdvection(), :LobattoLegendre, 0.09)
@test abs(mass_diff) < 10*eps()
@test ener_diff < 5.e-12
mass_diff, ener_diff = mass_and_energy_differences(ConstantLinearAdvection(), :GaussLegendre, 0.9)
@test abs(mass_diff) < 10*eps()
@test ener_diff < 5.e-8
mass_diff, ener_diff = mass_and_energy_differences(ConstantLinearAdvection(), :GaussLegendre, 0.09)
@test abs(mass_diff) < 10*eps()
@test ener_diff < 5.e-12
sol_serial, _  = compute_solution(ConstantLinearAdvection(), 4, 20, :LobattoLegendre, :EC, :EC, 0.9, UniformPeriodicMesh1D, Val{:serial}())
sol_threads, _ = compute_solution(ConstantLinearAdvection(), 4, 20, :LobattoLegendre, :EC, :EC, 0.9, UniformPeriodicMesh1D, Val{:threads}())
@test norm(sol_serial[end] - sol_threads[end], Inf) < 10*eps()
sol_serial, _  = compute_solution(ConstantLinearAdvection(), 4, 20, :LobattoLegendre, :EC, :EC, 0.9, UniformMesh1D, Val{:serial}())
sol_threads, _ = compute_solution(ConstantLinearAdvection(), 4, 20, :LobattoLegendre, :EC, :EC, 0.9, UniformMesh1D, Val{:threads}())
@test norm(sol_serial[end] - sol_threads[end], Inf) < 10*eps()
sol_serial, _  = compute_solution(ConstantLinearAdvection(), 4, 20, :GaussLegendre, :EC, :EC, 0.9, UniformPeriodicMesh1D, Val{:serial}())
sol_threads, _ = compute_solution(ConstantLinearAdvection(), 4, 20, :GaussLegendre, :EC, :EC, 0.9, UniformPeriodicMesh1D, Val{:threads}())
@test norm(sol_serial[end] - sol_threads[end], Inf) < 10*eps()
sol_serial, _  = compute_solution(ConstantLinearAdvection(), 4, 20, :GaussLegendre, :EC, :EC, 0.9, UniformMesh1D, Val{:serial}())
sol_threads, _ = compute_solution(ConstantLinearAdvection(), 4, 20, :GaussLegendre, :EC, :EC, 0.9, UniformMesh1D, Val{:threads}())
@test norm(sol_serial[end] - sol_threads[end], Inf) < 10*eps()

mass_diff, ener_diff = mass_and_energy_differences(Burgers(), :LobattoLegendre, 0.9)
@test abs(mass_diff) < 10*eps()
@test ener_diff < 5.e-8
mass_diff, ener_diff = mass_and_energy_differences(Burgers(), :LobattoLegendre, 0.09)
@test abs(mass_diff) < 10*eps()
@test ener_diff < 5.e-12
mass_diff, ener_diff = mass_and_energy_differences(Burgers(), :GaussLegendre, 0.9)
@test abs(mass_diff) < 10*eps()
@test ener_diff < 5.e-8
mass_diff, ener_diff = mass_and_energy_differences(Burgers(), :GaussLegendre, 0.09)
@test abs(mass_diff) < 10*eps()
@test ener_diff < 5.e-12
sol_serial, _  = compute_solution(Burgers(), 4, 20, :LobattoLegendre, :EC, :EC, 0.9, UniformPeriodicMesh1D, Val{:serial}())
sol_threads, _ = compute_solution(Burgers(), 4, 20, :LobattoLegendre, :EC, :EC, 0.9, UniformPeriodicMesh1D, Val{:threads}())
@test norm(sol_serial[end] - sol_threads[end], Inf) < 10*eps()
sol_serial, _  = compute_solution(Burgers(), 4, 20, :LobattoLegendre, :EC, :EC, 0.9, UniformMesh1D, Val{:serial}())
sol_threads, _ = compute_solution(Burgers(), 4, 20, :LobattoLegendre, :EC, :EC, 0.9, UniformMesh1D, Val{:threads}())
@test norm(sol_serial[end] - sol_threads[end], Inf) < 10*eps()
sol_serial, _  = compute_solution(Burgers(), 4, 20, :GaussLegendre, :EC, :EC, 0.9, UniformPeriodicMesh1D, Val{:serial}())
sol_threads, _ = compute_solution(Burgers(), 4, 20, :GaussLegendre, :EC, :EC, 0.9, UniformPeriodicMesh1D, Val{:threads}())
@test norm(sol_serial[end] - sol_threads[end], Inf) < 10*eps()
sol_serial, _  = compute_solution(Burgers(), 4, 20, :GaussLegendre, :EC, :EC, 0.9, UniformMesh1D, Val{:serial}())
sol_threads, _ = compute_solution(Burgers(), 4, 20, :GaussLegendre, :EC, :EC, 0.9, UniformMesh1D, Val{:threads}())
@test norm(sol_serial[end] - sol_threads[end], Inf) < 10*eps()

mass_diff, ener_diff = mass_and_energy_differences(Cubic(), :LobattoLegendre, 0.9)
@test abs(mass_diff) < 10*eps()
@test ener_diff < 5.e-8
mass_diff, ener_diff = mass_and_energy_differences(Cubic(), :LobattoLegendre, 0.09)
@test abs(mass_diff) < 10*eps()
@test ener_diff < 5.e-12
mass_diff, ener_diff = mass_and_energy_differences(Cubic(), :GaussLegendre, 0.9)
@test abs(mass_diff) < 10*eps()
@test ener_diff < 5.e-8
mass_diff, ener_diff = mass_and_energy_differences(Cubic(), :GaussLegendre, 0.09)
@test abs(mass_diff) < 10*eps()
@test ener_diff < 5.e-12
sol_serial, _  = compute_solution(Cubic(), 4, 20, :LobattoLegendre, :EC, :EC, 0.9, UniformPeriodicMesh1D, Val{:serial}())
sol_threads, _ = compute_solution(Cubic(), 4, 20, :LobattoLegendre, :EC, :EC, 0.9, UniformPeriodicMesh1D, Val{:threads}())
@test norm(sol_serial[end] - sol_threads[end], Inf) < 10*eps()
sol_serial, _  = compute_solution(Cubic(), 4, 20, :LobattoLegendre, :EC, :EC, 0.9, UniformMesh1D, Val{:serial}())
sol_threads, _ = compute_solution(Cubic(), 4, 20, :LobattoLegendre, :EC, :EC, 0.9, UniformMesh1D, Val{:threads}())
@test norm(sol_serial[end] - sol_threads[end], Inf) < 10*eps()
sol_serial, _  = compute_solution(Cubic(), 4, 20, :GaussLegendre, :EC, :EC, 0.9, UniformPeriodicMesh1D, Val{:serial}())
sol_threads, _ = compute_solution(Cubic(), 4, 20, :GaussLegendre, :EC, :EC, 0.9, UniformPeriodicMesh1D, Val{:threads}())
@test norm(sol_serial[end] - sol_threads[end], Inf) < 10*eps()
sol_serial, _  = compute_solution(Cubic(), 4, 20, :GaussLegendre, :EC, :EC, 0.9, UniformMesh1D, Val{:serial}())
sol_threads, _ = compute_solution(Cubic(), 4, 20, :GaussLegendre, :EC, :EC, 0.9, UniformMesh1D, Val{:threads}())
@test norm(sol_serial[end] - sol_threads[end], Inf) < 10*eps()

mass_diff, ener_diff = mass_and_energy_differences(Quartic(), :LobattoLegendre, 0.9)
@test abs(mass_diff) < 10*eps()
@test ener_diff < 5.e-8
mass_diff, ener_diff = mass_and_energy_differences(Quartic(), :LobattoLegendre, 0.09)
@test abs(mass_diff) < 10*eps()
@test ener_diff < 5.e-12
mass_diff, ener_diff = mass_and_energy_differences(Quartic(), :GaussLegendre, 0.9)
@test abs(mass_diff) < 10*eps()
@test ener_diff < 5.e-8
mass_diff, ener_diff = mass_and_energy_differences(Quartic(), :GaussLegendre, 0.09)
@test abs(mass_diff) < 10*eps()
@test ener_diff < 5.e-12
sol_serial, _  = compute_solution(Quartic(), 4, 20, :LobattoLegendre, :EC, :EC, 0.9, UniformPeriodicMesh1D, Val{:serial}())
sol_threads, _ = compute_solution(Quartic(), 4, 20, :LobattoLegendre, :EC, :EC, 0.9, UniformPeriodicMesh1D, Val{:threads}())
@test norm(sol_serial[end] - sol_threads[end], Inf) < 10*eps()
sol_serial, _  = compute_solution(Quartic(), 4, 20, :LobattoLegendre, :EC, :EC, 0.9, UniformMesh1D, Val{:serial}())
sol_threads, _ = compute_solution(Quartic(), 4, 20, :LobattoLegendre, :EC, :EC, 0.9, UniformMesh1D, Val{:threads}())
@test norm(sol_serial[end] - sol_threads[end], Inf) < 10*eps()
sol_serial, _  = compute_solution(Quartic(), 4, 20, :GaussLegendre, :EC, :EC, 0.9, UniformPeriodicMesh1D, Val{:serial}())
sol_threads, _ = compute_solution(Quartic(), 4, 20, :GaussLegendre, :EC, :EC, 0.9, UniformPeriodicMesh1D, Val{:threads}())
@test norm(sol_serial[end] - sol_threads[end], Inf) < 10*eps()
sol_serial, _  = compute_solution(Quartic(), 4, 20, :GaussLegendre, :EC, :EC, 0.9, UniformMesh1D, Val{:serial}())
sol_threads, _ = compute_solution(Quartic(), 4, 20, :GaussLegendre, :EC, :EC, 0.9, UniformMesh1D, Val{:threads}())
@test norm(sol_serial[end] - sol_threads[end], Inf) < 10*eps()

mass_diff, ener_diff = mass_and_energy_differences(Quintic(), :LobattoLegendre, 0.9)
@test abs(mass_diff) < 10*eps()
@test ener_diff < 5.e-8
mass_diff, ener_diff = mass_and_energy_differences(Quintic(), :LobattoLegendre, 0.09)
@test abs(mass_diff) < 10*eps()
@test ener_diff < 5.e-12
mass_diff, ener_diff = mass_and_energy_differences(Quintic(), :GaussLegendre, 0.9)
@test abs(mass_diff) < 10*eps()
@test ener_diff < 5.e-8
mass_diff, ener_diff = mass_and_energy_differences(Quintic(), :GaussLegendre, 0.09)
@test abs(mass_diff) < 10*eps()
@test ener_diff < 5.e-12
sol_serial, _  = compute_solution(Quintic(), 4, 20, :LobattoLegendre, :EC, :EC, 0.9, UniformPeriodicMesh1D, Val{:serial}())
sol_threads, _ = compute_solution(Quintic(), 4, 20, :LobattoLegendre, :EC, :EC, 0.9, UniformPeriodicMesh1D, Val{:threads}())
@test norm(sol_serial[end] - sol_threads[end], Inf) < 10*eps()
sol_serial, _  = compute_solution(Quintic(), 4, 20, :LobattoLegendre, :EC, :EC, 0.9, UniformMesh1D, Val{:serial}())
sol_threads, _ = compute_solution(Quintic(), 4, 20, :LobattoLegendre, :EC, :EC, 0.9, UniformMesh1D, Val{:threads}())
@test norm(sol_serial[end] - sol_threads[end], Inf) < 10*eps()
sol_serial, _  = compute_solution(Quintic(), 4, 20, :GaussLegendre, :EC, :EC, 0.9, UniformPeriodicMesh1D, Val{:serial}())
sol_threads, _ = compute_solution(Quintic(), 4, 20, :GaussLegendre, :EC, :EC, 0.9, UniformPeriodicMesh1D, Val{:threads}())
@test norm(sol_serial[end] - sol_threads[end], Inf) < 10*eps()
sol_serial, _  = compute_solution(Quintic(), 4, 20, :GaussLegendre, :EC, :EC, 0.9, UniformMesh1D, Val{:serial}())
sol_threads, _ = compute_solution(Quintic(), 4, 20, :GaussLegendre, :EC, :EC, 0.9, UniformMesh1D, Val{:threads}())
@test norm(sol_serial[end] - sol_threads[end], Inf) < 10*eps()

mass_diff, ener_diff = mass_and_energy_differences(Sextic(), :LobattoLegendre, 0.9)
@test abs(mass_diff) < 10*eps()
@test ener_diff < 5.e-8
mass_diff, ener_diff = mass_and_energy_differences(Sextic(), :LobattoLegendre, 0.09)
@test abs(mass_diff) < 10*eps()
@test ener_diff < 5.e-12
mass_diff, ener_diff = mass_and_energy_differences(Sextic(), :GaussLegendre, 0.9)
@test abs(mass_diff) < 10*eps()
@test ener_diff < 5.e-8
mass_diff, ener_diff = mass_and_energy_differences(Sextic(), :GaussLegendre, 0.09)
@test abs(mass_diff) < 10*eps()
@test ener_diff < 5.e-12
sol_serial, _  = compute_solution(Sextic(), 4, 20, :LobattoLegendre, :EC, :EC, 0.9, UniformPeriodicMesh1D, Val{:serial}())
sol_threads, _ = compute_solution(Sextic(), 4, 20, :LobattoLegendre, :EC, :EC, 0.9, UniformPeriodicMesh1D, Val{:threads}())
@test norm(sol_serial[end] - sol_threads[end], Inf) < 10*eps()
sol_serial, _  = compute_solution(Sextic(), 4, 20, :LobattoLegendre, :EC, :EC, 0.9, UniformMesh1D, Val{:serial}())
sol_threads, _ = compute_solution(Sextic(), 4, 20, :LobattoLegendre, :EC, :EC, 0.9, UniformMesh1D, Val{:threads}())
@test norm(sol_serial[end] - sol_threads[end], Inf) < 10*eps()
sol_serial, _  = compute_solution(Sextic(), 4, 20, :GaussLegendre, :EC, :EC, 0.9, UniformPeriodicMesh1D, Val{:serial}())
sol_threads, _ = compute_solution(Sextic(), 4, 20, :GaussLegendre, :EC, :EC, 0.9, UniformPeriodicMesh1D, Val{:threads}())
@test norm(sol_serial[end] - sol_threads[end], Inf) < 10*eps()
sol_serial, _  = compute_solution(Sextic(), 4, 20, :GaussLegendre, :EC, :EC, 0.9, UniformMesh1D, Val{:serial}())
sol_threads, _ = compute_solution(Sextic(), 4, 20, :GaussLegendre, :EC, :EC, 0.9, UniformMesh1D, Val{:threads}())
@test norm(sol_serial[end] - sol_threads[end], Inf) < 10*eps()

mass_diff, ener_diff = mass_and_energy_differences(Septic(), :LobattoLegendre, 0.9)
@test abs(mass_diff) < 10*eps()
@test ener_diff < 5.e-8
mass_diff, ener_diff = mass_and_energy_differences(Septic(), :LobattoLegendre, 0.09)
@test abs(mass_diff) < 10*eps()
@test ener_diff < 5.e-12
mass_diff, ener_diff = mass_and_energy_differences(Septic(), :GaussLegendre, 0.9)
@test abs(mass_diff) < 10*eps()
@test ener_diff < 5.e-8
mass_diff, ener_diff = mass_and_energy_differences(Septic(), :GaussLegendre, 0.09)
@test abs(mass_diff) < 10*eps()
@test ener_diff < 5.e-12
sol_serial, _  = compute_solution(Septic(), 4, 20, :LobattoLegendre, :EC, :EC, 0.9, UniformPeriodicMesh1D, Val{:serial}())
sol_threads, _ = compute_solution(Septic(), 4, 20, :LobattoLegendre, :EC, :EC, 0.9, UniformPeriodicMesh1D, Val{:threads}())
@test norm(sol_serial[end] - sol_threads[end], Inf) < 10*eps()
sol_serial, _  = compute_solution(Septic(), 4, 20, :LobattoLegendre, :EC, :EC, 0.9, UniformMesh1D, Val{:serial}())
sol_threads, _ = compute_solution(Septic(), 4, 20, :LobattoLegendre, :EC, :EC, 0.9, UniformMesh1D, Val{:threads}())
@test norm(sol_serial[end] - sol_threads[end], Inf) < 10*eps()
sol_serial, _  = compute_solution(Septic(), 4, 20, :GaussLegendre, :EC, :EC, 0.9, UniformPeriodicMesh1D, Val{:serial}())
sol_threads, _ = compute_solution(Septic(), 4, 20, :GaussLegendre, :EC, :EC, 0.9, UniformPeriodicMesh1D, Val{:threads}())
@test norm(sol_serial[end] - sol_threads[end], Inf) < 10*eps()
sol_serial, _  = compute_solution(Septic(), 4, 20, :GaussLegendre, :EC, :EC, 0.9, UniformMesh1D, Val{:serial}())
sol_threads, _ = compute_solution(Septic(), 4, 20, :GaussLegendre, :EC, :EC, 0.9, UniformMesh1D, Val{:threads}())
@test norm(sol_serial[end] - sol_threads[end], Inf) < 10*eps()

mass_diff, ener_diff = mass_and_energy_differences(Octic(), :LobattoLegendre, 0.9)
@test abs(mass_diff) < 10*eps()
@test ener_diff < 5.e-8
mass_diff, ener_diff = mass_and_energy_differences(Octic(), :LobattoLegendre, 0.09)
@test abs(mass_diff) < 10*eps()
@test ener_diff < 5.e-12
mass_diff, ener_diff = mass_and_energy_differences(Octic(), :GaussLegendre, 0.9)
@test abs(mass_diff) < 10*eps()
@test ener_diff < 5.e-8
mass_diff, ener_diff = mass_and_energy_differences(Octic(), :GaussLegendre, 0.09)
@test abs(mass_diff) < 10*eps()
@test ener_diff < 5.e-12
sol_serial, _  = compute_solution(Octic(), 4, 20, :LobattoLegendre, :EC, :EC, 0.9, UniformPeriodicMesh1D, Val{:serial}())
sol_threads, _ = compute_solution(Octic(), 4, 20, :LobattoLegendre, :EC, :EC, 0.9, UniformPeriodicMesh1D, Val{:threads}())
@test norm(sol_serial[end] - sol_threads[end], Inf) < 10*eps()
sol_serial, _  = compute_solution(Octic(), 4, 20, :LobattoLegendre, :EC, :EC, 0.9, UniformMesh1D, Val{:serial}())
sol_threads, _ = compute_solution(Octic(), 4, 20, :LobattoLegendre, :EC, :EC, 0.9, UniformMesh1D, Val{:threads}())
@test norm(sol_serial[end] - sol_threads[end], Inf) < 10*eps()
sol_serial, _  = compute_solution(Octic(), 4, 20, :GaussLegendre, :EC, :EC, 0.9, UniformPeriodicMesh1D, Val{:serial}())
sol_threads, _ = compute_solution(Octic(), 4, 20, :GaussLegendre, :EC, :EC, 0.9, UniformPeriodicMesh1D, Val{:threads}())
@test norm(sol_serial[end] - sol_threads[end], Inf) < 10*eps()
sol_serial, _  = compute_solution(Octic(), 4, 20, :GaussLegendre, :EC, :EC, 0.9, UniformMesh1D, Val{:serial}())
sol_threads, _ = compute_solution(Octic(), 4, 20, :GaussLegendre, :EC, :EC, 0.9, UniformMesh1D, Val{:threads}())
@test norm(sol_serial[end] - sol_threads[end], Inf) < 10*eps()


text:
using BinaryProvider

# Parse some basic command-line arguments
const verbose = "--verbose" in ARGS
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))

# These are the two binary objects we care about
products = Product[
    LibraryProduct(prefix, "libfastcluster", :libfastcluster)
]

# Download binaries from hosted location
bin_prefix = "https://github.com/jmboehm/FastclusterBuilder/releases/download/0.0.1"
download_info = Dict(
    Linux(:aarch64, :glibc) => ("$bin_prefix/libfastcluster.v0.0.1.aarch64-linux-gnu.tar.gz", "74f873d6120ca95d7a4ef500c1f17d8636c5f1165452e4efccb8b66db0b0d921"),
    Linux(:armv7l, :glibc)  => ("$bin_prefix/libfastcluster.v0.0.1.arm-linux-gnueabihf.tar.gz", "04708fc8f5e5cb9b5432a842d1c29f3a8739445f1613655b31350f7a461590b4"),
    Linux(:powerpc64le, :glibc) => ("$bin_prefix/libfastcluster.v0.0.1.powerpc64le-linux-gnu.tar.gz", "9c5ff18272a241df216885a45b76ec7b535ff0c9a0314ce7d8052ac7288f638f"),
    Linux(:i686, :glibc)    => ("$bin_prefix/libfastcluster.v0.0.1.i686-linux-gnu.tar.gz", "ca945df0fe0248c7e84ff0035201de832e46941a67c1deaff38ad7411986e9f6"),
    Linux(:x86_64, :glibc)  => ("$bin_prefix/libfastcluster.v0.0.1.x86_64-linux-gnu.tar.gz", "1762bc297a6f2bd9eafd3e3da2d9ba5e56b502b6834d1a861fc0eee8a0be77d5"),

    Linux(:aarch64, :musl)  => ("$bin_prefix/libfastcluster.v0.0.1.aarch64-linux-musl.tar.gz", "fab9df6f18f62550ea85690e3e3070de7fe9f9f6c991198953dae75d350ac878"),
    Linux(:armv7l, :musl)   => ("$bin_prefix/libfastcluster.v0.0.1.arm-linux-musleabihf.tar.gz", "9af5b1631a5e8ecf33dcf660d79270c7fe428ea6d3614e9d15a3c52f4bf16b70"),
    Linux(:i686, :musl)     => ("$bin_prefix/libfastcluster.v0.0.1.i686-linux-musl.tar.gz", "b167e77606216af5baaa58726b314b930944c1bb169c83dc93861eaa62a2adf8"),
    Linux(:x86_64, :musl)   => ("$bin_prefix/libfastcluster.v0.0.1.x86_64-linux-musl.tar.gz", "6d4f7d89656aeac60e889b470bd39986ede9d88d1ea351b02da5822f5dbbc55a"),

    FreeBSD(:x86_64)        => ("$bin_prefix/libfastcluster.v0.0.1.x86_64-unknown-freebsd11.1.tar.gz", "2d70eae69980f537abfe0395e5539ad46a0db3c5f402236931677829065ee55e"),
    MacOS(:x86_64)          => ("$bin_prefix/libfastcluster.v0.0.1.x86_64-apple-darwin14.tar.gz", "7154e6d81bf91df3e59e2d2e76243b62e098d4b45a8e4df1d5bf14d6d6857c06"),

    Windows(:i686)          => ("$bin_prefix/libfastcluster.v0.0.1.i686-w64-mingw32.tar.gz", "3926031d25e52f894ef89b5f65de4b9bc57aa09c80d2f6e9824353c5c6429c9d"),
    Windows(:x86_64)        => ("$bin_prefix/libfastcluster.v0.0.1.x86_64-w64-mingw32.tar.gz", "3926031d25e52f894ef89b5f65de4b9bc57aa09c80d2f6e9824353c5c6429c9d"),
)
# To get hashes:
# function url2hash(url::String)
#     path = download(url)
#     open(io-> bytes2hex(BinaryProvider.sha256(io)), path)
# end
# for i in keys(download_info)
#         @show i
#         @show url2hash(download_info[i][1])
# end

# First, check to see if we're all satisfied
if any(!satisfied(p; verbose=verbose) for p in products)
    try
        # Download and install binaries
        url, tarball_hash = choose_download(download_info)
        install(url, tarball_hash; prefix=prefix, force=true, verbose=true)
    catch e
        if typeof(e) <: ArgumentError
            error("Your platform $(Sys.MACHINE) is not supported by this package!")
        else
            rethrow(e)
        end
    end

    # Finally, write out a deps.jl file
    # I'm commenting this out, since our library seems to have some
    # issues with Libdl.dlopen(). But we don't really need it anyway.
    #write_deps_file(joinpath(@__DIR__, "deps.jl"), products)
end


text:
using AtomicGraphNets
using Test

@testset "AtomicGraphNets.jl" begin

    @testset "Layers" begin
        include("layer_tests.jl")
    end

    @testset "Models" begin
        include("model_tests.jl")
    end

    @testset "Examples" begin
        include("example_tests.jl")
    end

end


text:
using ParallelUtils
using Base.Test

@ensureworkers(5)

a = 1
sendto(2, a2=a)
@test getfrom(2, :a2) == 1

max_stars_repo_name: slundberg/ParallelUtils

max_issues_repo_name: slundberg/ParallelUtils

max_forks_repo_name: slundberg/ParallelUtils


text:
mutable struct MatrixDataIntegrator{yType, uType, lType, rType}
    Δy::yType
    u::uType
    left_factor::lType
    right_factor::rType
    sign::Int
end

function step!(integrator::MatrixDataIntegrator, dt, ::Bool=true)
    step!(integrator)
end

function step!(integrator::MatrixDataIntegrator)
    @unpack Δy, u, left_factor, right_factor, sign = integrator
    u .+= sign*left_factor'*Δy*right_factor
end

function set_u!(integrator::MatrixDataIntegrator, unew)
    integrator.u .= unew
end

function update_data!(x,y::Function,t::Number,dt::Number)
    x .= y(t+dt)
end

function update_data!(x,y::AbstractArray,t::Int,dt::Int)
    x .= deepcopy(y[t + dt])
end


text:
using WordCloud
using Test
using Random

include("test_render.jl")
include("test_textprocessing.jl")

@testset "WordCloud.jl" begin
    # @show pwd()
    # overall test
    wc = runexample(:random)
    @test getstate(wc) == :generate!
    @test isempty(WordCloud.outofbounds(wc.maskqtree, wc.qtrees))
    paint(wc)
    paint(wc, "test.jpg", background=outline(wc.mask, color=(1, 0, 0.2, 0.7), linewidth=2), ratio=0.5)
    paint(wc, "test.svg")
    @test getparameter(wc, :maskoccupying) == WordCloud.occupying(WordCloud.QTree.kernel(wc.maskqtree[1]), WordCloud.QTree.FULL)

    # placewords!
    placewords!(wc, style=:gathering)
    words = ["." for i in 1:500]
    weights = [1 for i in 1:length(words)]
    @test_throws ErrorException begin #no room
        wc = wordcloud(words, weights, maskshape=ellipse, masksize=(5,5), backgroundsize=(10,10), density=1000, angles=0)
        placewords!(wc)
    end

    # wordcloud factory
    wc = wordcloud(["singleword"=>12], mask=shape(box, 200, 150, 40, color=0.15), density=0.45, run=generate!) #singleword & Pair
    wc = wordcloud("giving a single word is ok. giving several words is ok too",
            mask=shape(squircle, 200, 150, color=0.15, rt=2.2), density=0.45, transparent=(1,1,1,0)) #String & small mask
    @test_throws AssertionError wordcloud(["1"],[2,3], density=0.1)|>generate! #length unmatch
    @test_throws AssertionError wordcloud(String[],Int[], density=0.1)|>generate! #empty inputs
    ##############no mask file
    wc = wordcloud(["test"], [1], maskcolor="green", outline=5)
    @test WordCloud.alpha(parsecolor(getbackgroundcolor(wc))) == 0
    wc = wordcloud(["test"], [1], backgroundcolor="blue", outline=5)
    @test parsecolor(getmaskcolor(wc)) == parsecolor("blue")
    wc = wordcloud(["test"], [1], maskcolor="green")
    @test getparameter(wc, :outline) == 0
    wc = wordcloud(["test"], [1], backgroundcolor="blue")
    @test getparameter(wc, :outline) == 0
    wc = wordcloud(["test"], [1], masksize=(100,100), outline=0)
    @test all(size(wc.mask) .> 105)
    wc = wordcloud(["test"], [1], masksize=(100,100), outline=30, padding=30)
    @test all(size(wc.mask) .> 200)
    ##############svg mask
    svgfile = "test.svg"
    wordcloud(["test"], [1], colors="#DE2910", mask=svgfile, maskcolor=:original)
    wordcloud(["test"], [1], mask=open(svgfile))
    wc = wordcloud(["test"], [1], mask=svgfile, backgroundcolor=0) #warning#can't edit the svg to remove original backgroundcolor,
    #so it's only work when the svgfile has a transparent background
    wc2 = wordcloud(["test"], [1], mask=open(svgfile), padding=20)
    wc3 = wordcloud(["test"], [1], mask=open(svgfile), padding=(10,-10))
    @test all(size(wc2.mask) .> size(wc.mask))
    @test all(size(wc2.svgmask) .> size(wc.svgmask))
    @test (size(wc3.svgmask) .> size(wc.svgmask)) == (true, false)
    wordcloud(["test"], [1], mask=open(svgfile), padding=10, backgroundcolor="red")#warning#
    ##############png mask
    pngfile = pkgdir(WordCloud)*"/res/heart_mask.png"
    wordcloud(["test"], [1], colors="#DE2910", mask=pngfile, maskcolor=:original)
    wc = wordcloud(["test"], [1], mask=open(pngfile), maskcolor="yellow", ratio=0.5)
    @test getbackgroundcolor(wc) in WordCloud.DEFAULTSYMBOLS
    wc = wordcloud(["test"], [1], mask=pngfile, maskcolor="green", outline=5)
    @test getbackgroundcolor(wc) in WordCloud.DEFAULTSYMBOLS
    @test getparameter(wc, :outline) == 5
    wc = wordcloud(["test"], [1], colors="#DE2910", mask=pngfile, backgroundcolor=0)
    @test getmaskcolor(wc) in WordCloud.DEFAULTSYMBOLS
    wc = wordcloud(["test"], [1], colors="#DE2910", mask=pngfile, maskcolor=1, backgroundcolor=0)
    @test parsecolor(getmaskcolor(wc)) == parsecolor(1)
    @test parsecolor(getbackgroundcolor(wc)) == parsecolor(0)
    wordcloud(["test"], [1], mask=pngfile, backgroundcolor=:maskcolor) #backgroundcolor=maskcolor=:default, didn't change anything
    wordcloud(["test"], [1], mask=pngfile, maskcolor="#faeef8", backgroundcolor=:maskcolor)
    wordcloud(["test"], [1], mask=pngfile, maskcolor="#faeef8", backgroundcolor=:maskcolor, outline=3)
    wordcloud(["test"], [1], mask=pngfile, backgroundcolor=:maskcolor, outline=3, smoothness=0.7)
    wc = wordcloud(["test"], [1], mask=open(pngfile), backgroundcolor=:auto)
    @test getmaskcolor(wc) == getbackgroundcolor(wc)
    wc = wordcloud(["test"], [1], mask=open(pngfile), maskcolor=:auto)
    @test getbackgroundcolor(wc) == :default
    @test size(wc.mask) == (572, 640)
    wc = wordcloud(["test"], [1], mask=pngfile, masksize=(200,200))
    @test size(wc.mask) == (200,200)
    wc = wordcloud(["test"], [1], mask=pngfile, ratio=0.3)
    @test all(size(wc.mask) .< 200)
    wc2 = wordcloud(["test"], [1], mask=pngfile, outline=50, ratio=0.3)
    @test size(wc2.mask) == size(wc.mask)
    wc3 = wordcloud(["test"], [1], mask=pngfile, padding=(-50, 50))
    @test (size(wc3.mask) .> (572, 640)) == (false, true)
    wc = wordcloud(["test"], [1], mask=pngfile, padding=100)
    @test all(size(wc.mask) .> 700)
    # get&set
    wc = wordcloud(
            processtext(open("../res/alice.txt"), stopwords=WordCloud.stopwords_en ∪ ["said"], maxnum=300),
            mask = "../res/alice_mask.png", maskcolor = "#faeef8", backgroundcolor = 0.97,
            colors = (WordCloud.colorschemes[:Set1_5].colors..., ),
            angles = (0, 90));
    rescale!(wc, 1.23)
    pin(wc, ["little", "know"]) do
        @test length(wc.words)==298
        setpositions!(wc, 1, (2,2))
        setpositions!(wc, [1, "Alice", "one"], (-1, -2))
        setpositions!(wc, [1, "Alice", "one"], [(10,10),(10,20),(21,2)])
        setpositions!(wc, "time", (0,0), type=setcenter!)
    end
    @test getpositions(wc, [1, "Alice", "one"])[3] == (21,2)
    @test WordCloud.QTree.kernelsize(wc.qtrees[WordCloud.index(wc, "time")]) == size(getimages(wc, "time"))
    @test .-reverse(size(getimages(wc, "time"))) .÷ 2 == getpositions(wc, ["time", getwords(wc, 9)])[1]
    w = getweights(wc, getwords(wc, [1,2]))
    setwords!(wc, [1,2], ["zz","yy"])
    @test getweights(wc, "zz") == w[1]
    setimages!(wc, [1,2], wc.imgs[[4,5]])
    setimages!(wc, 1, wc.imgs[[4,5]])
    setimages!(wc, 1, wc.imgs[4])
    setsvgimages!(wc, 1, wc.svgs[6])
    setsvgimages!(wc, 6, wc.svgs[6]) #the results of setsvgimages! and initword! may not be identical
    @test wc.imgs[1] == wc.imgs[6]

    for s = [:reset, :average, :clipping, :blending, :reset]
        recolor!(wc, style=s)
    end

    #utils
    wc.qtrees[1][1]|>imageof
    bg = getmask(wc)
    istrans = c->maximum(c[1:3])<128
    mask = WordCloud.imagemask(bg, istrans)
    s = showmask(bg, mask)
    @test all(bg[mask] .== s[mask])
    @test all(bg[.!mask] .!= s[.!mask])
end


text:
# embedding_map.jl

"A map between embedded spaces."
struct EmbeddingMap{S,T} <: AbstractMap{S,T}
    function EmbeddingMap{S,T}() where {S,T}
        @assert embedded(spacetype(S), spacetype(T))
        new{S,T}()
    end
end

embedding_map(::Type{S}, ::Type{T}) where {S,T} = EmbeddingMap{S,T}()
embedding_map(::Type{T}, ::Type{T}) where {T} = IdentityMap{T}()

(m::EmbeddingMap)(x) = applymap(m, x)

applymap(map::EmbeddingMap{S,T}, x::S) where {S,T} = convert_space(spacetype(T), x)

"A restriction map from a space to an embedded space."
struct RestrictionMap{S,T} <: AbstractMap{S,T}
    function RestrictionMap{S,T}() where {S,T}
        @assert embedded(spacetype(T), spacetype(S))
        new{S,T}()
    end
end

restriction_map(::Type{S}, ::Type{T}) where {S,T} = RestrictionMap{S,T}()
restriction_map(::Type{T}, ::Type{T}) where {T} = IdentityMap{T}()

applymap(map::RestrictionMap{S,T}, x::S) where {S,T} = restrict_space(spacetype(T), x)

left_inverse(map::EmbeddingMap{S,T}) where {S,T} = RestrictionMap{T,S}()
right_inverse(map::RestrictionMap{S,T}) where {S,T} = EmbeddingMap{T,S}()

"A map between isomorphic spaces."
struct IsomorphismMap{S,T} <: AbstractMap{S,T}
    function IsomorphismMap{S,T}() where {S,T}
        @assert embedded(spacetype(T), spacetype(S))
        new{S,T}()
    end
end

isomorphism_map(::Type{S}, ::Type{T}) where {S,T} = IsomorphismMap{S,T}()
isomorphism_map(::Type{T}, ::Type{T}) where {T} = IdentityMap{T}()

applymap(map::IsomorphismMap{S,T}, x::S) where {S,T} = convert_space(spacetype(T), x)

apply_inverse(map::IsomorphismMap{S,T}, y::T) where {S,T} = convert_space(spacetype(S), y)

inv(map::IsomorphismMap{S,T}) where {S,T} = IsomorphismMap{T,S}()


text:
@testset "read_graph" begin

    using DelimitedFiles

    @testset "simple" begin
        folder = "data/networks/simple"
        for fname in readdir(folder)
            file = joinpath(folder, fname)
            graph = read_graph(file, "simple")
            data, header = readdlm(file, '\t', Int64, header=true)
            @test length(graph.edges) == size(data, 1)
            @test length(graph.first_out) == length(graph.nodes) + 1
            @test length(graph.end_node) == length(graph.edges)
        end
    end

    @testset "distance" begin
        folder = "data/networks/distance"
        for fname in readdir(folder)
            file = joinpath(folder, fname)
            graph = read_graph(file, "distance")
            data, header = readdlm(file, '\t', Int64, header=true)
            @test length(graph.edges) == size(data, 1)
            @test length(graph.first_out) == length(graph.nodes) + 1
            @test length(graph.end_node) == length(graph.edges)
        end
    end

    @testset "capacity" begin
        folder = "data/networks/capacity"
        for fname in readdir(folder)
            file = joinpath(folder, fname)
            graph = read_graph(file, "capacity")
            data, header = readdlm(file, '\t', Int64, header=true)
            @test length(graph.edges) == size(data, 1)
            @test length(graph.first_out) == length(graph.nodes) + 1
            @test length(graph.end_node) == length(graph.edges)
        end
    end

    @testset "assignment" begin
        folder = "data/networks/assignment"
        for fname in readdir(folder)
            file = joinpath(folder, fname)
            graph = read_graph(file, "assignment")
            data, header = readdlm(file, '\t', Int64, header=true)
            @test length(graph.edges) == size(data, 1)
            @test length(graph.first_out) == length(graph.nodes) + 1
            @test length(graph.end_node) == length(graph.edges)
        end
    end

end


text:
using GLMakie
using JLD2

filename = "bathymetry-360x150-latitude-75.0.jld2"
file = jldopen(filename)
h = file["bathymetry"]
close(file)
h[h .> 0] .= NaN

fig = Figure(resolution=(1800, 800))
ax = Axis(fig[1, 1], xlabel="Longitude", ylabel="Latitude")
hm = heatmap!(ax, h)
Colorbar(fig[1, 2], hm, label="Elevation (m)")
display(fig)


text:
function minimumvertex!(g, distance)
    let i = 0, v = 0
        while i < countvertices(g)
            i += 1
            if getmark(g, i) == 0
                v = i
                break
            end
        end
        while i < countvertices(g)
            i += 1
            if getmark(g, i) == 0 && distance[i] < distance[v]
                v = i
            end
        end
        return v
    end
end


text:
using EchelleCCFs
using Test

    @testset "CCF shape constructors" begin
        Δv = 7e3
        c = EchelleCCFs.speed_of_light_mps
        Δz = Δv/c
        λ = 5000.0

        import EchelleCCFs: λ_min, λ_max, integrate
        @testset "Tophat mask" begin
            @test_nowarn TopHatCCFMask(Δv)
            m = TopHatCCFMask(Δv)
            @test_nowarn λ_min(m,λ)
            @test_nowarn λ_max(m,λ)
            @test_nowarn integrate(m, (λ_min(m,λ)-λ)*c/λ,  (λ_max(m,λ)-λ)*c/λ)
            @test integrate(m, (λ_min(m,λ)-λ)*c/λ,  (λ_max(m,λ)-λ)*c/λ) ≈ 1    rtol = 1e-3
        end
        @testset "Half Cos mask" begin
            @test_nowarn CosCCFMask(Δv)
            m = CosCCFMask(Δv)
            @test_nowarn λ_min(m,λ)
            @test_nowarn λ_max(m,λ)
            @test_nowarn integrate(m, (λ_min(m,λ)-λ)*c/λ,  (λ_max(m,λ)-λ)*c/λ)
            @test integrate(m, (λ_min(m,λ)-λ)*c/λ,  (λ_max(m,λ)-λ)*c/λ) ≈ 1    rtol = 1e-3
        end
        @testset "Gaussian mask" begin
            @test_nowarn GaussianCCFMask(Δv)
            m = GaussianCCFMask(Δv)
            @test_nowarn λ_min(m,λ)
            @test_nowarn λ_max(m,λ)
            @test_nowarn integrate(m, (λ_min(m,λ)-λ)*c/λ,  (λ_max(m,λ)-λ)*c/λ)
            @test integrate(m, (λ_min(m,λ)-λ)*c/λ,  (λ_max(m,λ)-λ)*c/λ) ≈ 1    rtol = 1e-3
        end
        @testset "SuperGaussian mask" begin
            @test_nowarn SuperGaussianCCFMask(Δv,1.1)
            m = SuperGaussianCCFMask(Δv,1.1)
            @test_nowarn λ_min(m,λ)
            @test_nowarn λ_max(m,λ)
            @test_nowarn integrate(m, (λ_min(m,λ)-λ)*c/λ,  (λ_max(m,λ)-λ)*c/λ)
            @test integrate(m, (λ_min(m,λ)-λ)*c/λ,  (λ_max(m,λ)-λ)*c/λ) ≈ 1    rtol = 1e-3
        end
        #=  TODO: Test before using
        @testset "Gaussian Mixture mask" begin
            @test_nowarn GaussianMixtureCCFMask([0.75, 0.25], [Δv, 0.5*Δv], 2*Δv)
            m = GaussianMixtureCCFMask([0.75, 0.25], [Δv, 0.5*Δv], 2*Δv)
            @test_nowarn λ_min(m,λ)
            @test_nowarn λ_max(m,λ)
            @test_nowarn integrate(m, (λ_min(m,λ)-λ)*c/λ,  (λ_max(m,λ)-λ)*c/λ)
            @test integrate(m, (λ_min(m,λ)-λ)*c/λ,  (λ_max(m,λ)-λ)*c/λ) ≈ 1    rtol = 1e-3
        end
        =#
    end

max_stars_repo_name: RvSpectML/EchelleCCFs

max_issues_repo_name: RvSpectML/EchelleCCFs

max_forks_repo_name: RvSpectML/EchelleCCFs


text:
import Pkg
Pkg.activate(".")
using Test
using PandaModels; const _PdM = PandaModels
import PowerModels; const _PM = PowerModels

_PM.silence()

pdm_path = joinpath(dirname(pathof(PandaModels)), "..")
data_path = joinpath(pdm_path, "test", "data")
cast_ts = joinpath(data_path, "cigre_with_timeseries.json")

##
pm2 = _PdM.load_pm_from_json(cast_ts)
param = _PdM.extract_params!(pm2)

##
pm = _PdM.load_pm_from_json(cast_ts)
_PdM.active_powermodels_silence!(pm)
pm = _PdM.check_powermodels_data!(pm)
# params = _PdM.extract_params!(pm)
# pm = delete!(pm, "user_defined_params")
model = _PdM.get_model(pm["pm_model"])
solver = _PdM.get_solver(pm)
mn = _PdM.set_pq_values_from_timeseries(pm)
# pm_mn = _PM.instantiate_model(mn, model, _PM.build_mn_opf, multinetwork=true)

result = _PdM._run_v_stab_ts(
    mn,
    model,
    solver,
    setting = Dict("output" => Dict("branch_flows" => true)),
    ext = _PdM.extract_params!(pm),
)

# using InfrastructureModels; const _IM = InfrastructureModels
#
# for (n, network) in _IM.nws(mn)
#         variable_bus_voltage(pm, nw=n)
#         variable_gen_power(pm, nw=n)
#         variable_branch_power(pm, nw=n)
#         variable_dcline_power(pm, nw=n)
#
#         constraint_model_voltage(pm, nw=n)
#
#         for i in ids(pm, :ref_buses, nw=n)
#             constraint_theta_ref(pm, i, nw=n)
#         end
#
#         for i in ids(pm, :bus, nw=n)
#             constraint_power_balance(pm, i, nw=n)
#         end
#
#         for i in ids(pm, :branch, nw=n)
#             constraint_ohms_yt_from(pm, i, nw=n)
#             constraint_ohms_yt_to(pm, i, nw=n)
#
#             constraint_voltage_angle_difference(pm, i, nw=n)
#
#             constraint_thermal_limit_from(pm, i, nw=n)
#             constraint_thermal_limit_to(pm, i, nw=n)
#         end
#
#         for i in ids(pm, :dcline, nw=n)
#             constraint_dcline_power_losses(pm, i, nw=n)
#         end
#     end

# _PdM.active_powermodels_silence!(pm)
# pm = _PdM.check_powermodels_data!(pm)
# model = _PdM.get_model(pm["pm_model"])
# solver = _PdM.get_solver(pm)
# result = _PdM._run_q_flex(
#     pm,
#     model,
#     solver,
#     setting = Dict("output" => Dict("branch_flows" => true)),
#     ext = _PdM.extract_params!(pm),
# )
# mn = set_pq_values_from_timeseries(pm)

# result = _run_mn_vd(
#     pm,
#     model,
#     solver,
#     setting = Dict("output" => Dict("branch_flows" => true)),
#     ext = extract_params!(pm),
# )

# result = run_pandamodels_q_flex(case_q_flex)
# import InfrastructureModels
# using JuMP
# pm = _PdM.load_pm_from_json(case_q_flex)
# # active_powermodels_silence!(pm)
# pm =  _PdM.check_powermodels_data!(pm)
# model =  _PdM.get_model(pm["pm_model"])
# solver =  _PdM.get_solver(pm)
# #
# # # ext = _PdM.extract_params!(pm)
# #
# # result = _run_q_flex(
# # pm,
# # model,
# # solver,
# # setting = Dict("output" => Dict("branch_flows" => true)),
# # ext = _PdM.extract_params!(pm),
# # )
#
#
# # pm = load_pm_from_json(json_path)
# # active_powermodels_silence!(pm)
# # pm = check_powermodels_data!(pm)
# # model = get_model(pm["pm_model"])
# # solver = get_solver(pm)
# ##  Load System Data
# # pm = load_pm_from_json(json_path)
# # solver = get_solver(pm["pm_solver"], pm["pm_nl_solver"], pm["pm_mip_solver"],
# # pm["pm_log_level"], pm["pm_time_limit"], pm["pm_nl_time_limit"], pm["pm_mip_time_limit"])
#
# # if haskey(pm, "user_defined_params")
#     # params = Dict{Symbol,Dict{String,Any}}()
#     # for key in keys(pm["user_defined_params"])
#     #     params[Symbol(key)] = pm["user_defined_params"][key]
#     # end
#     # data = delete!(pm, "user_defined_params")
# # end
# #
# if haskey(pm, "user_defined_params")
#     user_defined_params = pm["user_defined_params"]
#     data = delete!(pm, "user_defined_params")
# else
#     data = pm
# end
#
# # additional modification
# for (i, branch) in pm["branch"]
#     pm["branch"][i]["angmin"] = -6.28
#     pm["branch"][i]["angmax"] = 6.28
# end
#
# # use build_ref to filter out inactive components
# ref = PowerModels.build_ref(data)[:it][:pm][:nw][0]
#
# ###############################################################################
# # 1. Building the Optimal Power Flow Model
# ###############################################################################
# # Initialize a JuMP Optimization Model
# #-------------------------------------
# model = Model(solver)
# #nlp_optimizer = optimizer_with_attributes(Ipopt.Optimizer)
# #model = Model(nlp_optimizer)
# # Add Optimization and State Variables
# # ------------------------------------
#
# # Add voltage angles va for each bus
# @variable(model, va[i in keys(ref[:bus])])
# # note: [i in keys(ref[:bus])] adds one `va` variable for each bus in the network
#
# # Add voltage angles vm for each bus
# @variable(model, ref[:bus][i]["vmin"] <= vm[i in keys(ref[:bus])] <= ref[:bus][i]["vmax"], start=1.0)
# # note: this vairable also includes the voltage magnitude limits and a starting value
#
# # Add active power generation variable pg for each generator (including limits)
# @variable(model, ref[:gen][i]["pmin"] <= pg[i in keys(ref[:gen])] <= ref[:gen][i]["pmax"])
# # Add reactive power generation variable qg for each generator (including limits)
# @variable(model, ref[:gen][i]["qmin"] <= qg[i in keys(ref[:gen])] <= ref[:gen][i]["qmax"])
#
# # Add power flow variables p to represent the active power flow for each branch
# @variable(model, -ref[:branch][l]["rate_a"] <= p[(l,i,j) in ref[:arcs]] <= ref[:branch][l]["rate_a"])
# # Add power flow variables q to represent the reactive power flow for each branch
# @variable(model, -ref[:branch][l]["rate_a"] <= q[(l,i,j) in ref[:arcs]] <= ref[:branch][l]["rate_a"])
# # note: ref[:arcs] includes both the from (i,j) and the to (j,i) sides of a branch
#
# # Add Objective Function
# # ----------------------
#
# @objective(model, Min, sum((q[(content["element_index"],
#                             content["f_bus"],
#                             content["t_bus"])] - content["value"])^2
#                         for (i, content) in user_defined_params["setpoint_q"]))
#
# # Add Constraints
# # ---------------
#
# # Fix the voltage angle to zero at the reference bus ## why?
# for (i,bus) in ref[:ref_buses]
#     @constraint(model, va[i] == 0)
# end
#
# # Nodal power balance constraints
# for (i,bus) in ref[:bus]
#     # Build a list of the loads and shunt elements connected to the bus i
#     bus_loads = [ref[:load][l] for l in ref[:bus_loads][i]]
#     bus_shunts = [ref[:shunt][s] for s in ref[:bus_shunts][i]]
#
#     # Active power balance at node i
#     @constraint(model,
#         sum(p[a] for a in ref[:bus_arcs][i]) +                  # sum of active power flow on lines from bus i +
#         sum(p_dc[a_dc] for a_dc in ref[:bus_arcs_dc][i]) ==     # sum of active power flow on HVDC lines from bus i =
#         sum(pg[g] for g in ref[:bus_gens][i]) -                 # sum of active power generation at bus i -
#         sum(load["pd"] for load in bus_loads) -                 # sum of active load consumption at bus i -
#         sum(shunt["gs"] for shunt in bus_shunts)*vm[i]^2        # sum of active shunt element injections at bus i
#     )
#
#     # Reactive power balance at node i
#     @constraint(model,
#         sum(q[a] for a in ref[:bus_arcs][i]) +                  # sum of reactive power flow on lines from bus i +
#         sum(q_dc[a_dc] for a_dc in ref[:bus_arcs_dc][i]) ==     # sum of reactive power flow on HVDC lines from bus i =
#         sum(qg[g] for g in ref[:bus_gens][i]) -                 # sum of reactive power generation at bus i -
#         sum(load["qd"] for load in bus_loads) +                 # sum of reactive load consumption at bus i -
#         sum(shunt["bs"] for shunt in bus_shunts)*vm[i]^2        # sum of reactive shunt element injections at bus i
#     )
# end
#
# # Branch power flow physics and limit constraints
# for (i,branch) in ref[:branch]
#     # Build the from variable id of the i-th branch, which is a tuple given by (branch id, from bus, to bus)
#     f_idx = (i, branch["f_bus"], branch["t_bus"])
#     # Build the to variable id of the i-th branch, which is a tuple given by (branch id, to bus, from bus)
#     t_idx = (i, branch["t_bus"], branch["f_bus"])
#     # note: it is necessary to distinguish between the from and to sides of a branch due to power losses
#
#     p_fr = p[f_idx]                     # p_fr is a reference to the optimization variable p[f_idx]
#     q_fr = q[f_idx]                     # q_fr is a reference to the optimization variable q[f_idx]
#     p_to = p[t_idx]                     # p_to is a reference to the optimization variable p[t_idx]
#     q_to = q[t_idx]                     # q_to is a reference to the optimization variable q[t_idx]
#     # note: adding constraints to p_fr is equivalent to adding constraints to p[f_idx], and so on
#
#     vm_fr = vm[branch["f_bus"]]         # vm_fr is a reference to the optimization variable vm on the from side of the branch
#     vm_to = vm[branch["t_bus"]]         # vm_to is a reference to the optimization variable vm on the to side of the branch
#     va_fr = va[branch["f_bus"]]         # va_fr is a reference to the optimization variable va on the from side of the branch
#     va_to = va[branch["t_bus"]]         # va_fr is a reference to the optimization variable va on the to side of the branch
#
#     # Compute the branch parameters and transformer ratios from the data
#     g, b = PowerModels.calc_branch_y(branch)
#     tr, ti = PowerModels.calc_branch_t(branch)
#     g_fr = branch["g_fr"]
#     b_fr = branch["b_fr"]
#     g_to = branch["g_to"]
#     b_to = branch["b_to"]
#     tm = branch["tap"]^2
#     # note: tap is assumed to be 1.0 on non-transformer branches
#
#
#     # AC Power Flow Constraints
#
#     # From side of the branch flow
#     @NLconstraint(model, p_fr ==  (g+g_fr)/tm*vm_fr^2 + (-g*tr+b*ti)/tm*(vm_fr*vm_to*cos(va_fr-va_to)) + (-b*tr-g*ti)/tm*(vm_fr*vm_to*sin(va_fr-va_to)) )
#     @NLconstraint(model, q_fr == -(b+b_fr)/tm*vm_fr^2 - (-b*tr-g*ti)/tm*(vm_fr*vm_to*cos(va_fr-va_to)) + (-g*tr+b*ti)/tm*(vm_fr*vm_to*sin(va_fr-va_to)) )
#
#     # To side of the branch flow
#     @NLconstraint(model, p_to ==  (g+g_to)*vm_to^2 + (-g*tr-b*ti)/tm*(vm_to*vm_fr*cos(va_to-va_fr)) + (-b*tr+g*ti)/tm*(vm_to*vm_fr*sin(va_to-va_fr)) )
#     @NLconstraint(model, q_to == -(b+b_to)*vm_to^2 - (-b*tr+g*ti)/tm*(vm_to*vm_fr*cos(va_fr-va_to)) + (-g*tr-b*ti)/tm*(vm_to*vm_fr*sin(va_to-va_fr)) )
#
#     # Voltage angle difference limit
#     @constraint(model, va_fr - va_to <= branch["angmax"])
#     @constraint(model, va_fr - va_to >= branch["angmin"])
#
#     # Apparent power limit, from side and to side
#     @constraint(model, p_fr^2 + q_fr^2 <= branch["rate_a"]^2)
#     @constraint(model, p_to^2 + q_to^2 <= branch["rate_a"]^2)
# end
#
#
# ###############################################################################
# # 3. Solve the Optimal Power Flow Model and Review the Results
# ###############################################################################
#
# # Solve the optimization problem
# optimize!(model)
# ###############################################################################
# # 4. Create Result Dictionary such that the PowerModels Results can be used by pandapower
# ###############################################################################
# solution =  Dict{String,Any}()
# push!(solution, "baseMVA" => data["baseMVA"])
# push!(solution, "per_unit" => data["per_unit"])
# push!(solution, "gen" => data["gen"])
# push!(solution, "bus" => data["bus"])
# push!(solution, "branch" => data["branch"])
# push!(solution, "multinetwork" => false)
# push!(solution, "multiinfrastructrue" => false)
#
# for (i, gen) in solution["gen"]
#     index = gen["index"]
#     gen["qg"] = value(model[:qg][index])
#     gen["pg"] = value(model[:pg][index])
# end
#
# for (i, bus) in solution["bus"]
#     index = bus["index"]
#     bus["vm"] = value(model[:vm][index])
#     bus["va"] = value(model[:va][index])
# end
#
# for (i, branch) in solution["branch"]
#     index = branch["index"]
#     push!(branch, "qf" => value(model[:q][(index, branch["f_bus"], branch["t_bus"])]))
#     push!(branch, "qt" => value(model[:q][(index, branch["t_bus"], branch["f_bus"])]))
#     push!(branch, "pf" => value(model[:p][(index, branch["f_bus"], branch["t_bus"])]))
#     push!(branch, "pt" => value(model[:p][(index, branch["t_bus"], branch["f_bus"])]))
# end
# println("The solver termination status is")
#
# result = Dict{String,Any}(
#         "optimizer" => JuMP.solver_name(model),
#         "termination_status" => JuMP.termination_status(model),   # "LOCALLY_SOLVED",
#         "primal_status" => JuMP.primal_status(model),
#         "dual_status" => JuMP.dual_status(model),
#         "objective" => InfrastructureModels._guard_objective_value(model),
#         "objective_lb" => InfrastructureModels._guard_objective_bound(model),
#         "solve_time" => solve_time,
#         "solution" => solution)
#
# return result


text:
## Kernel types

# abstract kernel supertype, all others inherit from Kernel
abstract type AbstractKernel end

# all subtypes characterize similarity calculations
struct InverseDistanceKernel <: AbstractKernel end

struct GaussianKernel <: AbstractKernel
    σ::Float64
end
GaussianKernel() = GaussianKernel(1.0)

struct CosineKernel <: AbstractKernel end

struct CustomKernel <: AbstractKernel
    func::Function
end

## similarity computation for each kernel

# multiple dispatch over kernel types
@inline function similarity(k::InverseDistanceKernel, x::AbstractArray{T, N}, y::AbstractArray{T, N}) where {T,N}
    ret = zero(T)
    m = length(x)
    @inbounds @simd for k in 1:m
        ret += (x[k] - y[k])^2
    end
    return ret == 0.0 ? 0.0 : 1 / sqrt(ret)
end

@inline function similarity(k::GaussianKernel, x::AbstractArray{T, N}, y::AbstractArray{T, N}) where {T,N}
    ret = zero(T)
    m = length(x)
    @inbounds @simd for k in 1:m
        ret += (x[k] - y[k])^2
    end
    return ret == 0.0 ? 0.0 : exp(-ret /(2 * k.σ^2))
end

function similarity(k::CosineKernel, x::AbstractArray{T, N}, y::AbstractArray{T, N}) where {T,N}
    return (x'*y)/(norm(x)*norm(y))
end

function similarity(k::CustomKernel, x::AbstractArray{T, N}, y::AbstractArray{T, N}) where {T,N}
    return k.func(x,y)
end


text:
using Cairo
c  = CairoRGBSurface(512, 128);
cr = CairoContext(c);
save(cr);
set_source_rgb(cr, 0.8, 0.8, 0.8);  	# light gray
rectangle(cr, 0.0, 0.0, 512.0, 128.0); 	# background
fill(cr);
restore(cr);
save(cr);

x0=51.2;  y0=64.0;
x1=204.8; y1=115.4;
x2=307.2; y2=12.8;
x3=460.8; y3=64.0;

move_to (cr, x0, y0);
curve_to (cr, x1, y1, x2, y2, x3, y3);
set_line_width (cr, 10.0);
stroke (cr);
restore(cr);

move_to(cr, 12.0, 12.0);
set_source_rgb (cr, 0, 0, 0);
show_text(cr,"Stroke");
write_to_png(c,"stroke.png");

max_stars_repo_path: Module 3/Chapter07/cairo-stroke.jl

max_issues_repo_path: Module 3/Chapter07/cairo-stroke.jl

max_forks_repo_path: Module 3/Chapter07/cairo-stroke.jl


text:
module types

export Body, Moon, Command_Module, EarthMoonSystem

type Body{T}
    mass::T
    velocity::Vector{T}
    radius::T
    position::Vector{T}
end

typealias Moon Body

type Command_Module{T}
    mass::T
    velocity::Vector{T}
    radius::T
    position::Vector{T}
    positionE::Vector{T}
    positionH::Vector{T}
    velocityE::Vector{T}
    velocityH::Vector{T}
end

type EarthMoonSystem
    time::Float64
    earth::Body
    moon::Moon
    command_module::Command_Module
end

end

max_stars_repo_name: rodrigolece/grupo-fisica-compu

max_issues_repo_name: rodrigolece/grupo-fisica-compu

max_forks_repo_name: rodrigolece/grupo-fisica-compu