text:
# functions from GasChromatographyTools.jl:
##---begin-UI-functions-----------------------------------------------------------------------------
## functions defining PlutoUI widgets for Pluto notebooks
"""
    UI_Column(sp)

Construct a combined PlutoUI widget for the settings of the GC Column with then selectable stationary phases `sp`.

# UI fields
* ``L``: column length in m.
* ``d``: column diameter in mm.
* ``d_f``: film thickness in µm.
* stat. phase: stationary phase of the column
* Gas: mobile phase
"""
function UI_Column(sp; default=(10.0, 0.25, 0.25, 1, "He"))
		PlutoUI.combine() do Child
			@htl("""

Column settings

			L [m]: $(
				Child(NumberField(0.1:0.1:100.0; default=default[1]))
			)  d [mm]: $(
				Child(NumberField(0.01:0.01:1.00; default=default[2]))
			)  d_f [µm]: $(
				Child(NumberField(0.01:0.01:1.00; default=default[3]))
			)  stat. phase: $(
				Child(Select(sp; default=sp[default[4]]))
			)  Gas: $(
				Child(Select(["He", "H2", "N2"]; default=default[5]))
			)

			""")
	end
end

"""
    UI_Program(; default=("0 60 300 300 120", "40 40 170 300 300", "0 0 40 60 0", "-3 -3 -3 -3 -3", "18 18 58 98 98", "0 0 0 0 0"))

Construct a combined PlutoUI widget for the settings of the program of a GC
column with or without a thermal gradient (depending on the default tuple).

# With thermal gradient
For default as a **tuple of six strings** the folwing fields will be shown in the
widget:
* example default tupel: default=("0 60 300 300 120", "40 40 170 300 300", "0 0 40 60 0", "-3 -3 -3 -3 -3", "18 18 58 98 98", "0 0 0 0 0")
*`time steps`: the time steps after which duration the values of temperature,
inlet pressure, ΔT and α are achieved by linear interpolation (in s).
*`temperature steps`: the temperature steps (in °C).
*`ΔT steps`: the steps of the temperature difference (in °C) between column inlet
and outlet.
*`α steps`: the steps of the gradient profile (α = 0 ... linear change of
temperature along column, α < 0 ... concave exponential profile, α > 0 ...
convexe exponential profile).
*``p_{in}`` steps: the steps of the inlet pressure (in kPa(g))
*``p_{out}`` steps: the steps of the outlet pressure (in kPa(a))

# Without thermal gradient
For a default as a **tuple of four strings** the folwing fields will be shown in the
widget:
* example default tupel: default=("0 60 600 120", "40 40 300 300", "18 18 98 98", "vacuum")
*`time steps`: the time steps after which duration the values of temperature,
inlet pressure, ΔT and α are achieved by linear interpolation (in s).
*`temperature steps`: the temperature steps (in °C).
``p_{in}`` steps: the steps of the inlet pressure (in kPa(g))
`column outlet` selection of the outlet of the colum, "vacuum" (``p_{out} =
0.0`` kPa(a)) or "atmosphere" (``p_{out} = 101.3`` kPa(a)).

"""
function UI_Program(; default=("0 60 300 300 120", "40 40 170 300 300", "0 0 40 60 0", "-3 -3 -3 -3 -3", "18 18 58 98 98", "0 0 0 0 0"))
	if length(default) == 6
		PlutoUI.combine() do Child
			@htl("""

Program settings

			Note: Same number of entrys for every text field.

			$(
				Child(TextField((50,1); default=default[1]))
			) time steps [s]

			$(
				Child(TextField((50,1); default=default[2]))
			) temperature steps [°C]

			$(
				Child(TextField((50,1); default=default[3]))
			) ΔT steps [°C]

			$(
				Child(TextField((50,1); default=default[4]))
			) α steps

			$(
				Child(TextField((50,1); default=default[5]))
			) p_in steps [kPa(g)]

			$(
				Child(TextField((50,1); default=default[6]))
			) p_out steps [kPa(a)]

			""")
		end
	elseif length(default) == 4
		PlutoUI.combine() do Child
			@htl("""

Program settings

			Note: Same number of entrys for every text field.

			$(
				Child(TextField((50,1); default=default[1]))
			) time steps [s]

			$(
				Child(TextField((50,1); default=default[2]))
			) temperature steps [°C]

			$(
				Child(TextField((50,1); default=default[3]))
			) p_in steps [kPa(g)]

			$(
				Child(Select(["vacuum", "atmosphere"]; default=default[4]))
				) column outlet

			""")
		end
	end
end

"""
    UI_Substance(sol; default=(1:4, ))

Construct a combined PlutoUI widget for the settings of the substances separated
in the simulated GC system with the selectable substances `subs`.

Depending on the tupel of `default` the widget is setup.

For `default = (1:4, 0.0, 0.0)` the UI fields are:
* Select Substances: Selection of the substances, which will be simulated,
  default selection = 1st to 4th substance.
* Injection time: Start time (in s) of the simulation. The same for all selected
  substances. Default is 0.0 s.
* Injection width: Peak width (in s) of all selected substances at the time of
  injection. Default is 0.0 s.

For `default = (1:4,)` the UI fields are:
* Select Substances: Selection of the substances, which will be simulated,
  default selection = 1st to 4th substance.
"""
function UI_Substance(sol; default=(1:4,))
	if length(sol)>10
		select_size = 10
	else
		select_size = length(sol)
	end
	if length(default) == 3
		PlutoUI.combine() do Child
			@htl("""

Substance settings

			Select Substances: $(
				Child(MultiSelect(sol; default=sol[default[1]], size=select_size))
			)
			Injection time [s]: $(
				Child(NumberField(0.0:0.1:100.0; default=default[2]))
			) and Injection width [s]: $(
				Child(NumberField(0.00:0.01:10.0; default=default[3]))
			)
			""")
		end
	elseif length(default) == 1
		PlutoUI.combine() do Child
			@htl("""

Substance settings

			Select Substances: $(
				Child(MultiSelect(sol; default=sol[default[1]], size=select_size))
			)
			""")
		end
	end
end

"""
    UI_Options()

Construct a combined PlutoUI widget for the settings of the options for the simulation.
"""
function UI_Options()
	PlutoUI.combine() do Child
		@htl("""

Option settings

		abstol: 1e $(
			Child(NumberField(-10:1:-3; default=-8))
		) reltol: 1e $(
			Child(NumberField(-8:1:-2; default=-5))
		) Tcontrol: $(
			Child(Select(["inlet", "outlet"]; default="inlet"))
		)
		""")
	end
end

"""
	setting_prog(prog_values)

Translates the Program parameters from a tuple defined by a PlutoUI widget into
the structure GasChromatographySimulator.Program.
"""
function setting_prog(prog_values, L)
	# make different variants based on the size/composition of `prog_values``
	time_steps = parse.(Float64, split(prog_values[1]))
	temp_steps = parse.(Float64, split(prog_values[2]))
	pin_steps = parse.(Float64, split(prog_values[3])).*1000.0.+101300.0
	if prog_values[4] == "vacuum"
		pout_steps = zeros(length(time_steps))
	elseif prog_values[4] == "atmosphere"
		pout_steps = 101300.0.*ones(length(time_steps))
	end
	prog = GasChromatographySimulator.Program( 	time_steps,
												temp_steps,
												pin_steps,
												pout_steps,
												L
												)
	return prog
end
##---end-UI-functions-------------------------------------------------------------------------------


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

function colwise(D::LocalGeoData, ::Type{AnisoDistance}, i::Int, xj::Vector{Int})
	Qi, ix = qmat(rotation(D,i),magnitude(D,i)), coords(D,i)
	d = Vector{Float64}()
	for j in xj
		Qj = qmat(rotation(D,j),magnitude(D,j))
		Qij = (Qi+Qj)/2
		jx = coords(D,j)
		push!(d,Distances.evaluate(Mahalanobis(Symmetric(Qij)), ix, jx))
	end
	d
end

function colwise(D::LocalGeoData, ::Type{KernelVariogram}, i::Int, xj::Vector{Int})
	Qi, ix = qmat(rotation(D,i),magnitude(D,i)), centro(D,i)
	d = Vector{Float64}()
	for j in xj
		Qj, jx = qmat(rotation(D,j),magnitude(D,j)), centro(D,j)
		dγ = abs(sill(D.refvario) - kccov(D.refvario, ix, jx, Qi, Qj))
		push!(d,dγ)
	end
	d
end

function colwise(D::LocalGeoData, ::Type{GraphDistance}, i::Int, j::Ints)
	d = dijkstra_shortest_paths(D.graph,i).dists
	view(d,j)
end

function evaluate(D::LocalGeoData, ::Type{AnisoDistance}, i::Int, j::Int)
	Qi, Qj = [qmat(rotation(D,x),magnitude(D,x)) for x in (i,j)]
	ix, jx = [coords(D,x) for x in (i,j)]
	Qij = (Qi+Qj)/2
	Distances.evaluate(Mahalanobis(Symmetric(Qij)), ix, jx)
end

function evaluate(D::LocalGeoData, ::Type{KernelVariogram}, i::Int, j::Int)
	Qi, Qj = [qmat(rotation(D,x),magnitude(D,x)) for x in (i,j)]
	ix, jx = [centro(D,x) for x in (i,j)]
	abs(sill(D.refvario) - kccov(D.refvario, ix, jx, Qi, Qj))
end

evaluate(D::LocalGeoData, ::Type{GraphDistance}, i::Int, j::Int) = colwise(GraphDistance, D, i, j)[1]

ext: jl


text:
using KernelFunctions,Plots,JLD,Distributions

include("../algo/kernels.jl")
include("../algo/regularized_Picard.jl")
include("../algo/utils.jl")

"""
    estimate_Gaussian(s,n,sigma,lambda,tol,intensity,it_max,p)

Estimates likelihood and correlation kernels and plots several figures. No output.

# Arguments
- `s::Integer`: the number of Gaussian-DPP samples used for estimation (up to 10).
- `n:Integer`: the number of uniformly sampled point for approximating the Fredholm determinant
- `sigma:Float`: RKHS kernel bandwidth
- `lambda:Float`: regularization parameter
- `tol:Float`: relative objective variation desired for regularized Picard iteration
- `intensity:Integer=100`: choose Gaussian-DPP data-samples with intensity 100 (or 50)
- `it_max:Integer=100000`: largest number of iterations for Picard
- `p:Integer= 10000`: the number of uniformly sampled point for approximating correlation kernel

"""
function estimate_Gaussian(s::Int64,n::Int64,sigma::Float64,lambda::Float64,tol::Float64,intensity::Int64=100,it_max::Int64=100000,p::Int64 = 10000)

    # constant added on diagonal of matrices
    epsilon = 1e-10;

    # sample uniformly in [0,1]^2 n points
    Fredholm_sple = rand(Uniform(0,1), n,2);

    # load DPP samples
    total_sples, indices_Fredholm_sple, indices_DPP_sples = add_DPP_samples_to_Fredholm_samples(s,Fredholm_sple,intensity)

    # create full kernel matrix
    x = (total_sples)'/sigma;
    k = SqExponentialKernel();
    K = kernelmatrix(k, x) + epsilon *I ;

    #####################################################################################################
    ## estimate B: insample likelihood
    #####################################################################################################

    print("Solving discrete problem with regularized Picard ...\n")

    X = randn(size(K));
    B = X*X'+ UniformScaling(epsilon);# initial positive definite B

    use_inverse = false
    @time begin
        B, R, obj, i_stop = regularized_Picard(B, K, indices_DPP_sples, indices_Fredholm_sple, lambda, it_max ,tol,use_inverse)
    end

    # create repo for storing results
    foldername = create_repo(s,n,p,sigma,lambda,tol);

    print("Plotting ...\n")

    #####################################################################################################
    ## plot objectives
    #####################################################################################################
    i_start = 5; # skip first iterates
    plot_objectives(i_start,i_stop, obj);
    savefig(foldername*"/objectives.pdf")

    #####################################################################################################
    ## scatter plot of in-sample likelihood kernel
    #####################################################################################################

    # color is diagonal of likelihood matrix
    plot_insample_likelihood(B, R, indices_Fredholm_sple, total_sples);
    savefig(foldername*"/in-sample_likelihood.pdf")

    #####################################################################################################
    ## estimate out-of-sample likelihood kernel on grid
    #####################################################################################################

    # test samples on a grid n_side x n_side in [a,b]^2
    n_side = 30;
    a = 0.; b = 1.;
    n_test = n_side^2;

    test_sples = flat_square_2d_grid(n_test, a, b);

    # evaluate likelihood kernel on the grid
    GramA = likelihood_kernel_Gram(B,R,total_sples,test_sples,k,sigma);

    display(heatmap(GramA,title="Gram matrix of estimated likelihood kernel"));
    savefig(foldername*"/Gram_likelihood.pdf")

    #####################################################################################################
    ## estimate out-of-sample correlation kernel on the same grid
    #####################################################################################################

    ########## Approximate correlation kernel ##########

    # sampling uniformly p points in window [0,1]^2
    c_1 = 0.; c_2 = 1.;
    unif_sples_correlation = rand(Uniform(c_1,c_2), p,2);

    # use these samples to compute correlation kernel
    GramK = correlation_kernel_Gram(B,R,unif_sples_correlation,total_sples,test_sples,k,sigma);

    display(heatmap(GramK,title="Gram matrix of estimated correlation kernel"));
    savefig(foldername*"/Gram_correlation.pdf")

    heatmap_intensity_correlation_kernel(GramK,n_side)
    savefig(foldername*"/intensity_correlation.pdf")

    # select one point x_0 almost at center of grid and display heatmap of k(x,x_0)
    heatmap_slice_Gram_matrix_at_grid_center(GramK,n_side);
    savefig(foldername*"/slice_correlation.pdf")

    ########## Exact correlation kernel ##########

    alpha_0 = 0.05;
    rho_0 = 100;
    x_test_0 = (test_sples)'/(alpha_0/sqrt(2));
    # NB: srt(2) to match the kernel definition in genGaussianSpatstats.R

    k_0 = SqExponentialKernel();
    GramK_0 = rho_0 * kernelmatrix(k_0, x_test_0) + epsilon *I ; # positive definite

    # display heatmap of GramK_0
    display(heatmap(GramK_0,title="Gram matrix of exact correlation kernel"));
    savefig(foldername*"/Gram_correlation_exact.pdf")

    #####################################################################################################
    ## save outputs
    #####################################################################################################

    save_results_to_jld(B, R, n, total_sples, GramK, GramA, GramK_0, obj, i_stop, foldername)    ;

