text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct RoundRotorMachine <: Machine
        R::Float64
        Td0_p::Float64
        Td0_pp::Float64
        Tq0_p::Float64
        Tq0_pp::Float64
        Xd::Float64
        Xq::Float64
        Xd_p::Float64
        Xq_p::Float64
        Xd_pp::Float64
        Xl::Float64
        Se::Tuple{Float64, Float64}
        ext::Dict{String, Any}
        states::Vector{Symbol}
        n_states::Int64
        internal::InfrastructureSystemsInternal
    end

Parameters of 4-states round-rotor synchronous machine with quadratic/exponential saturation:
IEEE Std 1110 §5.3.2 (Model 2.2). GENROU or GENROE model in PSSE and PSLF.

# Arguments
- `R::Float64`: Armature resistance, validation range: (0, nothing)
- `Td0_p::Float64`: Time constant of transient d-axis voltage, validation range: (0, nothing)
- `Td0_pp::Float64`: Time constant of sub-transient d-axis voltage, validation range: (0, nothing)
- `Tq0_p::Float64`: Time constant of transient q-axis voltage, validation range: (0, nothing)
- `Tq0_pp::Float64`: Time constant of sub-transient q-axis voltage, validation range: (0, nothing)
- `Xd::Float64`: Reactance after EMF in d-axis per unit, validation range: (0, nothing)
- `Xq::Float64`: Reactance after EMF in q-axis per unit, validation range: (0, nothing)
- `Xd_p::Float64`: Transient reactance after EMF in d-axis per unit, validation range: (0, nothing)
- `Xq_p::Float64`: Transient reactance after EMF in q-axis per unit, validation range: (0, nothing)
- `Xd_pp::Float64`: Sub-Transient reactance after EMF in d-axis per unit. Note: Xd_pp = Xq_pp, validation range: (0, nothing)
- `Xl::Float64`: Stator leakage reactance, validation range: (0, nothing)
- `Se::Tuple{Float64, Float64}`: Saturation factor at 1 and 1.2 pu flux: S(1.0) = B(|ψ_pp|-A)^2
- `ext::Dict{String, Any}`
- `states::Vector{Symbol}`: The states are:
	eq_p: q-axis generator voltage behind the transient reactance,
	ed_p: d-axis generator voltage behind the transient reactance,
	ψ_kd: flux linkage in the first equivalent damping circuit in the d-axis,
	ψ_kq: flux linkage in the first equivalent damping circuit in the d-axis
- `n_states::Int64`
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct RoundRotorMachine <: Machine
    "Armature resistance"
    R::Float64
    "Time constant of transient d-axis voltage"
    Td0_p::Float64
    "Time constant of sub-transient d-axis voltage"
    Td0_pp::Float64
    "Time constant of transient q-axis voltage"
    Tq0_p::Float64
    "Time constant of sub-transient q-axis voltage"
    Tq0_pp::Float64
    "Reactance after EMF in d-axis per unit"
    Xd::Float64
    "Reactance after EMF in q-axis per unit"
    Xq::Float64
    "Transient reactance after EMF in d-axis per unit"
    Xd_p::Float64
    "Transient reactance after EMF in q-axis per unit"
    Xq_p::Float64
    "Sub-Transient reactance after EMF in d-axis per unit. Note: Xd_pp = Xq_pp"
    Xd_pp::Float64
    "Stator leakage reactance"
    Xl::Float64
    "Saturation factor at 1 and 1.2 pu flux: S(1.0) = B(|ψ_pp|-A)^2"
    Se::Tuple{Float64, Float64}
    ext::Dict{String, Any}
    "The states are:
	eq_p: q-axis generator voltage behind the transient reactance,
	ed_p: d-axis generator voltage behind the transient reactance,
	ψ_kd: flux linkage in the first equivalent damping circuit in the d-axis,
	ψ_kq: flux linkage in the first equivalent damping circuit in the d-axis"
    states::Vector{Symbol}
    n_states::Int64
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function RoundRotorMachine(R, Td0_p, Td0_pp, Tq0_p, Tq0_pp, Xd, Xq, Xd_p, Xq_p, Xd_pp, Xl, Se, ext=Dict{String, Any}(), )
    RoundRotorMachine(R, Td0_p, Td0_pp, Tq0_p, Tq0_pp, Xd, Xq, Xd_p, Xq_p, Xd_pp, Xl, Se, ext, [:eq_p, :ed_p, :ψ_kd, :ψ_kq], 4, InfrastructureSystemsInternal(), )
end

function RoundRotorMachine(; R, Td0_p, Td0_pp, Tq0_p, Tq0_pp, Xd, Xq, Xd_p, Xq_p, Xd_pp, Xl, Se, ext=Dict{String, Any}(), )
    RoundRotorMachine(R, Td0_p, Td0_pp, Tq0_p, Tq0_pp, Xd, Xq, Xd_p, Xq_p, Xd_pp, Xl, Se, ext, )
end

# Constructor for demo purposes; non-functional.
function RoundRotorMachine(::Nothing)
    RoundRotorMachine(;
        R=0,
        Td0_p=0,
        Td0_pp=0,
        Tq0_p=0,
        Tq0_pp=0,
        Xd=0,
        Xq=0,
        Xd_p=0,
        Xq_p=0,
        Xd_pp=0,
        Xl=0,
        Se=(0.0, 0.0),
        ext=Dict{String, Any}(),
    )
end

"""Get RoundRotorMachine R."""
get_R(value::RoundRotorMachine) = value.R
"""Get RoundRotorMachine Td0_p."""
get_Td0_p(value::RoundRotorMachine) = value.Td0_p
"""Get RoundRotorMachine Td0_pp."""
get_Td0_pp(value::RoundRotorMachine) = value.Td0_pp
"""Get RoundRotorMachine Tq0_p."""
get_Tq0_p(value::RoundRotorMachine) = value.Tq0_p
"""Get RoundRotorMachine Tq0_pp."""
get_Tq0_pp(value::RoundRotorMachine) = value.Tq0_pp
"""Get RoundRotorMachine Xd."""
get_Xd(value::RoundRotorMachine) = value.Xd
"""Get RoundRotorMachine Xq."""
get_Xq(value::RoundRotorMachine) = value.Xq
"""Get RoundRotorMachine Xd_p."""
get_Xd_p(value::RoundRotorMachine) = value.Xd_p
"""Get RoundRotorMachine Xq_p."""
get_Xq_p(value::RoundRotorMachine) = value.Xq_p
"""Get RoundRotorMachine Xd_pp."""
get_Xd_pp(value::RoundRotorMachine) = value.Xd_pp
"""Get RoundRotorMachine Xl."""
get_Xl(value::RoundRotorMachine) = value.Xl
"""Get RoundRotorMachine Se."""
get_Se(value::RoundRotorMachine) = value.Se
"""Get RoundRotorMachine ext."""
get_ext(value::RoundRotorMachine) = value.ext
"""Get RoundRotorMachine states."""
get_states(value::RoundRotorMachine) = value.states
"""Get RoundRotorMachine n_states."""
get_n_states(value::RoundRotorMachine) = value.n_states
"""Get RoundRotorMachine internal."""
get_internal(value::RoundRotorMachine) = value.internal

"""Set RoundRotorMachine R."""
set_R!(value::RoundRotorMachine, val::Float64) = value.R = val
"""Set RoundRotorMachine Td0_p."""
set_Td0_p!(value::RoundRotorMachine, val::Float64) = value.Td0_p = val
"""Set RoundRotorMachine Td0_pp."""
set_Td0_pp!(value::RoundRotorMachine, val::Float64) = value.Td0_pp = val
"""Set RoundRotorMachine Tq0_p."""
set_Tq0_p!(value::RoundRotorMachine, val::Float64) = value.Tq0_p = val
"""Set RoundRotorMachine Tq0_pp."""
set_Tq0_pp!(value::RoundRotorMachine, val::Float64) = value.Tq0_pp = val
"""Set RoundRotorMachine Xd."""
set_Xd!(value::RoundRotorMachine, val::Float64) = value.Xd = val
"""Set RoundRotorMachine Xq."""
set_Xq!(value::RoundRotorMachine, val::Float64) = value.Xq = val
"""Set RoundRotorMachine Xd_p."""
set_Xd_p!(value::RoundRotorMachine, val::Float64) = value.Xd_p = val
"""Set RoundRotorMachine Xq_p."""
set_Xq_p!(value::RoundRotorMachine, val::Float64) = value.Xq_p = val
"""Set RoundRotorMachine Xd_pp."""
set_Xd_pp!(value::RoundRotorMachine, val::Float64) = value.Xd_pp = val
"""Set RoundRotorMachine Xl."""
set_Xl!(value::RoundRotorMachine, val::Float64) = value.Xl = val
"""Set RoundRotorMachine Se."""
set_Se!(value::RoundRotorMachine, val::Tuple{Float64, Float64}) = value.Se = val
"""Set RoundRotorMachine ext."""
set_ext!(value::RoundRotorMachine, val::Dict{String, Any}) = value.ext = val
"""Set RoundRotorMachine states."""
set_states!(value::RoundRotorMachine, val::Vector{Symbol}) = value.states = val
"""Set RoundRotorMachine n_states."""
set_n_states!(value::RoundRotorMachine, val::Int64) = value.n_states = val
"""Set RoundRotorMachine internal."""
set_internal!(value::RoundRotorMachine, val::InfrastructureSystemsInternal) = value.internal = val


text:
isdefined(Base, :__precompile__) && __precompile__()

module InfrastructureSystems

import CSV
import DataFrames
import Dates
import JSON
import JSON2
import PrettyTables
import TimeSeries

using DocStringExtensions

@template (FUNCTIONS, METHODS) = """
                                 $(TYPEDSIGNATURES)
                                 $(DOCSTRING)
                                 """

# IS should not export any function since it can have name clashes with other packages.
# Do not add export statements.

# Every subtype must implement InfrastructureSystems.get_name() or have a field called name.
abstract type InfrastructureSystemsType end

"""
Return the internal forecast storage container or nothing, if the type doesn't store
forecasts.

Subtypes need to implement this method if they store forecasts.
"""
function get_forecasts(value::T) where {T <: InfrastructureSystemsType}
    return nothing
end

set_forecasts!(value::T) where {T <: InfrastructureSystemsType} = nothing

get_name(value::InfrastructureSystemsType) = value.name
set_name!(value::InfrastructureSystemsType, name) = value.name = name
get_internal(value::InfrastructureSystemsType) = value.internal

include("common.jl")
include("internal.jl")
include("utils/recorder_events.jl")
include("utils/flatten_iterator_wrapper.jl")
include("utils/generate_structs.jl")
include("utils/lazy_dict_from_iterator.jl")
include("utils/logging.jl")
include("utils/stdout_redirector.jl")
include("utils/utils.jl")
include("time_series_data.jl")
include("time_series_storage.jl")
include("hdf5_time_series_storage.jl")
include("in_memory_time_series_storage.jl")

include("forecasts.jl")
include("forecast_metadata.jl")
include("component.jl")
include("components.jl")
include("generated/includes.jl")
include("supplemental_constructors.jl")
include("forecast_parser.jl")
include("timeseries_formats.jl")
include("results.jl")
include("serialization.jl")
include("system_data.jl")
include("validation.jl")
include("utils/print.jl")
include("utils/test.jl")
include("units.jl")

end # module

ext: jl


text:
struct Template
	id::Int
  title::String
	date::String
	category::String
	rating::Int
	userid::Int
	locked::Any
	canread::String
	canwrite::String
	lastchange::DateTime
	up_item_id::Any
	has_attachment::Any
	next_step::Any
	category_id::Int
	color::String
	fullname::String
	recent_comment::Any
	has_comment::Int
	bookable::Int
end

# we need this to parse json strings as Julia structs
StructTypes.StructType(::Type{Template}) = StructTypes.Struct()


text:
#UpperBoundsUnimodal.jl

#####
#Function assumes that we have already translated to V_c space
###
#h 1d function that represents  E[h(V)] = 0
#H is a 1d function that represents int_[t, mode] h  (notice ordered limits)
#sep_fun(a, b, c, mode, pk, p_{k+1}) returns maximizers of
#	max_{t in [pk, p_{k+1}}  at + b t^2 + c H(t)
#returns val, tstar
function _vopp_unimodal(S, mode, h, H, sep_fun; N = 100, numCuts=100, TOL=1e-6, print_trace=false)
	m = Model(solver=GurobiSolver(OutputFlag = false))
	@variable(m, theta)
	@variable(m, lam[1:2])

	#build a grid containing mode
	ps = collect(range(0, stop=S, length=N+1)) #doc indexes from 0 to N
	if !(mode in ps)
		push!(ps, mode)
		sort!(ps)
	end
	N = length(ps) - 1
	jstar = findall(mode .== ps)[1]

	@variable(m, Q[1:N+1] >= 0)  #doc indexes from 0 to N
	@constraint(m, sum(Q) == 1.)

	#the constraint corresponding to t = m
	@constraint(m, theta + mode * lam[1] + h(mode) * lam[2] <=
				sum(ps[j] * Q[j] for j = 1:jstar)
				)

	@objective(m, Max, theta + (2 - mode) * lam[1])

	for iter = 1:numCuts
		status = solve(m)
		feasible = true

		print_trace && println("Starting Iteration: \t", iter, "\t", feasible)
		#t < mode constraints
		for k = 1:jstar-1
			a = sum(ps[j] * getvalue(Q[j]) for j = 1:k)
			a = a - getvalue(theta) + mode * getvalue(lam[1])
			b = -getvalue(lam[1])
			c =  getvalue(lam[2])

			val, tstar = sep_fun(a, b, c, mode, ps[k], ps[k+1])
			rhs = mode * sum(ps[j] * getvalue(Q[j]) for j = 1:k)
			rhs += sum(ps[j] * getvalue(Q[j]) * (mode - ps[j]) for j = (k+1):jstar)
			rhs -= mode * getvalue(theta)

			if val > rhs + TOL
				feasible = false
				print_trace && println("< iter $iter : \t $val \t $rhs")

				@constraint(m, (mode - tstar) * theta + (mode - tstar) * tstar * lam[1] +
								H(tstar) * lam[2] <= sum(ps[j] * Q[j] * (mode - tstar) for j = 1:k) +
												     sum(ps[j] * Q[j] * (mode - ps[j]) for j = (k+1):jstar)

					)
			end
		end
		print_trace && feasible && println("Iter $iter \t All t < m cnsts satisfied")

		#t > mode constraints
		for k = (jstar+1):N
			a = getvalue(theta) - mode * getvalue(lam[1])
			a -= sum(ps[j] * getvalue(Q[j]) for j = 1:k)
			b = getvalue(lam[1])
			c = -getvalue(lam[2])

			val, tstar = sep_fun(a, b, c, mode, ps[k], ps[k+1])
			rhs = mode * getvalue(theta)
			rhs -= mode * sum(ps[j] * getvalue(Q[j]) for j = 1:jstar)
			rhs -= sum(ps[j]^2 * getvalue(Q[j]) for j = (jstar + 1):k)

			if val > rhs + TOL
				feasible = false
				print_trace && println("> iter $iter : \t $val \t $rhs")

				@constraint(m, (tstar - mode) * theta + (tstar- mode) * tstar * lam[1] -
								H(tstar) * lam[2] <= sum(ps[j] * Q[j] * (tstar - mode) for j = 1:jstar) +
												     sum(ps[j] * Q[j] * (tstar - ps[j]) for j = (jstar+1):k)

					)
			end
		end
		print_trace && feasible && println("Iter $iter \t All t > m cnsts satisfied")

		if feasible
			if status != :Optimal
				println("Problem likely infeasible.  Dual Status:\t", status)
				return -1.0
			else #correctly solved
				return 1 / getobjectivevalue(m)
			end
		end
	end
	##Exhausted iteration count
	throw("Max Iteration Count Reached")
	return -1.
