text:
@defcomp n2ocycle begin
    e_globalN2Oemissions=Parameter(index=[time],unit="Mtonne/year")
    e_0globalN2Oemissions=Parameter(unit="Mtonne/year", default=11.22281348065454)
    c_N2Oconcentration=Variable(index=[time],unit="ppbv")
    pic_preindustconcN2O=Parameter(unit="ppbv", default=270.)
    exc_excessconcN2O=Variable(unit="ppbv")
    c0_N2Oconcbaseyr=Parameter(unit="ppbv", default=324.)
    re_remainN2O=Variable(index=[time],unit="Mtonne")
    re_remainN2Obase=Variable(unit="Mtonne")
    nte_natN2Oemissions=Variable(index=[time],unit="Mtonne/year")
    air_N2Ofractioninatm=Parameter(unit="%", default=100.)
    tea_N2Oemissionstoatm=Variable(index=[time],unit="Mtonne/year")
    teay_N2Oemissionstoatm=Variable(index=[time],unit="Mtonne/t")
    y_year=Parameter(index=[time],unit="year")
    y_year_0=Parameter(unit="year")
    res_N2Oatmlifetime=Parameter(unit="year", default=114.)
    den_N2Odensity=Parameter(unit="Mtonne/ppbv", default=7.8)
    stim_N2Oemissionfeedback=Parameter(unit="Mtonne/degreeC", default=0.)
    rtl_g0_baselandtemp=Parameter(unit="degreeC", default=1.1683981941310047)
    rtl_g_landtemperature=Parameter(index=[time],unit="degreeC")

    function run_timestep(p, v, d, t)
        #note that Hope (2009) states that Equations 1-12 for methane also apply to N2O

        if is_first(t)
            #eq.3 from Hope (2006) - natural emissions feedback, using global temperatures calculated in ClimateTemperature component
            nte_0=p.stim_N2Oemissionfeedback*p.rtl_g0_baselandtemp
            v.nte_natN2Oemissions[t]=p.stim_N2Oemissionfeedback*p.rtl_g0_baselandtemp
            #eq.6 from Hope (2006) - emissions to atmosphere depend on the sum of natural and anthropogenic emissions
            v.tea_N2Oemissionstoatm[t]=(p.e_globalN2Oemissions[t]+v.nte_natN2Oemissions[t])*p.air_N2Ofractioninatm/100
            tea_0=(p.e_0globalN2Oemissions+nte_0)*p.air_N2Ofractioninatm/100
            v.teay_N2Oemissionstoatm[t]=(v.tea_N2Oemissionstoatm[t]+tea_0)*(p.y_year[t]-p.y_year_0)/2
            #adapted from eq.1 in Hope(2006) - calculate excess concentration in base year
            v.exc_excessconcN2O=p.c0_N2Oconcbaseyr-p.pic_preindustconcN2O
            #Eq. 2 from Hope (2006) - base-year remaining emissions
            v.re_remainN2Obase=v.exc_excessconcN2O*p.den_N2Odensity
            v.re_remainN2O[t]=v.re_remainN2Obase*exp(-(p.y_year[t]-p.y_year_0)/p.res_N2Oatmlifetime)+
                v.teay_N2Oemissionstoatm[t]*p.res_N2Oatmlifetime*(1-exp(-(p.y_year[t]-p.y_year_0)/p.res_N2Oatmlifetime))/(p.y_year[t]-p.y_year_0)
        else
            #eq.3 from Hope (2006) - natural emissions (carbon cycle) feedback, using global temperatures calculated in ClimateTemperature component
            #Here assume still using area-weighted average regional temperatures (i.e. land temperatures) for natural emissions feedback
            v.nte_natN2Oemissions[t]=p.stim_N2Oemissionfeedback*p.rtl_g_landtemperature[t-1]
            #eq.6 from Hope (2006) - emissions to atmosphere depend on the sum of natural and anthropogenic emissions
            v.tea_N2Oemissionstoatm[t]=(p.e_globalN2Oemissions[t]+v.nte_natN2Oemissions[t])*p.air_N2Ofractioninatm/100
            #eq.7 from Hope (2006) - average emissions to atm over time period
            v.teay_N2Oemissionstoatm[t]=(v.tea_N2Oemissionstoatm[t]+v.tea_N2Oemissionstoatm[t-1])*(p.y_year[t]-p.y_year[t-1])/2
            #eq.10 from Hope (2006) - remaining emissions in atmosphere
            v.re_remainN2O[t]=v.re_remainN2O[t-1]*exp(-(p.y_year[t]-p.y_year[t-1])/p.res_N2Oatmlifetime)+
                v.teay_N2Oemissionstoatm[t]*p.res_N2Oatmlifetime*(1-exp(-(p.y_year[t]-p.y_year[t-1])/p.res_N2Oatmlifetime))/(p.y_year[t]-p.y_year[t-1])
        end

        #eq.11 from Hope(2006) - N2O concentration
        v.c_N2Oconcentration[t]=p.pic_preindustconcN2O+v.exc_excessconcN2O*v.re_remainN2O[t]/v.re_remainN2Obase

    end
end

max_stars_repo_name: lrennels/PAGE-ICE

max_issues_repo_name: lrennels/PAGE-ICE

max_forks_repo_name: lrennels/PAGE-ICE


text:
#=
solves a simple random linear optimization problem (LP):
    min  c'x
    s.t. Ax = b, x >= 0
=#

using SparseArrays

struct LinearOptNative{T <: Real} <: ExampleInstanceNative{T}
    m::Int
    n::Int
    nz_frac::Float64
end

function build(inst::LinearOptNative{T}) where {T <: Real}
    (m, n, nz_frac) = (inst.m, inst.n, inst.nz_frac)
    @assert 0 < nz_frac <= 1

    # generate random data
    A = (nz_frac >= 1) ? rand(T, m, n) : sprand(T, m, n, nz_frac)
    A .*= 10
    b = vec(sum(A, dims = 2))
    c = rand(T, n)
    G = Diagonal(-one(T) * I, n) # TODO uniformscaling
    h = zeros(T, n)
    cones = Cones.Cone{T}[Cones.Nonnegative{T}(n)]

    model = Models.Model{T}(c, A, b, G, h, cones)
    return model
end

ext: jl


text:
abstract type AbstractDomain end

abstract type RealDomain <: AbstractDomain end

struct RealNumbers <: RealDomain end

const ℝ = RealNumbers()

Base.minimum(::RealNumbers) = static(-Inf)
Base.maximum(::RealNumbers) = static(Inf)

Base.in(x, ::RealNumbers) = isreal(x)

Base.show(io::IO, ::typeof(ℝ)) = print(io, "ℝ")

struct BoundedReals{L,U} <: RealDomain
    lower :: L
    upper :: U
end

Base.in(x, b::BoundedReals) = b.lower ≤ x ≤ b.upper

export ℝ, ℝ₊, 𝕀, ℤ

const ℝ₊ = BoundedReals(static(0.0), static(Inf))
const 𝕀 = BoundedReals(static(0.0), static(1.0))

Base.minimum(b::BoundedReals) = b.lower
Base.maximum(b::BoundedReals) = b.upper

Base.show(io::IO, ::typeof(ℝ₊)) = print(io, "ℝ₊")
Base.show(io::IO, ::typeof(𝕀)) = print(io, "𝕀")

testvalue(::typeof(ℝ)) = 0.0
testvalue(::typeof(ℝ₊)) = 1.0
testvalue(::typeof(𝕀)) = 0.5

abstract type IntegerDomain <: AbstractDomain end

struct IntegerNumbers <: IntegerDomain end

Base.in(x, ::IntegerNumbers) = isinteger(x)

const ℤ = IntegerNumbers()

Base.show(io::IO, ::typeof(ℤ)) = print(io, "ℤ")

Base.minimum(::IntegerNumbers) = static(-Inf)
Base.maximum(::IntegerNumbers) = static(Inf)
struct BoundedInts{L,U} <: IntegerDomain
    lower :: L
    upper :: U
end

Base.in(x, b::BoundedInts) = x ∈ ℤ && b.lower ≤ x ≤ b.upper

Base.minimum(b::BoundedInts) = b.lower
Base.maximum(b::BoundedInts) = b.upper

function Base.show(io::IO, b::BoundedInts)
    io = IOContext(io, :compact => true)
    print(io, "ℤ[", b.lower, ":", b.upper, "]")
end

testvalue(b::BoundedInts) = min(b.lower, 0)

function Base.getindex(::typeof(ℤ), r::AbstractUnitRange)
    BoundedInts(extrema(r)...)
end


text:
using EnvironmentMigrators
using Test
using Pkg

function collect_files(path::AbstractString)
    flat_files = String[]
    for (root, dirs, files) in walkdir(path)
        for f in files
            push!(flat_files, joinpath(root,f))
        end
    end
    flat_files
end

@testset "EnvironmentMigrators.jl" begin
    Pkg.activate(mktempdir())
    Pkg.add("REPL")
    @test EnvironmentMigrators.list_shared_environments() isa Vector{String}
    @test EnvironmentMigrators.backup_current_environment() === nothing
    current_env = Base.active_project()
    current_dir = dirname(current_env)
    backup_dir = joinpath(current_dir, "backups")
    @test isfile(current_env)
    @test isdir(backup_dir)
    @test backup_dir |> readdir |> !isempty
    timestamp_dir = backup_dir |> readdir |> first |> x->joinpath(backup_dir, x)
    @test timestamp_dir |> readdir |> !isempty
    timetstamp_dir_contents = timestamp_dir |> readdir |> x->filter(s->endswith(s, "toml"), x)
    @test "Project.toml" in timetstamp_dir_contents
    @test "Manifest.toml" in timetstamp_dir_contents
    sleep(2)
    @test EnvironmentMigrators.migrate_selected_environment(Base.active_project()) === nothing
    sleep(2)
    t = mktempdir()
    @test EnvironmentMigrators.migrate_selected_environment(t; backup = true) === nothing
    t = mktempdir()
    @test EnvironmentMigrators.migrate_selected_environment(t; backup = false) === nothing
    t = mktempdir()
    touch(joinpath(t, "Project.toml"))
    touch(joinpath(t, "Manifest.toml"))
    @test EnvironmentMigrators.migrate_selected_environment(t; backup = true) === nothing
    t = mktempdir()
    touch(joinpath(t, "JuliaProject.toml"))
    touch(joinpath(t, "JuliaManifest.toml"))
    @test EnvironmentMigrators.migrate_selected_environment(t; backup = true) === nothing
    t = mktempdir()
    touch(joinpath(t, "JuliaProject.toml"))
    touch(joinpath(t, "JuliaManifest.toml"))
    Pkg.activate(t)
    @test EnvironmentMigrators.migrate_selected_environment(t; backup = true) === nothing
    backup_files = collect_files(joinpath(t, "backups"))
    @test backup_files |> x->filter(s->endswith(s, "JuliaProject.toml"), x) |> !isempty
    @test backup_files |> x->filter(s->endswith(s, "JuliaManifest.toml"), x) |> !isempty
end


text:
function joiner_lineofsight(jsonlinex, jsonliney)
  x = JSON.parse(jsonlinex)
  y = JSON.parse(jsonliney)
  key = ""
  value = ""
  if y[1] == x["flightid"]
    key = x["flightid"] * "_" * x["originairport"] * "_" * x["dstairport"]
    value = x["passengerid"]
  end
  ret = Dict(key => value)
  return ret
end

max_stars_repo_name: deyandyankov/nmr

max_issues_repo_name: deyandyankov/nmr

max_forks_repo_name: deyandyankov/nmr


text:
"""
    VariableLinearAdvectionNonperiodicSemidiscretisation

A semidiscretisation of the linear advection equation
    \$\\partial_t u(t,x) + \\partial_x ( a(x) u(t,x) ) = 0\$
with boundary conditions `left_bc(t)`, `right_bc(t)`.
"""
struct VariableLinearAdvectionNonperiodicSemidiscretisation{T,Derivative<:AbstractDerivativeOperator{T},
                                                            Dissipation,
                                                            SplitForm<:Union{Val{false}, Val{true}},
                                                            LeftBC, RightBC} <: AbstractSemidiscretisation
    derivative::Derivative
    dissipation::Dissipation
    a::Vector{T}
    tmp1::Vector{T}
    tmp2::Vector{T}
    split_form::SplitForm
    left_bc::LeftBC
    right_bc::RightBC

    function VariableLinearAdvectionNonperiodicSemidiscretisation(derivative::Derivative, dissipation::Dissipation, afunc, split_form::SplitForm, left_bc::LeftBC, right_bc::RightBC) where {T, Derivative<:AbstractDerivativeOperator{T}, Dissipation, SplitForm<:Union{Val{false}, Val{true}}, LeftBC, RightBC}
        if dissipation != nothing
            @argcheck size(derivative) == size(dissipation) DimensionMismatch
            @argcheck grid(derivative) == grid(dissipation) ArgumentError
        end
        N = size(derivative, 2)
        a = compute_coefficients(afunc, derivative)
        tmp1 = Array{T}(undef, N)
        tmp2 = Array{T}(undef, N)
        new{T,Derivative,Dissipation,SplitForm,LeftBC,RightBC}(derivative, dissipation, a, tmp1, tmp2, split_form, left_bc, right_bc)
    end
end

function Base.show(io::IO, semidisc::VariableLinearAdvectionNonperiodicSemidiscretisation)
    print(io, "Semidiscretisation of the linear advection equation\n")
    print(io, "  \$ \\partial_t u(t,x) + \\partial_x ( a(x) u(t,x) ) = 0 \$ \n")
    print(io, "with nonperiodic boundaries using")
    if semidisc.split_form == Val{true}()
        print(io, " a split form and: \n")
    else
        print(io, " no split form and: \n")
    end
    print(io, semidisc.derivative)
    print(io, semidisc.dissipation)
end

function godunov_flux_variablelinearadvection(uₗ::T, uᵣ::T, a::T) where {T<:Real}
    ifelse(a > 0, a*uₗ, a*uᵣ)
end

