text:
export 	images_10p,
	calculate_expected_triplets,
	calculate_required_images,
	calculate_images_ppi,
	calculate_incident_photons,
	minutues_to_measure,
	cor_ISC,
	cor_FSC,
	angle_between,
	angle_between_simple,
	euler,
	get_euler_angles,
	sphericalToCartesian,
	cartesianToSpherical,
	gaussian,
	gaussian_distribution,
	normal_distribution,
	replace_NaN!,
	random_rotation,
	random_rotation_step,
	SON_parametrized,
	SON_parametrized_derivative,
	random_doublet,
	random_triplet,
	serializeToFile,
	deserializeFromFile,
    gaussian_filter,
	generate_gaussian_kernel,
	calculate_ppi

images_10p = Int64[1280000,5120000,20480000,81920000,327680000,3276800000]
calculate_expected_triplets = (images,ppi::Int64) -> images * (ppi * (ppi-1) * (ppi-2) / 6.0)
calculate_required_images = (triplets,ppi::Int64) -> triplets / (ppi * (ppi-1) * (ppi-2) /6.0)
calculate_images_ppi = (ppi::Int64) -> map((img)->ceil(Int64,calculate_required_images(calculate_expected_triplets(img, 10), ppi)/8)*8, images_10p)
calculate_incident_photons = (ppi::Int64, incident_10p::Int64=550) -> round(Int64, (ppi / 10.0) * incident_10p)
calculate_ppi = (incident_photons::Int64, incident_10p::Int64=550) -> round(Int64, (incident_photons * 10.0) / incident_10p)

minutues_to_measure = (images, pulses_per_second::Int64=27000, hitrate::Float64=1.0) -> images / (pulses_per_second*60) / hitrate

###############################################################
#                     MATH
###############################################################

"""Pearson correlation between two vectors"""
function cor_ISC(a::Vector{Complex{Float64}}, b::Vector{Complex{Float64}})
    return cor(real(a), real(b))
end

"""Fourier shell correlation between two complex vectors"""
function cor_FSC(a::Vector{Complex{Float64}}, b::Vector{Complex{Float64}})
    na = norm(a)
    nb = norm(b)
    if na > eps() && nb > eps()
        return abs(sum(a.*conj(b))/(na*nb))
    else
        return 0.0
    end
end

"Angle between two vectors ranging from 0 to 2.0*pi"
function angle_between(p1::Vector{Float64}, p2::Vector{Float64})
    val = clamp(Base.dot(p1,p2)/ (norm(p1) * norm(p2)), -1.0, 1.0)
    if abs(val) > 0.0
        angle = acos(val)
        if cross(p1, p2)[3] > 0
            angle =  2.0*pi - angle
        end
        return angle
    else
        return 0.0
    end
end

"Angle between two vectors ranging from 0 to pi"
function angle_between_simple(p1::Vector{Float64}, p2::Vector{Float64})
    angle = acos( clamp(Base.dot(p1,p2)/ (norm(p1) * norm(p2)), -1.0, 1.0))
    return angle
end

"""
Calculate a three dimensional rotation matrix from the euler angles.

@param a: alpha, angle between the x-axis and the line of nodes
@param b: beta, angle between the z axis of the different coordinate systems
@param c: gamma, angle between the line of nodes and the X-axis
"""
function euler(a::Float64, b::Float64, c::Float64)

    ca, cb, cc = cos(a), cos(b), cos(c)
    sa, sb, sc = sin(a), sin(b), sin(c)

    @fastmath return Float64[[ cc * cb * ca - sc * sa, cc * cb * sa + sc * ca, -cc * sb] [-sc * cb * ca - cc * sa, -sc * cb * sa + cc * ca, sc * sb] [ sb * ca, sb * sa, cb ]]
end

"Calculates the Euler angles of a 3D rotation matrix"
function get_euler_angles(m)
    a = mod(atan2(m[3, 2], m[3, 1]), 2.0*pi)
    b = mod(atan2((m[3, 1] + m[3, 2]) / (cos(a) + sin(a)), m[3, 3]), 2.0*pi)
    c = mod(atan2(m[2, 3], -m[1, 3]), 2.0*pi)

    return a, b, c
end

"""Converts spherical coordinates to cartesian coordinates"""
function sphericalToCartesian(r::Float64, phi::Float64, theta::Float64)
    x = r*cos(phi)*sin(theta)
    y = r*sin(phi)*sin(theta)
    z = r*cos(theta)
    return Float64[x,y,z]
end

"""Converts a vector to spherical coordinates
Returns r, phi,theta
"""
function cartesianToSpherical(vec::Vector{Float64})
    r::Float64 = norm(vec)
    phi::Float64 = atan2(vec[2],vec[1])
    if phi < 0.0 phi = phi + 2.0*pi end
    theta::Float64 = r > 0.0 ? acos(clamp(vec[3]/r, -1.0, 1.0)) : 0.0
    return r,phi,theta
end

"Creates a gaussian distributed number with sigma and mu"
function gaussian(mu::Float64, sigma::Float64)
    return randn()*sigma + mu
end

"""Calcualtes the gaussian distribution wiht sigma and mu"""
function gaussian_distribution(x::Float64, mu::Float64, sigma::Float64)
    @assert sigma != 0.0
    s2 = 2.0 * sigma^2
    1/sqrt(pi*s2)*exp(-norm(x-mu)^2/s2)
end

"""Calcualtes the gaussian distribution wiht sigma and mu"""
function gaussian_distribution(x::Vector{Float64}, mu::Vector{Float64}, sigma::Vector{Float64})
    @assert sigma != 0.0
    return pdf(Distributions.MvNormal(mu, sigma), x)
end

"""Calcualtes the normal distribution wiht sigma=1.0 and mu=0.0"""
function normal_distribution(x::Vector)
    gaussian_distribution(x, 0.0, 1.0)
end

"""In a numerical Vector, replaces NaN with zero"""
function replace_nan!{T}(x::Array{T})
    for i = eachindex(x)
        if isnan(x[i])
            x[i] = zero(T)
        end
    end
end

"""In a numerical Vector, replaces NaN with zero"""
function replace_inf!{T}(x::Array{T})
    for i = eachindex(x)
        if isinf(x[i])
            x[i] = zero(T)
        end
    end
end
"""Calculates a random rotation in `dim`
Source: http://home.lu.lv/~sd20008/papers/essays/Random%20unitary%20[paper].pdf"""
function random_rotation(dim::Int64)
    Q,R = qr(randn(dim, dim))
    r = diag(R)
    r = r ./ abs.(r)
    L = diagm(r)
    res = Q*L
    if(det(res) > 0)
        return res
    else
        random_rotation(dim)
    end
end

"""returns a random rotation with a slight angle beta=pi/8
Taken from "Random rotations: characters and random walk on SO(N)"""
function random_rotation_step(b, beta=pi/8)
    R = eye(b)
    R[1,1] = cos(beta)
    R[1,2] = -sin(beta)
    R[2,1] = sin(beta)
    R[2,2] = cos(beta)

    C = random_rotation(b)
    B = inv(C)*R*C
    return B
end

"Calculates a rotation matrix in `dim` dimensions with dim*(dim-1)/2 rotation angles"
function SON_parametrized(dim, angles)
    R = eye(dim)
    assert(dim*(dim-1)/2 == length(angles))
    num = 1
    for i = 1:dim
        for j = i+1:dim
            M = eye(dim)
            c = cos(angles[num])
            s = sin(angles[num])
            M[i,i] = c
            M[j,j] = c
            M[i,j] = s
            M[j,i] = -s
            R *= M
            num += 1
        end
    end
    return R
end

"Calculates the derivative of R by the kth parameter"
function SON_parametrized_derivative(dim, angles, k)
    R = eye(dim)
    num = 1
    for i = 1:dim
        for j = i+1:dim
            c = cos(angles[num])
            s = sin(angles[num])
            if(num != k)
                M = eye(dim)
                M[i,i] = c
                M[j,j] = c
                M[i,j] = s
                M[j,i] = -s
                R *= M
            else
                M = zeros(dim,dim)
                M[i,i] = -s
                M[j,j] = -s
                M[i,j] = c
                M[j,i] = -c
                R *= M
            end
            num += 1
        end
    end
    return R
end

"""Creates random doublet"""
function random_doublet(K::Int64=8, half::Bool=true)
    k1 = rand(1:K)
    k2 = rand(1:K)
    if half
        if k1 >= k2
            return (k1, k2)
        else
            return random_doublet(K, half)
        end
    else
        return (k1,k2)
    end
end

"""Creates random triplet"""
function random_triplet(K=8, half::Bool=true)
    k1 = rand(1:K)
    k2 = rand(1:K)
    k3 = rand(1:K)
    if half
        if k1 >= k2 && k2 >= k3
            return (k1,k2,k3)
        else return random_triplet(K, half)
        end
    else
        return (k1,k2,k3)
    end
end

###############################################################
#                     FILES
###############################################################
using FileIO

"""Serializes the data to a file"""
function serializeToFile(filename::String, data)
	if contains(filename, ".jld")
		save("$(filename)", "data", data)
	else
		open(filename, "w+") do out
			serialize(out, data)
		end
	end
	println("Wrote to $filename.")
end

"""Deserializes a julia data object from a file using JLD.jl and fixing backward compatibility"""
function deserializeFromFile(filename::String)
	if contains(filename, ".jld")
		data = load("$(filename)")
		println("Read from $filename.")
		return data["data"]
	else
		open(filename, "r") do out
	    	data = deserialize(out)
	    end
		println("Read from $filename.")
	    return data
	end

	# try
	# 	return deserializeFromFile_raw(filename)
	# catch
		# try
	    # 	data = load("$(filename).jld")
		# 	return data["data"]
		# catch err
		# 	if isa(err, FileIO.UnknownFormat{FileIO.File{FileIO.DataFormat{:UNKNOWN}}}) || isa(err, FileIO.File{FileIO.DataFormat{:UNKNOWN}})
		# 		println("Trying to rewrite $filename into JLD format and reload")
		# 		run(`julia5 -e "using ThreePhotons5; using JLD; data = deserializeFromFile(\"$(filename)\"); save(\"$(filename).jld\", \"data\", data)"`)
		# 		try
		# 			data = load("$(filename).jld")
		# 			return data["data"]
		# 		catch err
		# 			println("Something failed loading $filename: $err")
		# 		end
		# 	else
		# 		try
		# 			data = load("$(filename).jld")
		# 			return data["data"]
		# 		catch err
		# 			println("Something failed loading $filename: $err")
		# 		end
		# 	end
		# end
	# end
end

function deserializeFromFile_raw(filename::String)
    out = open(filename, "r")
    data = deserialize(out)
    close(out)
    return data
end

###############################################################
#                     MISC
###############################################################

"""Say something via MacOS say function"""
function say(str::String)
    voices = ["Agnes",
    "Albert",
    "Alex",
    "Alice",
    "Alva",
    "Amelie",
    "Anna",
    "Bad",
    "Bahh",
    "Bells",
    "Boing",
    "Bruce",
    "Bubbles",
    "Carmit",
    "Cellos",
    "Damayanti",
    "Daniel",
    "Deranged",
    "Diego",
    "Ellen",
    "Fiona",
    "Fred",
    "Good",
    "Hysterical",
    "Ioana",
    "Joana",
    "Junior",
    "Kanya",
    "Karen",
    "Kathy",
    "Kyoko",
    "Laura",
    "Lekha",
    "Luciana",
    "Maged",
    "Mariska",
    "Mei-Jia",
    "Melina",
    "Milena",
    "Moira",
    "Monica",
    "Nora",
    "Paulina",
    "Pipe",
    "Princess",
    "Ralph",
    "Samantha",
    "Sara",
    "Satu",
    "Sin-ji",
    "Tessa",
    "Thomas",
    "Ting-Ting",
    "Trinoids",
    "Veena",
    "Vicki",
    "Victoria",
    "Whisper",
    "Xander",
    "Yelda",
    "Yuna",
    "Zarvox",
    "Zosia",
    "Zuzana"]
    try
        voice = voices[rand(1:length(voices))]
        # println(voice)
        run(`say $str --voice $voice`)
    end
end

function whisper(str)
    try
        run(`say $str --voice 'whisper'`)
    end
end

ϕ(x::Float64, σ::Float64) = exp(-x^2 / (2 * σ^2)) / (σ * sqrt(2 * π))
ϕ(x::Vector{Float64}, σ::Float64) = exp(-norm(x)^2 / (2 * σ^2)) / (σ * sqrt(2 * π))

function generate_gaussian_kernel(σ::Float64, kernel_size::Vector{Int64})
    center = map((dim) -> mod(dim,2) == 0 ? (dim/2)+1 : ceil(dim/2), kernel_size)
	if length(kernel_size) == 1
    	return Float64[ϕ(x-center[1],σ) for x=1:kernel_size[1]]
	elseif length(kernel_size) == 2
		return Float64[ϕ(Float64[x,y]-center,σ) for x=1:kernel_size[1], y=1:kernel_size[2]]
	end
end

function generate_gaussian_kernel(σ::Float64, kernel_size::Tuple)
	generate_gaussian_kernel(σ, Int64[t for t in kernel_size])
end

function generate_gaussian_kernel(σ::Float64, kernel_size::Int64)
	return generate_gaussian_kernel(σ, [kernel_size])
end

function gaussian_filter(data::Vector{Float64}, σ::Float64=1.0)
    kernel = generate_gaussian_kernel(σ, Base.size(data))
    fkern = fft(ifftshift(kernel))
    fdata = fft((data))
    return real((ifft(fkern.*fdata)))
end

function gaussian_filter(data::Matrix{Float64}, σ::Float64=1.0)
    kernel = generate_gaussian_kernel(σ, Base.size(data))
    fkern = fft(ifftshift(kernel))
    fdata = fft((data))
    return real((ifft(fkern.*fdata)))
end


text:
import MathOptInterface
const MOI = MathOptInterface
const MOIU = MOI.Utilities

import MathProgBase
const MPB = MathProgBase.SolverInterface

MOIU.@model(
    Model,
    (MOI.ZeroOne, MOI.Integer),
    (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan, MOI.Interval),
    (),
    (),
    (),
    (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction),
    (),
    (),
    true    # is_optimizer
)

const MOI_SCALAR_SETS = (
    MOI.EqualTo{Float64}, MOI.GreaterThan{Float64}, MOI.LessThan{Float64},
    MOI.Interval{Float64}
)