end

function vopp_ub_unimodal_MAD(mu, S, M, D, mode; numCuts=100, N=100, print_trace=false)
	#standardize
	c = mu * (1 - M)
	Sc = vopp.comp_Sc(S, M)
	mode_c = (mode - c) / (mu - c)

	#create subfunctions
	h(t) = M * abs(t - 1)/2 - D

	#H(t) = int_[t, mode_c] h
	function H(t)
		quadgk(h, t, 1, mode_c)[1]
	end

	#return maximizers of
	#max_{[pk, pk_p]}  at + bt^2 + c H(t)
	function sep_fun(a, b, c, mode, pk, pk_p)
		obj(t) = a * t + b * t^2 + c * H(t)
		val_k = obj(pk)
		val_k_p = obj(pk_p)

		tstar, val = 0., 0.
		if val_k > val_k_p
			tstar, val = pk, val_k
		else
			tstar, val = pk_p, val_k_p
		end

		t_less = (-2a - 2D + c * M) / (4b + c * M)
		if t_less <= 1 && pk <= t_less <= pk_p
			val_less = obj(t_less)
			if val_less >= val
				tstar, val = t_less, val_less
			end
		end

		t_more = (-2a - 2D - c * M) / (4b - c * M)
		if t_more >= 1 && pk <= t_more <= pk_p
			val_more = obj(t_more)
			if val_more >= val
				tstar, val = t_more, val_more
			end
		end
		return val, tstar
	end

	#Now kick it off!
	return _vopp_unimodal(Sc, mode_c, h, H, sep_fun, N=N, numCuts=numCuts, print_trace=print_trace)
end

function vopp_ub_unimodal_CV(mu, S, M, C, mode; N=100, numCuts=500, print_trace=false)
	#standardize
	c = mu * (1 - M)
	Sc = vopp.comp_Sc(S, M)
	mode_c = (mode - c) / (mu - c)

	#create subfunctions
	h(t) = M^2 * (t - 1)^2 - C^2

	#H(t) = int_[t, mode_c] h
	function H(t)
		quadgk(h, t, 1, mode_c)[1]
	end

	#return maximizers of
	#max_{[pk, pk_p]}  at + bt^2 + c H(t)
	function sep_fun(a, b, c, mode, pk, pk_p)
		obj(t) = a * t + b * t^2 + c * H(t)
		val_k = obj(pk)
		val_k_p = obj(pk_p)

		tstar, val = 0., 0.
		if val_k > val_k_p
			tstar, val = pk, val_k
		else
			tstar, val = pk_p, val_k_p
		end

		tp = b + c*M^2 + sqrt(b^2 + c *(a + 2b + c*C^2) * M^2)
		tp /= c*M^2

		if pk <= tp <= pk_p
			valp = obj(tp)
			if valp >= val
				tstar, val = tp, valp
			end
		end

		tm = b + c*M^2 - sqrt(b^2 + c *(a + 2b + c*C^2) * M^2)
		tm /= c*M^2

		if pk <= tm <= pk_p
			valm = obj(tm)
			if valm >= val
				tstar, val = tm, valm
			end
		end
		return val, tstar
	end

	#Now kick it off!
	return _vopp_unimodal(Sc, mode_c, h, H, sep_fun, N=N, numCuts=numCuts, print_trace=print_trace)
end

function vopp_ub_unimodal_GM(mu, S, M, B, mode; N=100, numCuts=500, print_trace=false)
	#standardize
	c = mu * (1 - M)
	Sc = vopp.comp_Sc(S, M)
	mode_c = (mode - c) / (mu - c)

	#create subfunctions
	h(t) = -log(M * t + 1 - M) + log(B/mu)

	#H(t) = int_[t, mode_c] h
	function H(t)
		quadgk(h, t, 1, mode_c)[1]
	end

	#return maximizers of
	#max_{[pk, pk_p]}  at + bt^2 + c H(t)
	function sep_fun(a, b, c, mode, pk, pk_p)
		obj(t) = a * t + b * t^2 + c * H(t)
		val_k = obj(pk)
		val_k_p = obj(pk_p)

		tstar, val = 0., 0.
		if val_k > val_k_p
			tstar, val = pk, val_k
		else
			tstar, val = pk_p, val_k_p
		end

		#handle degenerate case where c == 0
		if c == 0.
			tp = -a/2/b
			if pk <= tp <= pk_p
				valp = obj(tp)
				if valp >= val
					tstar, val = tp, valp
				end
			end
			return val, tstar
		end

		#else real work to be done
		inner = 2b/c * B/mu / M
    	inner *= exp( 2b/c * (1/M - 1) - a/c)

    	if inner < -1/MathConstants.e
    		tp = pk - 1
    	else
    		if inner < 0. ##two solutions.  Need to check -1 branch
				tp = 1- 1/M + c/2b * vopp.lambertw(inner)
				if pk <= tp <= pk_p
					valp = obj(tp)
					if valp >= val
						tstar, val = tp, valp
					end
				end
			end
			#always check principal branch
			tp= 1- 1/M + c/2b * vopp.lambertw0(inner)
		end

		if pk <= tp <= pk_p
			valp = obj(tp)
			if valp >= val
				tstar, val = tp, valp
			end
		end

		return val, tstar
	end

	#Now kick it off!
	return _vopp_unimodal(Sc, mode_c, h, H, sep_fun, N=N, numCuts=numCuts, print_trace=print_trace)
end

function vopp_ub_unimodal_IC(mu, S, M, phat, q, mode; N=100, numCuts=500, print_trace=false)
	#standardize
	c = mu * (1 - M)
	Sc = vopp.comp_Sc(S, M)
	mode_c = (mode - c) / (mu - c)

	#create subfunctions
	v0 = (phat - 1)/M + 1
	h(t) = t >= v0 ? 1. - q : -q

	#H(t) = int_[t, mode_c] h
	#lazy implementation.
	function H(t)
		quadgk(h, t, 1, mode_c)[1]
	end

	#return maximizers of
	#max_{[pk, pk_p]}  at + bt^2 + c H(t)
	function sep_fun(a, b, c, mode, pk, pk_p)
		obj(t) = a * t + b * t^2 + c * H(t)

		#check endpoints
		tstar, val = pk, obj(pk)
		if obj(pk_p) > val
			tstar = pk_p
			val = obj(pk_p)
		end

		#handle degenerate case where c == 0
		if c == 0.
			tp = -a/2/b
			if pk <= tp <= pk_p
				if obj(tp) >= val
					tstar, val = tp, obj(tp)
				end
			end
			return val, tstar
		end

		#else real work to be done
		#check the two possible critical points
		if b != 0
			t1 = c/2/b * (1 - q - a/c)
			if t1 >= v0 && (pk <= t1 < pk_p)
				if obj(t1) >= val
					val = obj(t1)
					tstar = t1
				end
			end

			t0 = -c/2/b * (q + a/c)
			if t0 < v0 && (pk <= t0 < pk_p)
				if obj(t0) >= val
					val = obj(t0)
					tstar = t0
				end
			end
		end
		return val, tstar
	end  #end definition of separator

	#Now kick it off!
	return _vopp_unimodal(Sc, mode_c, h, H, sep_fun, N=N, numCuts=numCuts, print_trace=print_trace)
end


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

print(xs...)   = print(stdout::IO, xs...)
println(xs...) = println(stdout::IO, xs...)
println(io::IO) = print(io, '\n')

function show end
function repr end

struct DevNull <: IO end
const devnull = DevNull()
write(::DevNull, ::UInt8) = 1
unsafe_write(::DevNull, ::Ptr{UInt8}, n::UInt)::Int = n
close(::DevNull) = nothing
wait_close(::DevNull) = wait()

let CoreIO = Union{Core.CoreSTDOUT, Core.CoreSTDERR}
    global write(io::CoreIO, x::UInt8) = Core.write(io, x)
    global unsafe_write(io::CoreIO, x::Ptr{UInt8}, nb::UInt) = Core.unsafe_write(io, x, nb)

    CoreIO = Union{CoreIO, DevNull}
    global read(::CoreIO, ::Type{UInt8}) = throw(EOFError())
    global isopen(::CoreIO) = true
    global isreadable(::CoreIO) = false
    global iswritable(::CoreIO) = true
    global flush(::CoreIO) = nothing
    global eof(::CoreIO) = true
    global wait_readnb(::CoreIO, nb::Int) = nothing
end

stdin = devnull
stdout = Core.stdout
stderr = Core.stderr

max_issues_repo_name: zhuangchaoxi/julia

max_forks_repo_name: zhuangchaoxi/julia


text:
__precompile__(true)

module Mocking

include("expr.jl")

export @patch, @mock, Patch, apply

const global GENERIC_ANONYMOUS = VERSION >= v"0.5-"

# When ENABLED is false the @mock macro is a noop.
global ENABLED = false
global PATCH_ENV = nothing

function __init__()
    # Attempt to detect when Mocking has been imported while running within Pkg.test()
    if isdefined(Base, :PROGRAM_FILE) && basename(Base.PROGRAM_FILE) == "runtests.jl"
        enable()
    end
end

function enable()
    ENABLED::Bool && return  # Abend early if enabled has already been set
    global ENABLED = true
    global PATCH_ENV = PatchEnv()

    # TODO: Support programatically disabling the use of the compilecache.
    opts = Base.JLOptions()
    if isdefined(opts, :use_compilecache) && Bool(opts.use_compilecache)
        warn("Mocking.jl will probably not work when compilecache is enabled. Please start Julia with `--compilecache=no`")
    end
end

immutable Patch
    signature::Expr
    body::Function
    modules::Set
    # translation::Dict

    function Patch(signature::Expr, body::Function, translation::Dict)
        trans = adjust_bindings(translation)
        absolute_binding!(signature.args[2:end], trans)  # TODO: Don't like that signature is modified
        modules = Set([v.args[1] for v in values(trans)])
        new(signature, body, modules)
    end
end

# TODO: Find non-eval way to determine module locations of Types
# evaling in the @patch scope seems to be problematic for pre-compliation
# first(methods(x)).sig.types[2:end]

# We can use the @patch macro to create a list of bindings used then pass that
# in as an array into Patch. At runtime the types and function names will be fully
# qualified

# We can support optional parameters and keywords by using generic functions on
# 0.4

function convert(::Type{Expr}, p::Patch)
    sig, body = p.signature, p.body
    params = call_parameters(sig)
    return :($sig = $body($(params...)))
end

macro patch(expr::Expr)
    if expr.head == :function
        name = expr.args[1].args[1]
        params = expr.args[1].args[2:end]
        body = expr.args[2]

    # Short-form function syntax
    elseif expr.head == :(=) && expr.args[1].head == :call
        name = expr.args[1].args[1]
        params = expr.args[1].args[2:end]
        body = expr.args[2]

    # Anonymous function syntax
    # elseif expr.head == :(->)
        # TODO: Determine how this could be supported
    else
        throw(ArgumentError("expression is not a function definition"))
    end

    # Determine the modules required for the parameter types
    bindings = unique(extract_bindings(params))

    signature = QuoteNode(Expr(:call, name, params...))

    # Need to evaluate the body of the function in the context of the `@patch` macro in
    # order to support closures.
    # func = Expr(:(->), Expr(:tuple, params...), body)
    func = Expr(:(=), Expr(:call, gensym(), params...), body)

    translations = []
    for b in bindings
        push!(translations, Expr(:(=>), QuoteNode(b), b))
    end

    return esc(:(Mocking.Patch( $signature, $func, Dict($(translations...)) )))
end

immutable PatchEnv
    mod::Module
    debug::Bool

    function PatchEnv(debug::Bool=false)
        # Be careful not to call this code during pre-compilation otherwise we'll see the
        # warning: "incremental compilation may be broken for this module"
        m = Core.eval(Mocking, :(module $(gensym()) end))
        new(m, debug)
    end
end

function PatchEnv(patches::Array{Patch}, debug::Bool=false)
    pe = PatchEnv(debug)
    apply!(pe, patches)
    return pe
end

function PatchEnv(patch::Patch, debug::Bool=false)
    pe = PatchEnv(debug)
    apply!(pe, patch)
    return pe
end

function apply!(pe::PatchEnv, p::Patch)
    for m in p.modules
        Core.eval(pe.mod, Expr(:import, splitbinding(m)...))
    end
    Core.eval(pe.mod, convert(Expr, p))
end

function apply!(pe::PatchEnv, patches::Array{Patch})
    for p in patches
        apply!(pe, p)
    end
end

function apply(body::Function, pe::PatchEnv)
    original_pe = get_active_env()
    set_active_env(pe)
    try
        return body()
    finally
        set_active_env(original_pe)
    end
end

function apply(body::Function, patches::Array{Patch}; debug::Bool=false)
    apply(body, PatchEnv(patches, debug))
end
function apply(body::Function, patch::Patch; debug::Bool=false)
    apply(body, PatchEnv(patch, debug))
end

function ismocked(pe::PatchEnv, func_name::Symbol, args::Tuple)
    if isdefined(pe.mod, func_name)
        func = Core.eval(pe.mod, func_name)
        types = map(typeof, tuple(args...))
        exists = method_exists(func, types)

        if pe.debug
            info("calling $func_name($(types...))")
            if exists
                m = first(methods(func, types))
                info("executing mocked function: $m")
            else
                m = first(methods(Core.eval(func_name), types))
                info("executing original function: $m")
            end
        end

        return exists
    end
    return false
end

set_active_env(pe::PatchEnv) = (global PATCH_ENV = pe)
get_active_env() = PATCH_ENV::PatchEnv

macro mock(expr)
    isa(expr, Expr) || error("argument is not an expression")
    expr.head == :call || error("expression is not a function call")
    ENABLED::Bool || return esc(expr)  # @mock is a no-op when Mocking is not ENABLED

    func = expr.args[1]
    func_name = QuoteNode(func)
    args = expr.args[2:end]

    env_var = gensym("env")
    args_var = gensym("args")

    result = quote
        local $env_var = Mocking.get_active_env()
        local $args_var = tuple($(args...))
        if Mocking.ismocked($env_var, $func_name, $args_var)
            $env_var.mod.$func($args_var...)
        else
            $func($args_var...)
        end
    end

    return esc(result)
end

end # module


text:
export gmm,
    shrink,
    estimate_component_prob,
    estimate_var_μ,
    rand_gmm

using Distributions
using StatsBase

