text:
using sem, Feather, ModelingToolkit, Statistics, LinearAlgebra,
    Optim, SparseArrays, Test, Zygote, LineSearches, ForwardDiff

## Observed Data #############################################################

growth_dat = Feather.read("test/comparisons/growth_dat.feather")
growth_dat_miss30 = Feather.read("test/comparisons/growth_dat_miss30.feather")
#growth_par = Feather.read("test/comparisons/growth_par.feather")

definition_dat = Feather.read("test/comparisons/definition_dat.feather")
definition_dat_unique =
    Feather.read("test/comparisons/definition_dat_unique.feather")
definition_par = Feather.read("test/comparisons/definition_par.feather")
definition_par_missing =
    Feather.read("test/comparisons/definition_par_missing.feather")
definition_par_missing_unique =
    Feather.read("test/comparisons/definition_par_missing_unique.feather")

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

semobserved = SemObsCommon(data = Matrix(growth_dat); meanstructure = true, rowwise = true)

diff_fin = SemFiniteDiff(BFGS(), Optim.Options())

## Model definition ##########################################################

@ModelingToolkit.variables x[1:7], m[1:2], load_t[1:4]

S =[x[1]  0     0     0     0     0
    0     x[2]  0     0     0     0
    0     0     x[3]  0     0     0
    0     0     0     x[4]  0     0
    0     0     0     0     x[5]  x[7]
    0     0     0     0     x[7]  x[6]]

F =[1.0 0 0 0 0 0
    0 1 0 0 0 0
    0 0 1 0 0 0
    0 0 0 1 0 0]

A =[0  0  0  0  1.0  load_t[1]
    0  0  0  0  1  load_t[2]
    0  0  0  0  1  load_t[3]
    0  0  0  0  1  load_t[4]
    0  0  0  0  0  0
    0  0  0  0  0  0]

M = [0.0, 0, 0, 0, m[1:2]...]

S = sparse(S)

#F
F = sparse(F)

#A
A = sparse(A)

data_def = Matrix(definition_dat)

start_val = [1.0, 1, 1, 1, 0.05, 0.05, 0.0, 1.0, 1.0]

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

imply = ImplySymbolicDefinition(
    A,
    S,
    F,
    M,
    [x[1:7]..., m[1:2]...],
    load_t,
    start_val,
    data_def
)

loss = Loss([SemML(semobserved, imply, 0.0, 0.0)])

model_fin = Sem(semobserved, imply, loss, diff_fin)

solution = sem_fit(model_fin)

par_order = [collect(1:5); 7; 6; 8; 9]

all(
    abs.(solution.minimizer .- definition_par.Estimate[par_order]
        ) .< 0.05*abs.(definition_par.Estimate[par_order]))

## missings ############################################################

data_miss = Matrix(growth_dat_miss30)

semobserved = SemObsMissing(data_miss)

data_def = data_def[sem.remove_all_missing(data_miss)[2], :]

imply = ImplySymbolicDefinition(
    A,
    S,
    F,
    M,
    [x[1:7]..., m[1:2]...],
    load_t,
    start_val,
    data_def
    )

loss = Loss([SemFIML(semobserved, imply, 0.0, 0.0)])

model_fin = Sem(semobserved, imply, loss, diff_fin)

solution = sem_fit(model_fin)

par_order = [collect(1:5); 7; 6; 8; 9]

all(
    abs.(solution.minimizer .- definition_par_missing.Estimate[par_order]
        ) .< 0.05*abs.(definition_par_missing.Estimate[par_order]))

## missings with unique defvar patterns #####################################

data_def_unique = Matrix(definition_dat_unique)
data_def_unique = data_def_unique[sem.remove_all_missing(data_miss)[2], :]

imply = ImplySymbolicDefinition(
    A,
    S,
    F,
    M,
    [x[1:7]..., m[1:2]...],
    load_t,
    start_val,
    data_def_unique
    )

loss = Loss([SemFIML(semobserved, imply, 0.0, 0.0)])

model_fin = Sem(semobserved, imply, loss, diff_fin)

solution = sem_fit(model_fin)

par_order = [collect(1:5); 7; 6; 8; 9]

all(
    abs.(solution.minimizer .- definition_par_missing_unique.Estimate[par_order]
    ) .< 0.05*abs.(definition_par_missing_unique.Estimate[par_order])
)


text:
@test GeneralizedMetropolisHastings._valuestuple(1,0) == (1,0)
@test GeneralizedMetropolisHastings._valuestuple(1,1) == (1,)
@test GeneralizedMetropolisHastings._valuestuple(1,2) == (1,2)
@test GeneralizedMetropolisHastings._valuestuple(2,0) == (2,0)
@test GeneralizedMetropolisHastings._valuestuple(2,1) == (2,)
@test GeneralizedMetropolisHastings._valuestuple(2,2) == (2,2)

@test ndims(zeros(GeneralizedMetropolisHastings._valuestuple(1,0))) == 2
@test ndims(zeros(GeneralizedMetropolisHastings._valuestuple(1,1))) == 1
@test ndims(zeros(GeneralizedMetropolisHastings._valuestuple(1,2))) == 2
@test ndims(zeros(GeneralizedMetropolisHastings._valuestuple(2,0))) == 2
@test ndims(zeros(GeneralizedMetropolisHastings._valuestuple(2,1))) == 1
@test ndims(zeros(GeneralizedMetropolisHastings._valuestuple(2,2))) == 2

@test GeneralizedMetropolisHastings._tensortuple(1,1,0) == (1,1,0)
@test GeneralizedMetropolisHastings._tensortuple(1,1,1) == (1,1)
@test GeneralizedMetropolisHastings._tensortuple(1,1,2) == (1,1,2)
@test GeneralizedMetropolisHastings._tensortuple(2,1,0) == (2,1,0)
@test GeneralizedMetropolisHastings._tensortuple(2,1,1) == (2,1)
@test GeneralizedMetropolisHastings._tensortuple(2,1,2) == (2,1,2)

@test ndims(zeros(GeneralizedMetropolisHastings._tensortuple(1,1,0))) == 3
@test ndims(zeros(GeneralizedMetropolisHastings._tensortuple(1,1,1))) == 2
@test ndims(zeros(GeneralizedMetropolisHastings._tensortuple(1,1,2))) == 3
@test ndims(zeros(GeneralizedMetropolisHastings._tensortuple(1,2,0))) == 3
@test ndims(zeros(GeneralizedMetropolisHastings._tensortuple(1,2,1))) == 2
@test ndims(zeros(GeneralizedMetropolisHastings._tensortuple(1,2,2))) == 3

v1 = [1.0 2.0 3.0;4.0 5.0 6.0]
ll1 = [1.0,2.0,3.0]
lp1 = [0.1,0.2,0.3]
gll1 = 2*v1
glp1 = 3*v1
tll1 = similar(eye(2),2,2,3) ; tll1[:,:,1] = eye(2) ; tll1[:,:,2] = eye(2) ; tll1[:,:,3] = eye(2) ; tll1
tlp1 = similar(eye(2),2,2,3) ; tlp1[:,:,1] = ones(2,2) ; tlp1[:,:,2] = ones(2,2) ; tlp1[:,:,3] = ones(2,2) ; tlp1
tt1 = similar(eye(2,3),2,4,3) ; tt1[:,:,1] = ones(2,4) ; tt1[:,:,2] = 2*ones(2,4) ; tt1[:,:,3] = 3*ones(2,4) ; tt1

b1 = samples(:base,2,3,Float64,Float64)
b2 = samples(:base,2,1,Int,Float64)
b3 = samples(:base,2,3,Float32,Float32)
b4 = samples(:base,2,1,Int16,Float32)

@test typeof(b1) <: BaseSample && eltype(b1.values) == Float64 && size(b1.values) == (2,3) && size(b1.loglikelihood) == (3,) && size(b1.logprior) == (3,)
@test typeof(b2) <: BaseSample && eltype(b2.values) == Int && size(b2.values) == (2,) && size(b2.loglikelihood) == (1,) && size(b2.logprior) == (1,)
@test typeof(b3) <: BaseSample && eltype(b3.values) == Float32 && size(b3.values) == (2,3) && size(b3.loglikelihood) == (3,) && size(b3.logprior) == (3,)
@test typeof(b4) <: BaseSample && eltype(b4.values) == Int16 && size(b4.values) == (2,) && size(b4.loglikelihood) == (1,) && size(b4.logprior) == (1,)
@test ndims(b1.values) == 2 && ndims(b2.values) == 1

copy!(b1.values,v1)
copy!(b1.loglikelihood,ll1)
copy!(b1.logprior,lp1)
bc = samples(:base,2,3,Float64,Float64)
br = samples(:base,2,1,Float64,Float64)
@test b1 != bc
copy!(bc,b1)
@test b1 == bc
copy!(br,1,b1,3)
@test br.values == b1.values[:,3] && br.loglikelihood[1] == b1.loglikelihood[3] && br.logprior[1] == b1.logprior[3]
bo = samples(:base,2,1,Float64,Float64)
copy!(bo,1,br,1)
@test bo.values == br.values && bo.loglikelihood == br.loglikelihood && bo.logprior == br.logprior

bs1 = similar(b1)
bs2 = similar(b1,2)

@test isa(bs1,BaseSample) && numparas(bs1) == numparas(b1) && numsamples(bs1) == numsamples(b1) && sampletype(bs1) == sampletype(b1) && calculationtype(bs1) == calculationtype(b1)
@test isa(bs2,BaseSample) && numparas(bs2) == numparas(b1) && numsamples(bs2) == 2 && sampletype(bs2) == sampletype(b1) && calculationtype(bs2) == calculationtype(b1)

bc1 = copy(b1)
bc2 = copy(b1,1)
bc3 = copy(b1,[3,2])
bc4 = similar(b1) ; copy!(bc4,1,bc2,1) ; copy!(bc4,[2,3],bc3,[2,1])

@test bc1 == b1
@test bc4 == b1
@test numsamples(copy(b1,[])) == 0

offset!(bc1,[1.0,2.0])
@test bc1.values == (b1.values + [1.0 1.0 1.0;2.0 2.0 2.0])

###############################################################################################################################################################
g1 = samples(:gradient,2,3,Float64,Float64)
g2 = samples(:gradient,2,1,Float32,Float32)

@test typeof(g1) <: GradientSample && size(g1.values) == (2,3) && size(g1.loglikelihood) == (3,) && size(g1.logprior) == (3,) && size(g1.gradloglikelihood) == (2,3) && size(g1.gradlogprior) == (2,3)
@test typeof(g2) <: GradientSample && size(g2.values) == (2,) && size(g2.loglikelihood) == (1,) && size(g2.logprior) == (1,) && size(g2.gradloglikelihood) == (2,) && size(g2.gradlogprior) == (2,)
@test ndims(g1.values) == 2 && ndims(g2.values) == 1

copy!(g1.values,v1)
copy!(g1.loglikelihood,ll1)
copy!(g1.logprior,lp1)
copy!(g1.gradloglikelihood,gll1)
copy!(g1.gradlogprior,glp1)
gc = samples(:gradient,2,3,Float64,Float64)
gr = samples(:gradient,2,1,Float64,Float64)
@test g1 != gc
copy!(gc,g1)
@test g1 == gc
copy!(gr,1,g1,2)
@test gr.values == g1.values[:,2] && gr.loglikelihood[1] == g1.loglikelihood[2] && gr.logprior[1] == g1.logprior[2] && gr.gradloglikelihood == g1.gradloglikelihood[:,2] && gr.gradlogprior == g1.gradlogprior[:,2]

gs1 = similar(g1)
gs2 = similar(g1,2)
gs3 = similar(g1,2,:base)

@test isa(gs1,GradientSample) && numparas(gs1) == numparas(g1) && numsamples(gs1) == numsamples(g1) && sampletype(gs1) == sampletype(g1) && calculationtype(gs1) == calculationtype(g1)
@test isa(gs2,GradientSample) && numparas(gs2) == numparas(g1) && numsamples(gs2) == 2 && sampletype(gs2) == sampletype(g1) && calculationtype(gs2) == calculationtype(g1)
@test isa(gs3,BaseSample) && numparas(gs3) == numparas(g1) && numsamples(gs3) == 2 && sampletype(gs3) == sampletype(g1) && calculationtype(gs3) == calculationtype(g1)

gc1 = copy(g1)
gc2 = copy(g1,1)
gc3 = copy(g1,[3,2])
gc4 = similar(g1) ; copy!(gc4,1,gc2,1) ; copy!(gc4,[2,3],gc3,[2,1])
gc5 = copy(g1,[3,2],:base)

@test gc1 == g1
@test gc4 == g1
@test isa(gc5,BaseSample) && g1.values[:,[3,2]] == gc5.values
@test numsamples(copy(g1,[])) == 0

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

t1 = samples(:tensor,2,3,Float32,Float32)
t2 = samples(:tensor,2,1,Float64,Float64)

@test typeof(t1) <: TensorSample && size(t1.values) == (2,3) && size(t1.loglikelihood) == (3,) && size(t1.logprior) == (3,) && size(t1.gradloglikelihood) == (2,3) && size(t1.gradlogprior) == (2,3) && size(t1.tensorloglikelihood) == (2,2,3) && size(t1.tensorlogprior) == (2,2,3)
@test typeof(t2) <: TensorSample && size(t2.values) == (2,) && size(t2.loglikelihood) == (1,) && size(t2.logprior) == (1,) && size(t2.gradloglikelihood) == (2,) && size(t2.gradlogprior) == (2,) && size(t2.tensorloglikelihood) == (2,2) && size(t2.tensorlogprior) == (2,2)
@test ndims(t1.values) == 2 && ndims(t2.values) == 1

copy!(t1.values,v1)
copy!(t1.loglikelihood,ll1)
copy!(t1.logprior,lp1)
copy!(t1.gradloglikelihood,gll1)
copy!(t1.gradlogprior,glp1)
copy!(t1.tensorloglikelihood,tll1)
copy!(t1.tensorlogprior,tlp1)
tc = samples(:tensor,2,3,Float32,Float32)
tr = samples(:tensor,2,1,Float32,Float32)
@test t1 != tc
copy!(tc,t1)
@test t1 == tc
copy!(tr,1,t1,2)
@test tr.values == t1.values[:,2] && tr.loglikelihood[1] == t1.loglikelihood[2] && tr.logprior[1] == t1.logprior[2]
@test tr.gradloglikelihood == t1.gradloglikelihood[:,2] && tr.gradlogprior == t1.gradlogprior[:,2] && tr.tensorloglikelihood == t1.tensorloglikelihood[:,:,2] && tr.tensorlogprior == t1.tensorlogprior[:,:,2]

ts1 = similar(t1)
ts2 = similar(t1,2)
ts3 = similar(t1,2,:base)

@test isa(ts1,TensorSample) && numparas(ts1) == numparas(t1) && numsamples(ts1) == numsamples(t1) && sampletype(ts1) == sampletype(t1) && calculationtype(ts1) == calculationtype(t1)
@test isa(ts2,TensorSample) && numparas(ts2) == numparas(t1) && numsamples(ts2) == 2 && sampletype(ts2) == sampletype(t1) && calculationtype(ts2) == calculationtype(t1)
@test isa(ts3,BaseSample) && numparas(ts3) == numparas(t1) && numsamples(ts3) == 2 && sampletype(ts3) == sampletype(t1) && calculationtype(ts3) == calculationtype(t1)