const MOI_SCALAR_FUNCTIONS = (
    MOI.ScalarAffineFunction{Float64}, MOI.ScalarQuadraticFunction{Float64}
)

"Struct to contain the MPB solution."
struct MPBSolution
    status::Symbol
    raw_status_string::String
    objective_value::Float64
    primal_solution::Dict{MOI.VariableIndex, Float64}
end

"""
    Optimizer(
        solver_command::String,
        options::Vector{String} = String[];
        filename::String = ""
    )

# Example

    Optimizer(Ipopt.amplexe, ["print_level=0"])
"""
function Optimizer(
        solver_command::String,
        options::Vector{String} = String[];
        filename::String = ""
)
    model = Model{Float64}()
    model.ext[:MPBSolver] = AmplNLSolver(solver_command, options, filename = filename)
    model.ext[:VariablePrimalStart] = Dict{MOI.VariableIndex, Float64}()
    return model
end

function MOI.empty!(model::Model{Float64})
    model.name = ""
    model.senseset = false
    model.sense = MOI.FEASIBILITY_SENSE
    model.objectiveset = false
    model.objective = MOI.ScalarAffineFunction{Float64}(MOI.ScalarAffineTerm{Float64}[], 0.0)
    model.num_variables_created = 0
    model.variable_indices = nothing
    empty!(model.single_variable_mask)
    empty!(model.lower_bound)
    empty!(model.upper_bound)
    empty!(model.var_to_name)
    model.name_to_var = nothing
    model.nextconstraintid = 0
    empty!(model.con_to_name)
    model.name_to_con = nothing
    empty!(model.constrmap)
    MOI.empty!(model.moi_scalaraffinefunction)
    MOI.empty!(model.moi_scalarquadraticfunction)
    solver = model.ext[:MPBSolver]
    empty!(model.ext)
    model.ext[:MPBSolver] = solver
    model.ext[:VariablePrimalStart] = Dict{MOI.VariableIndex, Float64}()
    return
end

Base.show(io::IO, ::Model) = print(io, "An AmplNLWriter model")

MOI.get(::Model, ::MOI.SolverName) = "AmplNLWriter"

set_to_bounds(set::MOI.LessThan) = (-Inf, set.upper)
set_to_bounds(set::MOI.GreaterThan) = (set.lower, Inf)
set_to_bounds(set::MOI.EqualTo) = (set.value, set.value)
set_to_bounds(set::MOI.Interval) = (set.lower, set.upper)
set_to_cat(set::MOI.ZeroOne) = :Bin
set_to_cat(set::MOI.Integer) = :Int

struct NLPEvaluator{T} <: MPB.AbstractNLPEvaluator
    inner::T
    variable_map::Dict{MOI.VariableIndex, Int}
    num_inner_con::Int
    has_nlobjective::Bool
    objective_expr::Union{Nothing, Expr}
    scalar_constraint_expr::Vector{Expr}
end

"""
MathProgBase expects expressions with variables denoted by `x[i]` for contiguous
`i`. However, JuMP 0.19 creates expressions with `x[MOI.VariableIndex(i)]`. So
we have to recursively walk the expression replacing instances of
MOI.VariableIndex by a corresponding integer.
"""
function replace_variableindex_by_int(variable_map, expr::Expr)
    for (i, arg) in enumerate(expr.args)
        expr.args[i] = replace_variableindex_by_int(variable_map, arg)
    end
    return expr
end
function replace_variableindex_by_int(variable_map, expr::MOI.VariableIndex)
    return variable_map[expr]
end
replace_variableindex_by_int(variable_map, expr) = expr

# In the next section we match up the MPB nonlinear functions to the MOI API.
# This is basically just a rename.
function MPB.initialize(d::NLPEvaluator, requested_features::Vector{Symbol})
    if d.inner !== nothing
        MOI.initialize(d.inner, requested_features)
    end
    return
end

function MPB.features_available(d::NLPEvaluator)
    if d.inner !== nothing
        return MOI.features_available(d.inner)
    else
        return [:ExprGraph]
    end
end

function MPB.obj_expr(d::NLPEvaluator)
    if d.has_nlobjective
        # If d.objective_expr === nothing, then the objective must be nonlinear.
        # Query it from the inner NLP evaluator.
        expr = MOI.objective_expr(d.inner)
        return replace_variableindex_by_int(d.variable_map, expr)
    elseif d.objective_expr !== nothing
        # d.objective_expr is a SingleVariable, ScalarAffineFunction, or a
        # ScalarQuadraticFunction from MOI. If it is unset, it will be `nothing`
        # (we enforce this when creating the NLPEvaluator in `optimize!`).
        return d.objective_expr
    else
        return :(0.0)
    end
end

function MPB.constr_expr(d::NLPEvaluator, i)
    # There are two types of constraints in the model:
    # i = 1, 2, ..., d.num_inner_con are the nonlinear constraints. We access
    # these from the inner NLP evaluator.
    # i = d.num_inner_con + 1, d.num_inner_con + 2, ..., N are the linear or
    # quadratic constraints added by MOI.
    if i <= d.num_inner_con
        expr = MOI.constraint_expr(d.inner, i)
        return replace_variableindex_by_int(d.variable_map, expr)
    else
        return d.scalar_constraint_expr[i - d.num_inner_con]
    end
end

# In the next section, we need to convert MOI functions and sets into expression
# graphs. First, we convert functions (SingleVariable, ScalarAffine, and
# ScalarQuadratic) into expression graphs.
function func_to_expr_graph(func::MOI.SingleVariable, variable_map)
    return Expr(:ref, :x, variable_map[func.variable])
end

function func_to_expr_graph(func::MOI.ScalarAffineFunction, variable_map)
    expr = Expr(:call, :+, func.constant)
    for term in func.terms
        push!(expr.args, Expr(:call, :*, term.coefficient,
            Expr(:ref, :x, variable_map[term.variable_index])
        ))
    end
    return expr
end

function func_to_expr_graph(func::MOI.ScalarQuadraticFunction, variable_map)
    expr = Expr(:call, :+, func.constant)
    for term in func.affine_terms
        push!(expr.args, Expr(:call, :*, term.coefficient,
            Expr(:ref, :x, variable_map[term.variable_index])
        ))
    end
    for term in func.quadratic_terms
        index_1 = variable_map[term.variable_index_1]
        index_2 = variable_map[term.variable_index_2]
        coef = term.coefficient
        # MOI defines quadratic as 1/2 x' Q x :(
        if index_1 == index_2
            coef *= 0.5
        end
        push!(expr.args, Expr(:call, :*, coef,
            Expr(:ref, :x, index_1),
            Expr(:ref, :x, index_2)
        ))
    end
    return expr
end

MOI.supports(::Model, ::MOI.NLPBlock) = true

function MOI.set(model::Model, ::MOI.NLPBlock, block)
    model.ext[:NLPBlock] = block
    return
end

MOI.get(model::Model, ::MOI.NLPBlock) = get(model.ext, :NLPBlock, nothing)

function MOI.supports(
    ::Model, ::MOI.VariablePrimalStart, ::Type{MOI.VariableIndex}
)
    return true
end

function MOI.set(
    model::Model,
    ::MOI.VariablePrimalStart,
    variable::MOI.VariableIndex,
    value::Union{Nothing, Real},
)
    if value === nothing
        delete!(model.ext[:VariablePrimalStart], variable)
    else
        model.ext[:VariablePrimalStart][variable] = Float64(value)
    end
    return
end

function MOI.get(
    model::Model, ::MOI.VariablePrimalStart, variable::MOI.VariableIndex
)
    return get(model.ext[:VariablePrimalStart], variable, nothing)
end

# Next, we need to combine a function `Expr` with a MOI set into a comparison.

function funcset_to_expr_graph(func::Expr, set::MOI.LessThan)
    return Expr(:call, :<=, func, set.upper)
end

function funcset_to_expr_graph(func::Expr, set::MOI.GreaterThan)
    return Expr(:call, :>=, func, set.lower)
end

function funcset_to_expr_graph(func::Expr, set::MOI.EqualTo)
    return Expr(:call, :(==), func, set.value)
end

function funcset_to_expr_graph(func::Expr, set::MOI.Interval)
    return Expr(:comparison, set.lower, :<=, func, :<=, set.upper)
end

# A helper function that converts a MOI function and MOI set into a comparison
# expression.

function moi_to_expr_graph(func, set, variable_map)
    func_expr = func_to_expr_graph(func, variable_map)
    return funcset_to_expr_graph(func_expr, set)
end

# Now we're ready to optimize model.
function MOI.optimize!(model::Model)
    # Extract the MathProgBase solver from the model. Recall it's a MOI
    # attribute which we're careful not to delete on calls to `empty!`.
    mpb_solver = model.ext[:MPBSolver]

    # Get the optimzation sense.
    opt_sense = MOI.get(model, MOI.ObjectiveSense())
    sense = opt_sense == MOI.MAX_SENSE ? :Max : :Min

    # Get the NLPBlock from the model.
    nlp_block = MOI.get(model, MOI.NLPBlock())

    # Extract the nonlinear constraint bounds. We're going to append to these
    # g_l and g_u vectors later.
    num_con = 0
    moi_nlp_evaluator = nlp_block !== nothing ? nlp_block.evaluator : nothing
    if nlp_block !== nothing
        num_con += length(nlp_block.constraint_bounds)
    end
    g_l = fill(-Inf, num_con)
    g_u = fill(Inf, num_con)
    if nlp_block !== nothing
        for (i, bound) in enumerate(nlp_block.constraint_bounds)
            g_l[i] = bound.lower
            g_u[i] = bound.upper
        end
    end

    # Intialize the variables. We need to form a mapping between the MOI
    # VariableIndex and an Int in order to replace instances of
    # `x[VariableIndex]` with `x[i]` in the expression graphs.
    variables = MOI.get(model, MOI.ListOfVariableIndices())
    num_var = length(variables)
    variable_map = Dict{MOI.VariableIndex, Int}()
    for (i, variable) in enumerate(variables)
        variable_map[variable] = i
    end

    # Extract variable bounds.
    x_l = fill(-Inf, num_var)
    x_u = fill(Inf, num_var)
    for set_type in MOI_SCALAR_SETS
        for c_ref in MOI.get(model,
            MOI.ListOfConstraintIndices{MOI.SingleVariable, set_type}())
            c_func = MOI.get(model, MOI.ConstraintFunction(), c_ref)
            c_set = MOI.get(model, MOI.ConstraintSet(), c_ref)
            v_index = variable_map[c_func.variable]
            lower, upper = set_to_bounds(c_set)
            # Note that there might be multiple bounds on the same variable
            # (e.g., LessThan and GreaterThan), so we should only update bounds
            # if they differ from the defaults.
            if lower > -Inf
                x_l[v_index] = lower
            end
            if upper < Inf
                x_u[v_index] = upper
            end
        end
    end

    # We have to convert all ScalarAffineFunction-in-Set constraints to an
    # expression graph.
    scalar_constraint_expr = Expr[]
    for func_type in MOI_SCALAR_FUNCTIONS
        for set_type in MOI_SCALAR_SETS
            for c_ref in MOI.get(model, MOI.ListOfConstraintIndices{
                    func_type, set_type}())
                c_func = MOI.get(model, MOI.ConstraintFunction(), c_ref)
                c_set = MOI.get(model, MOI.ConstraintSet(), c_ref)
                expr = moi_to_expr_graph(c_func, c_set, variable_map)
                push!(scalar_constraint_expr, expr)
                lower, upper = set_to_bounds(c_set)
                push!(g_l, lower)
                push!(g_u, upper)
            end
        end
    end

    # MOI objective
    obj_type = MOI.get(model, MOI.ObjectiveFunctionType())
    obj_func = MOI.get(model, MOI.ObjectiveFunction{obj_type}())
    obj_func_expr = func_to_expr_graph(obj_func, variable_map)
    if obj_func_expr == :(+ 0.0)
        obj_func_expr = nothing
    end

    # Build the nlp_evaluator
    nlp_evaluator = NLPEvaluator(
        moi_nlp_evaluator,
        variable_map,
        num_con,
        nlp_block === nothing ? false : nlp_block.has_objective,
        obj_func_expr,
        scalar_constraint_expr
    )

    # Create the MathProgBase model. Note that we pass `num_con` and the number
    # of linear constraints.
    mpb_model = MPB.NonlinearModel(mpb_solver)
    MPB.loadproblem!(mpb_model, num_var,
        num_con + length(scalar_constraint_expr), x_l, x_u, g_l, g_u, sense,
        nlp_evaluator)

    # Set any variables to :Bin if they are in ZeroOne and :Int if they are
    # Integer. The default is just :Cont.
    x_cat = fill(:Cont, num_var)
    for set_type in (MOI.ZeroOne, MOI.Integer)
        for c_ref in MOI.get(model,
            MOI.ListOfConstraintIndices{MOI.SingleVariable, set_type}())
            c_func = MOI.get(model, MOI.ConstraintFunction(), c_ref)
            c_set = MOI.get(model, MOI.ConstraintSet(), c_ref)
            v_index = variable_map[c_func.variable]
            x_cat[v_index] = set_to_cat(c_set)
        end
    end
    MPB.setvartype!(mpb_model, x_cat)

    # Set the VariablePrimalStart attributes for variables.
    variable_primal_start = fill(0.0, num_var)
    for (i, variable) in enumerate(variables)
        start_val = MOI.get(model, MOI.VariablePrimalStart(), variable)
        if start_val !== nothing
            variable_primal_start[i] = start_val
        end
    end
    MPB.setwarmstart!(mpb_model, variable_primal_start)

    # Solve the model!
    MPB.optimize!(mpb_model)

    # Extract and save the MathProgBase solution.
    primal_solution = Dict{MOI.VariableIndex, Float64}()
    for (variable, sol) in zip(variables, MPB.getsolution(mpb_model))
        primal_solution[variable] = sol
    end
    model.ext[:MPBModel] = mpb_model
    model.ext[:MPBSolutionAttribute] = MPBSolution(
        MPB.status(mpb_model),
        mpb_model.inner.solve_result,
        MPB.getobjval(mpb_model),
        primal_solution
    )
    return
end

# MOI accessors for the solution info.

function mpb_solution_attribute(model::Model)::Union{Nothing, MPBSolution}
    return get(model.ext, :MPBSolutionAttribute, nothing)
end