"""
fit Gaussian mixture model with all means set to 0.

A SNP contributing h2 to heritability has μ ~ N(0, var = N*h2)
and z|μ ~ N(μ, 1.0)

Here we decompose distribution of z into null component z ~ N(0; 1)
and components z ~ N(0; var = 1 + N*h2/M) where h2 is total heritability

Due to population stratification etc. the null component might have
variance > 1. Then z needs to be rescaled such that variance of null component ≈ 1

dp0: Dirichlet pseudo counts for null component

"""
function gmm(x, π0, σ0; ϵ = 1e-3, maxiter = 1000, dp0 = 0, verbose = false, forcenull = false)

    n_data = length(x)
    n_cluster = length(π0)

    # Dirichlet prior
    dp = zeros(Float64, 1, n_cluster)
    dp[1] = dp0

    π = π0[:]'
    σ = σ0[:]'
    π_old = π0[:]'
    σ_old = σ0[:]'
    pc = zeros(Float64, n_data, n_cluster)

    llold = -Inf
    converged = false
    n_iter = 0

    for iter = 1:maxiter
        n_iter += 1
        # E-step
        # P(c_j|x_j) propto P(x_j|c_j) * P(c_j)
        for i = 1:n_data
            pc[i,:] = pdf.(Normal.(0., σ), x[i]) .* π
        end # i
        pc .= pc ./ sum(pc; dims = 2)

        # M-step
        π .= (sum(pc; dims = 1) .+ dp) ./ (n_data + dp0)
        σ  .= sqrt.(sum(x.^2 .* pc; dims = 1) ./ sum(pc; dims = 1))
        if forcenull
            σ[1] = 1.0
        end

        # loglikelihood
        # sum_i log( sum_c P(x_i|c)*P(c) )
        ll = 0.0
        for i = 1:n_data
            ll += log(sum(pdf.(Normal.(0., σ), x[i]) .* π))
        end

        verbose && println("Iteration $iter : LogLikelihood $ll")

        if (maximum(abs.((π_old .- π) ./ π_old)) < ϵ) && (maximum(abs.((σ_old .- σ) ./ σ_old)) < ϵ)
            converged = true
            break
        end

        llold = ll
        π_old .= π
        σ_old = σ

    end

    (π = π, σ = σ, converged = converged, ll = llold, iter = n_iter)
end

"""
estimate component probability for z

pc[i, k] : probability that z[i] belongs to component (π[k], σ[k])
"""
function estimate_component_prob(z, π, σ)
    @assert ndims(π) == ndims(σ) == 2 "π and σ need to be row vectors"
    @assert ndims(z) == 1 "z needs to be a column vector"
    @assert size(π) == size(σ) "π and σ need to have same dimensions"
    n_data = length(z)
    n_cluster = size(π, 2)
    pc = zeros(Float64, n_data, n_cluster)
    for i = 1:n_data
        pc[i,:] = pdf.(Normal.(0., σ), z[i]) .* π
    end
    pc ./ sum(pc; dims = 2)
end

function rand_gmm(π, σ)
    K = length(π)
    k = sample(1:K, Weights(π[:]))
    rand(Normal(0.0, σ[k]), 1)[1]
end

"""
compute shrinkage factor for z

i.e. estimate μ|z  where z|μ ~ N(μ; sd = 1)
and prior for μ given by  P(μ) = sum_k π[k] N(0; var = σ[k]^2 - 1)

σ2_μ for first component is set to 0
"""
function shrink(z, π, σ)

    pc = estimate_component_prob(z, π, σ)
    σ2_μ = max.(σ .^2 .- 1.0, 0.0)
    σ2_μ[1] = 0.0

    pc * (σ2_μ ./ (1.0 .+ σ2_μ))'
end

"""
estimate var_μ = N h2 / M for each SNP

N  number of samples
h2 total heritability
M  number of SNPs

π: component probability as obtained from gmm
σ: standard deviation of component as obtained from gmm

1/var_μ needs to be added to diagonal of LD matrix
"""
function estimate_var_μ(z, π, σ)

    pc = estimate_component_prob(z, π, σ)
    σ2_μ = max.(σ .^2 .- 1.0, 0.0)

    pc * σ2_μ'
end

function test_gmm()

    x = vcat(rand(Normal(0, 1.5), 100),
             rand(Normal(0,1.), 800),
             rand(Normal(0, 5.), 100))

    σ0 = [1., 2., 3.]
    π0 = [0.9, 0.05, 0.05]

    res = gmm(x, π0, σ0)
end


text:
export solve

"""
```julia
solve(nls::AbstractNLS,
      θ_init::AbstractVector,
      bc::BoundConstraints,
      conf::Abstract_BC_Solver_Conf) -> Abstract_Solver_Result
```

Generic interface to solve a [`AbstractNLS`](@ref) problem **with
bound constraints**.

The used algorithm is defined through [`Abstract_BC_Solver_Conf`](@ref)
specializations.

The method returns a [`Abstract_BC_Solver_Result`](@ref) specialization.
"""
solve(nls::AbstractNLS,
      θ_init::AbstractVector,
      bc::BoundConstraints,
      conf::Abstract_BC_Solver_Conf) = @assert(false,"To implement")


text:
using MUMPSjInv
using Test
using LinearAlgebra
using SparseArrays

include("getDivGrad.jl");

A = getDivGrad(32,32,16);
n = size(A,1);

# REAL: test for single rhs
println("Test for real SPD matrix: one sparse random rhs");
rhs = sprandn(n,1,0.05);

x = solveMUMPS(A,rhs,1);
err  =  norm(A*x-Matrix(rhs)) / norm(Matrix(rhs));
@test err < 1e-14

# rhs as sparse vector instead of n X 1 sparse matrix
rhs = SparseVector(rhs)

x = solveMUMPS(A,rhs,1);
err  =  norm(A*x-Vector(rhs)) / norm(Vector(rhs));
@test err < 1e-14

# REAL: test for multiple rhs
println("Test for real SPD matrix: multiple sparse random rhs");
nrhs = 10;
rhs = sprandn(n,nrhs,0.03);

x = solveMUMPS(A,rhs,1);

err = zeros(nrhs)
for i=1:nrhs
        err[i] =  norm(A*x[:,i]-Vector(rhs[:,i])) / norm(Vector(rhs[:,i]));
end
@test eltype(x) == Float64
@test maximum(err) < 1e-14

println("Test for real SPD matrix: multiple sparse delta function rhs");
rhs = sparse(1.0*I,n,n)
rhs = rhs[:,1:nrhs]

x = solveMUMPS(A,rhs,1);

err = zeros(nrhs)
for i=1:nrhs
        err[i] =  norm(A*x[:,i]-Vector(rhs[:,i])) / norm(Vector(rhs[:,i]));
end
@test eltype(x) == Float64
@test maximum(err) < 1e-14

# COMPLEX: test for one rhs
println("Test for complex SPD matrix: one rhs");
r = rand(n);
A = A + im*spdiagm(0=>r);

rhs = sprandn(n,1,0.03) + im*sprandn(n,1,0.03)

x = solveMUMPS(A,rhs,2);

err=  norm(A*x-Matrix(rhs)) / norm(Matrix(rhs))
@test eltype(x) == ComplexF64
@test err < 1e-14

# sparse vector rhs
rhs = SparseVector(sprandn(n,1,0.03) + im*sprandn(n,1,0.03))

x = solveMUMPS(A,rhs,2);

err=  norm(A*x-Vector(rhs)) / norm(Vector(rhs))
@test eltype(x) == ComplexF64
@test err < 1e-14

# COMPLEX : test for multiple rhs
println("Test for complex SPD matrix: multiple sparse random rhs");
nrhs = 10;
rhs = sprandn(n,nrhs,0.03) + im*sprandn(n,nrhs,0.03);

x = solveMUMPS(A,rhs,2);
err = zeros(nrhs)
for i=1:nrhs
        err[i] =  norm(A*x[:,i]-Vector(rhs[:,i])) / norm(Vector(rhs[:,i]));
end
@test eltype(x) == ComplexF64
@test maximum(err) < 1e-14

println("DONE!")


text:
"""
    IntervalTopology([domain::IntervalDomain, ]faces::AbstractVector)

A 1D topology of elements, with element boundaries at `faces`.
"""
struct IntervalTopology{D <: IntervalDomain, V} <: AbstractTopology
    domain::D
    faces::V
end

IntervalTopology(faces::AbstractVector) =
    IntervalTopology(IntervalDomain(faces[1], faces[end]), faces)


text:
#  Copyright 2016, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.
#############################################################################
# JuMP
# An algebraic modelling langauge for Julia
# See http://github.com/JuliaOpt/JuMP.jl
#############################################################################
# test/nonlinear.jl
# Test general nonlinear
#############################################################################
using JuMP, FactCheck

# If solvers not loaded, load them (i.e running just these tests)
!isdefined(:nlp_solvers) && include("solvers.jl")

facts("[nonlinear] Test getvalue on arrays") do
    m = Model()
    @variable(m, x, start = π/2)
    @NLexpression(m, f1, sin(x))
    @NLexpression(m, f2, sin(2x))
    @NLexpression(m, f3[i=1:2], sin(i*x))

    @fact getvalue(f1) --> roughly(1, 1e-5)
    @fact getvalue(f2) --> roughly(0, 1e-5)

    @fact getvalue([f1, f2]) --> getvalue(f3)

    v = [1.0, 2.0]
    @NLparameter(m, vparam[i=1:2] == v[i])
    @fact getvalue(vparam) --> v
    v[1] = 3.0
    setvalue(vparam, v)
    @fact getvalue(vparam) --> v
end

facts("[nonlinear] Test HS071 solves correctly") do
for nlp_solver in nlp_solvers
context("With solver $(typeof(nlp_solver))") do
    # hs071
    # Polynomial objective and constraints
    # min x1 * x4 * (x1 + x2 + x3) + x3
    # st  x1 * x2 * x3 * x4 >= 25
    #     x1^2 + x2^2 + x3^2 + x4^2 = 40
    #     1 <= x1, x2, x3, x4 <= 5
    # Start at (1,5,5,1)
    # End at (1.000..., 4.743..., 3.821..., 1.379...)
    m = Model(solver=nlp_solver)
    initval = [1,5,5,1]
    @variable(m, 1 <= x[i=1:4] <= 5, start=initval[i])
    @NLobjective(m, Min, x[1]*x[4]*(x[1]+x[2]+x[3]) + x[3])
    @NLconstraint(m, x[1]*x[2]*x[3]*x[4] >= 25)
    @NLconstraint(m, sum{x[i]^2,i=1:4} == 40)
    @fact MathProgBase.numconstr(m) --> 2
    status = solve(m)

    @fact status --> :Optimal
    @fact getvalue(x)[:] --> roughly(
        [1.000000, 4.742999, 3.821150, 1.379408], 1e-5)
end; end; end

facts("[nonlinear] Test HS071 solves correctly (no macros)") do
for nlp_solver in nlp_solvers
context("With solver $(typeof(nlp_solver))") do
    # hs071
    # Polynomial objective and constraints
    # min x1 * x4 * (x1 + x2 + x3) + x3
    # st  x1 * x2 * x3 * x4 >= 25
    #     x1^2 + x2^2 + x3^2 + x4^2 = 40
    #     1 <= x1, x2, x3, x4 <= 5
    # Start at (1,5,5,1)
    # End at (1.000..., 4.743..., 3.821..., 1.379...)
    m = Model(solver=nlp_solver)
    initval = [1,5,5,1]
    @variable(m, 1 <= x[i=1:4] <= 5, start=initval[i])
    JuMP.setNLobjective(m, :Min, :($(x[1])*$(x[4])*($(x[1])+$(x[2])+$(x[3])) + $(x[3])))
    JuMP.addNLconstraint(m, :($(x[1])*$(x[2])*$(x[3])*$(x[4]) >= 25))
    JuMP.addNLconstraint(m, :($(x[1])^2+$(x[2])^2+$(x[3])^2+$(x[4])^2 == 40))
    @fact MathProgBase.numconstr(m) --> 2
    status = solve(m)

    @fact status --> :Optimal
    @fact getvalue(x)[:] --> roughly(
        [1.000000, 4.742999, 3.821150, 1.379408], 1e-5)
end; end; end

if VERSION >= v"0.5-dev+5475"
eval("""
facts("[nonlinear] Test HS071 solves correctly (generators)") do
    # hs071
    # Polynomial objective and constraints
    # min x1 * x4 * (x1 + x2 + x3) + x3
    # st  x1 * x2 * x3 * x4 >= 25
    #     x1^2 + x2^2 + x3^2 + x4^2 = 40
    #     1 <= x1, x2, x3, x4 <= 5
    # Start at (1,5,5,1)
    # End at (1.000..., 4.743..., 3.821..., 1.379...)
    m = Model()
    initval = [1,5,5,1]
    @variable(m, 1 <= x[i=1:4] <= 5, start=initval[i])
    @NLobjective(m, Min, x[1]*x[4]*(x[1]+x[2]+x[3]) + x[3])
    @NLconstraint(m, x[1]*x[2]*x[3]*x[4] >= 25)
    @NLconstraint(m, sum(x[i]^2 for i=1:4) == 40)
    @fact MathProgBase.numconstr(m) --> 2
    status = solve(m)

    @fact status --> :Optimal
    @fact getvalue(x)[:] --> roughly(
        [1.000000, 4.742999, 3.821150, 1.379408], 1e-5)

    @fact isexpr(macroexpand(:(@NLconstraint(m, sum(x[i]^2 for i=1:4) == 40))),:error) --> true
end""");end

facts("[nonlinear] Test HS071 solves correctly, epigraph") do
for nlp_solver in nlp_solvers
context("With solver $(typeof(nlp_solver))") do
        # hs071, with epigraph formulation
        # Linear objective, nonlinear constraints
        # min t
        # st  t >= x1 * x4 * (x1 + x2 + x3) + x3
        #     ...
        m = Model(solver=nlp_solver)
        start = [1.0, 5.0, 5.0, 1.0]
        @variable(m, 1 <= x[i=1:4] <= 5, start = start[i])
        @variable(m, t, start = 100)
        @objective(m, Min, t)
        @NLconstraint(m, t >= x[1]*x[4]*(x[1]+x[2]+x[3]) + x[3])
        @NLconstraint(m, x[1]*x[2]*x[3]*x[4] >= 25)
        @NLconstraint(m, sum{x[i]^2,i=1:4} == 40)
        status = solve(m)

        @fact status --> :Optimal
        @fact getvalue(x)[:] --> roughly(
            [1.000000, 4.742999, 3.821150, 1.379408], 1e-5)
end; end; end

facts("[nonlinear] Test ifelse") do
for nlp_solver in nlp_solvers
(contains("$(typeof(nlp_solver))", "OsilSolver") || contains("$(typeof(nlp_solver))", "NLoptSolver")) && continue
context("With solver $(typeof(nlp_solver))") do
        m = Model(solver=nlp_solver)
        @variable(m, x, start = 2)
        # minimizer at smooth point, solvers should be okay
        @NLobjective(m, Min, ifelse( x <= 1, x^2, x) )
        status = solve(m)

        @fact status --> :Optimal
        @fact getvalue(x) --> roughly(0.0, 1e-5)
end; end; end