@test numsamples(copy(t1,[])) == 0

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

a1 = samples(:tangent,2,3,Float64,Float64,4)
a2 = samples(:tangent,2,1,Float64,Float64,4)

@test typeof(a1) <: TangentTensorSample && size(a1.values) == (2,3) && size(a1.loglikelihood) == (3,) && size(a1.logprior) == (3,) && size(a1.gradloglikelihood) == (2,3) && size(a1.gradlogprior) == (2,3) && size(a1.tensorloglikelihood) == (2,2,3) && size(a1.tensorlogprior) == (2,2,3) && size(a1.tangentvectors) == (2,4,3)
@test typeof(a2) <: TangentTensorSample && size(a2.values) == (2,) && size(a2.loglikelihood) == (1,) && size(a2.logprior) == (1,) && size(a2.gradloglikelihood) == (2,) && size(a2.gradlogprior) == (2,) && size(a2.tensorloglikelihood) == (2,2) && size(a2.tensorlogprior) == (2,2) && size(a2.tangentvectors) == (2,4)
@test ndims(a1.values) == 2 && ndims(a2.values) == 1

copy!(a1.values,v1)
copy!(a1.loglikelihood,ll1)
copy!(a1.logprior,lp1)
copy!(a1.gradloglikelihood,gll1)
copy!(a1.gradlogprior,glp1)
copy!(a1.tensorloglikelihood,tll1)
copy!(a1.tensorlogprior,tlp1)
copy!(a1.tangentvectors,tt1)
ac = samples(:tangent,2,3,Float64,Float64,4)
ar = samples(:tangent,2,1,Float64,Float64,4)
@test a1 != ac
copy!(ac,a1)
@test a1 == ac
copy!(ar,1,a1,3)
@test ar.values == a1.values[:,3] && ar.loglikelihood[1] == a1.loglikelihood[3] && ar.logprior[1] == a1.logprior[3] && ar.tangentvectors == a1.tangentvectors[:,:,3]
@test ar.gradloglikelihood == a1.gradloglikelihood[:,3] && ar.gradlogprior == a1.gradlogprior[:,3] && ar.tensorloglikelihood == a1.tensorloglikelihood[:,:,3] && ar.tensorlogprior == a1.tensorlogprior[:,:,3]

as1 = similar(a1)
as2 = similar(a1,2)
as3 = similar(a1,2,:base)

@test isa(as1,TangentTensorSample) && numparas(as1) == numparas(a1) && numsamples(as1) == numsamples(a1) && sampletype(as1) == sampletype(a1) && calculationtype(as1) == calculationtype(a1)
@test isa(as2,TangentTensorSample) && numparas(as2) == numparas(a1) && numsamples(as2) == 2 && sampletype(as2) == sampletype(a1) && calculationtype(as2) == calculationtype(a1)
@test isa(as3,BaseSample) && numparas(as3) == numparas(a1) && numsamples(as3) == 2 && sampletype(as3) == sampletype(a1) && calculationtype(as3) == calculationtype(a1)

@test numparas(b1) == 2 && numsamples(b1) == 3
@test numparas(b2) == 2 && numsamples(b2) == 1
@test numparas(g1) == 2 && numsamples(g1) == 3
@test numparas(t1) == 2 && numsamples(t1) == 3
@test numparas(a1) == 2 && numsamples(a1) == 3 && numtangents(a1) == 4

@test numsamples(copy(a1,[])) == 0

println()
println("====================")
println("Test show() function")
show(b1)
show(g1)
show(t1)
show(a1)
show(b2)
show(g2)
show(t2)
show(a2)
println("End  show() function")
println("====================")
println()

nothing

ext: jl


text:
using Stations
s = Stations.load_station_info("espc2.info")
Stations.load_station_data(s, "espc2_1.csv")

t = Stations.obs_times(s)
o = Stations.observations(s, "RH")
o1 = o[t[1]]

println(o1)


text:
# Define the new domain structure
struct MyNewDomain <: InfiniteOpt.InfiniteScalarDomain
    attr1::Float64 # REPLACE WITH DESIRED TYPE
    attr2::Float64 # REPLACE WITH DESIRED TYPE
    # ADD MORE ATTRIBUTES AS NEEDED
    # constructor
    function MyNewDomain(attr1, attr2)
        # PUT APPROPRIATE CHECKS HERE AND OTHER NEEDED OPERATIONS
        return new(convert(Float64, attr1), convert(Float64, attr2))
    end
end

# Extend supports_in_domain
function InfiniteOpt.supports_in_domain(
    supports::Union{Real, Vector{<:Real}},
    domain::MyNewDomain
    )::Bool
    # DETERMINE IF SUPPORTS ARE IN THE DOMAIN OF `domain`
    in_domain = all(domain.attr1 .<= supports .<= domain.attr2) # REPLACE WITH ACTUAL CHECK
    return in_domain
end

# Extend generate_support_values if possible
function InfiniteOpt.generate_support_values(
    domain::MyNewDomain;
    num_supports::Int = 10,
    sig_digits::Int = 5
    )::Tuple{Vector{Float64}, DataType}
    # REPLACE BELOW WITH METHODS TO GENERATE `num_samples` with `sig_fig`
    supports = collect(range(domain.attr1, stop = domain.attr2, length = num_supports))
    return round.(supports, sigdigits = sig_digits), UniformGrid
end

# Extend JuMP.has_lower_bound (optional if the answer is always false)
function JuMP.has_lower_bound(domain::MyNewDomain)::Bool
    # INSERT NECESSARY CHECKS IF NEEDED
    has_bound = true # REPLACE WITH ACTUAL RESULT
    return has_bound
end

# Extend JuMP.lower_bound if appropriate
function JuMP.lower_bound(domain::MyNewDomain)
    return domain.attr1 # REPLACE WITH ACTUAL
end

# Extend JuMP.set_lower_bound if appropriate
function JuMP.set_lower_bound(domain::MyNewDomain, lower::Real)::MyNewDomain
    return MyNewDomain(lower, domain.attr2) # REPLACE WITH ACTUAL CONSTRUCTOR
end

# Extend JuMP.has_upper_bound (optional if the answer is always false)
function JuMP.has_upper_bound(domain::MyNewDomain)::Bool
    # INSERT NECESSARY CHECKS IF NEEDED
    has_bound = true # REPLACE WITH ACTUAL RESULT
    return has_bound
end

# Extend JuMP.lower_bound if appropriate
function JuMP.upper_bound(domain::MyNewDomain)
    return domain.attr2 # REPLACE WITH ACTUAL
end

# Extend JuMP.set_upper_bound if appropriate
function JuMP.set_upper_bound(domain::MyNewDomain, upper::Real)::MyNewDomain
    return MyNewDomain(domain.attr1, upper) # REPLACE WITH ACTUAL CONSTRUCTOR
end

# Extend InfiniteOpt.MeasureToolbox.generate_expect_data if wanted
function InfiniteOpt.MeasureToolbox.generate_expect_data(
    domain::MyNewDomain,
    pref::GeneralVariableRef,
    num_supports::Int;
    kwargs...
    )
    new_domain = IntervalDomain(domain.attr1, domain.attr2)
    return generate_expect_data(new_domain, pref, num_supports; kwargs...) # REPLACE WITH ACTUAL
end


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

using Markdown
using InteractiveUtils

# ╔═╡ e077ab16-22e5-11eb-15ed-5f0654bf7328
using Pkg, DrWatson

# ╔═╡ e75a9a06-22e5-11eb-2e36-8d4f62b830ed
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using RDatasets, StatsPlots
end

# ╔═╡ ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
md"## Listing4.26"

# ╔═╡ fa68607e-22e5-11eb-0558-c9a4d9f77426
iris = dataset("datasets", "iris")

# ╔═╡ 8894d350-2511-11eb-39d6-110ab1ec30db
@df iris andrewsplot(:Species, cols(1:4),
	line=(fill=[:blue :red :green]), legend=:topleft)

# ╔═╡ 475ff888-22e6-11eb-2354-09f8f40a8e12
md"## End of listing4.26"

# ╔═╡ Cell order:
# ╟─ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
# ╠═e077ab16-22e5-11eb-15ed-5f0654bf7328
# ╠═e75a9a06-22e5-11eb-2e36-8d4f62b830ed
# ╠═fa68607e-22e5-11eb-0558-c9a4d9f77426
# ╠═8894d350-2511-11eb-39d6-110ab1ec30db
# ╟─475ff888-22e6-11eb-2354-09f8f40a8e12


text:
simple_model = @stochastic_model begin
    @stage 1 begin
        @variable(model, x₁ >= 40)
        @variable(model, x₂ >= 20)
        @objective(model, Min, 100*x₁ + 150*x₂)
        @constraint(model, x₁ + x₂ <= 120)
    end
    @stage 2 begin
        @decision x₁ x₂
        @uncertain q₁ q₂ d₁ d₂
        @variable(model, 0 <= y₁ <= d₁)
        @variable(model, 0 <= y₂ <= d₂)
        @objective(model, Min, q₁*y₁ + q₂*y₂)
        @constraint(model, 6*y₁ + 10*y₂ <= 60*x₁)
        @constraint(model, 8*y₁ + 5*y₂ <= 80*x₂)
    end
end

s₁ = Scenario(q₁ = -24.0, q₂ = -28.0, d₁ = 500.0, d₂ = 100.0, probability = 0.4)
s₂ = Scenario(q₁ = -28.0, q₂ = -32.0, d₁ = 300.0, d₂ = 300.0, probability = 0.6)
simple = instantiate(simple_model, [s₁,s₂])

push!(problems, (simple, "Simple"))


text:
"""
set_constraints!(sdp_pb::SDP_Problem)

Set attributes `obj_keys`, `name_to_ctr`, `id_to_ctr` for a  `sdp_pb`
`SDP_Problem` already set with matrix linear and const contributions.
"""
function set_constraints!(sdp_pb::SDP_Problem)
    # Collecting constraint names
    ctr_names = SortedSet{SDP_CtrObjName}(collect(keys(sdp_pb.cst_ctr)))
    union!(ctr_names, [k[1] for k in keys(sdp_pb.matrices)])
    union!(ctr_names, [k[1] for k in keys(sdp_pb.linear)])

    # Extracting objective keys
    obj_keys = SortedSet{SDP_CtrObjName}()
    for ctr_name in ctr_names
        if (ctr_name[1:2] == ("1", "1"))
            push!(obj_keys, ctr_name)
            delete!(ctr_names, ctr_name)
        end
    end

    length(obj_keys) == 0 && error("No ctrkey matching objective key (\"1\", \"1\", .)")
    sdp_pb.obj_keys = obj_keys

    # Filling constraints attribute
    ctr_id = 1
    for ctr_name in ctr_names
        sdp_pb.name_to_ctr[ctr_name] = (ctr_id, "EQ", 0, 0)
        sdp_pb.id_to_ctr[ctr_id] = ctr_name
        ctr_id += 1
    end
end

"""
set_blocks!(sdp_pb::SDP_Problem)

Set string to id map for each sdp block given a `SDP_Problem` already set with
SDP matrix contributions.
"""
function set_blocks!(sdp_pb::SDP_Problem)
    blocktovarsset = SortedDict{String, SortedSet{String}}()

    for (ctr_name, blockname, var1, var2) in keys(sdp_pb.matrices)
        # sanity check
        @assert haskey(sdp_pb.name_to_sdpblock, blockname)

        !haskey(blocktovarsset, blockname) && (blocktovarsset[blockname]=SortedSet{String}())
        union!(blocktovarsset[blockname], [var1, var2])
    end

    for (blockname, vars) in blocktovarsset
        n = 0
        cur_block = sdp_pb.name_to_sdpblock[blockname]
        for var in vars
            n += 1
            cur_block.var_to_id[var] = n
        end
    end
end

"""
set_scalvars!(sdp_pb::SDP_Problem)

Set string to id map for each scalar variable given a `SDP_Problem` already set
with scalar linear contributions.
"""
function set_scalvars!(sdp_pb::SDP_Problem)
    for (ctr_name, var) in keys(sdp_pb.linear)

        if !haskey(sdp_pb.scalvar_to_id, var)
            sdp_pb.scalvar_to_id[var] = length(sdp_pb.scalvar_to_id) + 1
        end
    end
end


text:
export diff

import Base.diff

abstract Patch

immutable Insert <: Patch b end
immutable Overwrite{T <: Node} <: Patch b::T end
immutable Delete <: Patch end
immutable Reorder <: Patch moves::Dict end
immutable DictDiff <: Patch updates::Dict end

function key_idxs(ns)
    i = j = 1
    idxs = Dict()
    keys = Any[]
    for x in ns
        k = key(x)
        if k == nothing
            idxs[j] = i
            push!(keys, j)
            j += 1
        else
            idxs[k] = i
            push!(keys, k)
        end
        i += 1
    end
    idxs, keys
end

function diff!(a::AbstractVector, b::AbstractVector, index, patches; parentpatch=Patch[])
    if a === b return parentpatch end

    len_a = length(a)
    len_b = length(b)

    if len_a == 1 && len_b == 1
        diff!(a[1], b[1], index+1, patches)
        return parentpatch
    end

    a_key_idxs, a_keys = key_idxs(a)
    b_key_idxs, b_keys = key_idxs(b)

    ord = Dict()
    n = max(len_a, len_b)

    for i = 1:n
        index += 1
        if i <= len_a
            moved_to = get(b_key_idxs, a_keys[i], 0)
            if moved_to == 0
                patches[index] = Delete()
            else
                diff!(a[i], b[moved_to], index, patches)
                if i != moved_to
                    ord[i] = moved_to
                end
            end
            index += count(a[i])
        end
        if i <= len_b
            moved_from = get(a_key_idxs, b_keys[i], 0)
            if moved_from == 0
                # We don't need to save where this goes
                parentpatch = push!(parentpatch, Insert(b[i]))
            end
        end
    end
    if !isempty(ord)
        push!(parentpatch, Reorder(ord))
    end
end

diff!(a::TextNode, b::TextNode, index, patches) =
    a === b || a.text == b.text ?
        nothing :
        patches[index] = Patch[Overwrite(b)]

diff!(a::Elem, b::TextNode, index, patches) =
    patches[index] = Patch[Overwrite(b)]

diff!(a::TextNode, b::Elem, index, patches) =
    patches[index] = Patch[Overwrite(b)]

# If element's tag or ns is changed
diff!(a::Elem, b::Elem, index, patches) =
    patches[index] = Patch[Overwrite(b)]

function diff!{ns, tag}(a::Elem{ns, tag}, b::Elem{ns, tag}, index, patches)
    if a === b return patches end

    patch = get(patches, index, Patch[])

    proppatch = diff(properties(a), properties(b))
    if !is(proppatch, nothing)
        patch = push!(patch, DictDiff(proppatch))
    end

    diff!(a.children, b.children, index, patches, parentpatch=patch)
    if !isempty(patch)
        patches[index] = patch
    end
end

are_equal(a::AbstractArray, b::AbstractArray) = a === b || a == b
are_equal(a::AbstractString, b::Symbol) = a == string(b)
are_equal(a::Symbol, b::AbstractString) = string(a) == b
are_equal(a, b) = a == b