function MOI.get(model::Model, ::MOI.VariablePrimal, var::MOI.VariableIndex)
    mpb_solution = mpb_solution_attribute(model)
    if mpb_solution === nothing
        return
    end
    return mpb_solution.primal_solution[var]
end

function MOI.get(model::Model, ::MOI.ConstraintPrimal, idx::MOI.ConstraintIndex)
    return MOIU.get_fallback(model, MOI.ConstraintPrimal(), idx)
end

function MOI.get(model::Model, ::MOI.ObjectiveValue)
    mpb_solution = mpb_solution_attribute(model)
    if mpb_solution === nothing
        return
    end
    return mpb_solution.objective_value
end

function MOI.get(model::Model, ::MOI.TerminationStatus)
    mpb_solution = mpb_solution_attribute(model)
    if mpb_solution === nothing
        return MOI.OPTIMIZE_NOT_CALLED
    end
    status = mpb_solution.status
    if status == :Optimal
        # TODO(odow): this is not always the case. What if Ipopt solves a
        # convex problem?
        return MOI.LOCALLY_SOLVED
    elseif status == :Infeasible
        return MOI.INFEASIBLE
    elseif status == :Unbounded
        return MOI.DUAL_INFEASIBLE
    elseif status == :UserLimit
        return MOI.OTHER_LIMIT
    elseif status == :Error
        return MOI.OTHER_ERROR
    end
    return MOI.OTHER_ERROR
end

function MOI.get(model::Model, ::MOI.PrimalStatus)
    mpb_solution = mpb_solution_attribute(model)
    if mpb_solution === nothing
        return MOI.NO_SOLUTION
    end
    status = mpb_solution.status
    if status == :Optimal
        return MOI.FEASIBLE_POINT
    end
    return MOI.NO_SOLUTION
end

function MOI.get(::Model, ::MOI.DualStatus)
    return MOI.NO_SOLUTION
end

function MOI.get(model::Model, ::MOI.ResultCount)
    if MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
        return 1
    else
        return 0
    end
end

function MOI.get(model::Model, ::MOI.RawStatusString)
    mpb_solution = mpb_solution_attribute(model)
    return mpb_solution.raw_status_string
end

ext: jl


text:
using DataFrames
using GLM
using MixedModels
using Test

using MixedModels: dataset, likelihoodratiotest
using GLM: ProbitLink
using StatsModels: lrtest, isnested

include("modelcache.jl")

@testset "isnested" begin
    slp = dataset(:sleepstudy)

    # these tests don't actually depend on the models being fit,
    # so we just construct them

    # mismatched RE terms
    m1 = LinearMixedModel(@formula(reaction ~ 1 + days + (1+days|subj)), slp)
    m2 = LinearMixedModel(@formula(reaction ~ 1 + days + (0+days|subj)), slp)
    @test !isnested(m1, m2)

    # mismatched FE
    m1 = LinearMixedModel(@formula(reaction ~ 1 + days + (1|subj)), slp)
    m2 = LinearMixedModel(@formula(reaction ~ 0 + days + (1|subj)), slp)
    @test !isnested(m1, m2)

    # mismatched grouping vars
    kb07  = dataset(:kb07)
    m1 = LinearMixedModel(@formula(rt_trunc ~ 1 + (1|subj)), kb07)
    m2 = LinearMixedModel(@formula(rt_trunc ~ 1 + (1|item)), kb07)
    @test !isnested(m1, m2)

    # fixed-effects specification in REML and
    # conversion of internal ArgumentError into @error for StatsModels.isnested
    kb07  = dataset(:kb07)
    m1 = fit(MixedModel, @formula(rt_trunc ~ 1 + prec + (1|subj)), kb07, REML=true, progress=false)
    m2 = fit(MixedModel, @formula(rt_trunc ~ 1 + prec + (1+prec|subj)), kb07, REML=true, progress=false)
    @test isnested(m1, m2)
    m2 = fit(MixedModel, @formula(rt_trunc ~ 1 + (1+prec|subj)), kb07, REML=true, progress=false)
    @test !isnested(m1, m2)

end

@testset "likelihoodratio test" begin
    slp = dataset(:sleepstudy);

    fm0 = fit(MixedModel,@formula(reaction ~ 1 + (1+days|subj)),slp, progress=false);
    fm1 = fit(MixedModel,@formula(reaction ~ 1 + days + (1+days|subj)),slp, progress=false);
    lm0 = lm(@formula(reaction ~ 1), slp)
    lm1 = lm(@formula(reaction ~ 1 + days), slp)

    @test MixedModels._iscomparable(lm0, fm1)
    @test !MixedModels._iscomparable(lm1, fm0)

    lrt = likelihoodratiotest(fm0,fm1)

    @test [deviance(fm0), deviance(fm1)] == lrt.deviance
    @test deviance(fm0) - deviance(fm1) == only(lrt.tests.deviancediff)
    @test only(lrt.tests.dofdiff) == 1
    @test sum(map(length,lrt.tests)) == 3
    @test sum(map(length,lrt.pvalues)) == 1
    @test sum(map(length,lrt.models)) == 4
    @test length(lrt.formulae) == 2
    show(IOBuffer(),lrt);
    @test :pvalues in propertynames(lrt)

    lrt = likelihoodratiotest(lm1,fm1)
    @test lrt.deviance ≈ likelihoodratiotest(lm1.model,fm1).deviance
    @test lrt.dof == [3, 6]
    @test lrt.deviance ≈ -2 * loglikelihood.([lm1, fm1])
    shown = sprint(show, lrt)
    @test occursin("-2 logLik", shown)
    @test !occursin("deviance", shown)

    # non nested FE between non-mixed and mixed
    @test_throws ArgumentError likelihoodratiotest(lm1, fm0)

    # mix of REML and ML
    fm0 = fit(MixedModel,@formula(reaction ~ 1 + (1+days|subj)),slp, REML=true, progress=false);
    @test_throws ArgumentError likelihoodratiotest(fm0,fm1)
    @test_throws ArgumentError likelihoodratiotest(lm0,fm0)

    # differing FE with REML
    fm1 = fit(MixedModel,@formula(reaction ~ 1 + days + (1+days|subj)),slp, REML=true, progress=false);

    @test_throws ArgumentError likelihoodratiotest(fm0,fm1)
    contra = MixedModels.dataset(:contra);
    # glm doesn't like categorical responses, so we convert it to numeric ourselves
    # TODO: upstream fix
    cc = DataFrame(contra);
    cc.usenum = ifelse.(cc.use .== "Y", 1 , 0)
    gmf = glm(@formula(usenum ~ 1+age+urban+livch), cc, Bernoulli());
    gmf2 = glm(@formula(usenum ~ 1+age+abs2(age)+urban+livch), cc, Bernoulli());
    gm0 = fit(MixedModel, @formula(use ~ 1+age+urban+livch+(1|urban&dist)), contra, Bernoulli(), fast=true, progress=false);
    gm1 = fit(MixedModel, @formula(use ~ 1+age+abs2(age)+urban+livch+(1|urban&dist)), contra, Bernoulli(), fast=true, progress=false);

    lrt = likelihoodratiotest(gmf, gm1)
    @test [-2 * loglikelihood(gmf), deviance(gm1)] ≈ lrt.deviance
    @test -2 * loglikelihood(gmf) - deviance(gm1) ≈ only(lrt.tests.deviancediff)
    shown = sprint(show, lrt)
    @test !occursin("-2 logLik", shown)
    @test occursin("deviance", shown)

    lrt = likelihoodratiotest(gm0,gm1);
    @test [deviance(gm0), deviance(gm1)] == lrt.deviance
    @test deviance(gm0) - deviance(gm1) == only(lrt.tests.deviancediff)
    @test first(lrt.tests.dofdiff) == 1
    @test sum(length, lrt.tests) == 3
    @test sum(length, lrt.pvalues) == 1
    @test sum(length, lrt.models) == 4
    @test length(lrt.formulae) == 2

    # mismatched links
    gm_probit = fit(MixedModel, @formula(use ~ 1+age+urban+livch+(1|urban&dist)), contra, Bernoulli(), ProbitLink(), fast=true, progress=false);
    @test_throws ArgumentError likelihoodratiotest(gmf, gm_probit)
    @test_throws ArgumentError likelihoodratiotest(gm0, gm_probit)

    # mismatched families
    gm_poisson = fit(MixedModel, @formula(use ~ 1+age+urban+livch+(1|urban&dist)), contra, Poisson(), fast=true, progress=false);
    @test_throws ArgumentError likelihoodratiotest(gmf, gm_poisson)
    @test_throws ArgumentError likelihoodratiotest(gm0, gm_poisson)

    @test !MixedModels._iscomparable(lm0, gm0)
    @test !MixedModels._iscomparable(gmf, fm1)

    @test MixedModels._iscomparable(gmf, gm0)
    @test !MixedModels._iscomparable(gmf2, gm0)

    @test MixedModels._isnested(gmf.mm.m, gm0.X)
    @test !MixedModels._isnested(gmf2.mm.m, gm0.X)
    # this skips the linear term so that the model matrices
    # have the same column rank
    @test !MixedModels._isnested(gmf2.mm.m[:,Not(2)], gm0.X)
end


text:
"""
    LexicographicallyLessThan(row_dim::Int, column_dim::Int)

Ensures that each column of the matrix is lexicographically less than
the next column.

Formally, for two columns:

``\\{(x, y) \\in \\mathbb{R}^\\mathtt{column\\_dim} \\times \\mathbb{R}^\\mathtt{column\\_dim} | \\exists j \\in \\{1, 2 \\dots \\mathtt{column\\_dim}\\}: x_j < y_j, \\forall i < j, x_i = y_i \\}``.

Also called [`lex_less`](https://sofdem.github.io/gccat/gccat/Clex_less.html).

The matrix is encoded by stacking the columns, matching the behaviour of
Julia's `vec` function.
"""
struct LexicographicallyLessThan <: MOI.AbstractVectorSet
    row_dim::Int
    column_dim::Int
end

function LexicographicallyLessThan(column_dim::Int)
    # By default, only two columns.
    return LexicographicallyLessThan(2, column_dim)
end

MOI.dimension(set::LexicographicallyLessThan) = set.row_dim * set.column_dim

"""
    LexicographicallyGreaterThan(row_dim::Int, column_dim::Int)

Ensures that each column of the matrix is lexicographically greater than
the next column.

Formally, for two columns:

``\\{(x, y) \\in \\mathbb{R}^\\mathtt{column\\_dim} \\times \\mathbb{R}^\\mathtt{column\\_dim} | \exists j \\in \\{1, 2 \\dots \\mathtt{column\\_dim}\\}: x_j > y_j, \\forall i < j, x_i = y_i \\}``.

Also called [`lex_greater`](https://sofdem.github.io/gccat/gccat/Clex_greater.html).

The matrix is encoded by stacking the columns, matching the behaviour of
Julia's `vec` function.
"""
struct LexicographicallyGreaterThan <: MOI.AbstractVectorSet
    row_dim::Int
    column_dim::Int
end

function LexicographicallyGreaterThan(column_dim::Int)
    # By default, only two columns.
    return LexicographicallyGreaterThan(2, column_dim)
end

MOI.dimension(set::LexicographicallyGreaterThan) = set.row_dim * set.column_dim

# TODO: bridge to LexicographicallyLessThan.

"""
    DoublyLexicographicallyLessThan(dimension::Int)

Ensures that each column of the matrix is lexicographically less than
the next column, and that each row of the matrix is lexicographically less
than the next row.

Also called [`lex2`](https://sofdem.github.io/gccat/gccat/Clex2.html).

The matrix is encoded by stacking the columns, matching the behaviour of
Julia's `vec` function.
"""
struct DoublyLexicographicallyLessThan <: MOI.AbstractVectorSet
    row_dim::Int
    column_dim::Int
end

MOI.dimension(set::DoublyLexicographicallyLessThan) = set.row_dim * set.column_dim

"""
    DoublyLexicographicallyGreaterThan(dimension::Int)

Ensures that each column of the matrix is lexicographically greater than
the next column, and that each row of the matrix is lexicographically greater
than the next row.

The matrix is encoded by stacking the columns, matching the behaviour of
Julia's `vec` function.
"""
struct DoublyLexicographicallyGreaterThan <: MOI.AbstractVectorSet
    row_dim::Int
    column_dim::Int
end

MOI.dimension(set::DoublyLexicographicallyGreaterThan) = set.row_dim * set.column_dim

# TODO: bridge to LexicographicallyLessThan.

"""
    Sort(dimension::Int)

Ensures that the first `dimension` elements is a sorted copy of the next
`dimension` elements.

## Example

    [a, b, c, d] in Sort(2)
    # Enforces that:
    # - the first part is sorted: a <= b
    # - the first part corresponds to the second one:
    #     - either a = c and b = d
    #     - or a = d and b = c
"""
struct Sort <: MOI.AbstractVectorSet
    dimension::Int
end

MOI.dimension(set::Sort) = 2 * set.dimension

"""
    SortPermutation(dimension::Int)

Ensures that the first `dimension` elements is a sorted copy of the next
`dimension` elements.

The last `dimension` elements give a permutation to get from the original array
to its sorted version.

## Example

    [a, b, c, d, i, j] in SortPermutation(2)
    # Enforces that:
    # - the first part is sorted: a <= b
    # - the first part corresponds to the second one:
    #     - either a = c and b = d: in this case, i = 1 and j = 2
    #     - or a = d and b = c: in this case, i = 2 and j = 1
"""
struct SortPermutation <: MOI.AbstractVectorSet
    dimension::Int
end

MOI.dimension(set::SortPermutation) = 3 * set.dimension

"""
    MinimumAmong(dimension::Int)

Ensures that the first element is the minimum value among the next
`dimension` elements.

## Example

    [a, b, c] in MinimumAmong(2)
    # Enforces that a == min(b, c)
"""
struct MinimumAmong <: MOI.AbstractVectorSet
    dimension::Int
end

MOI.dimension(set::MinimumAmong) = 1 + set.dimension

"""
    ArgumentMinimumAmong(dimension::Int)

Ensures that the first element is the index of the minimum value among the
next `dimension` elements.

## Example

    [a, b, c] in ArgumentMinimumAmong(2)
    # Enforces that a == argmin(b, c)
    # I.e., if b < c, a = 1, if b > c, a = 2
"""
struct ArgumentMinimumAmong <: MOI.AbstractVectorSet
    dimension::Int
end

MOI.dimension(set::ArgumentMinimumAmong) = 1 + set.dimension