facts("[nonlinear] Accepting fixed variables") do
for nlp_solver in convex_nlp_solvers
for simplify in [true, false]
context("With solver $(typeof(nlp_solver)), simplify = $simplify") do
    m = Model(solver=nlp_solver, simplify_nonlinear_expressions=simplify)
    @variable(m, x == 0)
    @variable(m, y ≥ 0)
    @objective(m, Min, y)
    @NLconstraint(m, y ≥ x^2)
    EnableNLPResolve()
    for α in 1:4
        setvalue(x, α)
        solve(m)
        @fact getvalue(y) --> roughly(α^2, 1e-6)
    end
    DisableNLPResolve()
end; end; end; end

facts("[nonlinear] Test QP solve through NL pathway") do
for nlp_solver in convex_nlp_solvers
context("With solver $(typeof(nlp_solver))") do
    # Solve a problem with quadratic objective with linear
    # constraints, but force it to use the nonlinear code.
    m = Model(solver=nlp_solver)
    @variable(m, 0.5 <= x <=  2)
    @variable(m, 0.0 <= y <= 30)
    @NLparameter(m, param == 1.0)
    @objective(m, Min, (x+y)^2)
    @NLconstraint(m, x + y >= param)
    status = solve(m)

    @fact status --> :Optimal
    @fact m.objVal --> roughly(1.0, 1e-6)
    @NLexpression(m, lhs, x+y)
    @fact getvalue(x)+getvalue(y) --> roughly(1.0, 1e-6)
    @fact getvalue(lhs) --> roughly(1.0, 1e-6)

    setvalue(param,10)
    @fact m.internalModelLoaded --> true
    status = solve(m)
    @fact m.objVal --> roughly(10.0^2, 1e-6)
    @fact getvalue(x)+getvalue(y) --> roughly(10.0, 1e-6)
    @fact getvalue(lhs) --> roughly(10.0, 1e-6)

end; end; end

facts("[nonlinear] Test quad con solve through NL pathway") do
for nlp_solver in convex_nlp_solvers
context("With solver $(typeof(nlp_solver))") do
    # Solve a problem with linear objective with quadratic
    # constraints, but force it to use the nonlinear code.
    m = Model(solver=nlp_solver)
    @variable(m, -2 <= x <= 2)
    @variable(m, -2 <= y <= 2)
    @NLobjective(m, Min, x - y)
    @constraint(m, x + x^2 + x*y + y^2 <= 1)
    status = solve(m)

    @fact status --> :Optimal
    @fact getobjectivevalue(m) --> roughly(-1-4/sqrt(3), 1e-6)
    @fact getvalue(x) + getvalue(y) --> roughly(-1/3, 1e-3)
end; end; end

facts("[nonlinear] Test resolve with parameter") do
for nlp_solver in convex_nlp_solvers
for simplify in [true,false]
context("With solver $(typeof(nlp_solver)), simplify = $simplify") do
    m = Model(solver=nlp_solver, simplify_nonlinear_expressions=simplify)
    @variable(m, z)
    @NLparameter(m, x == 1.0)
    @NLobjective(m, Min, (z-x)^2)
    status = solve(m)
    @fact status --> :Optimal
    @fact getvalue(z) --> roughly(1.0, 1e-3)

    setvalue(x, 5.0)
    status = solve(m)
    @fact status --> :Optimal
    @fact getvalue(z) --> roughly(5.0, 1e-3)
end; end; end; end

facts("[nonlinear] Test two-sided nonlinear constraints") do
for nlp_solver in convex_nlp_solvers
context("With solver $(typeof(nlp_solver))") do
    m = Model(solver=nlp_solver)
    @variable(m, x)
    @NLobjective(m, Max, x)
    l = -1
    u = 1
    @NLconstraint(m, l <= x <= u)
    status = solve(m)

    @fact status --> :Optimal
    @fact getobjectivevalue(m) --> roughly(u, 1e-6)

    @NLobjective(m, Min, x)
    status = solve(m)

    @fact status --> :Optimal
    @fact getobjectivevalue(m) --> roughly(l, 1e-6)
end; end; end

facts("[nonlinear] Test two-sided nonlinear constraints (no macros)") do
for nlp_solver in convex_nlp_solvers
context("With solver $(typeof(nlp_solver))") do
    m = Model(solver=nlp_solver)
    @variable(m, x)
    JuMP.setNLobjective(m, :Max, x)
    l = -1
    u = 1
    JuMP.addNLconstraint(m, :($l <= $x <= $u))
    status = solve(m)

    @fact status --> :Optimal
    @fact getobjectivevalue(m) --> roughly(u, 1e-6)

    JuMP.setNLobjective(m, :Min, x)
    status = solve(m)

    @fact status --> :Optimal
    @fact getobjectivevalue(m) --> roughly(l, 1e-6)
end; end; end

facts("[nonlinear] Quadratic equality constraints") do
for nlp_solver in nlp_solvers
context("With solver $(typeof(nlp_solver))") do
    m = Model(solver=nlp_solver)
    @variable(m, 0 <= x[1:2] <= 1)
    @constraint(m, x[1]^2 + x[2]^2 == 1/2)
    @NLobjective(m, Max, x[1] - x[2])
    status = solve(m)

    @fact status --> :Optimal
    @fact getvalue(x) --> roughly([sqrt(1/2), 0], 1e-6)
end; end; end

facts("[nonlinear] Test mixed integer nonlinear problems") do
for minlp_solver in minlp_solvers
context("With solver $(typeof(minlp_solver))") do
    ## Solve test problem 1 (Synthesis of processing system) in
     # M. Duran & I.E. Grossmann, "An outer approximation algorithm for
     # a class of mixed integer nonlinear programs", Mathematical
     # Programming 36, pp. 307-339, 1986.  The problem also appears as
     # problem synthes1 in the MacMINLP test set.
    m = Model(solver=minlp_solver)
    x_U = [2,2,1]
    @variable(m, x_U[i] >= x[i=1:3] >= 0)
    @variable(m, y[4:6], Bin)
    @NLobjective(m, Min, 10 + 10*x[1] - 7*x[3] + 5*y[4] + 6*y[5] + 8*y[6] - 18*log(x[2]+1) - 19.2*log(x[1]-x[2]+1))
    @NLconstraints(m, begin
        0.8*log(x[2] + 1) + 0.96*log(x[1] - x[2] + 1) - 0.8*x[3] >= 0
        log(x[2] + 1) + 1.2*log(x[1] - x[2] + 1) - x[3] - 2*y[6] >= -2
        x[2] - x[1] <= 0
        x[2] - 2*y[4] <= 0
        x[1] - x[2] - 2*y[5] <= 0
        y[4] + y[5] <= 1
    end)
    status = solve(m)

    @fact status --> :Optimal
    @fact getobjectivevalue(m) --> roughly(6.00976, 1e-5)
    @fact getvalue(x)[:] --> roughly([1.30098, 0.0, 1.0], 1e-5)
    @fact getvalue(y)[:] --> roughly([0.0, 1.0, 0.0], 1e-5)
end; end; end

facts("[nonlinear] Test continuous relaxation of minlp test problem") do
for nlp_solver in nlp_solvers
context("With solver $(typeof(nlp_solver))") do
    ## Solve continuous relaxation of test problem 1 (Synthesis of processing system) in
     # M. Duran & I.E. Grossmann, "An outer approximation algorithm for
     # a class of mixed integer nonlinear programs", Mathematical
     # Programming 36, pp. 307-339, 1986.  The problem also appears as
     # problem synthes1 in the MacMINLP test set.
     # Introduce auxiliary nonnegative variable for the x[1]-x[2]+1 term
    m = Model(solver=nlp_solver)
    x_U = [2,2,1]
    @variable(m, x_U[i] >= x[i=1:3] >= 0)
    @variable(m, 1 >= y[4:6] >= 0)
    @variable(m, z >= 0, start=1)
    @NLobjective(m, Min, 10 + 10*x[1] - 7*x[3] + 5*y[4] + 6*y[5] + 8*y[6] - 18*log(x[2]+1) - 19.2*log(z))
    @NLconstraints(m, begin
        0.8*log(x[2] + 1) + 0.96*log(z) - 0.8*x[3] >= 0
        log(x[2] + 1) + 1.2*log(z) - x[3] - 2*y[6] >= -2
        x[2] - x[1] <= 0
        x[2] - 2*y[4] <= 0
        x[1] - x[2] - 2*y[5] <= 0
        y[4] + y[5] <= 1
        x[1] - x[2] + 1 == z
    end)
    status = solve(m)

    @fact status --> :Optimal
    @fact getobjectivevalue(m) --> roughly(0.7593, 5e-5)
    @fact getvalue(x)[:] --> roughly([1.1465, 0.54645, 1.0], 2e-4)
    @fact getvalue(y)[:] --> roughly([0.2732, 0.3, 0.0], 2e-4)
    @fact getvalue(z) --> roughly(1.6, 2e-4)
end; end; end

facts("[nonlinear] Test maximization objective") do
for nlp_solver in convex_nlp_solvers
context("With solver $(typeof(nlp_solver))") do
    # Solve a simple problem with a maximization objective
    m = Model(solver=nlp_solver)
    @variable(m, -2 <= x <= 2); setvalue(x, -1.8)
    @variable(m, -2 <= y <= 2); setvalue(y,  1.5)
    @NLobjective(m, Max, y - x)
    @constraint(m, x + x^2 + x*y + y^2 <= 1)

    @fact solve(m) --> :Optimal
    @fact getobjectivevalue(m) --> roughly(1+4/sqrt(3), 1e-6)
    @fact getvalue(x) + getvalue(y) --> roughly(-1/3, 1e-3)
end; end; end

facts("[nonlinear] Test maximization objective (embedded expressions)") do
for nlp_solver in convex_nlp_solvers
for simplify in [true,false]
context("With solver $(typeof(nlp_solver)), simplify = $simplify") do
    m = Model(solver=nlp_solver, simplify_nonlinear_expressions=simplify)
    @variable(m, -2 <= x <= 2); setvalue(x, -1.8)
    @variable(m, -2 <= y <= 2); setvalue(y,  1.5)
    @NLobjective(m, Max, y - x)
    @NLexpression(m, quadexpr, x + x^2 + x*y + y^2)
    @NLconstraint(m, quadexpr <= 1)

    @fact solve(m) --> :Optimal
    @fact getobjectivevalue(m) --> roughly(1+4/sqrt(3), 1e-6)
    @fact getvalue(x) + getvalue(y) --> roughly(-1/3, 1e-3)
    @fact getvalue(quadexpr) --> roughly(1, 1e-5)
    @NLexpression(m, quadexpr2, x + x^2 + x*y + y^2)
    @fact getvalue(quadexpr2) --> roughly(1, 1e-5)
end; end; end; end

facts("[nonlinear] Test infeasibility detection") do
for nlp_solver in convex_nlp_solvers
contains(string(typeof(nlp_solver)),"NLoptSolver") && continue
context("With solver $(typeof(nlp_solver))") do
    # (Attempt to) solve an infeasible problem
    m = Model(solver=nlp_solver)
    n = 10
    @variable(m, 0 <= x[i=1:n] <= 1)
    @NLobjective(m, Max, x[n])
    for i in 1:n-1
        @NLconstraint(m, x[i+1]-x[i] == 0.15)
    end
    @fact solve(m, suppress_warnings=true) --> :Infeasible
end; end; end

facts("[nonlinear] Test unboundedness detection") do
for nlp_solver in convex_nlp_solvers
contains(string(typeof(nlp_solver)),"NLoptSolver") && continue
context("With solver $(typeof(nlp_solver))") do
    # (Attempt to) solve an unbounded problem
    m = Model(solver=nlp_solver)
    @variable(m, x >= 0)
    @NLobjective(m, Max, x)
    @NLconstraint(m, x >= 5)
    @fact solve(m, suppress_warnings=true) --> :Unbounded
end; end; end

facts("[nonlinear] Test entropy maximization") do
for nlp_solver in convex_nlp_solvers
context("With solver $(typeof(nlp_solver))") do
    m = Model(solver=nlp_solver)
    N = 3
    @variable(m, x[1:N] >= 0, start = 1)
    @NLexpression(m, entropy[i=1:N], -x[i]*log(x[i]))
    @NLobjective(m, Max, sum{entropy[i], i = 1:N})
    @constraint(m, sum(x) == 1)

    @fact solve(m) --> :Optimal
    @fact norm(getvalue(x)[:] - [1/3,1/3,1/3]) --> roughly(0.0, 1e-4)
end; end; end

facts("[nonlinear] Test entropy maximization (reformulation)") do
for nlp_solver in convex_nlp_solvers
context("With solver $(typeof(nlp_solver))") do
    m = Model(solver=nlp_solver)
    idx = [1,2,3,4]
    @variable(m, x[idx] >= 0, start = 1)
    @variable(m, z[1:4], start = 0)
    @NLexpression(m, entropy[i=idx], -x[i]*log(x[i]))
    @NLobjective(m, Max, sum{z[i], i = 1:2} + sum{z[i]/2, i=3:4})
    @NLconstraint(m, z_constr1[i=1], z[i] <= entropy[i])
    @NLconstraint(m, z_constr1[i=2], z[i] <= entropy[i]) # duplicate expressions
    @NLconstraint(m, z_constr2[i=3:4], z[i] <= 2*entropy[i])
    @constraint(m, sum(x) == 1)

    @fact solve(m) --> :Optimal
    @fact norm([getvalue(x[i]) for i in idx] - [1/4,1/4,1/4,1/4]) --> roughly(0.0, 1e-4)
    @fact getvalue(entropy[1]) --> roughly(-(1/4)*log(1/4), 1e-4)
    @NLexpression(m, zexpr[i=1:4], z[i])
    @fact getvalue(zexpr[1]) --> roughly(-(1/4)*log(1/4), 1e-4)
end; end; end

facts("[nonlinear] Test derivatives of x^4, x < 0") do
for nlp_solver in convex_nlp_solvers
context("With solver $(typeof(nlp_solver))") do
    m = Model(solver=nlp_solver)
    @variable(m, x >= -1, start = -0.5)
    @NLobjective(m, Min, x^4)
    status = solve(m)

    @fact status --> :Optimal
    @fact getvalue(x) --> roughly(0.0, 1e-2)
end; end; end

facts("[nonlinear] Test nonlinear duals") do
for nlp_solver in nlp_solvers
for simplify in [true,false]
applicable(MathProgBase.getconstrduals, MathProgBase.NonlinearModel(nlp_solver)) || continue
context("With solver $(typeof(nlp_solver)), simplify = $simplify") do
    modA = Model(solver=nlp_solver, simplify_nonlinear_expressions=simplify)
    @variable(modA, x >= 0)
    @variable(modA, y <= 5)
    @variable(modA, 2 <= z <= 4)
    @variable(modA, 0 <= r[i=3:6] <= i)
    @NLobjective(modA, Min, -((x + y)/2.0 + 3.0)/3.0 - z - r[3])
    @constraint(modA, cons1, x+y >= 2)
    @constraint(modA, cons2, sum{r[i],i=3:5} <= (2 - x)/2.0)
    @NLconstraint(modA, cons3, 7.0*y <= z + r[6]/1.9)

    # Solution
    @fact solve(modA) --> :Optimal
    @fact getobjectivevalue(modA) --> roughly(-5.8446115, 1e-6)
    @fact getvalue(x)       --> roughly(0.9774436, 1e-6)
    @fact getvalue(y)       --> roughly(1.0225563, 1e-6)
    @fact getvalue(z)       --> roughly(4.0, 1e-6)
    @fact getvalue(r)[3]    --> roughly(0.5112781, 1e-6)
    @fact getvalue(r)[4]    --> roughly(0.0, 1e-6)
    @fact getvalue(r)[5]    --> roughly(0.0, 1e-6)
    @fact getvalue(r)[6]    --> roughly(6.0, 1e-6)

    # Reduced costs
    @fact getdual(x)    --> roughly( 0.0, 1e-6)
    @fact getdual(y)    --> roughly( 0.0, 1e-6)
    @fact getdual(z)    --> roughly(-1.0714286, 1e-6)
    @fact getdual(r)[3] --> roughly( 0.0, 1e-6)
    @fact getdual(r)[4] --> roughly(1.0, 1e-6)
    @fact getdual(r)[5] --> roughly(1.0, 1e-6)
    @fact getdual(r)[6] --> roughly(-0.03759398, 1e-6)

    # Row duals
    @fact getdual(cons1) --> roughly( 0.333333, 1e-6)
    @fact getdual(cons2) --> roughly(-1.0, 1e-6)
    @fact getdual(cons3) --> roughly(-0.0714286, 1e-6)