function diff(a::Associative, b::Associative)
    if a === b return nothing end

    updates = Dict()
    for (k, v) in a
        if haskey(b, k)
            if !are_equal(v, b[k])
                if isa(b[k], Associative) || isa(b[k], PropHook)
                    dictdiff = diff(v, b[k])
                    if !is(dictdiff, nothing)
                        updates[k] = dictdiff
                    end
                else
                    updates[k] = b[k]
                end
            end
        else
            # deleted
            if isa(v, Associative)
                updates[k] = diff(v, Dict())
            elseif isa(v, PropHook)
                updates[k] = diff(v, PropHook("", nothing))
            else
                updates[k] = nothing
            end
        end
    end
    for (k, v) in b
        if haskey(a, k) continue end
        updates[k] = v
    end
    if isempty(updates)
        return nothing
    else
        return updates
    end
end

function diff(a::Node, b::Node)
    patches = Dict()
    diff!(a, b, 0, patches)
    return patches
end


text:
# https://github.com/microsoft/Microsoft-MPI/blob/master/src/include/mpi.h

const MPI_Aint   = Int
const MPI_Offset = Int64
const MPI_Count  = Int64

const MPI_Comm = Cint
const MPI_COMM_NULL  = reinterpret(Cint, 0x04000000)
const MPI_COMM_SELF  = reinterpret(Cint, 0x44000001)
const MPI_COMM_WORLD = reinterpret(Cint, 0x44000000)

const MPI_Info = Cint
const MPI_INFO_NULL = reinterpret(Cint, 0x1c000000)

const MPI_Win = Cint
const MPI_WIN_NULL = reinterpret(Cint, 0x20000000)

const MPI_Op = Cint
const MPI_OP_NULL = reinterpret(Cint, 0x18000000)
const MPI_MAX     = reinterpret(Cint, 0x58000001)
const MPI_MIN     = reinterpret(Cint, 0x58000002)
const MPI_SUM     = reinterpret(Cint, 0x58000003)
const MPI_PROD    = reinterpret(Cint, 0x58000004)
const MPI_LAND    = reinterpret(Cint, 0x58000005)
const MPI_BAND    = reinterpret(Cint, 0x58000006)
const MPI_LOR     = reinterpret(Cint, 0x58000007)
const MPI_BOR     = reinterpret(Cint, 0x58000008)
const MPI_LXOR    = reinterpret(Cint, 0x58000009)
const MPI_BXOR    = reinterpret(Cint, 0x5800000a)
const MPI_MINLOC  = reinterpret(Cint, 0x5800000b)
const MPI_MAXLOC  = reinterpret(Cint, 0x5800000c)
const MPI_REPLACE = reinterpret(Cint, 0x5800000d)
const MPI_NO_OP   = reinterpret(Cint, 0x5800000e)

const MPI_Errhandler = Cint
const MPI_ERRORS_ARE_FATAL = Cint(0x54000000)
const MPI_ERRORS_RETURN = Cint(0x54000001)

const MPI_Request = Cint
const MPI_REQUEST_NULL = reinterpret(Cint, 0x2c000000)

const MPI_Datatype = Cint
const MPI_DATATYPE_NULL         = reinterpret(Cint, 0x0c000000)
const MPI_PACKED                = reinterpret(Cint, 0x4c00010f)
const MPI_CHAR                  = reinterpret(Cint, 0x4c000101)
const MPI_SIGNED_CHAR           = reinterpret(Cint, 0x4c000118)
const MPI_UNSIGNED_CHAR         = reinterpret(Cint, 0x4c000102)
const MPI_WCHAR                 = reinterpret(Cint, 0x4c00020e)
const MPI_BYTE                  = reinterpret(Cint, 0x4c00010d)
const MPI_SHORT                 = reinterpret(Cint, 0x4c000203)
const MPI_UNSIGNED_SHORT        = reinterpret(Cint, 0x4c000204)
const MPI_INT                   = reinterpret(Cint, 0x4c000405)
const MPI_UNSIGNED              = reinterpret(Cint, 0x4c000406)
const MPI_LONG                  = reinterpret(Cint, 0x4c000407)
const MPI_UNSIGNED_LONG         = reinterpret(Cint, 0x4c000408)
const MPI_FLOAT                 = reinterpret(Cint, 0x4c00040a)
const MPI_DOUBLE                = reinterpret(Cint, 0x4c00080b)

const MPI_INT8_T                = reinterpret(Cint, 0x4c000133)
const MPI_INT16_T               = reinterpret(Cint, 0x4c000234)
const MPI_INT32_T               = reinterpret(Cint, 0x4c000435)
const MPI_INT64_T               = reinterpret(Cint, 0x4c000836)
const MPI_UINT8_T               = reinterpret(Cint, 0x4c000137)
const MPI_UINT16_T              = reinterpret(Cint, 0x4c000238)
const MPI_UINT32_T              = reinterpret(Cint, 0x4c000439)
const MPI_UINT64_T              = reinterpret(Cint, 0x4c00083a)
const MPI_C_FLOAT_COMPLEX       = reinterpret(Cint, 0x4c000813)
const MPI_C_DOUBLE_COMPLEX      = reinterpret(Cint, 0x4c001014)

const MPI_File = Cint
const MPI_FILE_NULL = Cint(0)

const MPI_THREAD_SINGLE = Cint(0)
const MPI_THREAD_FUNNELED = Cint(1)
const MPI_THREAD_SERIALIZED = Cint(2)
const MPI_THREAD_MULTIPLE = Cint(3)
const MPI_PROC_NULL = Cint(-1)
const MPI_ANY_SOURCE = Cint(-2)
const MPI_ANY_TAG = Cint(-1)
const MPI_SUCCESS = Cint(0)
const MPI_UNDEFINED = Cint(-32766)
const MPI_LOCK_EXCLUSIVE = Cint(234)
const MPI_LOCK_SHARED = Cint(235)
const MPI_MAX_INFO_KEY = Cint(255)
const MPI_MAX_INFO_VAL = Cint(1024)
const MPI_MAX_OBJECT_NAME = Cint(128)
const MPI_TAG_UB = reinterpret(Cint, 0x64400001)
const MPI_COMM_TYPE_SHARED = Cint(1)
const MPI_ORDER_C = Cint(56)
const MPI_ORDER_FORTRAN = Cint(57)
const MPI_UNIVERSE_SIZE = reinterpret(Cint, 0x64400009)
const MPI_MAX_ERROR_STRING = Cint(512)
const MPI_SUCCESS = Cint(0)
const MPI_MODE_RDONLY = reinterpret(Cint, 0x00000002)
const MPI_MODE_RDWR = reinterpret(Cint, 0x00000008)
const MPI_MODE_WRONLY = reinterpret(Cint, 0x00000004)
const MPI_MODE_CREATE = reinterpret(Cint, 0x00000001)
const MPI_MODE_EXCL = reinterpret(Cint, 0x00000040)
const MPI_MODE_DELETE_ON_CLOSE = reinterpret(Cint, 0x00000010)
const MPI_MODE_UNIQUE_OPEN = reinterpret(Cint, 0x00000020)
const MPI_MODE_SEQUENTIAL = reinterpret(Cint, 0x00000100)
const MPI_MODE_APPEND = reinterpret(Cint, 0x00000080)
const MPI_MODE_NOCHECK    = Cint(1024)
const MPI_MODE_NOSTORE    = Cint(2048)
const MPI_MODE_NOPUT      = Cint(4096)
const MPI_MODE_NOPRECEDE  = Cint(8192)
const MPI_MODE_NOSUCCEED  = Cint(16384)

const MPI_SEEK_SET = Cint(600)
const MPI_SEEK_CUR = Cint(602)
const MPI_SEEK_END = Cint(604)

const MPI_IDENT = Cint(0)
const MPI_CONGRUENT = Cint(1)
const MPI_SIMILAR = Cint(2)
const MPI_UNEQUAL = Cint(3)

const MPI_BOTTOM = reinterpret(SentinelPtr, 0)
const MPI_IN_PLACE = reinterpret(SentinelPtr, -1)

const MPI_STATUS_IGNORE = reinterpret(SentinelPtr, 1)
const MPI_STATUSES_IGNORE = reinterpret(SentinelPtr, 1)

struct Status
    _pad1::Cint
    _pad2::Cint
    source::Cint
    tag::Cint
    error::Cint
end

const STATUS_EMPTY = Status(zero(Cint), zero(Cint), MPI_ANY_SOURCE, MPI_ANY_TAG, MPI_SUCCESS)


text:
using Conda

Conda.add("ase", channel = "conda-forge")
Conda.add("rdkit", channel = "conda-forge")
Conda.add("pymatgen", channel = "conda-forge")
if "mkl" in Conda._installed_packages()
    Conda.rm("mkl")
    Conda.add("mkl")
end


text:
using RasterIO
using Base.Test

function check_create(drivername::String, copy::Int=0)
    driver = RasterIO.driverbyname(drivername)
    if copy == 1
        createcopy = RasterIO.getmetadataitem(driver,"DCAP_CREATECOPY")
        check = createcopy == C_NULL ? false : true
    else
        create = RasterIO.getmetadataitem(driver,"DCAP_CREATE")
        check = create == C_NULL ? false : true
    end
    return check
end

# Check to make sure the drivers are properly registered
# GeoTiff and EHdr are built into GDAL, so this should work on all
# systems
@test check_create("GTiff",1)
@test check_create("GTiff")
@test check_create("EHdr",1)
@test check_create("EHdr")

# Open a tiny sample raster file and check that its data matches
# expectations
A = RasterIO.openraster("data/A.tif")
@test A.width == 3
@test A.height == 3
@test RasterIO.geotransform(A) == zeros(6)
@test RasterIO.wktprojection(A) == ""

# Write opened raster to file.
RasterIO.writeraster(A,"A","EHdr")

# Open the just-written raster and check that it matches the first
A2 = RasterIO.openraster("A")

@test A.nband == A2.nband
@test A.width == A2.width
@test A.height == A2.height
@test RasterIO.geotransform(A) == RasterIO.geotransform(A2)

RasterIO.closeraster(A2)

#Remove the written raster files so as not to clutter things up

for file in readdir()
    if file == "A" || file == "A.hdr"
        rm(file)
    end
end


text:
# # @__NAME__

# PREAMBLE

# PKG_SETUP

# ## Setup

using Pkg # hide
Pkg.activate(@__DIR__) # hide
Pkg.instantiate() # hide
using RigidBodyDynamics, StaticArrays, ForwardDiff
using Test, Random
Random.seed!(1); # to get repeatable results

# ## Jacobians with respect to $q$ and $v$ - the naive way

# First, we'll load our trusty double pendulum from a URDF:

srcdir = dirname(pathof(RigidBodyDynamics))
urdf = joinpath(srcdir, "..", "test", "urdf", "Acrobot.urdf")
mechanism = parse_urdf(urdf)

# Of course, we can create a `MechanismState` for the double pendulum, and compute its momentum in some random state:

float64state = MechanismState(mechanism)
rand!(float64state)
momentum(float64state)

# But now suppose we want the Jacobian of momentum with respect to the joint velocity vector $v$.
# We can do this using the `ForwardDiff.Dual` type and the `ForwardDiff.jacobian` function.
# The ForwardDiff package implements forward-mode [automatic differentiation](https://en.wikipedia.org/wiki/Automatic_differentiation).

# To use `ForwardDiff.jacobian` we'll create a function that maps `v` (as a `Vector`) to momentum (as a `Vector`):

q = configuration(float64state)
function momentum_vec(v::AbstractVector{T}) where T
    ## create a `MechanismState` that can handle the element type of `v` (which will be some `ForwardDiff.Dual`):
    state = MechanismState{T}(mechanism)

    ## set the state variables:
    set_configuration!(state, q)
    set_velocity!(state, v)

    ## return momentum converted to an `SVector` (as ForwardDiff expects an `AbstractVector`)
    Vector(SVector(momentum(state)))
end

# Let's first check that the function returns the same thing we got from `float64state`:

v = velocity(float64state)
@test momentum_vec(v) == SVector(momentum(float64state))

# That works, so now let's compute the Jacobian with `ForwardDiff`:

J = ForwardDiff.jacobian(momentum_vec, v)

# At this point we note that the matrix `J` is simply the momentum matrix (in world frame) of the `Mechanism`. In this case, RigidBodyDynamics.jl has a specialized algorithm for computing this matrix, so let's verify the results:

A = momentum_matrix(float64state)
@test J ≈ Array(A) atol = 1e-12

# Gradients with respect to $q$ can be computed in similar fashion.

# ## Improving performance

# Ignoring the fact that we have a specialized method available, let's look at the performance of using `ForwardDiff.jacobian`.

using BenchmarkTools
@benchmark ForwardDiff.jacobian($momentum_vec, $v)

# That's not great. Note all the allocations. We can do better by making the following modifications:
#
# 1. use an in-place version of the `jacobian` function, `ForwardDiff.jacobian!`
# 2. reimplement our `momentum_vec` function to be in-place as well
# 3. don't create a new `MechanismState` every time
#
# The third point is especially important; creating a `MechanismState` is expensive!
#
# Regarding the second point, we could also just stop converting momentum from a `StaticArrays.SVector` to a `Vector` to avoid allocations. However, the solution of making the function in-place also applies when the size of the output vector is not known statically (e.g., for `dynamics_bias!`).

# To facillitate reuse of `MechanismState`s while keeping the code nice and generic, we can use a `StateCache` object.
# `StateCache` is a container that stores `MechanismState`s of various types (associated with one `Mechanism`), and will ease the process of using `ForwardDiff`.
# Creating one is easy:

const statecache = StateCache(mechanism)

# `MechanismState`s of a given type can be accessed as follows (note that if a `MechanismState` of a certain type is already available, it will be reused):

float32state = statecache[Float32]
@test float32state === statecache[Float32]

# Now we'll use the `StateCache` to reimplement `momentum_vec`, making it in-place as well:

function momentum_vec!(out::AbstractVector, v::AbstractVector{T}) where T
    ## retrieve a `MechanismState` that can handle the element type of `v`:
    state = statecache[T]

    ## set the state variables:
    set_configuration!(state, q)
    set_velocity!(state, v)

    ## compute momentum and store it in `out`
    m = momentum(state)
    copyto!(out, SVector(m))
end

# Check that the in-place version works as expected on `Float64` inputs:

const out = zeros(6) # where we'll be storing our results
momentum_vec!(out, v)
@test out == SVector(momentum(float64state))

# And use `ForwardDiff.jacobian!` to compute the Jacobian:

const result = DiffResults.JacobianResult(out, v)
const config = ForwardDiff.JacobianConfig(momentum_vec!, out, v)
ForwardDiff.jacobian!(result, momentum_vec!, out, v, config)
J = DiffResults.jacobian(result)
@test J ≈ Array(A) atol = 1e-12

# Let's check the performance again:

@benchmark ForwardDiff.jacobian!($result, $momentum_vec!, $out, $v, $config)

# That's much better. Do note that the specialized algorithm is still faster:

q = copy(configuration(float64state))
@benchmark begin
    set_configuration!($float64state, $q)
    momentum_matrix!($A, $float64state)
end

# ## Time derivatives

