text:
#
# Me Learning Machine Learning
# Chapter 1: Julia and Knet
#
# In summary, I am (loosely) following the Knet tutorial
# (https://github.com/denizyuret/Knet.jl/blob/master/README.md#Tutorial-1)
# and adding lots of comments, including some of these noob-style "this line
# does x, y and z" comments (because I am a Julia noob :-) ).
#
# Leandro Motta Barros
#

using Knet

# ------------------------------------------------------------------------------
# Reading and pre-processing data
# ------------------------------------------------------------------------------

# This is based on the classic "Boston Housing" dataset, which I downloaded from
# https://archive.ics.uci.edu/ml/machine-learning-databases/housing/housing.data
# (see the file `housing.names` on the same server and directory for more info).
# This database contains 506 lines by 14 columns.

# Read the data into an array. By default, `readdlm` assumes that fields are
# separated by blanks -- quite handy, 'cause that's exactly what we have here.
rawdata = readdlm("housing.data") # [506 x 14]

# Take only the input (independent variables), which is everything but the 14th
# column. We also transpose the matrix, which at first seemed weird for me. It
# turns out that Julia stores arrays in column-major order, so having the data
# transposed will make it more efficient to iterate over the training examples.
x = rawdata[:,1:13]' # [13 x 506]

# Normalization of input variables. The operation itself is pretty standard,
# though the implementation might deserve a comment or two. Here (as in the
# previous line, indeed), there is a lot of Matlab/Octave taste. "Dot-operator"
# to execute element-wise operations. Second parameter to `mean` and `std` (`2`)
# meaning that the mean and standard deviation are to be made along the second
# dimension of the array. In other words, along all elements in a single column.
# Which, in still other words, is the same as saying that we are computing the
# mean and standard deviation among all values of each feature (recall that `x`
# was transposed in the previous line).
x = (x .- mean(x,2)) ./ std(x,2) # [13 x 506]

# Now take the expected output, which is on the 14th (and last) column in our
# data. Again, transpose.
y = rawdata[:,14:14]'

# ------------------------------------------------------------------------------
# Linear regression
# ------------------------------------------------------------------------------

# Initialize the weights. `randn` generates Gaussian-distributed random numbers
# with mean equals to zero and standard deviation of 1.0.

# TODO: Why this explicit array of `Any`? Why that additional zero?

w = Any[ 0.1 * randn(1,13), 0 ]

function train(w, data; lr = .1)
    for (x,y) in data
        dw = lossgradient(w, x, y)
        for i in 1:length(w)
            w[i] -= lr * dw[i]
        end
    end
    return w
end

show(w)

#print(x)


text:
@time @testset "Aqua (Auto QUality Assurance)" begin
    Aqua.test_all(
        PredictMDFuture;
        ambiguities       = false,
        stale_deps        = false,
        undefined_exports = false,
    )
end

ext: jl


text:
myFunction (a{T};
            x::Float64=sin(1))::Bool = false


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

using Markdown
using InteractiveUtils

# ╔═╡ ee0915fd-eacc-47ef-8ada-fd83909bae18
begin
	using PlutoUI
	using DataStructures
end

# ╔═╡ aa81ee70-e559-11eb-17be-9b561a8ae8db
md"""
# Chapter 9
#### Nicetak
 $(import Dates; Dates.format(Dates.today(), Dates.DateFormat("U d, Y")))
"""

# ╔═╡ d975985b-84e7-4d2c-968e-eb2f7bcff5c2
md"""
### 9.1
連結リストを用いてスタックとキューを実現してください.
"""

# ╔═╡ 585e300b-0844-4867-a252-b58fff71b63a
begin
	mutable struct MyStack{T}
		l::MutableLinkedList{T}

		function MyStack{T}() where T
			mst = new{T}()
			mst.l = MutableLinkedList{T}()

			return mst
		end
	end

	MyStack() = MyStack{Any}()

	function Base.push!(mst::MyStack, x)
		push!(mst.l, x)
		return mst
	end

	function Base.pop!(mst::MyStack)
		pop!(mst.l)
		return mst
	end

	function Base.show(io::IO, mst::MyStack)
    	print(io, typeof(mst), '(')
	    join(io, mst.l, ", ")
    	print(io, ')')
	end
end

# ╔═╡ be538624-5d9c-4849-ba32-2e0055764805
begin
	mutable struct MyQueue{T}
		l::MutableLinkedList{T}

		function MyQueue{T}() where T
			mq = new{T}()
			mq.l = MutableLinkedList{T}()

			return mq
		end
	end

	MyQueue() = MyQueue{Any}()

	function DataStructures.enqueue!(mq::MyQueue, x)
		push!(mq.l, x)
		return mq
	end

	function DataStructures.dequeue!(mq::MyQueue)
		popfirst!(mq.l)
		return mq
	end

	function Base.show(io::IO, mq::MyQueue)
    	print(io, typeof(mq), '(')
	    join(io, mq.l, ", ")
    	print(io, ')')
	end
end

# ╔═╡ 88bfe770-9756-4c79-be16-67c92798d700
with_terminal() do
	@show mst = MyStack{Int}()

	@show push!(mst, 3)
	@show push!(mst, 5)
	@show push!(mst, 7)

	@show pop!(mst)
	@show pop!(mst)

	@show push!(mst, 9)

end

# ╔═╡ ad3ad687-f324-42d3-b3ab-c1264cf1137d
with_terminal() do
	@show mq = MyQueue{Int}()

	@show enqueue!(mq, 3)
	@show enqueue!(mq, 5)
	@show enqueue!(mq, 7)

	@show dequeue!(mq)
	@show dequeue!(mq)

	@show enqueue!(mq, 9)

end

# ╔═╡ c54762ef-f758-470e-9329-7943e73f1057
md"""
## 9.2
逆ポーランド記法で記述された数式を入力として受け取って, その計算結果を出力するアルゴリズムを設計してください.
"""

# ╔═╡ c37576ea-fdab-49ae-8967-f525ca26cb0f
function rpn(s::String)
	st = Stack{Any}()
	for e in eval.(Meta.parse.(split(s)))
		if isa(e, Function)
			arg2 = pop!(st)
			arg1 = pop!(st)
			push!(st, e(arg1, arg2))
		else
			push!(st, e)
		end
	end

	return pop!(st)
end

# ╔═╡ f2be00f6-b8eb-4337-a83a-e73ce805bfdf
with_terminal() do
	ss = [
		"1 2 + 3 4 + +",
		"3 4 * 1 * 2 -",
		"2 3 * 1 4 * +"
	]
	for s in ss
		println("\"$s\" = $(rpn(s))")
	end
end

# ╔═╡ e507c176-96e0-48b0-9772-ac72519f1be0
md"""
### 9.3
"(()(())())(()())"のような'('と')'からなる, 長さ$2N$の文字列が与えられます. この文字列において括弧列が整合しているかどうかを判定し, さらに何文字目と何文字目の括弧が対応しているかを$N$組求める処理を,
$O(N)$で実行するアルゴリズムを設計してください.
"""

# ╔═╡ b7ce3689-afe5-4854-a4dd-7d937bab4fdc
function parenthesis(s::String)
	N = length(s) ÷ 2

	st = Stack{Int}()
	ans = zeros(Int, N, 2)
	idx = 1

	for i in eachindex(s)
		if s[i] == '('
			push!(st, i)
		else
			if isempty(st)
				return println("invalid")
			else
				j = pop!(st)
				ans[idx, :] = [j i]
				idx += 1
			end
		end
	end

	if isempty(st)
		println(ans)
	else
		println("invalid")
	end
end

# ╔═╡ e3c9bb00-6bc4-44f7-9744-6505f2a2797b
with_terminal() do
	ss = ["(((())())())", "(()", "())"]
	for s ∈ ss
		print("\"$s\" → ")
		parenthesis(s)
	end
end

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
DataStructures = "~0.18.10"
PlutoUI = "~0.7.9"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "727e463cfebd0c7b999bbf3e9e7e16f254b94193"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.34.0"