end; end; end; end

facts("[nonlinear] Test nonlinear duals (Max)") do
for nlp_solver in nlp_solvers
applicable(MathProgBase.getconstrduals, MathProgBase.NonlinearModel(nlp_solver)) || continue
context("With solver $(typeof(nlp_solver))") do
    modA = Model(solver=nlp_solver)
    @variable(modA, x >= 0)
    @variable(modA, y <= 5)
    @variable(modA, 2 <= z <= 4)
    @variable(modA, 0 <= r[i=3:6] <= i)
    @NLobjective(modA, Max, ((x + y)/2.0 + 3.0)/3.0 + z + r[3])
    @constraint(modA, cons1, x+y >= 2)
    @constraint(modA, cons2, sum{r[i],i=3:5} <= (2 - x)/2.0)
    cons3 = @NLconstraint(modA, 7.0*y <= z + r[6]/1.9)

    # Solution
    @fact solve(modA) --> :Optimal
    @fact getobjectivevalue(modA) --> roughly(5.8446115, 1e-6)
    @fact getvalue(x)       --> roughly(0.9774436, 1e-6)
    @fact getvalue(y)       --> roughly(1.0225563, 1e-6)
    @fact getvalue(z)       --> roughly(4.0, 1e-6)
    @fact getvalue(r)[3]    --> roughly(0.5112781, 1e-6)
    @fact getvalue(r)[4]    --> roughly(0.0, 1e-6)
    @fact getvalue(r)[5]    --> roughly(0.0, 1e-6)
    @fact getvalue(r)[6]    --> roughly(6.0, 1e-6)

    # Reduced costs
    @fact getdual(x)    --> roughly( 0.0, 1e-6)
    @fact getdual(y)    --> roughly( 0.0, 1e-6)
    @fact getdual(z)    --> roughly(1.0714286, 1e-6)
    @fact getdual(r)[3] --> roughly( 0.0, 1e-6)
    @fact getdual(r)[4] --> roughly(-1.0, 1e-6)
    @fact getdual(r)[5] --> roughly(-1.0, 1e-6)
    @fact getdual(r)[6] --> roughly(0.03759398, 1e-6)

    # Row duals
    @fact getdual(cons1) --> roughly(-0.333333, 1e-6)
    @fact getdual(cons2) --> roughly(1.0, 1e-6)
    @fact getdual(cons3) --> roughly(0.0714286, 1e-6)
end; end; end

facts("[nonlinear] Test changing objectives") do
for nlp_solver in nlp_solvers
context("With solver $(typeof(nlp_solver))") do
    m = Model(solver=nlp_solver)
    @variable(m, x >= 0)
    @variable(m, y >= 0)
    @objective(m, Max, x+y)
    @NLconstraint(m, x+2y <= 1)
    @fact solve(m) --> :Optimal
    @fact getvalue(x) --> roughly(1.0,1e-4)
    @fact getvalue(y) --> roughly(0.0,1e-4)
    @fact getobjectivevalue(m) --> roughly(1.0,1e-4)

    @objective(m, Max, 2x+y)
    @fact solve(m) --> :Optimal
    @fact getvalue(x) --> roughly(1.0,1e-4)
    @fact getvalue(y) --> roughly(0.0,1e-4)
    @fact getobjectivevalue(m) --> roughly(2.0,1e-4)
end; end; end

facts("[nonlinear] Test Hessian chunking code") do
for nlp_solver in nlp_solvers
context("With solver $(typeof(nlp_solver))") do
    m = Model(solver=nlp_solver)
    @variable(m, x[1:18] >= 1, start = 1.2)
    @NLobjective(m, Min, prod{x[i],i=1:18})
    @fact solve(m) --> :Optimal
    @fact getvalue(x) --> roughly(ones(18),1e-4)
end; end; end

if VERSION >= v"0.5-dev+5475"
eval("""
facts("[nonlinear] Test Hessian chunking code (generators)") do
    m = Model()
    @variable(m, x[1:18] >= 1, start = 1.2)
    @NLobjective(m, Min, prod(x[i] for i=1:18))
    @fact solve(m) --> :Optimal
    @fact getvalue(x) --> roughly(ones(18),1e-4)
end"""); end

#############################################################################
# Test that output is produced in correct MPB form
type DummyNLPSolver <: MathProgBase.AbstractMathProgSolver
end
type DummyNLPModel <: MathProgBase.AbstractNonlinearModel
end
MathProgBase.NonlinearModel(s::DummyNLPSolver) = DummyNLPModel()
function MathProgBase.loadproblem!(m::DummyNLPModel, numVar, numConstr, x_l, x_u, g_lb, g_ub, sense, d::MathProgBase.AbstractNLPEvaluator)
    MathProgBase.initialize(d, [:ExprGraph])
    objexpr = MathProgBase.obj_expr(d)
    facts("[nonlinear] Test NL MPB interface ($objexpr)") do
        @fact objexpr --> anyof(:(x[1]^x[2]), :(-1.0*x[1]+1.0*x[2]))
        @fact MathProgBase.isconstrlinear(d,1) --> true
        @fact MathProgBase.isconstrlinear(d,3) --> true
        @fact MathProgBase.constr_expr(d,1) --> :(2.0*x[1] + 1.0*x[2] <= 1.0)
        if VERSION >= v"0.5.0-dev+3231"
            @fact MathProgBase.constr_expr(d,2) --> :(2.0*x[1] + 1.0*x[2] <= -0.0)
        else
            @fact MathProgBase.constr_expr(d,2) --> :(2.0*x[1] + 1.0*x[2] <= 0.0)
        end
        @fact MathProgBase.constr_expr(d,3) --> :(-5.0 <= 2.0*x[1] + 1.0*x[2] <= 5.0)
        if numConstr > 3
            @fact MathProgBase.constr_expr(d,4) --> :(2.0*x[1]*x[1] + 1.0*x[2] + -2.0 >= 0)
            @fact MathProgBase.constr_expr(d,5) --> :(sin(x[1]) * cos(x[2]) - 5 == 0.0)
            @fact MathProgBase.constr_expr(d,6) --> :(1.0*x[1]^2 - 1.0 == 0.0)
            @fact MathProgBase.constr_expr(d,7) --> :(2.0*x[1]^2 - 2.0 == 0.0)
            @fact MathProgBase.constr_expr(d,8) --> :(-0.5 <= sin(x[1]) <= 0.5)
        end
    end
end
MathProgBase.setwarmstart!(m::DummyNLPModel,x) = nothing
MathProgBase.optimize!(m::DummyNLPModel) = nothing
MathProgBase.status(m::DummyNLPModel) = :Optimal
MathProgBase.getobjval(m::DummyNLPModel) = NaN
MathProgBase.getsolution(m::DummyNLPModel) = [1.0,1.0]
MathProgBase.setvartype!(m::DummyNLPModel,vartype) = @fact any(vartype .== :Fixed) --> false
function test_nl_mpb()
    m = Model(solver=DummyNLPSolver())
    @variable(m, x == 1)
    @variable(m, y, Bin)
    @objective(m, Min, -x+y)
    @constraint(m, 2x+y <= 1)
    @constraint(m, 2x+y <= 0)
    @constraint(m, -5 <= 2x+y <= 5)
    #solve(m) # FIXME maybe?
    lb,ub = JuMP.constraintbounds(m)
    @fact lb --> [-Inf,-Inf,-5.0]
    @fact ub --> [1.0,-0.0,5.0]

    @constraint(m, 2x^2+y >= 2)
    @NLconstraint(m, sin(x)*cos(y) == 5)
    @NLconstraint(m, nlconstr[i=1:2], i*x^2 == i)
    @NLconstraint(m, -0.5 <= sin(x) <= 0.5)
    solve(m)

    @NLobjective(m, Min, x^y)
    solve(m)

    lb,ub = JuMP.constraintbounds(m)
    @fact lb --> [-Inf,-Inf,-5.0,0.0,0.0,0.0,0.0,-0.5]
    @fact ub --> [1.0,-0.0,5.0,Inf,0.0,0.0,0.0,0.5]
end
test_nl_mpb()

facts("[nonlinear] Expression graph for linear problem") do
    m = Model()
    @variable(m, x)
    @constraint(m, 0 <= x <= 1)
    @objective(m, Max, x)
    d = JuMP.NLPEvaluator(m)
    MathProgBase.initialize(d, [:ExprGraph])
    @fact MathProgBase.obj_expr(d) --> :(+(1.0 * x[1]))
end

facts("[nonlinear] Expression graph for ifelse") do
    m = Model()
    @variable(m, x, start = 2)
    @NLobjective(m, Min, ifelse( x <= 1, x^2, x) )
    d = JuMP.NLPEvaluator(m)
    MathProgBase.initialize(d, [:ExprGraph])
    @fact MathProgBase.obj_expr(d) --> :(ifelse( x[1] <= 1, x[1]^2, x[1]))
end

facts("[nonlinear] Expression graphs for corner cases") do
    m = Model()
    @variable(m, x, start = 2)
    @constraint(m, 0 <= 1)
    @NLconstraint(m, x <= sum{0, i in []} + prod{1, i in []})
    d = JuMP.NLPEvaluator(m)
    MathProgBase.initialize(d, [:ExprGraph])
    @fact MathProgBase.constr_expr(d,1) --> :(0 <= 1.0)
    @fact MathProgBase.constr_expr(d,2) --> :(x[1] - (0 + 1) <= 0.0)
end

facts("[nonlinear] Hessians through MPB") do
    # Issue 435
    m = Model()
    @variable(m, a, start = 1)
    @variable(m, b, start = 2)
    @variable(m, c, start = 3)

    @NLexpression(m, foo, a * b + c^2)

    @NLobjective(m, Min, foo)
    d = JuMP.NLPEvaluator(m)
    MathProgBase.initialize(d, [:Hess])
    I,J = MathProgBase.hesslag_structure(d)
    V = zeros(length(I))
    MathProgBase.eval_hesslag(d, V, m.colVal, 1.0, Float64[])
    hess_raw = sparse(I,J,V)
    # Convert from lower triangular
    hess_sparse = hess_raw + hess_raw' - sparse(diagm(diag(hess_raw)))
    @fact hess_sparse --> roughly([0.0 1.0 0.0; 1.0 0.0 0.0; 0.0 0.0 2.0])

    # make sure we don't get NaNs in this case
    @NLobjective(m, Min, a * b + 3*c^2)
    d = JuMP.NLPEvaluator(m)
    MathProgBase.initialize(d, [:Hess])
    setvalue(c, -1.0)
    V = zeros(length(I))
    MathProgBase.eval_hesslag(d, V, m.colVal, 1.0, Float64[])
    hess_raw = sparse(I,J,V)
    hess_sparse = hess_raw + hess_raw' - sparse(diagm(diag(hess_raw)))
    @fact hess_sparse --> roughly([0.0 1.0 0.0; 1.0 0.0 0.0; 0.0 0.0 6.0])

    # Initialize again
    MathProgBase.initialize(d, [:Hess])
    V = zeros(length(I))
    MathProgBase.eval_hesslag(d, V, m.colVal, 1.0, Float64[])
    hess_raw = sparse(I,J,V)
    hess_sparse = hess_raw + hess_raw' - sparse(diagm(diag(hess_raw)))
    @fact hess_sparse --> roughly([0.0 1.0 0.0; 1.0 0.0 0.0; 0.0 0.0 6.0])
end

facts("[nonlinear] Hess-vec through MPB") do
    m = Model()
    @variable(m, a, start = 1)
    @variable(m, b, start = 2)
    @variable(m, c, start = 3)

    @NLobjective(m, Min, a*b + c^2)
    @constraint(m, c*b <= 1)
    @NLconstraint(m, a^2/2 <= 1)
    d = JuMP.NLPEvaluator(m)
    MathProgBase.initialize(d, [:HessVec])
    h = ones(3) # test that input values are overwritten
    v = [2.4,3.5,1.2]
    MathProgBase.eval_hesslag_prod(d, h, m.colVal, v, 1.0, [2.0,3.0])
    correct = [3.0 1.0 0.0; 1.0 0.0 2.0; 0.0 2.0 2.0]*v
    @fact h --> roughly(correct)
end

facts("[nonlinear] NaN corner case (#695)") do

    m = Model()
    x0 = 0.0
    y0 = 0.0

    @variable(m, x >= -1, start = 1.0)
    @variable(m, y, start = 2.0)

    @NLobjective(m, Min, (x - x0) /(sqrt(y0) + sqrt(y)))

    d = JuMP.NLPEvaluator(m)
    MathProgBase.initialize(d, [:HessVec])
    h = ones(2)
    v = [2.4,3.5]
    MathProgBase.eval_hesslag_prod(d, h, m.colVal, v, 1.0, Float64[])
    correct = [0.0 -1/(2*2^(3/2)); -1/(2*2^(3/2)) 3/(4*2^(5/2))]*v
    @fact h --> roughly(correct)
end

if length(convex_nlp_solvers) > 0
    facts("[nonlinear] Error on NLP resolve") do
        m = Model(solver=convex_nlp_solvers[1])
        @variable(m, x, start = 1)
        @NLobjective(m, Min, x^2)
        solve(m)
        setvalue(x, 2)
        @fact_throws ErrorException solve(m)
        EnableNLPResolve()
        status = solve(m)
        @fact status --> :Optimal
    end
end

mysquare(x) = x^2
function myf(x,y)
    return (x-1)^2+(y-2)^2
end