function (disc::VariableLinearAdvectionNonperiodicSemidiscretisation)(du, u, p, t)
    @unpack a, tmp1, tmp2, derivative, dissipation, split_form, left_bc, right_bc = disc
    @boundscheck begin
        @argcheck length(u) == length(a)
        @argcheck length(u) == length(tmp1)
        @argcheck length(u) == length(tmp2)
        @argcheck length(u) == length(du)
    end

    # volume terms
    if typeof(split_form) <: Val{true}
        m1_2 = -one(eltype(u)) / 2

        ## a * D * u
        mul!(tmp1, derivative, u)
        @. du = m1_2 * a * tmp1
        ## u * D * a
        mul!(tmp1, derivative, a)
        @. du += m1_2 * u * tmp1
        ## D * a*u
        @. tmp2 = a*u
        mul!(tmp1, derivative, tmp2)
        @. du += m1_2 * tmp1
    else
        @. tmp2 = a * u
        mul!(tmp1, derivative, tmp2)
        @. du = -tmp1
    end

    # dissipation
    if dissipation != nothing
        mul!(tmp1, dissipation, u)
        @. du += tmp1
    end

    # boundary conditions via Godunov's flux
    @inbounds fnum_left = godunov_flux_variablelinearadvection(left_bc(t), u[1], a[1])
    @inbounds du[1] += (fnum_left - a[1]*u[1]) / left_boundary_weight(derivative)
    @inbounds fnum_right = godunov_flux_variablelinearadvection(u[end], right_bc(t), a[end])
    @inbounds du[end] -= (fnum_right - a[end]*u[end]) / right_boundary_weight(derivative)

    nothing
end


text:
@testset "232.implement-queue-using-stacks.jl" begin
    q = MyQueue()
    @test myqueue_push!(q, 1) == [1]
    @test myqueue_push!(q, 2) == [1, 2]
    @test myqueue_peek(q) == 1
    @test myqueue_pop!(q) == 1
    @test myqueue_isempty(q) == false
end


text:
using MVDA, Plots, Random, StableRNGs, LaTeXStrings

# Note: only show samples used in fitting and selecting models

w, h = default(:size)
upscale = 1.0
font_main = Plots.font("Computer Modern", 24)
font_large = Plots.font("Computer Modern", 16)
font_small = Plots.font("Computer Modern", 12)
default()
default(
    markerstrokewidth=0.5,
    titlefont=font_main,
    guidefont=font_large,
    tickfont=font_small,
    legendfont=font_small,
    size=(1.5*w, 1.5*h),
    markersize=5,
    linewidth=3,
    dpi=200,
    thickness_scaling=1.0,
    left_margin=5Plots.mm,
    bottom_margin=5Plots.mm,
)
scalefontsizes(upscale)

make_figure = function(Y, X, is2D::Bool=true)
    if is2D
        scatter(X[:,1], X[:,2],
            group=Y,
            aspect_ratio=1,
            xlabel=L"feature $x_1$",
            ylabel=L"feature $x_2$",
            markersize=4,
            markerstrokewidth=1e-4,
        )
    else
        scatter(X[:,1], X[:,2], X[:,3],
            group=Y,
            aspect_ratio=1,
            xlabel=L"feature $x_1$",
            ylabel=L"feature $x_2$",
            zlabel=L"feature $x_3$",
            markersize=4,
            markerstrokewidth=1e-4,
        )
    end
end

draw_cloud_boundaries! = function() # assumes 3 classes
    xs = range(-2,2,length=11)
    plot!(xs, x -> sqrt(2)/2*x, color=:black, style=:dash, label=nothing)
    plot!(xs, x -> -sqrt(2)/2*x, color=:black, style=:dash, label=nothing)
    vline!([0.0], color=:black, style=:dash, label=nothing)
end

draw_circle! = function(r)
    xs = range(-r, r, length=100)
    ys = sqrt.(r^2 .- xs .^2)
    plot!(xs, ys, color=:black, style=:dash, label="", lw=2)
    plot!(xs, -ys, color=:black, style=:dash, label="", lw=2)
end

length(ARGS) < 1 && error("Script requires an argument to specify a target directory.")
dir = ARGS[1]

# VDA Simulation
n_cv, n_test = 500, 500
nsamples = n_cv + n_test
p = 50
c = 10
samples_per_class = nsamples ÷ c
rng = StableRNG(1234)

Y, X, _ = MVDA.simulate_ground_truth(p, c, samples_per_class, d=1.0, rng=rng, sigma=1.0)
fig_filename = joinpath(dir, "10clouds-1.png")
make_figure(Y[1:n_cv], X[1:n_cv,:], false)
title!(L"10clouds / $n=%$n_cv$ / $p=%$p$ / $c=%$c$ / $d=1$")
@info "Saving $(fig_filename)"
savefig(fig_filename)

Y, X, _ = MVDA.simulate_ground_truth(p, c, samples_per_class, d=1.0, rng=rng, sigma=5.0)
fig_filename = joinpath(dir, "10clouds-2.png")
make_figure(Y[1:n_cv], X[1:n_cv,:], false)
title!(L"10clouds / $n=%$n_cv$ / $p=%$p$ / $c=%$c$ / $d=1$")
@info "Saving $(fig_filename)"
savefig(fig_filename)

Y, X, _ = MVDA.simulate_ground_truth(p, c, samples_per_class, d=1.0, rng=rng, sigma=10.0)
fig_filename = joinpath(dir, "10clouds-3.png")
make_figure(Y[1:n_cv], X[1:n_cv,:], false)
title!(L"10clouds / $n=%$n_cv$ / $p=%$p$ / $c=%$c$ / $d=1$")
@info "Saving $(fig_filename)"
savefig(fig_filename)

Y, X, _ = MVDA.simulate_ground_truth(p, c, samples_per_class, d=3.0, rng=rng, sigma=1.0)
fig_filename = joinpath(dir, "10clouds-4.png")
make_figure(Y[1:n_cv], X[1:n_cv,:], false)
title!(L"10clouds / $n=%$n_cv$ / $p=%$p$ / $c=%$c$ / $d=3$")
@info "Saving $(fig_filename)"
savefig(fig_filename)

Y, X, _ = MVDA.simulate_ground_truth(p, c, samples_per_class, d=3.0, rng=rng, sigma=5.0)
fig_filename = joinpath(dir, "10clouds-5.png")
make_figure(Y[1:n_cv], X[1:n_cv,:], false)
title!(L"10clouds / $n=%$n_cv$ / $p=%$p$ / $c=%$c$ / $d=3$")
@info "Saving $(fig_filename)"
savefig(fig_filename)

Y, X, _ = MVDA.simulate_ground_truth(p, c, samples_per_class, d=3.0, rng=rng, sigma=10.0)
fig_filename = joinpath(dir, "10clouds-6.png")
make_figure(Y[1:n_cv], X[1:n_cv,:], false)
title!(L"10clouds / $n=%$n_cv$ / $p=%$p$ / $c=%$c$ / $d=3$")
@info "Saving $(fig_filename)"
savefig(fig_filename)

# Gaussian Clouds
n_cv, n_test = 250, 10^3
nsamples = n_cv + n_test
nclasses = 3
Y, X = MVDA.simulate_gaussian_clouds(nsamples, nclasses; sigma=0.25, rng=StableRNG(1903))
fig_filename = joinpath(dir, "clouds.png")
make_figure(Y[1:n_cv], X[1:n_cv,:])
draw_cloud_boundaries!()
title!(L"clouds / $n=%$n_cv$ / $p=2$ / $c=%$nclasses$")
@info "Saving $(fig_filename)"
savefig(fig_filename)