"""
    MaximumAmong(dimension::Int)

Ensures that the first element is the maximum value among the next
`dimension` elements.

## Example

    [a, b, c] in MaximumAmong(2)
    # Enforces that a == max(b, c)
"""
struct MaximumAmong <: MOI.AbstractVectorSet
    dimension::Int
end

MOI.dimension(set::MaximumAmong) = 1 + set.dimension

"""
    ArgumentMaximumAmong(dimension::Int)

Ensures that the first element is the index of the maximum value among the
next `dimension` elements.

## Example

    [a, b, c] in ArgumentMaximumAmong(2)
    # Enforces that a == argmax(b, c)
    # I.e., if b > c, a = 1, if b < c, a = 2
"""
struct ArgumentMaximumAmong <: MOI.AbstractVectorSet
    dimension::Int
end

MOI.dimension(set::ArgumentMaximumAmong) = 1 + set.dimension

"""
    Increasing(dimension::Int)

Ensures that the elements of the vector are in increasing order (<= operation).

## Example

    [a, b, c] in Increasing(3)
    # Enforces that a <= b <= c
"""
struct Increasing <: MOI.AbstractVectorSet
    dimension::Int
end

"""
    Decreasing(dimension::Int)

Ensures that the elements of the vector are in decreasing order (>= operation).

## Example

    [a, b, c] in Decreasing(3)
    # Enforces that a >= b >= c
"""
struct Decreasing <: MOI.AbstractVectorSet
    dimension::Int
end

# isbits types, nothing to copy
function copy(
    set::Union{
        LexicographicallyLessThan,
        LexicographicallyGreaterThan,
        DoublyLexicographicallyLessThan,
        DoublyLexicographicallyGreaterThan,
        Sort,
        SortPermutation,
        MinimumAmong,
        ArgumentMinimumAmong,
        MaximumAmong,
        ArgumentMaximumAmong,
        Increasing,
        Decreasing,
    },
)
    return set
end


text:
export fft_view_tuple

"""
    fft_view_tuple(x[, shift, dims])

"""
function fft_view_tuple(x, shift=true, dims=nothing; logarithmic=false)
    if isnothing(dims)
        dims = size(x)
    end
    if shift
        x_fft = fftshift(fft(x))
    else
        x_fft = fft(x)
    end

    x_fft_i = imag(x_fft)
    x_fft_r = real(x_fft)
    x_fft_abs = abs.(x_fft)

    if logarithmic
        l(x) = log(1 + abs.(x))
    else
        l = identity
    end

    return l.(x_fft_i), l.(x_fft_r), l.(x_fft_abs)
end

max_stars_repo_name: JuliaPhysics/PhysicalOptics

max_issues_repo_name: JuliaPhysics/PhysicalOptics

max_forks_repo_name: JuliaPhysics/PhysicalOptics


text:
#==========================================================================================
Conversions between DataTables and DataFrames
==========================================================================================#
"""
    convertArray(array_type, a)

Convert between `NullableArray`s and `DataArray`s.  The former is used by `DataTables` while
the latter is used by `DataFrames`.

Note that this is named `convertArray` rather than simply `convert` so as not to conflict
with the existing definitions for `convert`.
"""
function convertArray{T,N}(::Type{DataArrays.DataArray}, v::NullableArray{T,N})
    DataArrays.DataArray{T,N}(copy(v.values), BitArray(v.isnull))
end

function convertArray{T,N}(::Type{NullableArray}, v::DataArrays.DataArray{T,N})
    NullableArray{T,N}(copy(v.data), BitArray(v.na))
end

convertArray(::Type{DataArrays.DataArray}, v::Array) = DataArrays.DataArray(v)

convertArray(::Type{NullableArray}, v::Array) = NullableArray(v)

# TODO this will no longer be needed once Feather is updated
convertArray(::Type{NullableArray}, v::NullableArray) = copy(v)
export convertArray

function convert(::Type{DataFrames.DataFrame}, dt::DataTable)
    df = DataFrames.DataFrame()
    for name ∈ names(dt)
        df[name] = convertArray(DataArrays.DataArray, dt[name])
    end
    df
end

function convert(::Type{DataTable}, df::DataFrames.DataFrame)
    dt = DataTable()
    for name ∈ names(df)
        dt[name] = convertArray(NullableArray, df[name])
    end
    dt
end
export convert


text:
MODULES  = (:LaxWendroff,
            :Upwind,
            :Godunov,
            :SemiLagrangian,
            :PFC)

for mod_name in MODULES
    include(joinpath(dirname(@__FILE__),"..","..","src","VlasovSolver","$mod_name.jl"))
end

function test_1d_advection_step(mod_name, Δx, Δt, v, f₀, f₁, exp_norm_dev, args...; plot_needed=false, kwargs...)
    if plot_needed
        plt = plot(f₀, label = "initial")
        plot!(f₁, label = "expected")
    end

    f = similar(f₀)

    advect! = eval(:($mod_name)).generate_solver(f₀, f, v*Δt/Δx, args...; kwargs...)
    advect!()

    if plot_needed
        plot!(f, label = "$mod_name $args (constant velocity)")
    end

    println("$mod_name $args $(values(values(kwargs))) (constant velocity): $(Δx*norm(f - f₁))")
    @test Δx*norm(f - f₁) ≈ 0 atol=exp_norm_dev

    advect! = eval(:($mod_name)).generate_solver(f₀, f, args...; kwargs...)
    advect!(v*Δt/Δx)

    if plot_needed
        plot!(f, label = "$mod_name $args")
    end

    println("$mod_name $args $(values(values(kwargs))): $(Δx*norm(f - f₁))")
    @test Δx*norm(f - f₁) ≈ 0 atol=exp_norm_dev

    if plot_needed
        display(plt)
    end
end

@testset "Test 1D advection solvers" begin
    Δx = 0.01
    Δt = 0.8*Δx
    v = 1.0

    f₀ = [1.0 + 0.01*sin(2π*i*Δx) for i = 0:99]
    f₁ = [1.0 + 0.01*sin(2π*(i*Δx - v*Δt)) for i = 0:99]

    test_1d_advection_step(:Upwind, Δx, Δt,  v, f₀, f₁, 3e-7)
    test_1d_advection_step(:Upwind, Δx, Δt, -v, f₁, f₀, 3e-7)

    test_1d_advection_step(:LaxWendroff, Δx, Δt,  v, f₀, f₁, 1e-8)
    test_1d_advection_step(:LaxWendroff, Δx, Δt, -v, f₁, f₀, 1e-8)

    for riemann_solver in (:Riemann_constant, :Riemann_linear)
        test_1d_advection_step(:Godunov, Δx, Δt,  v, f₀, f₁, 1e-6, riemann_solver)
        test_1d_advection_step(:Godunov, Δx, Δt, -v, f₁, f₀, 1e-6, riemann_solver)
        if riemann_solver != :Riemann_constant
            for flux_limiter in (:VanLeer,)
                test_1d_advection_step(:Godunov, Δx, Δt,  v, f₀, f₁, 1e-6, riemann_solver; flux_limiter=flux_limiter)
                test_1d_advection_step(:Godunov, Δx, Δt, -v, f₁, f₀, 1e-6, riemann_solver; flux_limiter=flux_limiter)
            end
        end
    end

    test_1d_advection_step(:SemiLagrangian, Δx, Δt,  v, f₀, f₁, 3e-7; interpolation_order = :Linear)
    test_1d_advection_step(:SemiLagrangian, Δx, Δt, -v, f₁, f₀, 3e-7; interpolation_order = :Linear)

    test_1d_advection_step(:SemiLagrangian, Δx, Δt,  v, f₀, f₁, 2e-9; interpolation_order = :Quadratic)
    test_1d_advection_step(:SemiLagrangian, Δx, Δt, -v, f₁, f₀, 2e-9; interpolation_order = :Quadratic)

    test_1d_advection_step(:SemiLagrangian, Δx, Δt,  v, f₀, f₁, 2e-11; interpolation_order = :Cubic)
    test_1d_advection_step(:SemiLagrangian, Δx, Δt, -v, f₁, f₀, 2e-11; interpolation_order = :Cubic)

    test_1d_advection_step(:PFC, Δx, Δt,  v, f₀, f₁, 3e-8)
    test_1d_advection_step(:PFC, Δx, Δt, -v, f₁, f₀, 3e-8)

    f₀ = [exp(-((i*Δx - 0.5)/0.15)^2) for i = 0:100]
    f₁ = [exp(-((i*Δx - v*Δt - 0.5)/0.15)^2) for i = 0:100]

    test_1d_advection_step(:Upwind, Δx, Δt,  v, f₀, f₁, 1e-4)
    test_1d_advection_step(:Upwind, Δx, Δt, -v, f₁, f₀, 1e-4)

    test_1d_advection_step(:LaxWendroff, Δx, Δt,  v, f₀, f₁, 1e-5)
    test_1d_advection_step(:LaxWendroff, Δx, Δt, -v, f₁, f₀, 1e-5)

    test_1d_advection_step(:Godunov, Δx, Δt,  v, f₀, f₁, 1e-4, :Riemann_constant)
    test_1d_advection_step(:Godunov, Δx, Δt, -v, f₁, f₀, 1e-4, :Riemann_constant)

    test_1d_advection_step(:Godunov, Δx, Δt,  v, f₀, f₁, 1.1e-4, :Riemann_linear)
    test_1d_advection_step(:Godunov, Δx, Δt, -v, f₁, f₀, 1.1e-4, :Riemann_linear)

    test_1d_advection_step(:Godunov, Δx, Δt,  v, f₀, f₁, 1.2e-4, :Riemann_linear; flux_limiter=:VanLeer)
    test_1d_advection_step(:Godunov, Δx, Δt, -v, f₁, f₀, 1.2e-4, :Riemann_linear; flux_limiter=:VanLeer)

    test_1d_advection_step(:SemiLagrangian, Δx, Δt,  v, f₀, f₁, 3e-5; interpolation_order = :Linear)
    test_1d_advection_step(:SemiLagrangian, Δx, Δt, -v, f₁, f₀, 3e-5; interpolation_order = :Linear)

    test_1d_advection_step(:SemiLagrangian, Δx, Δt,  v, f₀, f₁, 4e-7; interpolation_order = :Quadratic)
    test_1d_advection_step(:SemiLagrangian, Δx, Δt, -v, f₁, f₀, 4e-7; interpolation_order = :Quadratic)

    test_1d_advection_step(:SemiLagrangian, Δx, Δt,  v, f₀, f₁, 4e-8; interpolation_order = :Cubic)
    test_1d_advection_step(:SemiLagrangian, Δx, Δt, -v, f₁, f₀, 4e-8; interpolation_order = :Cubic)

    test_1d_advection_step(:PFC, Δx, Δt,  v, f₀, f₁, 5e-6)
    test_1d_advection_step(:PFC, Δx, Δt, -v, f₁, f₀, 5e-6)

    f₀ = zeros(Float64, 100)
    for i = 40:50
        f₀[i] = 1.0
    end
    f₁ = zeros(Float64, 100)
    for i = 41:50
        f₁[i] = 1.0
    end
    f₁[40] = 0.2
    f₁[51] = 0.8
    f₂ = zeros(Float64, 100)
    for i = 40:49
        f₂[i] = 1.0
    end
    f₂[39] = 0.8
    f₂[50] = 0.2

    test_1d_advection_step(:Upwind, Δx, Δt,  v, f₀, f₁, 1e-18)
    test_1d_advection_step(:Upwind, Δx, Δt, -v, f₀, f₂, 1e-18)

    test_1d_advection_step(:LaxWendroff, Δx, Δt,  v, f₀, f₁, 1.6e-3)
    test_1d_advection_step(:LaxWendroff, Δx, Δt, -v, f₀, f₂, 1.6e-3)

    test_1d_advection_step(:Godunov, Δx, Δt,  v, f₀, f₁, 1e-18, :Riemann_constant)
    test_1d_advection_step(:Godunov, Δx, Δt, -v, f₀, f₂, 1e-18, :Riemann_constant)

    test_1d_advection_step(:Godunov, Δx, Δt,  v, f₀, f₁, 1e-2, :Riemann_linear)
    test_1d_advection_step(:Godunov, Δx, Δt, -v, f₀, f₂, 1e-2, :Riemann_linear)

    test_1d_advection_step(:Godunov, Δx, Δt,  v, f₀, f₁, 1e-2, :Riemann_linear; flux_limiter=:VanLeer)
    test_1d_advection_step(:Godunov, Δx, Δt, -v, f₀, f₂, 1e-2, :Riemann_linear; flux_limiter=:VanLeer)

    test_1d_advection_step(:SemiLagrangian, Δx, Δt,  v, f₀, f₁, 5e-17;interpolation_order = :Linear)
    test_1d_advection_step(:SemiLagrangian, Δx, Δt, -v, f₀, f₂, 5e-17; interpolation_order = :Linear)

    test_1d_advection_step(:SemiLagrangian, Δx, Δt,  v, f₀, f₁, 2e-3; interpolation_order = :Quadratic)
    test_1d_advection_step(:SemiLagrangian, Δx, Δt, -v, f₀, f₂, 2e-3; interpolation_order = :Quadratic)

    test_1d_advection_step(:SemiLagrangian, Δx, Δt,  v, f₀, f₁, 2e-3; interpolation_order = :Cubic)
    test_1d_advection_step(:SemiLagrangian, Δx, Δt, -v, f₀, f₂, 2e-3; interpolation_order = :Cubic)

    test_1d_advection_step(:PFC, Δx, Δt,  v, f₀, f₁, 1e-18)
    test_1d_advection_step(:PFC, Δx, Δt, -v, f₀, f₂, 1e-18)
end


text:
using Plots

include("../src/parzenwindow.jl")

# Data generation
data1 = ones(20)*2+randn(20)
data2 = ones(10)*6+randn(10)
data = vcat(data1, data2)
data = reshape(data, (1, size(data)[1]))

# Plot data
plt = plot(title = "Parzen window estimation (window size 'h')", xlabel="Data value", ylabel="Likelihood")
scatter!(data1, zeros(20), color="black", label="Data")
scatter!(data2, zeros(10), color="black", label = "")

# Plot parzen window function
colors = ["red", "green", "blue", "black"]
kernel(x) = k(x)
for step in 0.1:0.5:1.6
    Y = []
    model(x) = create_parzen_window(step, data, kernel, x)
    for i in minimum(data):0.01:maximum(data)
        append!(Y, model(Vector{Float64}([i])))
    end
    plot!(plt, minimum(data):0.01:maximum(data), Y, color=colors[Int((step+0.4)*2)], label=string("h = ", string(step)))