end


text:
using JLD2, FileIO, Test

fn = joinpath(tempdir(), "test.jld2")

# Test load macros
jldopen(fn, "w") do f
    write(f, "loadmacrotestvar1", ['a', 'b', 'c'])
    write(f, "loadmacrotestvar2", 1)
end

@eval begin # wrapped in eval since @load with no args needs file at compile time
function func1()
    @load $fn
    @test loadmacrotestvar1 == ['a', 'b', 'c']
    @test loadmacrotestvar2 == 1
end
end

func1()

function func2()
    @load fn loadmacrotestvar1 loadmacrotestvar2
    @test loadmacrotestvar1 == ['a', 'b', 'c']
    @test loadmacrotestvar2 == 1
end

func2()

@test !isdefined(@__MODULE__, :loadmacrotestvar1) # should not be in global scope
@test (@eval @load $fn) == [:loadmacrotestvar1, :loadmacrotestvar2]
@test loadmacrotestvar1 == ['a', 'b', 'c']
@test loadmacrotestvar2 == 1

# Test save macros
hello = "world"
@save fn hello
jldopen(fn, "r") do f
    @test read(f, "hello") == "world"
end

cmd = """
using JLD2
hello = "there"
@save $(repr(fn))
"""
run(`$(Base.julia_cmd()) -e $cmd`)
jldopen(fn, "r") do f
    @test read(f, "hello") == "there"
end

# Test Dict save/load
save(format"JLD2", fn, Dict("the"=>"quick", "brown"=>"fox", "stuff"=>reshape(1:4, (2, 2))))
@test load(fn) == Dict("the"=>"quick", "brown"=>"fox", "stuff"=>reshape(1:4, (2, 2)))

# Test load/save with pairs
save(format"JLD2", fn, "jumps", "over", "the", "lazy", "dog", reshape(1:4, (2, 2)))
@test load(fn, "jumps", "the", "dog") == ("over", "lazy", reshape(1:4, (2, 2)))
@test load(fn, "jumps") == "over"

jldopen(fn, "r+") do f
    @test !isempty(f)
    @test haskey(f, "dog")
    @test !haskey(f, "notdog")
    @test keys(f) == ["jumps", "the", "dog"]
    x1 = rand(UInt8, 1024^2+1024)
    x2 = rand(UInt8, 1024)
    f["x1"] = x1
    f["x2"] = x2
    @test f["x1"] == x1
    @test f["x2"] == x2
end

# Issue #19
save(fn, Dict("a"=>[1,2,3]))
io = open(fn)
@info("The next error message (involving \"loading nothing\") is a sign of normal operation")
@test_throws MethodError load(io)
close(io)

# Issue #33
d = Dict("params/p1" => 1,
         "params/p2" => 2.,
         "params/p3/p1" => 94,
         "data" => [[1,2,3], [4.,5.,6]])
save(fn, d)
@test load(fn) == d

# Issue #106
mutable struct MyMutableTest
    a::Int
    b::Vector{Int}
end
Base.getproperty(df::MyMutableTest, s::Symbol) =
    throw(ArgumentError("should not be called"))
Base.setproperty!(df::MyMutableTest, s::Symbol, x::Int) =
    throw(ArgumentError("should not be called"))
Base.isequal(x::MyMutableTest, y::MyMutableTest) =
    isequal(getfield(x, :a), getfield(y, :a)) && isequal(getfield(x, :b), getfield(y, :b))
mmtd = Dict("A" => MyMutableTest(1, [10]))
save(fn, mmtd)
@test isequal(load(fn), mmtd)


text:
@testset "WignerSeitz" begin
    # --- `wignerseitz` ---
    # hexagonal lattice (example from space group 147)
    Rs = SVector{3,Float64}.([[1,0,0], [-.5,√3/2,0], [0,0,1.25]])        # direct basis
    Gs = SVector{3,Float64}.([[2π,2π/√3,0], [0,4π/√3,0], [0,0,2π/1.25]]) # reciprocal basis
    cell = wignerseitz(Gs)
    @test length(cell) == 8
    @test length(vertices(cell)) == 12
    @test basis(cell) == Gs

    # ------------------------------------------------------------------------------------ #
    # for CI and testing generally, we cannot depend on a specific sorting/ordering
    # from `wignerseitz` (and ultimately, `PySpatial.Voronoi` or `PySpatial.ConvexHull`),
    # so to be robust to this, we compare output after sorting to a canonical order
    function canonicalize!(c::Cell)
        vs = vertices(c)
        fs = faces(c)
        # sort the vertices in some canonical order (we also use `round` to be robust to
        # floating point differences across architectures)
        I = sortperm(vs, by=v->round.(v, digits=5))
        permute!(vs, I)
        # relabel the face-indices to match the new vertex sorting
        replace!.(fs, (I .=> eachindex(I))...)
        # sort the faces in some canonical order; note that each face
        # must themselves also be sorted (because they are cyclical)
        map!(circshift_minimal_first, fs, fs)
        sort!(fs)

        return c
    end

    function circshift_minimal_first(v::Vector)
        _, idx = findmin(v)
        return circshift(v, length(v)-idx+1)
    end

    # actually canonicalize `cell`
    canonicalize!(cell)
    # ------------------------------------------------------------------------------------ #

    # test show/display (canonicalization is crucial here)
    cell_show_reference = """
    Cell{3} (8 faces, 12 vertices):
      verts: [-0.666667, 0.333333, -0.5]
             [-0.666667, 0.333333, 0.5]
             [-0.333333, -0.333333, -0.5]
             [-0.333333, -0.333333, 0.5]
             [-0.333333, 0.666667, -0.5]
             [-0.333333, 0.666667, 0.5]
             [0.333333, -0.666667, -0.5]
             [0.333333, -0.666667, 0.5]
             [0.333333, 0.333333, -0.5]
             [0.333333, 0.333333, 0.5]
             [0.666667, -0.333333, -0.5]
             [0.666667, -0.333333, 0.5]
      faces: [1, 2, 6, 5]
             [1, 3, 4, 2]
             [1, 5, 9, 11, 7, 3]
             [2, 4, 8, 12, 10, 6]
             [3, 7, 8, 4]
             [5, 6, 10, 9]
             [7, 11, 12, 8]
             [9, 10, 12, 11]
      basis: [6.283185, 3.627599, 0.0]
             [0.0, 7.255197, 0.0]
             [0.0, 0.0, 5.026548]"""

    cell_show_referenceᶜ = """
    Cell{3} (8 faces, 12 vertices):
      verts: [-4.18879, 0.0, -2.513274]
             [-4.18879, 0.0, 2.513274]
             [-2.094395, -3.627599, -2.513274]
             [-2.094395, -3.627599, 2.513274]
             [-2.094395, 3.627599, -2.513274]
             [-2.094395, 3.627599, 2.513274]
             [2.094395, -3.627599, -2.513274]
             [2.094395, -3.627599, 2.513274]
             [2.094395, 3.627599, -2.513274]
             [2.094395, 3.627599, 2.513274]
             [4.18879, 0.0, -2.513274]
             [4.18879, 0.0, 2.513274]
      faces: [1, 2, 6, 5]
             [1, 3, 4, 2]
             [1, 5, 9, 11, 7, 3]
             [2, 4, 8, 12, 10, 6]
             [3, 7, 8, 4]
             [5, 6, 10, 9]
             [7, 11, 12, 8]
             [9, 10, 12, 11]
      basis: [6.283185, 3.627599, 0.0]
             [0.0, 7.255197, 0.0]
             [0.0, 0.0, 5.026548]""" # cartesian representation
    test_show(cell, cell_show_reference)
    test_show(Brillouin.cartesianize(cell), cell_show_referenceᶜ)

    # test iteration of Cell struct
    @test cell[1] ≈ [[-2/3,  1/3, -1/2],
                     [-2/3,  1/3, 1/2],
                     [-1/3,  2/3, 1/2],
                     [-1/3,  2/3, -1/2]]
    @test cell[4] ≈ [[-2/3,  1/3, 1/2],
                     [-1/3, -1/3, 1/2],
                     [ 1/3, -2/3, 1/2],
                     [ 2/3, -1/3, 1/2],
                     [ 1/3,  1/3, 1/2],
                     [-1/3,  2/3, 1/2]]

    # test that everything works the same if we use ordinary vectors instead of SVectors
    @test wignerseitz(Gs) ≈ wignerseitz(collect.(Gs))

    # test error on out-of-bounds indexing error
    @test_throws BoundsError cell[9]

    # test :triangles output
    cell′ = wignerseitz(Gs; merge = true)
    @test all(v -> length(v) == 3, vertices(cell′))
end


text:
"""
# module Tokenize

# Examples

```jldoctest
julia>
```
"""
module Tokenize

using Compat: isnothing
using IterTools: takewhile
using Parameters: @with_kw

export Tokenizer,
    update_chars,
    lex,
    lex_name,
    lex_string,
    lex_numeric

const PUNCTUATION = raw"=+-*/\\()[]{},:;%&~<>?`|$#@"
const WHITESPACE = " \t\r\v\f"  # '\v' => '\x0b', '\f' => '\x0c' in Python

@with_kw mutable struct Tokenizer
    characters = nothing
    prior_char::Union{Nothing, AbstractChar} = nothing
    char::Union{Nothing, AbstractChar} = nothing
    idx::Int = 0
    prior_delim::Union{Nothing, AbstractChar} = nothing
    group_token::Union{Nothing, AbstractChar} = nothing  # Set to true if inside a namelist group
end  # struct Tokenizer

"""
    update_chars(tk::Tokenizer)

Update the current charters in the tokenizer.
"""
function update_chars(tk::Tokenizer)
    tk.prior_char, tk.char = tk.char, next(tk.characters, '\n')
    tk.idx += 1
end  # function update_chars

function lex(tk::Tokenizer, line)
    tokens = String[]
    tk.idx = 0   # Bogus value to ensure idx = 1 after first iteration
    tk.characters = Iterators.Stateful(line)  # An iterator generated by `line`
    update_chars(tk)

    while tk.char != '\n'  # NOTE: Cannot be "\n", which is a string!
        # Update namelist group status
        occursin(tk.char, raw"&$") && (tk.group_token = tk.char)

        if !isnothing(tk.group_token) && ((tk.group_token, tk.char) in (('&', '/'), ('$', '$')))
            # A namelist ends, the value cannot be the default value (`nothing`)
            # Because it is being compared below
            tk.group_token = '\0'
        end

        word = ""  # Initialize or clear `word` if exists
        # Ignore whitespace
        if occursin(tk.char, WHITESPACE)  # " \t\r\x0b\x0c"
            while occursin(tk.char, WHITESPACE)
                word *= tk.char  # Read one char to `word`
                update_chars(tk)  # Read the next char until meet a non-whitespace char
            end
        # Ignore comment
        elseif occursin(tk.char, raw"!#") || isnothing(tk.group_token)  # Comment line
            # Abort the iteration and build the comment token
            word = line[tk.idx:end]  # There is no '\n' at line end, no worry! Lines are already separated at line ends
            tk.char = '\n'  # NOTE: Cannot be "\n", which is a string!
        # Lex string
        elseif occursin(tk.char, raw"\"'") || !isnothing(tk.prior_delim)  # Meet a string
            word = lex_string(tk)
        # Lex variable
        elseif isletter(tk.char)  # Meet a variable
            word = lex_name(tk, line)
        # Meet a sign
        elseif occursin(tk.char, "+-")
            # Lookahead to check for IEEE value
            # tk.characters, lookahead = repeat([Iterators.Stateful(collect(tk.characters))], 2)  # FIXME:
            # ieee_val = join(takewhile(isletter, lookahead), "")
            # if lowercase(ieee_val) in ("inf", "infinity", "nan")  # Meet an nan/infinity
            #     word = tk.char * ieee_val
            #     tk.characters = lookahead
            #     tk.prior_char = ieee_val[end]
            #     tk.char = next(lookahead, '\n')
            # else
            #     word = lex_numeric(tk)  # Meet a number
            # end
            word = lex_numeric(tk)  # Meet a number
        # Meet a number
        elseif isdigit(tk.char)
            word = lex_numeric(tk)
        # Meet a dot
        elseif tk.char == '.'
            update_chars(tk)
            if isdigit(tk.char)
                frac = lex_numeric(tk)  # A fraction of a number
                word = '.' * frac
            else
                word = '.'  # If not followed by a number
                while isletter(tk.char)
                    word *= tk.char
                    update_chars(tk)
                end
                # A word containing `.` ends
                if tk.char == '.'
                    word *= tk.char
                    update_chars(tk)
                end
            end
        # Meet a punctuation
        elseif occursin(tk.char, PUNCTUATION)
            word = tk.char
            update_chars(tk)
        else
            # This should never happen
            error("This should never happen!")
        end
        push!(tokens, string(word))
    end  # while loop
    return tokens
end  # function lex

"""
    lex_name(tk::Tokenizer, line)

Tokenize a Fortran name, such as a variable or subroutine.
"""
function lex_name(tk::Tokenizer, line::AbstractString)
    endindex = tk.idx
    for char in line[tk.idx:end]
        !isalnum(char) && !occursin(char, raw"'\"_") && break
        endindex += 1
    end

    word = line[tk.idx:endindex - 1]  # Do not include non-alphanumeric and non-'"_ characters in `word`

    tk.idx = endindex - 1  # Do not include non-alphanumeric and non-'"_ characters
    # Update iterator, minus first character which was already read
    # Continue iterating from `length(word)` => drop the first `length(word) - 1` characters
    tk.characters = Iterators.Stateful(collect(Iterators.drop(tk.characters, length(word) - 1)))
    update_chars(tk)
    return word
end  # function lex_name

"""
    lex_string(tk::Tokenizer)

Tokenize a Fortran string.
"""
function lex_string(tk::Tokenizer)
    word = ""

    if !isnothing(tk.prior_delim)  # A previous quotation mark presents
        delim = tk.prior_delim  # Read until `delim`
        tk.prior_delim = nothing
    else
        delim = tk.char  # No previous quotation mark presents
        word *= tk.char  # Read this character
        update_chars(tk)
    end

    while true
        if tk.char == delim
            # Check for escaped delimiters
            update_chars(tk)
            if tk.char == delim
                word *= repeat(delim, 2)
                update_chars(tk)
            else
                word *= delim
                break
            end
        elseif tk.char == '\n'
            tk.prior_delim = delim
            break
        else
            word *= tk.char
            update_chars(tk)
        end
    end

    return word