if length(convex_nlp_solvers) > 0
    facts("[nonlinear] User-defined functions") do
        JuMP.register(:myf, 2, myf, autodiff=true)
        JuMP.register(:myf_2, 2, myf, (g,x,y) -> (g[1] = 2(x-1); g[2] = 2(y-2)))
        JuMP.register(:mysquare, 1, mysquare, autodiff=true)
        JuMP.register(:mysquare_2, 1, mysquare, x-> 2x, autodiff=true)
        JuMP.register(:mysquare_3, 1, mysquare, x-> 2x, x -> 2.0)

        m = Model(solver=convex_nlp_solvers[1])

        @variable(m, x[1:2] >= 0.5)
        @NLobjective(m, Min, myf(x[1],mysquare(x[2])))

        d = JuMP.NLPEvaluator(m)
        MathProgBase.initialize(d, [:Grad])
        gradout = zeros(2)
        xval = [1,sqrt(2.0)]
        @fact MathProgBase.eval_f(d, xval) --> roughly(0.0,1e-10)
        MathProgBase.eval_grad_f(d, gradout, xval)
        @fact gradout --> roughly([0.0,0.0],1e-10)

        @fact solve(m) --> :Optimal

        @fact getvalue(x) --> roughly(xval)

        @NLobjective(m, Min, myf_2(x[1],mysquare_2(x[2])))

        d = JuMP.NLPEvaluator(m)
        MathProgBase.initialize(d, [:Grad])
        gradout = zeros(2)
        xval = [1,sqrt(2.0)]
        @fact MathProgBase.eval_f(d, xval) --> roughly(0.0,1e-10)
        MathProgBase.eval_grad_f(d, gradout, xval)
        @fact gradout --> roughly([0.0,0.0],1e-10)
        setvalue(x[1],0.5)
        setvalue(x[2],0.5)
        @fact solve(m) --> :Optimal

        @fact getvalue(x) --> roughly(xval)

        # Test just univariate functions because hessians are disabled
        # if any multivariate functions are present.
        @NLobjective(m, Min, mysquare(x[1]-1) + mysquare_2(x[2]-2) + mysquare_3(x[1]))
        @fact solve(m) --> :Optimal
        @fact getvalue(x) --> roughly([0.5,2.0],1e-4)

    end
end


text:
# Julia code
# Copyright (c) 2022 the SLICOTMath.jl developers
# Portions extracted from SLICOT-Reference distribution:
# Copyright (c) 2002-2020 NICONET e.V.

function run_mc01rd(datfile, io=stdout)
    NIN = 5
    NOUT = 6
    DP1MAX = 10
    DP2MAX = 10
    DP3MAX = 10
    LENP3 = max(DP1MAX+DP2MAX,DP3MAX)+1
    P1 = Array{Float64,1}(undef, DP1MAX+1)
    P2 = Array{Float64,1}(undef, DP2MAX+1)
    P3 = Array{Float64,1}(undef, DP1MAX+DP2MAX+DP3MAX+1)
    f = open(datfile,"r")
    readline(f)
    vs = split(readline(f))
    DP1 = parse(BlasInt, vs[1])
if ( DP1<=-2 || DP1>DP1MAX )
else

    vs = String[]
    _isz = DP1+1
    while length(vs) < _isz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    P1[1:_isz] .= parsex.(Float64, vs)

    vs = split(readline(f))
    DP2 = parse(BlasInt, vs[1])
if ( DP2<=-2 || DP2>DP2MAX )
else

    vs = String[]
    _isz = DP2+1
    while length(vs) < _isz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    P2[1:_isz] .= parsex.(Float64, vs)

    vs = split(readline(f))
    DP3 = parse(BlasInt, vs[1])
if ( DP3<=-2 || DP3>DP3MAX )
else

    vs = String[]
    _isz = DP3+1
    while length(vs) < _isz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    P3[1:_isz] .= parsex.(Float64, vs)

end # if
    vs = split(readline(f))
    ALPHA = parse(Float64, replace(vs[1],'D'=>'E'))
# interp call 1

    DP3, INFO = SLICOT.mc01rd!(DP1, DP2, DP3, ALPHA, P1, P2, P3)
    println(io, "DP3 = $DP3")
    @test INFO == 0
    INFO == 0 || return

if ( INFO!=0 )
else
if ( DP3>=0 )
end # if
end # if
end # if
end # if
    close(f)
end # run_X()


text:
"""
Linear Kernel
"""
struct ArdLinearKernel{T, VT<:AbstractVector{T}} <: AbstractKernel
	lσb::VT
	lσv::VT
	c::VT
end
@functor ArdLinearKernel
function (LK::ArdLinearKernel{T, VT})(x; λ=T(1e-6)) where {T, VT}
	length(LK.c) == size(x, 1) && length(LK.lσv) == size(x, 1) || throw(DimensionMismatch("size of offset parameter should be the same as number of features"))
	n = size(x, 2)
	shifted_x = x .- LK.c
	scaled_shifted_x = @. shifted_x*exp(-LK.lσv)
	d = inner_prod_metric(scaled_shifted_x)
	σb_square = exp(2*LK.lσb[1])
	(@. d+σb_square) + Diagonal(λ*ones(n))
end
function (LK::ArdLinearKernel{T, VT})(x, xo) where {T, VT}
	shifted_x = x .- LK.c
	scaled_shifted_x = @. shifted_x*exp(-LK.lσv)
	shifted_xo = xo .- LK.c
	scaled_shifted_xo = @. shifted_xo*exp(-LK.lσv)
	d = inner_prod_metric(scaled_shifted_x, scaled_shifted_xo)
	σb_square = exp(2*LK.lσb[1])
	d .+ σb_square
end

function ArdLinearKernel(n_features::Int)
	lσb = rand(1)
	lσv = rand(n_features)
	c  = rand(n_features)
	ArdLinearKernel(lσb, lσv, c)
end

reset(LK::ArdLinearKernel) = ArdLinearKernel(length(LK.c))

function Base.show(io::IO, kernel::ArdLinearKernel)
	print(io, "ArdLinearKernel(")
	print(io, "lσv=", kernel.lσv, ", ", "c=", kernel.c, ", ", "lσb=", kernel.lσb)
	print(io, ")")
end

## iso kernel
struct IsoLinearKernel{T, VT<:AbstractVector{T}} <: AbstractKernel
	lσv::VT
end
@functor IsoLinearKernel
function (ILK::IsoLinearKernel{T, VT})(x; λ=T(1e-6)) where {T, VT}
	n = size(x, 2)
	scaled_x = x .* exp(-ILK.lσv[1])
	d = inner_prod_metric(scaled_x)
	d + Diagonal(λ*ones(n))
end
function (ILK::IsoLinearKernel{T, VT})(x, xo) where {T, VT}
	scaled_x = x .* exp(-ILK.lσv[1])
	scaled_xo = xo .* exp(-ILK.lσv[1])
	d = inner_prod_metric(scaled_x, scaled_xo)
	d
end

function IsoLinearKernel()
	lσv = rand(1)
	IsoLinearKernel(lσv)
end

reset(IKL::IsoLinearKernel) = IsoLinearKernel()

function Base.show(io::IO, kernel::IsoLinearKernel)
	print(io, "IsoLinearKernel(", "lσv=", kernel.lσv[1], ")")
end

max_issues_repo_name: HamletWantToCode/GPFlux

max_forks_repo_name: HamletWantToCode/GPFlux


text:
"""
Exact real arithmetic using a Boehm scaled-integer representation.

Reference:
Keith Briggs, Implementing exact real arithmetic in Python, C++ and C.
*Theor. Comp. Sci.* **351** (2006), 74--81

"""
module ExactReals

export ExactReal, value_range

include("exact_reals.jl")
include("arithmetic.jl")

end # module


text:
using StaticArrays

function aspectsOf(format::vk.VkFormat)::vk.VkImageAspectFlags
    if format == vk.VK_FORMAT_UNDEFINED
        return 0
    elseif format == vk.VK_FORMAT_R4G4_UNORM_PACK8 ||
            format == vk.VK_FORMAT_R4G4B4A4_UNORM_PACK16 ||
            format == vk.VK_FORMAT_B4G4R4A4_UNORM_PACK16 ||
            format == vk.VK_FORMAT_R5G6B5_UNORM_PACK16 ||
            format == vk.VK_FORMAT_B5G6R5_UNORM_PACK16 ||
            format == vk.VK_FORMAT_R5G5B5A1_UNORM_PACK16 ||
            format == vk.VK_FORMAT_B5G5R5A1_UNORM_PACK16 ||
            format == vk.VK_FORMAT_A1R5G5B5_UNORM_PACK16 ||
            format == vk.VK_FORMAT_R8_UNORM ||
            format == vk.VK_FORMAT_R8_SNORM ||
            format == vk.VK_FORMAT_R8_USCALED ||
            format == vk.VK_FORMAT_R8_SSCALED ||
            format == vk.VK_FORMAT_R8_UINT ||
            format == vk.VK_FORMAT_R8_SINT ||
            format == vk.VK_FORMAT_R8_SRGB ||
            format == vk.VK_FORMAT_R8G8_UNORM ||
            format == vk.VK_FORMAT_R8G8_SNORM ||
            format == vk.VK_FORMAT_R8G8_USCALED ||
            format == vk.VK_FORMAT_R8G8_SSCALED ||
            format == vk.VK_FORMAT_R8G8_UINT ||
            format == vk.VK_FORMAT_R8G8_SINT ||
            format == vk.VK_FORMAT_R8G8_SRGB ||
            format == vk.VK_FORMAT_R8G8B8_UNORM ||
            format == vk.VK_FORMAT_R8G8B8_SNORM ||
            format == vk.VK_FORMAT_R8G8B8_USCALED ||
            format == vk.VK_FORMAT_R8G8B8_SSCALED ||
            format == vk.VK_FORMAT_R8G8B8_UINT ||
            format == vk.VK_FORMAT_R8G8B8_SINT ||
            format == vk.VK_FORMAT_R8G8B8_SRGB ||
            format == vk.VK_FORMAT_B8G8R8_UNORM ||
            format == vk.VK_FORMAT_B8G8R8_SNORM ||
            format == vk.VK_FORMAT_B8G8R8_USCALED ||
            format == vk.VK_FORMAT_B8G8R8_SSCALED ||
            format == vk.VK_FORMAT_B8G8R8_UINT ||
            format == vk.VK_FORMAT_B8G8R8_SINT ||
            format == vk.VK_FORMAT_B8G8R8_SRGB ||
            format == vk.VK_FORMAT_R8G8B8A8_UNORM ||
            format == vk.VK_FORMAT_R8G8B8A8_SNORM  ||
            format == vk.VK_FORMAT_R8G8B8A8_USCALED ||
            format == vk.VK_FORMAT_R8G8B8A8_SSCALED ||
            format == vk.VK_FORMAT_R8G8B8A8_UINT ||
            format == vk.VK_FORMAT_R8G8B8A8_SINT ||
            format == vk.VK_FORMAT_R8G8B8A8_SRGB ||
            format == vk.VK_FORMAT_B8G8R8A8_UNORM ||
            format == vk.VK_FORMAT_B8G8R8A8_SNORM ||
            format == vk.VK_FORMAT_B8G8R8A8_USCALED ||
            format == vk.VK_FORMAT_B8G8R8A8_SSCALED ||
            format == vk.VK_FORMAT_B8G8R8A8_UINT ||
            format == vk.VK_FORMAT_B8G8R8A8_SINT ||
            format == vk.VK_FORMAT_B8G8R8A8_SRGB ||
            format == vk.VK_FORMAT_A8B8G8R8_UNORM_PACK32 ||
            format == vk.VK_FORMAT_A8B8G8R8_SNORM_PACK32 ||
            format == vk.VK_FORMAT_A8B8G8R8_USCALED_PACK32 ||
            format == vk.VK_FORMAT_A8B8G8R8_SSCALED_PACK32 ||
            format == vk.VK_FORMAT_A8B8G8R8_UINT_PACK32 ||
            format == vk.VK_FORMAT_A8B8G8R8_SINT_PACK32 ||
            format == vk.VK_FORMAT_A8B8G8R8_SRGB_PACK32 ||
            format == vk.VK_FORMAT_A2R10G10B10_UNORM_PACK32 ||
            format == vk.VK_FORMAT_A2R10G10B10_SNORM_PACK32 ||
            format == vk.VK_FORMAT_A2R10G10B10_USCALED_PACK32 ||
            format == vk.VK_FORMAT_A2R10G10B10_SSCALED_PACK32 ||
            format == vk.VK_FORMAT_A2R10G10B10_UINT_PACK32 ||
            format == vk.VK_FORMAT_A2R10G10B10_SINT_PACK32 ||
            format == vk.VK_FORMAT_A2B10G10R10_UNORM_PACK32 ||
            format == vk.VK_FORMAT_A2B10G10R10_SNORM_PACK32 ||
            format == vk.VK_FORMAT_A2B10G10R10_USCALED_PACK32 ||
            format == vk.VK_FORMAT_A2B10G10R10_SSCALED_PACK32 ||
            format == vk.VK_FORMAT_A2B10G10R10_UINT_PACK32 ||
            format == vk.VK_FORMAT_A2B10G10R10_SINT_PACK32 ||
            format == vk.VK_FORMAT_R16_UNORM ||
            format == vk.VK_FORMAT_R16_SNORM ||
            format == vk.VK_FORMAT_R16_USCALED ||
            format == vk.VK_FORMAT_R16_SSCALED ||
            format == vk.VK_FORMAT_R16_UINT ||
            format == vk.VK_FORMAT_R16_SINT ||
            format == vk.VK_FORMAT_R16_SFLOAT ||
            format == vk.VK_FORMAT_R16G16_UNORM ||
            format == vk.VK_FORMAT_R16G16_SNORM ||
            format == vk.VK_FORMAT_R16G16_USCALED ||
            format == vk.VK_FORMAT_R16G16_SSCALED ||
            format == vk.VK_FORMAT_R16G16_UINT ||
            format == vk.VK_FORMAT_R16G16_SINT ||
            format == vk.VK_FORMAT_R16G16_SFLOAT ||
            format == vk.VK_FORMAT_R16G16B16_UNORM ||
            format == vk.VK_FORMAT_R16G16B16_SNORM ||
            format == vk.VK_FORMAT_R16G16B16_USCALED ||
            format == vk.VK_FORMAT_R16G16B16_SSCALED ||
            format == vk.VK_FORMAT_R16G16B16_UINT ||
            format == vk.VK_FORMAT_R16G16B16_SINT ||
            format == vk.VK_FORMAT_R16G16B16_SFLOAT ||
            format == vk.VK_FORMAT_R16G16B16A16_UNORM ||
            format == vk.VK_FORMAT_R16G16B16A16_SNORM ||
            format == vk.VK_FORMAT_R16G16B16A16_USCALED ||
            format == vk.VK_FORMAT_R16G16B16A16_SSCALED ||
            format == vk.VK_FORMAT_R16G16B16A16_UINT ||
            format == vk.VK_FORMAT_R16G16B16A16_SINT ||
            format == vk.VK_FORMAT_R16G16B16A16_SFLOAT ||
            format == vk.VK_FORMAT_R32_UINT ||
            format == vk.VK_FORMAT_R32_SINT ||
            format == vk.VK_FORMAT_R32_SFLOAT ||
            format == vk.VK_FORMAT_R32G32_UINT ||
            format == vk.VK_FORMAT_R32G32_SINT ||
            format == vk.VK_FORMAT_R32G32_SFLOAT ||
            format == vk.VK_FORMAT_R32G32B32_UINT ||
            format == vk.VK_FORMAT_R32G32B32_SINT ||
            format == vk.VK_FORMAT_R32G32B32_SFLOAT ||
            format == vk.VK_FORMAT_R32G32B32A32_UINT ||
            format == vk.VK_FORMAT_R32G32B32A32_SINT ||
            format == vk.VK_FORMAT_R32G32B32A32_SFLOAT ||
            format == vk.VK_FORMAT_R64_UINT ||
            format == vk.VK_FORMAT_R64_SINT ||
            format == vk.VK_FORMAT_R64_SFLOAT ||
            format == vk.VK_FORMAT_R64G64_UINT ||
            format == vk.VK_FORMAT_R64G64_SINT ||
            format == vk.VK_FORMAT_R64G64_SFLOAT ||
            format == vk.VK_FORMAT_R64G64B64_UINT ||
            format == vk.VK_FORMAT_R64G64B64_SINT ||
            format == vk.VK_FORMAT_R64G64B64_SFLOAT ||
            format == vk.VK_FORMAT_R64G64B64A64_UINT ||
            format == vk.VK_FORMAT_R64G64B64A64_SINT ||
            format == vk.VK_FORMAT_R64G64B64A64_SFLOAT ||
            format == vk.VK_FORMAT_B10G11R11_UFLOAT_PACK32 ||
            format == vk.VK_FORMAT_E5B9G9R9_UFLOAT_PACK32
        return vk.VK_IMAGE_ASPECT_COLOR_BIT
    elseif format == vk.VK_FORMAT_D16_UNORM ||
            format == vk.VK_FORMAT_X8_D24_UNORM_PACK32 ||
            format == vk.VK_FORMAT_D32_SFLOAT
        return vk.VK_IMAGE_ASPECT_DEPTH_BIT
    elseif format == vk.VK_FORMAT_S8_UINT
        return vk.VK_IMAGE_ASPECT_STENCIL_BIT
    elseif format == vk.VK_FORMAT_D16_UNORM_S8_UINT ||
            format == vk.VK_FORMAT_D24_UNORM_S8_UINT ||
            format == vk.VK_FORMAT_D32_SFLOAT_S8_UINT
        return vk.VK_IMAGE_ASPECT_DEPTH_BIT | vk.VK_IMAGE_ASPECT_STENCIL_BIT
    elseif format == vk.VK_FORMAT_BC1_RGB_UNORM_BLOCK ||
            format == vk.VK_FORMAT_BC1_RGB_SRGB_BLOCK ||
            format == vk.VK_FORMAT_BC1_RGBA_UNORM_BLOCK ||
            format == vk.VK_FORMAT_BC1_RGBA_SRGB_BLOCK ||
            format == vk.VK_FORMAT_BC2_UNORM_BLOCK ||
            format == vk.VK_FORMAT_BC2_SRGB_BLOCK ||
            format == vk.VK_FORMAT_BC3_UNORM_BLOCK ||
            format == vk.VK_FORMAT_BC3_SRGB_BLOCK ||
            format == vk.VK_FORMAT_BC4_UNORM_BLOCK ||
            format == vk.VK_FORMAT_BC4_SNORM_BLOCK ||
            format == vk.VK_FORMAT_BC5_UNORM_BLOCK ||
            format == vk.VK_FORMAT_BC5_SNORM_BLOCK ||
            format == vk.VK_FORMAT_BC6H_UFLOAT_BLOCK ||
            format == vk.VK_FORMAT_BC6H_SFLOAT_BLOCK ||
            format == vk.VK_FORMAT_BC7_UNORM_BLOCK ||
            format == vk.VK_FORMAT_BC7_SRGB_BLOCK ||
            format == vk.VK_FORMAT_ETC2_R8G8B8_UNORM_BLOCK ||
            format == vk.VK_FORMAT_ETC2_R8G8B8_SRGB_BLOCK ||
            format == vk.VK_FORMAT_ETC2_R8G8B8A1_UNORM_BLOCK ||
            format == vk.VK_FORMAT_ETC2_R8G8B8A1_SRGB_BLOCK ||
            format == vk.VK_FORMAT_ETC2_R8G8B8A8_UNORM_BLOCK ||
            format == vk.VK_FORMAT_ETC2_R8G8B8A8_SRGB_BLOCK ||
            format == vk.VK_FORMAT_EAC_R11_UNORM_BLOCK ||
            format == vk.VK_FORMAT_EAC_R11_SNORM_BLOCK ||
            format == vk.VK_FORMAT_EAC_R11G11_UNORM_BLOCK ||
            format == vk.VK_FORMAT_EAC_R11G11_SNORM_BLOCK ||
            format == vk.VK_FORMAT_ASTC_4x4_UNORM_BLOCK ||
            format == vk.VK_FORMAT_ASTC_4x4_SRGB_BLOCK ||
            format == vk.VK_FORMAT_ASTC_5x4_UNORM_BLOCK ||
            format == vk.VK_FORMAT_ASTC_5x4_SRGB_BLOCK ||
            format == vk.VK_FORMAT_ASTC_5x5_UNORM_BLOCK ||
            format == vk.VK_FORMAT_ASTC_5x5_SRGB_BLOCK ||
            format == vk.VK_FORMAT_ASTC_6x5_UNORM_BLOCK ||
            format == vk.VK_FORMAT_ASTC_6x5_SRGB_BLOCK ||
            format == vk.VK_FORMAT_ASTC_6x6_UNORM_BLOCK ||
            format == vk.VK_FORMAT_ASTC_6x6_SRGB_BLOCK ||
            format == vk.VK_FORMAT_ASTC_8x5_UNORM_BLOCK ||
            format == vk.VK_FORMAT_ASTC_8x5_SRGB_BLOCK ||
            format == vk.VK_FORMAT_ASTC_8x6_UNORM_BLOCK ||
            format == vk.VK_FORMAT_ASTC_8x6_SRGB_BLOCK ||
            format == vk.VK_FORMAT_ASTC_8x8_UNORM_BLOCK ||
            format == vk.VK_FORMAT_ASTC_8x8_SRGB_BLOCK ||
            format == vk.VK_FORMAT_ASTC_10x5_UNORM_BLOCK ||
            format == vk.VK_FORMAT_ASTC_10x5_SRGB_BLOCK ||
            format == vk.VK_FORMAT_ASTC_10x6_UNORM_BLOCK ||
            format == vk.VK_FORMAT_ASTC_10x6_SRGB_BLOCK ||
            format == vk.VK_FORMAT_ASTC_10x8_UNORM_BLOCK ||
            format == vk.VK_FORMAT_ASTC_10x8_SRGB_BLOCK ||
            format == vk.VK_FORMAT_ASTC_10x10_UNORM_BLOCK ||
            format == vk.VK_FORMAT_ASTC_10x10_SRGB_BLOCK ||
            format == vk.VK_FORMAT_ASTC_12x10_UNORM_BLOCK ||
            format == vk.VK_FORMAT_ASTC_12x10_SRGB_BLOCK ||
            format == vk.VK_FORMAT_ASTC_12x12_UNORM_BLOCK ||
            format == vk.VK_FORMAT_ASTC_12x12_SRGB_BLOCK ||
            format == vk.VK_FORMAT_PVRTC1_2BPP_UNORM_BLOCK_IMG ||
            format == vk.VK_FORMAT_PVRTC1_4BPP_UNORM_BLOCK_IMG ||
            format == vk.VK_FORMAT_PVRTC2_2BPP_UNORM_BLOCK_IMG ||
            format == vk.VK_FORMAT_PVRTC2_4BPP_UNORM_BLOCK_IMG ||
            format == vk.VK_FORMAT_PVRTC1_2BPP_SRGB_BLOCK_IMG ||
            format == vk.VK_FORMAT_PVRTC1_4BPP_SRGB_BLOCK_IMG ||
            format == vk.VK_FORMAT_PVRTC2_2BPP_SRGB_BLOCK_IMG ||
            format == vk.VK_FORMAT_PVRTC2_4BPP_SRGB_BLOCK_IMG
        return vk.VK_IMAGE_ASPECT_COLOR_BIT
    else
        error("lava/common/FormatInfo lava::aspectsOf(): Unknown ", format, ", assuming color.")
    end
    return vk.VK_IMAGE_ASPECT_COLOR_BIT