# We can also use ForwardDiff to compute time derivatives. Let's verify that energy is conserved for the double pendulum in the absence of nonconservative forces (like damping). That is, we expect that the time derivative of the pendulum's total energy is zero when its state evolves according to the passive dynamics.

# Let's first compute the joint acceleration vector $\dot{v}$ using the passive dynamics:

dynamicsresult = DynamicsResult(mechanism)
set_configuration!(float64state, q)
set_velocity!(float64state, v)
dynamics!(dynamicsresult, float64state)
v̇ = dynamicsresult.v̇

# Now for the time derivative of total energy. ForwardDiff has a `derivative` function that can be used to take derivatives of functions that map a scalar to a scalar. But in this example, we'll instead use ForwardDiff's `Dual` type directly. `ForwardDiff.Dual` represents a (potentially multidimensional) dual number, i.e., a type that stores both the value of a function evaluated at a certain point, as well as the partial derivatives of the function, again evaluated at the same point. See the [ForwardDiff documentation](http://www.juliadiff.org/ForwardDiff.jl/stable/dev/how_it_works.html) for more information.

# We'll create a vector of `Dual`s representing the value and derivative of $q(t)$:

q̇ = v
q_dual = ForwardDiff.Dual.(q, q̇)

# **Note**: for the double pendulum, $\dot{q} = v$, but this is not the case in general for `Mechanism`s created using RigidBodyDynamics.jl. For example, the `QuaternionSpherical` joint type uses a unit quaternion to represent the joint configuration, but angular velocity (in body frame) to represent velocity. In general $\dot{q}$ can be computed from the velocity vector $v$ stored in a `MechanismState` using
#
# ```julia
# configuration_derivative(::MechanismState)
# ```
#
# or its in-place variant, `configuration_derivative!`.

# We'll do the same thing for $v(t)$:

v_dual = ForwardDiff.Dual.(v, v̇)

# Now we're ready to compute the total energy (kinetic + potential) using these `ForwardDiff.Dual` inputs. We'll use our `StateCache` again:

T = eltype(q_dual)
state = statecache[T]
set_configuration!(state, q_dual)
set_velocity!(state, v_dual)
energy_dual = kinetic_energy(state) + gravitational_potential_energy(state)

# Note that the result type of `energy_dual` is again a `ForwardDiff.Dual`. We can extract the energy and its time derivative (mechanical power) from `energy_dual` as follows:

energy = ForwardDiff.value(energy_dual)
partials = ForwardDiff.partials(energy_dual)
power = partials[1];

# So the total energy in the system is:

energy

# **Note**: the total energy is negative because the origin of the world frame is used as a reference for computing gravitational potential energy, i.e., the center of mass of the double pendulum is somewhere below this origin.

# And we can verify that, indeed, there is no power flow into or out of the system:

@test power ≈ 0 atol = 1e-14

max_stars_repo_path: examples/5. Derivatives and gradients using ForwardDiff/5. Derivatives and gradients using ForwardDiff.jl


text:
export isequal_sparse

function isequal_sparse(A_1::SparseMatrixCSC, A_2::SparseMatrixCSC)
    (   size(A_1)==size(A_2)
    && nonzeros(A_1)==nonzeros(A_2)
    && rowvals(A_1)==rowvals(A_2)
    && A_1.colptr==A_2.colptr)
end

if VERSION < v"0.4.0-dev"
    # TODO: remove these two for 0.4
    rowvals(S::SparseMatrixCSC) = S.rowval
    nzrange(S::SparseMatrixCSC, col::Integer) = S.colptr[col]:(S.colptr[col+1]-1)
end


text:
using OMEinsum
using Test

@testset "OMEinsum.jl" begin
    # Write your own tests here.
end


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

# Parse "GIT URLs" syntax (URLs and a scp-like syntax). For details see:
# https://git-scm.com/docs/git-clone#_git_urls_a_id_urls_a
const URL_REGEX = r"""
^(?:(?<scheme>ssh|git|https?)://)?+
(?:
    (?<user>.*?)
    (?:\:(?<password>.*?))?@
)?
(?<host>[A-Za-z0-9\-\.]+)
(?(<scheme>)
    # Only parse port when not using scp-like syntax
    (?:\:(?<port>\d+))?
    /?
    |
    :?
)
(?<path>
    # Require path to be preceeded by '/'. Alternatively, ':' when using scp-like syntax.
    (?<=(?(<scheme>)/|:))
    .*
)?
$
"""x

"""
    version() -> VersionNumber

Return the version of libgit2 in use, as a [`VersionNumber`](@ref man-version-number-literals).
"""
function version()
    major = Ref{Cint}(0)
    minor = Ref{Cint}(0)
    patch = Ref{Cint}(0)
    ccall((:git_libgit2_version, :libgit2), Cvoid,
          (Ref{Cint}, Ref{Cint}, Ref{Cint}), major, minor, patch)
    return VersionNumber(major[], minor[], patch[])
end
const VERSION = version()

"""
    isset(val::Integer, flag::Integer)

Test whether the bits of `val` indexed by `flag` are set (`1`) or unset (`0`).
"""
isset(val::Integer, flag::Integer) = (val & flag == flag)

"""
    reset(val::Integer, flag::Integer)

Unset the bits of `val` indexed by `flag`, returning them to `0`.
"""
reset(val::Integer, flag::Integer) = (val &= ~flag)

"""
    toggle(val::Integer, flag::Integer)

Flip the bits of `val` indexed by `flag`, so that if a bit is `0` it
will be `1` after the toggle, and vice-versa.
"""
toggle(val::Integer, flag::Integer) = (val |= flag)

"""
    features()

Return a list of git features the current version of libgit2 supports, such as
threading or using HTTPS or SSH.
"""
function features()
    feat = ccall((:git_libgit2_features, :libgit2), Cint, ())
    res = Consts.GIT_FEATURE[]
    for f in instances(Consts.GIT_FEATURE)
        isset(feat, Cuint(f)) && Base.push!(res, f)
    end
    return res
end

"""
    LibGit2.posixpath(path)

Standardise the path string `path` to use POSIX separators.
"""
function posixpath end
if Sys.iswindows()
    posixpath(path) = replace(path,'\\' => '/')
elseif Sys.isunix()
    posixpath(path) = path
end

"""
    LibGit2.git_url(; kwargs...) -> String

Create a string based upon the URL components provided. When the `scheme` keyword is not
provided the URL produced will use the alternative [scp-like syntax](https://git-scm.com/docs/git-clone#_git_urls_a_id_urls_a).

# Keywords

  * `scheme::AbstractString=""`: the URL scheme which identifies the protocol to be used.
    For HTTP use "http", SSH use "ssh", etc. When `scheme` is not provided the output format
    will be "ssh" but using the scp-like syntax.
  * `username::AbstractString=""`: the username to use in the output if provided.
  * `password::AbstractString=""`: the password to use in the output if provided.
  * `host::AbstractString=""`: the hostname to use in the output. A hostname is required to
    be specified.
  * `port::Union{AbstractString,Integer}=""`: the port number to use in the output if
    provided. Cannot be specified when using the scp-like syntax.
  * `path::AbstractString=""`: the path to use in the output if provided.

# Examples
```jldoctest
julia> LibGit2.git_url(username="git", host="github.com", path="JuliaLang/julia.git")
"git@github.com:JuliaLang/julia.git"

julia> LibGit2.git_url(scheme="https", host="github.com", path="/JuliaLang/julia.git")
"https://github.com/JuliaLang/julia.git"

julia> LibGit2.git_url(scheme="ssh", username="git", host="github.com", port=2222, path="JuliaLang/julia.git")
"ssh://git@github.com:2222/JuliaLang/julia.git"
```
"""
function git_url(;
        scheme::AbstractString="",
        username::AbstractString="",
        password::AbstractString="",
        host::AbstractString="",
        port::Union{AbstractString,Integer}="",
        path::AbstractString="")

    port_str = string(port)
    scp_syntax = isempty(scheme)

    isempty(host) && throw(ArgumentError("A host needs to be specified"))
    scp_syntax && !isempty(port_str) && throw(ArgumentError("Port cannot be specified when using scp-like syntax"))

    io = IOBuffer()
    !isempty(scheme) && print(io, scheme, "://")

    if !isempty(username) || !isempty(password)
        print(io, username)
        !isempty(password) && print(io, ':', password)
        print(io, '@')
    end

    print(io, host)
    !isempty(port_str) && print(io, ':', port_str)

    if !isempty(path)
        if scp_syntax
            print(io, ':')
        elseif !startswith(path, '/')
            print(io, '/')
        end
        print(io, path)
    end

    return String(take!(io))
end

function credential_identifier(scheme::AbstractString, host::AbstractString)
    string(isempty(scheme) ? "ssh" : scheme, "://", host)
end

function credential_identifier(url::AbstractString)
    m = match(URL_REGEX, url)
    scheme = something(m[:scheme], "")
    host = m[:host]
    credential_identifier(scheme, host)
end

max_stars_repo_name: perrutquist/julia

max_issues_repo_name: perrutquist/julia

max_forks_repo_name: perrutquist/julia


text:
# Unit testing for texpr

using Base.Test
using Devectorize

###########################################################
#
#   simple expression construction
#
###########################################################

# tnum

ex = texpr(1.5)
@test isa(ex, TNum{Float64})
@test ex.val == 1.5
@test ex == tnum(1.5)
@test eval(ju_expr(ex)) == 1.5

# tvar & tscalarvar

a = 12.5

ex = texpr(:a)
@test isa(ex, TVar)
@test ex.name == :a
@test ex == tvar(:a)
@test eval(ju_expr(ex)) == a

ex = tscalarvar(:a)
@test isa(ex, TScalarVar)
@test ex.name == :a
@test ex == tscalarvar(:a)
@test eval(ju_expr(ex)) == a

# tqvar

ex = texpr(:(a.b))
@test isa(ex, TQVar)
@test ex.form == :(a.b)
@test ex == tqvar(:(a.b))

type A1
    b
end
a = A1(2.3)
@test eval(ju_expr(ex)) == a.b

ex = texpr(:(a.b.c))
@test isa(ex, TQVar)
@test ex.form == :(a.b.c)
@test ex == tqvar(:(a.b.c))

type A2
    c
end
a = A1(A2(4.5))
@test eval(ju_expr(ex)) == a.b.c

###########################################################
#
#   scalar reference expression construction
#
###########################################################

a = [10 20 30; 40 50 60]

ex = texpr(:(a[1]))
@test isa(ex, TGeneralRef1)
@test ex.host == tvar(:a)
@test ex.i == 1
@test ex == tref(:(a[1]))
@test eval(ju_expr(ex)) == a[1]

ex = texpr(:(a[x]))
@test isa(ex, TGeneralRef1)
@test ex.host == tvar(:a)
@test ex.i == :x
@test ex == tref(:(a[x]))
x = 2
@test eval(ju_expr(ex)) == a[x]

ex = texpr(:(a[1,2]))
@test isa(ex, TGeneralRef2)
@test ex.host == tvar(:a)
@test ex.i == 1
@test ex.j == 2
@test ex == tref(:(a[1,2]))
@test eval(ju_expr(ex)) == a[1,2]

ex = texpr(:(a[x,2]))
@test isa(ex, TGeneralRef2)
@test ex.host == tvar(:a)
@test ex.i == :x
@test ex.j == 2
@test ex == tref(:(a[x,2]))
@test eval(ju_expr(ex)) == a[x,2]

ex = texpr(:(a[1,y]))
@test isa(ex, TGeneralRef2)
@test ex.host == tvar(:a)
@test ex.i == 1
@test ex.j == :y
@test ex == tref(:(a[1,y]))
y = 3
@test eval(ju_expr(ex)) == a[1,y]

ex = texpr(:(a[x,y]))
@test isa(ex, TGeneralRef2)
@test ex.host == tvar(:a)
@test ex.i == :x
@test ex.j == :y
@test ex == tref(:(a[x,y]))
@test eval(ju_expr(ex)) == a[x,y]

a = A1(a)

ex = texpr(:(a.b[2]))
@test isa(ex, TGeneralRef1)
@test ex.host == tqvar(:(a.b))
@test ex.i == 2
@test ex == tref(:(a.b[2]))
@test eval(ju_expr(ex)) == a.b[2]

ex = texpr(:(a.b[2,3]))
@test isa(ex, TGeneralRef2)
@test ex.host == tqvar(:(a.b))
@test ex.i == 2
@test ex.j == 3
@test ex == tref(:(a.b[2,3]))
@test eval(ju_expr(ex)) == a.b[2,3]

a = {"s" => 100}

ex = texpr(:(a["s"]))
@test isa(ex, TGeneralRef1)
@test ex.host == tvar(:a)
@test ex.i == "s"
@test ex == tref(:(a["s"]))
@test eval(ju_expr(ex)) == a["s"]

ex = texpr(:(a["s", "t"]))
@test isa(ex, TGeneralRef2)
@test ex.host == tvar(:a)
@test ex.i == "s"
@test ex.j == "t"
@test ex == tref(:(a["s", "t"]))

###########################################################
#
#   1D array reference expression construction
#
###########################################################

ex = texpr(:(a[:]))
@test isa(ex, TRef1D)
@test ex.host == tvar(:a)
@test ex.rgn == TColon()
@test ex == tref(:(a[:]))

ex = texpr(:(a[2:]))
@test isa(ex, TRef1D)
@test ex.host == tvar(:a)
@test ex.rgn == TInterval(2, nothing)
@test ex == tref(:(a[2:]))

ex = texpr(:(a[i0:]))
@test isa(ex, TRef1D)
@test ex.host == tvar(:a)
@test ex.rgn == TInterval(:i0, nothing)
@test ex == tref(:(a[i0:]))

ex = texpr(:(a[3:7]))
@test isa(ex, TRef1D)
@test ex.host == tvar(:a)
@test ex.rgn == TInterval(3, 7)
@test ex == tref(:(a[3:7]))

ex = texpr(:(a[i0:7]))
@test isa(ex, TRef1D)
@test ex.host == tvar(:a)
@test ex.rgn == TInterval(:i0, 7)
@test ex == tref(:(a[i0:7]))

ex = texpr(:(a[3:i1]))
@test isa(ex, TRef1D)
@test ex.host == tvar(:a)
@test ex.rgn == TInterval(3, :i1)
@test ex == tref(:(a[3:i1]))

ex = texpr(:(a[i0:i1]))
@test isa(ex, TRef1D)
@test ex.host == tvar(:a)
@test ex.rgn == TInterval(:i0, :i1)
@test ex == tref(:(a[i0:i1]))

ex = texpr(:(a.b[:]))
@test isa(ex, TRef1D)
@test ex.host == tqvar(:(a.b))
@test ex.rgn == TColon()
@test ex == tref(:(a.b[:]))

###########################################################
#
#   column reference expression construction
#
###########################################################

ex = texpr(:(a[:, 5]))
@test isa(ex, TRefCol)
@test ex.host == tvar(:a)
@test ex.rrgn == TColon()
@test ex.icol == 5
@test ex == tref(:(a[:, 5]))

ex = texpr(:(a[2:, 5]))
@test isa(ex, TRefCol)
@test ex.host == tvar(:a)
@test ex.rrgn == TInterval(2, nothing)
@test ex.icol == 5
@test ex == tref(:(a[2:, 5]))