end  # function lex_string

"""
    lex_numeric(tk::Tokenizer)

Tokenize a Fortran numerical value.
"""
function lex_numeric(tk::Tokenizer)
    word = ""
    frac = false

    if tk.char == '-'
        word *= tk.char  # `word == "-"``
        update_chars(tk)
    end

    # Read as long as `tk.char` is a digit, or not a dot
    while isdigit(tk.char) || (tk.char == '.' && !frac)
        # Only allow one decimal point
        if tk.char == '.'
            frac = true  # If meet '.', break the loop
        end
        word *= tk.char
        update_chars(tk)
    end

    # Check for float exponent
    if occursin(tk.char, "eEdD")
        word *= tk.char
        update_chars(tk)
    end

    if occursin(tk.char, "+-")
        word *= tk.char
        update_chars(tk)
    end

    while isdigit(tk.char)
        word *= tk.char
        update_chars(tk)
    end

    return word
end  # function lex_numeric

isalnum(c) = isletter(c) || isnumeric(c)

function next(iterable, default)
    x = iterate(iterable)
    isnothing(x) ? default : first(x)
end

end


text:
# NOTE: This code is copied from: https://github.com/JuliaPOMDP/LocalApproximationValueIteration.jl
module LocalApproximationPolicyEvaluation
    using Random
    using Printf
    using POMDPs
    using POMDPModelTools
    using POMDPSimulators
    using GridInterpolations
    using LocalFunctionApproximation
    using Distributions
    import POMDPs: Solver, solve, Policy, action, value
    using POMDPLinter: @POMDP_require, @req, @subreq, @warn_requirements

    export LocalPolicyEvalSolver, LocalPolicyEvalPolicy, action_and_prob

    mutable struct LocalPolicyEvalSolver{I<:LocalFunctionApproximator, RNG<:AbstractRNG} <: Solver
        interp::I # Will be copied over by value to each policy
        action_probability::Function # a function with signature action_probability(mdp, s, a)
        max_iterations::Int64 # max number of iterations
        belres::Float64 # the Bellman Residual
        verbose::Bool # Whether to print while solving or not
        rng::RNG # Seed if req'd
        is_mdp_generative::Bool # Whether to treat underlying MDP model as generative
        n_generative_samples::Int64 # If underlying model generative, how many samples to use
    end

    # Default constructor
    function LocalPolicyEvalSolver(interp::I, action_probability::Function;
                                                    max_iterations::Int64=100, belres::Float64=1e-3,
                                                    verbose::Bool=false, rng::RNG=Random.GLOBAL_RNG,
                                                    is_mdp_generative::Bool=false, n_generative_samples::Int64=0) where {I<:LocalFunctionApproximator, RNG<:AbstractRNG}
        return LocalPolicyEvalSolver(interp, action_probability, max_iterations, belres, verbose, rng, is_mdp_generative, n_generative_samples)
    end

    # Unparameterized constructor just for getting requirements
    function LocalPolicyEvalSolver()
        throw(ArgumentError("LocalPolicyEvalSolver needs a LocalFunctionApproximator object for construction!"))
    end

    # NOTE : We work directly with the value function
    # And extract actions at the end by using the interpolation object
    mutable struct LocalPolicyEvalPolicy{I<:LocalFunctionApproximator, RNG<:AbstractRNG} <: Policy
        interp::I # General approximator to be used in VI
        action_probability::Function # a function with signature action_probability(mdp, s, a)
        action_map::Vector # Maps the action index to the concrete action type
        mdp::Union{MDP,POMDP} # Uses the model for indexing in the action function
        is_mdp_generative::Bool # (Copied from solver.is_mdp_generative)
        n_generative_samples::Int64 # (Copied from solver.n_generative_samples)
        rng::RNG # (Copied from solver.rng)
    end

    # The policy can be created using the MDP and solver information
    # The policy's function approximation object (interp) is obtained by deep-copying over the
    # solver's interp object. The other policy parameters are also obtained from the solver
    function LocalPolicyEvalPolicy(mdp::Union{MDP,POMDP},
                                                    solver::LocalPolicyEvalSolver)
        return LocalPolicyEvalPolicy(deepcopy(solver.interp),solver.action_probability, ordered_actions(mdp),mdp,
                                                      solver.is_mdp_generative,solver.n_generative_samples,solver.rng)
    end

    @POMDP_require solve(solver::LocalPolicyEvalSolver, mdp::Union{MDP,POMDP}) begin

        P = typeof(mdp)
        S = statetype(P)
        A = actiontype(P)
        @req discount(::P)
        @subreq ordered_actions(mdp)

        @req actionindex(::P, ::A)
        @req actions(::P, ::S)
        as = actions(mdp)
        a = first(as)

        @req convert_s(::Type{S}, ::AbstractVector{Float64}, ::P)
        @req convert_s(::Type{Vector{Float64}}, ::S, ::P)

        # Have different requirements depending on whether solver MDP is generative or explicit
        if solver.is_mdp_generative
            @req gen(::DDNOut{(:sp, :r)}, ::P, ::S, ::A, ::typeof(solver.rng))
        else
            @req transition(::P, ::S, ::A)
            pts = get_all_interpolating_points(solver.interp)
            pt = first(pts)
            ss = POMDPs.convert_s(S, pt, mdp)
            dist = transition(mdp, ss, a)
            D = typeof(dist)
            @req support(::D)
        end

    end

    function POMDPs.solve(solver::LocalPolicyEvalSolver, mdp::Union{MDP,POMDP})

        @warn_requirements solve(solver, mdp)

        # Ensure that generative model has a non-zero number of samples
        if solver.is_mdp_generative
            @assert solver.n_generative_samples > 0
        end

        # Solver parameters
        max_iterations = solver.max_iterations
        belres = solver.belres
        discount_factor = discount(mdp)

        # Initialize the policy
        policy = LocalPolicyEvalPolicy(mdp,solver)

        total_time = 0.0
        iter_time = 0.0

        # Get attributes of interpolator
        # Since the policy object is created by the solver, it directly
        # modifies the value of the interpolator of the created policy
        num_interps = n_interpolating_points(policy.interp)
        interp_points = get_all_interpolating_points(policy.interp)
        interp_values = get_all_interpolating_values(policy.interp)

        # Obtain the vector of states by converting the corresponding
        # vector of interpolation points/samples to the state type
        # using the user-provided convert_s function
        S = statetype(typeof(mdp))
        interp_states = Vector{S}(undef, num_interps)
        for (i, pt) in enumerate(interp_points)
            interp_states[i] = POMDPs.convert_s(S, pt, mdp)
        end

        # State transition dictionary
        next_state_dict = Dict()

        # Outer loop for Value Iteration
        for i = 1 : max_iterations
            residual::Float64 = 0.0
            iter_time = @elapsed begin

            for (istate,s) in enumerate(interp_states)
                sub_aspace = actions(mdp, s)

                if isterminal(mdp, s)
                    interp_values[istate] = 0.0
                else
                    old_util = interp_values[istate]
                    total_util = 0

                    for a in sub_aspace
                        u::Float64 = 0.0

                        # Do bellman backup based on generative / explicit model
                        if solver.is_mdp_generative
                            # Generative Model
                            for j in 1:solver.n_generative_samples
                                sp_point, r, isTerm = 0, 0, false
                                if solver.n_generative_samples ==1 && haskey(next_state_dict, (s,a))
                                    sp_point, r, isTerm = next_state_dict[(s,a)]
                                else
                                    sp, r = gen(mdp, s, a, solver.rng)
                                    isTerm = isterminal(mdp, sp)
                                    sp_point = POMDPs.convert_s(Vector{Float64}, sp, mdp)
                                    solver.n_generative_samples == 1  && (next_state_dict[(s,a)] = (sp_point, r, isTerm))
                                end

                                u += r

                                # Only interpolate sp if it is non-terminal
                                if !isTerm
                                    u += discount_factor*compute_value(policy.interp, sp_point)
                                end
                            end
                            u = u / solver.n_generative_samples
                        else
                            # Explicit Model
                            dist = transition(mdp, s, a)
                            for (sp, p) in weighted_iterator(dist)
                                p == 0.0 ? continue : nothing
                                r = reward(mdp, s, a, sp)
                                u += p*r

                            # Only interpolate sp if it is non-terminal
                                if !isterminal(mdp, sp)
                                    sp_point = POMDPs.convert_s(Vector{Float64}, sp, mdp)
                                    u += p * (discount_factor*LocalFunctionApproximation.compute_value(policy.interp, sp_point))
                                end
                            end # next-states
                        end
                        total_util += solver.action_probability(mdp, s, a)*u
                    end #action

                    # Update this interpolant value
                    interp_values[istate] = total_util
                    util_diff = abs(total_util - old_util)
                    util_diff > residual ? (residual = util_diff) : nothing
                end
            end #state

            end #time
            total_time += iter_time
            solver.verbose ? @printf("[Iteration %-4d] residual: %10.3G | iteration runtime: %10.3f ms, (%10.3G s total)\n", i, residual, iter_time*1000.0, total_time) : nothing
            residual < belres ? break : nothing

        end #main
        return policy
    end

    function POMDPs.value(policy::LocalPolicyEvalPolicy, s::S) where S
        # Call the conversion function on the state to get the corresponding vector
        # That represents the point at which to interpolate the function
        s_point = POMDPs.convert_s(Vector{Float64}, s, policy.mdp)
        val = LocalFunctionApproximation.compute_value(policy.interp, s_point)
        return val
    end

    # Not explicitly stored in policy - extract from value function interpolation
    POMDPs.action(policy::LocalPolicyEvalPolicy, s, rng= Random.GLOBAL_RNG) = action_and_prob(policy, s, rng)[1]

    function Distributions.logpdf(policy::LocalPolicyEvalPolicy, s, a)
        mdp = policy.mdp
        us = [value(policy, s, a)*policy.action_probability(mdp, s, a) for a in actions(policy.mdp, s)]
        if sum(us) == 0
            us = [policy.action_probability(mdp, s, a) for a in actions(policy.mdp, s)]
        end
        p = value(policy, s, a)*policy.action_probability(mdp, s, a) / sum(us)
        return (p == 0) ? log(policy.action_probability(mdp, s, a)) : log(p)
    end

    Distributions.logpdf(policy::LocalPolicyEvalPolicy, h::SimHistory) = sum([logpdf(policy, s, a) for (s,a) in eachstep(h, (:s, :a))])

    function action_and_prob(policy::LocalPolicyEvalPolicy, s, rng= Random.GLOBAL_RNG)
        mdp = policy.mdp
        us = [value(policy, s, a)*policy.action_probability(mdp, s, a) for a in actions(policy.mdp, s)]
        if sum(us) == 0
            us = [policy.action_probability(mdp, s, a) for a in actions(policy.mdp, s)]
        end
        us /= sum(us)
        ai = rand(rng, Categorical(us))
        policy.action_map[ai], us[ai]
    end

    # GlobalApproximationFailureProbe the action-value for some state-action pair
    # This is also used in the above function
    function POMDPs.value(policy::LocalPolicyEvalPolicy, s::S, a::A) where {S,A}

        mdp = policy.mdp
        discount_factor = discount(mdp)

        u::Float64 = 0.0

        # As in solve(), do different things based on whether
        # mdp is generative or explicit
        if policy.is_mdp_generative
            for j in 1:policy.n_generative_samples
                sp, r = gen(mdp, s, a, policy.rng)
                sp_point = POMDPs.convert_s(Vector{Float64}, sp, mdp)
                u += r
                if !isterminal(mdp, sp)
                    u += discount_factor*LocalFunctionApproximation.compute_value(policy.interp, sp_point)
                end
            end
            u = u / policy.n_generative_samples
        else
            dist = transition(mdp, s, a)
            for (sp, p) in weighted_iterator(dist)
                p == 0.0 ? continue : nothing
                r = reward(mdp, s, a, sp)
                u += p*r

                # Only interpolate sp if it is non-terminal
                if !isterminal(mdp, sp)
                    sp_point = POMDPs.convert_s(Vector{Float64}, sp, mdp)
                    u += p*(discount_factor*LocalFunctionApproximation.compute_value(policy.interp, sp_point))
                end
            end
        end

        return u
    end

end


text:
function layoutable(::Type{Button}, fig_or_scene::FigureLike; bbox = nothing, kwargs...)

    scene = get_scene(fig_or_scene)

    default_attrs = default_attributes(Button, scene).attributes
    theme_attrs = subtheme(scene, :Button)
    attrs = merge!(merge!(Attributes(kwargs), theme_attrs), default_attrs)

    @extract attrs (padding, textsize, label, font, halign, valign, cornerradius,
        cornersegments, strokewidth, strokecolor, buttoncolor,
        labelcolor, labelcolor_hover, labelcolor_active,
        buttoncolor_active, buttoncolor_hover, clicks)

    decorations = Dict{Symbol, Any}()

    layoutobservables = LayoutObservables{Button}(attrs.width, attrs.height, attrs.tellwidth, attrs.tellheight,
        halign, valign, attrs.alignmode; suggestedbbox = bbox)

    textpos = Node(Point2f0(0, 0))

    subarea = lift(layoutobservables.computedbbox) do bbox
        round_to_IRect2D(bbox)
    end
    subscene = Scene(scene, subarea, camera=campixel!)

    # buttonrect is without the left bottom offset of the bbox
    buttonrect = lift(layoutobservables.computedbbox) do bbox
        BBox(0, width(bbox), 0, height(bbox))
    end

    on(buttonrect) do rect
        textpos[] = Point2f0(left(rect) + 0.5f0 * width(rect), bottom(rect) + 0.5f0 * height(rect))
    end

    roundedrectpoints = lift(roundedrectvertices, buttonrect, cornerradius, cornersegments)

    mousestate = Node(:out)

    bcolors = (; out = buttoncolor, active = buttoncolor_active, hover = buttoncolor_hover)
    bcolor = lift((s,_...)->bcolors[s][], Any, mousestate, values(bcolors)...)
    button = poly!(subscene, roundedrectpoints, strokewidth = strokewidth, strokecolor = strokecolor,
        color = bcolor, raw = true, inspectable = false)
    decorations[:button] = button

    lcolors = (; out = labelcolor, active = labelcolor_active, hover = labelcolor_hover)
    lcolor = lift((s,_...)->lcolors[s][], Any, mousestate, values(lcolors)...)
    labeltext = text!(subscene, label, position = textpos, textsize = textsize, font = font,
        color = lcolor, align = (:center, :center), raw = true, space = :data, inspectable = false)

    decorations[:label] = labeltext

    # move text in front of background to be sure it's not occluded
    translate!(labeltext, 0, 0, 1)

    onany(label, textsize, font, padding) do label, textsize, font, padding
        textbb = FRect2D(boundingbox(labeltext))
        autowidth = width(textbb) + padding[1] + padding[2]
        autoheight = height(textbb) + padding[3] + padding[4]
        layoutobservables.autosize[] = (autowidth, autoheight)
    end

    mouseevents = addmouseevents!(scene, layoutobservables.computedbbox)

    onmouseover(mouseevents) do _
        mousestate[] = :hover
        return Consume(false)
    end

    onmouseout(mouseevents) do _
        mousestate[] = :out
        return Consume(false)
    end

    onmouseleftup(mouseevents) do _
        mousestate[] = :hover
        return Consume(true)
    end

    onmouseleftdown(mouseevents) do _
        mousestate[] = :active
        clicks[] = clicks[] + 1
        return Consume(true)
    end

    label[] = label[]
    # trigger bbox
    layoutobservables.suggestedbbox[] = layoutobservables.suggestedbbox[]

    Button(fig_or_scene, layoutobservables, attrs, decorations)