end

const TYPE_TO_VK_FORMAT_MAP = Dict([
    (Vec{1, UInt8}, vk.VK_FORMAT_R8_UNORM),
    (Vec{2, UInt8}, vk.VK_FORMAT_R8G8_UNORM),
    (Vec{3, UInt8}, vk.VK_FORMAT_R8G8B8_UNORM),
    (Vec{4, UInt8}, vk.VK_FORMAT_R8G8B8A8_UNORM),
    (Vec{1, Int8}, vk.VK_FORMAT_R8_SNORM),
    (Vec{2, Int8}, vk.VK_FORMAT_R8G8_SNORM),
    (Vec{3, Int8}, vk.VK_FORMAT_R8G8B8_SNORM),
    (Vec{4, Int8}, vk.VK_FORMAT_R8G8B8A8_SNORM),
    (Vec{1, UInt32}, vk.VK_FORMAT_R32_UINT),
    (Vec{2, UInt32}, vk.VK_FORMAT_R32G32_UINT),
    (Vec{3, UInt32}, vk.VK_FORMAT_R32G32B32_UINT),
    (Vec{4, UInt32}, vk.VK_FORMAT_R32G32B32A32_UINT),
    (Vec{1, Int32}, vk.VK_FORMAT_R32_SINT),
    (Vec{2, Int32}, vk.VK_FORMAT_R32G32_SINT),
    (Vec{3, Int32}, vk.VK_FORMAT_R32G32B32_SINT),
    (Vec{4, Int32}, vk.VK_FORMAT_R32G32B32A32_SINT),
    (Vec{1, Float32}, vk.VK_FORMAT_R32_SFLOAT),
    (Vec{2, Float32}, vk.VK_FORMAT_R32G32_SFLOAT),
    (Vec{3, Float32}, vk.VK_FORMAT_R32G32B32_SFLOAT),
    (Vec{4, Float32}, vk.VK_FORMAT_R32G32B32A32_SFLOAT)
])
function vkTypeOfFormat(type)::vk.VkFormat
    return get(TYPE_TO_VK_FORMAT_MAP, type, vk.VK_FORMAT_UNDEFINED)
end


text:
function collatzLength(n)
  i = 1
  while n != 1
    if mod(n,2) == 0
      n = n / 2
    else
      n = 3n + 1
    end
    i = i + 1
  end
  i
end

max = 1
maxAt = 1
for i in 1:1000000
  if collatzLength(i) > max
    max = collatzLength(i)
    maxAt = i
  end
end

println(maxAt)

max_stars_repo_name: vmchale/julia-problem-solving

max_issues_repo_name: vmchale/julia-problem-solving

max_forks_repo_name: vmchale/julia-problem-solving


text:
# ---
# title: 459. Repeated Substring Pattern
# id: problem459
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: String
# link: <https://leetcode.com/problems/repeated-substring-pattern/description/>
# hidden: true
# ---
#
# Given a non-empty string check if it can be constructed by taking a substring
# of it and appending multiple copies of the substring together. You may assume
# the given string consists of lowercase English letters only and its length
# will not exceed 10000.
#
#
#
# **Example 1:**
#
#
#
#     Input: "abab"
#     Output: True
#     Explanation: It's the substring "ab" twice.
#
#
# **Example 2:**
#
#
#
#     Input: "aba"
#     Output: False
#
#
# **Example 3:**
#
#
#
#     Input: "abcabcabcabc"
#     Output: True
#     Explanation: It's the substring "abc" four times. (And the substring "abcabc" twice.)
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
"""
Returns x^n mod p
"""
function exp!(x::BigInt, n::BigInt, p::BigInt = big(1000007))
    ans = 1
    x = x % p

    while n > 0
        if ((n & 1) != 0)
            ans = (ans * x) % p
        end
        x = (x * x) % p
        n = n >> 1
    end
    ans
end

exp!(x::T, n::T, p::T = 1000007) where {T <: Integer} =
exp!(big(x), big(n), big(p))

"""
Modular multiplicative inverse

a: Number whose inverse is to be found
n: modular base
"""
function inverse(a::BigInt, n::BigInt)
    g, x, _ = gcdx(a, n)

    if (g != 1)
        error("The inverse doesn't exist")
    end
    (x % n + n) % n
end

inverse(a::T, n::T) where {T <: Integer} =
inverse(big(a), big(n))

"""
Seive of Eratosthenes

Returns primes less than n
"""
function seive(n::BigInt)
    seive_list = ones(Bool, n)

    for i in 2:trunc(BigInt, sqrt(n))
        if seive_list[i] == true
           for j in (i*i):i:n
                seive_list[j] = false
           end
        end
   end
    filter(x -> seive_list[x] == true, 2:n)
end

seive(n::T) where {T <: Integer} =
seive(big(n))

"""
Euler totient function of n

phi!(n) = number of positive Numbers co-prime with n.
n: Int
"""
function phi!(n::Number)
    result = n
    p = 2

    while p * p <= n
        if (n % p == 0)
            while n % p == 0
                n /= p
            end
            result -= result / p
        end
        p += 1
    end

    if (n > 1)
        result -= result / n
    end
    Int(result)
end

"""
Returns a random Number of size bit_size bits
"""
function get_random_int(bit_size::T) where {T <: Integer}
    bit_vector = bitrand(bit_size)
    sum([(2 ^ BigInt((i - 1))) * bit for (i, bit) in enumerate(reverse(bit_vector))])
end

"""
Returns the Jacobi symbol
If b is prime, it returns the Legendre Symbol
"""
function jacobi!(a::BigInt, b::BigInt)
    if (b <= 0)
        error("b must be >= 1")
    end
    if (b & 1 == 0)
        error("b must be odd")
    end

    ans = 1
    if (a < 0)
        a = -a
        if (b % 4 == 3)
            ans = -ans
        end
    end

    while a != 0
        while a % 2 == 0
            a /= 2
            if ((b % 8) in [3, 5])
                ans = -ans
            end
        end
        a, b = b, a
        if a % 4 == 3 && b % 4 == 3
            ans = -ans
        end
        a = a % b
    end

    if (b == 1)
        return (ans)
    end
    0
end

jacobi!(a::T, b::T) where {T <: Integer} =
jacobi!(big(a), big(b))

"""
Generates a large a prime number by incremental search
"""
function generate_large_prime(bit_size::T, primality_test::Function = miller_rabin) where {T <: Integer}
    p = get_random_int(bit_size)

    if (p & 1 != 1)
        p += 1
    end
    while primality_test(p) == false
        p += 2
    end
    p
end

"""
Easier to convert Int to bytes
"""
to_bytes(x) = reinterpret(UInt8, x)


text:
## This file includes data type and structure
struct ScenarioData
    contin
    load24
end

struct PenaltyParameter
    z
    x
end

struct PHmultiplier
    z
    x
end


text:
export computeTensorizedSP

function computeTensorizedSP(m::Int,
                             α::AbstractVector{<:AbstractVector{<:Real}},β::AbstractVector{<:AbstractVector{<:Real}},
                             nodes::AbstractVector{<:AbstractVector{<:Real}},weights::AbstractVector{<:AbstractVector{<:Real}},
                             ind::AbstractMatrix{<:Int};
                             issymmetric::BitArray=falses(length(α)))
    m < 1 && throw(DomainError(m, "`dimension` has to be positive"))
    l,p = size(ind)
    l -= 1
    T = spzeros(reduce(+,[l^i for i in 0:m]))
    for tensor_ind in collect(with_replacement_combinations(0:l, m))
        index = 1 + reduce(+,[idx*l^(m-i) for (i,idx) in enumerate(tensor_ind)])
        T[index] = computeSP(tensor_ind, α, β, nodes, weights, ind; issymmetric=issymmetric)
    end
    return T