end
plt

max_stars_repo_name: kozvojtex/AnomalyDetection

max_issues_repo_name: kozvojtex/AnomalyDetection

max_forks_repo_name: kozvojtex/AnomalyDetection


text:
module NeuralNets
using SparseArrays
using DataFrames
using LightGraphs
using MetaGraphs
#using Query
using ..Neurons
using ..Neurons.Segments.Synapses
using ..SWCs

using Query

export NeuralNet
const NeuralNet = MetaDiGraph

"""
    NeuralNet(neuronId2neuron::Dict{Int, Neuron}; neuronIdList::Vector{Int})
construct neuron graph.
Parameters:
    neuronId2neuron: the dictionary mapping ID to neuron
    neuronIdList: the vertex id was ordered according to this list
Return:
    net: the MetaDiGraph representing neural network.
"""
function NeuralNet(neuronList::Vector{Neuron{T}}, neuronIdList::Vector{Int}) where T
    neuronIdSet = Set{Int}(neuronIdList)
    net = NeuralNet(length(neuronIdList))

    # name the vertices
    neuronId2vertexId = Dict{Int,Int}()
    for (vertexId, neuronId) in enumerate(neuronIdList)
        neuronId2vertexId[neuronId] = vertexId
        neuron = neuronList[ vertexId ]
        set_props!(net, vertexId, Dict(:id=>neuronId, :skeleton=>neuron))
    end

    # add edges
    selfConnectNeuronIdSet = Set{Int}()
    for (i,neuronId) in enumerate(neuronIdList)
        neuron = neuronList[ i ]
        vertexId = neuronId2vertexId[ neuronId ]
        preSynapseList = Neurons.get_all_pre_synapse_list(neuron)
        for preSynapse in preSynapseList
            postNeuronId = Synapses.get_post_synaptic_segmentation_id(preSynapse)
            if postNeuronId == neuronId
                @warn("skipping self connection edge of neuron: ", neuronId)
                push!(selfConnectNeuronIdSet, neuronId)
            elseif postNeuronId in neuronIdSet
                # this is an effective synaptic edge
                postNeuronVertexId = neuronId2vertexId[ postNeuronId ]
                edge = (vertexId, postNeuronVertexId)
                add_edge!(net, edge)
                # get properties
                if has_prop(net, edge[1], edge[2], :synapses)
                    prop = get_prop(net, edge[1], edge[2], :synapses)
                    totalPSDSize = get_prop(net, edge[1], edge[2], :total_psd_size)
                else
                    prop = Vector{Synapse{T}}()
                    totalPSDSize = 0
                end
                push!(prop, preSynapse)
                set_prop!(net, edge[1], edge[2], :synapses, prop)
                set_prop!(net, edge[1], edge[2], :num_synapses, length(prop))
                totalPSDSize += Synapses.get_psd_size(preSynapse)
                set_prop!(net, edge[1], edge[2], :total_psd_size, totalPSDSize)
            end
        end
    end
    @warn("skipped self connection edges from $(length(selfConnectNeuronIdSet)) neurons.")
    net
end

function NeuralNet( syn::DataFrame; neuronDict::Dict{Int, Neuron}=Dict{Int,Neuron}(),
                                    cellIdList=get_cell_id_list(syn) )
    N = length(cellIdList)
    net = NeuralNet(N)
    # name the vertices
    for (i,cellId) in enumerate(cellIdList)
        if haskey(neuronDict, cellId)
            set_props!(net, i, Dict(:id=>cellId, :skeleton=>neuronDict[cellId]))
        else
            @warn("no skeleton attatched: ", cellId)
        end
    end

    # filter out the valid set
    cellIdSet = Set(cellIdList)
	syn =  @from i in syn begin
        @where 	i.presyn_segid in cellIdSet &&
                i.postsyn_segid in cellIdSet
        @select i
        @collect DataFrame
    end
    # construct a look up table
    cellId2vertexIdMap = Dict{Int, Int}()
    for (i,cellId) in enumerate(cellIdList)
        cellId2vertexIdMap[cellId] = i
    end
    # add edges
    numSelfConnection = 0
    for row in eachrow(syn)
        preId  = round(Int, row[:presyn_segid])
        postId = round(Int, row[:postsyn_segid])
        if preId == postId
            numSelfConnection += 1
            @warn("self connection synapse on neuron: $(preId), skipping the connection")
            continue
        end
        edge = ( cellId2vertexIdMap[preId], cellId2vertexIdMap[postId] )
        add_edge!(net, edge)
        # get properties
        if has_prop(net, edge[1], edge[2], :synapses)
            prop = get_prop(net, edge[1], edge[2], :synapses)
        else
            prop = Vector{Synapse}()
        end
        synapse = Synapse(row)
        push!(prop, synapse )
        set_prop!(net, edge[1], edge[2], :synapses, prop)
    end
    println("number of self connection synapses: $(numSelfConnection)")
    net
end

function get_neuron_id_list(self::NeuralNet)
    neuronIdList = Vector{Int}()
    sizehint!(neuronIdList, nv(self))

    for v in vertices(self)
        push!(neuronIdList, get_prop(self, v, :id))
    end
    neuronIdList
end

"""
    get_synapse_number_connectivity_matrix(self::NeuralNet)
Return:
    sparse matrix representing the connectivity matrix.
    The rows are presynase and the columns are post synapses.
    For example, [i,j] means the axon of neuron i innervate the dendrite of neuron j
"""
function get_synapse_number_connectivity_matrix(self::NeuralNet)
    N = nv(self)
    ret = spzeros(Int, N,N)
    for edge in edges(self)
        synapseList = get_prop(self, edge, :synapses)
        ret[src(edge), dst(edge)] = length(synapseList)
    end
    ret
end

"""
    get_psd_size_connectivity_matrix(self::NeuralNet)
Return:
    sparse matrix representing the connectivity matrix.
    The rows are presynapse and the columns are post synapses.
    For example [i,j] means the axon of neuron i innvervate the dendrite of neuron j.
    The value represent the total volume of psd
"""
function get_psd_size_connectivity_matrix(self::NeuralNet)
    N = nv(self)
    ret = spzeros(Int, N,N)
    for edge in edges(self)
        synapses = get_prop(self, edge, :synapses)
        totalSynSize = 0
        for synapse in synapses
            # note that the size was measured using voxel number
            totalSynSize += Synapses.get_psd_size(synapse)
        end
        ret[src(edge), dst(edge)] = totalSynSize
    end
    ret
end

"""
    labels2orders(label::Vector)

transform clustering label to orders

Parameters:
    labels: the assignment of each neuron. The label must start from 1 to N, where N is the maximum number of groups.
"""
function labels2orders(labels::Vector)
    orders = zero(labels)
    i = 0
    @inbounds for l in unique(labels)
        idxes = findall(x->x==l, labels)
        for idx in idxes
            i += 1
            orders[idx] = i
        end
    end
    return orders
end

"""
reorder_connectivity_matrix(connMatrix::Union{SparseMatrixCSC{T}, Matrix{T}}, orders::Vector)

reorder the rows and columns of the connectivity matrix for visualization.

Parameters:
    connMatrix: the connectivity matrix
    orders: the rearranged order of samples.

Return:
    a reordered connectivity matrix which is also a SparseMatrixCSC
"""
function reorder_matrix(connMatrix::Union{SparseMatrixCSC{T},Matrix{T}}, orders::Vector) where T

    # the reordered connectivity Matrix
    if isa(connMatrix, SparseMatrixCSC)
        ret = spzeros(T, size(connMatrix)...,)
    else
        ret = zeros(T, size(connMatrix))
    end

    @inbounds for (i, o1) in enumerate(orders)
        for (j, o2) in enumerate(orders)
            if connMatrix[o1,o2] !== zero(T)
                ret[i,j] = connMatrix[o1,o2]
            end
        end
    end
    return ret
end

function remove_singletons(self::NeuralNet)
    @error("unimplemented")
end

function get_nblast_similarity_matrix(self::NeuralNet)
    error("not implemented.")
end

end # end of module


text:
using Flux, Gym, Printf, Zygote
using Flux.Optimise: update!
using Statistics: mean
#using CuArrays

#Load game environment

env = make("Pendulum-v0")
reset!(env)
# ----------------------------- Parameters -------------------------------------

STATE_SIZE = length(state(env)) # returns state from obs space
ACTION_SIZE = 1#length(env.actions)
ACTION_BOUND = env._env.action_space.high[1]
MAX_REWARD = 0f0 # Max reward in a timestep
MAX_EP = 10
MAX_EP_LENGTH = 1000
SEQ_LEN = 4

# ------------------------------ Model Architecture ----------------------------

model = Chain(Dense(STATE_SIZE, 24, relu),
              Dense(24, 48, relu),
              Dense(48, ACTION_SIZE)) |> gpu

η = 3f-2

opt = ADAM(η)

loss(r) = Flux.mse(r, MAX_REWARD)

# ----------------------------- Helper Functions -------------------------------

function μEpisode(env::EnvWrapper)
    l = 0
    for frames ∈ 1:SEQ_LEN
        #render(env, ctx)
        #sleep(0.01)
        a = model(state(env))
        s, r, done, _ = step!(env, a)
        if trainable(env)
            l += loss(r)
        end

        game_over(env) && break
    end
    return l
end

function episode!(env::EnvWrapper)
  reset!(env)
  while !game_over(env)
    if trainable(env)
      grads = gradient(()->μEpisode(env), params(model))
      update!(opt, params(model), grads)
    else
      μEpisode(env)
    end
  end

  env.total_reward
end

# -------------------------------- Testing -------------------------------------

function test(env::EnvWrapper)
  score_mean = 0f0
  testmode!(env)
  for e=1:100
    total_reward = episode!(env)
    score_mean += total_reward / 100
  end
  testmode!(env, false)
  return score_mean
end

# ------------------------------ Training --------------------------------------

for e=1:MAX_EP
  total_reward = episode!(env)
  total_reward = @sprintf "%9.3f" total_reward
  print("Episode: $e | Score: $total_reward | ")
  score_mean = test(env)
  score_mean = @sprintf "%9.3f" score_mean
  println("Mean score over 100 test episodes: $score_mean")
end

max_stars_repo_name: jstrube/model-zoo

max_issues_repo_name: jstrube/model-zoo

max_forks_repo_name: jstrube/model-zoo


text:
export run_acdcscopf, run_acdcscopf_nocl, run_acdcscopf_nocl_PLdim

""
function run_acdcscopf(file::String, model_type::Type, solver; kwargs...)
    data = _PM.parse_file(file)
    PowerModelsACDC.process_additional_data!(data)
    return run_acdcscopf(data, model_type, solver; ref_extensions = [add_ref_dcgrid!], kwargs...)
end

""
function run_acdcscopf(data::Dict{String,Any}, model_type::Type, solver; kwargs...)
    return _PM.run_model(data, model_type, solver, post_acdcscopf; ref_extensions = [add_ref_dcgrid!], kwargs...)
    # return PowerModels.optimize_model!(pm, solver; solution_builder = get_solution_acdc)
end
""
function run_acdcscopf_nocl(file::String, model_type::Type, solver; kwargs...)
    data = _PM.parse_file(file)
    PowerModelsACDC.process_additional_data!(data)
    return run_acdcscopf(data, model_type, solver; ref_extensions = [add_ref_dcgrid!], kwargs...)
end

""
function run_acdcscopf_nocl(data::Dict{String,Any}, model_type::Type, solver; kwargs...)
    return _PM.run_model(data, model_type, solver, post_acdcscopf_nocl; ref_extensions = [add_ref_dcgrid!], kwargs...)
    # return PowerModels.optimize_model!(pm, solver; solution_builder = get_solution_acdc)
end

""
function run_acdcscopf_nocl_PLdim(file::String, model_type::Type, solver; kwargs...)
    data = _PM.parse_file(file)
    PowerModelsACDC.process_additional_data!(data)
    return run_acdcscopf(data, model_type, solver; ref_extensions = [add_ref_dcgrid!], kwargs...)
end

""
function run_acdcscopf_nocl_PLdim(data::Dict{String,Any}, model_type::Type, solver; kwargs...)
    return _PM.run_model(data, model_type, solver, post_acdcscopf_nocl_PLdim; ref_extensions = [add_ref_dcgrid!], kwargs...)
    # return PowerModels.optimize_model!(pm, solver; solution_builder = get_solution_acdc)
end

""
function post_acdcscopf(pm::_PM.AbstractPowerModel)

     for (n, networks) in pm.ref[:nw]
    _PM.variable_bus_voltage(pm; nw = n)
    _PM.variable_gen_power(pm; nw = n)
    _PM.variable_branch_power(pm; nw = n)

    variable_active_dcbranch_flow(pm; nw = n)
    variable_dcbranch_current(pm; nw = n)
    variable_dc_converter(pm; nw = n)
    variable_dcgrid_voltage_magnitude(pm; nw = n)
    variable_frequency_stab(pm; nw = n)
    end

    # for (n, networks) in pm.ref[:nw]
    objective_min_cost_OPF(pm)
    # end

    for (n, networks) in pm.ref[:nw]
    _PM.constraint_model_voltage(pm; nw = n)
    constraint_voltage_dc(pm; nw = n)

    for i in _PM.ids(pm, :ref_buses)
        _PM.constraint_theta_ref(pm, i; nw = n)
    end

    for i in _PM.ids(pm, :bus)
        constraint_power_balance_ac(pm, i; nw = n)
    end

    for i in _PM.ids(pm, :branch)
        _PM.constraint_ohms_yt_from(pm, i; nw = n)
        _PM.constraint_ohms_yt_to(pm, i; nw = n)
        _PM.constraint_voltage_angle_difference(pm, i; nw = n) #angle difference across transformer and reactor - useful for LPAC if available?
        _PM.constraint_thermal_limit_from(pm, i; nw = n)
        _PM.constraint_thermal_limit_to(pm, i; nw = n)
    end

    for i in _PM.ids(pm, :busdc)
        constraint_power_balance_dc(pm, i; nw = n)
    end
    for i in _PM.ids(pm, nw = n, :branchdc)
        constraint_ohms_dc_branch(pm, i; nw = n)
    end
    for i in _PM.ids(pm, :convdc)
        constraint_converter_losses(pm, i; nw = n)
        constraint_converter_current(pm, i; nw = n)
        constraint_conv_transformer(pm, i; nw = n)
        constraint_conv_reactor(pm, i; nw = n)
        constraint_conv_filter(pm, i; nw = n)
        if pm.ref[:nw][pm.cnw][:convdc][i]["islcc"] == 1
            constraint_conv_firing_angle(pm, i; nw = n)
        end
    end
 end
     constraint_frequency_stab_OPF(pm)