end


text:
# test/macros.jl
# Testing macros work correctly
#############################################################################
using JuMP, FactCheck
using Base.Test

# To ensure the tests work on Windows and Linux/OSX, we need
# to use the correct comparison operators
const leq = JuMP.repl[:leq]
const geq = JuMP.repl[:geq]
const  eq = JuMP.repl[:eq]
const Vert = JuMP.repl[:Vert]
const sub2 = JuMP.repl[:sub2]

facts("[macros] Check Julia expression parsing") do
    sumexpr = :(sum{x[i,j] * y[i,j], i = 1:N, j in 1:M; i != j})
    @fact sumexpr.head --> :curly
    @fact length(sumexpr.args) --> 5
    @fact sumexpr.args[1] --> :sum
    @fact sumexpr.args[2].head --> :parameters
    @fact sumexpr.args[3] --> :(x[i,j] * y[i,j])
    @fact sumexpr.args[4].head --> :(=)
    if VERSION < v"0.5.0-dev+3231"
        @fact sumexpr.args[5].head --> :in
    else
        @fact sumexpr.args[5].head --> :call
        @fact sumexpr.args[5].args[1] --> :in
    end

    sumexpr = :(sum{x[i,j] * y[i,j], i in 1:N, j = 1:M})
    @fact sumexpr.head --> :curly
    @fact length(sumexpr.args) --> 4
    @fact sumexpr.args[1] --> :sum
    @fact sumexpr.args[2] --> :(x[i,j] * y[i,j])
    if VERSION < v"0.5.0-dev+3231"
        @fact sumexpr.args[3].head --> :in
    else
        @fact sumexpr.args[3].head --> :call
        @fact sumexpr.args[3].args[1] --> :in
    end

    @fact sumexpr.args[4].head --> :(=)
end

# generator syntax only parses on 0.5
if VERSION >= v"0.5-dev+5475"
    eval("""
facts("[macros] Check Julia expression parsing (0.5)") do
    sumexpr = :(sum(x[i,j] * y[i,j] for i = 1:N, j in 1:M if i != j))
    @fact sumexpr.head --> :call
    @fact sumexpr.args[1] --> :sum
    @fact sumexpr.args[2].head --> :generator
    @fact sumexpr.args[2].args[1] --> :(x[i,j] * y[i,j])
    @fact sumexpr.args[2].args[2].head --> :filter
    @fact sumexpr.args[2].args[2].args[1] --> :(i != j)
    @fact sumexpr.args[2].args[2].args[2] --> :(i = 1:N)
    @fact sumexpr.args[2].args[2].args[3] --> :(j = 1:M)

    sumexpr = :(sum(x[i,j] * y[i,j] for i = 1:N, j in 1:M))
    @fact sumexpr.head --> :call
    @fact sumexpr.args[1] --> :sum
    @fact sumexpr.args[2].head --> :generator
    @fact sumexpr.args[2].args[1] --> :(x[i,j] * y[i,j])
    @fact sumexpr.args[2].args[2] --> :(i = 1:N)
    @fact sumexpr.args[2].args[3] --> :(j = 1:M)
end
"""); end

facts("[macros] Check @constraint basics") do
    m = Model()
    @variable(m, w)
    @variable(m, x)
    @variable(m, y)
    @variable(m, z)
    t = 10

    @constraint(m, 3x - y == 3.3(w + 2z) + 5)
    @fact string(m.linconstr[end]) --> "3 x - y - 3.3 w - 6.6 z $eq 5"
    @constraint(m, 3x - y == (w + 2z)*3.3 + 5)
    @fact string(m.linconstr[end]) --> "3 x - y - 3.3 w - 6.6 z $eq 5"
    @constraint(m, (x+y)/2 == 1)
    @fact string(m.linconstr[end]) --> "0.5 x + 0.5 y $eq 1"
    @constraint(m, -1 <= x-y <= t)
    @fact string(m.linconstr[end]) --> "-1 $leq x - y $leq 10"
    @constraint(m, -1 <= x+1 <= 1)
    @fact string(m.linconstr[end]) --> "-2 $leq x $leq 0"
    @constraint(m, -1 <= x <= 1)
    @fact string(m.linconstr[end]) --> "-1 $leq x $leq 1"
    @constraint(m, -1 <= x <= sum{0.5, i = 1:2})
    @fact string(m.linconstr[end]) --> "-1 $leq x $leq 1"
    @fact_throws @constraint(m, x <= t <= y)
    @fact macroexpand(:(@constraint(m, 1 >= x >= 0))).head --> :error
    @fact macroexpand(:(@constraint(1 <= x <= 2, foo=:bar))).head --> :error

    @expression(m, aff, 3x - y - 3.3(w + 2z) + 5)
    @fact string(aff) --> "3 x - y - 3.3 w - 6.6 z + 5"

    @constraint(m, 3 + 5*7 <= 0)
    @fact string(m.linconstr[end]) --> "0 $leq -38"

    @expression(m, qaff, (w+3)*(2x+1)+10)
    @fact string(qaff) --> "2 w*x + 6 x + w + 13"
end

facts("[macros] Checking @variable with reverse direction bounds") do
    m = Model()
    @variable(m, 3.2 >= x >= 1)
    @fact m.colLower --> [1.0]
    @fact m.colUpper --> [3.2]
end

facts("[macros] sum{}") do
    m = Model()
    @variable(m, x[1:3,1:3])
    @variable(m, y)
    C = [1 2 3; 4 5 6; 7 8 9]
    @constraint(m, sum{ C[i,j]*x[i,j], i in 1:2, j = 2:3 } <= 1)
    @fact string(m.linconstr[end]) --> "2 x[1,2] + 3 x[1,3] + 5 x[2,2] + 6 x[2,3] $leq 1"
    @constraint(m, sum{ C[i,j]*x[i,j], i = 1:3, j in 1:3; i != j} == y)
    @fact string(m.linconstr[end]) --> "2 x[1,2] + 3 x[1,3] + 4 x[2,1] + 6 x[2,3] + 7 x[3,1] + 8 x[3,2] - y $eq 0"

    @constraint(m, sum{ C[i,j]*x[i,j], i = 1:3, j = 1:i} == 0);
    @fact string(m.linconstr[end]) --> "x[1,1] + 4 x[2,1] + 5 x[2,2] + 7 x[3,1] + 8 x[3,2] + 9 x[3,3] $eq 0"

    @constraint(m, sum{ 0*x[i,1], i=1:3} == 0)
    @fact string(m.linconstr[end]) --> "0 $eq 0"

    @constraint(m, sum{ 0*x[i,1] + y, i=1:3} == 0)
    @fact string(m.linconstr[end]) --> "3 y $eq 0"

end

if VERSION >= v"0.5-dev+5475"
eval("""
facts("[macros] sum(generator)") do
    m = Model()
    @variable(m, x[1:3,1:3])
    @variable(m, y)
    C = [1 2 3; 4 5 6; 7 8 9]
    @constraint(m, sum( C[i,j]*x[i,j] for i in 1:2, j = 2:3 ) <= 1)
    @fact string(m.linconstr[end]) --> "2 x[1,2] + 3 x[1,3] + 5 x[2,2] + 6 x[2,3] $leq 1"
    @constraint(m, sum( C[i,j]*x[i,j] for i = 1:3, j in 1:3 if i != j) == y)
    @fact string(m.linconstr[end]) --> "2 x[1,2] + 3 x[1,3] + 4 x[2,1] + 6 x[2,3] + 7 x[3,1] + 8 x[3,2] - y $eq 0"

    @constraint(m, sum( C[i,j]*x[i,j] for i = 1:3, j = 1:i) == 0);
    @fact string(m.linconstr[end]) --> "x[1,1] + 4 x[2,1] + 5 x[2,2] + 7 x[3,1] + 8 x[3,2] + 9 x[3,3] $eq 0"

    @constraint(m, sum( 0*x[i,1] for i=1:3) == 0)
    @fact string(m.linconstr[end]) --> "0 $eq 0"

    @constraint(m, sum( 0*x[i,1] + y for i=1:3) == 0)
    @fact string(m.linconstr[end]) --> "3 y $eq 0"

    @fact isexpr(macroexpand(:(@constraint(m, sum( 0*x[i,1] + y for i=1:3 for j in 1:3) == 0))),:error) --> true

end"""); end

facts("[macros] Problem modification") do
    m = Model()
    @variable(m, x[1:3,1:3])
    C = [1 2 3; 4 5 6; 7 8 9]

    @constraint(m, sum{ x[i,j]*(C[i,j]-1), i in 1:3, j = 1:3; i != j} == 0)
    @fact string(m.linconstr[end]) --> "x[1,2] + 2 x[1,3] + 3 x[2,1] + 5 x[2,3] + 6 x[3,1] + 7 x[3,2] $eq 0"

    con = @constraint(m, sum{ C[i,j]*x[i,j], i = 1:3, j = 1:3; i != j} == 0)
    @fact string(m.linconstr[end]) --> "2 x[1,2] + 3 x[1,3] + 4 x[2,1] + 6 x[2,3] + 7 x[3,1] + 8 x[3,2] $eq 0"

    @variable(m, y, objective = 0, inconstraints = [con], coefficients = [-1.0])
    @fact string(m.linconstr[end]) --> "2 x[1,2] + 3 x[1,3] + 4 x[2,1] + 6 x[2,3] + 7 x[3,1] + 8 x[3,2] - y $eq 0"

    JuMP.setRHS(con, 3)
    @fact string(m.linconstr[end]) --> "2 x[1,2] + 3 x[1,3] + 4 x[2,1] + 6 x[2,3] + 7 x[3,1] + 8 x[3,2] - y $eq 3"
end

facts("[macros] Using pre-built affine is OK in macro") do
    m = Model()
    @variable(m, x)
    @variable(m, y)
    temp = x + 2y + 1
    @constraint(m, 3*temp - x - 2 >= 0)
    @fact string(m.linconstr[end]) --> "6 y + 2 x $geq -1"
    # More complex expression
    a = 1.0*x
    @constraint(m, (2+2)*((3+4)*(1+a)) == 0)
    @fact string(m.linconstr[end]) --> "28 x $eq -28"
    @fact string(a) --> "x"

    @fact string(@LinearConstraint(1 + 0*temp == 0)) --> "0 $eq -1"
end

facts("[macros] Test ranges in @variable") do
    m = Model()
    @variable(m, x[1:5])
    @variable(m, y[3:2:9])
    @variable(m, z[4:3:8])
    @variable(m, w[6:5])

    @fact length(x) --> 5
    @fact length(y) --> 4
    @fact length(z) --> 2
    @fact length(w) --> 0

    @fact x[end].col --> x[5].col
    @fact y[3].m --> y[5].m
    @fact y[3].m --> y[7].m
    @fact y[3].m --> y[9].m
    @fact_throws z[8].col  # KeyError
    @fact_throws w[end]  # BoundsError
end

facts("[macros] Unicode comparisons") do
    m = Model()
    @variable(m, 0 ≤ x ≤ 1)
    @variable(m, y ≥ 2)
    @variable(m, z ≤ 3)
    @fact m.colUpper --> [1.0, Inf,  3.0]
    @fact m.colLower --> [0.0, 2.0, -Inf]
    @constraint(m, 0 ≤ x + y ≤ 1)
    @constraint(m, x + z ≤ 2)
    @constraint(m, y + z ≥ 3)
    @constraint(m, y*z ≤ 1)
    @fact m.linconstr[1].lb --> 0.0
    @fact m.linconstr[1].ub --> 1.0
    @fact m.linconstr[2].lb --> -Inf
    @fact m.linconstr[2].ub --> 2.0
    @fact m.linconstr[3].lb --> 3.0
    @fact m.linconstr[3].ub --> Inf
    @fact m.quadconstr[1].sense --> :(<=)
end

facts("[macros] Three argument @constraint") do
    m = Model()
    @variable(m, x[1:5])
    @variable(m, y[2:2:6])

    @constraint(m, c, x[4] - y[4] == 1)
    @fact string(m.linconstr[c.idx]) --> "x[4] - y[4] $eq 1"

    @constraint(m, d[i in 1:5,j=6:-2:2], x[i] - y[j] == 2)
    @fact string(m.linconstr[d[4,4].idx]) --> "x[4] - y[4] $eq 2"

    @constraint(m, q[i=1:5], x[i]^2 == 1)
    @fact string(m.quadconstr[q[5].idx]) --> "x[5]² - 1 $eq 0"
end

facts("[macros] @constraints") do
    m = Model()
    @variable(m, x)
    @variable(m, y[1:3])

    @constraints(m, begin
        x + y[1] == 1
        ref[i=1:3], y[1] + y[i] >= i
    end)

    @fact string(m.linconstr[1]) --> "x + y[1] $eq 1"
    @fact string(m.linconstr[2]) --> "2 y[1] $geq 1"
    @fact string(m.linconstr[3]) --> "y[1] + y[2] $geq 2"
    @fact string(m.linconstr[4]) --> "y[1] + y[3] $geq 3"