end
function computeTensorizedSP(m::Int, op::AbstractVector, ind::AbstractMatrix{<:Int})
    α, β = coeffs(op)
    nodes, weights = nw(op)
    computeTensorizedSP(m, α, β, nodes, weights, ind; issymmetric=issymmetric.(op))
end

function computeTensorizedSP(m::Int, mop::MultiOrthoPoly)
    any([typeof(op.quad) == typeof(EmptyQuad()) for op in mop.uni]) && throw(InconsistencyError("at least one quadrature rule missing"))
    computeTensorizedSP(m, mop.uni, mop.ind)
end

function computeTensorizedSP(m::Int, op::AbstractOrthoPoly)
    typeof(op.quad) == typeof(EmptyQuad()) && throw(InconsistencyError("no quadrature rule provided"))
    computeTensorizedSP(m, [op], calculateMultiIndices(1, deg(op)))
end


text:
# Based on your calculations, the planned course doesn't seem to make any
# sense. You find the submarine manual and discover that the process is
# actually slightly more complicated.
#
# In addition to horizontal position and depth, you'll also need to track a
# third value, aim, which also starts at 0. The commands also mean something
# entirely different than you first thought:
#
# down X increases your aim by X units.
# up X decreases your aim by X units.
# forward X does two things:
# It increases your horizontal position by X units.
# It increases your depth by your aim multiplied by X.
#
# Again note that since you're on a submarine, down and up do the opposite of
# what you might expect: "down" means aiming in the positive direction.
#
# Now, the above example does something different:
#
# forward 5 adds 5 to your horizontal position, a total of 5. Because your aim is 0, your depth does not change.
# down 5 adds 5 to your aim, resulting in a value of 5.
# forward 8 adds 8 to your horizontal position, a total of 13. Because your aim is 5, your depth increases by 8*5=40.
# up 3 decreases your aim by 3, resulting in a value of 2.
# down 8 adds 8 to your aim, resulting in a value of 10.
# forward 2 adds 2 to your horizontal position, a total of 15. Because your aim is 10, your depth increases by 2*10=20 to a total of 60.
#
# After following these new instructions, you would have a horizontal position
# of 15 and a depth of 60. (Multiplying these produces 900.)
#
# Using this new interpretation of the commands, calculate the horizontal
# position and depth you would have after following the planned course. What do
# you get if you multiply your final horizontal position by your final depth?

function mapper(l)
  (dir, amt) = split(l, ' ')
  (dir, parse(Int, amt))
end
function reducer((d, h, a), (dir, amt))
  if dir == "down"
    (d, h, a + amt)
  elseif dir == "up"
    (d, h, a - amt)
  else
    (d + a * amt, h + amt, a)
  end
end

(depth, horz, _) = mapfoldl(mapper, reducer, eachline("input.txt"); init=(0, 0, 0))
println("Result: ", depth * horz)


text:
using Test
import GeoInterface
import ArchGDAL; const AG = ArchGDAL

@testset "Reproject a Geometry" begin
    @testset "Method 1" begin
        GDAL.__init__()
        source = AG.unsafe_importEPSG(2927); target = AG.unsafe_importEPSG(4326)
            transform = AG.unsafe_createcoordtrans(source, target)
                point = AG.unsafe_fromWKT("POINT (1120351.57 741921.42)")
                    @test AG.toWKT(point) == "POINT (1120351.57 741921.42)"
                    AG.transform!(point, transform)
                    @test GeoInterface.coordinates(point) ≈ [-122.598135, 47.348801]
                AG.destroy(point)
            AG.destroy(transform)
        AG.destroy(target); AG.destroy(source)
    end

    @testset "Method 2" begin
        GDAL.__init__()
        AG.importEPSG(2927) do source; AG.importEPSG(4326) do target
            AG.createcoordtrans(source, target) do transform
                AG.fromWKT("POINT (1120351.57 741921.42)") do point
                    @test AG.toWKT(point) == "POINT (1120351.57 741921.42)"
                    AG.transform!(point, transform)
                    @test GeoInterface.coordinates(point) ≈ [-122.598135, 47.348801]
        end end end end
    end
end

@testset "Get Projection" begin
    AG.registerdrivers() do
        AG.read("ospy/data1/sites.shp") do dataset
            layer = AG.getlayer(dataset, 0)
            spatialref = AG.getspatialref(layer)
            @test AG.toWKT(spatialref)[1:6] == "PROJCS"
            @test AG.toWKT(spatialref, false)[1:6] == "PROJCS"
            @test AG.toPROJ4(spatialref) == "+proj=utm +zone=12 +datum=WGS84 +units=m +no_defs "
            @test AG.toXML(spatialref)[1:17] == "<gml:ProjectedCRS"
            @test AG.toMICoordSys(spatialref) == "Earth Projection 8, 104, \"m\", -111, 0, 0.9996, 500000, 0"
            AG.nextfeature(layer) do feature
                @test AG.toPROJ4(AG.getspatialref(AG.getgeom(feature))) == "+proj=utm +zone=12 +datum=WGS84 +units=m +no_defs "
            end
        end
    end
    GDAL.__init__()
    AG.importEPSG(26912) do spatialref
        @test AG.toPROJ4(spatialref) == "+proj=utm +zone=12 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs "
        @test AG.toWKT(spatialref)[1:6] == "PROJCS"
        AG.morphtoESRI!(spatialref)
        @test AG.toPROJ4(spatialref) == "+proj=utm +zone=12 +ellps=GRS80 +units=m +no_defs "
    end
end


text:
using Flux, Metalhead
using Flux: onehotbatch, onecold, crossentropy, throttle
using Metalhead: trainimgs
using Base.Iterators: partition

# VGG16 and VGG19 models

vgg16() = Chain(
  Conv((3, 3), 3 => 64, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(64),
  Conv((3, 3), 64 => 64, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(64),
  x -> maxpool(x, (2, 2)),
  Conv((3, 3), 64 => 128, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(128),
  Conv((3, 3), 128 => 128, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(128),
  x -> maxpool(x, (2,2)),
  Conv((3, 3), 128 => 256, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(256),
  Conv((3, 3), 256 => 256, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(256),
  Conv((3, 3), 256 => 256, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(256),
  x -> maxpool(x, (2, 2)),
  Conv((3, 3), 256 => 512, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(512),
  Conv((3, 3), 512 => 512, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(512),
  Conv((3, 3), 512 => 512, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(512),
  x -> maxpool(x, (2, 2)),
  Conv((3, 3), 512 => 512, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(512),
  Conv((3, 3), 512 => 512, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(512),
  Conv((3, 3), 512 => 512, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(512),
  x -> maxpool(x, (2, 2)),
  x -> reshape(x, :, size(x, 4)),
  Dense(512, 4096, relu),
  Dropout(0.5),
  Dense(4096, 4096, relu),
  Dropout(0.5),
  Dense(4096, 10),
  softmax) |> gpu

vgg19() = Chain(
  Conv((3, 3), 3 => 64, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(64),
  Conv((3, 3), 64 => 64, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(64),
  x -> maxpool(x, (2, 2)),
  Conv((3, 3), 64 => 128, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(128),
  Conv((3, 3), 128 => 128, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(128),
  x -> maxpool(x, (2, 2)),
  Conv((3, 3), 128 => 256, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(256),
  Conv((3, 3), 256 => 256, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(256),
  Conv((3, 3), 256 => 256, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(256),
  Conv((3, 3), 256 => 256, relu, pad=(1, 1), stride=(1, 1)),
  x -> maxpool(x, (2, 2)),
  Conv((3, 3), 256 => 512, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(512),
  Conv((3, 3), 512 => 512, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(512),
  Conv((3, 3), 512 => 512, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(512),
  Conv((3, 3), 512 => 512, relu, pad=(1, 1), stride=(1, 1)),
  x -> maxpool(x, (2, 2)),
  Conv((3, 3), 512 => 512, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(512),
  Conv((3, 3), 512 => 512, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(512),
  Conv((3, 3), 512 => 512, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(512),
  Conv((3, 3), 512 => 512, relu, pad=(1, 1), stride=(1, 1)),
  x -> maxpool(x, (2, 2)),
  x -> reshape(x, :, size(x, 4)),
  Dense(512, 4096, relu),
  Dropout(0.5),
  Dense(4096, 4096, relu),
  Dropout(0.5),
  Dense(4096, 10),
  softmax) |> gpu

# Function to convert the RGB image to Float64 Arrays

getarray(X) = float.(permutedims(channelview(im), (2, 3, 1)))

# Fetching the train and validation data and getting them into proper shape

X = trainimgs(CIFAR10)
imgs = [getarray(X[i].img) for i in 1:50000]
labels = onehotbatch([X[i].ground_truth.class for i in 1:50000],1:10)
train = gpu.([(cat(4, imgs[i]...), labels[:,i]) for i in partition(1:49000, 1000)])
valset = collect(49001:50000)
valX = cat(4, imgs[valset]...) |> gpu
valY = labels[:, valset] |> gpu

# Defining the loss and accuracy functions

m = vgg16()

loss(x, y) = crossentropy(m(x), y)

accuracy(x, y) = mean(onecold(m(x), 1:10) .== onecold(y, 1:10))

# Defining the callback and the optimizer

evalcb = throttle(() -> @show(accuracy(valX, valY)), 10)

opt = ADAM(params(m))

# Starting to train models

Flux.train!(loss, train, opt, cb = evalcb)

# Fetch the test data from Metalhead and get it into proper shape.
# CIFAR-10 does not specify a validation set so valimgs fetch the testdata instead of testimgs

test = valimgs(CIFAR10)

testimgs = [getarray(test[i].img) for i in 1:10000]
testY = onehotbatch([test[i].ground_truth.class for i in 1:10000], 1:10) |> gpu
testX = cat(4, testimgs...) |> gpu

# Print the final accuracy

@show(accuracy(testX, testY))

max_stars_repo_name: lcy-seso/model-zoo

max_issues_repo_name: lcy-seso/model-zoo

max_forks_repo_name: lcy-seso/model-zoo


text:
"""
This module represents an explicit, direct sparse 2D finite-difference Poisson solver for heterogeneous media,
i.e. media having spatially varying (space-dependent) medium parameters.
The following functionality is currently available in this module:
* 'Forward' problem: given the source, and the conductivity distribution, solve for the electrical potential ψ:
	 Example problem: solve Aψ=∇⋅j; A=∇⋅(σ(x,z)∇), for ψ.
* 'Inverse' problem: given wave propagation related pore pressure, and mechanical medium properties, calculate the source term:
	 Example problem: solve Aψ=-∇⋅j; A=∇⋅([Q*k/η](x,z)∇B*P), for ∇⋅j;
* The Boundary Value Problem that is currently implemented assumes Neumann boundary conditions at all boundaries.
Developed by:\\\n
Niels Grobbe, Massachusetts Institute of Technology, USA.\\\n
In collaboration with: Aimé Fournier & Laurent Demanet, Massachusetts Institute of Technology, USA.\\\n
Date: October, 2017 \\\n
Contact: ngrobbe@gmail.com
"""
module Poisson

include("core.jl")
include("expt.jl")
include("expt_born.jl")
include("getprop.jl")
include("inversion.jl")

end # end module Poisson


text:
import Base: copy
export Clamp, constant, placeholder, @ensureVariables

"""
Description:

    A factor that clamps a variable to a constant value.

    f(out) = δ(out - value)

Interfaces:

    1. out

Construction:

    Clamp(out, value, id=:some_id)
    Clamp(value, id=:some_id)
"""
mutable struct Clamp{T<:VariateType} <: DeltaFactor
    id::Symbol
    interfaces::Vector{Interface}
    i::Dict{Symbol,Interface}
    value::Any

    # Fields for algorithm assembly
    dist_or_msg::Type # Define whether the clamp outputs a Message or ProbabilityDistribution
    buffer_id::Symbol # Specify the buffer id for an attached placeholder
    buffer_index::Int64 # Specify the buffer index when the placeholder attaches a vector

    Clamp(V::Type) = new{V}() # Instantiate an empty clamp without adding it to the graph

    function Clamp(out::Variable, value::Any; id=generateId(Clamp{variateType(value)}))
        self = new{variateType(value)}(id, Array{Interface}(undef, 1), Dict{Symbol,Interface}(), value)
        addNode!(currentGraph(), self)
        self.i[:out] = self.interfaces[1] = associate!(Interface(self), out)

        return self
    end
end

"""
Make a shallow copy of a Clamp without adding it to the graph.
The copied Clamp can be used as altered inboud during algorithm
assembly, without altering the original source Clamp.
"""
function copy(src::Clamp{V}) where V<:VariateType
    nd = Clamp(V)
    for field in fieldnames(Clamp)
        isdefined(src, field) && setfield!(nd, field, getfield(src, field))
    end

    return nd
end

variateType(::Number) = Univariate
variateType(::AbstractVector) = Multivariate
variateType(::AbstractMatrix) = MatrixVariate

"""
`constant` creates a `Variable` which is linked to a new `Clamp`,
and returns this variable.

    y = constant(3.0, id=:y)
"""
function constant(value::Any; id=generateId(Clamp{variateType(value)}))
    # This is basically an outer constructor for Clamp,
    # but the function returns the new variable linked to the Clamp.
    var = Variable(id=id)
    Clamp(var, value, id=id)

    return var
end

"""
`@ensureVariables(...)` casts all non-`Variable` arguments to `Variable` through constant(arg).
"""
macro ensureVariables(args...)
    checks = map(arg -> :(isa($arg, Variable) || ($arg = constant($arg))), args)
    expr   = reduce((current, item) -> :($current; $item), checks)
    return esc(expr)
end

"""
`placeholder(...)` creates a `Clamp` node and registers this node
as a data placeholder with the current graph.

    # Link variable y to buffer with id :y,
    # indicate that Clamp will hold Float64 values.
    placeholder(y, :y, datatype=Float64)

    # Link variable y to index 3 of buffer with id :y.
    # Specify the data type by passing a default value for the Clamp.
    placeholder(y, :y, index=3, default=0.0)

    # Indicate that the Clamp will hold an array of size `dims`,
    # with Float64 elements.
    placeholder(X, :X, datatype=Float64, dims=(3,2))
"""
function placeholder(   var::Variable,
                        buffer_id::Symbol;
                        index::Int=0,
                        default::Any=nothing,
                        datatype::DataType=Float64,
                        dims::Tuple=())

    # Build placeholder id
    constant_id = :placeholder_*var.id

    # Define a default value
    if default != nothing
        value = default
    elseif isempty(dims)
        value = zero(datatype)
    else
        value = zeros(datatype, dims)
    end

    node = Clamp(var, value, id=constant_id)
    current_graph.placeholders[node] = (buffer_id, index)

    return var
end

function placeholder(   buffer_id::Symbol;
                        index::Int=0,
                        default::Any=nothing,
                        datatype::DataType=Float64,
                        dims::Tuple=(),
                        var_id=buffer_id)

    var = Variable(id=var_id)
    return placeholder(var, buffer_id, index=index, default=default, datatype=datatype, dims=dims)
end