ex = texpr(:(a[i0:, 5]))
@test isa(ex, TRefCol)
@test ex.host == tvar(:a)
@test ex.rrgn == TInterval(:i0, nothing)
@test ex.icol == 5
@test ex == tref(:(a[i0:, 5]))

ex = texpr(:(a[3:7, 5]))
@test isa(ex, TRefCol)
@test ex.host == tvar(:a)
@test ex.rrgn == TInterval(3, 7)
@test ex.icol == 5
@test ex == tref(:(a[3:7, 5]))

ex = texpr(:(a[i0:7, 5]))
@test isa(ex, TRefCol)
@test ex.host == tvar(:a)
@test ex.rrgn == TInterval(:i0, 7)
@test ex.icol == 5
@test ex == tref(:(a[i0:7, 5]))

ex = texpr(:(a[3:i1, 5]))
@test isa(ex, TRefCol)
@test ex.host == tvar(:a)
@test ex.rrgn == TInterval(3, :i1)
@test ex.icol == 5
@test ex == tref(:(a[3:i1, 5]))

ex = texpr(:(a[i0:i1, 5]))
@test isa(ex, TRefCol)
@test ex.host == tvar(:a)
@test ex.rrgn == TInterval(:i0, :i1)
@test ex.icol == 5
@test ex == tref(:(a[i0:i1, 5]))

ex = texpr(:(a[:, j]))
@test isa(ex, TRefCol)
@test ex.host == tvar(:a)
@test ex.rrgn == TColon()
@test ex.icol == :j
@test ex == tref(:(a[:, j]))

ex = texpr(:(a[2:, j]))
@test isa(ex, TRefCol)
@test ex.host == tvar(:a)
@test ex.rrgn == TInterval(2, nothing)
@test ex.icol == :j
@test ex == tref(:(a[2:, j]))

ex = texpr(:(a[i0:, j]))
@test isa(ex, TRefCol)
@test ex.host == tvar(:a)
@test ex.rrgn == TInterval(:i0, nothing)
@test ex.icol == :j
@test ex == tref(:(a[i0:, j]))

ex = texpr(:(a[3:7, j]))
@test isa(ex, TRefCol)
@test ex.host == tvar(:a)
@test ex.rrgn == TInterval(3, 7)
@test ex.icol == :j
@test ex == tref(:(a[3:7, j]))

ex = texpr(:(a[i0:7, j]))
@test isa(ex, TRefCol)
@test ex.host == tvar(:a)
@test ex.rrgn == TInterval(:i0, 7)
@test ex.icol == :j
@test ex == tref(:(a[i0:7, j]))

ex = texpr(:(a[3:i1, j]))
@test isa(ex, TRefCol)
@test ex.host == tvar(:a)
@test ex.rrgn == TInterval(3, :i1)
@test ex.icol == :j
@test ex == tref(:(a[3:i1, j]))

ex = texpr(:(a[i0:i1, j]))
@test isa(ex, TRefCol)
@test ex.host == tvar(:a)
@test ex.rrgn == TInterval(:i0, :i1)
@test ex.icol == :j
@test ex == tref(:(a[i0:i1, j]))

ex = texpr(:(a.b[:,j]))
@test isa(ex, TRefCol)
@test ex.host == tqvar(:(a.b))
@test ex.rrgn == TColon()
@test ex.icol == :j
@test ex == tref(:(a.b[:,j]))

###########################################################
#
#   row reference expression construction
#
###########################################################

ex = texpr(:(a[6, :]))
@test isa(ex, TRefRow)
@test ex.host == tvar(:a)
@test ex.irow == 6
@test ex.crgn == TColon()
@test ex == tref(:(a[6, :]))

ex = texpr(:(a[6, 2:]))
@test isa(ex, TRefRow)
@test ex.host == tvar(:a)
@test ex.irow == 6
@test ex.crgn == TInterval(2, nothing)
@test ex == tref(:(a[6, 2:]))

ex = texpr(:(a[6, i0:]))
@test isa(ex, TRefRow)
@test ex.host == tvar(:a)
@test ex.irow == 6
@test ex.crgn == TInterval(:i0, nothing)
@test ex == tref(:(a[6, i0:]))

ex = texpr(:(a[6, 3:7]))
@test isa(ex, TRefRow)
@test ex.host == tvar(:a)
@test ex.irow == 6
@test ex.crgn == TInterval(3, 7)
@test ex == tref(:(a[6, 3:7]))

ex = texpr(:(a[6, i0:7]))
@test isa(ex, TRefRow)
@test ex.host == tvar(:a)
@test ex.irow == 6
@test ex.crgn == TInterval(:i0, 7)
@test ex == tref(:(a[6, i0:7]))

ex = texpr(:(a[6, 3:i1]))
@test isa(ex, TRefRow)
@test ex.host == tvar(:a)
@test ex.irow == 6
@test ex.crgn == TInterval(3, :i1)
@test ex == tref(:(a[6, 3:i1]))

ex = texpr(:(a[6, i0:i1]))
@test isa(ex, TRefRow)
@test ex.host == tvar(:a)
@test ex.irow == 6
@test ex.crgn == TInterval(:i0, :i1)
@test ex == tref(:(a[6, i0:i1]))

ex = texpr(:(a[i, :]))
@test isa(ex, TRefRow)
@test ex.host == tvar(:a)
@test ex.irow == :i
@test ex.crgn == TColon()
@test ex == tref(:(a[i, :]))

ex = texpr(:(a[i, 2:]))
@test isa(ex, TRefRow)
@test ex.host == tvar(:a)
@test ex.irow == :i
@test ex.crgn == TInterval(2, nothing)
@test ex == tref(:(a[i, 2:]))

ex = texpr(:(a[i, i0:]))
@test isa(ex, TRefRow)
@test ex.host == tvar(:a)
@test ex.irow == :i
@test ex.crgn == TInterval(:i0, nothing)
@test ex == tref(:(a[i, i0:]))

ex = texpr(:(a[i, 3:7]))
@test isa(ex, TRefRow)
@test ex.host == tvar(:a)
@test ex.irow == :i
@test ex.crgn == TInterval(3, 7)
@test ex == tref(:(a[i, 3:7]))

ex = texpr(:(a[i, i0:7]))
@test isa(ex, TRefRow)
@test ex.host == tvar(:a)
@test ex.irow == :i
@test ex.crgn == TInterval(:i0, 7)
@test ex == tref(:(a[i, i0:7]))

ex = texpr(:(a[i, 3:i1]))
@test isa(ex, TRefRow)
@test ex.host == tvar(:a)
@test ex.irow == :i
@test ex.crgn == TInterval(3, :i1)
@test ex == tref(:(a[i, 3:i1]))

ex = texpr(:(a[i, i0:i1]))
@test isa(ex, TRefRow)
@test ex.host == tvar(:a)
@test ex.irow == :i
@test ex.crgn == TInterval(:i0, :i1)
@test ex == tref(:(a[i, i0:i1]))

ex = texpr(:(a.b[i, :]))
@test isa(ex, TRefRow)
@test ex.host == tqvar(:(a.b))
@test ex.irow == :i
@test ex.crgn == TColon()
@test ex == tref(:(a.b[i, :]))

###########################################################
#
#   2D reference expression construction
#
###########################################################

ex = texpr(:(a[:, :]))
@test isa(ex, TRef2D)
@test ex.host == tvar(:a)
@test ex.rrgn == TColon()
@test ex.crgn == TColon()
@test ex == tref(:(a[:, :]))

ex = texpr(:(a[:, 2:]))
@test isa(ex, TRef2D)
@test ex.host == tvar(:a)
@test ex.rrgn == TColon()
@test ex.crgn == TInterval(2, nothing)
@test ex == tref(:(a[:, 2:]))

ex = texpr(:(a[:, j0:]))
@test isa(ex, TRef2D)
@test ex.host == tvar(:a)
@test ex.rrgn == TColon()
@test ex.crgn == TInterval(:j0, nothing)
@test ex == tref(:(a[:, j0:]))

ex = texpr(:(a[:, 3:7]))
@test isa(ex, TRef2D)
@test ex.host == tvar(:a)
@test ex.rrgn == TColon()
@test ex.crgn == TInterval(3, 7)
@test ex == tref(:(a[:, 3:7]))

ex = texpr(:(a[:, j0:7]))
@test isa(ex, TRef2D)
@test ex.host == tvar(:a)
@test ex.rrgn == TColon()
@test ex.crgn == TInterval(:j0, 7)
@test ex == tref(:(a[:, j0:7]))

ex = texpr(:(a[:, 3:j1]))
@test isa(ex, TRef2D)
@test ex.host == tvar(:a)
@test ex.rrgn == TColon()
@test ex.crgn == TInterval(3, :j1)
@test ex == tref(:(a[:, 3:j1]))

ex = texpr(:(a[:, j0:j1]))
@test isa(ex, TRef2D)
@test ex.host == tvar(:a)
@test ex.rrgn == TColon()
@test ex.crgn == TInterval(:j0, :j1)
@test ex == tref(:(a[:, j0:j1]))

ex = texpr(:(a[i0:, :]))
@test isa(ex, TRef2D)
@test ex.host == tvar(:a)
@test ex.rrgn == TInterval(:i0, nothing)
@test ex.crgn == TColon()
@test ex == tref(:(a[i0:, :]))

ex = texpr(:(a[i0:, 2:]))
@test isa(ex, TRef2D)
@test ex.host == tvar(:a)
@test ex.rrgn == TInterval(:i0, nothing)
@test ex.crgn == TInterval(2, nothing)
@test ex == tref(:(a[i0:, 2:]))

ex = texpr(:(a[i0:, j0:]))
@test isa(ex, TRef2D)
@test ex.host == tvar(:a)
@test ex.rrgn == TInterval(:i0, nothing)
@test ex.crgn == TInterval(:j0, nothing)
@test ex == tref(:(a[i0:, j0:]))

ex = texpr(:(a[i0:, 3:7]))
@test isa(ex, TRef2D)
@test ex.host == tvar(:a)
@test ex.rrgn == TInterval(:i0, nothing)
@test ex.crgn == TInterval(3, 7)
@test ex == tref(:(a[i0:, 3:7]))

ex = texpr(:(a[i0:, j0:7]))
@test isa(ex, TRef2D)
@test ex.host == tvar(:a)
@test ex.rrgn == TInterval(:i0, nothing)
@test ex.crgn == TInterval(:j0, 7)
@test ex == tref(:(a[i0:, j0:7]))

ex = texpr(:(a[i0:, 3:j1]))
@test isa(ex, TRef2D)
@test ex.host == tvar(:a)
@test ex.rrgn == TInterval(:i0, nothing)
@test ex.crgn == TInterval(3, :j1)
@test ex == tref(:(a[i0:, 3:j1]))

ex = texpr(:(a[i0:, j0:j1]))
@test isa(ex, TRef2D)
@test ex.host == tvar(:a)
@test ex.rrgn == TInterval(:i0, nothing)
@test ex.crgn == TInterval(:j0, :j1)
@test ex == tref(:(a[i0:, j0:j1]))

ex = texpr(:(a[i0:i1, :]))
@test isa(ex, TRef2D)
@test ex.host == tvar(:a)
@test ex.rrgn == TInterval(:i0, :i1)
@test ex.crgn == TColon()
@test ex == tref(:(a[i0:i1, :]))

ex = texpr(:(a[i0:i1, 2:]))
@test isa(ex, TRef2D)
@test ex.host == tvar(:a)
@test ex.rrgn == TInterval(:i0, :i1)
@test ex.crgn == TInterval(2, nothing)
@test ex == tref(:(a[i0:i1, 2:]))

ex = texpr(:(a[i0:i1, j0:]))
@test isa(ex, TRef2D)
@test ex.host == tvar(:a)
@test ex.rrgn == TInterval(:i0, :i1)
@test ex.crgn == TInterval(:j0, nothing)
@test ex == tref(:(a[i0:i1, j0:]))

ex = texpr(:(a[i0:i1, 3:7]))
@test isa(ex, TRef2D)
@test ex.host == tvar(:a)
@test ex.rrgn == TInterval(:i0, :i1)
@test ex.crgn == TInterval(3, 7)
@test ex == tref(:(a[i0:i1, 3:7]))

ex = texpr(:(a[i0:i1, j0:7]))
@test isa(ex, TRef2D)
@test ex.host == tvar(:a)
@test ex.rrgn == TInterval(:i0, :i1)
@test ex.crgn == TInterval(:j0, 7)
@test ex == tref(:(a[i0:i1, j0:7]))

ex = texpr(:(a[i0:i1, 3:j1]))
@test isa(ex, TRef2D)
@test ex.host == tvar(:a)
@test ex.rrgn == TInterval(:i0, :i1)
@test ex.crgn == TInterval(3, :j1)
@test ex == tref(:(a[i0:i1, 3:j1]))

ex = texpr(:(a[i0:i1, j0:j1]))
@test isa(ex, TRef2D)
@test ex.host == tvar(:a)
@test ex.rrgn == TInterval(:i0, :i1)
@test ex.crgn == TInterval(:j0, :j1)
@test ex == tref(:(a[i0:i1, j0:j1]))

ex = texpr(:(a.b[i0:i1, j0:j1]))
@test isa(ex, TRef2D)
@test ex.host == tqvar(:(a.b))
@test ex.rrgn == TInterval(:i0, :i1)
@test ex.crgn == TInterval(:j0, :j1)
@test ex == tref(:(a.b[i0:i1,j0:j1]))

###########################################################
#
#   Fun call expression construction
#
###########################################################

# maps

ex = texpr(:(sin(a)))
@test isa(ex, TMap)
@test ex.fun == :sin
@test length(ex.args) == 1
@test ex.args[1] == tvar(:a)
@test ex == tcall(:(sin(a)))

ex = texpr(:(a + b))
@test isa(ex, TMap)
@test ex.fun == :(+)
@test length(ex.args) == 2
@test ex.args[1] == tvar(:a)
@test ex.args[2] == tvar(:b)
@test ex == tcall(:(a + b))

ex = texpr(:(clamp(a, b, c)))
@test isa(ex, TMap)
@test ex.fun == :clamp
@test length(ex.args) == 3
@test ex.args[1] == tvar(:a)
@test ex.args[2] == tvar(:b)
@test ex.args[3] == tvar(:c)
@test ex == tcall(:(clamp(a, b, c)))

ex = texpr(:(a + b .* c))
@test isa(ex, TMap)
@test ex.fun == :+
@test length(ex.args) == 2
@test ex.args[1] == tvar(:a)
@test ex.args[2] == tcall(:(b .* c))

ex = texpr(:(a.x + max(b[2], c[:,i])))
@test isa(ex, TMap)
@test ex.fun == :+
@test length(ex.args) == 2
@test ex.args[1] == tqvar(:(a.x))
@test ex.args[2].fun == :max
@test ex.args[2].args[1] == tref(:(b[2]))
@test ex.args[2].args[2] == tref(:(c[:,i]))

# reductions

ex = texpr(:(sum(a)))
@test isa(ex, TReduc)
@test ex.fun == :sum
@test length(ex.args) == 1
@test ex.args[1] == tvar(:a)

ex = texpr(:(max(a)))
@test isa(ex, TReduc)
@test ex.fun == :max
@test length(ex.args) == 1
@test ex.args[1] == tvar(:a)