end

facts("[macros] @NLconstraints") do
    m = Model()
    @variable(m, 0 <= x <= 1)
    @variable(m, y[1:3])
    @objective(m, Max, x)

    @NLconstraints(m, begin
        ref[i=1:3], y[i] == 0
        x + y[1] * y[2] * y[3] <= 0.5
    end)

    @fact length(m.nlpdata.nlconstr) --> 4
    d = JuMP.NLPEvaluator(m)
    MathProgBase.initialize(d, [:ExprGraph])

    @fact MathProgBase.constr_expr(d,1) --> :(x[2] - 0.0 == 0.0)
    @fact MathProgBase.constr_expr(d,2) --> :(x[3] - 0.0 == 0.0)
    @fact MathProgBase.constr_expr(d,3) --> :(x[4] - 0.0 == 0.0)
    @fact MathProgBase.constr_expr(d,4) --> :((x[1] + x[2] * x[3] * x[4]) - 0.5 <= 0.0)

end

facts("[macros] Vectors in nonlinear expressions") do
    m = Model()
    @variable(m, x[1:3])
    @fact_throws ErrorException @NLobjective(m, Min, x)
    @fact_throws ErrorException @NLobjective(m, Min, [1,2,3])
end

facts("[macros] @objective with quadratic") do
    m = Model()
    @variable(m, x[1:5])
    @objective(m, Max, sum{i*x[i]*x[j], i=1:5, j=5:-1:1; isodd(i) && iseven(j)} + 2x[5])

    @fact string(m.obj) --> "x[1]*x[2] + 3 x[2]*x[3] + x[1]*x[4] + 3 x[3]*x[4] + 5 x[2]*x[5] + 5 x[4]*x[5] + 2 x[5]"
end

facts("[macros] @constraint with quadratic") do
    m = Model()
    @variable(m, x[1:5])

    @constraint(m, x[3]*x[1] + sum{x[i]*x[5-i+1], i=1:5; 2 <= i <= 4} + 4x[5] == 1)
    @fact string(m.quadconstr[end]) --> "x[1]*x[3] + x[3]² + 2 x[2]*x[4] + 4 x[5] - 1 $eq 0"

    @constraint(m, sum{sum{(x[i] - 2)*x[j],j=4:5},i=2:3} >= -3*x[2]*2*x[4])
    @fact string(m.quadconstr[end]) --> "7 x[2]*x[4] + x[3]*x[4] + x[2]*x[5] + x[3]*x[5] - 4 x[4] - 4 x[5] $geq 0"

    foo(x) = x
    @constraint(m, x[1] ≤ foo(x[1])^2)
    @fact string(m.quadconstr[end]) --> "-x[1]² + x[1] $leq 0"

    @constraint(m, sum{x[i],i=1:2}*sum{x[i],i=2:3} >= 0)
    @fact string(m.quadconstr[end]) --> "x[1]*x[2] + x[2]² + x[1]*x[3] + x[2]*x[3] $geq 0"
    @constraint(m, x[1]^2 + x[2]*x[3] >= 0)
    @fact string(m.quadconstr[end]) --> "x[1]² + x[2]*x[3] $geq 0"
    @constraint(m, x[1]^2 + (x[2]+3)*(x[3]-1) >= 0)
    @fact string(m.quadconstr[end]) --> "x[1]² + x[2]*x[3] + 3 x[3] - x[2] - 3 $geq 0"
    @constraint(m, sum{x[i],i=1:2}^2 >= 0)
    @fact string(m.quadconstr[end]) --> "x[1]² + 2 x[1]*x[2] + x[2]² $geq 0"

    myquadexpr = x[1]*x[2]
    @constraint(m, sum{i*myquadexpr + x[i], i=1:3} + sum{x[i] + myquadexpr*i, i=1:3} == 0)
    @fact string(m.quadconstr[end]) --> "12 x[1]*x[2] + 2 x[1] + 2 x[2] + 2 x[3] $eq 0"

    @constraint(m, (x[1] + x[2])*sum{ 0*x[i] + x[3], i=1:3} == 0)
    @fact string(m.quadconstr[end]) --> "3 x[1]*x[3] + 3 x[2]*x[3] $eq 0"

    @fact string(@QuadConstraint(1 + 0*myquadexpr == 0)) --> "1 $eq 0"

    @variable(m, y)
    @fact string(@QuadConstraint(1 + (2y)*y   == 0)) --> "2 y² + 1 $eq 0"
    @fact string(@QuadConstraint(1 +   y *y*2 == 0)) --> "2 y² + 1 $eq 0"
    z = 2y
    @fact string(@QuadConstraint(y*y + y*z == 0)) --> "3 y² $eq 0"
end

facts("[macros] Triangular indexing, iteration") do
    n = 10
    trimod = Model()
    @variable(trimod, x[i=1:n,j=i:n])
    @variable(trimod, y[i=3:2:7,j=-i])
    @fact MathProgBase.numvar(trimod) --> n*(n+1)/2 + 3
    S = Any[(i,i+2) for i in 1:5]
    @variable(trimod, z[(i,j)=S,k=i:j])
    @fact length(z.tupledict) --> 15
    @constraint(trimod, cref[i=1:n,j=i:n], x[i,j] + y[5,-5] == 1)
    @fact MathProgBase.numconstr(trimod) --> n*(n+1)/2

    cntr = zeros(Bool, n, n)
    for ((i,j),var) in zip(keys(x),values(x))
        @fact x[i,j] --> exactly(var)
        cntr[i,j] = true
    end
    for i in 1:n, j in 1:n
        @fact cntr[i,j] --> (j >= i)
    end
end

facts("[macros] Multidimensional indexing") do
    model = Model()
    I1 = 1:5
    I2 = 2:8
    I3 = 5:6
    @variable(model, x[1:5,2:8,5:6])
    coll = Int[]
    for v in values(x)
        push!(coll, v.col)
    end
    p = 1
    match = true
    for v in x.innerArray
        match &= (coll[p] == v.col)
        p += 1
    end
    @fact match --> true
end

facts("[macros] @expression") do
    model = Model()
    @variable(model, x[1:3,1:3])
    @expression(model, expr, sum{i*x[i,j] + j, i=1:3,j in 1:3})
    @fact string(expr) --> "x[1,1] + x[1,2] + x[1,3] + 2 x[2,1] + 2 x[2,2] + 2 x[2,3] + 3 x[3,1] + 3 x[3,2] + 3 x[3,3] + 18"

    @fact_throws @expression(model, blah[i=1:3], x[i,1]^2)

    @expression(model, y[i=1:2], sum{x[i,1]; i == 1})
    @fact string(y[1]) --> "x[1,1]"
    @fact string(y[2]) --> "0"

    # deprecated versions
    @expression(expr2, sum{i*x[i,j] + j, i=1:3,j in 1:3})
    @fact string(expr2) --> "x[1,1] + x[1,2] + x[1,3] + 2 x[2,1] + 2 x[2,2] + 2 x[2,3] + 3 x[3,1] + 3 x[3,2] + 3 x[3,3] + 18"
    expr2 = @expression(sum{i*x[i,j] + j, i=1:3,j in 1:3})
    @fact string(expr2) --> "x[1,1] + x[1,2] + x[1,3] + 2 x[2,1] + 2 x[2,2] + 2 x[2,3] + 3 x[3,1] + 3 x[3,2] + 3 x[3,3] + 18"
end

facts("[macros] Conditions in constraint indexing") do
    model = Model()
    @variable(model, x[1:10])
    @constraint(model, c1[i=1:9;isodd(i)], x[i] + x[i+1] <= 1)
    @NLconstraint(model, c2[i=["red","blue","green"], k=9:-2:2; (i == "red" && isodd(k)) || (k >=4 && (i == "blue" || i == "green"))], x[k]^3 <= 1)
    @fact length(model.linconstr) --> 5
    @fact string(model.linconstr[1]) --> "x[1] + x[2] $leq 1"
    @fact string(model.linconstr[2]) --> "x[3] + x[4] $leq 1"
    @fact string(model.linconstr[3]) --> "x[5] + x[6] $leq 1"
    @fact string(model.linconstr[4]) --> "x[7] + x[8] $leq 1"
    @fact string(model.linconstr[5]) --> "x[9] + x[10] $leq 1"
    @fact length(model.nlpdata.nlconstr) --> 10
end

facts("[macros] Test changes in condition parsing") do
    ex = :(x[12;3])
    @fact ex.head --> :typed_vcat
    @fact ex.args --> [:x, 12, 3]

    ex = :(x[i=1:3,j=S;isodd(i) && i+j>=2])
    @fact ex.head --> :typed_vcat
    @fact ex.args --> [:x,
                      Expr(:parameters, Expr(:&&, :(isodd(i)), :(i+j>=2))),
                      Expr(:(=), :i, :(1:3)),
                      Expr(:(=), :j, :S)]
end

facts("[macros] Norm parsing") do
    model = Model()
    @variable(model, x[1:2,1:2])
    @constraint(model, -2norm2{x[i,j], i in 1:2, j=1:2} + x[1,2] >= -1)
    @constraint(model, -2norm2{x[i,j], i=1:2, j in 1:2; iseven(i+j)} + x[1,2] >= -1)
    @constraint(model, 1 >= 2*norm2{x[i,1], i in 1:2})
    @fact string(model.socconstr[1]) --> "2.0 $Vert[x[1,1],x[1,2],x[2,1],x[2,2]]$Vert$sub2 $leq x[1,2] + 1"
    @fact string(model.socconstr[2]) --> "2.0 $Vert[x[1,1],x[2,2]]$Vert$sub2 $leq x[1,2] + 1"
    @fact string(model.socconstr[3]) --> "2.0 $Vert[x[1,1],x[2,1]]$Vert$sub2 $leq 1"
    @fact_throws @constraint(model, (x[1,1]+1)*norm2{x[i,j], i=1:2, j=1:2} + x[1,2] >= -1)
    @fact_throws @constraint(model, norm2{x[i,j], i=1:2, j=1:2} + x[1,2] >= -1)
end

if VERSION >= v"0.5-dev+5475"
eval("""
facts("[macros] Generator norm parsing") do
    model = Model()
    @variable(model, x[1:2,1:2])
    @constraint(model, -2norm2(x[i,j] for i in 1:2, j=1:2) + x[1,2] >= -1)
    @constraint(model, -2norm2(x[i,j] for i=1:2, j in 1:2 if iseven(i+j)) + x[1,2] >= -1)
    @constraint(model, 1 >= 2*norm2(x[i,1] for i in 1:2))
    @fact string(model.socconstr[1]) --> "2.0 $Vert[x[1,1],x[1,2],x[2,1],x[2,2]]$Vert$sub2 $leq x[1,2] + 1"
    @fact string(model.socconstr[2]) --> "2.0 $Vert[x[1,1],x[2,2]]$Vert$sub2 $leq x[1,2] + 1"
    @fact string(model.socconstr[3]) --> "2.0 $Vert[x[1,1],x[2,1]]$Vert$sub2 $leq 1"
    @fact_throws @constraint(model, (x[1,1]+1)*norm2(x[i,j] for i=1:2, j=1:2) + x[1,2] >= -1)
    @fact_throws @constraint(model, norm2(x[i,j] for i=1:2, j=1:2) + x[1,2] >= -1)
end""");end

facts("[macros] Extraneous terms in QuadExpr (#535)") do
    model = Model()
    @variable(model, x)
    @variable(model, y)
    @constraint(model, x*x <= y*y)
    @fact string(model.quadconstr[1]) --> "x² - y² $leq 0"
end

facts("[macros] Special-case binary multiplication in addtoexpr_reorder (#537)") do
    dual = Model()
    @variable(dual, α[1:3] <= 0)
    @variable(dual, γ >= 0)
    @constraint(dual, ones(3)*γ .<= α)
    @fact string(dual.linconstr[1]) --> "γ - α[1] $leq 0"
    @fact string(dual.linconstr[2]) --> "γ - α[2] $leq 0"
    @fact string(dual.linconstr[3]) --> "γ - α[3] $leq 0"
end

facts("[macros] Indices in macros don't leak out of scope (#582)") do
    m = Model()
    cnt = 4
    for i in 5:8
        @variable(m, x[i=1:3,j=1:3] ≤ i)
        cnt += 1
        @fact i --> cnt
    end
    cnt = 4
    for i in 5:8
        @variable(m, y[i=2:4,j=1:3] ≤ i)
        cnt += 1
        @fact i --> cnt
    end
    cnt = 4
    for i in 5:8
        @variable(m, z[i=[1:3;],j=1:3] ≤ i)
        cnt += 1
        @fact i --> cnt
    end
    @fact m.colUpper --> vcat(repeat([1.0,2.0,3.0], inner=[3], outer=[4]),
                              repeat([2.0,3.0,4.0], inner=[3], outer=[4]),
                              repeat([1.0,2.0,3.0], inner=[3], outer=[4]))
    @variable(m, x[i=1:3] ≤ i)
    cnt = 4
    for i in 5:8
        @constraint(m, sum{x[i], i=1, j=1, k=1} == 1)
        cnt += 1
        @fact i --> cnt
    end
    cnt = 4
    for i in 5:8
        @constraint(m, norm2{x[i], i=1, j=1, k=1} <= 1)
        cnt += 1
        @fact i --> cnt
    end
    cnt = 4
    for i in 5:8
        @constraint(m, c[i=1:3,j=1:3], x[i] == 1)
        cnt += 1
        @fact i --> cnt
    end
end

facts("[macros] Issue #621") do
    m = Model()
    @variable(m, x)

    q = x^2
    a = x+1
    con = @QuadConstraint(a+q*3 <= 0)
    @fact string(con) --> "3 x² + x + 1 $leq 0"
end