end

function post_acdcscopf_nocl(pm::_PM.AbstractPowerModel)

     for (n, networks) in pm.ref[:nw]
    _PM.variable_bus_voltage(pm; nw = n)
    _PM.variable_gen_power(pm; nw = n)
    _PM.variable_branch_power(pm; nw = n)

    variable_active_dcbranch_flow(pm; nw = n)
    variable_dcbranch_current(pm; nw = n)
    variable_dc_converter(pm; nw = n)
    variable_dcgrid_voltage_magnitude(pm; nw = n)
    variable_frequency_stab(pm; nw = n)
    end

    # for (n, networks) in pm.ref[:nw]
    objective_min_cost_OPF_nocl(pm)
    # end

    for (n, networks) in pm.ref[:nw]
    _PM.constraint_model_voltage(pm; nw = n)
    constraint_voltage_dc(pm; nw = n)

    for i in _PM.ids(pm, :ref_buses)
        _PM.constraint_theta_ref(pm, i; nw = n)
    end

    for i in _PM.ids(pm, :bus)
        constraint_power_balance_ac(pm, i; nw = n)
    end

    for i in _PM.ids(pm, :branch)
        _PM.constraint_ohms_yt_from(pm, i; nw = n)
        _PM.constraint_ohms_yt_to(pm, i; nw = n)
        _PM.constraint_voltage_angle_difference(pm, i; nw = n) #angle difference across transformer and reactor - useful for LPAC if available?
        _PM.constraint_thermal_limit_from(pm, i; nw = n)
        _PM.constraint_thermal_limit_to(pm, i; nw = n)
    end

    for i in _PM.ids(pm, :busdc)
        constraint_power_balance_dc(pm, i; nw = n)
    end
    for i in _PM.ids(pm, nw = n, :branchdc)
        constraint_ohms_dc_branch(pm, i; nw = n)
    end
    for i in _PM.ids(pm, :convdc)
        constraint_converter_losses(pm, i; nw = n)
        constraint_converter_current(pm, i; nw = n)
        constraint_conv_transformer(pm, i; nw = n)
        constraint_conv_reactor(pm, i; nw = n)
        constraint_conv_filter(pm, i; nw = n)
        if pm.ref[:nw][pm.cnw][:convdc][i]["islcc"] == 1
            constraint_conv_firing_angle(pm, i; nw = n)
        end
    end
 end
     constraint_frequency_stab_OPF(pm)
end

function post_acdcscopf_nocl_PLdim(pm::_PM.AbstractPowerModel)
    display("post_acdcscopf_nocl_PLdim")
     for (n, networks) in pm.ref[:nw]
    _PM.variable_bus_voltage(pm; nw = n)
    _PM.variable_gen_power(pm; nw = n)
    _PM.variable_branch_power(pm; nw = n)

    variable_active_dcbranch_flow(pm; nw = n)
    variable_dcbranch_current(pm; nw = n)
    variable_dc_converter(pm; nw = n)
    variable_dcgrid_voltage_magnitude(pm; nw = n)
    variable_frequency_stab(pm; nw = n)
    end

    # for (n, networks) in pm.ref[:nw]
    objective_min_cost_OPF_nocl_PLdim(pm)
    # end

    for (n, networks) in pm.ref[:nw]
    _PM.constraint_model_voltage(pm; nw = n)
    constraint_voltage_dc(pm; nw = n)

    for i in _PM.ids(pm, :ref_buses)
        _PM.constraint_theta_ref(pm, i; nw = n)
    end

    for i in _PM.ids(pm, :bus)
        constraint_power_balance_ac(pm, i; nw = n)
    end

    for i in _PM.ids(pm, :branch)
        _PM.constraint_ohms_yt_from(pm, i; nw = n)
        _PM.constraint_ohms_yt_to(pm, i; nw = n)
        _PM.constraint_voltage_angle_difference(pm, i; nw = n) #angle difference across transformer and reactor - useful for LPAC if available?
        _PM.constraint_thermal_limit_from(pm, i; nw = n)
        _PM.constraint_thermal_limit_to(pm, i; nw = n)
    end

    for i in _PM.ids(pm, :busdc)
        constraint_power_balance_dc(pm, i; nw = n)
    end
    for i in _PM.ids(pm, nw = n, :branchdc)
        constraint_ohms_dc_branch(pm, i; nw = n)
    end
    for i in _PM.ids(pm, :convdc)
        constraint_converter_losses(pm, i; nw = n)
        constraint_converter_current(pm, i; nw = n)
        constraint_conv_transformer(pm, i; nw = n)
        constraint_conv_reactor(pm, i; nw = n)
        constraint_conv_filter(pm, i; nw = n)
        if pm.ref[:nw][pm.cnw][:convdc][i]["islcc"] == 1
            constraint_conv_firing_angle(pm, i; nw = n)
        end
    end
 end
     constraint_frequency_stab_OPF_MP_PLdim(pm)
end


text:
"""
    Weak classifier `h`

structure representing a weak classifier which selects class for a measurement with given weight

# Fields
- `feature`: index in a feature vector defining by which part the decision is made
- `threshold`: separating value for given `feature`
- `bigger`: if set to ``true``, the feature value is compared by ``≥`` with the `threshold`, otherwise by `<`
- `weight`: how much this weak classifier contributes to overall decision
"""
struct WeakClassifier
    feature::Integer
    threshold::Real
    bigger::Bool
    weight::Real

    """
        WeakClassifier(feature::Integer, threshold::Real, bigger::Bool, error::Real)

    constructs weak classifer, which is meant to be saved into strong classifier

    how many data samples were misclassified is given by `error` argument, which determines the `weight` of the weak classifier
    """
    function WeakClassifier(feature::Integer, threshold::Real, bigger::Bool, error::Real)
        weight = log((1 - error) / error) / 2
        return new(feature, threshold, bigger, weight)
    end

    """
        WeakClassifier(feature::Integer, threshold::Real, bigger::Bool, error::Real)

    constructs weak classifer, which is meant for finding the best weak classifier on current data and does not incorporate the `weight calculation`
    """
    function WeakClassifier(feature::Integer, threshold::Real, bigger::Bool)
        return new(feature, threshold, bigger, 0.0)
    end
end

function Base.show(io::IO, h::WeakClassifier)
    print(io, "h: x[$(h.feature)] $(h.bigger ? ">=" : "<") $(h.threshold) [α = $(h.weight)]")
end

"""
    (h::WeakClassifier)(x::Vector{<:Real})

classifies one data sample `x` using weak classifier `h`

output is either `1`, if the sample satisfies the classifiers condition, or `-1`
"""
function (h::WeakClassifier)(x::Vector{<:Real})
    value = x[h.feature]
    if h.bigger
        return value >= h.threshold ? 1 : -1
    else
        return value < h.threshold ? 1 : -1
    end
end

"""
    Strong classifier `H`

structure representing adaboost classifier that selects class of a given sample using mulititude of weighted weak classifiers

# Fields
- `weights` - determine how important is correct classification of a given sample; initially set to 1/N
- `weaks` - set of weak classifiers used to making decisions
- `Z` - upper bound on classification error in individual steps of learning
- `samples` - number of data samples used for learning
"""
struct StrongClassifier
    weights::Vector{<:Real}
    weaks::Vector{WeakClassifier}
    Z::Vector{<:Real}
    samples::Integer

    """
        StrongClassifier(samples::Integer)

        constructs empty strong classifier with uniform initial data weights
    """
    function StrongClassifier(samples::Integer)
        initial = 1 / samples
        return new(fill(initial, samples), Vector{WeakClassifier}(undef, 0), zeros(0), samples)
    end
end

function Base.show(io::IO, H::StrongClassifier)
    println(io, "=== Strong classifier H ===")
    for i in 1:length(H.weaks)
        println(io, "$(H.weaks[i]) ... Z = $(H.Z[i])")
    end
end

"""
    (H::StrongClassifier)(x::Vector{<:Real})

classify given data sample `x` using strong classifier `H`

output is either `1` or `-1` based on weighted decisions of all weak classifiers
"""
function (H::StrongClassifier)(x::Vector{<:Real})
    return Integer(sign(sum([h.weight * h(x) for h in H.weaks])))
end

"""
    updateWeights!(H::StrongClassifier, h::WeakClassifier, X::Matrix{<:Real}, y::Vector{<:Integer})

update data weights based on correct labels `y` and those assigned by new weak classifier `h`
- increase weight for misclassified samples and decrease for correctly classified ones
- recompute error upper bound Z
"""
function updateWeights!(H::StrongClassifier, h::WeakClassifier, X::Matrix{<:Real}, y::Vector{<:Integer})
    newweights = [H.weights[i] * exp(-h.weight * y[i] * h(X[:, i])) for i in 1:size(X, 2)]
    Z = sum(newweights)
    push!(H.Z, Z * (length(H.Z) > 0 ? H.Z[end] : 1))
    H.weights .= newweights / Z
end

"""
    bestWeak(H::StrongClassifier, X::Matrix{<:Real}, Xsorted::Matrix{<:Real}, y::Vector{<:Integer})

select a weak classifier that misclassifies the least amount of samples in current state of data weights

`Xsorted` - matrix X with presorted rows to try different `thresholds`

returns optimal weak classifier and its classification error
"""
function bestWeak(H::StrongClassifier, X::Matrix{<:Real}, Xsorted::Matrix{<:Real}, y::Vector{<:Integer})
    dim, n = size(X)
    opterror = Inf64
    local opth

    for feature in 2:dim
        x = Xsorted[feature, :]
        for threshold in x
            improved = false
            h = WeakClassifier(feature, threshold, true)
            comp = [h(X[:, i]) for i in 1:n]
            diff = comp .!= y
            err = dot(diff, H.weights)
            if err < opterror
                opterror = err
                opth = h
                improved = true
            end

            h = WeakClassifier(feature, threshold, false)
            comp = [h(X[:, i]) for i in 1:n]
            diff = comp .!= y
            err = dot(diff, H.weights)
            if err < opterror
                opterror = err
                opth = h
                improved = true
            end

            if !improved
                break
            end
        end
    end

    return WeakClassifier(opth.feature, opth.threshold, opth.bigger, opterror), opterror
end