ex = texpr(:(min(a)))
@test isa(ex, TReduc)
@test ex.fun == :min
@test length(ex.args) == 1
@test ex.args[1] == tvar(:a)

ex = texpr(:(mean(a)))
@test isa(ex, TReduc)
@test ex.fun == :mean
@test length(ex.args) == 1
@test ex.args[1] == tvar(:a)

ex = texpr(:(dot(a, b)))
@test isa(ex, TReduc)
@test ex.fun == :dot
@test length(ex.args) == 2
@test ex.args[1] == tvar(:a)
@test ex.args[2] == tvar(:b)

ex = texpr(:(sum(a, 1)))
@test isa(ex, TColwiseReduc)
@test ex.fun == :sum
@test length(ex.args) == 1
@test ex.args[1] == tvar(:a)

ex = texpr(:(max(a, (), 1)))
@test isa(ex, TColwiseReduc)
@test ex.fun == :max
@test length(ex.args) == 1
@test ex.args[1] == tvar(:a)

ex = texpr(:(min(a, (), 1)))
@test isa(ex, TColwiseReduc)
@test ex.fun == :min
@test length(ex.args) == 1
@test ex.args[1] == tvar(:a)

ex = texpr(:(mean(a, 1)))
@test isa(ex, TColwiseReduc)
@test ex.fun == :mean
@test length(ex.args) == 1
@test ex.args[1] == tvar(:a)

ex = texpr(:(sum(a, 2)))
@test isa(ex, TRowwiseReduc)
@test ex.fun == :sum
@test length(ex.args) == 1
@test ex.args[1] == tvar(:a)

ex = texpr(:(max(a, (), 2)))
@test isa(ex, TRowwiseReduc)
@test ex.fun == :max
@test length(ex.args) == 1
@test ex.args[1] == tvar(:a)

ex = texpr(:(min(a, (), 2)))
@test isa(ex, TRowwiseReduc)
@test ex.fun == :min
@test length(ex.args) == 1
@test ex.args[1] == tvar(:a)

ex = texpr(:(mean(a, 2)))
@test isa(ex, TRowwiseReduc)
@test ex.fun == :mean
@test length(ex.args) == 1
@test ex.args[1] == tvar(:a)

ex = texpr(:(sum(a[1,2] + b .* c[:,:])))
@test isa(ex, TReduc)
@test ex.fun == :sum
@test length(ex.args) == 1
@test ex.args[1] == tcall(:(a[1,2] + b .* c[:,:]))

ex = texpr(:(sum(a[1,2] + b .* c[:,:], 1)))
@test isa(ex, TColwiseReduc)
@test ex.fun == :sum
@test length(ex.args) == 1
@test ex.args[1] == tcall(:(a[1,2] + b .* c[:,:]))

ex = texpr(:(sum(a[1,2] + b .* c[:,:], 2)))
@test isa(ex, TRowwiseReduc)
@test ex.fun == :sum
@test length(ex.args) == 1
@test ex.args[1] == tcall(:(a[1,2] + b .* c[:,:]))

###########################################################
#
#   Assignments and blocks
#
###########################################################

# assignment

ex = texpr(:(a = b))
@test isa(ex, TAssign)
@test ex.lhs == tvar(:a)
@test ex.rhs == tvar(:b)
@test is_trivial_assignment(ex)
@test ex == tassign(:(a = b))

ex = texpr(:(r.a = b))
@test isa(ex, TAssign)
@test ex.lhs == tqvar(:(r.a))
@test ex.rhs == tvar(:b)
@test is_trivial_assignment(ex)
@test ex == tassign(:(r.a = b))

ex = texpr(:(r.a = s.b))
@test isa(ex, TAssign)
@test ex.lhs == tqvar(:(r.a))
@test ex.rhs == tqvar(:(s.b))
@test is_trivial_assignment(ex)
@test ex == tassign(:(r.a = s.b))

ex = texpr(:(a = b + c))
@test isa(ex, TAssign)
@test ex.lhs == tvar(:a)
@test ex.rhs == tcall(:(b + c))
@test !is_trivial_assignment(ex)
@test ex == tassign(:(a = b + c))

ex = texpr(:(r.a[:,1] = sin(a) + b[:]))
@test isa(ex, TAssign)
@test ex.lhs == tref(:(r.a[:,1]))
@test ex.rhs == tcall(:(sin(a) + b[:]))
@test !is_trivial_assignment(ex)
@test ex == tassign(:(r.a[:,1] = sin(a) + b[:]))

# op-assignment

ex = texpr(:(r += a))
@test isa(ex, TAssign)
@test ex.lhs == tvar(:r)
@test ex.rhs == tcall(:(r + a))
@test ex == tassign(:(r = r + a))

ex = texpr(:(r.x -= a.y))
@test isa(ex, TAssign)
@test ex.lhs == tqvar(:(r.x))
@test ex.rhs == tcall(:(r.x - a.y))
@test ex == tassign(:(r.x = r.x - a.y))

ex = texpr(:(r .*= a + b))
@test isa(ex, TAssign)
@test ex.lhs == tvar(:r)
@test ex.rhs == tcall(:(r .* (a + b)))
@test ex == tassign(:(r = r .* (a + b)))

ex = texpr(:(r ./= a + b))
@test isa(ex, TAssign)
@test ex.lhs == tvar(:r)
@test ex.rhs == tcall(:(r ./ (a + b)))
@test ex == tassign(:(r = r ./ (a + b)))

# tblock

blk = quote
    a = b + 1
    c.x[:] = a[:,j]
    c .*= 2
end
ex = texpr(blk)
@test isa(ex, TBlock)
@test length(ex.stmts) == 3
@test ex.stmts[1] == tassign(:(a = b + 1))
@test ex.stmts[2] == tassign(:(c.x[:] = a[:,j]))
@test ex.stmts[3] == tassign(:(c = c .* 2))
@test ex == tblock(blk)


text:
using RSDSGE
using JLD
using Base.Test

# Run the tests of the model constructor
@testset "Model setup" begin
    include("modelsetup.jl")
end

# Check model solving works properly
@testset "Solving" begin
    include("solve.jl")
end

# Check model utils
@testset "Model utils" begin
    include("modelutils.jl")
end


text:
using Statistics
using Oceananigans.BuoyancyModels: g_Earth

using Oceananigans.Models.HydrostaticFreeSurfaceModels:
    ImplicitFreeSurface,
    FreeSurface,
    implicit_free_surface_step!,
    implicit_free_surface_linear_operation!

function run_implicit_free_surface_solver_tests(arch, grid)

    Δt = 900
    Nx = grid.Nx
    Ny = grid.Ny

    # Create a model
    model = HydrostaticFreeSurfaceModel(architecture = arch,
                                        grid = grid,
                                        momentum_advection = nothing,
                                        free_surface=ImplicitFreeSurface())

    # Create a divergent velocity
    u, v, w = model.velocities
    imid = Int(floor(grid.Nx / 2)) + 1
    jmid = Int(floor(grid.Ny / 2)) + 1
    CUDA.@allowscalar u[imid, jmid, 1] = 1

    implicit_free_surface_step!(model.free_surface, Event(device(arch)), model, Δt, 1.5)

    # Extract right hand side "truth"
    right_hand_side = model.free_surface.implicit_step_right_hand_side

    # Compute left hand side "solution"
    g = g_Earth
    η = model.free_surface.η
    ∫ᶻ_Axᶠᶜᶜ = model.free_surface.vertically_integrated_lateral_face_areas.xᶠᶜᶜ
    ∫ᶻ_Ayᶜᶠᶜ = model.free_surface.vertically_integrated_lateral_face_areas.yᶜᶠᶜ

    left_hand_side = ReducedField(Center, Center, Nothing, arch, grid; dims=3)
    implicit_free_surface_linear_operation!(left_hand_side, η, ∫ᶻ_Axᶠᶜᶜ, ∫ᶻ_Ayᶜᶠᶜ, g, Δt)

    # Compare
    extrema_tolerance = 1e-10
    std_tolerance = 1e-10

    CUDA.@allowscalar begin
        @test minimum(abs, interior(left_hand_side) .- interior(right_hand_side)) < extrema_tolerance
        @test maximum(abs, interior(left_hand_side) .- interior(right_hand_side)) < extrema_tolerance
        @test std(interior(left_hand_side) .- interior(right_hand_side)) < std_tolerance
    end

    return nothing
end

@testset "Implicit free surface solver tests" begin
    for arch in archs

        rectilinear_grid = RegularRectilinearGrid(size = (128, 1, 5),
                                                  x = (0, 1000kilometers), y = (0, 1), z = (-400, 0),
                                                  topology = (Bounded, Periodic, Bounded))

        lat_lon_grid = RegularLatitudeLongitudeGrid(size = (90, 90, 5),
                                                    longitude = (-30, 30),
                                                    latitude = (15, 75),
                                                    z = (-4000, 0))

        for grid in (rectilinear_grid, lat_lon_grid)
            @info "Testing implicit free surface solver [$(typeof(arch)), $(typeof(grid).name.wrapper)]..."
            run_implicit_free_surface_solver_tests(arch, grid)
        end
    end
end


text:
include("args.jl")
using DelimitedFiles

function sw_p_to_lambda_den(sw, p)
    sw = tf.reshape(sw, (1, m, n, 1))
    p = tf.reshape(p, (1, m, n, 1))
    sw = tf.image.resize_bilinear(sw, (nz, nx))
    p = tf.image.resize_bilinear(p, (nz, nx))
    sw = cast(sw, Float64)
    p = cast(p, Float64)
    sw = squeeze(sw)
    p = squeeze(p)
    # tran_lambda, tran_den = Gassman(sw)
    # tran_lambda, tran_den = RockLinear(sw) # test linear relationship
    tran_lambda, tran_den = Patchy(sw)
    return tran_lambda, tran_den
end

if !isdir("figures_summary")
  mkdir("figures_summary")
end

iter = 100
Prj_names = "CO2/flow_fit_results";
K_name = "/K$iter.txt"
K = readdlm(Prj_names*K_name)

tfCtxTrue = tfCtxGen(m,n,h,NT,Δt,Z,X,ρw,ρo,μw,μo,K,g,ϕ,qw,qo, sw0, true)
out_sw_true, out_p_true = imseq(tfCtxTrue)
lambdas = Array{PyObject}(undef, n_survey)
dens = Array{PyObject}(undef, n_survey)
for i = 1:n_survey
    sw = out_sw_true[survey_indices[i]]
    p = out_p_true[survey_indices[i]]
    lambdas[i], dens[i] = sw_p_to_lambda_den(sw, p)
end

sess = Session();init(sess);

vps = Array{PyObject}(undef, n_survey)
for i=1:n_survey
  vps[i] = sqrt((lambdas[i] + 2.0 * tf_shear_sat1[i])/dens[i])
end
V = run(sess, vps);
S = run(sess, out_sw_true);
P = run(sess, out_p_true);

z_inj = (9-1)*h + h/2.0
x_inj = (3-1)*h + h/2.0
z_prod = (9-1)*h + h/2.0
x_prod = (28-1)*h + h/2.0

rc("axes", titlesize=30)
rc("axes", labelsize=30)
rc("xtick", labelsize=28)
rc("ytick", labelsize=28)
rc("legend", fontsize=30)
fig1,axs = subplots(3,3, figsize=[30,15], sharex=true, sharey=true)
ims = Array{Any}(undef, 9)
for iPrj = 1:3
  for jPrj = 1:3
    ims[(iPrj-1)*3+jPrj] = axs[iPrj,jPrj].imshow(V[(iPrj-1)*3+jPrj], extent=[0,n*h,m*h,0], vmin=3350, vmax=3500);
    axs[iPrj,jPrj].title.set_text("Snapshot $((iPrj-1)*3+jPrj)")
    if jPrj == 1 || jPrj == 1
      axs[iPrj,jPrj].set_ylabel("Depth (m)")
    end
    if iPrj == 3 || iPrj == 3
      axs[iPrj,jPrj].set_xlabel("Distance (m)")
    end
    # cb = fig1.colorbar(ims[(iPrj-1)*3+jPrj], ax=axs[iPrj,jPrj])
    # cb.set_label("Vp")
    axs[iPrj,jPrj].scatter(x_inj, z_inj, c="r", marker=">", s=128)
    axs[iPrj,jPrj].scatter(x_prod, z_prod, c="r", marker="<", s=128)
  end
end
fig1.subplots_adjust(wspace=0.02, hspace=0.18)
cbar_ax = fig1.add_axes([0.91, 0.08, 0.01, 0.82])
cb1 = fig1.colorbar(ims[1], cax=cbar_ax)
cb1.set_label("Vp (m/s)")
savefig("figures_summary/predicted_Vp_evo_sep.pdf",bbox_inches="tight",pad_inches = 0);

fig2,axs = subplots(3,3, figsize=[30,15], sharex=true, sharey=true)
ims = Array{Any}(undef, 9)
for iPrj = 1:3
  for jPrj = 1:3
    ims[(iPrj-1)*3+jPrj] = axs[iPrj,jPrj].imshow(S[survey_indices[(iPrj-1)*3+jPrj], :, :], extent=[0,n*h,m*h,0], vmin=0.0, vmax=0.6);
    axs[iPrj,jPrj].title.set_text("Snapshot $((iPrj-1)*3+jPrj)")
    if jPrj == 1 || jPrj == 1
      axs[iPrj,jPrj].set_ylabel("Depth (m)")
    end
    if iPrj == 3 || iPrj == 3
      axs[iPrj,jPrj].set_xlabel("Distance (m)")
    end
    # if iPrj ==2 && jPrj == 3
    # cb = fig2.colorbar(ims[(iPrj-1)*3+jPrj], ax=axs[iPrj,jPrj])
    # cb.set_label("Saturation")
    axs[iPrj,jPrj].scatter(x_inj, z_inj, c="r", marker=">", s=128)
    axs[iPrj,jPrj].scatter(x_prod, z_prod, c="r", marker="<", s=128)
  end
end
# fig2.subplots_adjust(wspace=0.04, hspace=0.042)
fig2.subplots_adjust(wspace=0.02, hspace=0.18)
cbar_ax = fig2.add_axes([0.91, 0.08, 0.01, 0.82])
cb2 = fig2.colorbar(ims[1], cax=cbar_ax)
cb2.set_label("Saturation")
savefig("figures_summary/predicted_Saturation_evo_sep.pdf",bbox_inches="tight",pad_inches = 0);

# fig3,axs = subplots(3,3, figsize=[30,15], sharex=true, sharey=true)
# ims = Array{Any}(undef, 9)
# for iPrj = 1:3
#   for jPrj = 1:3
#     ims[(iPrj-1)*3+jPrj] = axs[iPrj,jPrj].imshow(P[survey_indices[(iPrj-1)*3+jPrj], :, :]*1.4504e-04, extent=[0,n*h,m*h,0], vmin=-2500.0, vmax=500);
#     axs[iPrj,jPrj].title.set_text("Snapshot $((iPrj-1)*3+jPrj)")
#     if jPrj == 1 || jPrj == 1
#       axs[iPrj,jPrj].set_ylabel("Depth (m)")
#     end
#     if iPrj == 3 || iPrj == 3
#       axs[iPrj,jPrj].set_xlabel("Distance (m)")
#     end
#     # if iPrj ==2 && jPrj == 3
#     # cb = fig2.colorbar(ims[(iPrj-1)*3+jPrj], ax=axs[iPrj,jPrj])
#     # cb.set_label("Saturation")
#     axs[iPrj,jPrj].scatter(x_inj, z_inj, c="r", marker=">")
#     axs[iPrj,jPrj].scatter(x_prod, z_prod, c="r", marker="<")
#   end
# end
# # fig2.subplots_adjust(wspace=0.04, hspace=0.042)
# fig3.subplots_adjust(wspace=0.02, hspace=0.18)
# cbar_ax = fig3.add_axes([0.91, 0.08, 0.01, 0.82])
# cb3 = fig3.colorbar(ims[1], cax=cbar_ax)
# cb3.set_label("Potential (psi)")
# savefig("figures_summary/Potential_evo_patchy_true.pdf",bbox_inches="tight",pad_inches = 0);