facts("[macros] @variables and @constraints") do
    m = Model()
    @variables m begin
        0 ≤ x[i=1:2] ≤ i
        y ≥ 2, Int, (start = 0.7)
        z ≤ 3, (start=10)
        q, (Bin, start=0.5)
    end
    @constraints m begin
        0 ≤ x[1] + y ≤ 1
        x[1] + z ≤ 2
        y + z ≥ 3
        y*z ≤ 1
    end
    @fact m.colUpper --> [1.0, 2.0, Inf,  3.0, 1.0]
    @fact m.colLower --> [0.0, 0.0, 2.0, -Inf, 0.0]
    @fact m.colCat --> [:Cont, :Cont, :Int, :Cont, :Bin]
    @fact getvalue(y) --> 0.7
    @fact getvalue(z) --> 10
    @fact getvalue(q) --> 0.5
    @fact m.linconstr[1].lb --> 0.0
    @fact m.linconstr[1].ub --> 1.0
    @fact m.linconstr[2].lb --> -Inf
    @fact m.linconstr[2].ub --> 2.0
    @fact m.linconstr[3].lb --> 3.0
    @fact m.linconstr[3].ub --> Inf
    @fact m.quadconstr[1].sense --> :(<=)
end

facts("[macros] @expressions and @NLexpressions") do
    m = Model()
    @variable(m, x)

    @expressions(m, begin
        myex[i=1:2], x + i
        myex2,       x + 3
    end)
    setvalue(x, 1)
    @fact getvalue(myex[1]) --> 2
    @fact getvalue(myex[2]) --> 3
    @fact getvalue(myex2)   --> 4
    @fact_throws getvalue(myex[3])

    @NLexpressions(m, begin
        nlex,         x^2
        nlex2[i=1:2], x^i + sqrt(x)
    end)
    setvalue(x, 2)
    @fact getvalue(nlex) --> 4
    @fact getvalue(nlex2[1]) --> 2^1 + sqrt(2)
    @fact getvalue(nlex2[2]) --> 2^2 + sqrt(2)
    @fact_throws getvalue(nlex2[3])
end

facts("[macros] No bare symbols in constraint macros") do
    m = Model()
    @variable(m, x)
    @fact macroexpand(:(@constraint(m, x))).head --> :error
    @fact macroexpand(:(@constraint(m, :foo))).head --> :error
    @fact macroexpand(:(@SDconstraint(m, x))).head --> :error
    @fact macroexpand(:(@SDconstraint(m, :foo))).head --> :error
    @fact macroexpand(:(@NLconstraint(m, x))).head --> :error
    @fact macroexpand(:(@NLconstraint(m, :foo))).head --> :error
end

facts("[macros] LB/UB kwargs") do
    m = Model()
    @variable(m, a, lowerbound=0)
    @variable(m, b, lowerbound=0, upperbound=1)
    @variable(m, c, upperbound=1, lowerbound=0)
    @variable(m, d, upperbound=1)
    @variable(m, e >= 0, upperbound=1)
    @variable(m, f <= 1, lowerbound=0)
    @fact m.colLower --> [0,0,0,-Inf,0,0]
    @fact m.colUpper --> [Inf,1,1,1,1,1]

    @fact macroexpand(:(@variable(m, g >= 0, lowerbound=1))).head --> :error
    @fact macroexpand(:(@variable(m, h <= 1, upperbound=1))).head --> :error
    @fact macroexpand(:(@variable(m, 0 <= i <= 1, lowerbound=1))).head --> :error
    @fact macroexpand(:(@variable(m, 0 <= j <= 1, upperbound=1))).head --> :error
end

facts("[macros] Anonymous versions of macros") do
    m = Model()
    x = @variable(m, [1:3], lowerbound=1.0)
    y = @variable(m, [2:3,1:3], upperbound=1.0, lowerbound=0.0, Bin)
    z = @variable(m, [[:red,:blue]], Int)
    w = @variable(m, [i=1:4,j=1:4;isodd(i+j)], SemiCont)
    # v = @variable(m, [i=1:3,j=1:3], Symmetric, lowerbound = eye(3)[i,j])
    u = @variable(m, [1:4,1:4], SDP)
    @fact macroexpand(:(@variable(m, [1:3] <= 1))).head --> :error

    @fact getlowerbound(x[1]) --> 1.0
    @fact getupperbound(x[1]) --> Inf
    @fact getcategory(x[1]) --> :Cont
    @fact getlowerbound(y[2,1]) --> 0.0
    @fact getupperbound(y[2,1]) --> 1.0
    @fact getcategory(y[2,1]) --> :Bin
    @fact getlowerbound(z[:red]) --> -Inf
    @fact getupperbound(z[:red]) --> Inf
    @fact getcategory(z[:red]) --> :Int
    @fact getlowerbound(w[1,2]) --> -Inf
    @fact getupperbound(w[1,2]) --> Inf
    @fact getcategory(w[1,2]) --> :SemiCont
    # @fact getlowerbound(v[1,2]) --> 0.0
    # @fact getupperbound(v[1,2]) --> Inf
    # @fact getcategory(v[1,2]) --> :Cont
    @fact getlowerbound(u[1,2]) --> -Inf
    @fact getupperbound(u[1,2]) --> Inf
    @fact getcategory(u[1,2]) --> :Cont

    c = @constraint(m, [i=1:3], x[i] <= z[:red])
    d = @NLconstraint(m, [i=1:3], x[i]^3 == 1)
    e = @NLexpression(m, [i=2:3], y[2,i]^3)
    f = @expression(m, [i=[:red,:blue]], u[1,2] + 2z[i])

    # not sure how else to test this
    @fact c[1] --> c[1]
    @fact d[1] --> d[1]
    @fact e[2] --> e[2]
    @fact f[:red] --> f[:red]
end

facts("[macros] Colons in index sets") do
    m = Model()
    S = [:]
    @fact_throws ErrorException @variable(m, x[S])
end

facts("[macros] getconstraint") do
    m = Model()
    @variable(m, x)
    _c1 = @constraint(m, c1, x == 1)
    _c2 = @constraint(m, c2[i=1:3], x == i)
    _c3 = @NLconstraint(m, c3, x^3 == 1)
    _c4 = @NLconstraint(m, c4[i=1:3], x^i == 1)
    @fact _c1 --> getconstraint(m, :c1)
    @fact _c2 --> getconstraint(m, :c2)
    @fact _c3 --> getconstraint(m, :c3)
    @fact _c4 --> getconstraint(m, :c4)
end

facts("[macros] Anonymous singleton variables") do
    m = Model()
    x = @variable(m)
    y = @variable(m, lowerbound=0, upperbound=1)
    @fact x --> Variable(m, 1)
    @fact y --> Variable(m, 2)
end


text:
function LoveNumbers(arg0::Vector{Vector{jdouble}}, arg1::Vector{Vector{jdouble}}, arg2::Vector{Vector{jdouble}})
    return LoveNumbers((Vector{Vector{jdouble}}, Vector{Vector{jdouble}}, Vector{Vector{jdouble}}), arg0, arg1, arg2)
end

function get_imaginary(obj::LoveNumbers, arg0::jint, arg1::jint)
    return jcall(obj, "getImaginary", jdouble, (jint, jint), arg0, arg1)
end

function get_plus(obj::LoveNumbers, arg0::jint, arg1::jint)
    return jcall(obj, "getPlus", jdouble, (jint, jint), arg0, arg1)
end

function get_real(obj::LoveNumbers, arg0::jint, arg1::jint)
    return jcall(obj, "getReal", jdouble, (jint, jint), arg0, arg1)
end

function get_size(obj::LoveNumbers)
    return jcall(obj, "getSize", jint, ())
end


text:
"""
    Script to generate the scatter plot associated with:
    - benchmark_randV_randE.jl
"""
using Pkg
Pkg.activate(".")
using LTMSim
using SimpleHypergraphs
using PyPlot
using Serialization
using Statistics
using LaTeXStrings

project_path = dirname(pathof(LTMSim))

# y-values = Time (s)
fname = "benchmark_randV_randE.data"
data_path = joinpath(project_path, "..", "res", "journal", fname)

# x-values = target set size
fname_tss = "randV_randE_with_noOpt.data"
data_path_tss = joinpath(project_path, "..", "res", "journal", fname_tss)

res_path = joinpath(project_path, "..", "res", "journal", "plot", "benchmark")

hg_files = readdir(joinpath(project_path, "..", "data", "hgs"))
hg_names = [split(file, ".")[1] for file in hg_files]

hg_names[8] = "music-rev"
hg_names[10] = "restaurants-rev"
hg_names[11] = "bars-rev"

#
# DATA
#
data = deserialize(data_path)
data_tss = deserialize(data_path_tss)
hgs = [prune_hypergraph!(hg_load(joinpath(project_path, "..", "data", "hgs", hg_file))) for hg_file in hg_files]

#
# PLOTTING INFO
#
algorithms = [
    "BinarySearch(H)",
    "BinarySearch(H)-noOpt",
    "Greedy(H)",
    "Greedy(H)-noOpt",
    "Greedy([H]₂)",
    "Greedy([H]₂)-noOpt",
    "SubTSS(H)",
    "SubTSS(H)-noOpt"
]

labels_dict = Dict{String, String}(
    "BinarySearch(H)" => "StaticGreedy",
    "BinarySearch(H)-noOpt" => "StaticGreedy-noOpt",
    "Greedy(H)" => "DynamicGreedy",
    "Greedy(H)-noOpt" => "DynamicGreedy-noOpt",
    "Greedy([H]₂)" => L"DynamicGreedy_{[H]_2}",
    "Greedy([H]₂)-noOpt" => L"DynamicGreedy_{[H]_2}-noOpt",
    "SubTSS(H)" => "SubTSS",
    "SubTSS(H)-noOpt" => "SubTSS-noOpt"
)

#
# SCATTER PLOT
#
colorz = ["#2C7BB6", "#D7191C", "#FF8900", "#33CC33", "#cc0099", "#4d4dff", "#008080", "#2C7BB6"]
markers = ["o", "^", "*", "d", "p", ".", "h", "x"]

for (hg_index, hg_name) in enumerate(hg_names)
    println(hg_name)

    labels = []

    clf()
    plt.figure(figsize=(7,5))

    for (algo_index, algo) in enumerate(algorithms)
        println(algo_index, " ", algo, " ")

        # get execution time (mean)
        y = mean(data[algo][hg_index])

        # get target set size (%)
        x = mean(data_tss[algo][hg_index] / nhv(hgs[hg_index]))

        push!(
            labels,
            labels_dict[algo]
        )

        plt.scatter(x, y, c=colorz[algo_index], marker=markers[algo_index], label=labels_dict[algo])
    end

    plt.legend()

    plt.title("$hg_name")

    plt.xticks(fontsize="x-large") #range(0, length(ticks) * 2, step=2),
    plt.yticks(fontsize="x-large")

    plt.xlim(0, 1)
    xlabel("Seed set size / n", fontstyle = "italic", fontsize="large", labelpad=10) #, fontweight="semibold"

    plt.ylim(bottom=0)
    ylabel("Execution time (s)", fontstyle = "italic", fontsize="large", labelpad=10) #fontweight="semibold",

    plt.tight_layout()

    PyPlot.savefig(joinpath(res_path, "$hg_name.png"))
end


text:
@testset "Flux variability summary" begin
    model = load_model(model_paths["e_coli_core.json"])

    sol = flux_variability_analysis_dict(
        model,
        Tulip.Optimizer;
        bounds = objective_bounds(0.90),
        modifications = [change_optimizer_attribute("IPM_IterationsLimit", 2000)],
    )

    fr = flux_variability_summary(sol)
    @test isapprox(
        fr.biomass_fluxes["BIOMASS_Ecoli_core_w_GAM"][1],
        0.7865293520891825;
        atol = TEST_TOLERANCE,
    )
    @test isapprox(
        fr.exchange_fluxes["EX_for_e"][2],
        11.322324494491848;
        atol = TEST_TOLERANCE,
    )
end


text:
# Unit testing for element-wise expressions on scalar backend

import Devectorize
import Devectorize.@devec
import Devectorize.@inspect_devec
import Devectorize.dump_devec
import Devectorize.sqr
import Devectorize.blend

using Base.Test

a = [1., 2., 3., 4., 5., 6., 7., 8.]
b = [3., 4., 5., 6., 8., 7., 6., 5.]
c = [9., 8., 7., 6., 4., 2., 3., 1.]
cv = 12.0
cv_a = [12., 12., 12., 12., 12., 12., 12., 12.]

abc = [a b c]

#################################################
#
#   basic element-wise maps
#
#################################################

# unary

@devec r = -a
@test isequal(r, -a)

r[:] = 0.
@devec r[:] = -a
@test isequal(r, -a)

# binary

@devec r = a + b
@test isequal(r, a + b)

r[:] = 0.
@devec r[:] = a + b
@test isequal(r, a + b)

@devec r = a + 1.
@test isequal(r, a + 1.)

r[:] = 0.
@devec r[:] = a + 1.
@test isequal(r, a + 1.)

@devec r = 1. + b
@test isequal(r, 1. + b)

r[:] = 0.
@devec r[:] = 1. + b
@test isequal(r, 1. + b)

@devec r[:] = 1. + 2.
@test isequal(r, fill(3., size(a)))

# ternary

@devec r = +(a, b, c)
@test isequal(r, a + b + c)

r[:] = 0.
@devec r[:] = +(a, b, c)
@test isequal(r, a + b + c)

@devec r = +(1., b, c)
@test isequal(r, 1. + b + c)

r[:] = 0.
@devec r[:] = +(1., b, c)
@test isequal(r, 1. + b + c)

@devec r = +(a, 2., c)
@test isequal(r, a + 2. + c)

r[:] = 0.
@devec r[:] = +(a, 2., c)
@test isequal(r, a + 2. + c)

@devec r = +(a, b, 3.)
@test isequal(r, a + b + 3.)

r[:] = 0.
@devec r[:] = +(a, b, 3.)
@test isequal(r, a + b + 3.)

@devec r = +(1., 2., c)
@test isequal(r, 1. + 2. + c)

r[:] = 0.
@devec r[:] = +(1., 2., c)
@test isequal(r, 1. + 2. + c)

@devec r = +(1., b, 3.)
@test isequal(r, 1. + b + 3.)

r[:] = 0.
@devec r[:] = +(1., b, 3.)
@test isequal(r, 1. + b + 3.)

@devec r = +(a, 2., 3.)
@test isequal(r, a + 2. + 3.)

r[:] = 0.
@devec r[:] = +(a, 2., 3.)
@test isequal(r, a + 2. + 3.)

@devec r[:] = +(1., 2., 3.)
@test isequal(r, fill(6., size(a)))

# comparison and blend

@devec r = blend(a + 2 .>= b, a + 2, b)
@test isequal(r, max(a+2, b))

r[:] = 0.
@devec r[:] = blend(a + 2 .>= b, a + 2, b)
@test isequal(r, max(a+2, b))