[[DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "7d9d316f04214f7efdbb6398d545446e246eff02"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.10"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[LinearAlgebra]]
deps = ["Libdl"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "438d35d2d95ae2c5e8780b330592b6de8494e779"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.0.3"

[[Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[PlutoUI]]
deps = ["Base64", "Dates", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "Suppressor"]
git-tree-sha1 = "44e225d5837e2a2345e69a1d1e01ac2443ff9fcb"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.9"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[Suppressor]]
git-tree-sha1 = "a819d77f31f83e5792a76081eee1ea6342ab8787"
uuid = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
version = "0.2.0"

[[TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╟─aa81ee70-e559-11eb-17be-9b561a8ae8db
# ╠═ee0915fd-eacc-47ef-8ada-fd83909bae18
# ╟─d975985b-84e7-4d2c-968e-eb2f7bcff5c2
# ╠═585e300b-0844-4867-a252-b58fff71b63a
# ╠═be538624-5d9c-4849-ba32-2e0055764805
# ╟─88bfe770-9756-4c79-be16-67c92798d700
# ╟─ad3ad687-f324-42d3-b3ab-c1264cf1137d
# ╟─c54762ef-f758-470e-9329-7943e73f1057
# ╠═c37576ea-fdab-49ae-8967-f525ca26cb0f
# ╟─f2be00f6-b8eb-4337-a83a-e73ce805bfdf
# ╟─e507c176-96e0-48b0-9772-ac72519f1be0
# ╠═b7ce3689-afe5-4854-a4dd-7d937bab4fdc
# ╟─e3c9bb00-6bc4-44f7-9744-6505f2a2797b
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
"""
An SCOPF formulation conforming to the ARPA-e GOC Challenge 1 specification.
A DC power flow approximation is used. Power balance and line flow constraints
are strictly enforced in the first stage.  Contingency branch flow constraints
are enforced by PTDF cuts and penalized based on a conservative linear
approximation of the formulation's specification.

This formulation is used in conjunction with the contingency filters that
generate PTDF cuts.
"""
function run_scopf_cuts_soft(file, model_constructor, solver; kwargs...)
    return run_model(file, model_constructor, solver, build_scopf_cuts_soft; kwargs...)
end

""
function build_scopf_cuts_soft(pm::AbstractPowerModel)
    PowerModels.variable_voltage(pm)
    PowerModels.variable_generation(pm)
    PowerModels.variable_branch_flow(pm)

    variable_branch_contigency_flow_violation(pm)
    variable_gen_contigency_flow_violation(pm)
    variable_gen_contigency_capacity_violation(pm)

    for i in ids(pm, :bus)
        expression_bus_generation(pm, i)
        expression_bus_withdrawal(pm, i)
    end

    PowerModels.constraint_model_voltage(pm)

    for i in ids(pm, :ref_buses)
        PowerModels.constraint_theta_ref(pm, i)
    end

    for i in ids(pm, :bus)
        PowerModels.constraint_power_balance(pm, i)
    end

    for i in ids(pm, :branch)
        PowerModels.constraint_ohms_yt_from(pm, i)
        PowerModels.constraint_ohms_yt_to(pm, i)

        PowerModels.constraint_voltage_angle_difference(pm, i)

        PowerModels.constraint_thermal_limit_from(pm, i)
        PowerModels.constraint_thermal_limit_to(pm, i)
    end

    for (i,cut) in enumerate(ref(pm, :branch_flow_cuts))
        constraint_branch_contingency_ptdf_thermal_limit_from_soft(pm, i)
        constraint_branch_contingency_ptdf_thermal_limit_to_soft(pm, i)
    end

    bus_withdrawal = var(pm, :bus_wdp)

    for (i,cut) in enumerate(ref(pm, :gen_flow_cuts))
        branch = ref(pm, :branch, cut.branch_id)
        gen = ref(pm, :gen, cut.gen_id)
        gen_bus = ref(pm, :bus, gen["gen_bus"])
        gen_set = ref(pm, :area_gens)[gen_bus["area"]]
        alpha_total = sum(gen["alpha"] for (i,gen) in ref(pm, :gen) if gen["index"] != cut.gen_id && i in gen_set)

        cont_bus_injection = Dict{Int,Any}()
        for (i, bus) in ref(pm, :bus)
            inj = 0.0
            for g in ref(pm, :bus_gens, i)
                if g != cut.gen_id
                    if g in gen_set
                        inj += var(pm, :pg, g) + gen["alpha"]*var(pm, :pg, cut.gen_id)/alpha_total
                    else
                        inj += var(pm, :pg, g)
                    end
                end
            end
            cont_bus_injection[i] = inj
        end

        #rate = branch["rate_a"]
        rate = branch["rate_c"]
        @constraint(pm.model,  sum( weight*(cont_bus_injection[bus_id] - bus_withdrawal[bus_id]) for (bus_id, weight) in cut.bus_injection) <= rate + var(pm, :gen_cont_flow_vio, i))
        @constraint(pm.model, -sum( weight*(cont_bus_injection[bus_id] - bus_withdrawal[bus_id]) for (bus_id, weight) in cut.bus_injection) <= rate + var(pm, :gen_cont_flow_vio, i))
    end

    for (i,gen_cont) in enumerate(ref(pm, :gen_contingencies))
        #println(gen_cont)
        gen = ref(pm, :gen, gen_cont.idx)
        gen_bus = ref(pm, :bus, gen["gen_bus"])
        gen_set = ref(pm, :area_gens)[gen_bus["area"]]
        response_gens = Dict(g => ref(pm, :gen, g) for g in gen_set if g != gen_cont.idx)

        # factor of 1.2 accounts for losses in a DC model
        #@constraint(pm.model, sum(gen["pmax"] - var(pm, :pg, g) for (g,gen) in response_gens) >= 1.2*var(pm, :pg, gen_cont.idx))
        @constraint(pm.model, var(pm, :gen_cont_cap_vio, i) + sum(gen["pmax"] - var(pm, :pg, g) for (g,gen) in response_gens) >= var(pm, :pg, gen_cont.idx))
        #@constraint(pm.model, sum(gen["pmin"] - var(pm, :pg, g) for (g,gen) in response_gens) <= var(pm, :pg, gen_cont.idx))
    end

    ##### Setup Objective #####
    objective_variable_pg_cost(pm)
    # explicit network id needed because of conductor-less
    pg_cost = var(pm, pm.cnw, :pg_cost)
    branch_cont_flow_vio = var(pm, pm.cnw, :branch_cont_flow_vio)
    gen_cont_flow_vio = var(pm, pm.cnw, :gen_cont_flow_vio)
    gen_cont_cap_vio = var(pm, pm.cnw, :gen_cont_cap_vio)

    @objective(pm.model, Min,
        sum( pg_cost[i] for (i,gen) in ref(pm, :gen) ) +
        sum( 5e5*branch_cont_flow_vio[i] for i in 1:length(ref(pm, :branch_flow_cuts)) ) +
        sum( 5e5*gen_cont_flow_vio[i] for i in 1:length(ref(pm, :gen_flow_cuts)) ) +
        sum( 5e5*gen_cont_cap_vio[i] for i in 1:length(ref(pm, :gen_contingencies)) )
    )
end

"a variant of `run_scopf_cuts_dc_soft` with a different generator response function"
function run_scopf_cuts_soft_bpv(file, model_constructor, solver; kwargs...)
    return run_model(file, model_constructor, solver, build_scopf_cuts_soft_bpv; kwargs...)
end

""
function build_scopf_cuts_soft_bpv(pm::AbstractPowerModel)
    PowerModels.variable_voltage(pm)
    PowerModels.variable_generation(pm)
    PowerModels.variable_branch_flow(pm)

    variable_branch_contigency_flow_violation(pm)
    variable_gen_contigency_flow_violation(pm)

    for i in ids(pm, :bus)
        expression_bus_generation(pm, i)
        expression_bus_withdrawal(pm, i)
    end

    PowerModels.constraint_model_voltage(pm)

    for i in ids(pm, :ref_buses)
        PowerModels.constraint_theta_ref(pm, i)
    end

    for i in ids(pm, :bus)
        PowerModels.constraint_power_balance(pm, i)
    end

    for i in ids(pm, :branch)
        PowerModels.constraint_ohms_yt_from(pm, i)
        PowerModels.constraint_ohms_yt_to(pm, i)

        PowerModels.constraint_voltage_angle_difference(pm, i)

        PowerModels.constraint_thermal_limit_from(pm, i)
        PowerModels.constraint_thermal_limit_to(pm, i)
    end

    for (i,cut) in enumerate(ref(pm, :branch_flow_cuts))
        constraint_branch_contingency_ptdf_thermal_limit_from_soft(pm, i)
        constraint_branch_contingency_ptdf_thermal_limit_to_soft(pm, i)
    end

    bus_withdrawal = var(pm, :bus_wdp)

    for (i,cut) in enumerate(ref(pm, :gen_flow_cuts))
        branch = ref(pm, :branch, cut.branch_id)

        bus_num = length(ref(pm, :bus))
        cont_bus_injection = Dict{Int,Any}()
        for (b,bus) in ref(pm, :bus)
            inj = 0.0
            for g in ref(pm, :bus_gens, b)
                if g != cut.gen_id
                    inj += var(pm, :pg, g)
                end
            end
            cont_bus_injection[b] = inj + var(pm, :pg, cut.gen_id)/bus_num
        end

        #rate = branch["rate_a"]
        rate = branch["rate_c"]
        @constraint(pm.model,  sum( weight*(cont_bus_injection[bus_id] - bus_withdrawal[bus_id]) for (bus_id, weight) in cut.bus_injection) <= rate + var(pm, :gen_cont_flow_vio, i))
        @constraint(pm.model, -sum( weight*(cont_bus_injection[bus_id] - bus_withdrawal[bus_id]) for (bus_id, weight) in cut.bus_injection) <= rate + var(pm, :gen_cont_flow_vio, i))
        #@constraint(pm.model,  sum( weight*(cont_bus_injection[bus_id] - bus_withdrawal[bus_id]) for (bus_id, weight) in cut.bus_injection) <= rate)
        #@constraint(pm.model, -sum( weight*(cont_bus_injection[bus_id] - bus_withdrawal[bus_id]) for (bus_id, weight) in cut.bus_injection) <= rate)
    end

    ##### Setup Objective #####
    objective_variable_pg_cost(pm)
    # explicit network id needed because of conductor-less
    pg_cost = var(pm, pm.cnw, :pg_cost)
    branch_cont_flow_vio = var(pm, pm.cnw, :branch_cont_flow_vio)
    gen_cont_flow_vio = var(pm, pm.cnw, :gen_cont_flow_vio)

    @objective(pm.model, Min,
        sum( pg_cost[i] for (i,gen) in ref(pm, :gen) ) +
        sum( 5e5*branch_cont_flow_vio[i] for i in 1:length(ref(pm, :branch_flow_cuts)) ) +
        sum( 5e5*gen_cont_flow_vio[i] for i in 1:length(ref(pm, :gen_flow_cuts)) )
    )
end


text:
module RBM
using Winston

export dream
export showWeights

function xp(x)
    return 1.0 ./ (1+exp(-x))
end

function vtohprob{T<:Real,P<:Integer}(rbmw::Matrix{T},vstate::Matrix{P})
    return xp(rbmw * vstate)
end

function htovprob{T<:Real,P<:Integer}(rbmw::Matrix{T},hstate::Matrix{P})
    return xp(rbmw' * hstate)
end

function gradient{T<:Real,P<:Integer}(vstate::Matrix{P},hstate::Matrix{T})
    hstate * vstate'/size(vstate)[2]
end

function bernoulli{T<:Real}(v::Vector{T})
    return map((x)->x>rand()?1:0,v)
end

function bernoulli{T<:Real}(m::Matrix{T})
    return map((x)->x>rand()?1:0,m)
end

function cd{T<:Real}(rbmw::Matrix{T},data::Matrix{T})
    vdata=bernoulli(data)
    h0=bernoulli(vtohprob(rbmw,vdata))
    vh0=gradient(vdata,h0)
    v1=bernoulli(htovprob(rbmw,h0))
    h1=bernoulli(vtohprob(rbmw,v1))
    vh1=gradient(v1,h1)
    return vh0-vh1
end

function showWeights{T<:Real}(W::Matrix{T},row,col)
    r,c = size(W)
    plots = []
    for i in 1:r
        plots=vcat(plots,Winston.imagesc((reshape(W[i,:],28,28))'))
    end

    mgrid=Winston.Table(row,col)
    c=1
    for i = 1:row
        for j=1:col
            mgrid[i,j] = plots[c]
            c+=1
        end
    end
    display(mgrid)
end

function dream{T<:Real}(rbmw::Matrix{T},data::Matrix{T})
    vdata=bernoulli(data)
    h0=bernoulli(vtohprob(rbmw,vdata))
    v1=bernoulli(htovprob(rbmw,h0))

    rw = size(vdata)[2]
    mgrid = Winston.Table(1,2)
    for i in 1:rw
        p1=Winston.imagesc(reshape(vdata[:,i],28,28)')
        p2=Winston.imagesc(reshape(v1[:,i],28,28)')
        mgrid[1,1]=p1
        mgrid[1,2]=p2
        display(mgrid)
    end
    return v1
end

function rbm{T<:Real,P<:Integer}(nhid::P,data::Matrix{T},lr::T,niter::P,mbatchsz::P,mom::T)
    psz,nsz=size(data)
    momentum=mom
    if nsz % mbatchsz != 0
        println("test cases must be multiple of batch size")
        return
    end
    model=(rand(nhid,psz)*2-1)*0.1
    momSpeed=zeros(nhid,psz)

    startMini=1
    for iter in 1:niter
        if mod(iter,50)==0
            println("iter:",iter," ",sqrt(mean(map((x)->x*x,model))))
        end
        mbatch=data[:,startMini:(startMini+mbatchsz-1)]
        startMini=(startMini+mbatchsz) % nsz
        gradient=cd(model,mbatch)
        momSpeed=momentum*momSpeed+gradient
        model=model+momSpeed*lr
    end
    return model
end
end

using RBM
#trainData=readdlm("./mnist_train.csv",',',Float64)
#trainLab=trainData[:,1]
#trainData=trainData[:,2:size(trainData)[2]]/255.0
#trainData=trainData'

testData=readdlm("./mnist_test.csv",',',Float64)
testLab=testData[:,1]
testData=testData[:,2:size(testData)[2]]/255.0
testData=testData'

hidSz=30
mbatchsz=100
#niter=20000
niter = 500
lr=0.01
mom=0.9

#W=RBM.rbm(hidSz,trainData,lr,niter,mbatchsz,mom)

W=readcsv("./weighnobias.csv")

RBM.dream(W,testData[:,1:50]);


text:
# integration segment (a,b), estimated integral I, and estimated error E
struct Segment{TX,TI,TE}
    a::TX
    b::TX
    I::TI
    E::TE
end
Base.@pure Base.promote_rule(::Type{Segment{TX,TI,TE}}, ::Type{Segment{TX′,TI′,TE′}}) where {TX,TI,TE,TX′,TI′,TE′} =
    Segment{promote_type(TX,TX′), promote_type(TI,TI′), promote_type(TE,TE′)}
Base.convert(::Type{T}, s::Segment) where {T<:Segment} = T(s.a,s.b,s.I,s.E)
Base.isless(i::Segment, j::Segment) = isless(i.E, j.E)

# Internal routine: approximately integrate f(x) over the interval (a,b)
# by evaluating the integration rule (x,w,gw). Return a Segment.
function evalrule(f, a,b, x,w,gw, nrm)
    # Ik and Ig are integrals via Kronrod and Gauss rules, respectively
    s = convert(eltype(x), 0.5) * (b-a)
    n1 = 1 - (length(x) & 1) # 0 if even order, 1 if odd order
    # unroll first iterationof loop to get correct type of Ik and Ig
    fg = f(a + (1+x[2])*s) + f(a + (1-x[2])*s)
    fk = f(a + (1+x[1])*s) + f(a + (1-x[1])*s)
    Ig = fg * gw[1]
    Ik = fg * w[2] + fk * w[1]
    for i = 2:length(gw)-n1
        fg = f(a + (1+x[2i])*s) + f(a + (1-x[2i])*s)
        fk = f(a + (1+x[2i-1])*s) + f(a + (1-x[2i-1])*s)
        Ig += fg * gw[i]
        Ik += fg * w[2i] + fk * w[2i-1]
    end
    if n1 == 0 # even: Gauss rule does not include x == 0
        Ik += f(a + s) * w[end]
    else # odd: don't count x==0 twice in Gauss rule
        f0 = f(a + s)
        Ig += f0 * gw[end]
        Ik += f0 * w[end] +
            (f(a + (1+x[end-1])*s) + f(a + (1-x[end-1])*s)) * w[end-1]
    end
    Ik_s, Ig_s = Ik * s, Ig * s # new variable since this may change the type
    E = nrm(Ik_s - Ig_s)
    if isnan(E) || isinf(E)
        throw(DomainError(a+s, "integrand produced $E in the interval ($a, $b)"))
    end
    return Segment(oftype(s, a), oftype(s, b), Ik_s, E)
end


text:
using QuakeML, Test
import Dates

@testset "Parsing" begin
    @testset "Local parsing" begin
        @test QuakeML.local_parse(String, "xyzx̂") == "xyzx̂"
        @test QuakeML.local_parse(String, SubString("φabc")) == "φabc"
        @test QuakeML.local_parse(SubString, "ABC ♪") == "ABC ♪"
        @test QuakeML.local_parse(Float64, "-1.0") == -1.0 == parse(Float64, "-1.0")
        @test QuakeML.local_tryparse(Int, "12") == 12
        @test QuakeML.local_tryparse(Int, "1.23") === nothing
        @test QuakeML.local_tryparse(String, "⊫") == "⊫"
    end

    @testset "Dates" begin
        let
            # Too much precision
            @test QuakeML.readstring("""
                <?xml version="1.0" encoding="UTF-8"?>
                <quakeml xmlns="http://quakeml.org/xmlns/quakeml/1.0">
                  <eventParameters publicID="smi:TEST/test">
                    <creationInfo>
                      <agencyID>ISC</agencyID>
                      <creationTime>2020-03-11T11:45:58.1234567890</creationTime>
                    </creationInfo>
                  </eventParameters>
                </quakeml>
                """).creation_info.creation_time == Dates.DateTime(2020, 03, 11, 11, 45, 58, 123)

            # Time zone ahead of UTC and too much precision
            @test QuakeML.readstring("""
                <?xml version="1.0" encoding="UTF-8"?>
                <quakeml xmlns="http://quakeml.org/xmlns/quakeml/1.0">
                  <eventParameters publicID="smi:TEST/test">
                    <creationInfo>
                      <agencyID>ISC</agencyID>
                      <creationTime>2000-01-01T12:34:00.123456789-12:34</creationTime>
                    </creationInfo>
                  </eventParameters>
                </quakeml>
                """).creation_info.creation_time == Dates.DateTime(2000, 1, 1, 0, 0, 0, 123)
            # UTC specified
            @test QuakeML.readstring("""
                <?xml version="1.0" encoding="UTF-8"?>
                <quakeml xmlns="http://quakeml.org/xmlns/quakeml/1.0">
                  <eventParameters publicID="smi:TEST/test">
                    <creationInfo>
                      <agencyID>ISC</agencyID>
                      <creationTime>2020-03-11T11:45:58.1234567890Z</creationTime>
                    </creationInfo>
                  </eventParameters>
                </quakeml>
                """).creation_info.creation_time == Dates.DateTime(2020, 03, 11, 11, 45, 58, 123)

        end
    end

    @testset "Version" begin
        @test_logs (:warn,
          "document is StationXML version 1.3.0; only v1.2 data will be read") QuakeML.readstring("""
                <?xml version="1.0" encoding="UTF-8"?>
                <quakeml xmlns="http://quakeml.org/xmlns/quakeml/1.3">
                  <eventParameters publicID="smi:TEST/test">
                  </eventParameters>
                </quakeml>
                """)
    end
end

max_stars_repo_name: anowacki/QuakeML


text:
"""
GenX: An Configurable Capacity Expansion Model
Copyright (C) 2021,  Massachusetts Institute of Technology
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
A complete copy of the GNU General Public License v2 (GPLv2) is available
in LICENSE.txt.  Users uncompressing this from an archive may not have
received this license file.  If not, see <http://www.gnu.org/licenses/>.
"""

@doc raw"""
	investment_charge(EP::Model, inputs::Dict)

This function defines the expressions and constraints keeping track of total available storage charge capacity across all resources as well as constraints on capacity retirements. The function also adds investment and fixed O\&M related costs related to charge capacity to the objective function.

The total capacity of each resource is defined as the sum of the existing capacity plus the newly invested capacity minus any retired capacity.

```math
\begin{aligned}
& \Delta^{total,charge}_{y,z} =(\overline{\Delta^{charge}_{y,z}}+\Omega^{charge}_{y,z}-\Delta^{charge}_{y,z}) \forall y \in \mathcal{O}^{asym}, z \in \mathcal{Z}
\end{aligned}
```

One cannot retire more capacity than existing capacity.
```math
\begin{aligned}
&\Delta^{charge}_{y,z} \leq \overline{\Delta^{charge}_{y,z}}
	\hspace{4 cm}  \forall y \in \mathcal{O}^{asym}, z \in \mathcal{Z}
\end{aligned}
```

For resources where $\overline{\Omega_{y,z}^{charge}}$ and $\underline{\Omega_{y,z}^{charge}}$ is defined, then we impose constraints on minimum and maximum power capacity.
```math
\begin{aligned}
& \Delta^{total,charge}_{y,z} \leq \overline{\Omega}^{charge}_{y,z}
	\hspace{4 cm}  \forall y \in \mathcal{O}^{asym}, z \in \mathcal{Z} \\
& \Delta^{total,charge}_{y,z}  \geq \underline{\Omega}^{charge}_{y,z}
	\hspace{4 cm}  \forall y \in \mathcal{O}^{asym}, z \in \mathcal{Z}
\end{aligned}
```

In addition, this function adds investment and fixed O&M related costs related to charge capacity to the objective function:
```math
\begin{aligned}
& 	\sum_{y \in \mathcal{O}^{asym} } \sum_{z \in \mathcal{Z}}
	\left( (\pi^{INVEST,charge}_{y,z} \times    \Omega^{charge}_{y,z})
	+ (\pi^{FOM,charge}_{y,z} \times  \Delta^{total,charge}_{y,z})\right)
\end{aligned}
```
"""
function investment_charge(EP::Model, inputs::Dict)

	println("Charge Investment Module")

	dfGen = inputs["dfGen"]

	G = inputs["G"]     # Number of resources (generators, storage, DR, and DERs)

	STOR_ASYMMETRIC = inputs["STOR_ASYMMETRIC"] # Set of storage resources with asymmetric (separte) charge/discharge capacity components

	NEW_CAP_CHARGE = inputs["NEW_CAP_CHARGE"] # Set of asymmetric charge/discharge storage resources eligible for new charge capacity
	RET_CAP_CHARGE = inputs["RET_CAP_CHARGE"] # Set of asymmetric charge/discharge storage resources eligible for charge capacity retirements

	### Variables ###

	## Storage capacity built and retired for storage resources with independent charge and discharge power capacities (STOR=2)

	# New installed charge capacity of resource "y"
	@variable(EP, vCAPCHARGE[y in NEW_CAP_CHARGE] >= 0)

	# Retired charge capacity of resource "y" from existing capacity
	@variable(EP, vRETCAPCHARGE[y in RET_CAP_CHARGE] >= 0)

	### Expressions ###

	@expression(EP, eTotalCapCharge[y in STOR_ASYMMETRIC],
		if (y in intersect(NEW_CAP_CHARGE, RET_CAP_CHARGE))
			dfGen[!,:Existing_Charge_Cap_MW][y] + EP[:vCAPCHARGE][y] - EP[:vRETCAPCHARGE][y]
		elseif (y in setdiff(NEW_CAP_CHARGE, RET_CAP_CHARGE))
			dfGen[!,:Existing_Charge_Cap_MW][y] + EP[:vCAPCHARGE][y]
		elseif (y in setdiff(RET_CAP_CHARGE, NEW_CAP_CHARGE))
			dfGen[!,:Existing_Charge_Cap_MW][y] - EP[:vRETCAPCHARGE][y]
		else
			dfGen[!,:Existing_Charge_Cap_MW][y] + EP[:vZERO]
		end
	)

	## Objective Function Expressions ##

	# Fixed costs for resource "y" = annuitized investment cost plus fixed O&M costs
	# If resource is not eligible for new charge capacity, fixed costs are only O&M costs
	@expression(EP, eCFixCharge[y in STOR_ASYMMETRIC],
		if y in NEW_CAP_CHARGE # Resources eligible for new charge capacity
			dfGen[!,:Inv_Cost_Charge_per_MWyr][y]*vCAPCHARGE[y] + dfGen[!,:Fixed_OM_Cost_Charge_per_MWyr][y]*eTotalCapCharge[y]
		else
			dfGen[!,:Fixed_OM_Cost_Charge_per_MWyr][y]*eTotalCapCharge[y]
		end
	)

	# Sum individual resource contributions to fixed costs to get total fixed costs
	@expression(EP, eTotalCFixCharge, sum(EP[:eCFixCharge][y] for y in STOR_ASYMMETRIC))

	# Add term to objective function expression
	EP[:eObj] += eTotalCFixCharge

	### Constratints ###

	## Constraints on retirements and capacity additions
	#Cannot retire more charge capacity than existing charge capacity
 	@constraint(EP, cMaxRetCharge[y in RET_CAP_CHARGE], vRETCAPCHARGE[y] <= dfGen[!,:Existing_Charge_Cap_MW][y])

  	#Constraints on new built capacity

	# Constraint on maximum charge capacity (if applicable) [set input to -1 if no constraint on maximum charge capacity]
	# DEV NOTE: This constraint may be violated in some cases where Existing_Charge_Cap_MW is >= Max_Charge_Cap_MWh and lead to infeasabilty
	@constraint(EP, cMaxCapCharge[y in intersect(dfGen[!,:Max_Charge_Cap_MW].>0, STOR_ASYMMETRIC)], eTotalCapCharge[y] <= dfGen[!,:Max_Charge_Cap_MW][y])

	# Constraint on minimum charge capacity (if applicable) [set input to -1 if no constraint on minimum charge capacity]
	# DEV NOTE: This constraint may be violated in some cases where Existing_Charge_Cap_MW is <= Min_Charge_Cap_MWh and lead to infeasabilty
	@constraint(EP, cMinCapCharge[y in intersect(dfGen[!,:Min_Charge_Cap_MW].>0, STOR_ASYMMETRIC)], eTotalCapCharge[y] >= dfGen[!,:Min_Charge_Cap_MW][y])

	return EP
end

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
using Documenter, LibCEED, LinearAlgebra

makedocs(sitename="LibCEED.jl Docs",
         format=Documenter.HTML(prettyurls=false),
         pages=[
             "Home" => "index.md",
             "Ceed Objects" => [
                "Ceed.md",
                "CeedVector.md",
                "ElemRestriction.md",
                "Basis.md",
                "QFunction.md",
                "Operator.md",
             ],
             "Utilities" => [
                "Misc.md",
                "Globals.md",
                "Quadrature.md",
             ],
             "LibCEED.md",
             "C.md",
             "UserQFunctions.md",
             "Examples.md",
         ])

deploydocs(
    repo="github.com/CEED/libCEED-julia-docs.git",
    devbranch="main",
    push_preview=true,
)

max_stars_repo_name: AdelekeBankole/libCEED

max_issues_repo_name: AdelekeBankole/libCEED

max_forks_repo_name: AdelekeBankole/libCEED


text:
@testset "DC Power Flow Models Monitored Line Flow Constraints and Static Unbounded" begin
    system = PSB.build_system(PSITestSystems, "c_sys5_ml")
    limits = PSY.get_flow_limits(PSY.get_component(MonitoredLine, system, "1"))
    for model in [DCPPowerModel, StandardPTDFModel]
        template = get_thermal_dispatch_template_network(
            NetworkModel(model; PTDF=PSY.PTDF(system)),
        )
        model_m = DecisionModel(template, system; optimizer=HiGHS_optimizer)
        @test build!(model_m; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT
        @test check_variable_bounded(model_m, FlowActivePowerVariable, MonitoredLine)
        @test check_variable_unbounded(model_m, FlowActivePowerVariable, Line)

        @test solve!(model_m) == RunStatus.SUCCESSFUL
        @test check_flow_variable_values(
            model_m,
            FlowActivePowerVariable,
            MonitoredLine,
            "1",
            limits.from_to,
        )
    end
end

@testset "DC Power Flow Models Monitored Line Flow Constraints and Static with Bounds" begin
    system = PSB.build_system(PSITestSystems, "c_sys5_ml")
    set_rate!(PSY.get_component(Line, system, "2"), 1.5)
    for model in [DCPPowerModel, StandardPTDFModel]
        template = get_thermal_dispatch_template_network(
            NetworkModel(model; PTDF=PSY.PTDF(system)),
        )
        set_device_model!(template, DeviceModel(Line, StaticBranch))
        set_device_model!(template, DeviceModel(MonitoredLine, StaticBranchUnbounded))
        model_m = DecisionModel(template, system; optimizer=HiGHS_optimizer)
        @test build!(model_m; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT

        @test check_variable_unbounded(model_m, FlowActivePowerVariable, MonitoredLine)
        # Broken
        # @test check_variable_bounded(model_m, FlowActivePowerVariable, Line)

        @test solve!(model_m) == RunStatus.SUCCESSFUL
        @test check_flow_variable_values(model_m, FlowActivePowerVariable, Line, "2", 1.5)
    end
end

@testset "AC Power Flow Monitored Line Flow Constraints" begin
    system = PSB.build_system(PSITestSystems, "c_sys5_ml")
    limits = PSY.get_flow_limits(PSY.get_component(MonitoredLine, system, "1"))
    template = get_thermal_dispatch_template_network(ACPPowerModel)
    model_m = DecisionModel(template, system; optimizer=ipopt_optimizer)
    @test build!(model_m; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT

    @test check_variable_bounded(model_m, FlowActivePowerFromToVariable, MonitoredLine)
    @test check_variable_unbounded(model_m, FlowReactivePowerFromToVariable, MonitoredLine)

    @test solve!(model_m) == RunStatus.SUCCESSFUL
    @test check_flow_variable_values(
        model_m,
        FlowActivePowerFromToVariable,
        FlowReactivePowerFromToVariable,
        MonitoredLine,
        "1",
        0.0,
        limits.from_to,
    )
end

@testset "DC Power Flow Models for HVDCLine with with Line Flow Constraints, TapTransformer & Transformer2W Unbounded" begin
    ratelimit_constraint_keys = [
        PSI.ConstraintKey(RateLimitConstraint, Transformer2W, "ub"),
        PSI.ConstraintKey(RateLimitConstraint, Transformer2W, "lb"),
        PSI.ConstraintKey(RateLimitConstraint, TapTransformer, "ub"),
        PSI.ConstraintKey(RateLimitConstraint, TapTransformer, "lb"),
    ]

    system = PSB.build_system(PSITestSystems, "c_sys14_dc")
    hvdc_line = PSY.get_component(HVDCLine, system, "DCLine3")
    limits_from = PSY.get_active_power_limits_from(hvdc_line)
    limits_to = PSY.get_active_power_limits_to(hvdc_line)
    limits_min = min(limits_from.min, limits_to.min)
    limits_max = min(limits_from.max, limits_to.max)

    tap_transformer = PSY.get_component(TapTransformer, system, "Trans3")
    rate_limit = PSY.get_rate(tap_transformer)

    transformer = PSY.get_component(Transformer2W, system, "Trans4")
    rate_limit2w = PSY.get_rate(tap_transformer)

    for model in [DCPPowerModel, StandardPTDFModel],
        hvdc_model in [HVDCDispatch, HVDCLossless]

        template =
            get_template_dispatch_with_network(NetworkModel(model; PTDF=PSY.PTDF(system)))
        set_device_model!(template, HVDCLine, hvdc_model)
        model_m = DecisionModel(template, system; optimizer=ipopt_optimizer)
        @test build!(model_m; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT

        @test check_variable_bounded(model_m, FlowActivePowerVariable, HVDCLine)
        @test check_variable_unbounded(model_m, FlowActivePowerVariable, TapTransformer)
        @test check_variable_unbounded(model_m, FlowActivePowerVariable, Transformer2W)

        psi_constraint_test(model_m, ratelimit_constraint_keys)

        @test solve!(model_m) == RunStatus.SUCCESSFUL

        @test check_flow_variable_values(
            model_m,
            FlowActivePowerVariable,
            HVDCLine,
            "DCLine3",
            limits_min,
            limits_max,
        )
        @test check_flow_variable_values(
            model_m,
            FlowActivePowerVariable,
            TapTransformer,
            "Trans3",
            -rate_limit,
            rate_limit,
        )
        @test check_flow_variable_values(
            model_m,
            FlowActivePowerVariable,
            Transformer2W,
            "Trans4",
            -rate_limit2w,
            rate_limit2w,
        )
    end
end

@testset "DC Power Flow Models for Unbounded HVDCLine, and StaticBranchBounds for TapTransformer & Transformer2W" begin
    system = PSB.build_system(PSITestSystems, "c_sys14_dc")
    hvdc_line = PSY.get_component(HVDCLine, system, "DCLine3")
    limits_from = PSY.get_active_power_limits_from(hvdc_line)
    limits_to = PSY.get_active_power_limits_to(hvdc_line)
    limits_min = min(limits_from.min, limits_to.min)
    limits_max = min(limits_from.max, limits_to.max)

    tap_transformer = PSY.get_component(TapTransformer, system, "Trans3")
    rate_limit = PSY.get_rate(tap_transformer)

    transformer = PSY.get_component(Transformer2W, system, "Trans4")
    rate_limit2w = PSY.get_rate(tap_transformer)

    for model in [DCPPowerModel, StandardPTDFModel]
        template =
            get_template_dispatch_with_network(NetworkModel(model; PTDF=PSY.PTDF(system)))
        set_device_model!(template, DeviceModel(HVDCLine, HVDCUnbounded))
        set_device_model!(template, DeviceModel(TapTransformer, StaticBranchBounds))
        set_device_model!(template, DeviceModel(Transformer2W, StaticBranchBounds))
        model_m = DecisionModel(template, system; optimizer=ipopt_optimizer)
        @test build!(model_m; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT

        if model == DCPPowerModel
            @test check_variable_unbounded(model_m, FlowActivePowerVariable, HVDCLine)
            @test check_variable_bounded(model_m, FlowActivePowerVariable, TapTransformer)
            @test check_variable_bounded(model_m, FlowActivePowerVariable, TapTransformer)
        end

        @test solve!(model_m) == RunStatus.SUCCESSFUL

        @test check_flow_variable_values(
            model_m,
            FlowActivePowerVariable,
            HVDCLine,
            "DCLine3",
            limits_min,
            limits_max,
        )
        @test check_flow_variable_values(
            model_m,
            FlowActivePowerVariable,
            TapTransformer,
            "Trans3",
            -rate_limit,
            rate_limit,
        )
        @test check_flow_variable_values(
            model_m,
            FlowActivePowerVariable,
            Transformer2W,
            "Trans4",
            -rate_limit2w,
            rate_limit2w,
        )
    end
end

@testset "AC Power Flow Models for HVDCLine Flow Constraints and TapTransformer & Transformer2W Unbounded" begin
    ratelimit_constraint_keys = [
        PSI.ConstraintKey(RateLimitConstraintFromTo, Transformer2W),
        PSI.ConstraintKey(RateLimitConstraintToFrom, Transformer2W),
        PSI.ConstraintKey(RateLimitConstraintFromTo, TapTransformer),
        PSI.ConstraintKey(RateLimitConstraintToFrom, TapTransformer),
        PSI.ConstraintKey(FlowRateConstraintFromTo, HVDCLine),
        PSI.ConstraintKey(FlowRateConstraintToFrom, HVDCLine),
        PSI.ConstraintKey(HVDCPowerBalance, HVDCLine),
    ]

    system = PSB.build_system(PSITestSystems, "c_sys14_dc")

    hvdc_line = PSY.get_component(HVDCLine, system, "DCLine3")
    limits_from = PSY.get_active_power_limits_from(hvdc_line)
    limits_to = PSY.get_active_power_limits_to(hvdc_line)
    limits_min = min(limits_from.min, limits_to.min)
    limits_max = min(limits_from.max, limits_to.max)

    tap_transformer = PSY.get_component(TapTransformer, system, "Trans3")
    rate_limit = PSY.get_rate(tap_transformer)

    transformer = PSY.get_component(Transformer2W, system, "Trans4")
    rate_limit2w = PSY.get_rate(tap_transformer)

    template = get_template_dispatch_with_network(ACPPowerModel)
    set_device_model!(template, DeviceModel(HVDCLine, HVDCDispatch))
    model_m = DecisionModel(template, system; optimizer=ipopt_optimizer)
    @test build!(model_m; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT

    check_variable_bounded(model_m, FlowReactivePowerToFromVariable, HVDCLine)
    check_variable_bounded(model_m, FlowActivePowerVariable, HVDCLine)
    @test check_variable_bounded(model_m, FlowActivePowerFromToVariable, TapTransformer)
    @test check_variable_unbounded(model_m, FlowReactivePowerFromToVariable, TapTransformer)
    @test check_variable_bounded(model_m, FlowActivePowerToFromVariable, Transformer2W)
    @test check_variable_unbounded(model_m, FlowReactivePowerToFromVariable, Transformer2W)

    psi_constraint_test(model_m, ratelimit_constraint_keys)

    @test solve!(model_m) == RunStatus.SUCCESSFUL

    @test check_flow_variable_values(
        model_m,
        FlowActivePowerVariable,
        FlowReactivePowerToFromVariable,
        HVDCLine,
        "DCLine3",
        limits_max,
    )
    @test check_flow_variable_values(
        model_m,
        FlowActivePowerFromToVariable,
        FlowReactivePowerFromToVariable,
        TapTransformer,
        "Trans3",
        rate_limit,
    )
    @test check_flow_variable_values(
        model_m,
        FlowActivePowerToFromVariable,
        FlowReactivePowerToFromVariable,
        Transformer2W,
        "Trans4",
        rate_limit2w,
    )
end

@testset "DC Power Flow Models for HVDCLine Dispatch and TapTransformer & Transformer2W Unbounded" begin
    ratelimit_constraint_keys = [
        PSI.ConstraintKey(RateLimitConstraint, Transformer2W, "ub"),
        PSI.ConstraintKey(RateLimitConstraint, Line, "ub"),
        PSI.ConstraintKey(RateLimitConstraint, Line, "lb"),
        PSI.ConstraintKey(RateLimitConstraint, TapTransformer, "ub"),
        PSI.ConstraintKey(RateLimitConstraint, Transformer2W, "lb"),
        PSI.ConstraintKey(RateLimitConstraint, TapTransformer, "lb"),
        PSI.ConstraintKey(FlowRateConstraint, HVDCLine),
    ]

    system = PSB.build_system(PSITestSystems, "c_sys14_dc")

    hvdc_line = PSY.get_component(HVDCLine, system, "DCLine3")
    limits_from = PSY.get_active_power_limits_from(hvdc_line)
    limits_to = PSY.get_active_power_limits_to(hvdc_line)
    limits_min = min(limits_from.min, limits_to.min)
    limits_max = min(limits_from.max, limits_to.max)

    tap_transformer = PSY.get_component(TapTransformer, system, "Trans3")
    rate_limit = PSY.get_rate(tap_transformer)

    transformer = PSY.get_component(Transformer2W, system, "Trans4")
    rate_limit2w = PSY.get_rate(tap_transformer)

    template = get_template_dispatch_with_network(
        NetworkModel(StandardPTDFModel; PTDF=PSY.PTDF(system)),
    )
    set_device_model!(template, DeviceModel(HVDCLine, HVDCDispatch))
    model_m = DecisionModel(template, system; optimizer=ipopt_optimizer)
    @test build!(model_m; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT

    @test check_variable_bounded(model_m, FlowActivePowerVariable, HVDCLine)
    @test !check_variable_bounded(model_m, FlowActivePowerVariable, TapTransformer)
    @test !check_variable_bounded(model_m, FlowActivePowerVariable, Transformer2W)
    @test check_variable_unbounded(model_m, FlowActivePowerVariable, Line)

    psi_constraint_test(model_m, ratelimit_constraint_keys)

    @test solve!(model_m) == RunStatus.SUCCESSFUL

    @test check_flow_variable_values(
        model_m,
        FlowActivePowerVariable,
        HVDCLine,
        "DCLine3",
        limits_max,
    )
    @test check_flow_variable_values(
        model_m,
        FlowActivePowerVariable,
        TapTransformer,
        "Trans3",
        rate_limit,
    )
    @test check_flow_variable_values(
        model_m,
        FlowActivePowerVariable,
        Transformer2W,
        "Trans4",
        rate_limit2w,
    )
end


text:
__precompile__(true)
"""
API Tools package

Copyright 2018-2019 Gandalf Software, Inc., Scott P. Jones

Licensed under MIT License, see LICENSE.md

(@def macro "stolen" from DiffEqBase.jl/src/util.jl :-) )
"""
module ModuleInterfaceTools

const debug = Ref(false)
const showeval = Ref(false)

const V6_COMPAT = VERSION < v"0.7-"
const BIG_ENDIAN = (ENDIAN_BOM == 0x01020304)

_stdout() = stdout
_stderr() = stderr

Base.parse(::Type{Expr}, args...; kwargs...) =
    Meta.parse(args...; kwargs...)

export @api, V6_COMPAT, BIG_ENDIAN

_api_def(name, definition) =
    quote
        macro $(esc(name))()
            esc($(Expr(:quote, definition)))
        end
    end

const SymSet = Set{Symbol}

abstract type AbstractAPI end

struct TMP_API <: AbstractAPI
    mod::Module
    base::SymSet
    public::SymSet
    develop::SymSet
    public!::SymSet
    develop!::SymSet
    modules::SymSet

    TMP_API(mod::Module) = new(mod, SymSet(), SymSet(), SymSet(), SymSet(), SymSet(), SymSet())
end

const SymList = Tuple{Vararg{Symbol}}

struct API <: AbstractAPI
    mod::Module
    base::SymList
    public::SymList
    develop::SymList
    public!::SymList
    develop!::SymList
    modules::SymList
end

API(api::TMP_API) =
    API(api.mod, SymList(api.base), SymList(api.public), SymList(api.develop),
        SymList(api.public!), SymList(api.develop!), SymList(api.modules))

function Base.show(io::IO, api::AbstractAPI)
    println(io, "ModuleInterfaceTools.API: ", api.mod)
    for fld in (:base, :public, :develop, :public!, :develop!, :modules)
        syms = getfield(api, fld)
        isempty(syms) && continue
        print(fld, ":")
        for s in syms
            print(" ", s)
        end
        println()
        println()
    end
end

function m_eval(mod, expr)
    try
        showeval[] && println("m_eval($mod, $expr)")
        Core.eval(mod, expr)
    catch ex
        println("m_eval($mod, $expr)");
        println(sprint(showerror, ex, catch_backtrace()))
        #rethrow(ex)
    end
end

"""
@api <cmd> [<symbols>...]

 * freeze                # use at end of module to freeze API

 * list     <modules>... # list API(s) of given modules (or current if none given)

 * use      <modules>... # use, without importing (i.e. can't extend)
 * use!     <modules>... # use, without importing (i.e. can't extend), "export"
 * test     <modules>... # using public and develop APIs, for testing purposes
 * extend   <modules>... # for development, imports api & dev, use api & dev definitions
 * extend!  <modules>... # for development, imports api & dev, use api & dev definitions, "export"
 * reexport <modules>... # export public definitions from those modules

 * base     <names...>   # Add functions from Base that are part of the API (extendible)
 * base!    <names...>   # Add functions from Base or define them if not in Base
 * public   <names...>   # Add other symbols that are part of the public API (structs, consts)
 * public!  <names...>   # Add functions that are part of the public API (extendible)
 * develop  <names...>   # Add other symbols that are part of the development API
 * develop! <names...>   # Add functions that are part of the development API (extendible)
 * define!  <names...>   # Define functions to be extended, public API
 * defdev!  <names...>   # Define functions to be extended, develop API
 * modules  <names...>   # Add submodule names that are part of the API

 * path     <paths...>  # Add paths to LOAD_PATH

 * def <name> <expr>    # Same as the @def macro, creates a macro with the given name

"""
macro api(cmd::Symbol)
    mod = __module__
    cmd == :list   ? _api_list(mod) :
    cmd == :freeze ? _api_freeze(mod) :
    cmd == :test   ? _api_test(mod) :
    error("@api unrecognized command: $cmd")
end

function _api_display(mod, nam)
    if isdefined(mod, nam) && (api = m_eval(mod, nam)) !== nothing
        show(api);
    else
        println("Exported from $mod:")
        syms = names(mod)
        if !isempty(syms)
            print(fld, ":")
            for s in syms
                print(" ", s)
            end
        end
    end
    println()
end

_api_list(mod::Module) = (_api_display(mod, :__api__) ; _api_display(mod, :__tmp_api__))

function _api_freeze(mod::Module)
    ex = :( const __api__ = ModuleInterfaceTools.API(__tmp_api__) ; __tmp_api__ = nothing )
    isdefined(mod, :__tmp_api__) && m_eval(mod, :( __tmp_api__ !== nothing ) ) && m_eval(mod, ex)
    nothing
end

function _api_path(curmod, exprs)
    for exp in exprs
        if isa(exp, Expr) || isa(exp, Symbol)
            str = m_eval(curmod, exp)
        elseif isa(exp, String)
            str = exp
        else
            error("@api path: syntax error $exp")
        end
        m_eval(curmod, :( push!(LOAD_PATH, $str) ))
    end
    nothing
end

const _cmduse = (:use, :use!, :test, :extend, :extend!, :reexport, :list)
const _cmdadd =
    (:modules, :public, :develop, :public!, :develop!, :base, :base!, :define!, :defdev!)

_ff(lst, val) = (ret = findfirst(isequal(val), lst); ret === nothing ? 0 : ret)

function _add_def!(curmod, grp, exp)
    debug[] && print("_add_def!($curmod, $grp, $exp::$(typeof(exp))")
    if isa(exp, Symbol)
        sym = exp
    elseif isa(exp, AbstractString)
        sym = Symbol(exp)
    else
        error("@api $grp: syntax error $exp")
    end
    if grp == :base! && isdefined(Base, sym)
        m_eval(curmod, :(import Base.$sym ))
        m_eval(curmod, :(push!(__tmp_api__.base, $(QuoteNode(sym)))))
        return
    end
    m_eval(curmod, :(function $sym end))
    m_eval(curmod, (grp == :defdev!
                    ? :(push!(__tmp_api__.develop!, $(QuoteNode(sym))))
                    : :(push!(__tmp_api__.public!,  $(QuoteNode(sym))))))
end

function push_args!(symbols, lst, grp)
    for ex in lst
        if isa(ex, Expr) && ex.head == :tuple
            push_args!(symbols, ex.args)
        elseif isa(ex, Symbol)
            push!(symbols, ex)
        elseif isa(ex, AbstractString)
            push!(symbols, Symbol(ex))
        else
            error("@api $grp: syntax error $ex")
        end
    end
end

"""Initialize the temp api variable for this module"""
_init_api(curmod) =
    isdefined(curmod, :__tmp_api__) ||
        m_eval(curmod, :( global __tmp_api__ = ModuleInterfaceTools.TMP_API($curmod)))

"""Add symbols"""
function _add_symbols(curmod, grp, exprs)
    if debug[]
        print("_add_symbols($curmod, $grp, $exprs)")
        isdefined(curmod, :__tmp_api__) && print(" => ", m_eval(curmod, :__tmp_api__))
        println()
    end
    _init_api(curmod)
    if grp == :base! || grp == :define! || grp == :defdev!
        for ex in exprs
            if isa(ex, Expr) && ex.head == :tuple
                for sym in ex.args
                    _add_def!(curmod, grp, sym)
                end
            else
                _add_def!(curmod, grp, ex)
            end
        end
    else
        symbols = SymSet()
        for ex in exprs
            if isa(ex, Expr) && ex.head == :tuple
                push_args!(symbols, ex.args, grp)
            elseif isa(ex, Symbol)
                push!(symbols, ex)
            elseif isa(ex, AbstractString)
                push!(symbols, Symbol(ex))
            else
                error("@api $grp: syntax error $ex")
            end
        end
        if grp == :base
            for sym in symbols
                m_eval(curmod, :( import Base.$sym ))
            end
        end
        for sym in symbols
            m_eval(curmod, :( push!(__tmp_api__.$grp, $(QuoteNode(sym)) )))
        end
    end
    debug[] && println("after add symbols: ", m_eval(curmod, :__tmp_api__))
    nothing
end

has_api(mod) = isdefined(mod, :__api__)
get_api(curmod, mod) = m_eval(curmod, :( $mod.__api__ ))

function _api_extend(curmod, modules, cpy::Bool)
    for nam in modules
        mod = m_eval(curmod, nam)
        if has_api(mod)
            api = get_api(curmod, mod)
            _do_list(curmod, cpy, :import, Base, :Base, :base,  api)
            _do_list(curmod, cpy, :import, mod, nam, :public!,  api)
            _do_list(curmod, cpy, :import, mod, nam, :develop!, api)
            _do_list(curmod, cpy, :using,  mod, nam, :public,   api)
            _do_list(curmod, cpy, :using,  mod, nam, :develop,  api)
        else
            _do_list(curmod, cpy, :import, mod, nam, :public!,  names(mod))
        end
    end
    nothing
end

function _api_use(curmod, modules, cpy::Bool)
    for nam in modules
        mod = m_eval(curmod, nam)
        if has_api(mod)
            api = get_api(curmod, mod)
            _do_list(curmod, cpy, :using, mod, nam, :public,  api)
            _do_list(curmod, cpy, :using, mod, nam, :public!, api)
        else
            _do_list(curmod, cpy, :using, mod, nam, :public!, names(mod))
        end
    end
    nothing
end

function _api_reexport(curmod, modules)
    for nam in modules
        mod = m_eval(curmod, nam)
        if has_api(mod)
            api = get_api(curmod, mod)
            l1, l2, l3 = getfield(api, :modules), getfield(api, :public), getfield(api, :public!)
            if !(isempty(l1) && isempty(l2) && isempty(l3))
                debug[] && println("Reexport: api = $api:$l1,$l2,$l3")
                m_eval(curmod, Expr( :export, l1..., l2..., l3... ))
            end
        end
    end
    nothing
end

function _api_list(curmod, modules)
    for nam in modules
        _api_list(m_eval(curmod, nam))
    end
    nothing
end

_api_test(mod) = m_eval(mod, :(using Test))

function _api(curmod::Module, cmd::Symbol, exprs)
    cmd == :def && return _api_def(exprs...)
    cmd == :path && return _api_path(curmod, exprs)

    ind = _ff(_cmdadd, cmd)
    ind == 0 || return _add_symbols(curmod, cmd, exprs)

    _ff(_cmduse, cmd) == 0 && error("Syntax error: @api $cmd $exprs")

    debug[] && print("_api($curmod, $cmd, $exprs)")

    # Be nice and set up standard Test
    cmd == :test && _api_test(curmod)

    modules = SymSet()
    for ex in exprs
        if isa(ex, Expr) && ex.head == :tuple
            # Some of these might not just be modules
            # might have module(symbols, !syms, sym => other), need to add support for that
            for sym in ex.args
                if isa(sym, Symbol)
                    push!(modules, sym)
                    m_eval(curmod, :(import $sym))
                else
                    println("Not a symbol: $sym");
                    dump(sym);
                end
            end
        elseif isa(ex, Symbol)
            push!(modules, ex)
            m_eval(curmod, :(import $ex))
        else
            error("@api $cmd: syntax error $ex")
        end
    end
    debug[] && println(" => $modules")

    cmd == :reexport && return _api_reexport(curmod, modules)
    cmd == :list && return _api_list(curmod, modules)

    cpy = (cmd == :use!) || (cmd == :extend!)
    cpy && _init_api(curmod)

    for nam in modules
        mod = m_eval(curmod, nam)
        if has_api(mod)
            api = get_api(curmod, mod)
            for sym in getfield(api, :modules)
                if isdefined(mod, sym)
                    m_eval(curmod, :(using $nam.$sym))
                    cpy && m_eval(curmod, :( push!(__tmp_api__.modules, $(QuoteNode(sym)) )))
                else
                    println(_stderr(), "Warning: Exported symbol $sym is not defined in $nam")
                end
            end
        end
    end

    ((cmd == :use || cmd == :use!)
     ? _api_use(curmod, modules, cpy)
     : _api_extend(curmod, modules, cpy))
end

function makecmd(cmd, nam, sym)
    (cmd == :using
     ? Expr(cmd, Expr(:(:), Expr(:., nam), Expr(:., sym)))
     : Expr(cmd, Expr(:., nam, sym)))
end

_do_list(curmod, cpy, cmd, mod, nam, grp, api::API) =
    _do_list(curmod, cpy, cmd, mod, nam, grp, getfield(api, grp))

function _do_list(curmod, cpy, cmd, mod, nam, grp, lst)
    debug[] && println("_do_list($curmod, $cpy, $cmd, $mod, $nam, $grp, $lst)")
    for sym in lst
        if isdefined(mod, sym)
            m_eval(curmod, makecmd(cmd, nam, sym))
            cpy && m_eval(curmod, :( push!(__tmp_api__.$grp, $(QuoteNode(sym)) )))
        else
            println(_stderr(), "Warning: Exported symbol $sym is not defined in $nam")
        end
    end
end

macro api(cmd::Symbol, exprs...)
    _api(__module__, cmd, exprs)
end

end # module ModuleInterfaceTools


text:
using Test
using LinearAlgebra # dot

# using Einsum
using Tullio: @einsum

@testset "Test that vars in Main aren't overwritten by einsum" begin
    i = -1
    y = randn(10)
    @einsum x[i] := y[i]
    @test i == -1
end

@testset "Test that B is overwritten by := operator" begin
    B = randn(10, 10)
    A = randn(5, 10)
    @einsum B[i, j] := A[i, j] # this should run without a problem
    @test size(B) == size(A)
end

@testset "CP decomposition test case" begin

    # preallocated test case
    A = zeros(5, 6, 7)
    B = similar(A)
    C = similar(A)

    X = randn(5, 2)
    Y = randn(6, 2)
    Z = randn(7, 2)

    @einsum A[i, j, k] = X[i, r] * Y[j, r] * Z[k, r]
    # @einsimd B[i, j, k] = X[i, r] * Y[j, r] * Z[k, r]
    # @vielsum C[i, j, k] = X[i, r] * Y[j, r] * Z[k, r]

    for i = 1:5
        for j = 1:6
            for k = 1:7
                s = 0.0
                for r = 1:2
                    s += X[i, r] * Y[j, r] * Z[k, r]
                end
                @test isapprox(A[i, j, k], s)
                # @test isapprox(B[i, j, k], s)
                # @test isapprox(C[i, j, k], s)
            end
        end
    end

    # without preallocation
    @einsum A2[i, j, k] := X[i, r] * Y[j, r] * Z[k, r]
    @test isapprox(A, A2)

end

@testset "Interesting test case, can throw an error that local vars are declared twice." begin
    A = zeros(5, 6, 7)
    X = randn(5, 2)
    Y = randn(6, 2)
    Z = randn(7, 2)
    if true
        @einsum A[i, j, k] = X[i, r] * Y[j, r] * Z[k, r]
    else
        @einsum A[i, j, k] = X[i, r] * Y[j, r] * Z[k, r]
    end
end

@testset "From #21: local `T` does not interfer with internal T" begin
    function test(x::Vector{T}, y::Vector{T}) where T
        @einsum z := x[i] * y[i]
        return z
    end
    @test_nowarn test(rand(3), rand(3))
end

@testset "From #20: local `s` does not interfere with internal s" begin
    x = rand(2, 3)
    @test_nowarn @einsum y[i] := x[i, s]
end

@testset "At one point this threw an error because the lhs had no indices/arguments" begin
    x = randn(10)
    y = randn(10)
    @einsum k := x[i] * y[i]
    @test isapprox(k, dot(x, y))
end

@testset "Elementwise multiplication (this should create nested loops with no no summation.)" begin
    x = randn(10)
    y = randn(10)
    @einsum k[i] := x[i] * y[i]
    # @einsimd k2[i] := x[i] * y[i]
    # @vielsum k3[i] := x[i] * y[i]
    @test isapprox(k, x .* y)
    # @test isapprox(k2, x .* y)
    # @test isapprox(k3, x .* y)
end

@testset "Transpose a block matrix" begin
    z = [rand(2, 2) for i = 1:2, j = 1:2]
    @einsum t[i, j] := transpose(z[j, i])
    @test isapprox(z[1, 1],  t[1, 1]')
    @test isapprox(z[2, 2],  t[2, 2]')
    @test isapprox(z[1, 2],  t[2, 1]')
    @test isapprox(z[2, 1],  t[1, 2]')
end

@testset "Mapping functions" begin
    A = randn(10, 10)
    @einsum B[i, j] := exp(A[i, j])
    @test isapprox(exp.(A), B)
end

@testset "Example from numpy" begin
    A = reshape(collect(1:25), 5, 5)
    @einsum B[i] := A[i, i]
    @test all(B .== [1, 7, 13, 19, 25])

    # @vielsum C[i] := A[i, i]
    # @test all(C .== [1, 7, 13, 19, 25])

end

@testset "Adding a scalar, to-done" begin

    A = collect(reshape(1:12, 3,4))
    @einsum A[i, j] = A[i,j] + 50
    @test A == collect(50 .+ reshape(1:12, 3,4))

end

@testset "Test in-place operations" begin
    A = randn(5, 6, 7)
    B = randn(5, 6, 7)
    A1 = copy(A)
    B1 = copy(B)

    X = randn(5, 2)
    Y = randn(6, 2)
    Z = randn(7, 2)

    @einsum A[i, j, k] += X[i, r] * Y[j, r] * Z[k, r]
    # @einsimd B[i, j, k] += X[i, r] * Y[j, r] * Z[k, r]

    for i = 1:5
        for j = 1:6
            for k = 1:7
                s = 0.0
                for r = 1:2
                    s += X[i, r] * Y[j, r] * Z[k, r]
                end
                @test isapprox(A[i, j, k], A1[i, j, k] + s)
                # @test isapprox(B[i, j, k], B1[i, j, k] + s)
            end
        end
    end

end

@testset "scalar += dot" begin

    x = randn(10)
    y = randn(10)
    k0 = randn()
    k = k0
    @einsum k += x[i] * y[i]
    @test isapprox(k, k0 + dot(x, y))

end

@testset "test *= operator" begin

    A = randn(5, 6, 7)
    B = randn(5, 6, 7)
    A1 = copy(A)
    B1 = copy(B)

    X = randn(5, 2)
    Y = randn(6, 2)
    Z = randn(7, 2)

    @einsum A[i, j, k]  *= X[i, r] * Y[j, r] * Z[k, r]
    # @einsimd B[i, j, k] *= X[i, r] * Y[j, r] * Z[k, r]

    for i = 1:5
        for j = 1:6
            for k = 1:7
                s = 0.0
                for r = 1:2
                    s += X[i, r] * Y[j, r] * Z[k, r]
                end
                @test isapprox(A[i, j, k], A1[i, j, k] * s)
                # @test isapprox(B[i, j, k], B1[i, j, k] * s)
            end
        end
    end

    x = randn(10)
    y = randn(10)
    k0 = randn()
    k = k0
    @einsum k *= x[i] * y[i]
    @test isapprox(k, k0 * dot(x, y))
end

@testset "Test offsets" begin
    X = randn(10)

    # without preallocation
    @einsum A[i] := X[i + 5]
    @test_broken size(A) == (5,)  # here Tullio returns an OffsetArray
    @test_broken all(A .== X[6:end])

    # with preallocation
    B = zeros(10)
    @einsum B[i] = X[i + 5]
    @test size(B) == (10,)
    @test all(B[1:5] .== X[6:end])
end

@testset "Test symbolic offsets" begin
    offset = 5
    X = randn(10)

    # without preallocation
    # @einsum A[i] := X[i + :offset] # error on 1.0
    @einsum A[i] := X[i + $offset]  # here Tullio returns an OffsetArray
    @test_broken size(A) == (5,)
    @test_broken all(A .== X[6:end])

    # with preallocation
    B = zeros(10)
    # @einsum B[i] = X[i + :offset] # error on 1.0
    @einsum B[i] = X[i + $offset]
    @test size(B) == (10,)
    @test all(B[1:5] .== X[6:end])
end

@testset "Test adding/subtracting constants" begin
    k = 5
    X = randn(10)

    # without preallocation
    @einsum A[i] := X[i] + k
    @einsum B[i] := X[i] - k
    @test isapprox(A, X .+ k)
    @test isapprox(B, X .- k)

    @einsum A[i] := X[i] + $k  # Tullio prefers $k, it becomes a function argument
    @einsum B[i] := X[i] - $k
    @test isapprox(A, X .+ k)
    @test isapprox(B, X .- k)

    # with preallocation
    C, D = zeros(10), zeros(10)
    @einsum C[i] = X[i] + $k
    @einsum D[i] = X[i] - $k
    @test isapprox(C, X .+ k)
    @test isapprox(D, X .- k)
end

@testset "Test multiplying/dividing constants" begin
    k = 5
    X = randn(10)

    # without preallocation
    @einsum A[i] := X[i] * k
    @einsum B[i] := X[i] / k
    @test isapprox(A, X .* k)
    @test isapprox(B, X ./ k)

    @einsum A[i] := X[i] * $k  # Tullio prefers $k, it becomes a function argument
    @einsum B[i] := X[i] / $k
    @test isapprox(A, X .* k)
    @test isapprox(B, X ./ k)

    # with preallocation
    C, D = zeros(10), zeros(10)
    @einsum C[i] = X[i] * k
    @einsum D[i] = X[i] / k
    @test isapprox(C, X .* k)
    @test isapprox(D, X ./ k)

    @einsum C[i] = X[i] * $k  # Tullio prefers $k, it becomes a function argument
    @einsum D[i] = X[i] / $k
    @test isapprox(C, X .* k)
    @test isapprox(D, X ./ k)
end

@testset "Test indexing with a constant" begin
    A = randn(10, 2)
    j = 2
    # @einsum B[i] := A[i, :j] # error on Julia 1.0, i.e. this broke at some point in Einsum.jl
    @einsum B[i] := A[i, $j]   # Tullio's notation, here the $ is not optional!
    @test all(B .== A[:, j])
    @einsum C[i] := A[i, 1]
    @test all(C .== A[:, 1])

    D = zeros(10, 3)
    # @einsum D[i, 1] = A[i, :j]
    @einsum D[i, 1] = A[i, $j]
    @test isapprox(D[:, 1], A[:, j])
    # @einsum D[i, :j] = A[i, :j]
    @einsum D[i, $j] = A[i, $j]
    @test isapprox(D[:, j], A[:, j])
end

@testset "Better type inference on allocating arrays" begin
    B1 = ones(Int, 5)
    B2 = ones(Float32, 5)
    B3 = ones(5)
    C = randn(5)
    @einsum A1[i, j] := B1[i] * C[j]
    @einsum A2[i, j] := B2[i] * C[j]
    @einsum A3[i, j] := B3[i] * C[j]

    @test eltype(A1) == Float64
    @test eltype(A2) == Float64
    @test eltype(A3) == Float64
    @test isapprox(A1, A3)
    @test isapprox(A2, A3)
end

@testset "Scalar output, issue #37" begin

    a = [1 0; 0 1]
    @einsum b := a[i,i]
    @einsum c[] := a[i,i]

    @test b == c[] == 2
    @test b isa Int
    @test c isa Array{Int,0}

end

#========== some extra things from issues not tests ==========#

@testset "shifts, issue 6" begin
    # https://github.com/ahwillia/Einsum.jl/issues/6
    # discussion of things which there were attempts to make work

    A = zeros(10);
    X = randn(10);
    Y = randn(10);
    @einsum A[j] = X[j]*Y[j+3]
    @test isapprox(A, X.*[Y[4:end];zeros(3)])

    @einsum A2[j] := X[j+3]
    @test axes(A2,1) == -2:7

    offset = 3
    @einsum A3[j] := X[j+$offset]
    @test axes(A3,1) == -2:7

    @einsum A5[i] := X[i-5]
    @test size(A5) == (10,)
    @test_broken all(A5[6:end] .== X[1:5]) # not true in Tullio's conventions

    @einsum A6[i] := X[i+3]*X[i-3]
    # @test size(A) == (7,) # that can't make sense
    # @test isapprox(A[4:7], X[7:end].*X[1:4])

end
@testset "shifts, issue 12" begin
    # https://github.com/ahwillia/Einsum.jl/pull/12

    B = collect(1:10)

    # "produce errors?" -- no!
    A = zeros(5)
    @einsum A[i] = B[i+5]
    intersect(axes(A, 1), axes(B, 1) .- 5) # legal indices 1:5
    @test A[1] == 6

    A = zeros(5)
    @einsum A[i] = B[i-5] #
    intersect(axes(A, 1), axes(B, 1) .+ 5) # empty range, so nothing changed
    @test all(A .== 0)

    # "legal?" -- yes!
    A = zeros(10)
    @einsum A[i] = B[i+5]
    @test A == [6, 7, 8, 9, 10, 0, 0, 0, 0, 0]

    A = zeros(10)
    @einsum A[i] = B[i-5]
    @test A == [0, 0, 0, 0, 0, 1, 2, 3, 4, 5]

#=
    # "Side note. This is silly, but would be kind of cool to support this:"

    B = collect(1:10)
    A = zeros(10)
    @einsum A[i] = B[i>>5]
    A == [6, 7, 8, 9, 10, 1, 2, 3, 4, 5]

    @einsum A[i] = B[i<<2]
    A == [3, 4, 5, 6, 7, 8, 9, 10, 1, 2]

    # Cyclic indices would be neat but this isn't the notation
=#
end

@tullio avx=true


text:
function LazyLoadedGravityFields(arg0::DataProvidersManager, arg1::TimeScale)
    return LazyLoadedGravityFields((DataProvidersManager, TimeScale), arg0, arg1)
end

function add_default_ocean_tides_readers(obj::LazyLoadedGravityFields)
    return jcall(obj, "addDefaultOceanTidesReaders", void, ())
end

function add_default_potential_coefficients_readers(obj::LazyLoadedGravityFields)
    return jcall(obj, "addDefaultPotentialCoefficientsReaders", void, ())
end

function add_ocean_tides_reader(obj::LazyLoadedGravityFields, arg0::OceanTidesReader)
    return jcall(obj, "addOceanTidesReader", void, (OceanTidesReader,), arg0)
end

function add_potential_coefficients_reader(obj::LazyLoadedGravityFields, arg0::PotentialCoefficientsReader)
    return jcall(obj, "addPotentialCoefficientsReader", void, (PotentialCoefficientsReader,), arg0)
end

function clear_ocean_tides_readers(obj::LazyLoadedGravityFields)
    return jcall(obj, "clearOceanTidesReaders", void, ())
end

function clear_potential_coefficients_readers(obj::LazyLoadedGravityFields)
    return jcall(obj, "clearPotentialCoefficientsReaders", void, ())
end

function configure_ocean_load_deformation_coefficients(obj::LazyLoadedGravityFields, arg0::OceanLoadDeformationCoefficients)
    return jcall(obj, "configureOceanLoadDeformationCoefficients", void, (OceanLoadDeformationCoefficients,), arg0)
end

function get_constant_normalized_provider(obj::LazyLoadedGravityFields, arg0::jint, arg1::jint)
    return jcall(obj, "getConstantNormalizedProvider", NormalizedSphericalHarmonicsProvider, (jint, jint), arg0, arg1)
end

function get_constant_unnormalized_provider(obj::LazyLoadedGravityFields, arg0::jint, arg1::jint)
    return jcall(obj, "getConstantUnnormalizedProvider", UnnormalizedSphericalHarmonicsProvider, (jint, jint), arg0, arg1)
end

function get_normalized_provider(obj::LazyLoadedGravityFields, arg0::jint, arg1::jint)
    return jcall(obj, "getNormalizedProvider", NormalizedSphericalHarmonicsProvider, (jint, jint), arg0, arg1)
end

function get_ocean_load_deformation_coefficients(obj::LazyLoadedGravityFields)
    return jcall(obj, "getOceanLoadDeformationCoefficients", OceanLoadDeformationCoefficients, ())
end

function get_ocean_tides_waves(obj::LazyLoadedGravityFields, arg0::jint, arg1::jint)
    return jcall(obj, "getOceanTidesWaves", List, (jint, jint), arg0, arg1)
end

function get_unnormalized_provider(obj::LazyLoadedGravityFields, arg0::jint, arg1::jint)
    return jcall(obj, "getUnnormalizedProvider", UnnormalizedSphericalHarmonicsProvider, (jint, jint), arg0, arg1)
end

function read_gravity_field(obj::LazyLoadedGravityFields, arg0::jint, arg1::jint)
    return jcall(obj, "readGravityField", PotentialCoefficientsReader, (jint, jint), arg0, arg1)
end


text:
"""
    MultiPathfinderResult

Container for results of multi-path Pathfinder.

# Fields
- `input`: User-provided input object, e.g. either `logp`, `(logp, ∇logp)`, `optim_fun`,
    `optim_prob`, or another object.
- `optimizer`: Optimizer used for maximizing the log-density
- `rng`: Pseudorandom number generator that was used for sampling
- `optim_prob::SciMLBase.OptimizationProblem`: Otimization problem used for
    optimization
- `logp`: Log-density function
- `fit_distribution::Distributions.MixtureModel`: uniformly-weighted mixture of ELBO-
    maximizing multivariate normal distributions from each run.
- `draws::AbstractMatrix{<:Real}`: draws from `fit_distribution` with size `(dim, ndraws)`,
    potentially resampled using importance resampling to be closer to the target
    distribution.
- `draw_component_ids::Vector{Int}`: component id of each draw in `draws`.
- `fit_distribution_transformed`: `fit_distribution` transformed to the same space as the
    user-supplied target distribution. This is only different from `fit_distribution` when
    integrating with other packages, and its type depends on the type of `input`.
- `draws_transformed`: `draws` transformed to be draws from `fit_distribution_transformed`.
- `pathfinder_results::Vector{<:PathfinderResult}`: results of each single-path Pathfinder
    run.
- `psis_result::Union{Nothing,<:PSIS.PSISResult}`: If importance resampling was used, the
    result of Pareto-smoothed importance resampling. `psis_result.pareto_shape` also
    diagnoses whether `draws` can be used to compute estimates from the target distribution.
    See [`PSIS.PSISResult`](https://psis.julia.arviz.org/stable/api/#PSIS.PSISResult) for
    details
"""
struct MultiPathfinderResult{I,O,R,OF,LP,FD,D,C,FDT,DT,PFR,PR}
    input::I
    optimizer::O
    rng::R
    optim_fun::OF
    logp::LP
    fit_distribution::FD
    draws::D
    draw_component_ids::C
    fit_distribution_transformed::FDT
    draws_transformed::DT
    pathfinder_results::PFR
    psis_result::PR
end

function Base.show(io::IO, ::MIME"text/plain", result::MultiPathfinderResult)
    println(io, "Multi-path Pathfinder result")
    println(io, "  runs: $(length(result.pathfinder_results))")
    print(io, "  draws: $(size(result.draws, 2))")
    psis_result = result.psis_result
    if psis_result !== nothing
        println(io)
        k = psis_result.pareto_shape
        assessment = if k > 1
            "very bad"
        elseif k > 0.7
            "bad"
        elseif k > 0.5
            "ok"
        else
            "good"
        end
        print(io, "  Pareto shape diagnostic: ", round(k; digits=2), " ($assessment)")
    end
    return nothing
end

"""
    multipathfinder(logp, ndraws; kwargs...)
    multipathfinder(logp, ∇logp, ndraws; kwargs...)
    multipathfinder(fun::SciMLBase.OptimizationFunction, ndraws; kwargs...)

Run [`pathfinder`](@ref) multiple times to fit a multivariate normal mixture model.

For `nruns=length(init)`, `nruns` parallel runs of pathfinder produce `nruns` multivariate
normal approximations ``q_k = q(\\phi | \\mu_k, \\Sigma_k)`` of the posterior. These are
combined to a mixture model ``q`` with uniform weights.

``q`` is augmented with the component index to generate random samples, that is, elements
``(k, \\phi)`` are drawn from the augmented mixture model
```math
\\tilde{q}(\\phi, k | \\mu, \\Sigma) = K^{-1} q(\\phi | \\mu_k, \\Sigma_k),
```
where ``k`` is a component index, and ``K=`` `nruns`. These draws are then resampled with
replacement. Discarding ``k`` from the samples would reproduce draws from ``q``.

If `importance=true`, then Pareto smoothed importance resampling is used, so that the
resulting draws better approximate draws from the target distribution ``p`` instead of
``q``. This also prints a warning message if the importance weighted draws are unsuitable
for approximating expectations with respect to ``p``.

# Arguments
- `logp`: a callable that computes the log-density of the target distribution.
- `∇logp`: a callable that computes the gradient of `logp`. If not provided, `logp` is
    automatically differentiated using the backend specified in `ad_backend`.
- `fun::SciMLBase.OptimizationFunction`: an optimization function that represents
    `-logp(x)` with its gradient. It must have the necessary features (e.g. a Hessian
    function) for the chosen optimization algorithm. For details, see
    [GalacticOptim.jl: OptimizationFunction](https://galacticoptim.sciml.ai/stable/API/optimization_function/).
- `ndraws::Int`: number of approximate draws to return

# Keywords
- `init`: iterator of length `nruns` of initial points of length `dim` from which each
    single-path Pathfinder run will begin. `length(init)` must be implemented. If `init` is
    not provided, `dim` and `nruns` must be.
- `nruns::Int`: number of runs of Pathfinder to perform. Ignored if `init` is provided.
- `ad_backend=AD.ForwardDiffBackend()`: AbstractDifferentiation.jl AD backend used to
    differentiate `logp` if `∇logp` is not provided.
- `ndraws_per_run::Int`: The number of draws to take for each component before resampling.
    Defaults to a number such that `ndraws_per_run * nruns > ndraws`.
- `importance::Bool=true`: Perform Pareto smoothed importance resampling of draws.
- `rng::AbstractRNG=Random.GLOBAL_RNG`: Pseudorandom number generator. It is recommended to
    use a parallelization-friendly PRNG like the default PRNG on Julia 1.7 and up.
- `executor::Transducers.Executor`: Transducers.jl executor that determines if and how
    to run the single-path runs in parallel. If `rng` is known to be thread-safe, the
    default is `Transducers.PreferParallel()` (parallel executation, defaulting to
    multi-threading). Otherwise, it is `Transducers.SequentialEx()` (no parallelization).
- `executor_per_run::Transducers.Executor=Transducers.SequentialEx()`: Transducers.jl
    executor used within each run to parallelize PRNG calls. Defaults to no parallelization.
    See [`pathfinder`](@ref) for a description.
- `kwargs...` : Remaining keywords are forwarded to [`pathfinder`](@ref).

# Returns
- [`MultiPathfinderResult`](@ref)
"""
function multipathfinder end

function multipathfinder(
    logp, ndraws::Int; ad_backend=AD.ForwardDiffBackend(), input=logp, kwargs...
)
    return multipathfinder(build_optim_function(logp; ad_backend), ndraws; input, kwargs...)
end
function multipathfinder(
    logp,
    ∇logp,
    ndraws::Int;
    ad_backend=AD.ForwardDiffBackend(),
    input=(logp, ∇logp),
    kwargs...,
)
    return multipathfinder(
        build_optim_function(logp, ∇logp; ad_backend), ndraws; input, kwargs...
    )
end
function multipathfinder(
    optim_fun::SciMLBase.OptimizationFunction,
    ndraws::Int;
    init=nothing,
    input=optim_fun,
    nruns::Int=init === nothing ? -1 : length(init),
    ndraws_elbo::Int=DEFAULT_NDRAWS_ELBO,
    ndraws_per_run::Int=max(ndraws_elbo, cld(ndraws, max(nruns, 1))),
    rng::Random.AbstractRNG=Random.GLOBAL_RNG,
    history_length::Int=DEFAULT_HISTORY_LENGTH,
    optimizer=default_optimizer(history_length),
    executor::Transducers.Executor=_default_executor(rng),
    executor_per_run=Transducers.SequentialEx(),
    importance::Bool=true,
    kwargs...,
)
    if optim_fun.grad === nothing || optim_fun.grad isa Bool
        throw(ArgumentError("optimization function must define a gradient function."))
    end
    if init === nothing
        nruns > 0 || throw(
            ArgumentError("A positive `nruns` must be set or `init` must be provided.")
        )
        _init = fill(init, nruns)
    else
        _init = init
    end
    if ndraws > ndraws_per_run * nruns
        @warn "More draws requested than total number of draws across replicas. Draws will not be unique."
    end
    logp(x) = -optim_fun.f(x, nothing)

    # run pathfinder independently from each starting point
    trans = Transducers.Map() do (init_i)
        return pathfinder(
            optim_fun;
            rng,
            history_length,
            optimizer,
            ndraws=ndraws_per_run,
            init=init_i,
            executor=executor_per_run,
            ndraws_elbo,
            kwargs...,
        )
    end
    iter_sp = Transducers.withprogress(_init; interval=1e-3) |> trans
    pathfinder_results = Folds.collect(iter_sp, executor)
    fit_distributions =
        pathfinder_results |> Transducers.Map(x -> x.fit_distribution) |> collect
    draws_all = reduce(hcat, pathfinder_results |> Transducers.Map(x -> x.draws))

    # draw samples from augmented mixture model
    inds = axes(draws_all, 2)
    sample_inds, psis_result = if importance
        log_densities_fit =
            pathfinder_results |>
            Transducers.MapCat() do x
                return Distributions.logpdf(x.fit_distribution, x.draws)
            end |>
            collect
        iter_logp = eachcol(draws_all) |> Transducers.Map(logp)
        log_densities_target = Folds.collect(iter_logp, executor)
        log_densities_ratios = log_densities_target - log_densities_fit
        resample(rng, inds, log_densities_ratios, ndraws)
    else
        resample(rng, inds, ndraws), nothing
    end

    fit_distribution = Distributions.MixtureModel(fit_distributions)
    draws = draws_all[:, sample_inds]

    # get component ids (k) of draws in ϕ
    draw_component_ids = cld.(sample_inds, ndraws_per_run)

    # placeholders
    fit_distribution_transformed = fit_distribution
    draws_transformed = draws

    return MultiPathfinderResult(
        input,
        optimizer,
        rng,
        optim_fun,
        logp,
        fit_distribution,
        draws,
        draw_component_ids,
        fit_distribution_transformed,
        draws_transformed,
        pathfinder_results,
        psis_result,
    )
end


text:
export Dirac

struct Dirac{X} <: AbstractMeasure
    supp :: X
end

Dirac(x::X) = Dirac{X}(x)

sampletype(μ::Dirac{X}) = X

function (μ::Dirac{X})(s) where {X}
    μ.supp ∈ s && return 1
    return 0
end

function logdensity(μ::Dirac{X}, x::X)
    μ.supp ∈ s && return 0.0
    return -Inf
end

Base.rand(μ::Dirac) = μ.supp


text:
function Quaternion(arg0::Vector{jdouble})
    return Quaternion((Vector{jdouble},), arg0)
end

function Quaternion(arg0::jdouble, arg1::Vector{jdouble})
    return Quaternion((jdouble, Vector{jdouble}), arg0, arg1)
end

function Quaternion(arg0::jdouble, arg1::jdouble, arg2::jdouble, arg3::jdouble)
    return Quaternion((jdouble, jdouble, jdouble, jdouble), arg0, arg1, arg2, arg3)
end

function add(::Type{Quaternion}, arg0::Quaternion, arg1::Quaternion)
    return jcall(Quaternion, "add", Quaternion, (Quaternion, Quaternion), arg0, arg1)
end

function add(obj::Quaternion, arg0::Quaternion)
    return jcall(obj, "add", Quaternion, (Quaternion,), arg0)
end

function dot_product(::Type{Quaternion}, arg0::Quaternion, arg1::Quaternion)
    return jcall(Quaternion, "dotProduct", jdouble, (Quaternion, Quaternion), arg0, arg1)
end

function dot_product(obj::Quaternion, arg0::Quaternion)
    return jcall(obj, "dotProduct", jdouble, (Quaternion,), arg0)
end

function equals(obj::Quaternion, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function equals(obj::Quaternion, arg0::Quaternion, arg1::jdouble)
    return jcall(obj, "equals", jboolean, (Quaternion, jdouble), arg0, arg1)
end

function get_conjugate(obj::Quaternion)
    return jcall(obj, "getConjugate", Quaternion, ())
end

function get_inverse(obj::Quaternion)
    return jcall(obj, "getInverse", Quaternion, ())
end

function get_norm(obj::Quaternion)
    return jcall(obj, "getNorm", jdouble, ())
end

function get_positive_polar_form(obj::Quaternion)
    return jcall(obj, "getPositivePolarForm", Quaternion, ())
end

function get_q0(obj::Quaternion)
    return jcall(obj, "getQ0", jdouble, ())
end

function get_q1(obj::Quaternion)
    return jcall(obj, "getQ1", jdouble, ())
end

function get_q2(obj::Quaternion)
    return jcall(obj, "getQ2", jdouble, ())
end

function get_q3(obj::Quaternion)
    return jcall(obj, "getQ3", jdouble, ())
end

function get_scalar_part(obj::Quaternion)
    return jcall(obj, "getScalarPart", jdouble, ())
end

function get_vector_part(obj::Quaternion)
    return jcall(obj, "getVectorPart", Vector{jdouble}, ())
end

function hash_code(obj::Quaternion)
    return jcall(obj, "hashCode", jint, ())
end

function is_pure_quaternion(obj::Quaternion, arg0::jdouble)
    return jcall(obj, "isPureQuaternion", jboolean, (jdouble,), arg0)
end

function is_unit_quaternion(obj::Quaternion, arg0::jdouble)
    return jcall(obj, "isUnitQuaternion", jboolean, (jdouble,), arg0)
end

function multiply(::Type{Quaternion}, arg0::Quaternion, arg1::Quaternion)
    return jcall(Quaternion, "multiply", Quaternion, (Quaternion, Quaternion), arg0, arg1)
end

function multiply(obj::Quaternion, arg0::Quaternion)
    return jcall(obj, "multiply", Quaternion, (Quaternion,), arg0)
end

function multiply(obj::Quaternion, arg0::jdouble)
    return jcall(obj, "multiply", Quaternion, (jdouble,), arg0)
end

function normalize(obj::Quaternion)
    return jcall(obj, "normalize", Quaternion, ())
end

function subtract(::Type{Quaternion}, arg0::Quaternion, arg1::Quaternion)
    return jcall(Quaternion, "subtract", Quaternion, (Quaternion, Quaternion), arg0, arg1)
end

function subtract(obj::Quaternion, arg0::Quaternion)
    return jcall(obj, "subtract", Quaternion, (Quaternion,), arg0)
end

function to_string(obj::Quaternion)
    return jcall(obj, "toString", JString, ())
end


text:
# Computes graph Stein discrepancy based on the affine objective and a greedy
# 2-spanner of the complete graph on sample points x(i,:) with edges weighted by
# ||x(i,:) - x(i2,:)||
#
# For each j in {1,...,d}, returns a solution g_j to
#   max_{g_j}
#   sum_i q(i) gobjcoeff(i,j) * g_j(x(i,:)) + <gradgobjcoeff(i,:,j), grad g_j(x(i,:)))>
#   subject to, for each spanner edge (i, i2),
#    |g_j(x(i,:))| <= c_{1}, ||grad g_j(x(i,:))||^* <= c_{2},
#    |g_j(x(i,:)) - g_j(x(i2,:))| <= c_{2} ||x(i,:) - x(i2,:)||,
#    ||grad g_j(x(i,:)) - grad g_j(x(i2,:))||^*
#      <= c_3 ||x(i,:) - x(i2,:)||
#    |g_j(x(i,:)) - g_j(x(i2,:)) - <grad g_j(x(i2,:)), x(i,:) - x(i2,:)>|
#      <= (c_3/2) (||x(i,:) - x(i2,:)||)^2
#    |g_j(x(i,:)) - g_j(x(i2,:)) - <grad g_j(x(i,:)), x(i,:) - x(i2,:)>|
#      <= (c_3/2) (||x(i,:) - x(i2,:)||)^2
# where the norm ||.|| is the L1 norm and ||.||^* is its dual.
#
# Args:
#  sample - SteinDiscrete object representing a sample
#  gobjcoefficients - a n x d matrix so jth column has coefficients of
#      the objective for gj
#  gradgobjcoefficients - a n x d x d tensor so [:,:,j] represents
#      the coefficients of grad gj on each row
#  solver - optimization program solver supported by JuMP. Examples:

function affine_graph_discrepancy(sample::SteinDiscrete,
                                  gobjcoefficients::Array{Float64,2},
                                  gradgobjcoefficients::Array{Float64,3},
                                  solver::AbstractMathProgSolver,
                                  operator::AbstractString;
                                  c1::Float64=1.0,
                                  c2::Float64=1.0,
                                  c3::Float64=1.0,
                                  supportlowerbounds=nothing,
                                  supportupperbounds=nothing)

    ## Extract inputs
    points = sample.support
    weights = sample.weights
    n = length(weights)
    d = size(points,2)

    ## reset the support bounds if needed
    if supportlowerbounds == nothing
        supportlowerbounds = repmat([-Inf], d)
    elseif isa(supportlowerbounds, Number)
        supportlowerbounds = rempat([supportlowerbounds], d)
    end
    if supportupperbounds == nothing
        supportupperbounds = repmat([Inf], d)
    elseif isa(supportupperbounds, Number)
        supportupperbounds = rempat([supportupperbounds], d)
    end

    ## Find spanner edge set
    println("[Computing spanner edges]")
    tic(); edges = getspanneredges(points); edgetime = toc()
    println("\tusing $(size(edges,1)) of $(binomial(n,2)) edges");
    if n > 1
        # Compute differences between points connected by an edge
        diffs = points[edges[:,1],:] - points[edges[:,2],:]
        # Compute L1 distances between points connected by an edge
        distances = sum(abs.(diffs),2)
        scaled_squared_distances = (c3/2)*distances.^2
    end

    ## Prepare return values
    objval = Array{Float64}(1, d)
    gopt = Array{Float64}(n, d)
    gradgopt = Array{Float64}(n, d, d)
    solvetime = Array{Float64}(1, d)

    # Distance cutoff for enforcing Lipschitz function constraints
    lipfunccutoff = 2*c1/c2
    # Distance cutoff for enforcing Lipschitz gradient constraints
    lipgradcutoff = 2*c2/c3
    # Distance cutoff for enforcing Taylor compatibility constraints
    taylorcutoff = 4*c2/c3

    ## Solve a different problem for each sample coordinate
    println("[Solving optimization program]")
    for j = 1:d
        tic()
        ## Define optimization problem
        m = Model(solver=solver)
        # Define optimization variables
        @variable(m, -c1 <= g[i=1:n] <= c1)
        @variable(m, -c2 <= gradg[i=1:n,k=1:d] <= c2)
        # Define the optimization objective
        gobj = AffExpr(g[1:n], vec(gobjcoefficients[:,j]), 0.0)
        gradgobj = AffExpr()
        for k = 1:d
            gradgobj += AffExpr(gradg[1:n,k], vec(gradgobjcoefficients[1:n,k,j]), 0.0)
        end
        @objective(m, Max, gobj + gradgobj)
        # Find finite limits of support in dimension j
        limits = [
            supportlowerbounds[j],
            supportupperbounds[j]
        ]
        limits = filter(isfinite, limits)
        # Add boundary constraints if needed
        for i = 1:n, bj = limits
            slackij = points[i,j] - bj
            # Add constraints to ensure gj can vanish on boundary
            # whenever abs(slackij) < lipfunccutoff
            # (otherwise constraints will never be active)
            if abs(slackij) < lipfunccutoff
                @constraint(m, g[i] >= -c2 * abs(slackij))
                @constraint(m, g[i] <= c2 * abs(slackij))
            end
            # Add constraints to ensure \grad gj can vanish in non-j dimension
            # whenever abs(slackij) < lipgradcutoff
            # (otherwise constraints will never be active)
            if abs(slackij) < lipgradcutoff
                constrained_dims = filter(x -> x != j, 1:d)
                for k = constrained_dims
                    @constraint(m, gradg[i,k] >= -c3 * abs(slackij))
                    @constraint(m, gradg[i,k] <= c3 * abs(slackij))
                end
            end
            # Add \grad gj constraints on jth dimension
            # whenever abs(slackij) < taylorcutoff
            # (otherwise constraints will never be active)
            if abs(slackij) < taylorcutoff
                @constraint(m, g[i] - slackij * gradg[i,j] <= (c3/2) * slackij^2)
                @constraint(m, g[i] - slackij * gradg[i,j] >= (-c3/2) * slackij^2)
            end
        end
        # Add pairwise constraints
        if n > 1
            for i = 1:length(distances)
                v1 = edges[i,1]; v2 = edges[i,2]
                # Add Lipschitz function constraints
                # whenever distance < lipgradcutoff
                # (otherwise constraints will never be active)
                if distances[i] < lipfunccutoff
                    @constraint(m, g[v1] - g[v2] <= c2 * distances[i])
                    @constraint(m, g[v1] - g[v2] >= -c2 * distances[i])
                end
                # Add Lipschitz gradient constraints
                # whenever distances[i] < lipgradcutoff
                # (otherwise constraints will never be active)
                if distances[i] < lipgradcutoff
                    for k = 1:d
                        @constraint(m, gradg[v1,k] - gradg[v2,k] <= c3*distances[i])
                        @constraint(m, gradg[v1,k] - gradg[v2,k] >= -c3*distances[i])
                    end
                end
                # Add Taylor compatibility constraints relating g and gradg
                # whenever distance < taylorcutoff
                # (otherwise constraints will never be active)
                if distances[i] < taylorcutoff
                    expr1 = AffExpr(vec(gradg[v1,:]), vec(diffs[i,:]), 0.0)
                    @constraint(m, g[v1] - g[v2] - expr1 <= scaled_squared_distances[i])
                    @constraint(m, g[v1] - g[v2] - expr1 >= -scaled_squared_distances[i])
                    expr2 = AffExpr(vec(gradg[v2,:]), vec(diffs[i,:]), 0.0)
                    @constraint(m, g[v1] - g[v2] - expr2 <= scaled_squared_distances[i])
                    @constraint(m, g[v1] - g[v2] - expr2 >= -scaled_squared_distances[i])
                end
            end
        end
        # Solve the problem
        @time status = JuMP.solve(m)
        solvetime[j] = toc()
        # Package the results
        objval[j] = getobjectivevalue(m)
        gopt[:,j] = getvalue(g)[1:n]
        gradgopt[:,:,j] = getvalue(gradg)[1:n, 1:d]
    end
    # Compute (T_Pg)(x) = {g part} + {grad g part}
    operatorgopt = sum(gopt .* gobjcoefficients,2)
    for j = 1:d
        operatorgopt += sum(gradgopt[:,:,j] .* gradgobjcoefficients[:,:,j],2)
    end

    AffineDiscrepancyResult(
        points,
        weights,
        edges,
        objval,
        gopt,
        gradgopt,
        operatorgopt,
        edgetime,
        solvetime,
        operator
    );
end


text:
# uniqueness checks
abstract type Uniqueness end

struct AllUniqueTrait <: Uniqueness end
const AllUnique = AllUniqueTrait()

struct NotUniqueTrait <: Uniqueness end
const NotUnique = NotUniqueTrait()

struct UnkownUniqueTrait <: Uniqueness end
const UnkownUnique = UnkownUniqueTrait()

Uniqueness(::T) where {T} = Uniqueness(T)
Uniqueness(::Type{T}) where {T} = UnkownUnique
Uniqueness(::Type{T}) where {T<:AbstractRange} = AllUnique

function check_index_uniqueness(idx::AbstractIndex, u::Uniqueness=NotUnique)
    return _check_index_uniqueness(keys(idx), u)
end
function check_index_uniqueness(idx, u::Uniqueness=NotUnique)
    return _check_index_uniqueness(idx, u)
end

_check_index_uniqueness(idx, ::AllUniqueTrait) = nothing
function _check_index_uniqueness(idx, ::UnkownUniqueTrait)
    return allunique(idx) ? nothing : _check_index_uniqueness(idx, NotUnique)
end
function _check_index_uniqueness(idx, ::NotUniqueTrait)
    error("all keys of each index must be unique.")
end

# length checks
abstract type AbstractLengthCheck end

struct LengthCheckedTrait <: AbstractLengthCheck end
const LengthChecked = LengthCheckedTrait()

struct LengthNotCheckedTrait <: AbstractLengthCheck end
const LengthNotChecked = LengthNotCheckedTrait()

function check_index_length(axs, idx, lc::AbstractLengthCheck=LengthNotChecked)
    return _check_length(axs, idx, lc)
end

_check_length(axs, idx, ::LengthCheckedTrait) = nothing

function _check_length(axs, idx, ::LengthNotCheckedTrait)
    if length(axs) == length(idx)
        return nothing
    else
        error("Length of parent axes and index must be of equal length, got an
               axis of length $(length(axs)) and index of length $(length(idx)).")
    end
end

# dimensions check
function check_dimensions(parent::AbstractArray{T,N1}, indices::Tuple{Vararg{Any,N2}}) where {T,N1,N2}
    N1 === N2 || error("Parent array has $N1 dimensions but there are $N2 indices.")
    return nothing
end

# all checks
function check_index_params(parent, indices, uc::Uniqueness, lc::AbstractLengthCheck)
    check_dimensions(parent, indices)        # ensure dimensions correspond to indices
    for (axs,idx) in zip(axes(parent),indices)
        check_index_length(axs, idx, lc)     # ensure length of index match axes
        check_index_uniqueness(idx, uc)      # ensure all keys are unique
    end
    return nothing
end


text:
function printgrid()
    dashed = " -"^4*" "
    plusdash = "+"*dashed*"+"*dashed*"+"
    verticals = ("|"*" "^9)^2*"|"
    s = (plusdash * "\n" * (verticals * "\n")^4)^2*plusdash * "\n"
    println(s)
end