# iter = 100
# Prj_names = ["CO2", "CO2_1src", "CO2_2surveys", "CO2_6surveys"]
# K_name = "/K$iter.txt"

# fig,axs = subplots(2,2, figsize=[18,8], sharex=true, sharey=true)
# for iPrj = 1:2
#   for jPrj = 1:2
#     # println(ax)
#     A = readdlm(Prj_names[(iPrj-1)*2 + jPrj] * K_name)
#     im = axs[iPrj,jPrj].imshow(A, extent=[0,n*h,m*h,0]);
#     if jPrj == 1 || jPrj == 1
#       axs[iPrj,jPrj].set_ylabel("Depth (m)")
#     end
#     if iPrj == 2 || iPrj == 2
#       axs[iPrj,jPrj].set_xlabel("Distance (m)")
#     end
#     axs[iPrj,jPrj].text(-0.1,1.1,string("(" * Char((iPrj-1)*2 + jPrj+'a'-1) * ")"),transform=axs[iPrj,jPrj].transAxes,size=12,weight="bold")
#   end
# end
# fig.subplots_adjust(bottom=0.1, top=0.9, left=0.1, right=0.9,
#                     wspace=0.1, hspace=0.2)
# cb_ax = fig.add_axes([0.93, 0.1, 0.02, 0.8])
# cbar = fig.colorbar(im, cax=cb_ax)

# cb = fig.colorbar()
# clim([20, 120])
# cb.set_label("Permeability (md)")

# fig = figure()
# ax = fig.add_subplot(111)    # The big subplot
# ax1 = fig.add_subplot(211)
# ax2 = fig.add_subplot(212)

# # Turn off axis lines and ticks of the big subplot
# ax.spines["top"].set_color("none")
# ax.spines["bottom"].set_color("none")
# ax.spines["left"].set_color("none")
# ax.spines["right"].set_color("none")
# ax.tick_params(labelcolor="w", top="off", bottom="off", left="off", right="off")

# # Set common labels
# ax.set_xlabel("common xlabel")
# ax.set_ylabel("common ylabel")

# ax1.set_title('ax1 title')
# ax2.set_title('ax2 title')


text:
# Portions extracted from SLICOT-Reference distribution:
# Copyright (c) 2002-2020 NICONET e.V.