@devec r = blend(a + 2 .<= b, a + 2, b)
@test isequal(r, min(a + 2, b))

r[:] = 0.
@devec r[:] = blend(a + 2 .<= b, a + 2, b)
@test isequal(r, min(a + 2, b))

#################################################
#
#   interaction with qualified names
#
#################################################

type S
    a::Array
    b::Array
    c::Array
end

s = S(a, b, c)

@devec r = sqr(s.a)
@test isequal(r, sqr(s.a))

@devec r = s.a + s.b .* s.c
@test isequal(r, (s.a + s.b .* s.c))

@devec r = +(s.a, s.b, s.c)
@test isequal(r, s.a + s.b + s.c)

#################################################
#
#   interaction with references
#
#################################################

@devec r = sqr(abc[:,1])
@test isequal(r, sqr(a))

@devec r = sqr(abc[:,:])
@test isequal(r, sqr(abc))

@devec r = abc[:,:] .* cv
@test isequal(r, abc .* cv)

j = 3
@devec r = abc[:,1] + b[:] .* abc[:,j]
@test isequal(r, a + b .* c)

r = zeros(8, 3)
rr = zeros(size(r))
@devec r[:,2:3] = abc[:, 1:2] + cv
rr[:,2:3] = abc[:, 1:2] + cv
@test isequal(r, rr)

#################################################
#
#   op-assignment
#
#################################################

r = a
@devec r += b
@test isequal(r, a + b)

r0 = r
@devec r[:] += c
@test r === r0
@test isequal(r, a + b + c)

r = a
@devec r .*= 2.
@test isequal(r, a * 2)

@devec r -= (c + 3)
@test isequal(r, a * 2 - (c + 3))

#################################################
#
#   blocks
#
#################################################

@devec begin
    x = a + b .* c
    y = sqr(a) + c
end
@test isequal(x, a + b .* c)
@test isequal(y, sqr(a) + c)


text:
using ArgParse
using NBody

integrators = Dict(
  "euler"     => NBody.forward_step,
  "leapfrog"  => NBody.leapfrog_step,
  "rk2"       => NBody.rk2_step,
  "rk4"       => NBody.rk4_step,
  "hermite"   => NBody.hermite_step
)

function parse_commandline()
  s = ArgParseSettings()

  @add_arg_table s begin
    "--total_duration", "-t"
      help = "Duration of the integration"
      arg_type = Float64
      default = 10.0

    "--output_interval", "-o"
      help = "Snapshot output interval"
      arg_type = Float64
      default = 1.0

    "--diagnostics_interval", "-e"
      help = "Diagnostics output interval"
      arg_type = Float64
      default = 1.0

    "--step_size", "-d"
      help = "Integration time step"
      arg_type = Float64
      default = 0.001

    "--step_size_control", "-c"
      help = "Integration time step scaling control"
      arg_type = Float64
      default = 0.01

    "--integrator", "-i"
      help = "Integration method"
      arg_type = String
      default = "rk4"
      range_tester = x -> haskey(integrators, x)

    "--dynamic_step", "-q"
      help = "Enable dynamic step size"
      action = :store_true

    "--exact_time"
      help = "Force output to occur at exact times"
      action = :store_true

    "--extra_diagnostics", "-x"
      help = "Extra diagnostics"
      action = :store_true

    "--softening_length", "-s"
      help = "Softening ength"
      arg_type = Float64
      default = 0.0

    "--no_init_out", "-n"
      help = "Don't output initial snapshot"
      action = :store_false
  end

  return parse_args(s)
end

function main(args)

  nb = read_nbody()

  # Should get moved in with rest of args
  nb.soften_len = args["softening_length"]

  method = integrators[args["integrator"]]

  final_snapshot = args["output_interval"] > args["total_duration"]

  ea = EvolveArgs(
    method,
    args["softening_length"],
    args["step_size"],
    args["step_size_control"],
    args["total_duration"],
    args["output_interval"],
    args["diagnostics_interval"],
    !args["no_init_out"],
    args["extra_diagnostics"],
    args["dynamic_step"],
    args["exact_time"],
    final_snapshot
  )

  write_info(args, ea)

  evolve(nb, ea)
end

function write_info(args, ea)
  println(STDERR, "A simple N-Body code")

  d = Dict(
    "Integration method"              => ea.integ_method,
    "Integration time step"           => ea.dt,
    "Diagnostics output interval"     => ea.dt_stats,
    "Snapshot output interval"        => ea.dt_output,
    "Duration of the integration"     => ea.time_end,
    "Softening length"                => ea.soften_len
    )
  for (k,v) in d
    @printf(STDERR, "%-30s %-20s\n", k, v)
  end
end

parsed_args = parse_commandline()

main(parsed_args)


text:
using Logging

# Enable @debug logging
function set_log_level(level)
    global_logger(ConsoleLogger(stderr, level))
end

log_debug() = set_log_level(Logging.Debug)
log_info()  = set_log_level(Logging.Info)

"""
    interpolate(values, ts=10)

Linearly interpolate annual values from a vector of values
for timestep of length `ts`.
"""
function interpolate(oldvalues::Vector{T}, ts::Int=10) where T <: Union{Float64, Int}
    count = length(oldvalues)
    newvalues::Vector{Float64} = zeros((count-1) * ts + 1)
    fracs = collect(0.0:1/ts:1.0)

    for i = 1:count - 1
        first = oldvalues[i]
        last  = oldvalues[i+1]
        diff  = last - first

        first_idx = (i - 1) * ts + 1
        end_idx   = first_idx + ts - 1
        newvalues[first_idx:end_idx] = first .+ diff * fracs[1:end-1]
    end

    newvalues[end] = oldvalues[end]
    return newvalues
end

# MacroTools has a "prettify", so we have to import to "extend"
# even though our function is unrelated. This seems unfortunate.
import MacroTools.prettify

"""
    MacroTools.prettify(s::String)

Accepts a camelcase or snakecase string, and makes it human-readable
e.g. camelCase -> Camel Case; snake_case -> Snake Case
"""
function MacroTools.prettify(s::String)
    s = replace(s, r"_" => s" ")
    s = replace(s, r"([a-z])([A-Z])" =>  s"\1 \2")
    s = replace(s, r"([A-Z]+)([A-Z])" => s"\1 \2")        # handle case of consecutive caps by splitting last from rest

    # Capitalize the first letter of each word
    s_arr = split(s)

    for (i, word) in enumerate(s_arr)
        s_arr[i] = "$(uppercase(word[1]))$(word[2:length(word)])"
    end

    # Return our string
    return join(s_arr, " ")
end

prettify(s::Symbol) = prettify(string(s))


text:
export CoupledSimulation, step!, run!

"""
    CoupledSimulation

An abstract type representing a coupled simulation.
"""
abstract type CoupledSimulation end

"""
    run!(::CoupledSimulation)

A simple outer timestepping loop for coupled system runs.

This will be formalized when the run! functionality for component
models is implemented so to have a consistent interface.
"""
function run!(sim::CoupledSimulation)
    clock = sim.clock
    while !stop_time_exceeded(clock)
        step!(sim, clock.dt)
        tick!(clock)
    end
end

"""
    step!(sim, dt)

Advances a simulation by `dt`.

Note that `dt` is not necessarily the simulation's timestep length;
a simuation could take several shorter steps that total to `dt`.
"""
function step!(sim, dt) end

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine


text:
# Minimum variance analysis

export mva

"""
    mva(Bx, By, Bz; verbose=false)

Perform minimum variance analysis to vector components defined in orthogonal
coordinates `Bx`, `By` and `Bz`.
If λ₁ ≥ λ₂ ≥ λ₃ are 3 eigenvalues of the constructed matrix M, then a good
indicator of nice fitting LMN coordinate system should have λ₂/λ₃ > 5. Set
`verbose=true` to turn on the check.
"""
function mva(Bx, By, Bz; verbose=false)

   B̄1 = mean(Bx)
   B̄2 = mean(By)
   B̄3 = mean(Bz)
   B̄11= mean(Bx.*Bx) - B̄1*B̄1
   B̄22= mean(By.*By) - B̄2*B̄2
   B̄33= mean(Bz.*Bz) - B̄3*B̄3
   B̄12= mean(Bx.*By) - B̄1*B̄2
   B̄23= mean(By.*Bz) - B̄2*B̄3
   B̄31= mean(Bz.*Bx) - B̄3*B̄1
   # Construct the matrix
   M = [B̄11 B̄12 B̄31; B̄12 B̄22 B̄23; B̄31 B̄23 B̄33]

   # Compute the eigen values and ratios (descending order)
   F = eigen(M, sortby = x -> -abs(x))

   if verbose
      println("Eigenvalues:", F.values)
      println("Eigenvectors:")
      println(F.vectors[:,1])
      println(F.vectors[:,2])
      println(F.vectors[:,3])
      r = F.values[2] / F.values[3]
      println("Ratio of intermediate variance to minimum variance = ",
         round(r, digits=3))
      if r ≥ 5
         println("Seems to be a proper MVA attempt!")
      else
         @warn "Take the MVA result with a grain of salt!"
      end
   end
   F
end


text:
# using Distributions
# using BayesInference

export ldensweight_mat, getEy, qnorm_mix, getQuant, getlogdens_EY;