"""
    boost(data::Dataset, labels::Labels; limit::Integer = 100)

select at most `limit` weak classifiers for given dataset `data` and its correct `labels`
"""
function boost(data::Dataset, labels::Labels; limit::Integer = 100)
    X = data()
    y = labels()

    sorted = vcat([sort(row) for row in eachrow(X)]'...)

    n = data.passengers

    H = StrongClassifier(n)

    for _ in 1:limit
        h, error = bestWeak(H, X, sorted, y)

        if error >= 0.5
            break
        end

        push!(H.weaks, h)

        updateWeights!(H, h, X, y)

    end

    return H
end

"""
    classify(data::Dataset, H::StrongClassifier)

classify samples in `data` dataset with strong classifier `H`
"""
function classify(data::Dataset, H::StrongClassifier)
    return Labels(data, [H(data()[:, i]) for i in 1:data.passengers])
end


text:
###############################################################################
#
# Save figures as GIF
#
###############################################################################
"""
    save_gif!(figs::Array, fps::Int, gif::String)
    save_gif!(imgs::Array{Array,1}, fps::Int, fn::String)

Save an array of figures as GIF, given
- `figs` An array of figures
- `fps` Frame per second
- `fn` File name of the GIF to save
- `imgs` An array of array (loaded figures)
"""
function save_gif!(figs::Array{Figure,1}, fps::Int, fn::String)
    # load images to memory
    images = Array[];
    for fig in figs
        fig.savefig("temp_gif.png");
        img = load("temp_gif.png");
        push!(images, img);
    end
    rm("temp_gif.png");

    save_gif!(images, fps, fn);

    return nothing
end

function save_gif!(imgs::Array{Array,1}, fps::Int, fn::String)
    gif = cat(imgs...; dims=3);
    save(fn, gif; fps=fps);

    return nothing
end


text:
#=
These functions are used to display a ROOTFile is a tree-like fashion
by using `AbstractTrees` printing functions. We customize what the children
of ROOTFile and a TTree is, and how to print the final `node`.
=#
struct TKeyNode
    name::AbstractString
    classname::AbstractString
end
function children(f::T) where T <: Union{ROOTFile,ROOTDirectory}
    # display TTrees recursively
    # subsequent TTrees with duplicate fName will be skipped
    # since TKey cycle number is guaranteed to be decreasing
    # then all TKeys in the file which are not for a TTree
    seen = Set{String}()
    ch = Vector{Union{TTree,TKeyNode,ROOTDirectory}}()
    for k in keys(f)
        try
            obj = f[k]
            obj isa TTree || continue
            obj.fName ∈ seen && continue
            push!(ch, obj)
            push!(seen, obj.fName)
        catch
        end
    end
    tkeys = T === ROOTFile ? f.directory.keys : f.keys
    for tkey in tkeys
        kn = TKeyNode(tkey.fName, tkey.fClassName)
        kn.classname == "TTree" && continue
        if kn.classname == "TDirectory"
            push!(ch, f[tkey.fName])
        else
            push!(ch, kn)
        end
    end
    ch
end
function children(t::TTree)
    ks = keys(t)
    if length(ks) < 2
        return ks
    elseif length(ks) > 7
        return vcat(first(ks, 3), "⋮", ks[end-2:end])
    else
        return ks
    end
end
printnode(io::IO, t::TTree) = print(io, "$(t.fName) (TTree)")
printnode(io::IO, f::ROOTFile) = print(io, f.filename)
printnode(io::IO, f::ROOTDirectory) = print(io, "$(f.name) (TDirectory)")
printnode(io::IO, k::TKeyNode) = print(io, "$(k.name) ($(k.classname))")

Base.show(tree::LazyTree; kwargs...) = _show(stdout, tree; crop=:both, kwargs...)
Base.show(io::IO, tree::LazyTree; kwargs...) = _show(io, tree; kwargs...)
Base.show(io::IO, ::MIME"text/plain", tree::LazyTree) = _show(io, tree)
function _show(io::IO, tree::LazyTree; kwargs...)
    _hs = _make_header(tree)
    _ds = displaysize(io)
    PrettyTables.pretty_table(
        io,
        innertable(tree);
        header=_hs,
        alignment=:l,
        vlines=[1],
        hlines=[:header],
        crop_num_lines_at_beginning=2,
        row_number_alignment=:l,
        row_number_column_title="Row",
        show_row_number=true,
        compact_printing=false,
        formatters=(v, i, j) -> _treeformat(v, _ds[2] ÷ min(8, length(_hs[1]))),
        display_size=(min(_ds[1], 40), min(_ds[2], 160)),
        kwargs...
    )
    nothing
end

function Base.show(io::IO, ::MIME"text/plain", br::LazyBranch)
    print(io, summary(br))
    println(io, ": ")
    if length(br) < 200
        Base.print_array(IOContext(io, :limit => true), br[:])
    else
        head = @async br[1:100]
        tail = @async br[end-99:end]
        wait(head)
        wait(tail)
        Base.print_array(IOContext(io, :limit => true), Vcat(head.result, tail.result))
    end
    nothing
end

# stop crazy stracktrace
function Base.show(io::IO,
    ::Type{<:LazyTree{<:UnROOT.TypedTables.Table{NamedTuple{Ns, Vs}}}}) where {T, Ns, Vs}
    elip = length(Ns) > 5 ? "..." : ""
    println(io, "LazyTree with $(length(Ns)) branches:")
    println(io, join(first(Ns, 5), ", "), elip)
end

function Base.show(io::IO, ::MIME"text/html", tree::LazyTree)
    _hs = _make_header(tree)
    maxrows = 10
    maxcols = 30
    nrow = length(tree)
    t = @view innertable(tree)[1:min(maxrows,nrow)]
    ncol = length(Tables.columns(t))
    withcommas(value) = reverse(join(join.(Iterators.partition(reverse(string(value)),3)),","))
    write(io, "<p>")
    write(io, "$(withcommas(nrow)) rows × $(ncol) columns")
    if (nrow > maxrows) && (ncol > maxcols)
        write(io, " (omitted printing of $(withcommas(nrow-maxrows)) rows and $(ncol-maxcols) columns)")
    elseif (nrow > maxrows)
        write(io, " (omitted printing of $(withcommas(nrow-maxrows)) rows)")
    elseif (ncol > maxcols)
        write(io, " (omitted printing of $(ncol-maxcols) columns)")
    end
    write(io, "</p>")
    PrettyTables.pretty_table(
        io,
        t;
        header=_hs,
        alignment=:l,
        row_number_column_title="",
        show_row_number=true,
        compact_printing=false,
        filters_col     = ((_,i) -> i <= maxcols,),
        formatters=(v, i, j) -> _treeformat(v, 100),
        tf = PrettyTables.HTMLTableFormat(css = """th { color: #000; background-color: #fff; }"""),
        backend=Val(:html),
    )
    nothing
end
_symtup2str(symtup, trunc=15) = collect(first.(string.(symtup), trunc))
function _make_header(t)
    pn = propertynames(t)
    header = _symtup2str(pn)
    subheader = _symtup2str(Tables.columntype.(Ref(innertable(t)), pn))
    (header, subheader)
end
function _treeformat(val, trunc)
    s = if val isa AbstractArray{T} where T<:Integer
        T = Int
        replace(string(T.(val)), string(T)=>"")
    elseif val isa AbstractArray{T} where T<:AbstractFloat
        T = eltype(val)
        replace(string(round.(T.(val); sigdigits=3)), string(T)=>"")
    else
        string(val)
    end
    first(s, trunc)
end


text:
module Sesam

using ProgressMeter
using CSV
using DataFrames
using Dates
using JuMP
using Gurobi

println("Using Sesam module.")

# utility functions
include("util.jl")

# power plants
include("powerplants.jl")

# renewables
include("renewables.jl")

# storages
include("storages.jl")

# nodes
include("nodes.jl")

# lines
include("lines.jl")

# Powerflow helpers
include("powerflow.jl")

### models

# economic dispatch
include("economic_dispatch.jl")

# congestion management
include("congestion_management.jl")

# congestion management with PtG
include("congestion_management_ptg.jl")

# helper functions
include("helpers.jl")

export
    CSV,
    DataFrames,
    JuMP,
    Gurobi,
    ProgressMeter,
    ElectronDisplay,
    Nodes,
    Lines,
    linespmax,
    bprime,
    PowerPlants,
    Renewables,
    Storages,
    EconomicDispatch,
    EconomicDispatch80,
    CongestionManagement,
    CongestionManagementPHS80,
    CongestionManagementPtG,
    hoursInWeek,
    hoursInDay,
    hoursInHalfDay,
    has,
    getKeyVector,
    dayslicer

end


text:
using ReverseDiffSparse
using Compat, Compat.Test, Compat.LinearAlgebra
using MathProgBase

ex = :(sin(x[1]^2) + cos(x[2]*4)/5-2.0)

nd,const_values = expr_to_nodedata(ex)
adj = adjmat(nd)

#@show nd

storage = zeros(length(nd))
partials_storage = zeros(length(nd))
reverse_storage = zeros(length(nd))

x = [2.0,3.0]
#@show x
fval = forward_eval(storage,partials_storage,nd,adj,const_values,[],x,[])
true_val = sin(x[1]^2) + cos(x[2]*4)/5 -2.0
@test isapprox(fval,true_val)

grad = zeros(2)
reverse_eval(reverse_storage,partials_storage,nd,adj)
reverse_extract(grad,reverse_storage,nd,adj,[],1.0)

true_grad = [2*x[1]*cos(x[1]^2), -4*sin(x[2]*4)/5]
@test isapprox(grad,true_grad)

# subexpressions

nd_outer = [NodeData(SUBEXPRESSION,1,-1)]
li = list_subexpressions(nd_outer)
@test li == [1]
li,li_individual = order_subexpressions(Vector{NodeData}[nd_outer], Vector{NodeData}[nd])
@test li == [1]
@test li_individual[1] == [1]

nd_outer2 = [NodeData(CALL,operator_to_id[:+],-1),NodeData(SUBEXPRESSION,2,1),NodeData(SUBEXPRESSION,1,1)]
li = list_subexpressions(nd_outer2)
@test li == [1,2]
li,li_individual = order_subexpressions(Vector{NodeData}[nd_outer2], Vector{NodeData}[nd,nd_outer])
@test li == [1,2]
@test li_individual[1] == [1,2]

adj_outer = adjmat(nd_outer)
outer_storage = zeros(1)
outer_storage_partials = zeros(1)
fval = forward_eval(outer_storage,outer_storage_partials,nd_outer,adj_outer,[],[],x,[fval])
@test isapprox(fval,true_val)

outer_reverse_storage = zeros(1)
fill!(grad,0.0)
subexpr_output = zeros(1)
reverse_eval(outer_reverse_storage,outer_storage_partials,nd_outer,adj_outer)
reverse_extract(grad,outer_reverse_storage,nd_outer,adj_outer,subexpr_output,1.0)
@assert subexpr_output[1] == 1.0
reverse_eval(reverse_storage,partials_storage,nd,adj)
reverse_extract(grad,reverse_storage,nd,adj,[],subexpr_output[1])
@test isapprox(grad,true_grad)

ex = :((1/x[1])^x[2]-x[3])

nd,const_values = expr_to_nodedata(ex)
adj = adjmat(nd)

storage = zeros(length(nd))
partials_storage = zeros(length(nd))
reverse_storage = zeros(length(nd))

x = [2.5,3.5,1.0]
#@show x
fval = forward_eval(storage,partials_storage,nd,adj,const_values,[],x,[])
true_val = (1/x[1])^x[2]-x[3]
@test isapprox(fval,true_val)

grad = zeros(3)
reverse_eval(reverse_storage,partials_storage,nd,adj)
reverse_extract(grad,reverse_storage,nd,adj,[],1.0)

true_grad = [-x[2]*x[1]^(-x[2]-1), -((1/x[1])^x[2])*log(x[1]),-1]
@test isapprox(grad,true_grad)

# logical expressions
ex = :(x[1] > 0.5 && x[1] < 0.9)
nd,const_values = expr_to_nodedata(ex)
adj = adjmat(nd)
storage = zeros(length(nd))
partials_storage = zeros(length(nd))
fval = forward_eval(storage,partials_storage,nd,adj,const_values,[],[1.5],[])
@test fval == 0
fval = forward_eval(storage,partials_storage,nd,adj,const_values,[],[0.6],[])
@test fval == 1

ex = :(ifelse(x[1] >= 0.5 || x[1] <= 0.1,x[1],5))
nd,const_values = expr_to_nodedata(ex)
adj = adjmat(nd)
storage = zeros(length(nd))
partials_storage = zeros(length(nd))
reverse_storage = zeros(length(nd))
grad = zeros(1)
fval = forward_eval(storage,partials_storage,nd,adj,const_values,[],[1.5],[])
@test fval == 1.5
reverse_eval(reverse_storage,partials_storage,nd,adj)
reverse_extract(grad,reverse_storage,nd,adj,[],1.0)
@test grad[1] == 1

fval = forward_eval(storage,partials_storage,nd,adj,const_values,[],[-0.1],[])
@test fval == -0.1
fill!(grad,0)
reverse_eval(reverse_storage,partials_storage,nd,adj)
reverse_extract(grad,reverse_storage,nd,adj,[],1.0)
@test grad[1] == 1

fval = forward_eval(storage,partials_storage,nd,adj,const_values,[],[0.2],[])
@test fval == 5
fill!(grad,0)
reverse_eval(reverse_storage,partials_storage,nd,adj)
reverse_extract(grad,reverse_storage,nd,adj,[],1.0)
@test grad[1] == 0

ex = :(sqrt(x[1]))
nd,const_values = expr_to_nodedata(ex)
adj = adjmat(nd)
storage = zeros(length(nd))
partials_storage = zeros(length(nd))
reverse_storage = zeros(length(nd))
grad = zeros(1)
fval = forward_eval(storage,partials_storage,nd,adj,const_values,[],[-1.5],[])
@test isnan(fval)
reverse_eval(reverse_storage,partials_storage,nd,adj)
reverse_extract(grad,reverse_storage,nd,adj,[],1.0)
@test isnan(grad[1])

# parameters
nd = [NodeData(PARAMETER,1,-1)]
adj = adjmat(nd)
storage = zeros(length(nd))
partials_storage = zeros(length(nd))
fval = forward_eval(storage,partials_storage,nd,adj,[],[105.2],[-0.1],[])
@test fval == 105.2

# abs
ex = :(abs(x[1]))

nd,const_values = expr_to_nodedata(ex)
adj = adjmat(nd)

storage = zeros(length(nd))
partials_storage = zeros(length(nd))
reverse_storage = zeros(length(nd))
fval = forward_eval(storage,partials_storage,nd,adj,const_values,[],[2.0],[])
@test fval == 2.0
grad = zeros(1)
reverse_eval(reverse_storage,partials_storage,nd,adj)
reverse_extract(grad,reverse_storage,nd,adj,[],1.0)
@test grad[1] == 1.0

fval = forward_eval(storage,partials_storage,nd,adj,const_values,[],[-2.0],[])
@test fval == 2.0
grad = zeros(1)
reverse_eval(reverse_storage,partials_storage,nd,adj)
reverse_extract(grad,reverse_storage,nd,adj,[],1.0)
@test grad[1] == -1.0

# https://github.com/JuliaOpt/JuMP.jl/issues/855
ex = :(ifelse(x[1]<=3.0, (x[1]-2.0)^2, 2*log(x[1]-2.0)+1.0))

nd,const_values = expr_to_nodedata(ex)
adj = adjmat(nd)

storage = zeros(length(nd))
partials_storage = zeros(length(nd))
reverse_storage = zeros(length(nd))
fval = forward_eval(storage,partials_storage,nd,adj,const_values,[],[-1.0],[])
@test fval == 9.0
grad = zeros(1)
reverse_eval(reverse_storage,partials_storage,nd,adj)
reverse_extract(grad,reverse_storage,nd,adj,[],1.0)
@test grad[1] == -6.0

fval = forward_eval(storage,partials_storage,nd,adj,const_values,[],[2.0],[])
@test fval == 0.0
grad = zeros(1)
reverse_eval(reverse_storage,partials_storage,nd,adj)
reverse_extract(grad,reverse_storage,nd,adj,[],1.0)
@test grad[1] == 0.0

function test_linearity(ex,testval,IJ = [],indices=[])
    nd,const_values = expr_to_nodedata(ex)
    adj = adjmat(nd)
    linearity = classify_linearity(nd,adj,[])
    @test linearity[1] == testval
    idxset = Coloring.IndexedSet(100)
    edgelist = compute_hessian_sparsity(nd,adj,linearity,idxset,Array{Set{Tuple{Int,Int}}}(undef, 0), Array{Vector{Int}}(undef, 0))
    if linearity[1] != NONLINEAR
        @test length(edgelist) == 0
    elseif length(IJ) > 0
        @test IJ == edgelist
    end
    if length(indices) > 0
        ix = sort(collect(compute_gradient_sparsity(nd)))
        compute_gradient_sparsity!(idxset, nd)
        ix2 = sort(collect(idxset))
        empty!(idxset)
        @test ix == indices
        @test ix == ix2
    end
end

test_linearity(:(sin(x[1]^2) + cos(x[2]*4)-2.0), NONLINEAR, Set([(2,2),(1,1)]), [1,2])
test_linearity(:(3*4*(x[1]+x[2])), LINEAR)
test_linearity(:(x[3]*x[2]), NONLINEAR, Set([(3,2),(3,3),(2,2)]),[2,3])
test_linearity(:(3+4), CONSTANT)
test_linearity(:(sin(3)+x[1]), LINEAR)
test_linearity(:(cos(x[3])*sin(3)+x[1]), NONLINEAR, Set([(3,3)]),[1,3])
test_linearity(:(x[1]-3x[2]), LINEAR)
test_linearity(:(-x[1]), LINEAR)
test_linearity(:(+x[1]), LINEAR)
test_linearity(:(x[1]^x[2]), NONLINEAR, Set([(2,2),(1,1),(2,1)]))
test_linearity(:(x[1]/3+x[2]), LINEAR)
test_linearity(:(3/(x[1]*x[2])), NONLINEAR, Set([(2,2),(1,1),(2,1)]))
test_linearity(:(1/(x[1]+3)), NONLINEAR)
test_linearity(:(ifelse(x[1] <= 1,x[1],x[2])), PIECEWISE_LINEAR, Set([]))
test_linearity(:(ifelse(x[1] <= 1,x[1]^2,x[2])), NONLINEAR, Set([(1,1)]))
test_linearity(:(ifelse(1 <= 1,2,3)), CONSTANT)
test_linearity(:(1/ifelse(x[1] < 1, x[1],0)), NONLINEAR, Set([(1,1)]))

# eliminating fixed variables and constants
ex = :(sin(x[1]^2) + cos(x[2]*(2*2))/5-2.0)
nd,const_values = expr_to_nodedata(ex)
adj = adjmat(nd)
storage = zeros(length(nd))
partials_storage = zeros(length(nd))
x = [2.0,3.0]
fval = forward_eval(storage,partials_storage,nd,adj,const_values,[],x,[])
linearity = classify_linearity(nd,adj,[],[false,true])
new_nd = simplify_constants(storage,nd,adj,const_values,linearity)
@test length(new_nd) < length(nd)
new_adj = adjmat(new_nd)
true_val = sin(x[1]^2) + cos(x[2]*4)/5 -2.0
x[2] = 100 # this shouldn't affect the answer because we said x[2] is fixed
fval = forward_eval(storage,partials_storage,new_nd,new_adj,const_values,[],x,[])
@test isapprox(fval,true_val)
# all variables fixed
linearity = classify_linearity(nd,adj,[],[true,true])
new_nd = simplify_constants(storage,nd,adj,const_values,linearity)
@test length(new_nd) == 1

# user-defined functions
#Φ(x,y) = 1/3(y)^3 - 2x^2
# c(x) = cos(x)
struct ΦEvaluator <: MathProgBase.AbstractNLPEvaluator
end
function MathProgBase.eval_f(::ΦEvaluator,x)
    @assert length(x) == 2
    return (1/3)*x[2]^3-2x[1]^2
end
function MathProgBase.eval_grad_f(::ΦEvaluator,grad,x)
    grad[1] = -4x[1]
    grad[2] = x[2]^2
end
r = ReverseDiffSparse.UserOperatorRegistry()
register_multivariate_operator!(r,:Φ,ΦEvaluator())
register_univariate_operator!(r,:c,cos,x->-sin(x),x->-cos(x))
Φ(x,y) = MathProgBase.eval_f(ΦEvaluator(),[x,y])
ex = :(Φ(x[2],x[1]-1)*c(x[3]))
nd,const_values = expr_to_nodedata(ex,r)
@test ReverseDiffSparse.has_user_multivariate_operators(nd)
adj = adjmat(nd)
storage = zeros(length(nd))
partials_storage = zeros(length(nd))
reverse_storage = zeros(length(nd))
x = [2.0,3.0,4.0]
fval = forward_eval(storage,partials_storage,nd,adj,const_values,[],x,[],zeros(2),zeros(2),user_operators=r)
true_val = Φ(x[2],x[1]-1)*cos(x[3])
@test isapprox(fval,true_val)
grad = zeros(3)
reverse_eval(reverse_storage,partials_storage,nd,adj)
reverse_extract(grad,reverse_storage,nd,adj,[],1.0)

true_grad = [cos(x[3])*(x[1]-1)^2, -4cos(x[3])*x[2], -sin(x[3])*Φ(x[2],x[1]-1)]
@test isapprox(grad,true_grad)

using DualNumbers
using ForwardDiff

# dual forward test
function dualforward(ex, x; ignore_nan=false)
    nd,const_values = expr_to_nodedata(ex)
    adj = adjmat(nd)
    forward_storage = zeros(length(nd))
    partials_storage = zeros(length(nd))
    reverse_storage = zeros(length(nd))

    fval = forward_eval(forward_storage,partials_storage,nd,adj,const_values,[],x,[])
    grad = zeros(length(x))
    reverse_eval(reverse_storage,partials_storage,nd,adj)
    reverse_extract(grad,reverse_storage,nd,adj,[],1.0)

    zero_ϵ = zero(ForwardDiff.Partials{1,Float64})
    forward_storage_ϵ = fill(zero_ϵ,length(nd))
    partials_storage_ϵ = fill(zero_ϵ,length(nd))
    x_values_ϵ = fill(ForwardDiff.Partials((1.0,)),length(x))
    reverse_storage_ϵ = fill(zero_ϵ,length(nd))
    output_ϵ = fill(zero_ϵ,length(x))
    fval_ϵ = forward_eval_ϵ(forward_storage,forward_storage_ϵ,partials_storage,partials_storage_ϵ,nd,adj,x_values_ϵ,[])
    reverse_eval_ϵ(output_ϵ,reverse_storage,reverse_storage_ϵ,partials_storage,partials_storage_ϵ,nd,adj,[],[],2.0,zero_ϵ)
    @test isapprox(fval_ϵ[1], dot(grad,ones(length(x))))

    # compare with running dual numbers
    forward_dual_storage = zeros(DualNumbers.Dual{Float64},length(nd))
    partials_dual_storage = zeros(DualNumbers.Dual{Float64},length(nd))
    output_dual_storage = zeros(DualNumbers.Dual{Float64},length(x))
    reverse_dual_storage = zeros(DualNumbers.Dual{Float64},length(nd))
    x_dual = [DualNumbers.Dual(x[i],1.0) for i in 1:length(x)]
    fval = forward_eval(forward_dual_storage,partials_dual_storage,nd,adj,const_values,[],x_dual,[])
    reverse_eval(reverse_dual_storage,partials_dual_storage,nd,adj)
    reverse_extract(output_dual_storage,reverse_dual_storage,nd,adj,[],DualNumbers.Dual(2.0))
    for k in 1:length(nd)
        @test isapprox(epsilon(forward_dual_storage[k]), forward_storage_ϵ[k][1])
        if !(isnan(epsilon(partials_dual_storage[k])) && ignore_nan)
            @test isapprox(epsilon(partials_dual_storage[k]), partials_storage_ϵ[k][1])
        else
            @test !isnan(forward_storage_ϵ[k][1])
        end
        if !(isnan(epsilon(reverse_dual_storage[k])) && ignore_nan)
            @test isapprox(epsilon(reverse_dual_storage[k]), reverse_storage_ϵ[k][1]/2)
        else
            @test !isnan(reverse_storage_ϵ[k][1])
        end
    end
    for k in 1:length(x)
        if !(isnan(epsilon(output_dual_storage[k])) && ignore_nan)
            @test isapprox(epsilon(output_dual_storage[k]), output_ϵ[k][1])
        else
            @test !isnan(output_ϵ[k][1])
        end
    end
end

dualforward(:(sin(x[1]^2) + cos(x[2]*4)/5-2.0),[1.0,2.0])
dualforward(:(sin(x[1]^x[2]) + cos(x[2]*4)/5-2.0),[1.0,2.0])
dualforward(:(sin(x[1]^3) + cos(x[1]*x[2]*4)/5-2.0),[1.0,0.0])
dualforward(:(x[1]*x[2]),[3.427139283036299e-206,1.0], ignore_nan=true)

include("test_coloring.jl")


text:
import CollectiveSpins

# Band structure (Omega) and collective decay rate (Gamma) for 1D chain for x-axis (along chain)
# and in the yz-plane (plane perpendicular to the atom chain) polarized atoms.
# See Asenjo-Garcia et al 10.1103/PhysRevX.7.031024 Figs. 1b and 1c.

# Define chain parameters
a = 0.2   # spin-spin distance
polarization_par = [1, 0, 0]
polarization_perp = [0, 1, 1]

k_max = 10000
k_list = [iii*pi/(a*k_max) for iii=-k_max:k_max]

band_structure_par = []
band_structure_perp = []
decay_par = []
decay_perp = []
for k in k_list
    append!(band_structure_par, CollectiveSpins.Omega_k_chain(k, a, polarization_par))
    append!(band_structure_perp, CollectiveSpins.Omega_k_chain(k, a, polarization_perp))
    append!(decay_par, CollectiveSpins.Gamma_k_chain(k, a, polarization_par))
    append!(decay_perp, CollectiveSpins.Gamma_k_chain(k, a, polarization_perp))
end

using PyPlot

subplot(211)
plot(k_list*a/pi, band_structure_par)
plot(k_list*a/pi, band_structure_perp, color="red")
xlabel("k d/pi")
ylabel("Omega(k)/gamma0")

subplot(212)
plot(k_list*a/pi, decay_par)
plot(k_list*a/pi, decay_perp, color="red")
xlabel("k d/pi")
ylabel("Gamma(k)/gamma0")

tight_layout()
savefig("band_structure_decay.svg")


text:
"""
    Diagnostic(first_byte, last_byte; [error="msg" | warning="msg"])

A diagnostic message, referring to the source code byte range
first_byte:last_byte, with a `warning` or `error` message.

Messages should be concise, matter-of-fact and not include decorations:

* Concise: "Show don't tell". Where possible, let's show the user what's wrong
  by annotating their original source code via the byte range.
* Matter-of-fact: Admonishing the user isn't helpful. Let's gently show them
  what's wrong instead, using a neutral tone.
* Decorations: Capitalization, punctuation and diagnostic class ("error" /
  "warning") should be omitted. These decorations will be added by the
  formatting code.

TODO: At some point we should enhance Diagnostic to allow multiple sub-ranges
for better annotation. Let's follow the excellent precedent set by Rust's
[rustc_errors::Diagnostic](https://doc.rust-lang.org/stable/nightly-rustc/rustc_errors/struct.Diagnostic.html).

TODO: We should cater for extended descriptions containing multiple sentences
via a diagnostic code which can be used to look up detailed information. Again,
Rust does this well.
"""
struct Diagnostic
    first_byte::Int
    last_byte::Int
    level::Symbol
    message::String
end

function Diagnostic(first_byte, last_byte; error=nothing, warning=nothing)
    message = !isnothing(error)   ? error :
              !isnothing(warning) ? warning :
              error("No message in diagnostic")
    level = !isnothing(error) ? :error : :warning
    Diagnostic(first_byte, last_byte, level, message)
end

first_byte(d::Diagnostic) = d.first_byte
last_byte(d::Diagnostic)  = d.last_byte
is_error(d::Diagnostic)   = d.level == :error

function show_diagnostic(io::IO, diagnostic::Diagnostic, source::SourceFile)
    color,prefix = diagnostic.level == :error   ? (:light_red, "Error")      :
                   diagnostic.level == :warning ? (:light_yellow, "Warning") :
                   diagnostic.level == :note    ? (:light_blue, "Note")      :
                   (:normal, "Info")
    line, col = source_location(source, first_byte(diagnostic))
    linecol = "$line:$col"
    if !isnothing(source.filename)
        locstr = "$(source.filename):$linecol"
        if get(io, :color, false)
            # Also add hyperlinks in color terminals
            url = "file://$(abspath(source.filename))#$linecol"
            locstr = "\e]8;;$url\e\\$locstr\e]8;;\e\\"
        end
    else
        locstr = "line $linecol"
    end
    print(io, prefix, ": ")
    printstyled(io, diagnostic.message, color=color)
    printstyled(io, "\n", "@ $locstr", color=:light_black)
    print(io, "\n")

    p = first_byte(diagnostic)
    q = last_byte(diagnostic)
    text = sourcetext(source)
    if q < p || (p == q && source[p] == '\n')
        # An empty or invisible range!  We expand it symmetrically to make it
        # visible.
        p = max(firstindex(text), prevind(text, p))
        q = min(lastindex(text), nextind(text, q))
    end

    # p and q mark the start and end of the diagnostic range. For context,
    # buffer these out to the surrouding lines.
    a,b = source_line_range(source, p, context_lines_before=2, context_lines_after=1)
    c,d = source_line_range(source, q, context_lines_before=1, context_lines_after=2)

    hicol = (100,40,40)

    # TODO: show line numbers on left

    print(io, source[a:prevind(text, p)])
    # There's two situations, either
    if b >= c
        # The diagnostic range is compact and we show the whole thing
        # a...............
        # .....p...q......
        # ...............b
        _printstyled(io, source[p:q]; color=hicol)
    else
        # Or large and we trucate the code to show only the region around the
        # start and end of the error.
        # a...............
        # .....p..........
        # ...............b
        # (snip)
        # c...............
        # .....q..........
        # ...............d
        _printstyled(io, source[p:b]; color=hicol)
        println(io, "…")
        _printstyled(io, source[c:q]; color=hicol)
    end
    print(io, source[nextind(text,q):d])
    println(io)
end

function show_diagnostics(io::IO, diagnostics::AbstractVector{Diagnostic}, source::SourceFile)
    for d in diagnostics
        show_diagnostic(io, d, source)
    end
end

function show_diagnostics(io::IO, diagnostics::AbstractVector{Diagnostic}, text::AbstractString)
    if !isempty(diagnostics)
        show_diagnostics(io, diagnostics, SourceFile(text))
    end
end

function any_error(diagnostics::AbstractVector{Diagnostic})
    any(is_error(d) for d in diagnostics)
end


text:
using Documenter, RegistryActivity

makedocs(
    modules = [RegistryActivity],
    sitename = "RegistryActivity.jl",
)

deploydocs(
    repo = "github.com/JuliaEcosystem/RegistryActivity.jl.git",
    devbranch = "main",
)


text:
module DifferentiableStateSpaceModels

using Logging
using MatrixEquations
using TensorCast
using DistributionsAD
using ChainRulesCore
using DocStringExtensions
using Distributions
using Parameters
using MacroTools
using LinearAlgebra
using Zygote
using NLsolve
using Symbolics
using SymbolicUtils
using Latexify
using LaTeXStrings
using StructArrays
using RecursiveFactorization # Or perhaps MKL will dominate this on all hardware.

# Flip off when not debugging
using Infiltrator

export make_perturbation_model, default_model_cache_location, PerturbationModel, SolverCache, PerturbationSolverSettings, @include_example_module, @make_and_include_perturbation_model, generate_perturbation, generate_perturbation_derivatives!

export solve # will be replaced by SciML

include("utils.jl")
include("symbolic_utils.jl")
include("make_perturbation_model.jl")
include("types.jl")
include("generate_perturbation.jl")
include("generate_perturbation_derivatives.jl")
include("sequence.jl")
include("Examples/Examples.jl")

end # module