function run_mb04ed(datfile, io=stdout)
    NIN = 5
    NOUT = 6
    NMAX = 60
    LDB = NMAX÷2
    LDFG = NMAX÷2
    LDQ = NMAX
    LDU1 = NMAX÷2
    LDU2 = NMAX÷2
    LDWORK = 3*NMAX^2÷2 + max( NMAX, 24 ) + 3
    LDZ = NMAX
    LIWORK = NMAX + 9
    Z = Array{Float64,2}(undef, LDZ,NMAX)
    B = Array{Float64,2}(undef, LDB,NMAX÷2)
    FG = Array{Float64,2}(undef, LDFG,NMAX÷2+1)
    Q = Array{Float64,2}(undef, LDQ,NMAX)
    U1 = Array{Float64,2}(undef, LDU1,NMAX÷2)
    U2 = Array{Float64,2}(undef, LDU2,NMAX÷2)
    ALPHAR = Array{Float64,1}(undef, NMAX÷2)
    ALPHAI = Array{Float64,1}(undef, NMAX÷2)
    BETA = Array{Float64,1}(undef, NMAX÷2)
    IWORK = Array{BlasInt,1}(undef, LIWORK)
    DWORK = Array{Float64,1}(undef, LDWORK)
    f = open(datfile,"r")
    readline(f)
    vs = split(readline(f))
    JOB =  vs[1][1]
    COMPQ =  vs[2][1]
    COMPU =  vs[3][1]
    N = parse(BlasInt, vs[4])

    vs = String[]
    _isz,_jsz = (N,N)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       Z[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end

    vs = String[]
    _isz,_jsz = (N÷2,N÷2)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       B[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end

    vs = String[]
    _isz,_jsz = (N÷2,N÷2+1)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       FG[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end
if ( LSAME( COMPU, 'U' ) )

    vs = String[]
    _isz,_jsz = (N÷2,N÷2)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       U1[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end

    vs = String[]
    _isz,_jsz = (N÷2,N÷2)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       U2[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end
end # if
if ( N<0 || N>NMAX || mod( N, 2 )!=0 )
else
# interp call 1

    INFO = SLICOT.mb04ed!(JOB, COMPQ, COMPU, N, Z, B, FG, Q, U1, U2, ALPHAR, ALPHAI, BETA, LIWORK)
    @test INFO == 0
    INFO == 0 || return

if ( INFO!=0 )
else
# interp output 1
    println(io,"Z:")
    _nc = N
    _nr = N
    show(io,"text/plain",Z[1:_nr,1:_nc])
    println(io,)

# interp output 2
    println(io,"B:")
    _nc = N÷2
    _nr = N÷2
    show(io,"text/plain",B[1:_nr,1:_nc])
    println(io,)

# interp output 3
    println(io,"FG:")
    _nc = N÷2+1
    _nr = N÷2
    show(io,"text/plain",FG[1:_nr,1:_nc])
    println(io,)

# interp output 4
    println(io,"ALPHAR:")
    _nr = N÷2
    show(io,"text/plain",ALPHAR[1:_nr])
    println(io,)

# interp output 5
    println(io,"ALPHAI:")
    _nr = N÷2
    show(io,"text/plain",ALPHAI[1:_nr])
    println(io,)

# interp output 6
    println(io,"BETA:")
    _nr = N÷2
    show(io,"text/plain",BETA[1:_nr])
    println(io,)

# interp output 7
    println(io,"Q:")
    _nc = N
    _nr = N
    show(io,"text/plain",Q[1:_nr,1:_nc])
    println(io,)

if ( !LSAME( COMPU, 'N' ) )
# interp output 8
    println(io,"U1:")
    _nc = N÷2
    _nr = N÷2
    show(io,"text/plain",U1[1:_nr,1:_nc])
    println(io,)

# interp output 9
    println(io,"U2:")
    _nc = N÷2
    _nr = N÷2
    show(io,"text/plain",U2[1:_nr,1:_nc])
    println(io,)

end # if
end # if
end # if
    close(f)
end # run_X()


text:
using Decs
using Test

const d = [
    Dec(2, -1)
    Dec(2, -1)
    Dec(100, -4)
    Dec(1512, -2)
    Dec(-3, -2)
    Dec(-4, -6)
]

const testdecs =
[
 (Dec(1, -3),    "0.001",   0.001,       Float64),
 (Dec(1523, -2), "15.23",   15.23,       Float64),
 (Dec(543, 0),   "543",     543,         UInt),
 (Dec(-345, 0),  "-345",    -345,        Int),
 (Dec(123, 0),   "123",     123,         Int32),
 (Dec(-32, 0),   "-32",     -32,         Int8),
 (Dec(2001, 2),  "200100",  200100,      BigInt),
 (Dec(123, -2),  "1.23",    1.23,        Float64)
]

@testset "Decs" begin

@testset "Dec constructor" begin

@test isa(d, Vector{Dec})

end

@testset "Conversions" begin

@testset "String/Number to Dec" begin
    # Note: comparison to binary floats, where the value may not be representable exactly
    # may cause problems in the tests
    @testset "Direct" begin
        @test dec"0.01" == Dec(0.01) == Dec(1, -2)
        @test dec".001" == Dec(.001) == Dec(1, -3)
        @test dec"15.23" == Dec(15.23) == Dec(1523, -2)
        @test dec"543" == Dec(543) == Dec(543, 0)
        @test dec"-345" == Dec(-345) == Dec(-345, 0)
        @test dec"000123" == Dec(000123) == Dec(123, 0)
        @test dec"-00032" == Dec(-00032) == Dec(-32, 0)
        @test dec"200100" == Dec(200100) == Dec(2001, 2)
        @test dec"-.123" == Dec(-.123) == Dec(-123, -3)
        @test dec"1.23000" == Dec(1.23000) == Dec(123, -2)
        @test dec"4734.612" == Dec(4734.612) == Dec(4734612, -3)
        @test dec"541724.2" == Dec(541724.2) == Dec(5417242, -1)
        @test dec"2.5e6" == Dec(2.5e6) == Dec(25, 5)
        @test dec"2.385350e8" == Dec(2.385350e8) == Dec(238535, 3)
        @test dec"12.3e-4" == Dec(12.3e-4) == Dec(123, -5)

        @test dec"-12.3e4" == Dec(-12.3e4) == Dec(-123, 3)

        @test dec"-12.3e-4" == Dec(-12.3e-4) == Dec(-123, -5)

        @test dec"0.1234567891" == Dec(0.1234567891) == Dec(1234567891, -10)
        @test dec"0.12345678912" == Dec(0.12345678912) == Dec(12345678912, -11)
    end

    @testset "Using `dec`" begin
        @test dec("1.0") == Dec(1, 0)
        @test dec(8.1) == Dec(81, -1)
        @test dec.(Float64.(d)) == d
    end
end

@testset "Array{<:Number} to Array{Dec}" begin
    @test Dec.([0.1 0.2 0.3]) == [Dec(0.1) Dec(0.2) Dec(0.3)]
end

@testset "Dec to String" begin
    for (d, s, v, T) in testdecs
        @test string(d) == s
        @test T(d) == v
        #@test number(d) == v
    end
    @test string(Float64(Dec(543, 0))) == "543.0"
    #@test string(number(Dec(543, 0)))  == "543"
    #@test string(number(Dec(543, -1))) == "54.3"
end

@testset "Float32" begin
    # Note that 0.01f0 != 0.01
    d = Dec(1, -2)
    @test string(d) == "0.01"
    @test Float32(d) == 0.01f0
    #@test number(d) == 0.01
end

@testset "BigFloat" begin
    d = Dec(-123, -3)
    @test string(d) == "-0.123"
    @test BigFloat(d) == big"-0.123"
    #@test number(d) == -0.123
end

end # constructor

@testset "Normalization" begin

@test Dec(-151100, -4) == Dec(-1511, -2)
@test Dec(100100, -5) == Dec(1001, -3)
#@test normalize(Dec(-151100, -4)) == Dec(-1511, -2)
#@test normalize(Dec(100100, -5)) == Dec(1001, -3)
@test dec"3.0"    == Dec(3, 0)
@test dec"3.0"    == Dec(30, -1)
@test dec"3.1400" == Dec(314, -2)
@test dec"1234"   == Dec(1234, 0)

end # Normalization

@testset "Arithmetic" begin

@testset "Addition" begin
    @test Dec(0.1) + 0.2 == 0.1 + Dec(0.2) == Dec(0.1) + Dec(0.2) == Dec(0.3)
    @test Dec.([0.1 0.2]) .+ [0.3 0.1] == Dec.([0.4 0.3])
    @test Dec(2147483646) + Dec(1) == Dec(2147483647)
    @test Dec(-3, -2) + dec"0.2523410412138103" == Dec(2223410412138103, -16)
end

@testset "Subtraction" begin
    @test Dec(0.3) - 0.1 == 0.3 - Dec(0.1)
    @test 0.3 - Dec(0.1) == Dec(0.3) - Dec(0.1)
    @test Dec(0.3) - Dec(0.1) == Dec(0.2)
    @test Dec.([0.3 0.1]) .- [0.1 0.5] == Dec.([0.2 -0.4])
end

@testset "Negation" begin
    @test -Dec.([0.3 0.2]) == [-Dec(0.3) -Dec(0.2)]
    @test -Dec(0.3) == zero(Dec) - Dec(0.3)
    @test iszero(dec(12.1) - dec(12.1))
end

@testset "Multiplication" begin
    @test Dec(12.21) * Dec(2.12) == Dec(0, -4, 258852)
    @test Dec(12.2112543) * Dec(2.121352) == Dec(0, -13, 259043687318136)
    @test Dec(0.2) * 0.1 == 0.2 * Dec(0.1)
    @test 0.2 * Dec(0.1) == Dec(0.02)
    @test Dec(12.34) * 0.1234 == 12.34 * Dec(0.1234)
    @test 12.34 * Dec(0.1234) == Dec(1.522756)
    @test Dec(0.21084210) * -2 == -2 * Dec(0.21084210)
    @test -2 * Dec(0.21084210) == Dec(-0.4216842)
    @test Dec(0, 2, -1) * 0.0 == zero(Dec)
    @test Dec.([0.3, 0.6]) .* 5 == [Dec(0.3)*5, Dec(0.6)*5]
    @test one(Dec) * 1 == Dec(1, 0)
end

@testset "Inversion" begin
    @test inv(Dec(1, -1))   == Dec(1, 1)
    @test inv(Dec(1, 1))    == Dec(1, -1)
    @test inv(Dec(-2, -1))  == Dec(-5, 0)
    @test inv(Dec(-5, 0))   == Dec(-2, -1)
    @test inv(Dec(2, -2))   == Dec(5, 1)
    @test inv(Dec(5, 1))   == Dec(2, -2)
    @test inv(Dec(-4, -1))  == Dec(-25, -1)
    @test inv(Dec(-25, -1)) == Dec(-4, -1)
end

@testset "Division" begin
    @test Dec(0.2) / Dec(0.1) == Dec(2)
    @test Dec(0.3) / Dec(0.1) == Dec(3, 0)
    @test [Dec(0.3) / Dec(0.1), Dec(0.6) / Dec(0.1)] == [Dec(0.3), Dec(0.6)] ./ Dec(0.1)
    @test [Dec(0.3) / 0.1, Dec(0.6) / 0.1] == [Dec(0.3), Dec(0.6)] ./ 0.1
end

end # Arithmetic

@testset "Equality" begin

@testset "isequal" begin
    @test isequal(Dec(0, -3, 2), Dec(0, -3, 2))
    @test !isequal(Dec(0, -3, 2), Dec(0, 3, 2))
    @test isequal(Dec(0, -3, 2), 0.002)
    @test isequal(Dec(1, 0, 2), -2)
    @test !isequal(Dec(1, 0, 2), 2)
    @test !isequal(Dec(1, -1, 0), Dec(0, 0, 0))
end

@testset "==" begin
    @test Dec(2, -3) == Dec(2, -3)
    @test Dec(2, -3) != Dec(2, 3)
    @test Dec(2, -3) == 0.002

    @test -2 == Dec(-2, 0)
    @test 2 != Dec(-2, 0)

    @test Dec(-2, 0) == -2
    @test Dec(-2, 0) != 2

    bf_pi = BigFloat(pi)
    @test Dec(bf_pi) == bf_pi
    @test bf_pi == Dec(bf_pi)

    bi = big"4608230166434464229556241992703"
    @test Dec(bi) == bi
    @test bi == Dec(bi)

    @test dec(12.1) == dec(12.1)

    # Test negative zero
    @test Dec(1, -1, 0) == Dec(0)
end

@testset "<" begin
    @test Dec(-1, 1)      < Dec(1, 1)
    @test !(Dec(1, 1)     < Dec(-1, 1))
    @test Dec(2, -3)      < Dec(2, 3)
    @test !(Dec(2, 3)     < Dec(2, -3))
    @test !(dec(12.1)     < dec(12.1))

    # Tests with negative zero
    @test Dec(1, 1, 1)    < Dec(1, 1, 0)
    @test !(Dec(1, 1, 0)  < Dec(1, 1, 1))
    @test !(Dec(1, -1, 0) < Dec(0))
    @test !(Dec(0)        < Dec(1, -1, 0))
end

end # comparisons

@testset "Rounding" begin

@test round(Dec(7.123456), digits=0) == Dec(7)
@test round(Dec(7.123456), digits=2) == Dec(7.12)
@test round(Dec(7.123456), digits=3) == Dec(7.123)
@test round(Dec(7.123456), digits=5) == Dec(7.12346)
@test round(Dec(7.123456), digits=6) == Dec(7.123456)

@test round.(Dec.([0.1111, 0.2222, 0.8888]), digits=2) == Dec.([0.11, 0.22, 0.89])

function tet()
    a = dec"1.0000001"
    for i = 1:27
        a *= a
    end
    return a
end

# set DIGITS = 20 (aaljuffali's example)
#@test tet() == Dec(0, -20, 67453047074102193157641340)

end # Rounding
end # Decs


text:
Base.include(Main, "scopf.jl")

# to obtain a list of off-lines run:  grep "optimal" summaries/case30_SCOPF_summary | awk '{print $2;}' | tr ':\n' ' '

function main()
  if length(ARGS) < 1
    println("Usage: julia scopf_main.jl case_name lines_indexes")
    println("Cases are in 'data' directory: case9 case30 case57 case118 case300 case1354pegase case2383wp case2736sp case2737sop case2746wop case2869pegase case3012wp  case3120sp case3375wp case9241pegase")
    return
  end

  opfdata = opf_loaddata(ARGS[1])

  lines_off=Array{Line}(undef, length(ARGS)-1)
  for l in 1:length(lines_off)
    lines_off[l] = opfdata.lines[parse(Int,ARGS[l+1])]
  end
  scopfdata = SCOPFData(opfdata,lines_off)
  @assert length(lines_off) == length(ARGS)-1
  scopfmodel = scopf_model(scopfdata)
  opfmodel,status = scopf_solve(scopfmodel,scopfdata)
  if status==:Optimal
    scopf_outputAll(opfmodel,scopfdata)
  end
end

main()


text:
module YaoQX

using YaoAPI

include("register.jl")

end


text:
""""Macro which adds NVTX range only if CUDA is functional"""
macro nvtx_range(label, ex)
    if CUDA.functional()
        return quote
            NVTX.@range $(esc(label)) $(esc(ex))
        end
    else
        return quote
            $(esc(ex))
        end
    end
end


text:
#=
Copyright (c) 2015, Intel Corporation

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

    * Redistributions of source code must retain the above copyright notice,
      this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of Intel Corporation nor the names of its contributors
      may be used to endorse or promote products derived from this software
      without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=#

module StructureAnalysisStructOnly

    using Sparso: Sym, Symexpr, TypedExprNode
    using ..SymbolicAnalysis

    function symbolize(e :: Symexpr, tp :: Type, unique)
        if unique == true
            sym = new_symbol("unknown")
            return MiddleSymbol(sym)
        end

        if isa(e, Sym)
            if tp <: SparseMatrixCSC
                s = MiddleSymbol(:true)
            else
                dump(tp)
                assert(0)
            end
            return s
        else
            assert(0)
        end
    end

    function preprocess(property_proxies, symbol_info)
        predefined = Dict{Sym, AbstractSymbol}()
        for (s, v) in property_proxies
            if isa(v.structure_only, MiddleSymbol)
                assert(v.structure_only.value == :true)
                predefined[s] = v.structure_only
            end
        end
        predefined
    end

    function postprocess(res, property_proxies, symbol_info)
        for (s, v) in res
            assert(isa(s, Sym))
            if isa(v, MiddleSymbol) && v.value == :true
                property_proxies[s].structure_only = v
            end
        end
    end

    function set_true_action(e)
        e.svalue = MiddleSymbol(:true)
    end

    function pass_a1_action(e)
        e.svalue = e.args[1].svalue
    end

    function assign_action(e)
        e.args[1].svalue = e.args[2].svalue
    end

    function set_elem_action(e)
        val = e.args[2].raw_expr
        if !isa(val, Int) || (val != 0 && val != 1)
            e.args[1].svalue = MiddleSymbol(:false)
        end
    end

    const transfer_rules = (
        ((:(=), Any, Any), assign_action),

        ((:call, GlobalRef(Main, :spones), Any), set_true_action),
        ((:call, GlobalRef(Main, :zeros), Any, Any), set_true_action),
        ((:call, GlobalRef(Main, :sparse), Any), pass_a1_action),
        ((:call, GlobalRef(Main, :scale), SparseMatrixCSC, Any), pass_a1_action),
        ((:call, GlobalRef(Main, :setindex!), Array, Any, Any, Any), set_elem_action),

        ((:call, TypedExprNode(Function, :call, TopNode(:apply_type), GlobalRef(Main, :SparseMatrixCSC), GlobalRef(Main, :Float64), GlobalRef(Main, :Int32)), Any), pass_a1_action),
        ((:call, TypedExprNode(Function, :call, TopNode(:apply_type), GlobalRef(Main, :SparseMatrixCSC), GlobalRef(Main, :Cdouble), GlobalRef(Main, :Cint)), Any), pass_a1_action),
    )

    const pass_info = ("StructOnly", transfer_rules, preprocess, postprocess, symbolize)
end

max_stars_repo_name: IntelLabs/Sparso

max_issues_repo_name: IntelLabs/Sparso

max_forks_repo_name: IntelLabs/Sparso


text:
using Qsosed, Test

@testset "Test corona properties" begin
    # comparison with Xspec model

    bh = BlackHole(1e8 * M_SUN, 0.5, 0.0)
    corona =
        Corona(bh, hard_xray_fraction = 0.02, electron_energy = 100, reprocessing = true)
    Ledd = compute_eddington_luminosity(bh)

    computed_radius = compute_corona_radius(corona)
    @test computed_radius ≈ 8.9008961989970832 rtol = 1e-3

    temp_corona = compute_disk_temperature(bh, computed_radius)
    @test temp_corona ≈ 107102.98423666063 rtol = 1e-1

    dissip_lumin = compute_corona_dissipated_luminosity(corona) / Ledd
    @test dissip_lumin ≈ 2.0183853632627917e-2 rtol = 1e-2

    corona_luminosity = compute_corona_luminosity(corona) / Ledd
    @test corona_luminosity ≈ 3.6390444796974425e-2 rtol = 1e-1

    corona_photon_index = compute_corona_photon_index(corona)
    @test corona_photon_index ≈ 2.2826826043284854 rtol = 1e-2

    # Second black hole

    bh = BlackHole(1e6 * M_SUN, 0.025, 0.0)
    @test gravity_radius(bh) ≈ 1161.0642639605840 rtol = 1e-2
    corona =
        Corona(bh, hard_xray_fraction = 0.02, electron_energy = 100, reprocessing = true)
    Ledd = compute_eddington_luminosity(bh)

    computed_radius = compute_corona_radius(corona)
    @test computed_radius ≈ 89.731710842601927  rtol = 1e-2

    temp_corona = compute_disk_temperature(bh, computed_radius)
    @test temp_corona ≈ 47643.787634112188 rtol = 1e-1

    dissip_lumin = compute_corona_dissipated_luminosity(corona) / Ledd
    @test dissip_lumin ≈ 2.0001856311471888e-2 rtol = 1e-2

    corona_luminosity = compute_corona_luminosity(corona) / Ledd
    @test corona_luminosity ≈ 2.0621184497821025e-2 rtol = 1e-1

    corona_photon_index = compute_corona_photon_index(corona)
    @test corona_photon_index ≈ 1.6483970267857697  rtol = 1e-1

    warm = Warm(corona, electron_energy = 0.2, photon_index = 2.5)
    energy_range = 10 .^ range(-5, 3, length = 500);
    corona_sed = compute_corona_spectral_luminosity(corona, warm, energy_range)
    integ_lumin = sum(corona_sed[2:end] .* diff(energy_range))
    target_lumin = compute_corona_luminosity(corona)
    @test integ_lumin ≈ target_lumin rtol = 1e-5

end


text:
function amino_acid_sequence(dna)
    amino_acids = ""

    # Your code goes here!

    return amino_acids
end

max_stars_repo_name: basimr/lovelace

max_issues_repo_name: basimr/project-lovelace

max_forks_repo_name: basimr/project-lovelace


text:
using Revise
using JuMP
using GLPK
using ModelGraphs

include("DualDecompositionSolver.jl")
using .DualDecompositionSolver

DDSolver = DualDecompositionSolver

m2 = Model(with_optimizer(GLPK.Optimizer))
@variable(m2, xs[i in 1:2])
@variable(m2, y[i in 1:2] >= 2)
@constraint(m2, y[1] + y[2] <= 5)
@constraint(m2, xs[1] + xs[2] + y[1] + y[2] <= 10)
@objective(m2, Min, xs[1] + xs[2])

m3 = Model(with_optimizer(GLPK.Optimizer))
@variable(m3, xs[i in 1:2])
@variable(m3, y[i in 1:2] >= 0)
@constraint(m3, y[1] + y[2] <= 1)
@constraint(m3, xs[1] + xs[2] + y[1] + y[2] <= 10)
@objective(m3, Min, y[1] + y[2])

## Model Graph
graph = ModelGraph()
JuMP.set_optimizer(graph.mastermodel,with_optimizer(GLPK.Optimizer))

@linkvariable(graph, xm[i in 1:2] >= 2)
@masterconstraint(graph, xm[1] + xm[2] <= 10)
#@objective(graph.mastermodel, Max, 16*xm[1].vref + 10*xm[2].vref)

n1 = add_node!(graph,m2)
n2 = add_node!(graph,m3)

for i = 1:2
    link_variables!(graph[:xm][i],n1[:xs][i])
    link_variables!(graph[:xm][i],n2[:xs][i])
end

# #verify with full problem
glpk = with_optimizer(GLPK.Optimizer)
m,ref_map = aggregate(graph)
optimize!(m,glpk)

println()
println("Pure glpk solution:")
println(JuMP.value.(all_variables(m)))

dd_model = DDSolver.DDModel(graph)

for subproblem in dd_model.subproblems
    JuMP.set_optimizer(subproblem,glpk)
end

DDSolver.optimize!(dd_model)
#print solution
println("Dual decomposition solution: ")
i = 1
for subproblem in dd_model.subproblems
    println("Subproblem $i")
    println(value.(all_variables(subproblem)))
    global i += 1
end

# link_eq_matrix,b_eq,link_eq_variables = DDSolver.prepare_link_matrix(link_eq_constraints,sub_var_map)
# link_ineq_matrix,b_ineq,link_ineq_variables, = prepare_link_matrix(link_ineq_constraints,sub_var_map)
# link_var_matrix,link_var_variables = prepare_link_matrix(link_variables,master_var_map)

#link_var_matrix,link_var_variables,link_var_map = DDSolver.prepare_link_matrix(linkvariables)


text:
module Embeddings

using Statistics: norm
using DataDeps
using AutoHashEquals

export load_embeddings, language_files
export Word2Vec, GloVe, FastText_Text

abstract type EmbeddingSystem{LANG} end
include("fasttext.jl")
include("glove.jl")
include("word2vec.jl")

function __init__()
    for T in [Word2Vec, GloVe, FastText]
        init(T)
    end
end

@auto_hash_equals struct EmbeddingTable{M<:AbstractMatrix, A<:AbstractVector}
    embeddings::M
    vocab::A
end

const supported_languages_files = Dict{DataType, Vector{String}}()

load_datadep(datadep_name) = @datadep_str datadep_name # makes sure to get  __FILE__ right

default_file(::Type{T}, file_number) where T = language_files(T)[file_number] |> load_datadep

function language_files(::Type{T}) where {T<:EmbeddingSystem{LANG}} where LANG
    get!(supported_languages_files,T) do
        String[]
    end
end

# Catch the unionall where no type-param given
language_files(::Type{T}) where T =  language_files(T{:en})

"""

    load_embeddings(EmbeddingSystem, [embedding_file|default_file_number])
    load_embeddings(EmbeddingSystem{:lang}, [embedding_file|default_file_number])

Loaded the embeddings from a embedding file.
The embeddings should be of the type given by the Embedding system.

If the `embedding file` is not provided, a default embedding file will be used.
(It will be automatically installed if required).
EmbeddingSystems have a language type parameter.
For example `FastText_Text{:fr}` or `Word2Vec{:en}`, if that language parameter is not given it defaults to English.
(I am sorry for the poor state of the NLP field that many embedding formats are only available pretrained in English.)
Using this the correct  default embedding file will be installed for that language.
For some languages and embedding systems there are multiple possible files.
You can check the list of them using for example `language_files(FastText_Text{:de})`.
The first is nominally the most popular, but if you want to default to another you can do so by setting the `default_file_num`.

### Keyword Arguments:

 - `max_vocab_size` an integer, it specifies the maximum number of words to load (most formats are sorted by frequency so this keeps the the most common words). Default is to keep all of them
 - `keep_words=Set()` if a non-empty set of words is provided, then only word embeddings for words from that list are loaded. Otherwise (default) all words are loaded.

### Returns an `Embeddings` object.
This has 2 fields.

 - `embeddings` is a matrix, each column is the embedding for a word.
 - `vocab` is a vector of strings, ordered as per the columns of `embeddings`, such that the first string in vocab is the first column of `embeddings` etc

We do not include a method for getting the index of a column from a word.
This is trivial to define in code (`vocab2ind(vocab)=Dict(word=>ii for (ii,word) in enumerate(vocab))`),
and you might like to be doing this in a more consistant way, e.g using [MLLabelUtils.jl](https://github.com/JuliaML/MLLabelUtils.jl),
or you might like to build a much faster Dict solution on top of [InternedStrings.jl](https://github.com/JuliaString/InternedStrings.jl)

```

"""
function load_embeddings(::Type{T},
        default_file_number::Int=1;
        max_vocab_size=typemax(Int),
        keep_words=Set()) where T<:EmbeddingSystem

    embedding_file = default_file(T, default_file_number)
    load_embeddings(T, embedding_file; max_vocab_size=max_vocab_size, keep_words=keep_words)
end

function load_embeddings(::Type{T},
        embedding_file::AbstractString;
        max_vocab_size=typemax(Int),
        keep_words=Set()) where T<:EmbeddingSystem

    EmbeddingTable(_load_embeddings(T, embedding_file, max_vocab_size,  keep_words)...)
end

end