function ldensweight_mat(X_pred::Array{T,2}, sims::Union{Array{Dict{Symbol,Any},1}, Array{Any,1}},
    γδc::Union{Float64, Array{T, 1}, Nothing}=Inf) where T <: Real

    # fill(1.0e6, size(sims[1][:β_y])[2]) # default γδc under variance inflation method

    nsim = length(sims)
    H, K = size(sims[1][:β_y])
    npred, Kx = size(X_pred)

    K == Kx || throw(error("X_pred dimensions not aligned with simuations."))
    useγ = haskey(sims[1], :γ)
    if useγ
        γglobal = ( typeof(sims[1][:γ]) == BitArray{1} )
        γlocal = ( typeof(sims[1][:γ]) == BitArray{2} )
    end

    out = Array{T,3}(undef, nsim, npred, H)

    if K > 1
        for ii = 1:nsim

            if useγ

                if γglobal

                    if γδc == Inf  # subset method
                        γindx = findall(sims[ii][:γ])
                        nγ = length( γindx )

                        if nγ == 0
                            lNX = zeros(Float64, npred, H)
                        elseif nγ == 1
                            δγ_x = δ_x_modify_γ(sims[ii][:δ_x], sims[ii][:γ], γδc)
                            lNX = lNXmat(X_pred[:,γindx],
                                         sims[ii][:μ_x][:,γindx], δγ_x) # npred by H matrix
                        elseif nγ > 1 && !isnothing(sims[ii][:β_x]) # indicates Σx_type is :full
                            βγ_x, δγ_x = βδ_x_modify_γ(sims[ii][:β_x],
                                        sims[ii][:δ_x],sims[ii][:γ], γδc)

                            lNX = lNXmat(X_pred[:,γindx], sims[ii][:μ_x][:,γindx],
                                         βγ_x, δγ_x) # npred by H matrix
                        elseif nγ > 1 && isnothing(sims[ii][:β_x]) # indicates Σx_type is :diag
                            lNX = lNXmat_Σdiag(X_pred[:,γindx], sims[ii][:μ_x][:,γindx], sims[ii][:δ_x][:,γindx]) # npred by H matrix
                        end

                    elseif isnothing(γδc) # integration method
                        γindx = findall(sims[ii][:γ])
                        nγ = length( γindx )

                        if nγ == 0
                            lNX = zeros(Float64, npred, H)
                        elseif nγ == 1 && !isnothing(sims[ii][:β_x]) # indicates Σx_type is :full
                            σ2xs = [ sqfChol_to_Σ( [ sims[ii][:β_x][k][h,:] for k = 1:(K-1) ], sims[ii][:δ_x][h,:] ).mat[γindx, γindx][1] for h = 1:H ]
                            lNX = lNXmat(X_pred[:,γindx], sims[ii][:μ_x][:,γindx], σ2xs) # npred by H matrix
                        elseif nγ == 1 && isnothing(sims[ii][:β_x]) # indicates Σx_type is :diag
                            σ2xs = deepcopy( sims[ii][:δ_x][:,γindx] )
                            lNX = lNXmat(X_pred[:,γindx], sims[ii][:μ_x][:,γindx], σ2xs) # npred by H matrix
                        elseif nγ > 1 && !isnothing(sims[ii][:β_x]) # indicates Σx_type is :full
                            Σxs = [ PDMat( sqfChol_to_Σ( [ sims[ii][:β_x][k][h,:] for k = 1:(K-1) ], sims[ii][:δ_x][h,:] ).mat[γindx, γindx] ) for h = 1:H ]
                            lNX = hcat( [ logpdf(MultivariateNormal(sims[ii][:μ_x][h,γindx], Σxs[h]), Matrix(X_pred[:,γindx]')) for h = 1:H ]... ) # n by H matrix
                        elseif nγ > 1 && isnothing(sims[ii][:β_x]) # indicates Σx_type is :diag
                            lNX = lNXmat_Σdiag(X_pred[:,γindx], sims[ii][:μ_x][:,γindx], sims[ii][:δ_x][:,γindx]) # npred by H matrix
                        end

                    else  # variance-inflation method
                        !isnothing(sims[ii][:β_x]) || throw("Variance-inflation variable selection not implemented for diagonal Σx.")
                        βγ_x, δγ_x = βδ_x_modify_γ(sims[ii][:β_x],
                                    sims[ii][:δ_x], sims[ii][:γ], γδc)

                        lNX = lNXmat(X_pred, sims[ii][:μ_x],
                                [ βγ_x[k] for k = 1:(K-1) ],
                                δγ_x) # npred by H
                    end

                elseif γlocal

                    if γδc == Inf  # subset method

                        lNX = zeros(Float64, npred, H)

                        for h = 1:H

                            γ_h = sims[ii][:γ][h,:]
                            γindx = findall(γ_h)
                            nγ = length( γindx )

                            if nγ == 0
                                # do nothing because column h is already zeros
                            elseif nγ == 1

                                lNX[:,h] = logpdf.(Normal(sims[ii][:μ_x][h,γindx[1]], sqrt(sims[ii][:δ_x][h,γindx[1]])), vec(X_pred[:,γindx]))

                            elseif nγ > 1 && !isnothing(sims[ii][:β_x]) # indicates Σx_type is :full

                                βγ_x_h, δγ_x_h = βδ_x_h_modify_γ( [ sims[ii][:β_x][j][h,:] for j = 1:length(sims[ii][:β_x]) ],
                                    sims[ii][:δ_x][h,:], γ_h, γδc) # either variance-inflated or subset

                                lNX[:,h] = lNX_sqfChol(Matrix(X_pred[:,γindx]'), sims[ii][:μ_x][h,γindx], βγ_x_h, δγ_x_h, true)

                            elseif nγ > 1 && isnothing(sims[ii][:β_x]) # indicates Σx_type is :diag

                                lNX[:,h] = lNX_Σdiag( Matrix(X_pred[:,γindx]'), sims[ii][:μ_x][h,γindx], sims[ii][:δ_x][h,γindx] )

                            end

                        end

                    elseif isnothing(γδc) # integration method

                        lNX = zeros(Float64, npred, H)

                        for h = 1:H

                            γ_h = sims[ii][:γ][h,:]
                            γindx = findall(γ_h)
                            nγ = length( γindx )

                            if nγ == 0
                               # do nothing because column h is already zeros
                            elseif nγ == 1 && !isnothing(sims[ii][:β_x]) # indicates Σx_type is :full

                                σ2x = sqfChol_to_Σ( [ sims[ii][:β_x][j][h,:] for j = 1:length(sims[ii][:β_x]) ], sims[ii][:δ_x][h,:] ).mat[γindx, γindx][1]
                                lNX[:,h] = logpdf.(Normal(sims[ii][:μ_x][h,γindx[1]], sqrt(σ2x)), vec(X_pred[:,γindx[1]]))

                            elseif nγ == 1 && isnothing(sims[ii][:β_x]) # indicates Σx_type is :diag

                                σ2x = deepcopy( sims[ii][:δ_x][h,γindx[1]] )
                                lNX[:,h] = logpdf.(Normal(sims[ii][:μ_x][h,γindx[1]], sqrt(σ2x)), vec(X_pred[:,γindx[1]]))

                            elseif nγ > 1 && !isnothing(sims[ii][:β_x]) # indicates Σx_type is :full

                                Σx = PDMat( sqfChol_to_Σ( [ sims[ii][:β_x][j][h,:] for j = 1:length(sims[ii][:β_x]) ], sims[ii][:δ_x][h,:] ).mat[γindx, γindx] )
                                lNX[:,h] = logpdf( MultivariateNormal(sims[ii][:μ_x][h,γindx], Σx), Matrix(X_pred[:,γindx]') )

                            elseif nγ > 1 && isnothing(sims[ii][:β_x]) # indicates Σx_type is :diag

                                lNX[:,h] = lNX_Σdiag( Matrix(X_pred[:,γindx]'), sims[ii][:μ_x][h,γindx], sims[ii][:δ_x][h,γindx] )

                            end

                        end

                    else  # variance-inflation method
                        throw("Variance-inflation variable selection not implemented for local variable selection.")
                    end

                end

            else
                β_x = deepcopy(sims[ii][:β_x])
                δ_x = deepcopy(sims[ii][:δ_x])

                if isnothing(β_x) # indicates Σx_type is :diag
                    lNX = lNXmat_Σdiag(X_pred, sims[ii][:μ_x], δ_x) # npred by H
                else # indicates Σx_type is :full
                    lNX = lNXmat(X_pred, sims[ii][:μ_x],
                        [ β_x[k] for k = 1:(K-1) ],
                        δ_x) # npred by H
                end
            end

            lωNX_mat = broadcast(+, permutedims(sims[ii][:lω]), lNX) # npred by H
            lωNX_vec = vec( BayesInference.logsumexp(lωNX_mat, 2) ) # npred vec

            out[ii,:,:] = broadcast(-, lωNX_mat, lωNX_vec) # normalized, npred by H
        end
    else  # K = 1
        for ii = 1:nsim

            if useγ

                if γglobal

                    if γδc == Inf || isnothing(γδc) # subset or integration method

                        if sims[ii][:γ][1]
                            lNX = lNXmat(vec(X_pred), vec(sims[ii][:μ_x]), vec(sims[ii][:δ_x])) # npred by H matrix
                        else
                            lNX = zeros(Float64, npred, H)
                        end

                    else  # variance-inflation method

                        δγ_x = δ_x_modify_γ(sims[ii][:δ_x], sims[ii][:γ], γδc)
                        lNX = lNXmat(vec(X_pred), vec(sims[ii][:μ_x]), vec(δγ_x))

                    end

                elseif γlocal

                    if γδc == Inf || isnothing(γδc) # subset or integration method

                        lNX = zeros(Float64, npred, H)

                        for h = 1:H

                            if sims[ii][:γ][h,1]
                                σ2x = deepcopy(sims[ii][:δ_x][h,1])
                                lNX[:,h] = logpdf.( Normal(sims[ii][:μ_x][h,1], sqrt(σ2x)), vec(X_pred) )
                            else
                                # do nothing, the column is already zeros
                            end

                        end

                    else  # variance-inflation method
                        throw("Variance-inflation variable selection not implemented for local variable selection.")
                    end

                end

            else
                lNX = lNXmat(vec(X_pred), vec(sims[ii][:μ_x]), vec(sims[ii][:δ_x])) # npred by H
            end

            lωNX_mat = broadcast(+, permutedims(sims[ii][:lω]), lNX) # npred by H
            lωNX_vec = vec( BayesInference.logsumexp(lωNX_mat, 2) ) # npred vec

            out[ii,:,:] = broadcast(-, lωNX_mat, lωNX_vec) # normalized, npred by H
        end
    end

    return out
end

# npred = 100
# X_pred = hcat(collect(range(Float32(-2.0), length=npred, stop=Float32(2.0))))
#
# ldw = ldensweight_mat(X_pred, sims)
#
# dw = exp.(ldw)
# sum(dw[1,1,:])
#
# mean_dw = mean(dw, dims=1)
# h = 6
# plot([scatter(x=X_pred[:,1], y=mean_dw[1,:,h], mode="scatter")])

function getEy(X_pred::Array{T,2}, dw_mat::Array{T,3}, sims::Union{Array{Dict{Symbol,Any},1}, Array{Any,1}}) where T <: Real
    nsim, npred, H = size(dw_mat)
    nsim2 = length(sims)
    H2, K = size(sims[1][:β_y])
    npred3, K3 = size(X_pred)

    nsim == nsim2 || throw(error("Dimension mismatch."))
    H == H2 || throw(error("Dimension mismatch."))
    npred == npred3 || throw(error("Dimension mismatch."))
    K == K3 || throw(error("Dimension mismatch."))

    useγ = haskey(sims[1], :γ)
    if useγ
        γglobal = ( typeof(sims[1][:γ]) == BitArray{1} )
        γlocal = ( typeof(sims[1][:γ]) == BitArray{2} )
    end

    out = zeros(T, nsim, npred)

    for ii = 1:nsim
        if useγ
            if γglobal
                βγ_y = β_y_modify_γ(sims[ii][:β_y], sims[ii][:γ])
            elseif γlocal
                βγ_y = sims[ii][:β_y] .* sims[ii][:γ]
            end
        else
            βγ_y = deepcopy(sims[ii][:β_y])
        end
        for j = 1:npred
            for h = 1:H
                Ey_h = sims[ii][:μ_y][h] - sum( βγ_y[h,:] .* (X_pred[j,:] - sims[ii][:μ_x][h,:]) )
                out[ii,j] += dw_mat[ii, j, h] * Ey_h
            end
        end
    end

    out
end

function qnorm_mix(q::Float64, μ::Vector{Float64}, σ::Vector{Float64}, w::Vector{Float64})
    all(w .>= 0.0) || throw(error("Weights must be non-negative"))
    q > 0.0 && q < 1.0 || throw(error("Quantile must be between 0 and 1."))

    sum(w) ≈ 1.0 ? nothing : w = w ./ sum(w)

    d = Distributions.Normal()

    if q >= 0.5
        z_max = Distributions.quantile(d, q)
        z_min = -z_max

        x_max = maximum( σ .* z_max .+ μ )
        x_min = minimum( σ .* z_min .+ μ )
    else
        z_min = Distributions.quantile(d, q)
        z_max = -z_min

        x_max = maximum( σ .* z_max .+ μ )
        x_min = minimum( σ .* z_min .+ μ )
    end

    if abs(q - w'cdf.(d, ( x_max .- μ ) ./ σ)) < 1.0e-4
        out = x_max
    elseif abs(q - w'cdf.(d, ( x_min .- μ ) ./ σ)) < 1.0e-4
        out = x_min
    elseif x_max > x_min
        f(x) = ( q - w'cdf.(d, ( x .- μ ) ./ σ) )
        out = Roots.find_zero(f, (x_min, x_max))
    end

    return out
end

function getQuant(q::Float64, X_pred::Array{T,2}, dw_mat::Array{T,3}, sims::Union{Array{Dict{Symbol,Any},1}, Array{Any,1}}) where T <: Real
    nsim, npred, H = size(dw_mat)
    nsim2 = length(sims)
    H2, K = size(sims[1][:β_y])
    npred3, K3 = size(X_pred)

    q > 0.0 && q < 1.0 || throw(error("Quantile must be between 0 and 1."))
    nsim == nsim2 || throw(error("Dimension mismatch."))
    H == H2 || throw(error("Dimension mismatch."))
    npred == npred3 || throw(error("Dimension mismatch."))
    K == K3 || throw(error("Dimension mismatch."))

    useγ = haskey(sims[1], :γ)
    if useγ
        γglobal = ( typeof(sims[1][:γ]) == BitArray{1} )
        γlocal = ( typeof(sims[1][:γ]) == BitArray{2} )
    end

    out = zeros(T, nsim, npred)

    for ii = 1:nsim
        if useγ
            if γglobal
                βγ_y = β_y_modify_γ(sims[ii][:β_y], sims[ii][:γ])
            elseif γlocal
                βγ_y = sims[ii][:β_y] .* sims[ii][:γ]
            end
        else
            βγ_y = deepcopy(sims[ii][:β_y])
        end

        σ = sqrt.(sims[ii][:δ_y])

        for j = 1:npred

            μ = [ sims[ii][:μ_y][h] - sum( βγ_y[h,:] .* (X_pred[j,:] - sims[ii][:μ_x][h,:]) ) for h = 1:H ]
            w = deepcopy(dw_mat[ii, j, :])

            try
                out[ii,j] = qnorm_mix(q, μ, σ, w)
            catch
                println("No root found, simid $(ii) Xpred indx $(j)")
                println("μ=$(μ)")
                println("σ=$(σ)")
                println("w=$(w)")
                qnorm_mix(q, μ, σ, w) # for debugging
            end
        end
    end

    out
end

# Ey = getEy(X_pred, dw, sims)

# ii = 2
# plot([scatter(x=X_pred[:,1], y=Ey[ii,:], mode="scatter")])

# mean_Ey = mean(Ey, dims=1)
# p10_Ey = [quantile(Ey[:,j], 0.1) for j = 1:size(Ey,2)]
# p50_Ey = [quantile(Ey[:,j], 0.5) for j = 1:size(Ey,2)]
# p90_Ey = [quantile(Ey[:,j], 0.9) for j = 1:size(Ey,2)]
#
# trace1 = [scatter(x=X_pred[:,1], y=mean_Ey[1,:], mode="scatter")]
# trace2 = [scatter(x=X[:,1], y=y, mode="markers")]
# data = [trace1, trace2]
# Plotly.plot(trace2)
#
# trace0 = [scatter(x=X_pred[:,1], y=p10_Ey, mode="scatter")]
# Plotly.plot(trace0)
#
# npred = 5
# X_pred = hcat(collect(range(Float32(-2.0), length=npred, stop=Float32(2.0))))
# ngrid = 100
# y_grid = collect(range(Float32(-2.0), length=ngrid, stop=Float32(2.0)))

function getlogdens_EY(X_pred::Array{T,2}, y_grid::Array{T,1},
    ldw_mat::Array{T,3}, sims::Union{Array{Dict{Symbol,Any},1}, Array{Any,1}}) where T <: Real

    ngrid = length(y_grid)
    nsim, npred, H = size(ldw_mat)
    nsim2 = length(sims)
    H2, K = size(sims[1][:β_y])
    npred3, K3 = size(X_pred)

    nsim == nsim2 || throw(error("Dimension mismatch."))
    H == H2 || throw(error("Dimension mismatch."))
    npred == npred3 || throw(error("Dimension mismatch."))
    K == K3 || throw(error("Dimension mismatch."))

    useγ = haskey(sims[1], :γ)
    if useγ
        γglobal = ( typeof(sims[1][:γ]) == BitArray{1} )
        γlocal = ( typeof(sims[1][:γ]) == BitArray{2} )
    end

    ldens0 = Array{T,4}(undef, (nsim, npred, ngrid, H))
    Ey = zeros(T, nsim, npred)

    for ii = 1:nsim
        if useγ
            if γglobal
                βγ_y = β_y_modify_γ(sims[ii][:β_y], sims[ii][:γ])
            elseif γlocal
                βγ_y = sims[ii][:β_y] .* sims[ii][:γ]
            end
        else
            βγ_y = deepcopy(sims[ii][:β_y])
        end
        for j = 1:npred
            for h = 1:H
                Ey_h = sims[ii][:μ_y][h] - sum( βγ_y[h,:] .* (X_pred[j,:] - sims[ii][:μ_x][h,:]) )
                Ey[ii,j] += exp(ldw_mat[ii, j, h]) * Ey_h
                for jj = 1:ngrid
                    ldens0[ii,j,jj,h] = ldw_mat[ii, j, h] + logpdf(Normal(Ey_h, sqrt(sims[ii][:δ_y][h])), y_grid[jj])
                end
            end
        end
    end

    ldens = reshape(BayesInference.logsumexp(ldens0,4), (nsim, npred, ngrid))

    return ldens, Ey
end

# ldw = ldensweight_mat(X_pred, sims)
# ldens, Ey = getlogdens_EY(X_pred, y_grid, ldw, sims)
# dens = exp.(ldens)
# mean_dens = mean(dens, dims=1)
#
# xind = 4
# println(X_pred[xind,:])
#
# trace0 = [scatter(x=y_grid, y=mean_dens[1,xind,:], mode="scatter")]
# Plotly.plot(trace0)
#
# plot(dens[:,xind,25])
# plot(dens[:,xind,75])