# Nested Circles
n_cv, n_test = 250, 10^3
nsamples = n_cv + n_test
nclasses = 3
Y, X = MVDA.simulate_nested_circles(nsamples, nclasses; p=8//10, rng=StableRNG(1903))
fig_filename = joinpath(dir, "circles.png")
make_figure(Y[1:n_cv], X[1:n_cv,:])
title!(L"clouds / $n=%$n_cv$ / $p=2$ / $c=%$nclasses$")
foreach(r -> draw_circle!(sqrt(r)), 1:nclasses)
@info "Saving $(fig_filename)"
savefig(fig_filename)

# Waveform
n_cv, n_test = 375, 10^3
nsamples = n_cv + n_test
nfeatures = 21
Y, X = MVDA.simulate_waveform(nsamples, nfeatures; rng=StableRNG(1903))
fig_filename = joinpath(dir, "waveform.png")
Y, X = Y[1:n_cv], X[1:n_cv,:]
fig = plot(layout=@layout[a{1e-6h}; grid(3,1)])
title!(L"waveform / $n=%$n_cv$ / $p=21$ / $c=3$", subplot=1)
plot!(framestyle=:none, xticks=nothing, yticks=nothing, subplot=1)
@views for c in 1:3
    idx = findall(isequal(c), Y)
    foreach(i -> plot!(axes(X, 2), X[i,:], color=c, label=nothing, title="Class $c", alpha=0.25, subplot=c+1), idx)
end
ylabel!(L"feature $x_{j}$", subplot=3)
xlabel!(L"index $j$", subplot=4)
xticks!(1:nfeatures)
@info "Saving $(fig_filename)"
savefig(fig_filename)


text:
global tested_problem_types = [:multiobjective]

@testset "CoverEnumerator problem_type=$problem_type" for problem_type in tested_problem_types # FIXME use weights
    @testset "[a b c abc], mask=[a b]" begin
        sm = SetMosaic([Set([:a]), Set([:b]), Set([:c]), Set([:a, :b, :c])])
        sm_ab = mask(sm, [Set([:a, :b])], min_nmasked=1)

        # low penality to select sets, high probability to miss active element, so select a and b, and abc
        cover_params1 = CoverParams(sel_tax=-log(0.5), uncovered_factor=0.5)
        cover_problem = MultiobjCoverProblem(sm_ab, cover_params1)
        @test nvars(cover_problem) == 3
        @test aggscore([1.0, 1.0, 0.0], cover_problem) < aggscore([0.0, 0.0, 0.0], cover_problem)
        @test aggscore([0.0, 0.0, 1.0], cover_problem) < aggscore([0.0, 0.0, 0.0], cover_problem)
        @test aggscore([1.0, 1.0, 0.0], cover_problem) < aggscore([1.0, 1.0, 1.0], cover_problem)
        @test aggscore([0.0, 0.0, 1.0], cover_problem) < aggscore([1.0, 1.0, 1.0], cover_problem)
        cover_coll1a = collect(sm_ab, cover_params1, CoverEnumerationParams(max_set_score=10.0), problem_type=problem_type)
        @test length(cover_coll1a) == 2
        cover_coll1b = collect(sm_ab, cover_params1, CoverEnumerationParams(max_set_score=0.4), problem_type=problem_type)
        @test length(cover_coll1b) == 1
        cover_coll1c = collect(sm_ab, cover_params1, CoverEnumerationParams(max_set_score=-1.0), problem_type=problem_type)
        @test length(cover_coll1c) == 0

        # higher penalty to select sets, high probability to miss active element, so select abc
        cover_params2 = CoverParams(sel_tax=-log(0.75))
        cover_coll2a = collect(sm_ab, cover_params2, CoverEnumerationParams(max_set_score=10.0), problem_type=problem_type)
        if problem_type==:quadratic # FIXME workfor all
            @test length(cover_coll2a) == 1
        end
        cover_coll2b = collect(sm_ab, cover_params2, CoverEnumerationParams(max_set_score=-1.0), problem_type=problem_type)
        @test length(cover_coll2b) == 0

        # higher prior probability to select sets, lower probability to miss active element, so select a and b, then abc
        cover_params3 = CoverParams(sel_tax=0.0)
        cover_coll3 = collect(sm_ab, cover_params3, CoverEnumerationParams(max_set_score=10.0), problem_type=problem_type)
        @test length(cover_coll3) == 2 # FIXME use a and b weights
    end

    @testset "[a b c abc de], mask=[a b]" begin
        sm = SetMosaic([Set([:a]), Set([:b]), Set([:c]), Set([:a, :b, :c])],
                       Set([:a, :b, :c, :d, :e]))

        # higher prior probability to select sets, lower probability to miss active element, so select a and b, then abc
        cover_coll = collect(mask(sm, [Set([:a, :b])], min_nmasked=1), CoverParams(sel_tax=0.0),
                             CoverEnumerationParams(max_set_score=10.0), problem_type=problem_type, ϵ=0.01)
        @test length(cover_coll) == 2
        @test best_aggscore(cover_coll.results[1]) <= best_aggscore(cover_coll.results[2])
    end

    @testset "DataFrame(CoverCollection)" begin
        sm = SetMosaic([Set([:a]), Set([:b]), Set([:c]), Set([:a, :b, :c])],
                       Set([:a, :b, :c, :d, :e]))
        sm_ab = mask(sm, Dict(:X => Set([:a, :b])), min_nmasked=1)

        # higher prior probability to select sets, lower probability to miss active element, so select a and b, then a b c
        cover_coll = collect(sm_ab, CoverParams(sel_tax=0.0),
                             CoverEnumerationParams(max_set_score=10.0),
                             problem_type=problem_type, ϵ=0.01)

        df = DataFrame(cover_coll)
        @test size(df, 1) == 3
        @test df.cover_ix == [1, 1, 2]
        @test df.set_id == [1, 2, 4]
        @test df.experiment_ix == [1, 1, 1]
        @test df.experiment_id == [:X, :X, :X]
        @test df == DataFrame(cover_coll, min_nmasked=1)

        df2 = DataFrame(cover_coll, min_nmasked=0)
        @test size(df2, 1) == 3
        @test df2.cover_ix == [1, 1, 2]
        @test df2.set_id == [1, 2, 4]
    end

    @testset "multimask: [abd bcd c d abcde cde], mask=[abc be]" begin # FIXME take weights into account
        sm = SetMosaic(Dict("abd" => Set([:a, :b, :d]), "bcd" => Set([:b, :c, :d]),
                            "c" => Set([:c]), "d" => Set([:d]),
                            "abcde" => Set([:a, :b, :c, :d, :e]),
                            "cdef" => Set([:c, :d, :e, :f])),
                        Set([:a, :b, :c, :d, :e, :f]))
        sm_abc_be = mask(sm, Dict(:X => Set([:a, :b, :c]), :Y => Set([:b, :e])), min_nmasked=1)

        # higher prior probability to select sets, high overlap penalty,
        # high penalty for covering unmasked, so select [abd c], [bcd] and [abcde]
        cover_coll = collect(sm_abc_be, CoverParams(setXset_factor=1.0, covered_factor=0.05, uncovered_factor=0.1, sel_tax=-log(0.9)),
                             CoverEnumerationParams(max_set_score=0.0),
                             problem_type==:quadratic ?
                                OESC.QuadraticOptimizerParams() :
                                OESC.MultiobjOptimizerParams(ϵ=0.01),
                             false)

        df = DataFrame(cover_coll, best_only=true)
        @test size(df, 1) == 8
        @test df.cover_ix == [1, 1, 1, 1, 2, 2, 3, 3]
        @test df.experiment_ix == [1, 2, 1, 2, 1, 2, 1, 2]
        @test df.experiment_id == [:X, :Y, :X, :Y, :X, :Y, :X, :Y]
        @test df.set_id == ["c", "c", "abd", "abd", "abcde", "abcde", "bcd", "bcd"]
    end

    @testset "multiexperiment weighteed: [abd bcd c d abcde cde]" begin # FIXME take weights into account
        sm = SetMosaic(Dict("abd" => Set([:a, :b, :d]), "bcd" => Set([:b, :c, :d]),
                            "c" => Set([:c]), "d" => Set([:d]),
                            "abcde" => Set([:a, :b, :c, :d, :e]),
                            "cdef" => Set([:c, :d, :e, :f])),
                        Set([:a, :b, :c, :d, :e, :f]))
        sm_abc_be = assignweights(sm, Dict(:X => Dict("abd" => -5.0, "bcd" => -5.0,
                                                      "c" => -2.0, "abcd" => -6.0,
                                                      "cdef" => -3.0),
                                           :Y => Dict("abd" => -1.0, "bcd" => -1.0,
                                                      "abcde" => -4.0, "cdef" => -2.0)))

        # higher prior probability to select sets, high overlap penalty,
        # high penalty for covering unmasked, so select [abd c], [bcd] and [abcde]
        cover_coll = collect(sm_abc_be, CoverParams(setXset_factor=1.0, covered_factor=0.05, uncovered_factor=0.1, sel_tax=-log(0.9)),
                             CoverEnumerationParams(max_set_score=0.0),
                             problem_type==:quadratic ?
                                OESC.QuadraticOptimizerParams() :
                                OESC.MultiobjOptimizerParams(ϵ=0.01),
                             false)

        df = DataFrame(cover_coll, best_only=true)
        @test size(df, 1) == 8
        @test df.cover_ix == [1, 1, 1, 1, 2, 2, 3, 3]
        @test df.experiment_ix == [1, 2, 1, 2, 1, 2, 1, 2]
        @test df.experiment_id == [:X, :Y, :X, :Y, :X, :Y, :X, :Y]
        @test df.set_id == ["c", "c", "abd", "abd", "abcde", "abcde", "bcd", "bcd"]
    end
end


text:
using SharedArrays
using Interpolations: BSplineInterpolation, tweight, scale
using Random
# using BSplineExtensions
# using Interpolations
# using Test

@testset "bspline funs" begin

    x = SharedArray{Float64}(100, 100, 10, 10; pids=[1,])
    endsz = (map(i -> size(x,i), 3:4)...,)
    rand!(x)

    x2210 = copy(x)
    x2211 = copy(x)
    x1111 = copy(x)

    b1 = BSpline(Linear())
    b2 = BSpline(Quadratic(InPlace()))
    b0 = NoInterp()

    it2210 = (b2, b2, b1, b0)
    it2211 = (b2, b2, b1, b1)
    it1111 = (b1, b1, b1, b1)

    # interpolate a bunch of stuff in place
    itp_2210 = interpolate!(x2210, it2210, OnCell())
    itp_2211 = interpolate!(x2211, it2211, OnCell())
    itp_1111 = interpolate!(x1111, it1111, OnCell())

    @test itp_2210.coefs === x2210
    @test itp_2211.coefs === x2211
    @test itp_1111.coefs === x1111

    @test !all(x2210 .== x)
    @test !all(x2211 .== x)
    @test  all(x1111 .== x)

    # store these coefs
    coefs2210 = copy(x2210)
    coefs2211 = copy(x2211)
    coefs1111 = copy(x1111)

    # restore original x
    x2210 .= x
    x2211 .= x
    x1111 .= x

    # interpolate w/ homemade algorithm and see what we get!
    myitp_2210 = BSplineInterpolation(tweight(x2210), x2210, it2210, OnCell(), Val{0}())
    myitp_2211 = BSplineInterpolation(tweight(x2211), x2211, it2211, OnCell(), Val{0}())
    myitp_1111 = BSplineInterpolation(tweight(x1111), x1111, it1111, OnCell(), Val{0}())

    @test typeof(myitp_2210) == typeof(itp_2210)
    @test typeof(myitp_2211) == typeof(itp_2211)
    @test typeof(myitp_1111) == typeof(itp_1111)

    @test myitp_2210.coefs === x2210
    @test myitp_2211.coefs === x2211
    @test myitp_1111.coefs === x1111

    serial_prefilterByView!(x2210, myitp_2210, endsz...)
    serial_prefilterByView!(x2211, myitp_2211, endsz...)
    serial_prefilterByView!(x1111, myitp_1111, endsz...)

    @test all(myitp_2210.coefs .== coefs2210)
    @test all(myitp_2211.coefs .== coefs2211)
    @test all(myitp_1111.coefs .== coefs1111)

    gradient_d_impl(3, typeof(itp_2210))
    # gradient_impl(typeof(itp_2210))

    function test_grad_d(itp::BSplineInterpolation, x::Number...)
        g = Interpolations.gradient(itp, x...)
        for d = 1:length(g)
            @test g[d] == gradient_d(Val{d}, itp, x...)
        end
    end

    xs = 1.2, 2.2, 3.2, Int(4)
    # @code_warntype gradient_d(Val{1}, itp_2210, xs...)
    gradient!(Vector{Float64}(undef, 3), itp_2210, xs...)

    Interpolations.gradient(itp_2210, xs...)
    test_grad_d(itp_2210, xs...)
    test_grad_d(itp_1111, xs...)
    test_grad_d(itp_2211, xs...)

    @test_throws DomainError    gradient_d(Val{5}, itp_1111, 1:4...)
    @test_throws ErrorException gradient_d(Val{4}, itp_2210, 1:4...)

    # ------------------------ scale & gradient ---------------------

    rng1 = range(-30, stop=30, length=100)
    rng2 = range(20, stop=50, length=100)
    rng3 = range(30, stop=39, length=10)
    rng4 = 1:10

    sitp_2210 = scale(itp_2210, rng1, rng2, rng3, rng4)

    Random.seed!(1234)
    z0 = rand(10,10)
    za = copy(z0)
    zb = copy(z0')
    @test all(za .== zb')

    itpa = interpolate(za, (BSpline(Linear()), NoInterp()), OnGrid())
    itpb = interpolate(zb, (NoInterp(), BSpline(Linear())), OnGrid())
    @test all(itpa.coefs .== itpb.coefs')

    rng = range(1.0, stop=19.0, length=10)
    sitpa = scale(itpa, rng, 1:10)
    sitpb = scale(itpb, 1:10, rng)
    @test sitpa[2.0, 3] == sitpb[3,2.0]
    @test sitpa[2.0, 3] != sitpb[3,2.1]
    @test itpa[2.0, 3] == sitpa[3.0, 3]

    @test Interpolations.gradient(itpa, 2.0, 3)/step(rng) == Interpolations.gradient(sitpa, 3.0, 3)
    @test Interpolations.gradient(itpa, 2.0, 3)/step(rng) == Interpolations.gradient(itpb, 3, 2.0)/step(rng)
    # @test Interpolations.gradient(itpb, 3, 2.0)/step(rng) == Interpolations.gradient(sitpb, 3, 3.0)

    # @test Interpolations.gradient(sitpb, 3, 3.0) == Interpolations.gradient(sitpa, 3.0, 3)
    @test ShaleDrillingModel.fixedgradient(sitpa, 3.0, 3) == ShaleDrillingModel.fixedgradient(sitpb, 3, 3.0)
    @test gradient_d(Val{1}, sitpa, 3.0, 3) ==  ShaleDrillingModel.fixedgradient(sitpa, 3.0, 3)[1]

    println("done!")

    # Interpolations.gradient(sitpb, 3, 2.0)

    # @show mygradient(sitpa, 2.0, 3) , mygradient(sitpb, 3, 2.0)
    # @test Interpolations.gradient(sitpa, 2.0, 3) ==  Interpolations.gradient(sitpb, 3, 2.0)

    # @test Interpolations.gradient(sitpa, 2.0, 3) ==  Interpolations.gradient(sitpb, 3, 2.0)

    # it2201 = (b2, b2, b0, b1)
    # x2201 = copy(x)
    # itp_2201 = interpolate!(x2201, it2201, OnCell())
    # sitp_2201 = scale(itp_2201, rng1, rng2, rng4, rng3)
    #
    # function f(sitp::ScaledInterpolation{T,N,ITPT,IT,GT}) where {T,N,ITPT,IT<:DimSpec{InterpolationType},GT<:DimSpec{GridType}}
    #     @show IT.parameters
    #     @show length(IT.parameters)
    #     @show interp_types = length(IT.parameters) == N ? IT.parameters : tuple([IT.parameters[1] for _ in 1:N]...)
    #     @show interp_dimens = map(it -> interp_types[it] != NoInterp, 1:N)
    #     @show interp_indices = map(i -> interp_dimens[i] ? :(coordlookup(sitp.ranges[$i], xs[$i])) : :(xs[$i]), 1:N)
    # end
    #
    # mygradient(sitp_2201, 1.2, 2.2, Int(4), 3.2)
    # mygradient(sitp_2210, 1.2, 2.2, 3.2, Int(4))

    # f(sitp_2201)
end


text:
using Gtk
using GtkBuilderAid
using LibCURL

export start_intro

@GtkBuilderAid function_name(build_intro) begin

@guarded function language_toggled(cell, path_ptr, built)
  store = GAccessor.object(built, "programming_languages")
  path_str = pointer_to_string(path_ptr)
  iter = Gtk.iter_from_string_index(store, path_str)
  store[iter, 2] = Cint(!store[iter, 2])
  nothing
end

@guarded function submit(widget, built)
  store = GAccessor.object(built, "programming_languages")
  known_languages = []
  color_iter = Ref(GtkTreeIter())
  color = if ccall(
      (:gtk_combo_box_get_active_iter, Gtk.libgtk),
      Cint,
      (Ptr{GObject}, Ptr{GtkTreeIter}),
      GAccessor.object(built, "color_selection"),
      color_iter) != 0
    color_store = GAccessor.object(built, "colors")
    GtkTreeModel(color_store)[color_iter[], 1]
  else
    "DislikesColor"
  end

  model = GtkTreeModel(store)
  iter = Gtk.mutable(GtkTreeIter)
  if Gtk.get_iter_first(model, iter)
    last = true
    while last
      if store[iter, 2]
        push!(known_languages, store[iter, 1])
      end
      last = next(model, iter)
    end
  end

  curl = curl_easy_init()

  curl_easy_setopt(curl, CURLOPT_URL, string(
    "http://hackrva.org/cgi-bin/collect.py?languages=",
    join(known_languages, ","),
    "&color=",
    color))
  curl_easy_setopt(curl, CURLOPT_FOLLOWLOCATION, 1)

#   c_curl_read_cb = cfunction(curl_read_cb, Csize_t, (Ptr{Void}, Csize_t, Csize_t, Ptr{Void}))
#   curl_easy_setopt(curl, CURLOPT_READFUNCTION, c_curl_read_cb)

  println("Chosen color:")
  println(color)
  println("I know:")
  println(known_languages)
  destroy(GAccessor.object(built, "main_window"))
  # Blocks execution
  curl_easy_perform(curl)
end

end

@guarded function activate_cb(app_ptr, udata)

  app = GObject(app_ptr)
  built = @GtkBuilder(filename=joinpath(Pkg.dir("HackRvaJulia"), "resources", "introduction.glade"))

  build_intro(built, built)

  win = GAccessor.object(built, "main_window")
  push!(app, win)
  showall(win)

  nothing
end

"""
Starts a GTK+ GUI for introducing people to the power of julia. This
demonstrates that even complicated GUI applications are possible rather
easily. It doesn't even take that many lines of code either.
"""
function start_intro()
  println("Preventing a GTK related Heisenbug by Printing!")
  app = @GtkApplication("com.matt5sean3.juliademo.introduction", 0)
  signal_connect(activate_cb, app, "activate", Void, (), false, ())
  run(app)
end


text:
using LightGraphsFlows
using Test
import LightGraphs
import SimpleTraits
const lg = LightGraphs

using LinearAlgebra: diag

const testdir = dirname(@__FILE__)

testgraphs(g) = [g, lg.Graph{UInt8}(g), lg.Graph{Int16}(g)]
testdigraphs(g) = [g, lg.DiGraph{UInt8}(g), lg.DiGraph{Int16}(g)]

for t in [
        "edmonds_karp",
        "dinic",
        "boykov_kolmogorov",
        "push_relabel",
        "maximum_flow",
        "multiroute_flow",
        "mincost",
        "mincost_dep",  #Please remove this once depreciated mincost_flow iterface is removed
        "mincut",
        ]
    tp = joinpath(testdir, "$(t).jl")
    include(tp)
end


text:
using ContinuedFractions

cf = ContinuedFraction(sqrt(2))
@assert quotients(cf)[1] == 1
for i in 2:15
	@assert quotients(cf)[i] == 2
end
@assert abs(collect(convergents(cf))[end] - sqrt(2)) < 10e-16

cf = ContinuedFraction(sqrt(3))
@assert quotients(cf)[1] == 1
for i in 2:15
	if mod(i, 2) == 0
		@assert quotients(cf)[i] == 1
	else
		@assert quotients(cf)[i] == 2
	end
end
@assert abs(collect(convergents(cf))[end] - sqrt(3)) < 10e-16

cf = ContinuedFraction(1e)
@assert quotients(cf)[1] == 2
for i in 2:15
	if mod(i, 3) == 0
		@assert quotients(cf)[i] == fld(i, 3) * 2
	else
		@assert quotients(cf)[i] == 1
	end
end


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

using Base.Test

f = Float16(2.)
g = Float16(1.)

@test f >= g
@test f > g
@test g < f
@test g <= g
@test all([g g] .< [f f])
@test all([g g] .<= [f f])
@test all([f f] .> [g g])
@test all([f f] .>= [g g])
@test isless(g, f)
@test !isless(f, g)

@test convert(Bool,Float16(0.0)) == false
@test convert(Bool,Float16(1.0)) == true
@test_throws InexactError convert(Bool,Float16(0.1))

@test convert(Int128,Float16(-2.0)) == Int128(-2)
@test convert(UInt128,Float16(2.0)) == UInt128(2)

# convert(::Type{Int128},  x::Float16)
@test convert(Int128, Float16(1.0)) === Int128(1.0)
@test convert(Int128, Float16(-1.0)) === Int128(-1.0)
@test_throws InexactError convert(Int128, Float16(3.5))

# convert(::Type{UInt128}, x::Float16)
@test convert(UInt128, Float16(1.0)) === UInt128(1.0)
@test_throws InexactError convert(UInt128, Float16(3.5))
@test_throws InexactError convert(UInt128, Float16(-1))

@test Float16(0.5f0)^2 ≈ Float16(0.5f0^2)
@test round(Int,Float16(0.5f0)) == round(Int,0.5f0)
@test trunc(Int,Float16(0.9f0)) === trunc(Int,0.9f0) === 0
@test floor(Int,Float16(0.9f0)) === floor(Int,0.9f0) === 0
@test trunc(Int,Float16(1)) === 1
@test floor(Int,Float16(1)) === 1
@test ceil(Int,Float16(0.1f0)) === ceil(Int,0.1f0) === 1
@test ceil(Int,Float16(0)) === ceil(Int,0) === 0
@test round(Float16(0.1f0)) == round(0.1f0) == 0
@test round(Float16(0.9f0)) == round(0.9f0) == 1
@test trunc(Float16(0.9f0)) == trunc(0.9f0) == 0
@test floor(Float16(0.9f0)) == floor(0.9f0) == 0
@test trunc(Float16(1)) === Float16(1)
@test floor(Float16(1)) === Float16(1)
@test ceil(Float16(0.1)) == ceil(0.1)
@test ceil(Float16(0.9)) == ceil(0.9)
@test unsafe_trunc(UInt8, Float16(3)) === 0x03
@test unsafe_trunc(Int16, Float16(3)) === Int16(3)
@test unsafe_trunc(UInt128, Float16(3)) === UInt128(3)
@test unsafe_trunc(Int128, Float16(3)) === Int128(3)
@test unsafe_trunc(Int16, NaN16) === Int16(0)  #18771
@test fma(Float16(0.1),Float16(0.9),Float16(0.5)) ≈ fma(0.1,0.9,0.5)
@test muladd(Float16(0.1),Float16(0.9),Float16(0.5)) ≈ muladd(0.1,0.9,0.5)
@test convert(Int128,Float16(-1.0)) == Int128(-1)
@test convert(UInt128,Float16(5.0)) == UInt128(5)

@test -f === Float16(-2.)

@test f+g === Float16(3f0)
@test f-g === Float16(1f0)
@test f*g === Float16(2f0)
@test f/g === Float16(2f0)
@test f^g === Float16(2f0)
@test f^1 === Float16(2f0)
@test f^-g === Float16(0.5f0)

@test f + 2 === Float16(4f0)
@test f - 2 === Float16(0f0)
@test f*2 === Float16(4f0)
@test f/2 === Float16(1f0)
@test f + 2. === 4.
@test f - 2. === 0.
@test f*2. === 4.
@test f/2. === 1.

@test sin(f) ≈ sin(2f0)

@test isnan(NaN16)
@test isnan(-NaN16)
@test !isnan(Inf16)
@test !isnan(-Inf16)
@test !isnan(Float16(2.6))
@test NaN16 != NaN16
@test isequal(NaN16, NaN16)
@test repr(NaN16) == "NaN16"
@test sprint(showcompact, NaN16) == "NaN"

@test isinf(Inf16)
@test isinf(-Inf16)
@test !isinf(NaN16)
@test !isinf(-NaN16)
@test !isinf(Float16(2.6))
@test Inf16 == Inf16
@test Inf16 != -Inf16
@test -Inf16 < Inf16
@test isequal(Inf16, Inf16)
@test repr(Inf16) == "Inf16"
@test sprint(showcompact, Inf16) == "Inf"

@test repr(Float16(44099)) == "Float16(4.41e4)"

for z1 in (Float16(0.0), Float16(-0.0)), z2 in (Float16(0.0), Float16(-0.0))
    @test z1 == z2
    @test isequal(z1, z1)
    @test z1 === z1
    for elty in (Float32, Float64)
        z3 = convert(elty, z2)
        @test z1==z3
    end
end

@test Float16(2.5) == Float16(2.5)
@test Float16(2.5) != Float16(2.6)
@test isequal(Float16(0.0), Float16(0.0))
@test !isequal(Float16(-0.0), Float16(0.0))
@test !isequal(Float16(0.0), Float16(-0.0))

@test isnan(reinterpret(Float16,0x7c01))
@test !isinf(reinterpret(Float16,0x7c01))

@test nextfloat(Inf16) === Inf16
@test prevfloat(-Inf16) === -Inf16

# rounding in conversions
let f
    for f32 in [.3325f0, -.3325f0]
        f16 = Float16(f32)
        # need to round away from 0. make sure we picked closest number.
        @test abs(f32 - f16) < abs(f32 - nextfloat(f16))
        @test abs(f32 - f16) < abs(f32 - prevfloat(f16))
    end
    # halfway between and last bit is 1
    f = reinterpret(Float32,                           0b00111110101010100011000000000000)
    @test Float32(Float16(f)) === reinterpret(Float32, 0b00111110101010100100000000000000)
    # halfway between and last bit is 0
    f = reinterpret(Float32,                           0b00111110101010100001000000000000)
    @test Float32(Float16(f)) === reinterpret(Float32, 0b00111110101010100000000000000000)
end

# issue #5948
@test string(reinterpret(Float16, 0x7bff)) == "6.55e4"

@test log10(Float16(100)) == Float16(2.0)

#  #9939 (and #9897)
@test rationalize(Float16(0.1)) == 1//10

# issue #17148
@test rem(Float16(1.2), Float16(one(1.2))) == 0.20019531f0

# no domain error is thrown for negative values
@test cbrt(Float16(-1.0)) == -1.0

max_stars_repo_name: visr/julia

max_issues_repo_name: visr/julia

max_forks_repo_name: visr/julia


text:
import Base.convert
using OpenMDAO
using LinearAlgebra: norm

struct CCBladeResidualComp
    num_nodes
    num_radial
    af
    B
    turbine
    debug_print
    inputs
    outputs
    partials
end

convert(::Type{CCBladeResidualComp}, po::PyCall.PyObject) = CCBladeResidualComp(
    po.num_nodes, po.num_radial, po.af, po.B, po.turbine, po.debug_print,
    po.inputs, po.outputs, po.partials)

function CCBladeResidualComp(; num_nodes, num_radial, af, B, turbine, debug_print)

    # Putting these two lines at the top level of this file doesn't work:
    # get_rows_cols will be a "Null" PyObject, and then the first call to
    # get_rows_cols will give a nasty segfault. I should just rewrite
    # get_rows_cols in Julia.
    pyccblade = PyCall.pyimport("ccblade.ccblade_jl")
    get_rows_cols = pyccblade.get_rows_cols

    input_data = [
        VarData("r", [1, num_radial], 1., "m")
        VarData("chord", [1, num_radial], 1., "m")
        VarData("theta", [1, num_radial], 1., "rad")
        VarData("Vx", [num_nodes, num_radial], 1., "m/s")
        VarData("Vy", [num_nodes, num_radial], 1., "m/s")
        VarData("rho", [num_nodes, 1], 1., "kg/m**3")
        VarData("mu", [num_nodes, 1], 1., "N/m**2*s")
        VarData("asound", [num_nodes, 1], 1., "m/s")
        VarData("pitch", [num_nodes, 1], 1., "rad")
        VarData("Rhub", [1, 1], 1., "m")
        VarData("Rtip", [1, 1], 1., "m")
        VarData("precone", [1, 1], 1., "rad")]

    output_data = [
        VarData("phi", [num_nodes, num_radial], 1., "rad"),
        VarData("Np", [num_nodes, num_radial], 1., "N/m"),
        VarData("Tp", [num_nodes, num_radial], 1., "N/m"),
        VarData("a", [num_nodes, num_radial], 1.),
        VarData("ap", [num_nodes, num_radial], 1.),
        VarData("u", [num_nodes, num_radial], 1., "m/s"),
        VarData("v", [num_nodes, num_radial], 1., "m/s"),
        VarData("alpha", [num_nodes, num_radial], 1., "rad"),
        VarData("W", [num_nodes, num_radial], 1., "m/s"),
        VarData("cl", [num_nodes, num_radial], 1.),
        VarData("cd", [num_nodes, num_radial], 1.),
        VarData("cn", [num_nodes, num_radial], 1.),
        VarData("ct", [num_nodes, num_radial], 1.),
        VarData("F", [num_nodes, num_radial], 1.),
        VarData("G", [num_nodes, num_radial], 1.)]

    partials_data = Array{PartialsData, 1}()

    rows, cols = get_rows_cols(
        of_shape=(num_nodes, num_radial), of_ss="ij",
        wrt_shape=(1,), wrt_ss="k")
    of_names = ["phi", "Np", "Tp", "a", "ap", "u", "v", "alpha", "W", "cl", "cd", "cn", "ct", "F", "G"]
    for name in of_names
        push!(partials_data, PartialsData(name, "Rhub", rows=rows, cols=cols))
        push!(partials_data, PartialsData(name, "Rtip", rows=rows, cols=cols))
        push!(partials_data, PartialsData(name, "precone", rows=rows, cols=cols))
    end

    rows, cols = get_rows_cols(
        of_shape=(num_nodes, num_radial), of_ss="ij",
        wrt_shape=(num_radial,), wrt_ss="j")
    for name in of_names
        push!(partials_data, PartialsData(name, "r", rows=rows, cols=cols))
        push!(partials_data, PartialsData(name, "chord", rows=rows, cols=cols))
        push!(partials_data, PartialsData(name, "theta", rows=rows, cols=cols))
    end

    rows, cols = get_rows_cols(
        of_shape=(num_nodes, num_radial), of_ss="ij",
        wrt_shape=(num_nodes,), wrt_ss="i")
    for name in of_names
        push!(partials_data, PartialsData(name, "rho", rows=rows, cols=cols))
        push!(partials_data, PartialsData(name, "mu", rows=rows, cols=cols))
        push!(partials_data, PartialsData(name, "asound", rows=rows, cols=cols))
        push!(partials_data, PartialsData(name, "pitch", rows=rows, cols=cols))
    end

    rows, cols = get_rows_cols(
        of_shape=(num_nodes, num_radial), of_ss="ij",
        wrt_shape=(num_nodes, num_radial), wrt_ss="ij")
    for name in of_names
        push!(partials_data, PartialsData(name, "Vx", rows=rows, cols=cols))
        push!(partials_data, PartialsData(name, "Vy", rows=rows, cols=cols))
        push!(partials_data, PartialsData(name, "phi", rows=rows, cols=cols))
    end

    for name in ["Np", "Tp", "a", "ap", "u", "v", "alpha", "W", "cl", "cd", "cn", "ct", "F", "G"]
        push!(partials_data, PartialsData(name, name, rows=rows, cols=cols, val=1.0))
    end

    return CCBladeResidualComp(num_nodes, num_radial, af, B, turbine, debug_print, input_data, output_data, partials_data)
end

function OpenMDAO.apply_nonlinear!(self::CCBladeResidualComp, inputs, outputs, residuals)
    num_nodes = self.num_nodes
    num_radial = self.num_radial

    # Outputs.
    n = 1
    out = Outputs(
        outputs["Np"][n, :],
        outputs["Tp"][n, :],
        outputs["a"][n, :],
        outputs["ap"][n, :],
        outputs["u"][n, :],
        outputs["v"][n, :],
        outputs["phi"][n, :],
        outputs["alpha"][n, :],
        outputs["W"][n, :],
        outputs["cl"][n, :],
        outputs["cd"][n, :],
        outputs["cn"][n, :],
        outputs["ct"][n, :],
        outputs["F"][n, :],
        outputs["G"][n, :])

    phis = outputs["phi"]

    for n in 1:num_nodes

        # Rotor parameters.
        rotor = Rotor(
            inputs["r"][n, :],
            inputs["chord"][n, :],
            inputs["theta"][n, :],
            fill(self.af, num_radial),
            inputs["Rhub"][1],
            inputs["Rtip"][1],
            self.B,
            self.turbine,
            inputs["precone"][1])

        # Operating point parameters.
        op_point = OperatingPoint(
             inputs["Vx"][n, :],
             inputs["Vy"][n, :],
             inputs["pitch"][n],
             inputs["rho"][n],
             inputs["mu"][n],
             inputs["asound"][n])

        for idx in 1:num_radial
            residuals["phi"][n, idx] = residual!(phis[n, idx], rotor, op_point, out, idx, true)
        end

        @. residuals["Np"][n, :] = outputs["Np"][n, :] - out.Np
        @. residuals["Tp"][n, :] = outputs["Tp"][n, :] - out.Tp
        @. residuals["a"][n, :] = outputs["a"][n, :] - out.a
        @. residuals["ap"][n, :] = outputs["ap"][n, :] - out.ap
        @. residuals["u"][n, :] = outputs["u"][n, :] - out.u
        @. residuals["v"][n, :] = outputs["v"][n, :] - out.v
        @. residuals["alpha"][n, :] = outputs["alpha"][n, :] - out.alpha
        @. residuals["W"][n, :] = outputs["W"][n, :] - out.W
        @. residuals["cl"][n, :] = outputs["cl"][n, :] - out.cl
        @. residuals["cd"][n, :] = outputs["cd"][n, :] - out.cd
        @. residuals["cn"][n, :] = outputs["cn"][n, :] - out.cn
        @. residuals["ct"][n, :] = outputs["ct"][n, :] - out.ct
        @. residuals["F"][n, :] = outputs["F"][n, :] - out.F
        @. residuals["G"][n, :] = outputs["G"][n, :] - out.G

    end

end

function OpenMDAO.linearize!(self::CCBladeResidualComp, inputs, outputs, partials)
    num_nodes = self.num_nodes
    num_radial = self.num_radial

    for n in 1:num_nodes

        # Rotor parameters.
        rotor = Rotor(
            inputs["r"][n, :],
            inputs["chord"][n, :],
            inputs["theta"][n, :],
            fill(self.af, num_radial),
            inputs["Rhub"][1],
            inputs["Rtip"][1],
            self.B,
            self.turbine,
            inputs["precone"][1])

        # Operating point parameters.
        op_point = OperatingPoint(
             inputs["Vx"][n, :],
             inputs["Vy"][n, :],
             inputs["pitch"][n],
             inputs["rho"][n],
             inputs["mu"][n],
             inputs["asound"][n])

        for idx in 1:num_radial
            phi = outputs["phi"][n, idx]
            derivs = residual_partials(phi, rotor, op_point, outputs, idx)
            for (k, v) in derivs
                p = reshape(partials["phi", String(k)], num_nodes, num_radial)
                p[n, idx] = v
            end

            derivs = output_partials(phi, rotor, op_point, outputs, idx)
            for ((of, wrt), v) in derivs
                p = reshape(partials[String(of), String(wrt)], num_nodes, num_radial)
                p[n, idx] = -v
            end
        end
    end

end

function OpenMDAO.solve_nonlinear!(self::CCBladeResidualComp, inputs, outputs)
    num_nodes = self.num_nodes
    num_radial = self.num_radial

    for n in 1:num_nodes

        # Rotor parameters.
        rotor = Rotor(
            inputs["r"][n, :],
            inputs["chord"][n, :],
            inputs["theta"][n, :],
            fill(self.af, num_radial),
            inputs["Rhub"][1],
            inputs["Rtip"][1],
            self.B,
            self.turbine,
            inputs["precone"][1])

        # Operating point parameters.
        op_point = OperatingPoint(
             inputs["Vx"][n, :],
             inputs["Vy"][n, :],
             inputs["pitch"][n],
             inputs["rho"][n],
             inputs["mu"][n],
             inputs["asound"][n])

        # Do it.
        out = solve(rotor, op_point)

        # Set the outputs.
        @. outputs["phi"][n, :] = out.phi
        @. outputs["Np"][n, :] = out.Np
        @. outputs["Tp"][n, :] = out.Tp
        @. outputs["a"][n, :] = out.a
        @. outputs["ap"][n, :] = out.ap
        @. outputs["u"][n, :] = out.u
        @. outputs["v"][n, :] = out.v
        @. outputs["alpha"][n, :] = out.alpha
        @. outputs["W"][n, :] = out.W
        @. outputs["cl"][n, :] = out.cl
        @. outputs["cd"][n, :] = out.cd
        @. outputs["cn"][n, :] = out.cn
        @. outputs["ct"][n, :] = out.ct
        @. outputs["F"][n, :] = out.F
        @. outputs["G"][n, :] = out.G

    end

end


text:
using ChangesOfVariables, InverseFunctions, ArraysOfArrays, Statistics
using Optimisers
using Plots

using EuclidianNormalizingFlows: CenterStretch, JohnsonTrafo, HouseholderTrafo
using EuclidianNormalizingFlows: optimize_whitening

f_fwd_true =
    HouseholderTrafo([1.0, 0.3]) ∘
    CenterStretch([4.0, 4.1], [2.0, 2.1], [3.0, 3.1]) # ∘
    # JohnsonTrafo([10.0, 11.0], [3.5, 3.6], [10.0, 11.0], [1.0, 1.1])

XW = randn(2, 10^5)
histogram2d(XW[1,:], XW[2,:], nbins = 100)

X = f_fwd_true(XW)
histogram2d(X[1,:], X[2,:], nbins = 100)

initial_trafo =
    # inverse(JohnsonTrafo(rand(2), rand(2), rand(2), rand(2))) ∘
    inverse(CenterStretch([0.0, 0.0], [1.0, 1.0], [0.0, 0.0])) ∘
    inverse(HouseholderTrafo(randn(2)))

optimizer = ADAGrad()

smpls = nestedview(X)
nbatches = 1000
nepochs = 10

r = optimize_whitening(smpls, initial_trafo, optimizer, nbatches = nbatches, nepochs = nepochs)

XW2 = r.result(X)
histogram2d(XW2[1,:], XW2[2,:], nbins = 100)

cov(XW2, dims = 2)

using EuclidianNormalizingFlows: mvnormal_negll_trafo, mvnormal_negll_trafograd
batches = collect(Iterators.partition(smpls, round(Int, length(smpls)/nbatches)))
plot(r.negll_history, label = "batch negative log-likelihood")
ref_neg_ll = mean([mvnormal_negll_trafo(inverse(f_fwd_true), flatview(b)) for b in batches])
hline!([ref_neg_ll], label = "initial")
ref_neg_ll = mean([mvnormal_negll_trafo(initial_trafo, flatview(b)) for b in batches])
hline!([ref_neg_ll], label = "target")


text:
using MPIReco

@testset "sparse single- and multi-channel in-memory reconstruction" begin
  bSF = MPIFile("./data/HeadScanner/systemMatrix.mdf")
  b = MPIFile("./data/HeadScanner/measurement.mdf")

  numPeriodAverages = 65
  numPatches = div(acqNumPeriodsPerFrame(bSF), numPeriodAverages)
  bgCorrection = false

  # MP Reco
  numPeriodGrouping = 1
  maxMixingOrder = -1

  @time c1 = reconstruction(bSF, b, frames=1:10, numAverages=10,
           numPeriodGrouping=numPeriodGrouping, numPeriodAverages=numPeriodAverages,
           minFreq=30e3, maxFreq=600e3, maxMixingOrder=maxMixingOrder,
           lambd=0.01, iterations=100, bgCorrectionInternal=bgCorrection,
           spectralLeakageCorrection=false)

  exportImage("./img/Cartesian1.png", arraydata(c1[1,:,:,1,1]))

  # SP Reco
  numPeriodGrouping = numPatches
  maxMixingOrder = 12

  @time c2 = reconstruction(bSF, b, frames=1:10, numAverages=10,
           numPeriodGrouping=numPeriodGrouping, numPeriodAverages=numPeriodAverages,
           minFreq=30e3, maxFreq=600e3, maxMixingOrder=maxMixingOrder,
           lambd=0.01, iterations=100, bgCorrectionInternal=bgCorrection,
           spectralLeakageCorrection=false)

  exportImage("./img/Cartesian2.png", arraydata(c2[1,:,:,1,1]))

  ####  Low Level ####

  N = calibSize(bSF)

  ## MP
  numPeriodGrouping = 1
  maxMixingOrder = -1

  freq = filterFrequencies(bSF, numPeriodGrouping=numPeriodGrouping, numPeriodAverages=numPeriodAverages,
             minFreq=30e3, maxFreq=600e3,  maxMixingOrder=maxMixingOrder)

  S, grid = getSF(bSF, freq, numPeriodAverages=numPeriodAverages, bgCorrection=bgCorrection,
           numPeriodGrouping=numPeriodGrouping, spectralLeakageCorrection=false);

  u = getMeasurementsFD(b, frequencies=freq, bgCorrection=bgCorrection, frames=1:10, numAverages=10,
               numPeriodGrouping=numPeriodGrouping, numPeriodAverages=numPeriodAverages,
               spectralLeakageCorrection=false)

  c3 = reshape(reconstruction(transpose(S), u, lambd=0.01, iterations=100), N[1], N[2])

  exportImage("./img/Cartesian3.png", c3)

  ## SP
  numPeriodGrouping = numPatches
  maxMixingOrder = 12

  freq = filterFrequencies(bSF, numPeriodGrouping=numPeriodGrouping, numPeriodAverages=numPeriodAverages,
             minFreq=30e3, maxFreq=600e3,  maxMixingOrder=maxMixingOrder)

  S, grid = getSF(bSF, freq, numPeriodAverages=numPeriodAverages, bgCorrection=bgCorrection,
           numPeriodGrouping=numPeriodGrouping, spectralLeakageCorrection=false);

  u = getMeasurementsFD(b, frequencies=freq, bgCorrection=bgCorrection, frames=1:10, numAverages=10,
               numPeriodGrouping=numPeriodGrouping, numPeriodAverages=numPeriodAverages,
               spectralLeakageCorrection=false)

  c4 = reshape(reconstruction(transpose(S), u, lambd=0.01, iterations=100), N[1], N[2])

  exportImage("./img/Cartesian4.png", c4)

  ####  Multi Color ####

  @time c5 = reconstruction(MultiContrastFile([bSF,bSF]), b, frames=1:10, numAverages=10,
           numPeriodGrouping=numPeriodGrouping, numPeriodAverages=numPeriodAverages,
           minFreq=30e3, maxFreq=600e3, maxMixingOrder=maxMixingOrder,
           lambd=0.01, iterations=100, bgCorrectionInternal=bgCorrection,
           spectralLeakageCorrection=false)

  exportImage("./img/Cartesian5.png", arraydata(c5[1,:,:,1,1]))

end


text:
const _get_gpt_config = Dict(
    "OpenAIftlm" => Dict(
        :host => :openai,
        :vocab => 40000,
        :layer => 12,
        :hidden => 768,
        :head => 12,
        :bpe => :bpe,
        :name => "OpenAIftlm",
        :items => (:gpt_model, :bpe, :vocab, :tokenizer)
    )
)

function gpt_init()
    for (model, url, cksum) ∈ zip(
        ("OpenAIftlm",
         ),
        ("https://drive.google.com/uc?export=download&id=1jZ4wSrR7rGXprKz1a3O8PpEQ2o-nnMEA",
         ),
        ("e372ab1f42fff84efc6e97d30d0aadd823706e9b6d4084a8c474ba68c494e22c",
         )
    )
        config = _get_gpt_config[model]
        head = config[:head]
        hidden = config[:hidden]
        layer = config[:layer]
        vocab = config[:vocab]
        bpe = config[:bpe] == :bpe ? " bpe" : ""
        band = config[:host] == :openai ? "OpenAI released " : ""
        register(DataDep(
            "GPT-$model",
            """
            $(band)gpt model with $(layer)-layer, $(head)-head, $(hidden)-hidden, $(vocab)$(bpe) vocabularies
            """,
            url,
            cksum;
            fetch_method = download_gdrive
        ))
    end
end


text:
# Basic list

abstract List{T}

type ListNil{T} <: List{T} end

copy{T}(l::ListNil{T}) = l

type ListNode{T} <: List{T}
    head::T
    tail::List{T}
end

copy{T}(l::ListNode{T}) = ListNode{T}(l.head, l.tail)

head(l::ListNode) = l.head
tail(l::ListNode) = l.tail

# iterator
function start(l::List)
    l
end

function next(::List, l::List)
    (l.head, l.tail)
end

function done(::List, l::List)
    typeof(l) <: ListNil
end

function cons{T}(value, l::List{T})
    ListNode{T}(value, l)
end

function length{T}(l::List{T})
    n = 0
    while typeof(l) != ListNil{T}
        n += 1
        l = l.tail
    end
    n
end

function cat{T}(a::List{T}, b::List{T})
    if a === nothing
        b
    else
        a = copy(a)
        u = a
        while !is(u.tail, nothing)
            u.tail = copy(u.tail)
            u = u.tail
        end
        u.tail = b
        a
    end
end

function list{T}(xs::AbstractArray{T})
    if length(xs) == 0
        return nothing
    else
        u = v = ListNode(xs[1], ListNil{T}())
        i = 2
        while i <= length(xs)
            v.tail = ListNode{T}(xs[i], ListNil{T}())
            v = v.tail
            i += 1
        end
        u
    end
end

function show{T}(io::IO, a::List{T})
    print(io, "List([")
    while typeof(a) != ListNil{T}
        print(io, a.head)
        if typeof(a.tail) != ListNil{T}
            print(io, ", ")
        end
        a = a.tail
    end
    print(io, "])")
end


text:
#=

   Sequence waiting times.

   From
   Gunnar Blom, Lars Holst, Dennis Sandell:
   "Problems and Snapshots from the World of Probability"
   Page 4ff, Problem 1.4 Patterns I, (Part Problem 1)

   Problem: Waiting time for a pattern (coin flipping)

   What is the waiting time for the patterns[0,1] in a coin
   flipping sequence? For the pattern [1,1]?
   Answer:
     * waiting time for [0,1] = 4
     * waiting time for [1,1] = 6

   This model verifies this result:
   - [0,1]: 4.0031
   - [1,1]: 6.0300

   Mean length of some other sequences:
   - [0,1,0]: 9.9668
   - [1,1,1]: 14.1199

   Cf ~/webppl/sequence_waiting_times_1.wppl

=#
using Turing, StatsPlots
include("jl_utils.jl")

@model function sequence_waiting_times_1(pattern)
    patternLen = length(pattern)

    function flipSequence(a)
        len = length(a)
        if len >= patternLen && a[len-patternLen+1:end] == pattern
            return a
        else
            return flipSequence(push!(a,rand(Bernoulli(0.5))*1))
        end
    end

    a = flipSequence([])
    len ~ Dirac(length(a))

end

function run_seqence_waiting_times_1(pattern)
    println("\npattern: $pattern")
    model = sequence_waiting_times_1(pattern)

    # chns = sample(model, Prior(), 10_000)
    # chns = sample(model, MH(), 10_000)
    chns = sample(model, PG(5), 10_000)
    # chns = sample(model, SMC(), 10_000)
    # chns = sample(model, IS(), 10_000)

    display(chns)
    # display(plot(chns))

    # show_var_dist_pct(chns,:len)
end

run_seqence_waiting_times_1([0,1])
run_seqence_waiting_times_1([1,1])
run_seqence_waiting_times_1([0,1,0])
run_seqence_waiting_times_1([1,1,1])

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
using LengthChannels
using Test

@testset "LengthChannels.jl" begin
    l,b = 10, 4

    @testset "Typed Constructors" begin
        @info "Testing Typed Constructors"

        lc = LengthChannel{Int}(l, b) do ch
            for i = 1:100
                put!(ch, i)
            end
        end

        @test fullsizeof(lc) == l*sizeof(1)

        @test eltype(lc) <: Int

        @test length(lc) == l
        cc = collect(lc)
        @test length(cc) == l
        @test cc == 1:l
        @test isopen(lc)

        lc = LengthChannel{Int}(l, b, spawn=true) do ch
            for i = 1:100
                put!(ch, i)
            end
        end

        @test length(lc) == l
        cc = collect(lc)
        @test length(cc) == l
        @test cc == 1:l
        @test isopen(lc)

    end

    @testset "Untyped constructors" begin
        @info "Testing Untyped constructors"

        lc = LengthChannel(l, b) do ch
            for i = 1:100
                put!(ch, i)
            end
        end

        @test eltype(lc) == Any
        @test length(lc) == l
        cc = collect(lc)
        @test length(cc) == l
        @test cc == 1:l
        @test isopen(lc)

        lc = LengthChannel(l, b)
        @test length(lc) == l
        put!(lc, 1)
        @test length(lc) == l

        @test take!(lc) == 1

        @test !isready(lc)

    end

    @testset "Autoclose" begin
        @info "Testing Autoclose"

        lc = LengthChannel{Int}(l, b, autoclose=true) do ch
            for i = 1:100
                put!(ch, i)
            end
        end

        @test eltype(lc) <: Int

        @test length(lc) == l
        cc = collect(lc)
        @test length(cc) == l
        @test cc == 1:l
        @test !isopen(lc)

        lc = LengthChannel{Int}(l, b, autoclose=true, spawn=true) do ch
            for i = 1:100
                put!(ch, i)
            end
        end

        @test length(lc) == l
        cc = collect(lc)
        @test length(cc) == l
        @test cc == 1:l
        @test !isopen(lc)

    end

    @testset "Manual interaction" begin
        @info "Testing Manual interaction"

        lc = LengthChannel{Int}(l, b)
        @test length(lc) == l
        put!(lc, 1)
        @test length(lc) == l

        @test take!(lc) == 1

        @test !isready(lc)

    end

    @testset "Iteration" begin
        @info "Testing Iteration"

        lc = LengthChannel{Int}(l, b) do ch
            for i = 1:100
                put!(ch, i)
            end
        end

        c = 0
        for e in lc
            c += 1
            @test e == c
        end
        @test c == l

    end

    @testset "Channel closes early" begin
        @info "Testing Channel closes early"

        lc = LengthChannel{Int}(l, b) do ch
            for i = 1:5
                put!(ch, i)
            end
        end

        @test length(lc) == l
        if VERSION > v"1.3.100"
            @test_throws InvalidStateException collect(lc)
            @test !isopen(lc)
        end

        lc = LengthChannel{Int}(l, b) do ch
            for i = 1:5
                put!(ch, i)
            end
        end

        for e in 1:5
            @test e == first(lc)
        end
        sleep(0.1) # it seems to take a while to close it
        VERSION > v"1.3.100" && @test !isopen(lc)

    end

    @testset "Wrapping LengthChannel" begin
        @info "Testing Wrapping LengthChannel"

        lc = LengthChannel{Int}(l, b) do ch
            for i = 1:100
                put!(ch, i)
            end
        end
        wlc = LengthChannel{Float64}(x->Float64(x^2), lc)

        @test eltype(wlc) <: Float64

        @test length(wlc) == l
        cc = collect(wlc)
        @test length(cc) == l
        @test cc == (1:l).^2
        @test isopen(wlc)

        lc = LengthChannel{Int}(l, b) do ch
            for i = 1:100
                put!(ch, i)
            end
        end
        f(x) = x^2
        f(::Type) = Float64
        wlc = LengthChannel(f, lc)
        @test eltype(wlc) <: Float64

        @test length(wlc) == l
        cc = collect(wlc)
        @test length(cc) == l
        @test cc == (1:l).^2
        @test isopen(wlc)
    end
end


text:
using Distributed

function problem(x1::Float64, x2::Float64)
    x1^2 + x2^2
end

const Maximization = >
const Minimization = <

mutable struct Target
    problem::Function
    best_than::Function
    xdefzone::Array{Tuple{Float64,Float64},1}

    global_best::Vector{Float64}     # 截止目前，所有粒子获得的最优解
    global_value::Float64

    Target(p::Function, compare::Function, n::Int64, xdef::Array{Tuple{Float64,Float64},1})  =
            new(p, compare, xdef, Vector{Float64}(undef, n), ifelse(compare == Maximization, -Inf, Inf))
end

mutable struct SwarmParams
    population::Int64      # 群体数量
    iteration::Int64       # 迭代次数
    c1::Float64             # 学习因子c1
    c2::Float64             # 学习因子c2
    w::Float64              # 惯性因子

    SwarmParams(p::Int64, i::Int64, c1::Float64, c2::Float64, w::Float64) = new(p, i, c1, c2, w)
end

mutable struct Particle
    position::Vector{Float64}       # 粒子的当前位置
    velocity::Vector{Float64}       # 粒子的当前速度

    local_best_position::Vector{Float64}     # 截止目前，该粒子获得的最优解
    local_best_value::Float64

    Particle(n::Int64) = new(zeros(Float64,n), zeros(Float64,n), zeros(Float64,n), 0.0)
end

function update_particle(p::Particle, sc::SwarmParams, t::Target)
    p.velocity = sc.w .* p.velocity + sc.c1 * rand() .* (p.local_best_position .- p.position)+ sc.c2 * rand() .* (t.global_best .- p.position)

    @sync @distributed for i in 1:length(p.velocity)
        if p.velocity[i] > target.xdefzone[i][2] || p.velocity[i] < target.xdefzone[i][1]
            p.velocity[i] = ifelse(rand() > 0.5, -1.0, 1.0) * rand()
        end
    end

    p.position = p.position .+ p.velocity

    @sync @distributed for i in 1:length(p.position)
        if p.position[i] > target.xdefzone[i][2] || p.position[i] < target.xdefzone[i][1]
            p.position[i] = rand() .* (target.xdefzone[i][2]-target.xdefzone[i][1]) .+ target.xdefzone[i][1]
        end
    end
end

function optimize(target::Target, sc::SwarmParams)
    @show target.xdefzone

    dimension = length(target.global_best)

    swarm = Vector{Particle}()
    for i in 1:sc.population
        push!(swarm,Particle(dimension))
    end

    @sync @distributed for p in swarm
        p.velocity = rand(dimension)
        p.position = rand(dimension)

        for i in 1:length(p.position)
            p.position[i] = p.position[i] .* (target.xdefzone[i][2]-target.xdefzone[i][1]) .+ target.xdefzone[i][1]
        end

        p.local_best_position = p.position

        p.local_best_value = target.problem(p.position...)
    end

    println(">>> ---- starting ----")

    for i in 1:sc.iteration
        @show i
        @sync @distributed for p in swarm
            value = target.problem(p.position...)
            @show value

            if target.best_than(value, p.local_best_value)
                p.local_best_value = value
                p.local_best_position = p.position
            end
        end

        @sync @distributed for p in swarm
            if target.best_than(p.local_best_value, target.global_value)
                target.global_value = p.local_best_value
                target.global_best = p.position
            end
        end

        @sync @distributed for p in swarm
            update_particle(p, sc, target)
        end
    end
end

target = Target(problem, <, 2, [(-10.0, 10.0), (-10.0, 10.0)])
sc = SwarmParams(20, 500, 2.0, 2.0, 1.0)

optimize(target, sc)

println(target.global_value)
println(target.global_best)


text:
export get_residue,
       get_sourceside_wavefield,
       velocity_gradient

include("velocity_gradient.jl")


text:
# Functions to Test `solveBendersMilp`
# Edward J. Xu, edxu96@outlook.com
# Aug 11, 2019

"Test 1: mat_a and mat_b being column vector"
function do_test_1()
    n_x = 1
    n_y = 1
    vec_min_y = hcat([0])
    vec_max_y = hcat([10])
    vec_c = hcat([5])
    vec_f = hcat([-3])
    vec_b = hcat([4; 0; -13])
    mat_a = hcat([1; 2; 1])
    mat_b = hcat([2; -1; -3])
    solveBendersMilp(n_x, n_y, vec_min_y, vec_max_y, vec_c, vec_f, vec_b, mat_a, mat_b)
end

"Test 2: mat_a and mat_b being matrix"
function do_test_2()
    n_x = 2
    n_y = 2
    vec_min_y = hcat([0; 0])
    vec_max_y = hcat([10; 10])
    vec_c = hcat([5; 3])
    vec_f = hcat([-3; 1])
    vec_b = hcat([4; 0; -13])
    mat_a = [1 3; 2 1; 1 -5]
    mat_b = [2 -4; -1 2; -3 1]
    solveBendersMilp(n_x, n_y, vec_min_y, vec_max_y, vec_c, vec_f, vec_b, mat_a, mat_b)
end

"Test 3: sub and ray problem"
function do_test_3()
    n_x = 2
    n_y = 2
    vec_min_y = hcat([0; 0])
    vec_max_y = hcat([2; 2])
    vec_c = hcat([2; 6])
    vec_f = hcat([2; 3])
    vec_b = hcat([5; 4])
    mat_a = [-1 2; 1 -3]
    mat_b = [3 -1; 2 2]
    solveBendersMilp(n_x, n_y, vec_min_y, vec_max_y, vec_c, vec_f, vec_b, mat_a, mat_b)
end

function test_benders()
    do_test_1()
    do_test_2()
    do_test_3()
    return 0
end


text:
using AlphaZero

include("games.jl")
GAME = get(ENV, "GAME", "tictactoe")
SelectedGame = GAME_MODULE[GAME]
using .SelectedGame: Game, Training

const SESSION_DIR = joinpath("sessions", GAME)

session = Session(
  Game, Training.Network{Game},
  Training.params, Training.netparams, benchmark=Training.benchmark,
  dir=SESSION_DIR, autosave=false, save_intermediate=false)

start_explorer(session)


text:
function detect_anagrams(subject, candidates)
    sorted(s) = sort(collect(lowercase(s)))

    subject = lowercase(subject)
    subject_sorted = sorted(subject)

    function check(candidate)
        candidate = lowercase(candidate)
        (candidate != subject) && (sorted(candidate) == subject_sorted)
    end

    filter(check, candidates)
end

max_stars_repo_name: twolodzko/Learning

max_issues_repo_name: twolodzko/Learning

max_forks_repo_name: twolodzko/Learning


text:
using SparseArrays
using GalerkinSparseGrids

# The exact coefficients for f(x_1, ..., x_D) = 1 in the modal basis
function get_one_modal(D::Int, k::Int, n::Int)
    len = get_size(Val(D), k, n)
    return [i==1 ? 1.0 : 0.0 for i in 1:len]
end

# The exact coefficients of the function f(x_1, ..., x_D) = x_i
function get_xi_modal(D::Int, i::Int, k::Int, n::Int)
    len = get_size(Val(1), k, n)
    coeffs1_1D::Array{Float64, 1} = [i==1 ? 1 : 0 for i in 1:len] # Coeffs for 1
    coeffs2_1D::Array{Float64, 1} = [i==2 ? 1/(sqrt(3)) : 0 for i in 1:len]
    coeff_array = [d==i ? coeffs2_1D : coeffs1_1D for d in 1:D]
    return tensor_construct(D, k, n, coeff_array)
end

function get_xi2_modal(D::Int, i::Int, k::Int, n::Int)
    len = get_size(Val(1), k, n)
    # Coefficients for 1:
    coeffs1_1D::Array{Float64, 1} = [i==1 ? 1 : 0 for i in 1:len]
    # Coefficients for (v-1/2)^2:
    coeffs2_1D::Array{Float64, 1} = [i==3 ? 2/(3*sqrt(5)) : 0 for i in 1:len]
    coeffs2_1D += 1/3 * coeffs1_1D
    coeff_array = [d==i ? coeffs2_1D : coeffs1_1D for d in 1:D]
    return tensor_construct(D, k, n, coeff_array)
end

# The exact coefficients for f(x) = |x|^2 in the modal basis
function get_r2_modal(D::Int, k::Int, n::Int)
    coeff_array = [get_xi2_modal(D, i, k, n) for i in 1:D]

    return sum(coeff_array)
end

# The exact coefficients for f(x_1, ..., x_D) = x_i = x_i in the point basis
# ^* insofar as the transformation matrices are exact
function get_xi_point(D::Int, i::Int, k::Int, n::Int,
                        m2n::SparseMatrixCSC{Float64, Int},
                        n2p::SparseMatrixCSC{Float64, Int})
    return n2p * (m2n * get_xi_modal(D, i, k, n))
end

# The exact coefficients for f(x_1, ..., x_D) = |x|^2 in the point basis
# ^* insofar as the transformation matrices are exact
function get_r2_point(D::Int, k::Int, n::Int,
                        m2n::SparseMatrixCSC{Float64, Int},
                        n2p::SparseMatrixCSC{Float64, Int})
    v = n2p * (m2n * get_r2_modal(D, k, n))
    v[v .< 0] .= 0
    return v
end


text:
module KATE

using Flux
using Flux: onehot, argmax, chunk, batchseq, throttle, crossentropy, sigmoid, treelike, glorot_uniform
using Flux.Tracker: data
using StatsBase: wsample
using Base.Iterators: partition

struct KCompetetive{F,S,T}
  σ::F
  W::S
  b::T
  α::Float64 # boost
end

function KCompetetive(in::Integer, k::Integer, σ = tanh, α::Float64 = 6.26;
               initW = glorot_uniform, initb = zeros)
  if k > in
    throw(warn("Warning: k (out) should not be larger than dim: $in, found: $k, using $in"))
    k = in
  end
  return KCompetetive(σ, param(initW(k, in)), param(initb(k)), α)
end

treelike(KCompetetive)

function (a::KCompetetive)(x)
  W, b, σ, α = a.W, a.b, a.σ, a.α

  ps = Vector{Tuple{Int64,Float64}}()
  ns = Vector{Tuple{Int64,Float64}}()

  for (i, activation) = enumerate(data(x))
    if activation >= 0
        push!(ps, (i, activation))
    else
        push!(ns, (i, activation))
    end
  end
  ps = sort(ps, by=last)
  ns = sort(ns, by=last, rev=true)
  P = length(ps)
  N = length(ns)
  k = first(size(a.W))

  z = data(x)

  p = P-Int(k/2)
  if p > 0
    Epos = sum(map(last, ps[1:p]))
    # @show Epos
    for i = p+1:P
        # positive winners
        z[first(ps[i])] += α * Epos
    end
    for i = 1:p
        # positive losers
        z[first(ps[i])] = 0.0
    end
  end

  n = N-Int(k/2)
  if n > 0
    Eneg = sum(map(last, ns[1:n]))
    # @show Eneg
    for i = n+1:N
        # negative winners
        z[first(ns[i])] += α * Eneg
    end
    for i = 1:n
        # negative losers
        z[first(ns[i])] = 0.0
    end
  end

  #  result = σ.(W*x .+ b)
  #  @show result
  #  result
  σ.(W*x .+ b)
end

function Base.show(io::IO, l::KCompetetive)
  print(io, "KCompetetive(", size(l.W, 2), ", ", size(l.W, 1))
  l.σ == identity || print(io, ", ", l.σ)
  print(io, ")")
end

function non_empty_string(s)
    s != ""
end

function count_words(text::Array{String})
    unique_words = unique(text)
    word_counts = Dict{String,Int64}(map(word -> (word, 0), unique_words))
    map(word -> word_counts[word] += 1, text)
    word_counts
end

function nl(word::String, wc::Dict{String,Int64}, V::Int)
    n = wc[word]
    logwc = log(1+n)
    logwc/V*logwc
end

function normalize_log(text::Array{String}, word_counts::Dict{String,Int64})
    V = length(word_counts)
    map(word -> nl(word,word_counts,V), text)
end

function transform_text_to_input(src::String, limit::Int=1000)
    text_raw = readstring(src)
    text = map(s -> String(s),
        filter(non_empty_string, split(text_raw, r"[^\wäÄöÖüÜ&]+")))
    text = text[1:limit]
    word_counts = count_words(text)
    input = normalize_log(text, word_counts)
    return input, word_counts
end

function get_similar_words(model, query_id, vocab; topn=10)
    # @show vocab
    W = model[1].W.data
    W = W/norm(W)
    query = W[query_id]
    # @show query
    score = query*(W')
    # @show score
    # @show size(score), typeof(score)
    # @show score[:,1]
    # @show size(score[:,1]), typeof(score[:,1])
    vidx = sort(score[:,1])[1:topn]
    # @show vidx
    # weights = unitmatrix(weights) # normalize
    # query = weights[query_id]
    # score = query.dot(weights.T)
    # vidx = score.argsort()[::-1][:topn]

    return [vocab[idx] for idx in vidx]
end

end # module KATE


text:
function pkg_dir_simplecontainergenerator()::String
    this_file = @__FILE__ # PACKAGE_ROOT/src/pkg_dir.jl
    src_directory = dirname(this_file) # PACKAGE_ROOT/src
    package_root = dirname(src_directory) # PACKAGE_ROOT
    return package_root
end


text:
@doc raw"""
    PoissonLikelihood(λ::Real)->PoissonLikelihood

## Arguments
- `λ::Real` : Maximal Poisson rate

---

[Poisson Likelihood](https://en.wikipedia.org/wiki/Poisson_distribution) where a Poisson distribution is defined at every point in space (careful, it's different from continous Poisson processes).
```math
    p(y|f) = \text{Poisson}(y|\lambda \sigma(f))
```
Where ``\sigma`` is the logistic function.
Augmentation details will be released at some point (open an issue if you want to see them)
"""
function PoissonLikelihood(λ::Real)
    return PoissonLikelihood(ScaledLogistic(λ))
end

struct ScaledLogistic{T} <: AbstractLink
    λ::Vector{T}
end

ScaledLogistic(λ::Real) = ScaledLogistic([λ])

(l::ScaledLogistic)(f::Real) = l.λ[1] * logistic(f)

function implemented(
    ::PoissonLikelihood{<:ScaledLogistic}, ::Union{<:AnalyticVI,<:GibbsSampling}
)
    return true
end

function (l::PoissonLikelihood)(y::Real, f::Real)
    return pdf(l(f), y)
end

function Distributions.loglikelihood(l::PoissonLikelihood, y::Real, f::Real)
    return logpdf(l(f), y)
end

function Base.show(io::IO, l::PoissonLikelihood{<:ScaledLogistic})
    return print(io, "Poisson Likelihood (λ = $(l.invlink.λ[1]))")
end

function compute_proba(
    l::PoissonLikelihood, μ::AbstractVector{T}, σ²::AbstractVector{<:Real}
) where {T}
    N = length(μ)
    pred = zeros(T, N)
    sig_pred = zeros(T, N)
    for i in 1:N
        x = pred_nodes .* sqrt.(max(σ²[i], zero(σ²[i]))) .+ μ[i]
        pred[i] = dot(pred_weights, l.invlink.(x))
        sig_pred[i] = dot(pred_weights, l.invlink.(x) .^ 2) - pred[i]^2
    end
    return pred, sig_pred
end

### Local Updates ###
function init_local_vars(::PoissonLikelihood, batchsize::Int, T::DataType=Float64)
    return (; c=rand(T, batchsize), θ=zeros(T, batchsize), γ=rand(T, batchsize))
end

function local_updates!(
    local_vars,
    l::PoissonLikelihood{<:ScaledLogistic},
    y::AbstractVector,
    μ::AbstractVector,
    diagΣ::AbstractVector,
)
    λ = only(l.invlink.λ)
    map!(sqrt_expec_square, local_vars.c, μ, diagΣ)
    map!(local_vars.γ, μ, local_vars.c) do μ, c
        λ * safe_expcosh(-μ / 2, c / 2) / 2
    end
    map!(local_vars.θ, y, local_vars.γ, local_vars.c) do y, γ, c
        (y + γ) / c * tanh(c / 2)
    end
    l.invlink.λ .= sum(y) / sum(expectation.(logistic, μ, diagΣ))
    return local_vars
end

function sample_local!(
    local_vars, l::PoissonLikelihood{<:ScaledLogistic}, y::AbstractVector, f::AbstractVector
)
    map!(rand ∘ l, local_vars.γ, f) # sample n
    map!(local_vars.θ, y, local_vars.γ, f) do y, γ, f
        rand(PolyaGamma(y + Int(γ), abs(f))) # Sample ω
    end
    return local_vars
end

### Global Updates ###

@inline function ∇E_μ(
    ::PoissonLikelihood{<:ScaledLogistic}, ::AOptimizer, y::AbstractVector, state
)
    return ((y - state.γ) / 2,)
end
@inline function ∇E_Σ(
    ::PoissonLikelihood{<:ScaledLogistic}, ::AOptimizer, y::AbstractVector, state
)
    return (state.θ / 2,)
end

## ELBO Section ##
function expec_loglikelihood(
    l::PoissonLikelihood{<:ScaledLogistic},
    ::AnalyticVI,
    y,
    μ::AbstractVector,
    Σ::AbstractVector,
    state,
)
    tot = (dot(μ, (y - state.γ)) - dot(state.θ, abs2.(μ)) - dot(state.θ, Σ)) / 2
    tot += Zygote.@ignore(
        sum(y * log(l.invlink.λ[1])) - sum(logfactorial, y) - logtwo * sum((y + state.γ))
    )
    return tot
end

function AugmentedKL(l::PoissonLikelihood{<:ScaledLogistic}, state, y)
    return PoissonKL(l, state) + PolyaGammaKL(l, state, y)
end

function PoissonKL(l::PoissonLikelihood{<:ScaledLogistic}, state)
    return PoissonKL(state.γ, only(l.invlink.λ))
end

function PolyaGammaKL(::PoissonLikelihood{<:ScaledLogistic}, state, y)
    return PolyaGammaKL(y + state.γ, state.c, state.θ)
end


text:
# MIT License
#
# Copyright (c) 2020 Joshua E Gentry

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

struct Output
    chaining::CHAINING
    block::BLOCK
    block_len::UInt32
    counter::UInt64
    flags::UInt32
end

function chaining_value(self::Output)::CHAINING
    return compress(
        self.chaining,
        self.block,
        self.block_len,
        self.counter,
        self.flags
    )
end

function hash(
            self::Output,
            out::AbstractVector{UInt8}
        )::Nothing

    if length(out) == 32
        out = reinterpret(UInt32, out)
        # ---------------------------------------------------------------------------------
        # Normal case, produce a 256bit hash.
        cv = compress(
            self.chaining,
            self.block,
            self.block_len,
            UInt64(0),
            self.flags | ROOT
        )
        vstore(cv[1], out, 1)
        vstore(cv[2], out, 5)
    else
        # ---------------------------------------------------------------------------------
        # Extended case, produce however many bits' the caller wants.
        todo    = length(out)
        offset  = 1
        counter = UInt64(0)

        # Handle the full 64 byte chunks.
        while todo >= 64
            todo -= 64
            compress_xof(
                self.chaining,
                self.block,
                self.block_len,
                counter,
                self.flags | ROOT,
                out, offset
            )
            offset  += 64
            counter += 1
        end

        # Handle the last few bytes.
        if todo > 0
            buffer = Vector{UInt8}(undef, 64)
            compress_xof(
                self.chaining,
                self.block,
                self.block_len,
                counter,
                self.flags | ROOT,
                buffer, 1
            )
            copyto!(out, offset, buffer, 1, todo)
        end
    end
    nothing
end


text:
using JLD
using ProgressMeter
using Distributions
using StatsBase
include("make_SQ_model.jl")
include("utilities.jl")
include("self_confidence.jl")

function get_net_input_vals(net::Dict{Symbol,Any},inpts::Array{Symbol},sch::ML.Schema;return_scaled::Bool=false)
    info("# Getting inputs")
    ins = Array{Float64}(length(inpts),1)
    scaled_ins = Array{Float64}(length(inpts),1)
    for (i,n) in zip(collect(1:length(inpts)),inpts)
        ins[i] = net[:training_data][n]
        if return_scaled
            scaled_ins[i] = (ins[i]-mean(sch[n]))/std(sch[n])
        end
    end
    if return_scaled
        return ins, scaled_ins
    else
        return ins
    end
end

function get_surrogate(fname::String,epocs::Int64,fldr::String,net_type::String,
                       inpts::Array{Symbol},trusted_fname::String)
    # make an sq model if it doesn't exist already
    if !(any([contains(x,fname*"_") for x in readdir(fldr)]))
        info("No nn file exists, making one now")
        make_sq_model(net_type,inpts,num_epoc=epocs,trusted_fname=trusted_fname)
    end

    info("Loading xQ Model")
    param_files = searchdir(fldr,fname,".params")
    num_epocs = parse(split(match(r"-\d+",param_files[1]).match,"-")[2])

    fn = joinpath(fldr,fname)
    SQmodel = load_network(fn,num_epocs,fn*"_SQmodel.jld")

    return SQmodel
end

function calculate_sq(experiment_dict)
    trusted_cond = 0
    for c in keys(experiment_dict[:conditions])
        if experiment_dict[:conditions][c][:role] == :trusted
            trusted_cond = c
        end
    end
    for expr in keys(experiment_dict[:conditions])
        net_type = experiment_dict[:name]
        inpts = experiment_dict[:xQ][:inpts]
        epocs = experiment_dict[:xQ][:epocs]
        println("Processing: $(experiment_dict[:conditions][expr])")

        println("running '$expr' data")

        ###### Define the variables for the problem
        inputs = Dict()
        for i in inpts
            inputs[i] = "ML.Continuous"
        end

        log_fname = "$(net_type)_$(make_label_from_keys(inputs))"
        log_loc = experiment_dict[:xQ][:nn_loc]
        tfname = experiment_dict[:conditions][trusted_cond][:fname]
        tfold = experiment_dict[:conditions][trusted_cond][:fldr]
        cfname = experiment_dict[:conditions][expr][:fname]
        cfold = experiment_dict[:conditions][expr][:fldr]
        trusted_path = joinpath(tfold,tfname*".csv")

        ###### get surrogate model
        SQmodel = get_surrogate(log_fname,epocs,log_loc,net_type,inpts,trusted_path)
        info("restoring limits")
        limits = restore_eng_units(SQmodel.range,SQmodel.output_sch)
        println("limits: $limits")
        println("out_sch: $(SQmodel.output_sch)")

        ###### calculate xQ for each net
        df = JLD.load(joinpath(cfold,cfname*".jld"))
        d = df["problem_dict"]
        g = df["global_rwd_range"]
        info("Processing Networks")
        p = Progress(length(d),barglyphs=BarGlyphs("[=> ]"),color=:white)
        j = 0
        for net in d
            net_num = net[1]
            net = net[2]
            e_start = net[:evader_start]
            p_start = net[:pursuer_start]
            ins, scaled_ins = get_net_input_vals(net,inpts,SQmodel.input_sch,return_scaled=true)

            info("getting predictions")
            println("ins: $ins, $(typeof(ins))")
            _notused, R_star = SQ_predict(SQmodel,scaled_ins,use_eng_units=true)
            R_star_μ = R_star[:X3_1][1]
            R_star_σ = R_star[:X3_2][1]

            R_μ = net[:training_data][:X3_1] # mean of rwd dist
            R_σ = net[:training_data][:X3_2] # std of rwd dist

            info("R_star_μ = $R_star_μ, R_star_σ = $R_star_σ")
            info("R_μ = $R_μ, R_σ = $R_σ")
            # make sure sigma isn't zero or negative (i.e. model could have bad prediction, or data could be all an identical number)
            # is R_star <= 0. ? if yes R_star_σ = 1. : otherwise do nothing
            R_star_σ <= 0. ? R_star_σ = 1. : nothing
            R_σ <= 0. ? R_σ = 1. : nothing

            x_Q, x_Q_Dict = X3(Normal(R_μ,R_σ),Normal(R_star_μ,R_star_σ),global_rwd_range=g,return_raw_sq=true)

            net[:solver_quality] = Dict{Symbol,Any}()
            net[:solver_quality][:x_QDict] = deepcopy(x_Q_Dict)
            net[:solver_quality][:R_star] = R_star
            net[:solver_quality][:X_Q] = x_Q
            j += 1
            ProgressMeter.next!(p; showvalues = [(:network,"$j of $(length(d))")])
        end
        # re-write original file with xQ added to each network
        out_fname = joinpath("logs",net_type*"_$(expr)_"*make_label_from_keys(inputs)*"_solver.jld")
        info("# writing to $out_fname")
        jldopen(out_fname, "w") do file
            write(file, "problem_dict", d)
            write(file, "global_rwd_range", g)
        end
    end
end

max_stars_repo_name: COHRINT/FaMSeC

max_issues_repo_name: COHRINT/FaMSeC

max_forks_repo_name: COHRINT/FaMSeC


text:
module FiniteFloats

export AbstractFinite,
       Finite64, Finite32, Finite16,
       finite,
       typemaxneg, typeminneg,
       square, cube

import Base: hash, promote_rule, convert,
    string, show,
    typemax, typemin, floatmax, floatmin,
    significand, exponent, precision,
    unsigned, signed, float,
    prevfloat, nextfloat, isequal, isless,
    (==), (!=), (<), (<=), (>=), (>),
    (+), (-), (*), (/), (^),
    inv, div, rem, fld, mod, cld,
    round, trunc, ceil, floor,
    abs, signbit, copysign, flipsign, sign,
    frexp, ldexp, modf,
    min, max, minmax,
    clamp, sqrt, cbrt, hypot,
    exp, expm1, exp2, exp10,
    log, log1p, log2, log10,
    rad2deg, deg2rad, mod2pi, rem2pi,
    sin, cos, tan, csc, sec, cot,
    asin, acos, atan, acsc, asec, acot,
    sinh, cosh, tanh, csch, sech, coth,
    asinh, acosh, atanh, acsch, asech, acoth,
    sincos, sinc, sinpi, cospi,
    sind, cosd, tand, cscd, secd, cotd,
    asind, acosd, atand, acscd, asecd, acotd

import LinearAlgebra
import LinearAlgebra: adjoint, adjoint!, axpby!, axpy!,
          bunchkaufman, bunchkaufman!, cholesky, cholesky!,
          cond, condskeel, copy_transpose!,
          copyto!, cross, det, diag, diagind, diagm, dot,
          eigen, eigen!, eigmax, eigmin, eigvals,
          eigvals!, eigvecs, factorize, fill!,
          hessenberg, hessenberg!,
          isdiag, ishermitian, isposdef, isposdef!,
          issuccess, issymmetric, istril, istriu, kron, ldiv!,
          ldlt, ldlt!, lmul!, logabsdet, logdet,
          lowrankdowndate, lowrankdowndate!, lowrankupdate,
          lowrankupdate!, lq, lq!, lu, lu!,
          lyap, mul!, norm, normalize, normalize!, nullspace,
          opnorm, ordschur, ordschur!, pinv, qr, qr!,
          rank, rdiv!, rmul!, schur, schur!,
          svd, svd!, svdvals, svdvals!, sylvester,
          tr, transpose, transpose!,
          tril, tril!, triu, triu!

include("types.jl")

include("linearalgebra_types.jl")
include("linearalgebra.jl")

end # FiniteFloats


text:
# loss functions
function loss(w, srnn, visual, x, y, batchsizes; o=Dict(), values=[])
    finetune = get(o, :finetune, false)
    atype = typeof(AutoGrad.getval(w["wsoft"]))
    visual = convert(atype, visual)
    if finetune && haskey(w, "wcnn")
        visual = vgg16(w["wcnn"], visual; o=o)
    end
    lossval, nwords = decoder(w, srnn, visual, x, y, batchsizes; o=o)
    push!(values, AutoGrad.getval(lossval), AutoGrad.getval(nwords))
    return lossval/nwords
end

# loss gradient functions
lossgradient = grad(loss)

# loss function for decoder network
function decoder(w, srnn, vis, x, y, batchsizes=nothing; o=Dict())
    atype = typeof(AutoGrad.getval(w["wsoft"]))

    # set dropouts
    vembdrop = get(o, :vembdrop, 0.0)
    wembdrop = get(o, :wembdrop, 0.0)
    softdrop = get(o, :softdrop, 0.0)
    fc7drop  = get(o, :fc7drop, 0.0)

    # visual features
    vis = dropout(vis, fc7drop)
    vemb = w["vemb"] * vis
    vemb = dropout(vemb, vembdrop)

    wemb = w["wemb"][:,x]
    y1, h1, c1 = rnnforw(srnn, w["rnn"], vemb; hy=true, cy=true)
    yy, hy = rnnforw(srnn, w["rnn"], wemb, h1, c1; batchSizes=batchsizes, hy=true)
    y2 = reshape(yy,size(yy,1),size(yy,2)*size(yy,3))
    ypred = w["wsoft"] * dropout(y2, softdrop) .+ w["bsoft"]
    # @show size(ypred), size(y)
    nll(ypred, y; average=false), length(y)
end

# generate
function generate(w, s, vis, vocab; maxlen=20, beamsize=1)
    atype = typeof(AutoGrad.getval(w["wdec"])) <: KnetArray ? KnetArray : Array
    wcnn  = get(w, "wcnn", nothing)
    vis = convert(atype, vis)
    if wcnn != nothing
        vis = vgg19(wcnn, vis)
        vis = transpose(vis)
    end

    x = vis * w["vemb"]
    (s[1], s[2]) = lstm(w["wdec"], w["bdec"], s[1], s[2], x)

    # language generation with (sentence, state, probability) array
    sentences = Any[(Any[SOS],s,0.0)]
    while true
        changed = false
        for i = 1:beamsize
            # get current sentence
            curr = shift!(sentences)
            sentence, st, prob = curr

            # get last word
            word = sentence[end]
            if word == EOS || length(sentence) >= maxlen
                push!(sentences, curr)
                continue
            end

            # get probabilities
            x = w["wemb"][word2index(vocab,word),:]
            x = reshape(x, 1, length(x))
            (st[1], st[2]) = lstm(w["wdec"], w["bdec"], st[1], st[2], x)
            ypred = logp(st[1] * w["wsoft"] .+ w["bsoft"], 2)
            ypred = convert(Array{Float32}, ypred)[:]

            # add most probable predictions to array
            maxinds = sortperm(ypred, rev=true)
            for j = 1:beamsize
                ind = maxinds[j]
                new_word = index2word(vocab, ind)
                new_sentence = copy(sentence)
                new_state = copy(st)
                new_probability = prob + ypred[ind]
                push!(new_sentence, new_word)
                push!(sentences, (new_sentence, new_state, new_probability))
            end
            changed = true

            # skip first loop
            if word == SOS
                break
            end
        end

        orders = sortperm(map(x -> x[3], sentences), rev=true)
        sentences = sentences[orders[1:beamsize]]

        if !changed
            break
        end
    end

    sentence = first(sentences)[1]
    if sentence[end] == EOS
        pop!(sentence)
    end
    push!(sentence, ".")
    output = join(filter(word -> word != UNK, sentence[2:end]), " ")
end

max_stars_repo_name: ilkerkesen/kusrp

max_issues_repo_name: ilkerkesen/kusrp

max_forks_repo_name: ilkerkesen/kusrp


text:
using MacroTools: postwalk

macro pair_with_wasm(xs)
  postwalk(xs) do x
    !(x isa Expr) ? :($x, $(Symbol(x, :_wasm))) : x
  end
end

macro wast_str(wast)
  return wast |> WebAssembly.parsebrackets |> WebAssembly.parse
end

relu_ifelse(x) = ifelse(x > 0, x, 0)
relu_ifelse_wasm = wast"""
(func $#relu_Int64  (param i64) (result i64)
  (i64.const 0)
  (get_local 0)
  (get_local 0)
  (i64.const 0)
  (i64.lt_s)
  (select)
  (return))
"""

relu_ternary(x) = x > 0 ? x : 0
relu_ternary_wasm = wast"""
(func $#relu_ternary_Int64  (param i64) (result i64)
  (i64.const 0)
  (get_local 0)
  (i64.lt_s)
  (if
    (then
      (get_local 0)
      (return)))
  (i64.const 0)
  (return))
"""

relu_if_then_else = relu_ifelse
relu_if_then_else_wasm = wast"""
(func $#relu_if_then_else_Int64  (param i64) (result i64)
  (i64.const 0)
  (get_local 0)
  (i64.lt_s)
  (if
    (then
      (get_local 0))
    (else
      (i64.const 0)))
  (return))
"""

pow = function (x, n)
  r = 1
  while n > 0
    r *= x
    n -= 1
  end
  return r
end
pow_wasm = wast"""
(func $#pow_Int64_Int64  (param i64) (param i64) (result i64)
  (local i64) (local i64)
  (get_local 1)
  (set_local 3)
  (i64.const 1)
  (set_local 2)
  (block
    (loop
      (i64.const 0)
      (get_local 3)
      (i64.lt_s)
      (i32.eqz)
      (br_if 1)
      (get_local 2)
      (get_local 0)
      (i64.mul)
      (set_local 2)
      (get_local 3)
      (i64.const 1)
      (i64.sub)
      (set_local 3)
      (br 0)))
  (get_local 2)
  (return))
"""

addTwo(x, y) = x + y
addTwo_wasm = wast"""
(func $addTwo (param i32) (param i32) (result i32)
  (get_local 0)
  (get_local 1)
  (i32.add))
"""

tests = @pair_with_wasm [ relu_ifelse
                        , relu_ternary
                        , relu_if_then_else
                        , pow
                        , addTwo
                        ]

fib(x) = x <= 1 ? 1 : fib(x - 1) + fib(x - 2)
this(x) = pow(x + 1, x - 1)

function rand_test_wasm(f, wasm_f, n_tests = 50, max = 100)
 for i in 1:n_tests
   args = [rand(WebAssembly.jltype(typ)) % max for typ in wasm_f.params]
   WebAssembly.interpretwasm(wasm_f, Dict(), args)[1] != f(args...) && return false
 end
 return true
end

function rand_test_module(fs, m, n_tests = 50, max = 10)
  wasm_fs = interpret_module(m)
  for i in eachindex(fs)
    for j in 1:n_tests
      args = [rand(WebAssembly.jltype(typ)) % max for typ in m.funcs[i].params]
      wasm_fs[i](args...)[1] != fs[i](args...) && return false
    end
  end
  return true
end

@testset "Parse-Interpret" begin

relu_wasm = relu_ifelse_wasm

relu_wasm_expected = Func(Symbol("#relu_Int64"), [i64], [i64], [], Block([Const(0), Local(0), Local(0), Const(0), Op(i64, :lt_s), Select(), Return()]))
@test relu_wasm.body.body == relu_wasm_expected.body.body
@test relu_wasm.params == relu_wasm_expected.params
@test relu_wasm.result == relu_wasm_expected.result
@test relu_wasm.name == relu_wasm_expected.name

@test rand_test_wasm(relu_ifelse, relu_wasm)

# Bulk form for functions
root = "test/wast/functions/"

for test in tests
  @test rand_test_wasm(test[1], test[2])
end

# Sort of test module parsing
m = wast"""
(module
  (func $addTwo (param i32) (param i32) (result i32)
    (get_local 0)
    (get_local 1)
    (i32.add))
  (export "addTwo" (func $addTwo)))
"""
@test m.exports == [Export(:addTwo, :addTwo, :func)]
expected_func = Func(Symbol("addTwo"), [i32, i32], [i32], [], Block([Local(0), Local(1), Op(i32, :add)]))
@test m.funcs[1].body.body == expected_func.body.body
@test m.funcs[1].params == expected_func.params
@test m.funcs[1].result == expected_func.result
@test m.funcs[1].name == expected_func.name

m2 = wast"""
(module
  (export "this" (func $#this_Int64))
  (export "pow" (func $#pow_Int64_Int64))
  (export "fib" (func $#fib_Int64))
  (memory 1)
  (func $#fib_Int64 (param i64) (result i64)
    (get_local 0)
    (i64.const 1)
    (i64.le_s)
    (if
      (then
        (i64.const 1)
        (return)))
    (get_local 0)
    (i64.const 1)
    (i64.sub)
    (call $#fib_Int64)
    (get_local 0)
    (i64.const 2)
    (i64.sub)
    (call $#fib_Int64)
    (i64.add)
    (return))
  (func $#this_Int64 (param i64) (result i64)
    (get_local 0)
    (i64.const 1)
    (i64.add)
    (get_local 0)
    (i64.const 1)
    (i64.sub)
    (call $#pow_Int64_Int64)
    (return))
  (func $#pow_Int64_Int64 (param i64) (param i64) (result i64)
  (local i64) (local i64)
    (get_local 1)
    (set_local 3)
    (i64.const 1)
    (set_local 2)
    (block
      (loop
        (i64.const 0)
        (get_local 3)
        (i64.lt_s)
        (i32.eqz)
        (br_if 1)
        (get_local 2)
        (get_local 0)
        (i64.mul)
        (set_local 2)
        (get_local 3)
        (i64.const 1)
        (i64.sub)
        (set_local 3)
        (br 0)))
    (get_local 2)
    (return)))
"""
@test m2.exports == [Export(:this, Symbol("#this_Int64"), :func), Export(:pow, Symbol("#pow_Int64_Int64"), :func), Export(:fib, Symbol("#fib_Int64"), :func)]
@test rand_test_module([fib, this, pow], m2)

end