text:
mutable struct abstract_var
end

mutable struct context
end

mutable struct equtree
end

mutable struct equnode
end

mutable struct reshop_model
end

mutable struct reshop_var
end

mutable struct reshop_equ
end

mutable struct reshop_sp_matrix
end

mutable struct mathprgm
end

mutable struct ovf_definition
end

mutable struct reshop_options
end

mutable struct equil
end

mutable struct filter_ops
end

CONE_NONE = Cint(0)          # Unset/non-existent */
CONE_R_PLUS = Cuint(1)       # Non-negative real \f$\mathbb{R}_+\f$ */
CONE_R_MINUS = Cuint(2)      # Non-positive real \f$\mathbb{R}_-\f$  */
CONE_R = Cuint(3)            # Real \f$\mathbb{R}\f$ */
CONE_0 = Cuint(4)            # Zero singleton */
CONE_POLYHEDRAL = Cuint(5)   # Polyhedral cone */
CONE_SOC = Cuint(6)          # Second Order cone */
CONE_RSOC = Cuint(7)         # Rotated Second Order cone */
CONE_EXP = Cuint(8)          # Exponential cone */
CONE_DEXP = Cuint(9)         # Dual Exponential cone */
CONE_POWER = Cuint(10)       # Power cone */
CONE_DPOWER = Cuint(11)      # Dual Power cone */
__CONES_LEN = Cuint(12)

mutable struct option
	name::Vector{Cchar}
	typ::Cint
	value::Cdouble
end

# TODO with latest Julia
if VERSION < v"0.7"
	iswin = is_windows()
else
	iswin = Sys.iswindows()
end

#const libreshop = iswin ? "reshop" : "libreshop"
const term_str = iswin ? "\r\n" : "\n"

macro chk_index_or_size(idx)
quote
    if $(esc(idx)) >= reshop_valid_index_max
        error("the index/size returned by ReSHOP is not valid: upper limit is $(reshop_valid_index_max), value is ", $(esc(idx)))
    end
end
end

function ctx_add_lin_var(ctx::Ptr{context}, eidx, vidx, coeff::Cdouble)
	res = ccall((:rhp_equ_addvar, libreshop), Cint, (Ptr{context}, Cint, Cint, Cdouble), ctx, eidx, vidx, coeff)
	res != 0 && error("return code $res from ReSHOP")
end

function rhp_equ_add_linear(ctx::Ptr{context}, eidx, avar::Ptr{abstract_var}, coeffs::Vector{Cdouble})
	res = ccall((:rhp_equ_addlin, libreshop), Cint,
		    (Ptr{context}, Cint, Ptr{abstract_var}, Ref{Cdouble}),
			ctx, eidx, avar, coeffs)
	res != 0 && error("return code $res from ReSHOP")
end

function rhp_equ_add_linear_chk(ctx::Ptr{context}, eidx, avar::Ptr{abstract_var}, coeffs::Vector{Cdouble})
	res = ccall((:rhp_equ_addlinchk, libreshop), Cint,
		    (Ptr{context}, Cint, Ptr{abstract_var}, Ref{Cdouble}),
			ctx, eidx, avar, coeffs)
	res != 0 && error("return code $res from ReSHOP")
end

function ctx_alloc()
	ctx = ccall((:ctx_alloc, libreshop), Ptr{context}, (Cuint,), 2)
	if ctx == Ptr{context}(C_NULL)
		error("Could not allocate a context")
	end
	return ctx
end
function ctx_create(n, m)
  ctx = ctx_alloc()
	res = ccall((:model_reserve_eqns, libreshop), Cint, (Ptr{context}, Ptr{Cvoid}, Cuint), ctx, C_NULL, m);
	res != 0 && error("return code $res from ReSHOP")
	res = ccall((:model_reserve_vars, libreshop), Cint, (Ptr{context}, Ptr{Cvoid}, Cuint), ctx, C_NULL, n);
	res != 0 && error("return code $res from ReSHOP")
	res = ccall((:ctx_resize, libreshop), Cint, (Ptr{context}, Cuint, Cuint), ctx, n, m)
	res != 0 && error("return code $res from ReSHOP")
	return ctx
end

function ctx_getvar(ctx, idx)
	return ccall((:ctx_getvar, libreshop), Ptr{reshop_var}, (Ptr{context}, Cint), ctx, idx)
end

function ctx_m(ctx)
	return ccall((:ctx_m, libreshop), Cuint, (Ptr{context},), ctx)
end

function ctx_n(ctx)
	return ccall((:ctx_n, libreshop), Cuint, (Ptr{context},), ctx)
end

function ctx_numvar(ctx::Ptr{context})
	if ctx == Ptr{context}(C_NULL)
		return 0
	end
	res = ccall((:model_total_n, libreshop), Csize_t, (Ptr{context},), ctx)
	return res
end

function ctx_numequ(ctx::Ptr{context})
	if ctx == Ptr{context}(C_NULL)
		return 0
	end
	res = ccall((:model_total_m, libreshop), Csize_t, (Ptr{context},), ctx)
	@chk_index_or_size(res)
	return res
end

function ctx_setequname(ctx, eidx, name::String)
	res = ccall((:myo_set_equname, libreshop), Cint, (Ptr{context}, Cint, Cstring), ctx, eidx, name)
	res != 0 && error("return code $res from ReSHOP")
end

function ctx_setvarnames(ctx, names::Vector{String})
	res = ccall((:myo_set_varnames, libreshop), Cint, (Ptr{context}, Ptr{Ptr{Cchar}}, Cuint), ctx, names, length(names))
	res != 0 && error("return code $res from ReSHOP")
end

function ctx_setvarname(ctx, vidx, name::String)
	res = ccall((:myo_set_varname, libreshop), Cint, (Ptr{context}, Cint, Cstring), ctx, vidx, name)
	res != 0 && error("return code $res from ReSHOP")
end

function ctx_getvarname(ctx, vidx)
	res = ccall((:myo_get_varname, libreshop), Cstring, (Ptr{context}, Cint), ctx, vidx)
	res != C_NULL || return ""
	return unsafe_string(res)
end

function hack_last_vidx(ctx::Ptr{context})
	res = ccall((:model_total_n, libreshop), Csize_t, (Ptr{context},), ctx) - 1
	@chk_index_or_size(res)
	return res
end

function hack_exportempinfo(ctx_gms::Ptr{context}, mdl::Ptr{reshop_model})
	res = ccall((:hack_exportempinfo, libreshop), Cint, (Ptr{context}, Ptr{reshop_model}), ctx_gms, mdl)
	res != 0 && error("return code $res from ReSHOP")
end

function hack_solver_log()
	ccall((:hack_solver_log, libreshop), Cvoid, ())
end

function ctx_setvarval(ctx::Ptr{context}, idx, val::Cdouble)
	res = ccall((:ctx_setvarlone, libreshop), Cint, (Ptr{context}, Cint, Cdouble), ctx, idx, val)
	res != 0 && error("return code $res from ReSHOP")
end

function ctx_setvarlb(ctx::Ptr{context}, idx, val::Cdouble)
	res = ccall((:ctx_setvarlb, libreshop), Cint, (Ptr{context}, Cint, Cdouble), ctx, idx, val)
	res != 0 && error("return code $res from ReSHOP")
end

function ctx_setvarub(ctx::Ptr{context}, idx, val::Cdouble)
	res = ccall((:ctx_setvarub, libreshop), Cint, (Ptr{context}, Cint, Cdouble), ctx, idx, val)
	res != 0 && error("return code $res from ReSHOP")
end

function ctx_setvarfx(ctx::Ptr{context}, idx, val::Cdouble)
	ctx_setvarub(ctx, idx, val)
	ctx_setvarlb(ctx, idx, val)
end

function ctx_getvarbounds(ctx::Ptr{context}, idx)
	lb = Ref{Cdouble}(NaN)
	ub = Ref{Cdouble}(NaN)
	res = ccall((:ctx_getvarbounds, libreshop), Cint, (Ptr{context}, Cint, Ref{Cdouble}, Ref{Cdouble}), ctx, idx, lb, ub)
	res != 0 && error("return code $res from ReSHOP")
	return (lb.x, ub.x)
end

function ctx_getvarbstat(ctx::Ptr{context}, idx)
	val = Ref{Cint}(-1)
	res = ccall((:ctx_getvarstatone, libreshop), Cint, (Ptr{context}, Cint, Ref{Cint}), ctx, idx, val)
	res != 0 && error("return code $res from ReSHOP")
	return val.x
end

function ctx_getvarmult(ctx::Ptr{context}, idx)
	val = Ref{Cdouble}(NaN)
	res = ccall((:ctx_getvarmone, libreshop), Cint, (Ptr{context}, Cint, Ref{Cdouble}), ctx, idx, val)
	res != 0 && error("return code $res from ReSHOP")
	return val.x
end

#function ctx_getvarname(ctx::Ptr{context}, idx)
#	len = 1024
#	buf = Vector{UInt8}(undef, len)
#	res = ccall((:ctx_getvarname, libreshop), Cint, (Ptr{context}, Cint, Ptr{UInt8}, Cuint), ctx, idx, buf, len)
#	res != 0 && error("return code $res from ReSHOP")
#	e = findfirst(isequal(0), buf)
#	return String(buf[1:e])
#end

function ctx_getvarbyname(ctx::Ptr{context}, name::String)
	val = Ref{Cint}(-1)
	res = ccall((:ctx_getvarbyname, libreshop), Cint, (Ptr{context}, Cstring, Ptr{Cint}), ctx, name, val)
	if res == 5
	  return -2
	end
	return val.x
end

function ctx_getequbyname(ctx::Ptr{context}, name::String)
	val = Ref{Cint}(-1)
	res = ccall((:ctx_getequbyname, libreshop), Cint, (Ptr{context}, Cstring, Ptr{Cint}), ctx, name, val)
	if res == 5
	  return -2
	end
	return val.x
end

function reshop_get_equtype(ctx::Ptr{context}, idx)
        typ = Ref{Cuint}(1000)
	cone = Ref{Cuint}(1000)
	res = ccall((:ctx_getequtype, libreshop), Cint, (Ptr{context}, Cint, Ref{Cuint}, Ref{Cuint}), ctx, idx, typ, cone)
	res != 0 && error("return code $res from ReSHOP")
	return (typ.x, cone.x)
end

function reshop_getvartype(ctx::Ptr{context}, idx)
	typ = Ref{Cuint}(1000)
	res = ccall((:ctx_getvartype, libreshop), Cint, (Ptr{context}, Cint, Ptr{Cuint}), ctx, idx, typ)
	res != 0 && error("return code $res from ReSHOP")
	return typ.x
end

function ctx_getvarval(ctx::Ptr{context}, idx)
	val = Ref{Cdouble}(NaN)
	res = ccall((:ctx_getvarlone, libreshop), Cint, (Ptr{context}, Cint, Ref{Cdouble}), ctx, idx, val)
	res != 0 && error("return code $res from ReSHOP")
	return val.x
end

function ctx_getequ(ctx, idx)
	equ = ccall((:ctx_getequ, libreshop), Ptr{reshop_equ}, (Ptr{context}, Cint), ctx, idx)
end

function ctx_getequval(ctx, idx)
	val = Ref{Cdouble}(NaN)
	res = ccall((:ctx_getequlone, libreshop), Cint, (Ptr{context}, Cint, Ref{Cdouble}), ctx, idx, val)
	res != 0 && error("return code $res from ReSHOP")
	return val.x
end

function ctx_getrhs(ctx::Ptr{context}, idx)
	val = Ref{Cdouble}(NaN)
	res = ccall((:ctx_getrhs, libreshop), Cint, (Ptr{context}, Cint, Ref{Cdouble}), ctx, idx, val)
	res != 0 && error("return code $res from ReSHOP")
	return val.x
end

function ctx_getequmult(ctx::Ptr{context}, idx)
	val = Ref{Cdouble}(NaN)
	res = ccall((:ctx_getequmone, libreshop), Cint, (Ptr{context}, Cint, Ref{Cdouble}), ctx, idx, val)
	res != 0 && error("return code $res from ReSHOP")
	return val.x
end

function ctx_getequbstat(ctx::Ptr{context}, idx)
	val = Ref{Cint}(-1)
	res = ccall((:ctx_getequstatone, libreshop), Cint, (Ptr{context}, Cint, Ref{Cint}), ctx, idx, val)
	res != 0 && error("return code $res from ReSHOP")
	return val.x
end

function ctx_getequname(ctx::Ptr{context}, eidx)
	res = ccall((:myo_get_equname, libreshop), Cstring, (Ptr{context}, Cint), ctx, eidx)
	res != C_NULL || return ""
	return unsafe_string(res)
end

function ctx_getequbounds(ctx::Ptr{context}, idx)
	lb = Ref{Cdouble}(NaN)
	ub = Ref{Cdouble}(NaN)
	res = ccall((:ctx_getequbounds, libreshop), Cint, (Ptr{context}, Cint, Ref{Cdouble}, Ref{Cdouble}), ctx, idx, lb, ub)
	res != 0 && error("return code $res from ReSHOP")
	return (lb.x, ub.x)
end

function emp_init(mdl::Ptr{reshop_model})
	res = ccall((:reshop_alloc_emp, libreshop), Cint, (Ptr{reshop_model},), mdl)
	res != 0 && error("return code $res from ReSHOP")
end

function emp_create_equil(max_mp)
	return ccall((:mp_equil_alloc, libreshop), Ptr{equil}, (Cuint,), max_mp)
end

function emp_add_mp_mp(mp_parent, mp)
	res = ccall((:mathprgm_addmp, libreshop), Cint, (Ptr{mathprgm}, Ptr{mathprgm}), mp_parent, mp)
	res != 0 && error("return code $res from ReSHOP")
end

function emp_add_mp_equil(mp_parent, mpe)
	res = ccall((:mathprgm_addequil, libreshop), Cint, (Ptr{mathprgm}, Ptr{equil}), mp_parent, mpe)
	res != 0 && error("return code $res from ReSHOP")
end

function emp_add_equil(mdl::Ptr{reshop_model})
	mpe = Ref{Ptr{equil}}(C_NULL)
	res = ccall((:reshop_add_equil, libreshop), Cint, (Ptr{reshop_model}, Ref{Ptr{equil}}), mdl, mpe)
	res != 0 && error("return code $res from ReSHOP")

	return mpe.x
end

function emp_equil_add(mpe, mp)
	res = ccall((:mp_equil_add, libreshop), Cint, (Ptr{equil}, Ptr{mathprgm}), mpe, mp)
	res != 0 && error("return code $res from ReSHOP")
end

function emp_mp_ensure(mdl::Ptr{reshop_model}, nb)
	res = ccall((:reshop_ensure_mp, libreshop), Cint, (Ptr{reshop_model}, Cuint), mdl, nb)
	res != 0 && error("return code $res from ReSHOP")
end

function emp_mp_alloc(mdl::Ptr{reshop_model})
	mp = ccall((:mathprgm_alloc, libreshop), Ptr{mathprgm}, (Ptr{reshop_model},), mdl)
	return mp
end

function emp_mp_start(mp, typ)
	res = ccall((:mathprgm_addstart, libreshop), Cint, (Ptr{mathprgm}, Cuint), mp, typ)
	res != 0 && error("return code $res from ReSHOP")
end

function emp_mp_objdir(mp, sense)
	res = ccall((:mathprgm_addobjdir, libreshop), Cint, (Ptr{mathprgm}, Cint), mp, sense_to_reshop[sense])
	res != 0 && error("return code $res from ReSHOP")
end

function emp_mp_objequ(mp, idx)
	res = ccall((:mathprgm_addobjequ, libreshop), Cint, (Ptr{mathprgm}, Cint), mp, idx)
	res != 0 && error("return code $res from ReSHOP")
end

function emp_mp_objvar(mp, idx)
	res = ccall((:mathprgm_addobjvar, libreshop), Cint, (Ptr{mathprgm}, Cint), mp, idx)
	res != 0 && error("return code $res from ReSHOP")
end

function emp_mp_var(mp, idx)
	res = ccall((:mathprgm_addvar, libreshop), Cint, (Ptr{mathprgm}, Cint), mp, idx)
	res != 0 && error("return code $res from ReSHOP")
end

function emp_mp_constraint(mp, idx)
	res = ccall((:mathprgm_addconstraint, libreshop), Cint, (Ptr{mathprgm}, Cint), mp, idx)
	res != 0 && error("return code $res from ReSHOP")
end

function emp_mp_vipair(mp, eidx, vidx)
	res = ccall((:mathprgm_addvipair, libreshop), Cint, (Ptr{mathprgm}, Cint, Cint), mp, eidx, vidx)
	res != 0 && error("return code $res from ReSHOP")
end

function emp_mp_getobjvar(mp)
	return ccall((:mathprgm_getobjvar, libreshop), Cint, (Ptr{mathprgm},), mp)
end

function emp_mp_getobjdir(mp)
	return ccall((:mathprgm_getobjvar, libreshop), Cint, (Ptr{mathprgm},), mp)
end

function emp_mp_print(mp, ctx)
	ccall((:mathprgm_print, libreshop), Cint, (Ptr{mathprgm}, Ptr{context}), mp, ctx)
end

function reshop_postprocess(mdl_solver::Ptr{reshop_model})
	res = ccall((:rhp_postprocess, libreshop), Cint, (Ptr{reshop_model},), mdl_solver)
	res != 0 && error("return code $res from ReSHOP")
end

function emp_set_root(mdl::Ptr{reshop_model}, mpe::Ptr{equil})
	res = ccall((:rhp_emproot_setmpe, libreshop), Cint, (Ptr{reshop_model}, Ptr{equil}), mdl, mpe)
	res != 0 && error("return code $res from ReSHOP")
end

function emp_set_root(mdl::Ptr{reshop_model}, mp::Ptr{mathprgm})
	res = ccall((:rhp_emproot_setmp, libreshop), Cint, (Ptr{reshop_model}, Ptr{mathprgm}), mdl, mp)
	res != 0 && error("return code $res from ReSHOP")
end

function emp_ovf(mdl, name, ovf_vidx, args)
	ovf_def = Ref{Ptr{ovf_definition}}(C_NULL)
	args_var = reshop_avar(length(args), args)
	res = ccall((:rhp_ovf_add, libreshop), Cint, (Ptr{reshop_model}, Cstring, Cint, Ptr{abstract_var}, Ref{Ptr{ovf_definition}}),
							                              mdl, name, ovf_vidx, args_var, ovf_def)
	res != 0 && error("return code $res from ReSHOP")
	reshop_avar_free(args_var)
	return ovf_def.x
end

function emp_ovf_param(ovf_def, param_name, scalar::Number)
	res = ccall((:rhp_ovf_param_add_scalar, libreshop), Cint, (Ptr{ovf_definition}, Cstring, Cdouble), ovf_def, param_name, scalar)
	res != 0 && error("return code $res from ReSHOP")
end

function emp_ovf_param(ovf_def, param_name, arr::Vector)
  arrC = Vector{Cdouble}(arr)
	res = ccall((:rhp_ovf_param_add_vector, libreshop), Cint, (Ptr{ovf_definition}, Cstring, Cuint, Ptr{Cdouble}), ovf_def, param_name, length(arr), arrC)
	res != 0 && error("return code $res from ReSHOP")
end

function emp_ovf_check(ovf_def)
	res = ccall((:rhp_ovf_check, libreshop), Cint, (Ptr{ovf_definition},), ovf_def)
end

function equtree_var(ctx, tree, node, idx, coeff)
	res = ccall((:equtree_var, libreshop), Cint, (Ptr{context}, Ptr{equtree}, Ref{Ref{Ptr{equnode}}}, Cint, Cdouble),
		ctx,
		tree,
		node,
		idx,
		coeff)
	res != 0 && error("return code $res from ReSHOP")
end

function equtree_cst(ctx, tree, node, value)
	res = ccall((:equtree_cst, libreshop), Cint, (Ptr{context}, Ptr{equtree}, Ref{Ref{Ptr{equnode}}}, Cdouble),
		ctx,
		tree,
		node,
		value)
	res != 0 && error("return code $res from ReSHOP")
end

function equtree_arithm(tree, node, opcode, nb)
	res = ccall((:equtree_arithm, libreshop), Cint, (Ptr{equtree}, Ref{Ref{Ptr{equnode}}}, Cuint, Cuint),
		tree,
		node,
		opcode,
		nb)
	res != 0 && error("return code $res from ReSHOP")
end

function equtree_call(ctx, tree, node, fndata)
	res = ccall((:equtree_call, libreshop), Cint, (Ptr{context}, Ptr{equtree}, Ref{Ref{Ptr{equnode}}}, Cuint, Cuint),
		ctx,
		tree,
		node,
		fndata[1],
		fndata[2]
		)
	res != 0 && error("return code $res from ReSHOP")
end

function equtree_get_root_addr(tree::Ptr{equtree}, node::Ref{Ref{Ptr{equnode}}})
	res = ccall((:equtree_get_root_addr, libreshop), Cint, (Ptr{equtree}, Ref{Ref{Ptr{equnode}}}),
		tree,
		node)
	res != 0 && error("return code $res from ReSHOP")
end

function equtree_umin(ctx, tree, node)
	res = ccall((:equtree_umin, libreshop), Cint, (Ptr{equtree}, Ref{Ref{Ptr{equnode}}}),
		tree,
		node)
	res != 0 && error("return code $res from ReSHOP")
end

function equnode_get_child_addr(node::Ptr{equnode}, i::Int)
	child = Ref{Ref{Ptr{equnode}}}(C_NULL)
	res = ccall((:equnode_get_child_addr, libreshop), Cint, (Ptr{equnode}, Ref{Ref{Ptr{equnode}}}, Cuint),
		node,
		child,
		i)
	res != 0 && error("return code $res from ReSHOP")
	return child
end

function equnode_deref(node::Ref{Ref{Ptr{equnode}}})
	return ccall((:p2deref, libreshop), Ptr{equnode}, (Ref{Ref{Ptr{equnode}}},), node)
end

function ctx_get_solvername(ctx::Ptr{context})
	str = "?"^256
	res = ccall((:gams_getsolverstr, libreshop), Cint, (Ptr{context}, Cstring, Cuint), ctx, str, length(str))
	res != 0 && error("return code $res from ReSHOP")
	return str
end

function rhp_add_free_var(ctx::Ptr{context}, nb, avar::Ptr{abstract_var}=Ptr{abstract_var}(C_NULL))
	res = ccall((:model_add_free_vars, libreshop), Cint, (Ptr{context}, Cuint, Ptr{abstract_var}), ctx, nb, avar)
	res != 0 && error("return code $res from ReSHOP")
end

function rhp_add_neg_var(ctx::Ptr{context}, nb, avar::Ptr{abstract_var}=Ptr{abstract_var}(C_NULL))
	res = ccall((:model_add_neg_vars, libreshop), Cint, (Ptr{context}, Cuint, Ptr{abstract_var}), ctx, nb, avar)
	res != 0 && error("return code $res from ReSHOP")
end

function rhp_add_pos_var(ctx::Ptr{context}, nb, avar::Ptr{abstract_var}=Ptr{abstract_var}(C_NULL))
	res = ccall((:model_add_pos_vars, libreshop), Cint, (Ptr{context}, Cuint, Ptr{abstract_var}), ctx, nb, avar)
	res != 0 && error("return code $res from ReSHOP")
end

function reshop_alloc(ctx::Ptr{context})
	return ccall((:reshop_alloc, libreshop), Ptr{reshop_model}, (Ptr{context},), ctx)
end

function reshop_avar()
  return ccall((:avar_alloc, libreshop), Ptr{abstract_var}, ())
end

function reshop_avar(size, indices)
	indicesC = Vector{Cint}(indices)
	return ccall((:avar_alloc_list, libreshop), Ptr{abstract_var}, (Cuint, Ptr{Cint}), size, indicesC)
end

function rhp_avar_set(avar::Ptr{abstract_var}, indices::Integer)
	indicesC = Vector{Cint}([indices])
	res = ccall((:avar_set_list, libreshop), Cint, (Ptr{abstract_var}, Cuint, Ptr{Cint}), avar, length(indicesC), indicesC)
	res != 0 && error("return code $res from ReSHOP")
end

function rhp_avar_set(avar::Ptr{abstract_var}, indices::AbstractVector{Int32})
	indicesC = Vector{Cint}(indices)
	res = ccall((:avar_set_list, libreshop), Cint, (Ptr{abstract_var}, Cuint, Ptr{Cint}), avar, length(indices), indicesC)
	res != 0 && error("return code $res from ReSHOP")
end

function reshop_avar_free(avar::Ptr{abstract_var})
	return ccall((:avar_free, libreshop), Cvoid, (Ptr{abstract_var},), avar)
end

function reshop_avar_get(avar::Ptr{abstract_var}, i::UInt32)
	vidx = Ref{Cint}(-1)
	res = ccall((:avar_get, libreshop), Cint, (Ptr{abstract_var}, Cuint, Ref{Cint}), avar, i, vidx)
	res != 0 && error("return code $res from ReSHOP")
	return vidx.x
end

function reshop_avar_size(avar::Ptr{abstract_var})
	res = ccall((:avar_size, libreshop), Csize_t, (Ptr{abstract_var},), avar)
	@chk_index_or_size(res)
	return res
end

function reshop_decl_eqn(ctx::Ptr{context})
	minn = Ref{Cint}(-1)
	res = ccall((:rhp_add_equation, libreshop), Cint, (Ptr{context}, Ref{Cint}), ctx, minn)
	res != 0 && error("return code $res from ReSHOP")
end

function rhp_equ_add_lin_tree(ctx::Ptr{context}, eidx, qvals, avar::Ptr{abstract_var}, coeff)
	equ = ctx_getequ(ctx, eidx)
	res = ccall((:equ_add_lin_tree, libreshop), Cint, (Ptr{context}, Ptr{reshop_equ}, Ref{Cdouble}, Ptr{abstract_var}, Cdouble),
							ctx, equ, qvals, avar, coeff)
	res != 0 && error("return code $res from ReSHOP")
end

function rhp_equ_add_quadratic(ctx::Ptr{context}, eidx, mat::Ptr{reshop_sp_matrix}, avar::Ptr{abstract_var}, coeff)
	equ = ctx_getequ(ctx, eidx)
	res = ccall((:equ_add_quadratic, libreshop), Cint, (Ptr{context}, Ptr{reshop_equ}, Ptr{reshop_sp_matrix}, Ptr{abstract_var}, Cdouble),
							ctx, equ, mat, avar, coeff)
	res != 0 && error("return code $res from ReSHOP")
end

function rhp_equ_add_quadratic(ctx::Ptr{context}, eidx, vidx1, vidx2, coeffs)
	vidx1C = Vector{Cint}(vidx1)
	vidx2C = Vector{Cint}(vidx2)
	@assert length(vidx1) == length(vidx2) == length(coeffs)
	res = ccall((:rhp_equ_addquadabsolute, libreshop), Cint,
		    (Ptr{context}, Cint, Csize_t, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}, Cdouble),
		    ctx, eidx, length(vidx1C), vidx1C, vidx2C, coeffs, 1.)
	res != 0 && error("return code $res from ReSHOP")
end

function reshop_free(mdl::Ptr{reshop_model})
	return ccall((:reshop_free, libreshop), Cvoid, (Ref{Ptr{reshop_model}},), mdl)
end

function reshop_mat_coo(ridx, cidx, vals)
	# beaware of dragons! --xhub
	ridxC = Vector{Cint}(ridx)
	cidxC = Vector{Cint}(cidx)
	mat = ccall((:rhp_mat_triplet, libreshop), Ptr{reshop_sp_matrix}, (Cuint, Cuint, Csize_t, Ptr{Cint}, Ptr{Cint}, Ref{Cdouble}),
							0, 0, length(vals), ridxC, cidxC, vals)
	return mat
end

function reshop_mat_free(mat)
	ccall((:rhp_mat_free, libreshop), Cvoid, (Ptr{reshop_sp_matrix},), mat)
end

function rhp_set_objeqn(ctx::Ptr{context}, eidx)
	res = ccall((:model_setobjequ, libreshop), Cint, (Ptr{context}, Cint), ctx, eidx)
	res != 0 && error("return code $res from ReSHOP")
end

function rhp_set_objvar(ctx::Ptr{context}, vidx)
	res = ccall((:ctx_setobjvar, libreshop), Cint, (Ptr{context}, Cint), ctx, vidx)
	res != 0 && error("return code $res from ReSHOP")
end

function ctx_getobjequ(ctx::Ptr{context})
	tmpCint = Ref{Cint}(0)
	res = ccall((:ctx_getobjequ, libreshop), Cint, (Ptr{context}, Ptr{Cint}), ctx, tmpCint)
	res != 0 && error("return code $res from ReSHOP")
	return tmpCint.x
end

function ctx_getobjvar(ctx::Ptr{context})
	tmpCint = Ref{Cint}(0)
	res = ccall((:ctx_getobjvar, libreshop), Cint, (Ptr{context}, Ptr{Cint}), ctx, tmpCint)
	res != 0 && error("return code $res from ReSHOP")
	return tmpCint.x
end

function reshop_set_rhs(ctx::Ptr{context}, idx, val)
	res = ccall((:ctx_setrhs, libreshop), Cint, (Ptr{context}, Cint, Cdouble), ctx, idx, val)
	res != 0 && error("return code $res from ReSHOP")
end

function reshop_set_equtype(ctx::Ptr{context}, idx, cone)
	res = ccall((:ctx_setequtype, libreshop), Cint, (Ptr{context}, Cint, Cuint, Cuint), ctx, idx, 2, cone)
	res != 0 && error("return code $res from ReSHOP")
end

function reshop_set_vartype(ctx::Ptr{context}, idx, typ)
	res = ccall((:ctx_setvartype, libreshop), Cint, (Ptr{context}, Cint, Cuint), ctx, idx, typ)
	res != 0 && error("return code $res from ReSHOP")
end

function reshop_set_vartype(ctx::Ptr{context}, avar::Ptr{abstract_var}, ::MOI.SOS1{Float64}, weights::AbstractVector{Cdouble})
  res = ccall((:rhp_set_var_sos1, libreshop), Cint, (Ptr{context}, Ptr{abstract_var}, Ptr{Cdouble}), ctx, avar, weights)
	res != 0 && error("return code $res from ReSHOP")
end

function reshop_set_vartype(ctx::Ptr{context}, avar::Ptr{abstract_var}, ::MOI.SOS2{Float64}, weights::AbstractVector{Cdouble})
  res = ccall((:rhp_set_var_sos2, libreshop), Cint, (Ptr{context}, Ptr{abstract_var}, Ptr{Cdouble}), ctx, avar, weights)
	res != 0 && error("return code $res from ReSHOP")
end

function reshop_get_treedata(ctx, eidx::Int)
	tree = ccall((:myo_getequtree, libreshop), Ptr{equtree}, (Ptr{context}, Cint), ctx, eidx)
	node = Ref{Ref{Ptr{equnode}}}(C_NULL)
	res = equtree_get_root_addr(tree, node)
	res != 0 && error("return code $res from ReSHOP")
	return (tree, node)
end

function rhp_add_equ(ctx::Ptr{context})
	return rhp_add_equs(ctx, 1)
end

function rhp_add_equs(ctx::Ptr{context}, nb)
	eidx_start = Int(ctx_numequ(ctx))
	res = ccall((:myo_addinit_equs, libreshop), Cint, (Ptr{context}, Cuint), ctx, nb)
	res != 0 && error("return code $res from ReSHOP")
	return eidx_start
end

function rhp_add_var(ctx::Ptr{context}, avar::Ptr{abstract_var})
	return rhp_add_var(ctx, 1, avar)
end

function rhp_add_var(ctx::Ptr{context}, nb, avar::Ptr{abstract_var})
	res = ccall((:rhp_add_vars, libreshop), Cint, (Ptr{context}, Cuint, Ptr{abstract_var}), ctx, nb, avar)
	res != 0 && error("return code $res from ReSHOP")
end

function rhp_get_solvestat(ctx::Ptr{context})
	tmpCint = Ref{Cint}(0)
	res = ccall((:ctx_getsolvestat, libreshop), Cint, (Ptr{context}, Ref{Cint}), ctx, tmpCint)
	res != 0 && error("return code $res from ReSHOP")
	return tmpCint.x
end

function rhp_get_modelstat(ctx::Ptr{context})
	tmpCint = Ref{Cint}(0)
	res = ccall((:ctx_getmodelstat, libreshop), Cint, (Ptr{context}, Ref{Cint}), ctx, tmpCint)
	res != 0 && error("return code $res from ReSHOP")
	return tmpCint.x
end

function reshop_options_alloc()
	return ccall((:hack_options_alloc, libreshop), Ptr{reshop_options}, ())
end

function reshop_options_dealloc(o::Ptr{reshop_options})
	ccall((:hack_options_dealloc, libreshop), Cvoid, (Ptr{reshop_options},), o)
end

function reshop_option_get(opt::Ptr{reshop_options}, k::String, val::Ref{Cdouble})
	res = ccall((:option_get_d, libreshop), Cint, (Ptr{reshop_options}, Cstring, Ptr{Cdouble}), opt, k, val)
	res != 0 && error("return code $res from ReSHOP")
	return val.x
end

function reshop_option_get(opt::Ptr{reshop_options}, k::String, val::Ref{Cstring})
	res = ccall((:option_get_i, libreshop), Cint, (Ptr{reshop_options}, Cstring, Ptr{Cstring}), opt, k, val)
	res != 0 && error("return code $res from ReSHOP")
	return val.x
end

function reshop_option_get(opt::Ptr{reshop_options}, k::String, val::Ref{Cint})
	res = ccall((:option_get_i, libreshop), Cint, (Ptr{reshop_options}, Cstring, Ptr{Cint}), opt, k, val)
	res != 0 && error("return code $res from ReSHOP")
	return val.x
end

function reshop_option_get(opt::Ptr{reshop_options}, k::String, val::Ref{Cuint})
	res = ccall((:option_get_i, libreshop), Cint, (Ptr{reshop_options}, Cstring, Ptr{Cuint}), opt, k, val)
	res != 0 && error("return code $res from ReSHOP")
	return val.x
end

function reshop_option_get(opt::Ptr{reshop_options}, k::String)
	res = ccall((:option_get_type, libreshop), Cint, (Ptr{reshop_options}, Cstring), opt, k)
	if haskey(reshop_option_type, res)
		fn, type = reshop_option_type[res]
		val = Ref{type}()
		reshop_option_get(opt, k, val)
		return val.x
	else
		error("unknown option named $k, ReSHOP code is $res")
	end
end

function reshop_option_set(opt::Ptr{reshop_options}, k::String, v::Bool)
	return ccall((:option_set_b, libreshop), Cint, (Ptr{reshop_options}, Cstring, Cint), opt, k, v)
end

function reshop_option_set(opt::Ptr{reshop_options}, k::String, v::Integer)
	return ccall((:option_set_i, libreshop), Cint, (Ptr{reshop_options}, Cstring, Cint), opt, k, v)
end

function reshop_option_set(opt::Ptr{reshop_options}, k::String, v::Cdouble)
	return ccall((:option_set_d, libreshop), Cint, (Ptr{reshop_options}, Cstring, Cdouble), opt, k, v)
end

function reshop_option_set(opt::Ptr{reshop_options}, k::String, v::String)
	return ccall((:option_set_str, libreshop), Cint, (Ptr{reshop_options}, Cstring, Cstring), opt, k, v)
end

function reshop_set_modeltype(ctx::Ptr{context}, model_type)
    res = ccall((:ctx_setmodeltype, libreshop), Cint, (Ptr{context}, Cint), ctx, model_type)
    res != 0 && error("return code $res from ReSHOP")
end

function reshop_set_objsense(ctx::Ptr{context}, sense)
    res = ccall((:ctx_setobjsense, libreshop), Cint, (Ptr{context}, Cint), ctx, sense_to_reshop[sense])
    res != 0 && error("return code $res from ReSHOP")
end

function rhp_set_option(ctx::Ptr{context}, k::String, v::Bool)
	res = ccall((:rhp_set_option_b, libreshop), Cint, (Ptr{context}, Cstring, Cuchar), ctx, k, v)
	res != 0 && error("return code $res from ReSHOP")
end

function rhp_set_option(ctx::Ptr{context}, k::String, v::Cint)
	res = ccall((:rhp_set_option_i, libreshop), Cint, (Ptr{context}, Cstring, Cint), ctx, k, v)
	res != 0 && error("return code $res from ReSHOP")
end

function rhp_set_option(ctx::Ptr{context}, k::String, v::Cdouble)
	res = ccall((:rhp_set_option_d, libreshop), Cint, (Ptr{context}, Cstring, Cdouble), ctx, k, v)
	res != 0 && error("return code $res from ReSHOP")
end

function rhp_set_option(ctx::Ptr{context}, k::String, v::Union{Cstring, String})
	res = ccall((:rhp_set_option_s, libreshop), Cint, (Ptr{context}, Cstring, Cstring), ctx, k, v)
	res != 0 && error("return code $res from ReSHOP")
end

function ctx_dealloc(o::Ptr{context})
	ccall((:ctx_dealloc,  libreshop), Cvoid, (Ptr{context},), o)
end

function rhp_delete_var(ctx::Ptr{context}, vidx)
	res = ccall((:rhp_delete_var,  libreshop), Cint, (Ptr{context}, Cint), ctx, vidx)
	res != 0 && error("return code $res from ReSHOP")
end

function rhp_delete_eqn(ctx::Ptr{context}, eidx)
	res = ccall((:rhp_delete_eqn,  libreshop), Cint, (Ptr{context}, Cint), ctx, eidx)
	res != 0 && error("return code $res from ReSHOP")
end

function rhp_get_nb_bin_var(ctx::Ptr{context})
	res = ccall((:rhp_get_nb_var_bin, libreshop), Csize_t, (Ptr{context},), ctx)
	@chk_index_or_size(res)
	return res
end

function rhp_get_nb_int_var(ctx::Ptr{context})
	res = ccall((:rhp_get_nb_var_int, libreshop), Csize_t, (Ptr{context},), ctx)
	@chk_index_or_size(res)
	return res
end

function rhp_get_nb_var(ctx::Ptr{context}, ::Type{MOI.GreaterThan{Float64}})
	res = ccall((:rhp_get_nb_var_lb, libreshop), Csize_t, (Ptr{context},), ctx)
	@chk_index_or_size(res)
	return res
end

function rhp_get_nb_var(ctx::Ptr{context}, ::Type{MOI.LessThan{Float64}})
	res = ccall((:rhp_get_nb_var_ub, libreshop), Csize_t, (Ptr{context},), ctx)
	@chk_index_or_size(res)
	return res
end

function rhp_get_nb_var(ctx::Ptr{context}, ::Type{MOI.EqualTo{Float64}})
	res = ccall((:rhp_get_nb_var_fx, libreshop), Csize_t, (Ptr{context},), ctx)
	@chk_index_or_size(res)
	return res
end

function rhp_get_nb_var(ctx::Ptr{context}, ::Type{MOI.Interval{Float64}})
	res = ccall((:rhp_get_nb_var_interval, libreshop), Csize_t, (Ptr{context},), ctx)
	@chk_index_or_size(res)
	return res
end

function rhp_get_nb_var(ctx::Ptr{context}, ::Type{MOI.SOS1{Float64}})
  res = ccall((:rhp_get_nb_var_sos1, libreshop), Csize_t, (Ptr{context},), ctx)
  @chk_index_or_size(res)
  return res
end

function rhp_get_nb_var(ctx::Ptr{context}, ::Type{MOI.SOS2{Float64}})
  res = ccall((:rhp_get_nb_var_sos2, libreshop), Csize_t, (Ptr{context},), ctx)
  @chk_index_or_size(res)
  return res
end

function rhp_get_nb_lequ(ctx::Ptr{context}, ::Type{MOI.GreaterThan{Float64}})
	res = ccall((:rhp_get_nb_lequ_ge, libreshop), Csize_t, (Ptr{context},), ctx)
	@chk_index_or_size(res)
	return res
end

function rhp_get_nb_lequ(ctx::Ptr{context}, ::Type{MOI.LessThan{Float64}})
	res = ccall((:rhp_get_nb_lequ_le, libreshop), Csize_t, (Ptr{context},), ctx)
	@chk_index_or_size(res)
	return res
end

function rhp_get_nb_lequ(ctx::Ptr{context}, ::Type{MOI.EqualTo{Float64}})
	res = ccall((:rhp_get_nb_lequ_eq, libreshop), Csize_t, (Ptr{context},), ctx)
	@chk_index_or_size(res)
	return res
end

function rhp_get_sos_group(ctx::Ptr{context}, vidx, ::Type{MOI.SOS1{Float64}})
	grps = Ref{Ref{Cuint}}(C_NULL)
	res = ccall((:rhp_get_var_sos1, libreshop), Cint, (Ptr{context}, Cint, Ref{Ptr{Cuint}}), ctx, vidx, grps)
	res != 0 && error("return code $res from ReSHOP")
	return grps.x
end

function rhp_get_sos_group(ctx::Ptr{context}, vidx, ::Type{MOI.SOS2{Float64}})
	grps = Ref{Ptr{Cuint}}(C_NULL)
	res = ccall((:rhp_get_var_sos2, libreshop), Cint, (Ptr{context}, Cint, Ref{Ptr{Cuint}}), ctx, vidx, grps)
	res != 0 && error("return code $res from ReSHOP")
	return grps.x
end

function rhp_is_var_valid(ctx::Ptr{context}, vidx)
	return ccall((:rhp_is_var_valid, libreshop), Cint, (Ptr{context}, Cint,), ctx, vidx) == 1 ? true : false
end

function rhp_is_equ_valid(ctx::Ptr{context}, eidx)
	return ccall((:rhp_is_equ_valid, libreshop), Cint, (Ptr{context}, Cint,), ctx, eidx)  == 1 ? true : false
end


text:
""" Draw wiring diagrams using Compose.jl.
"""
module ComposeWiringDiagrams
export ComposePicture, to_composejl, layout_to_composejl

using Compat
using LinearAlgebra: dot
import Compose
const C = Compose

using ...WiringDiagrams
using ..WiringDiagramLayouts
using ..WiringDiagramLayouts: AbstractVector2D, Vector2D, BoxLayout, box_label,
  position, size, lower_corner, upper_corner, normal, tangent, wire_points

# Data types
############

""" A Compose context together with a given width and height.

We need this type because contexts have no notion of size or aspect ratio, but
wiring diagram layouts have fixed aspect ratios.
"""
struct ComposePicture
  context::Compose.Context
  width::Compose.Measure
  height::Compose.Measure
end

function Base.show(io::IO, ::MIME"text/html", pic::ComposePicture)
  jsmode = C.default_jsmode
  C.draw(C.SVGJS(io, pic.width, pic.height, false, jsmode=jsmode), pic.context)
end
function Base.show(io::IO, ::MIME"image/svg+xml", pic::ComposePicture)
  C.draw(C.SVG(io, pic.width, pic.height, false), pic.context)
end

const ComposeProperties = AbstractVector{<:Compose.Property}

""" Internal data type for configurable options of Compose.jl wiring diagrams.
"""
Base.@kwdef struct ComposeOptions
  base_unit::Compose.Measure = 4*C.mm
  props::AbstractDict = Dict()
  box_renderer::Function = render_box
  background_color::Union{AbstractString,C.Colorant,Nothing} = nothing
  rounded_boxes::Bool = true
end

# Wiring diagrams
#################

""" Draw a wiring diagram in Compose.jl.
"""
function to_composejl(args...; kw...)::ComposePicture
  layout_kw = filter(p -> first(p) ∉ fieldnames(ComposeOptions), kw)
  diagram = layout_diagram(args...; layout_kw...)

  compose_kw = filter(p -> first(p) ∈ fieldnames(ComposeOptions), kw)
  layout_to_composejl(diagram; compose_kw...)
end

""" Draw a wiring diagram in Compose.jl using the given layout.
"""
function layout_to_composejl(diagram::WiringDiagram; props=Dict(), kw...)
  props = merge(default_props, props)
  opts = ComposeOptions(; props=props, kw...)
  context = layout_to_context(diagram, opts)
  if !isnothing(opts.background_color)
    context = C.compose(C.context(),
      context,
      (C.context(), C.rectangle(), C.fill(opts.background_color)))
  end
  ComposePicture(context, (opts.base_unit .* size(diagram))...)
end

function layout_to_context(diagram::WiringDiagram, opts::ComposeOptions)
  # The origin of the Compose.jl coordinate system is in the top-left corner,
  # while the origin of wiring diagram layout is at the diagram center, so we
  # need to translate the coordinates using the `UnitBox`.
  units = C.UnitBox(-size(diagram)/2..., size(diagram)...)
  C.compose(C.context(units=units, tag=:diagram),
    C.compose(C.context(tag=:boxes),
      map(boxes(diagram)) do box
        layout_to_context(box, opts)
      end...
    ),
    C.compose(C.context(tag=:wires),
      map(wires(diagram)) do wire
        layout_to_context(diagram, wire, opts)
      end...
    ),
  )
end

function layout_to_context(box::Box, opts::ComposeOptions)
  C.compose(C.context(lower_corner(box)..., size(box)...,
                      units=C.UnitBox(), tag=:box),
    opts.box_renderer(box.value, opts)
  )
end

function layout_to_context(diagram::WiringDiagram, wire::Wire, opts::ComposeOptions)
  src, tgt = wire.source, wire.target
  src_pos, tgt_pos = position(diagram, src), position(diagram, tgt)
  points = [ src_pos ]
  prev = (src_pos, normal(diagram, src))
  for point in wire_points(wire)
    p, v = position(point), tangent(point)
    append!(points, tangents_to_controls(prev..., -v, p))
    push!(points, p)
    prev = (p, v)
  end
  append!(points, tangents_to_controls(prev..., normal(diagram, tgt), tgt_pos))
  push!(points, tgt_pos)

  C.compose(C.context(tag=:wire),
    piecewise_curve(map(Tuple, points)),
    opts.props[:wire]...
  )
end

""" Control points for Bezier curve from endpoints and unit tangent vectors.
"""
function tangents_to_controls(p1::AbstractVector2D, v1::AbstractVector2D,
    v2::AbstractVector2D, p2::AbstractVector2D; weight::Float64=0.5)
  v = p2 - p1
  (p1 + weight * dot(v,v1) * v1, p2 - weight * dot(v,v2) * v2)
end

""" Draw an atomic box in Compose.jl.
"""
function render_box(layout::BoxLayout, opts::ComposeOptions)
  render_box(Val(layout.shape), layout, opts)
end
function render_box(::Val{:rectangle}, layout::BoxLayout, opts::ComposeOptions)
  form = opts.rounded_boxes ? rounded_rectangle() : C.rectangle()
  labeled_box(form, layout, opts)
end
render_box(::Val{:circle}, layout::BoxLayout, opts::ComposeOptions) =
  labeled_box(C.circle(), layout, opts)
render_box(::Val{:ellipse}, layout::BoxLayout, opts::ComposeOptions) =
  labeled_box(C.ellipse(), layout, opts)
render_box(::Val{:junction}, layout::BoxLayout, opts::ComposeOptions) =
  C.compose(C.context(), C.circle(), box_props(layout, opts)...)
render_box(::Val{:invisible}, ::BoxLayout, ::ComposeOptions) = C.context()

function labeled_box(form::C.ComposeNode, layout::BoxLayout, opts::ComposeOptions)
  labeled_form(form, box_label(MIME("text/plain"), layout.value),
    form_props=box_props(layout, opts), text_props=opts.props[:text])
end

""" Get Compose.jl properties for box.
"""
function box_props(layout::BoxLayout, opts::ComposeOptions)
  style = if isnothing(layout.style)
    layout.shape in (:junction, :invisible) ? layout.shape : :box
  else
    layout.style
  end
  opts.props[style]
end

# Compose.jl forms
##################

""" Draw a form with centered text label in Compose.
"""
function labeled_form(form::C.ComposeNode, label::String; rounded::Bool=true,
    form_props::ComposeProperties=[], text_props::ComposeProperties=[])
  C.compose(C.context(),
    (C.context(order=1), form, form_props...),
    (C.context(order=2),
     C.text(0.5, 0.5, label, C.hcenter, C.vcenter),
     text_props...),
  )
end

""" Draw a rounded rectangle in Compose.
"""
function rounded_rectangle(args...; corner_radius::Compose.Measure=C.mm)
  w, h, px = Compose.w, Compose.h, Compose.px
  r = corner_radius
  C.compose(C.context(args...),
    C.compose(C.context(order=1),
      # Layer 1: fill.
      # Since we cannot use a stroke here, we add a single pixel fudge factor to
      # prevent a small but visible gap between the different fill regions.
      C.rectangle(r-1px, 0h, 1w-2r+2px, 1h),
      C.rectangle(0w, r-1px, 1w, 1h-2r+2px),
      C.sector(r,    r,    r, π,    3π/2),
      C.sector(1w-r, r,    r, 3π/2, 0),
      C.sector(r,    1h-r, r, π/2,  π),
      C.sector(1w-r, 1h-r, r, 0,    π/2),
      C.stroke(nothing),
    ),
    C.compose(C.context(order=2),
      # Layer 2: stroke.
      C.line([
        [ (r, 0h), (1w-r, 0h) ],
        [ (0w, r), (0w, 1h-r) ],
        [ (r, 1h), (1w-r, 1h) ],
        [ (1w, r), (1w, 1h-r) ],
      ]),
      C.arc(r,    r,    r, π,    3π/2),
      C.arc(1w-r, r,    r, 3π/2, 0),
      C.arc(r,    1h-r, r, π/2,  π),
      C.arc(1w-r, 1h-r, r, 0,    π/2),
    ),
  )
end

""" Draw a continuous piecewise cubic Bezier curve in Compose.

The points are given by a vector [p1, c1, c2, p2, c3, c4, p3, ...] of length
3n+1 where n >=1 is the number of Bezier curves.
"""
function piecewise_curve(points::Vector{<:Tuple})
  @assert length(points) % 3 == 1
  C.compose(C.context(),
    (C.curve(points[i:i+3]...) for i in range(1, length(points)-1, step=3))...
  )
end

# Compose.jl properties
#######################

const default_props = Dict{Symbol,Vector{Compose.Property}}(
  :box => [ C.stroke("black"), C.fill(nothing) ],
  :junction => [ C.stroke("black") ],
  :variant_junction => [ C.stroke("black"), C.fill(nothing) ],
  :wire => [ C.stroke("black") ],
  :text => [ C.fill("black") ],
)

end

ext: jl


text:
macro intrinsic(ex)
    esc(intrinsic_m(ex))
end

macro intrinsic(n::Int, ex)
    esc(intrinsic_m(ex, n))
end

function intrinsic_m(ex, n::Int=1)
    if ex isa Symbol
        def = JLStruct(;name=Symbol(ex, :Gate), supertype=IntrinsicRoutine)
        binding_name = ex
        binding = :(Core.@__doc__ const $ex = $(def.name)())
    else
        name, args, kw, whereparams, rettype = split_function_head(ex)
        kw === nothing || error("cannot have kwargs in intrinsic operation")
        rettype === nothing || error("cannot have rettype in intrinsic operation")

        typevars=whereparams
        fields = JLField[]
        for each in args
            field = @match each begin
                :($fname::$type) => JLField(;name=fname, type)
                ::Symbol => JLField(;name=each)
                _ => error("cannot handle $each")
            end
            push!(fields, field)
        end
        def = JLStruct(;name, fields, supertype=IntrinsicRoutine, typevars)
        binding_name = name
        binding = nothing
    end

    quote
        export $(def.name)
        Core.@__doc__ $(codegen_ast(def))
        $(codegen_helpers(def, n, binding_name))
        $binding
    end
end

nqubits(::Type{T}) where {T <: Routine} = error("YaoHIR.nqubits is not defined for $T")
isintrinsic(::Type{T}) where T = false

function codegen_helpers(def::JLStruct, n::Int, name::Symbol)
    quote
        $YaoHIR.nqubits(::Type{<:$(def.name)}) = $n
        $YaoHIR.isintrinsic(::Type{<:$(def.name)}) = true
        $YaoHIR.routine_name(::Type{<:$(def.name)}) = $(QuoteNode(name))
        $MLStyle.@as_record $(def.name)
    end
end

@intrinsic X
@intrinsic Y
@intrinsic Z
@intrinsic H
@intrinsic S
@intrinsic T
@intrinsic SWAP
@intrinsic shift(θ::T) where {T}
@intrinsic Rx(θ::T) where {T}
@intrinsic Ry(θ::T) where {T}
@intrinsic Rz(θ::T) where {T}
@intrinsic UGate(α::T, β::T, γ::T) where {T}


text:
export get_historical_data, get_exchange_volume, get_symbol_volume, get_symbol_volume_single_exchange
"""
    get_historical_data(duration::String, fsym::String, tsym::String; tryConversion::Bool=true,
                                e::String="all", aggregate::Int=1, aggregatePredictableTimePeriods::Bool=true,
                                limit::Int=30, allData::Bool=false, toTs::String="", explainPath::Bool=false, extraParams::String="", sign::Bool=false)

Get open, high, low, close, volumefrom and volumeto from the daily/hourly/minutewise historical data.
The values are based on 00:00 GMT time.
If e=CCCAGG and tryConversion=true, it attempts conversion through BTC or ETH to determine the best possible path.
The conversion type and symbol used are appended per historical point.
If you want to get all the available historical data, you can use limit=2000 and keep going back in time using the toTs param.
You can then keep requesting batches using: &limit=2000&toTs={the earliest timestamp received}.

1. tryConversion: If set to false, it will try to get only direct trading values. This parameter is only valid for e=CCCAGG value [ Default - true]
2. fsym [Required]: The cryptocurrency symbol of interest [ Min length - 1] [ Max length - 30]
3. tsym [Required]: The currency symbol to convert into [ Min length - 1] [ Max length - 30]
4. e: The exchange to obtain data from [ Min length - 2] [ Max length - 30] [ Default - CCCAGG]
5. aggregate: Time period to aggregate the data over (for daily it's days, for hourly it's hours and for minute histo it's minutes) [ Min - 1] [ Max - 30] [ Default - 1]
6. aggregatePredictableTimePeriods: Only used when the aggregate param is also in use. If false it will aggregate based on the current time.If the param is false and the time you make the call is 1pm - 2pm, with aggregate 2, it will create the time slots: ... 9am, 11am, 1pm.If the param is false and the time you make the call is 2pm - 3pm, with aggregate 2, it will create the time slots: ... 10am, 12am, 2pm.If the param is true (default) and the time you make the call is 1pm - 2pm, with aggregate 2, it will create the time slots: ... 8am, 10am, 12pm.If the param is true (default) and the time you make the call is 2pm - 3pm, with aggregate 2, it will create the time slots: ... 10am, 12am, 2pm. [ Default - true]
7. limit: The number of data points to return. If limit * aggregate > 2000 we reduce the limit param on our side. So a limit of 1000 and an aggerate of 4 would only return 2000 (max points) / 4 (aggregation size) = 500 total points + current one so 501. [ Min - 1] [ Max - 2000] [ Default - 30]
8. allData: Returns all data (only available on histo day) [ Default - false]
9. toTs: Returns historical data before that timestamp. If you want to get all the available historical data, you can use limit=2000 and keep going back in time using the toTs param. You can then keep requesting batches using: &limit=2000&toTs={the earliest timestamp received}
10. explainPath: If set to true, each point calculated will return the available options it used to make the calculation path decision. This is intended for calculation verification purposes, please note that manually recalculating the returned data point values from this data may not match exactly, this is due to levels of caching in some circumstances. [ Default - false]
11. extraParams: The name of your application (we recommend you send it) [ Min length - 1] [ Max length - 2000] [ Default - NotAvailable]
12. sign: If set to true, the server will sign the requests (by default we don't sign them), this is useful for usage in smart contracts [ Default - false]
"""
function get_historical_data(duration::String, fsym::String, tsym::String; tryConversion::Bool=true,
                                e::String="all", aggregate::Int=1, aggregatePredictableTimePeriods::Bool=true,
                                limit::Int=30, allData::Bool=false, toTs::String="", explainPath::Bool=false, extraParams::String="", sign::Bool=false)
    if duration in ["day", "minute", "hour"]
        url = build_historical_data_url(duration, fsym=fsym, tsym=tsym, e=e, tryConversion=tryConversion,
                          aggregate=aggregate, aggregatePredictableTimePeriods=aggregatePredictableTimePeriods,
                          limit=limit, allData=allData, toTs=toTs, explainPath=explainPath, extraParams=extraParams, sign=sign)
    # RAISE ERROR...
    end

    try
        str = make_API_call(url)
        jobj = JSON.Parser.parse(str)
        if jobj["Response"] == "Success"
            datadict = jobj["Data"]["Data"]
            df = DataFrame(time=Int[], datetime=Dates.DateTime[], open=Float64[], close=Float64[], high=Float64[], low=Float64[], volumefrom=Float64[], volumeto=Float64[])

            nrows = size(datadict)[1]
            for i in 1:nrows
                push!(df, datadict[i], cols=:union)
            end
            df[:,:datetime] = unix2datetime.(df[:,:time])
            return df
        end
    catch e
        return e
    end
end

"""
    get_exchange_volume(duration::String, tsym::String; e::String="all", aggregate::Int=1, aggregatePredictableTimePeriods::Bool=true, limit::Int=30,toTs::String="", extraParams::String="", sign::Bool=false)

Get total volume from the daily historical exchange data. The values are based on 00:00 GMT time. We store the data in BTC and we multiply by the BTC-tsym value. If you want to get all the available historical data, you can use limit=2000 and keep going back in time using the toTs param. You can then keep requesting batches using: &limit=2000&toTs={the earliest timestamp received}.

1. tsym [Required]: The currency symbol to convert into [ Min length - 1] [ Max length - 30]
2. e: The exchange to obtain data from [ Min length - 2] [ Max length - 30] [ Default - CCCAGG]
3. aggregate: Time period to aggregate the data over (for daily it's days, for hourly it's hours and for minute histo it's minutes) [ Min - 1] [ Max - 30] [ Default - 1]
4. aggregatePredictableTimePeriods: Only used when the aggregate param is also in use. If false it will aggregate based on the current time.If the param is false and the time you make the call is 1pm - 2pm, with aggregate 2, it will create the time slots: ... 9am, 11am, 1pm.If the param is false and the time you make the call is 2pm - 3pm, with aggregate 2, it will create the time slots: ... 10am, 12am, 2pm.If the param is true (default) and the time you make the call is 1pm - 2pm, with aggregate 2, it will create the time slots: ... 8am, 10am, 12pm.If the param is true (default) and the time you make the call is 2pm - 3pm, with aggregate 2, it will create the time slots: ... 10am, 12am, 2pm. [ Default - true]
5. limit: The number of data points to return. If limit * aggregate > 2000 we reduce the limit param on our side. So a limit of 1000 and an aggerate of 4 would only return 2000 (max points) / 4 (aggregation size) = 500 total points + current one so 501. [ Min - 1] [ Max - 2000] [ Default - 30]
6. toTs: Returns historical data before that timestamp. If you want to get all the available historical data, you can use limit=2000 and keep going back in time using the toTs param. You can then keep requesting batches using: &limit=2000&toTs={the earliest timestamp received}
7. extraParams: The name of your application (we recommend you send it) [ Min length - 1] [ Max length - 2000] [ Default - NotAvailable]
8. sign: If set to true, the server will sign the requests (by default we don't sign them), this is useful for usage in smart contracts [ Default - false]
"""
function get_exchange_volume(duration::String, tsym::String; e::String="all", aggregate::Int=1, aggregatePredictableTimePeriods::Bool=true, limit::Int=30,
                             toTs::String="", extraParams::String="", sign::Bool=false)
    if duration in ["day", "hour"]
        url = build_exchange_vol_url(duration, tsym=tsym, e=e,
                            aggregate=aggregate, aggregatePredictableTimePeriods=aggregatePredictableTimePeriods,
                            limit=limit, toTs=toTs, extraParams=extraParams, sign=sign)
    # RAISE ERROR...
    end

    try
        str = make_API_call(url)
        jobj = JSON.Parser.parse(str)
        datadict = jobj["Data"]
        df = DataFrame(time=Int[], datetime=Dates.DateTime[], volume=Float64[])

        nrows = size(datadict)[1]
        for i in 1:nrows
            push!(df, datadict[i], cols=:union)
        end
        df[:,:datetime] = unix2datetime.(df[:,:time])
        return df
    catch e
        return e
    end
end

"""
    get_symbol_volume(duration::String, fsym::String, tsym::String; aggregate::Int=1, limit::Int=30, toTs::String="", extraParams::String="", sign::Bool=false)

Get total volumes from the daily historical symbol volume data. The volumes are stored in fsym (base) and are converted using the matching historical fsym - tsym conversion rates. For example if you are asking for the fsym BTC and tsym USD, the volume returned is the total BTC that was part of all trades (base and quote) converted using the BTC-USD exchange rate
The value is based on 00:00 GMT time and is calculated as the sum of all trades where the symbol/asset was involved either as the from/base or as the to/quote/counter. The top_tier_volume_quote and base are totals based on exchanges included in our exchange benchmark that are ranked as AA, A or B. cccagg_volume_quote and _base are based on the exchanges included in our CCCAGG average. total_volume_quote and _base are the sum of all trades across all markets where the asset/coin was either the fsym/base or the tsym/quote/counter
If you want to get all the available historical data, you can use limit=2000 and keep going back in time using the toTs param. You can then keep requesting batches using: &limit=2000&toTs={the earliest timestamp received}.

1. fsym [Required]: The cryptocurrency symbol of interest [ Min length - 1] [ Max length - 30]
2. tsym [Required]: The currency symbol to convert into [ Min length - 1] [ Max length - 30]
3. aggregate: Time period to aggregate the data over (for daily it's days, for hourly it's hours and for minute histo it's minutes) [ Min - 1] [ Max - 30] [ Default - 1]
4. limit: The number of data points to return. If limit * aggregate > 2000 we reduce the limit param on our side. So a limit of 1000 and an aggerate of 4 would only return 2000 (max points) / 4 (aggregation size) = 500 total points + current one so 501. [ Min - 1] [ Max - 2000] [ Default - 30]
5. toTs: Returns historical data before that timestamp. If you want to get all the available historical data, you can use limit=2000 and keep going back in time using the toTs param. You can then keep requesting batches using: &limit=2000&toTs={the earliest timestamp received}
6. extraParams: The name of your application (we recommend you send it) [ Min length - 1] [ Max length - 2000] [ Default - NotAvailable]
7. sign: If set to true, the server will sign the requests (by default we don't sign them), this is useful for usage in smart contracts [ Default - false]
"""
function get_symbol_volume(duration::String, fsym::String, tsym::String; aggregate::Int=1, limit::Int=30,
                             toTs::String="", extraParams::String="", sign::Bool=false)
    if duration in ["day", "hour"]
        url = build_symbol_vol_url(duration, fsym=fsym, tsym=tsym,
                            aggregate=aggregate,
                            limit=limit, toTs=toTs, extraParams=extraParams, sign=sign)
    # RAISE ERROR...
    end

    try
        str = make_API_call(url)
        jobj = JSON.Parser.parse(str)
        datadict = jobj["Data"]
        df = DataFrame()

        nrows = size(datadict)[1]
        for i in 1:nrows
            push!(df, datadict[i], cols=:union)
        end
        df[:,:datetime] = unix2datetime.(df[:,:time])
        return df
    catch e
        return e
    end
end

"""
    get_symbol_volume_single_exchange(duration::String, fsym::String, tsym::String; e::String="all", aggregate::Int=1, limit::Int=30, toTs::String="", extraParams::String="", sign::Bool=false)
Get the daily from and to volume for a symbol on an exchange.
The value is based on 00:00 GMT time and is calculated as the sum of all trades where the symbol/asset was involved either as the from/base or as the to/quote/counter on the requested exchange.
If you want to get all the available historical data, you can use limit=2000 and keep going back in time using the toTs param. You can then keep requesting batches using: &limit=2000&toTs={the earliest timestamp received}.

1. fsym [Required]: The cryptocurrency symbol of interest [ Min length - 1] [ Max length - 30]
2. tsym [Required]: The currency symbol to convert into [ Min length - 1] [ Max length - 30]
3. e: The exchange to obtain data from [ Min length - 2] [ Max length - 30] [ Default - CCCAGG]
4. aggregate: Time period to aggregate the data over (for daily it's days, for hourly it's hours and for minute histo it's minutes) [ Min - 1] [ Max - 30] [ Default - 1]
5. limit: The number of data points to return. If limit * aggregate > 2000 we reduce the limit param on our side. So a limit of 1000 and an aggerate of 4 would only return 2000 (max points) / 4 (aggregation size) = 500 total points + current one so 501. [ Min - 1] [ Max - 2000] [ Default - 168]
6. toTs: Returns historical data before that timestamp. If you want to get all the available historical data, you can use limit=2000 and keep going back in time using the toTs param. You can then keep requesting batches using: &limit=2000&toTs={the earliest timestamp received}
7. extraParams: The name of your application (we recommend you send it) [ Min length - 1] [ Max length - 2000] [ Default - NotAvailable]
8. sign: If set to true, the server will sign the requests (by default we don't sign them), this is useful for usage in smart contracts [ Default - false]
"""
function get_symbol_volume_single_exchange(duration::String, fsym::String, tsym::String; e::String="all", aggregate::Int=1, limit::Int=30, toTs::String="", extraParams::String="", sign::Bool=false)
    if duration in ["day", "hour"]
        url = build_symbol_vol_single_exchange_url(duration, fsym=fsym, tsym=tsym, e=e,
                            aggregate=aggregate,
                            limit=limit, toTs=toTs, extraParams=extraParams, sign=sign)
    # RAISE ERROR...
    end

    try
        str = make_API_call(url)
        jobj = JSON.Parser.parse(str)
        datadict = jobj["Data"]
        df = DataFrame()

        nrows = size(datadict)[1]
        for i in 1:nrows
            push!(df, datadict[i], cols=:union)
        end
        df[:,:datetime] = unix2datetime.(df[:,:time])
        return df
    catch e
        return e
    end
end


text:
using jInv.Mesh
using jInv.Utils
using Base.Test
using OED

r = 4;
l = 3;
Mimg  = getRegularMesh([-1 1 -1 1.],[4,4])
Mgrid = getRegularMesh([-1.2 2.4 -1.6 1.4],[8,8])

Ftrue = randn(Mimg.nc,10)
OEDparam = getOEDTypeBParam(r,l,Ftrue,Mimg,Mgrid,noiseLevel=0.)
display(OEDparam)

T = getTransmissionMatrix(OEDparam);

(Nt,nt,nrhst) = getProblemSize(OEDparam)
@test getNumberOfDesignParameters(OEDparam)==l

println("--- test getM and derivatives ---")
v   = randn(Mimg.nc)
p0  = rand(getNumberOfDesignParameters(OEDparam))
pass, = checkDerivative(p-> getM(p,OEDparam)*v, p-> getdpMx(p,v,OEDparam), p0,out=true)
@test pass

println("--- test getM^T and derivatives ---")
v   = randn(Nt)
p0  = rand(getNumberOfDesignParameters(OEDparam))
pass, = checkDerivative(p-> getM(p,OEDparam)'*v, p->getdpMTx(p,v,OEDparam), p0,out=false)
@test pass

println("--- test getD and derivatives ---")
pass, = checkDerivative(p->getD(p,OEDparam)[:,1], p->getdpD(p,1,OEDparam),p0,out=false)
@test pass


text:
using VLConstraintBasedModelGenerationUtilities

# test phrase -
reaction_phrase = "3*M_a_c+8*M_b_c+M_d_c"

# test:
species = "M_gtp_c"
stm_coeff = extract_stochiometric_coefficient_from_phrase(reaction_phrase, species)


text:
# Metagraphs files are simply JLD2 files.

"""
    struct MGFormat <: AbstractGraphFormat end

You can save `MetaGraph`s in a `MGFormat`, currently based on `JLD2`.
"""
struct MGFormat <: AbstractGraphFormat end

"""
    struct DOTFormat <: AbstractGraphFormat end

If all metadata types support `pairs` or are `nothing`, you can save `MetaGraph`s in `DOTFormat`.
"""
struct DOTFormat <: AbstractGraphFormat end

function loadmg(fn::AbstractString)
    @load fn g
    return g
end

function savemg(fn::AbstractString, g::MetaGraph)
    @save fn g
    return 1
end

loadgraph(fn::AbstractString, gname::String, ::MGFormat) = loadmg(fn)
savegraph(fn::AbstractString, g::MetaGraph) = savemg(fn, g)

function show_meta_list(io::IO, meta)
    if meta !== nothing && length(meta) > 0
        next = false
        write(io, " [")
        for (key, value) in meta
            if next
                write(io, ", ")
            else
                next = true
            end
            write(io, key)
            write(io, " = ")
            show(io, value)
        end
        write(io, "]")
    end
end

function savedot(io::IO, g::MetaGraph)
    gprops = g.gprops
    metaindex = g.metaindex

    if is_directed(g)
        write(io, "digraph G {\n")
        dash = "->"
    else
        write(io, "graph T {\n")
        dash = "--"
    end

    if gprops !== nothing
        for (key, value) in pairs(gprops)
            write(io, "    ")
            write(io, key)
            write(io, " = ")
            show(io, value)
            write(io, '\n')
        end
    end

    for label in keys(g.vprops)
        write(io, "    ")
        write(io, label)
        show_meta_list(io, g[label])
        write(io, '\n')
    end

    for (label_1, label_2) in keys(g.eprops)
        write(io, "    ")
        write(io, label_1)
        write(io, ' ')
        write(io, dash)
        write(io, ' ')
        write(io, label_2)
        show_meta_list(io, g.eprops[arrange(g, label_1, label_2)])
        write(io, "\n")
    end
    write(io, "}\n")
end

function savegraph(fn::AbstractString, g::MetaGraph, ::DOTFormat)
    open(fn, "w") do fp
        savedot(fp, g)
    end
end


text:
export MPS,
       sample,
       sample!,
       leftLim,
       prime!,
       primelinks!,
       simlinks!,
       inner,
       productMPS,
       randomMPS,
       replaceBond!,
       rightLim,
       maxLinkDim,
       linkindex,
       siteindex,
       siteinds

mutable struct MPS
  N_::Int
  A_::Vector{ITensor}
  llim_::Int
  rlim_::Int

  MPS() = new(0,Vector{ITensor}(),0,0)

  MPS(N::Int) = MPS(N,Vector{ITensor}(undef,N),0,N+1)

  function MPS(N::Int,
               A::Vector{<:ITensor},
               llim::Int=0,
               rlim::Int=N+1)
    new(N,A,llim,rlim)
  end
end

function MPS(sites)
  N = length(sites)
  v = Vector{ITensor}(undef, N)
  l = [Index(1, "Link,l=$ii") for ii=1:N-1]
  for ii in eachindex(sites)
    s = sites[ii]
    if ii == 1
      v[ii] = ITensor(l[ii], s)
    elseif ii == N
      v[ii] = ITensor(l[ii-1], s)
    else
      v[ii] = ITensor(l[ii-1],s,l[ii])
    end
  end
  return MPS(N,v,0,N+1)
end

length(m::MPS) = m.N_
tensors(m::MPS) = m.A_
leftLim(m::MPS) = m.llim_
rightLim(m::MPS) = m.rlim_

function setLeftLim!(m::MPS,new_ll::Int)
  m.llim_ = new_ll
end

function setRightLim!(m::MPS,new_rl::Int)
  m.rlim_ = new_rl
end

isOrtho(m::MPS) = (leftLim(m)+1 == rightLim(m)-1)

function orthoCenter(m::MPS)
  !isOrtho(m) && error("MPS has no well-defined orthogonality center")
  return leftLim(m)+1
end

getindex(M::MPS, n::Integer) = getindex(tensors(M),n)

function setindex!(M::MPS,T::ITensor,n::Integer)
  (n <= leftLim(M)) && setLeftLim!(M,n-1)
  (n >= rightLim(M)) && setRightLim!(M,n+1)
  setindex!(tensors(M),T,n)
end

copy(m::MPS) = MPS(m.N_,copy(tensors(m)),m.llim_,m.rlim_)
similar(m::MPS) = MPS(m.N_, similar(tensors(m)), 0, m.N_)

eachindex(m::MPS) = 1:length(m)

function show(io::IO, M::MPS)
  print(io,"MPS")
  (length(M) > 0) && print(io,"\n")
  for (i, A) ∈ enumerate(tensors(M))
    println(io,"$i  $(inds(A))")
  end
end

function randomMPS(sites)
  M = MPS(sites)
  for i in eachindex(sites)
    randn!(M[i])
    normalize!(M[i])
  end
  M.llim_ = 1
  M.rlim_ = length(M)
  return M
end

function productMPS(ivals::Vector{IndexVal})
  N = length(ivals)
  As = Vector{ITensor}(undef,N)
  links  = Vector{Index}(undef,N)
  for n=1:N
    s = ind(ivals[n])
    links[n] = Index(1,"Link,l=$n")
    if n == 1
      A = ITensor(s,links[n])
      A[ivals[n],links[n](1)] = 1.0
    elseif n == N
      A = ITensor(links[n-1],s)
      A[links[n-1](1),ivals[n]] = 1.0
    else
      A = ITensor(links[n-1],s,links[n])
      A[links[n-1](1),ivals[n],links[n](1)] = 1.0
    end
    As[n] = A
  end
  return MPS(N,As,0,2)
end

function productMPS(sites,
                    states)
  if length(sites) != length(states)
    throw(DimensionMismatch("Number of sites and and initial states don't match"))
  end
  ivals = [state(sites[n],states[n]) for n=1:length(sites)]
  return productMPS(ivals)
end

function linkindex(M::MPS,j::Integer)
  N = length(M)
  j ≥ length(M) && error("No link index to the right of site $j (length of MPS is $N)")
  li = commonindex(M[j],M[j+1])
  if isnothing(li)
    error("linkindex: no MPS link index at link $j")
  end
  return li
end

function siteindex(M::MPS,j::Integer)
  N = length(M)
  if j == 1
    si = uniqueindex(M[j],M[j+1])
  elseif j == N
    si = uniqueindex(M[j],M[j-1])
  else
    si = uniqueindex(M[j],(M[j-1],M[j+1]))
  end
  return si
end

function siteinds(M::MPS)
  return [siteindex(M,j) for j in 1:length(M)]
end

function replacesites!(M::MPS,sites)
  for j in eachindex(M)
    sj = siteindex(M,j)
    replaceindex!(M[j],sj,sites[j])
  end
  return M
end

function inner(M1::MPS, M2::MPS)::Number
  N = length(M1)
  if length(M2) != N
      throw(DimensionMismatch("inner: mismatched lengths $N and $(length(M2))"))
  end
  M1dag = dag(M1)
  simlinks!(M1dag)
  O = M1dag[1]*M2[1]
  for j in eachindex(M1)[2:end]
    O = (O*M1dag[j])*M2[j]
  end
  return O[]
end

function replaceBond!(M::MPS,
                      b::Int,
                      phi::ITensor;
                      kwargs...)
  FU,FV = factorize(phi,inds(M[b]); which_factorization="automatic",
                        tags=tags(linkindex(M,b)), kwargs...)
  M[b]   = FU
  M[b+1] = FV

  dir = get(kwargs,:dir,"center")
  if dir=="fromright"
    M.llim_ ==b && (M.llim_ -= 1)
    M.rlim_ == b+2 && (M.rlim_ -= 1)
  elseif dir=="fromleft"
    M.llim_ == b-1 && (M.llim_ += 1)
    M.rlim_ == b+1 && (M.rlim_ += 1)
  end
end

"""
    sample!(m::MPS)

Given a normalized MPS m, returns a `Vector{Int}`
of `length(m)` corresponding to one sample
of the probability distribution defined by
squaring the components of the tensor
that the MPS represents. If the MPS does
not have an orthogonality center,
orthogonalize!(m,1) will be called before
computing the sample.
"""
function sample!(m::MPS)
  orthogonalize!(m,1)
  return sample(m)
end

"""
    sample(m::MPS)

Given a normalized MPS m with `orthoCenter(m)==1`,
returns a `Vector{Int}` of `length(m)`
corresponding to one sample of the
probability distribution defined by
squaring the components of the tensor
that the MPS represents
"""
function sample(m::MPS)
  N = length(m)

  if orthoCenter(m) != 1
    error("sample: MPS m must have orthoCenter(m)==1")
  end
  if abs(1.0-norm(m[1])) > 1E-8
    error("sample: MPS is not normalized, norm=$(norm(m[1]))")
  end

  result = zeros(Int,N)
  A = m[1]

  for j=1:N
    s = siteindex(m,j)
    d = dim(s)
    # Compute the probability of each state
    # one-by-one and stop when the random
    # number r is below the total prob so far
    pdisc = 0.0
    r = rand()
    # Will need n,An, and pn below
    n = 1
    An = ITensor()
    pn = 0.0
    while n <= d
      projn = ITensor(s)
      projn[s[n]] = 1.0
      An = A*projn
      pn = scalar(dag(An)*An)
      pdisc += pn
      (r < pdisc) && break
      n += 1
    end

    result[j] = n

    if j < N
      A = m[j+1]*An
      A *= (1.0/sqrt(pn))
    end
  end
  return result
end

@doc """
inner(psi::MPS, phi::MPS)

Compute <psi|phi>
""" inner


text:
"""
function to find MPS approximation of exp(A)v using the Arnoldi method with
numv Krylov space vectors
"""

function Arnoldi_expm(A, v, numv)

    rhov = Array{Array{Array{Complex{Float64}, 3}, 1}, 1}(numv)
    vprod = Array{Array{Array{Complex{Float64}, 3}, 1}, 1}(numv)
    N = complex(zeros(numv, numv))
    H = complex(zeros(numv, numv))

    rhov[1] = v
    vprod[1] = applyMPOtoMPS(rhov[1], A)
    H[1] = scal_prod(rhov[1], vprod[1])
    N[1] = scal_prod(rhov[1], rhov[1])

    coef = complex(zeros(numv))

    for jj = 2:numv

        coef[1:jj - 1] = [-H[ii, jj - 1]/N[ii, ii] for ii in 1:jj - 1]
        rhov[jj] = compress_sum_var(v, [1.0, coef[1:jj - 1]...],
                                                  [vprod[jj - 1], rhov[1:jj - 1]...], 4)
        vprod[jj] = applyMPOtoMPS(rhov[jj], A)
        H[jj, jj] = scal_prod(rhov[jj], vprod[jj])
        N[jj, jj] = scal_prod(rhov[jj], rhov[jj])

        for ll = 1:jj - 1

            H[jj, ll] = scal_prod(rhov[jj], vprod[ll])
            H[ll, jj] = scal_prod(rhov[ll], vprod[jj])
            N[jj, ll] = scal_prod(rhov[jj], rhov[ll])
            N[ll, jj] = conj(N[jj, ll])

        end

    end
    println(H)
    println(diag(N))
    coef = expm(N\H)[:, 1]
    println(coef)
    return compress_sum_var(v, coef, rhov, 4)

end

function Arnoldi_expm_norm(A, v, numv)

    rhov = Array{Array{Array{Complex{Float64}, 3}, 1}, 1}(numv)
    vprod = Array{Array{Array{Complex{Float64}, 3}, 1}, 1}(numv)
    N = complex(eye(numv))
    H = complex(zeros(numv, numv))

    rhov[1] = deepcopy(v)
    v_norm = normalize_lo!(rhov[1])
    vprod[1] = applyMPOtoMPS(rhov[1], A)
    H[1] = scal_prod(rhov[1], vprod[1])

    coef = complex(zeros(numv))

    for jj = 2:numv

        coef[1:jj - 1] = [-H[ii, jj - 1]/N[ii, ii] for ii in 1:jj - 1]
        rhov[jj] = compress_sum_var(v, [1.0, coef[1:jj - 1]...],
                                                  [vprod[jj - 1], rhov[1:jj - 1]...], 1)
        no = normalize_lo!(rhov[jj])
        println(no)
        vprod[jj] = applyMPOtoMPS(rhov[jj], A)
        H[jj, jj] = scal_prod(rhov[jj], vprod[jj])

        for ll = 1:jj - 1

            H[jj, ll] = scal_prod(rhov[jj], vprod[ll])
            H[ll, jj] = scal_prod(rhov[ll], vprod[jj])
            N[jj, ll] = scal_prod(rhov[jj], rhov[ll])
            N[ll, jj] = conj(N[jj, ll])

        end

    end
    coef = sqrt(v_norm)*expm(N\H)[:, 1]
    println(coef)
    return compress_sum_var(v, coef, rhov, 1)

end

function Arnoldi_expm!(A, rhov, numv, env, envop)

    N = complex(eye(numv))
    H = complex(zeros(numv, numv))

    H[1, 1] = scal_op_prod(rhov[1], A, rhov[1])

    for jj = 2:numv

        compress_sum_var_apply_H!(rhov[jj], envop,
                                    [- H[1:(jj - 1), jj - 1]..., 1.0],
                                        [rhov[1:(jj - 1)]..., rhov[jj - 1]], A, 1)
        normalize_lo!(rhov[jj])
        H[jj, jj] = scal_op_prod(rhov[jj], A, rhov[jj])

        for ll = 1:jj - 1

            H[jj, ll] = scal_op_prod(rhov[jj], A, rhov[ll])
            H[ll, jj] = scal_op_prod(rhov[ll], A, rhov[jj])
            N[jj, ll] = scal_prod(rhov[jj], rhov[ll])
            N[ll, jj] = conj(N[jj, ll])

        end

    end
    coef = expm(N\H)[:, 1]
    println(abs.(coef))
    rhov[1] = compress_sum_var(rhov[1], coef, rhov, 1)
    normalize_lo!(rhov[1])
    nothing

end
"""
function to evolve MPS using the 4th order Runga Kutta method
"""

function RK4stp(dt,H,psi,k0)

k1 = compress_var(k0,applyMPOtoMPS(psi,H[1]),1)
k2 = compress_var(k1,applyMPOtoMPS(
    compress_sum_var(psi,[1,dt/2],[psi,k1],1),H[2]),1)
k3 = compress_var(k2,applyMPOtoMPS(
    compress_sum_var(psi,[1,dt/2],[psi,k2],1),H[2]),1)
k4 = compress_var(k3,applyMPOtoMPS(
    compress_sum_var(psi,[1,dt],[psi,k3],1),H[3]),1)

psi = compress_sum_var(psi,[1,dt/6,dt/3,dt/3,dt/6],[psi,k1,k2,k3,k4],1)
return psi,k4

end

"""
function to evolve MPS using the 4th order Runga Kutta method
"""

function RK4stp_alt(dt,H,psi)

psi1 = compress_sum_var(psi,[1,dt/2],[psi,applyMPOtoMPS(psi,H[1])],1)
psi2 = compress_sum_var(psi1,[1,dt/2],[psi,applyMPOtoMPS(psi1,H[2])],1)
psi3 = compress_sum_var(psi2,[1,dt],[psi,applyMPOtoMPS(psi2,H[2])],1)
psi = compress_sum_var(psi3,[-1/3,1/3,2/3,1/3,dt/6],[psi,psi1,psi2,psi3,
                       applyMPOtoMPS(psi3,H[3])],1)

end

"""
function to evolve MPS using the 4th order Runga Kutta method
"""

function RK4stp_apply_H(dt,H,psi)

psi1 = compress_var_apply_H(psi,psi,H[1],1)
psi2 = compress_sum_var_apply_H(psi1,[1,dt/2],[psi,psi1],H[2],1)
psi3 = compress_sum_var_apply_H(psi2,[1,dt],[psi,psi2],H[2],1)
psi = compress_sum_var_apply_H(psi3,[-1/3,1/3,2/3,1/3],[psi,psi1,psi2,
                               psi3],H[3],1)

end

function RK4stp_apply_H_half(dt, H, psi)

psi1 = compress_var_apply_H_left(psi, psi, H[1])
psi2 = compress_sum_var_apply_H_right(psi1, [1, dt/2], [psi, psi1], H[2])
psi3 = compress_sum_var_apply_H_left(psi2, [1, dt], [psi, psi2], H[2])
compress_sum_var_apply_H_right(psi3, [-1/3, 1/3, 2/3, 1/3],
                                           [psi, psi1, psi2, psi3], H[3])

end

function RK4stp_apply_H_half!(dt, H, psi, psi1, psi2, psi3, env1, env2, env3, envop)

compress_var_apply_H_left!(psi1, envop, psi, psi, H[1])
compress_sum_var_apply_H_right!(psi2, [env1], envop, psi1, [1, dt/2], [psi, psi1], H[2])
compress_sum_var_apply_H_left!(psi3, [env1], envop, psi2, [1, dt], [psi, psi2], H[2])
compress_sum_var_apply_H_right!(psi, [env1, env2, env3], envop, psi3, [-1/3, 1/3, 2/3, 1/3],
                                [psi, psi1, psi2, psi3], H[3])

end

function RK4stp_apply_H_half!(dt, H, psi0, psi, psi1, psi2, psi3, env1, env2, env3, envop)

compress_var_apply_H_left!(psi1, envop, psi, psi, H[1])
compress_sum_var_apply_H_right!(psi2, [env1], envop, psi1, [1, dt/2], [psi, psi1], H[2])
compress_sum_var_apply_H_left!(psi3, [env1], envop, psi2, [1, dt], [psi, psi2], H[2])
compress_sum_var_apply_H_right!(psi0, [env1, env2, env3], envop, psi3, [-1/3, 1/3, 2/3, 1/3],
                                [psi, psi1, psi2, psi3], H[3])

end

function RK4stp_apply_H_half_two_site(dt,H,psi,svd_tol,dims,dlim)

psi1 = compress_var_apply_H_two_site_left(psi,psi,H[1],1,svd_tol,dims,dlim)
psi2 = compress_sum_var_apply_H_right(psi1,[1,dt/2],[psi,psi1],H[2],1)
psi3 = compress_sum_var_apply_H_left(psi2,[1,dt],[psi,psi2],H[2],1)
psi = compress_sum_var_apply_H_right(psi3,[-1/3,1/3,2/3,1/3],[psi,psi1,psi2,
                               psi3],H[3],1)

end

function RK4stp_apply_H_half_full_two_site(dt,H,psi,svd_tol,dims,dlim)

psi1 = compress_var_apply_H_two_site_left(psi,psi,H[1],1,svd_tol,dims,dlim)
psi2 = compress_sum_var_apply_H_two_site_right(psi1,[1,dt/2],[psi,psi1],H[2],1,svd_tol,dims,dlim)
psi3 = compress_sum_var_apply_H_two_site_left(psi2,[1,dt],[psi,psi2],H[2],1,svd_tol,dims,dlim)
psi = compress_sum_var_apply_H_two_site_right(psi3,[-1/3,1/3,2/3,1/3],[psi,psi1,psi2,
                               psi3],H[3],1,svd_tol,dims,dlim)

end

function RK4stp_apply_H_half_two_site_2(dt,H,psi,svd_tol,dims,dlim)

psi1 = compress_var_apply_H_left(psi,psi,H[1],1)
psi2 = compress_sum_var_apply_H_right(psi1,[1,dt/2],[psi,psi1],H[2],1)
psi3 = compress_sum_var_apply_H_left(psi2,[1,dt],[psi,psi2],H[2],1)
psi = compress_sum_var_apply_H_two_site_right(psi3,[-1/3,1/3,2/3,1/3],[psi,psi1,psi2,
                               psi3],H[3],1,svd_tol,dims,dlim)

end

"""
function to evolve MPS A (must be in left canonical form)
in time acording to the trapezoidal method
"""

function Evolve_trap(A,LF,LB,maxsweeps,nor)

	na = length(A)
	At = deepcopy(A)
	FL = Array{Array}(na)
	FR = Array{Array}(na)
	SL = Array{Array}(na)
	SR = Array{Array}(na)
	XY = applyMPOtoMPO(LF,Conj_MPO(LB))
	XX = applyMPOtoMPO(LB,Conj_MPO(LB))

	FL[1] = ones(1,1,1)
	SL[1] = ones(1,1,1)
	for jj = 1:na-1
		FL[jj+1] = update_le(At[jj],XY[jj],A[jj],FL[jj])
		SL[jj+1] = update_le(At[jj],XX[jj],At[jj],SL[jj])
	end
	FR[na] = ones(1,1,1)
	SR[na] = ones(1,1,1)

	for sweep = 1:maxsweeps

		for jj = na:-1:2
		    HH = reshape(O_prod_LR(SL[jj],XX[jj],SR[jj]),
		        size(SL[jj],1)*size(XX[jj],2)*size(SR[jj],1),
		        size(SL[jj],1)*size(XX[jj],2)*size(SR[jj],1))
		    HE = reshape(A_prod_LOR(FL[jj],A[jj],XY[jj],FR[jj]),
		        size(SL[jj],1)*size(XX[jj],2)*size(SR[jj],1),1)
		    At[jj] = reshape(HH\HE,size(SL[jj],1),size(XX[jj],2),size(SR[jj],1))
		    Q,R = qr(reshape(At[jj],size(At[jj],1),size(At[jj],3)*size(At[jj],2))')
		    At[jj] = reshape(Q',size(Q,2),size(At[jj],2),size(At[jj],3))
		    FR[jj-1] = update_re(At[jj],XY[jj],A[jj],FR[jj])
		    SR[jj-1] = update_re(At[jj],XX[jj],At[jj],SR[jj])
		end
	    HH = reshape(O_prod_LR(SL[1],XX[1],SR[1]),
	        size(SL[1],1)*size(XX[1],2)*size(SR[1],1),
	        size(SL[1],1)*size(XX[1],2)*size(SR[1],1))
	    HE = reshape(A_prod_LOR(FL[1],A[1],XY[1],FR[1]),
	        size(SL[1],1)*size(XX[1],2)*size(SR[1],1),1)
	    At[1] = reshape(HH\HE,size(SL[1],1),size(XX[1],2),size(SR[1],1));

		for jj = 1:na-1;
		    HH = reshape(O_prod_LR(SL[jj],XX[jj],SR[jj]),
		        size(SL[jj],1)*size(XX[jj],2)*size(SR[jj],1),
		        size(SL[jj],1)*size(XX[jj],2)*size(SR[jj],1))
		    HE = reshape(A_prod_LOR(FL[jj],A[jj],XY[jj],FR[jj]),
		        size(SL[jj],1)*size(XX[jj],2)*size(SR[jj],1),1)
		    At[jj] = reshape(HH\HE,size(SL[jj],1),size(XX[jj],2),size(SR[jj],1));
		    Q,R = qr(reshape(At[jj],size(At[jj],1)*size(At[jj],2),size(At[jj],3)));
		    At[jj] = reshape(Q,size(At[jj],1),size(At[jj],2),size(Q,2))
		    FL[jj+1] = update_le(At[jj],XY[jj],A[jj],FL[jj]);
		    SL[jj+1] = update_le(At[jj],XX[jj],At[jj],SL[jj]);
		end
	    HH = reshape(O_prod_LR(SL[na],XX[na],SR[na]),
	        size(SL[na],1)*size(XX[na],2)*size(SR[na],1),
	        size(SL[na],1)*size(XX[na],2)*size(SR[na],1))
	    HE = reshape(A_prod_LOR(FL[na],A[na],XY[na],FR[na]),
	        size(SL[na],1)*size(XX[na],2)*size(SR[na],1),1)
	    At[na] = reshape(HH\HE,size(SL[na],1),size(XX[na],2),size(SR[na],1))

	end

	no = sum(abs2(At[na][:]))

	if nor
	    At[na] = At[na]./sqrt(no)
	end

	return At, no

end


text:
# Optical flow `.flo` fileformat¹² as used by various optical flow
# benchmarks³⁴.
#
# Implementation notes:
#   - to match behaviour for image data the flow vector field is transformed to
#     julia col-major (i.e. height is first dimension)
#   - missing data (i.e. any component `> 1f9`) is represented using the
#     Missing type in julia
#   - as satisfied per spec `floatmax(Float32)` is used to save missing data
#
# [1] https://vision.middlebury.edu/flow/code/flow-code/README.txt
# [2] https://doi.org/10.1007/s11263-010-0390-2
# [3] https://vision.middlebury.edu/flow/
# [4] http://sintel.is.tue.mpg.de/

using FileIO: File, Stream, @format_str, skipmagic, magic

const iodims = 2
const T = Float32
const T_ = Union{Missing,T}
const nodata = floatmax(T)

function load(f::File{format"FLO"})
    open(f) do s
        skipmagic(s)
        load(s)
    end
end

function load(s::Stream{format"FLO"})
    sz = ntuple(iodims) do k
        read(s, Int32) |> ltoh
    end

    data = Array{T_}(undef, iodims, sz...)
    for I in eachindex(data)
        data[I] = read(s, T) |> ltoh
    end

    for J in CartesianIndices(sz)
        Jv = view(data, :, J)
        # spec wording is misleading: infering from the ground truth files
        # it should be "or" instead of "either or"
        if any(x -> abs(x) > 1f9, Jv)
            Jv .= missing
        end
    end
    # row- to col-major
    reverse!(data, dims = 1)
    return permutedims(data, (1, 3, 2, (4:iodims + 1)...))
end

function save(f::File{format"FLO"}, data::AbstractArray)
    open(f, "w") do s
        write(s, magic(format"FLO"))
        save(s, data)
    end
end

function save(s::Stream{format"FLO"}, data::AbstractArray)
    size(data, 1) != ndims(data) - 1 &&
        throw(ArgumentError("invalid dimensions $(size(data)) for flow field"))

    # col- to row-major
    data = permutedims(data, (1, 3, 2, (4:iodims + 1)...))
    reverse!(data, dims = 1)

    for k in 1:iodims
        write(s, Int32(size(data, k + 1)) |> htol)
    end

    for I in eachindex(data)
        v = convert(T, coalesce(data[I], nodata))
        write(s, v |> htol)
    end
end


text:
# domain spec.jl

abstract type ARDomain end

abstract type Invoice <: ARDomain end

abstract type Structure <: ARDomain  end

abstract type BodyItem <: Structure end

abstract type Payment <: ARDomain  end

"""
    meta(i::Invoice)

Returns the meta data of an invoice.

# example
```
julia> using AppliSales

julia> using AppliAR

julia> orders = AppliSales.process();

julia> AppliAR.process(orders);

julia> unpaid_invoices = AppliAR.retrieve_unpaid_invoices();

julia> m = AppliAR.meta(unpaid_invoices[1])
MetaInvoice("9715406426271665630", "LS", 2020-07-14T16:56:33.194, "€", 1.0)

julia> orderid = order_id(m)
"9715406426271665630"

julia> trainingid = training_id(m)
"LS"

julia> date(m)
2020-07-14T16:56:33

julia> currency(m)
"€"

julia> currency_ratio(m)
1.0
```
"""
function meta end

"""
    header(i::Invoice)

Returns the header of an invoice.

example
```
julia> using AppliSales

julia> using AppliAR

julia> orders = AppliSales.process();

julia> AppliAR.process(orders);

julia> unpaid_invoices = AppliAR.retrieve_unpaid_invoices();

julia> h = header(unpaid_invoices[1])
Header("A1001", "Scrooge Investment Bank", "1180 Seven Seas Dr", "FL 32830", "Lake Buena Vista", "USA", "PO-456", "Scrooge McDuck", "scrooge@duckcity.com")

julia> invoice_nbr(h)
"A1001"

julia> name(h)
"Scrooge Investment Bank"

julia> address(h)
"1180 Seven Seas Dr"

julia> postal_code(h)
"FL 32830"

julia> city(h)
"Lake Buena Vista"

julia> country(h)
"USA"

julia> order_ref(h)
"PO-456"

julia> name_contact(h)
"Scrooge McDuck

julia> email_contact(h)
"scrooge@duckcity.com"
```
"""
function header end

"""
    body(i::Invoice)

Returns the body of an invoice.

example
```
julia> using AppliSales

julia> using AppliAR

julia> orders = AppliSales.process();

julia> AppliAR.process(orders);

julia> unpaid_invoices = AppliAR.retrieve_unpaid_invoices();

julia> b = body(unpaid_invoices[1])
OpentrainingItem("Learn Smiling", 2019-08-30T00:00:00, 1000.0, ["Scrooge McDuck"], 0.21)

julia> name_training(b)
"Learn Smiling"

julia> date(b)
2019-08-30T00:00:00

julia> price_per_student(b)
1000.0

julia> students(b)
1-element Array{String,1}:
"Scrooge McDuck"

julia> vat_perc(b)
0.21
```
"""
function body end

"""
    id(i::Invoice)

Returns the id of an invoice.

example
```
julia> using AppliSales

julia> using AppliAR

julia> orders = AppliSales.process();

julia> AppliAR.process(orders);

julia> unpaid_invoices = AppliAR.retrieve_unpaid_invoices();

julia> id(unpaid_invoices[1])
"A1001"
```
"""
function id end

"""
    stm(i::PaidInvoice)

Returns the Bankstatement of a PaidInvoice.
"""
function stm end

max_stars_repo_name: rbontekoe/AppliAR

max_issues_repo_name: rbontekoe/AppliAR

max_forks_repo_name: rbontekoe/AppliAR


text:
include("activationsGPU.jl")
include("NeuralNetwork.jl")
using CuArrays, Statistics

"""
    initialize_parameters(layer_dims::CuArray{Int})::Dict{String, CuArray{Float64}}

Create and initialize a parameters dict to store parameters of a L-layer feedforward neural network.
The first argument to this function is an array representing number of units in each layer.
The second argument is a placeholder to infer the type of Arrays to initialize - either a normal array or a GPU based CuArray.

Parameters Wi (weights) and bi (biases) are initialized for every layer other than the input layer.
Weights are intialized randomly using Xavier initialization method (to normalize them) and
biases are initialized to zeros.

# Examples
```jldoctest
julia> parameters = initialize_parameters([5 10 1], CuArray([0]))
Dict{String,CuArray{Float32,N,P} where P where N} with 4 entries:
  "W2" => Float32[-0.00542962 0.00721382 … 0.185659 -0.355868]
  "W1" => Float32[0.679642 -0.365945 … 0.145359 -0.47302; 0.200724 -0.339901 … 1.11603 0.551288; … ; -0.93578 -0.0143422 … 0.433826 -0.425189; -0.70129 -0.077…
  "b2" => Float32[0.0]
  "b1" => Float32[0.0; 0.0; … ; 0.0; 0.0]

julia> for (key, value) in parameters
           println(key, " ", size(value))
       end
W2 (1, 10)
W1 (10, 5)
b2 (1, 1)
b1 (10, 1)
```
"""
function initialize_parameters(layer_dims::Array{Int}, Y::CuArray)::Dict{String, CuArray{Float32}}
    parameters = Dict{String, CuArray{Float32}}()
    for i = 2:length(layer_dims)
        parameters[string("W", i-1)] = CuArrays.randn(layer_dims[i], layer_dims[i-1]) ./ sqrt(layer_dims[i-1]) # Xavier initialization
        parameters[string("b", i-1)] = CuArrays.zeros(layer_dims[i], 1)
    end
    return parameters
end

function forward_prop(X::CuArray{Float32}, parameters::Dict{String, CuArray{Float32}}, activations::Tuple)::Tuple{CuArray{Float32}, Dict{String, CuArray{Float32}}}
    num_layers = length(parameters) ÷ 2

    caches = Dict{String, CuArray{Float32}}()

    caches[string("A", 0)] = X

    Ai = Nothing
    for i = 1:num_layers
        Zi = parameters[string("W", i)] * caches[string("A", i-1)] .+ parameters[string("b", i)] # * is equivivalent to np.dot
        Ai = activations[i].(Zi)
        caches[string("Z", i)] = Zi
        caches[string("A", i)] = Ai
    end

    return Ai, caches
end

function cost_binary(Y::CuArray{Float32}, Ŷ::CuArray{Float32})::Float32
    @assert length(Y) == length(Ŷ)
    m = length(Y)

    cost = - sum(Y .* log.(Ŷ) .+ (1 .- Y) .* log.(1 .- Ŷ)) / m
    # println("cost is ", string(cost, Y[1:30], Ŷ[1:30]))
    return cost
end

function backward_prop(Y::CuArray{Float32}, Ŷ::CuArray{Float32}, parameters::Dict{String, CuArray{Float32}}, caches::Dict{String, CuArray{Float32}}, layer_dims::Array{Int}, activations::Tuple)::Dict{String, CuArray{Float32}}

    num_layers = length(layer_dims)
    @assert length(Y) == length(Ŷ)
    m = size(Y)[2]
    # println("Number of examples(m) ", m)

    # activations_orig = activations
    # activations = []
    # for activation in activations_orig
    #     push!(activations, @eval ($(Symbol("$activation", "_back"))))
    # end
    # println(activations)

    dA = (.- Y ./ Ŷ .+ (1 .- Y) ./ (1 .- Ŷ))
    if all(isnan.(dA))
        println("dA was NaN!")
        dA = randn(Float32)
    end
    # println("dA: ", dA)

    grads = Dict{String, CuArray{Float32}}()

    for l in num_layers-1:-1:1
        dZ = dA .* activations[l].(caches[string("Z", l)])
        grads[string("dw", l)] = 1/m .* (dZ * transpose(caches[string("A", l-1)]))
        grads[string("db", l)] = 1/m .* sum(dZ, dims=2)
        dA = transpose(parameters[string("W", l)]) * dZ
    end

    return grads
end

function update_parameters_old(parameters::Dict{String, CuArray{Float32}}, grads::Dict{String, CuArray{Float32}}, layer_dims::Array{Int}, learning_rate::Number)::Dict{String, CuArray{Float32}}
    num_layers = length(layer_dims)
    for l = 1:num_layers-1
        parameters[string("W", l)] -= learning_rate .* grads[string("dw", l)]
        parameters[string("b", l)] -= learning_rate .* grads[string("db", l)]
        # if l > 1
        #     println("dw ", l, grads[string("dw", l)])
        #     println("db ", l, grads[string("db", l)])
        # end
    end
    return parameters
end

function update_params_kernel!(params, lr, grads, m, n)
    tx = (blockIdx().x-1) * blockDim().x + threadIdx().x
    ty = (blockIdx().y-1) * blockDim().y + threadIdx().y
    if tx <= m && ty <= n
        @inbounds params[tx, ty] = params[tx, ty] - lr * grads[tx, ty]
    end
    return nothing
end

function get_threads_blocks(m, n)
    nthreads = 16
    threads = (min(nthreads, m), min(nthreads, n))
    blocks = (cld(m, threads[1]), cld(n, threads[2]))
    return threads, blocks
end

function update_parameters(parameters::Dict{String, CuArray{Float32}}, grads::Dict{String, CuArray{Float32}}, layer_dims::Array{Int}, learning_rate::Number)::Dict{String, CuArray{Float32}}
    num_layers = length(layer_dims)
    for l = 1:num_layers-1
        m,n = size(parameters[string("W", l)])
        threads, blocks = get_threads_blocks(m, n)
        @cuda blocks=blocks threads=threads update_params_kernel!(parameters[string("W", l)], learning_rate, grads[string("dw", l)], m, n)

        m,n = size(parameters[string("b", l)])
        threads, blocks = get_threads_blocks(m, n)
        @cuda blocks=blocks threads=threads update_params_kernel!(parameters[string("b", l)], learning_rate, grads[string("db", l)], m, n)
    end
    synchronize()
    return parameters
end

function reshape_Y(Y::CuArray)
    Y = convert(CuArray{Float32, ndims(Y)}, Y)
    Y = reshape(Y, 1, :)
    return Y
end


text:
# jim.jl

using Plots
using LaTeXStrings
using MosaicViews
using FFTViews

# global default key/values
jim_def = Dict([
	:aspect_ratio => :equal,
	:clim => [],
	:color => :grays,
	:ncol => 0,
	:padval => 0,
	:mosaic_npad => 1,
	:title => "",
	:xlabel => "",
	:ylabel => "",
	:fft0 => false,
	:yflip => nothing, # defer to minimum value of y - see default below
	:abswarn => true, # warn when taking abs of complex images?
	])

"""
`jim(z, ...)`

jiffy image display of `x` using `heatmap`

in
* `z` image, can be 2D or higher, if higher then it uses `mosaicviews`

option
* `aspect_ratio` for heatmap; default `:equal`
* `clim` for heatmap; default `(minimum(z),maximum(z))`
* `color` colormap; default `:grays`
* `ncol` for mosaicview for 3D and higher arrays; default `0` does auto select
* `padval` padding value for mosaic view; default `minimum(z)`
* `mosaic_npad` # of pixel padding for mosaic view; default 1
* `fft0` if true use FFTView to display; default false
* `title` for heatmap; default `""`
* `xlabel` for heatmap; default `""`
* `ylabel` for heatmap; default `""`
* `yflip` for heatmap; default `true` if minimum(y) >= 0
* `x` for x axis; default `1:size(z,1)`
* `y` for y axis; default `1:size(z,2)`
* `xtick` for heatmap; default `[minimum(x),maximum(x)]`
* `ytick` for heatmap; default `[minimum(y),maximum(y)]`

out
* returns plot handle

2019-02-23 Jeff Fessler, University of Michigan
"""
function jim(z;
	aspect_ratio = jim_def[:aspect_ratio],
	clim = jim_def[:clim],
	color = jim_def[:color],
	ncol::Integer = jim_def[:ncol],
	padval = jim_def[:padval],
	mosaic_npad::Integer = jim_def[:mosaic_npad],
	title = jim_def[:title],
	xlabel = jim_def[:xlabel],
	ylabel = jim_def[:ylabel],
	fft0::Bool = jim_def[:fft0],
	x = fft0 ? Int.(-size(z,1)/2:size(z,1)/2-1) : (1:size(z,1)),
	y = fft0 ? Int.(-size(z,2)/2:size(z,2)/2-1) : (1:size(z,2)),
	xtick = (minimum(x) < 0 && maximum(x) > 0) ?
		 [minimum(x),0,maximum(x)] : [minimum(x),maximum(x)],
	ytick = (minimum(y) < 0 && maximum(y) > 0) ?
		 [minimum(y),0,maximum(y)] : [minimum(y),maximum(y)],
	yflip::Bool = jim_def[:yflip] == nothing ? minimum(y) >= 0 : jim_def[:yflip],
	abswarn::Bool = jim_def[:abswarn],
	)

	if !isreal(z)
		if abswarn
			@warn "magnitude"
		end
		z = abs.(z)
	end

	if isempty(clim) # must wait until after possible abs() to do this
		clim = (minimum(z), maximum(z))
	end

	if isempty(padval) # must wait until after possible abs() to do this
		padval = minimum(z)
	end

	xy = (x, y)
	if ndims(z) > 2
		if ncol == 0
			ncol = Int(floor(sqrt(prod(size(z)[3:end]))))
		end
		z = mosaicview(z, padval, ncol=ncol, npad=mosaic_npad)
		xy = () # no x,y for mosaic
	elseif fft0
		z = FFTView(z)[x,y]
	end

heatmap(xy..., z', transpose=false,
	aspect_ratio=aspect_ratio,
	clim=clim,
	color=color,
	title=title,
	yflip=yflip,
	xlabel=xlabel,
	ylabel=ylabel,
	xtick=xtick,
	ytick=ytick)

end # jim

"""
`jim(z, title; kwargs...)`
"""
function jim(z, title::Union{String,LaTeXString}; kwargs...)
	return jim(z; title=title, kwargs...)
end

"""
`jim(x, y, z; kwargs...)`
"""
function jim(x, y, z; kwargs...)
	return jim(z; x=x, y=y, kwargs...)
end

"""
`jim(x, y, z, title; kwargs...)`
"""
function jim(x, y, z, title::Union{String,LaTeXString}; kwargs...)
	return jim(z; x=x, y=y, title=title, kwargs...)
end

# show docstring if user calls it with no arguments
function jim()
	@doc jim
end

#"""
#`jim(abswarn=false)`
#"""
#function jim(; abswarn::Bool=jim_state_abswarn)
#	global jim_state_abswarn = abswarn
#end

"""
`jim(key, value)` set default value for one of the keys
"""
function jim(key::Symbol, value)
	global jim_def
	if !haskey(jim_def, key)
		error("no key $key")
	end
	jim_def[key] = value
end

"""
`jim(:test)`

`jim(:keys)` return default keys

`jim(:defs)` return Dict of default keys / vals
"""
function jim(test::Symbol)
	global jim_def
	if test == :keys
		return keys(jim_def)
	end
	if test == :defs
		return jim_def
	end
	@assert test == :test
	jim(:abswarn, false)
	jim(ones(4,3), title="test2")
	jim(ones(4,3,5), title="test3")
	jim(1:4, 5:9, zeros(4,5), title="test3")
	jim(zeros(4,5), x=1:4, y=5:9, title="test3")
	jim(x=1:4, y=5:9, rand(4,5), title="test4")
	true
end


text:
struct Network

    lsizes::Tuple{Vararg{Int}};
    b::Tuple{Vararg{Vector{Float64}}};
    w::Tuple{Vararg{Matrix{Float64}}};
    z::Tuple{Vararg{Vector{Float64}}};
    a::Tuple{Vararg{Vector{Float64}}};
    depth::Int;
    inlength::Int;
    outlength::Int;

    function Network(lsizes::Tuple{Vararg{Int}})
        @assert length(lsizes) > 0 && all(lsizes .> 0)

        #Initializing the bias randomly
        #For the first layer the bias is not usable.
        b = Vector{Vector{Float64}}();
        push!(b,zeros(1));
        for l in 2:length(lsizes)
            push!(b,randn(lsizes[l]))
        end

        #Initializing the ws
        #For the first layer the weight is not usable.
        w = Vector{Matrix{Float64}}();
        push!(w, zeros(1,1));
        for l in 2:length(lsizes)
            push!(w,randn(lsizes[l],lsizes[l - 1]))
        end

        #The network store the values z and a of the last feedforward
        z = Vector{Vector{Float64}}();
        a = Vector{Vector{Float64}}();
        for l in 1:length(lsizes)
            push!(z, zeros(lsizes[l]));
            push!(a, zeros(lsizes[l]));
        end

        #Constructing
        new(lsizes, tuple(b...), tuple(w...), tuple(z...), tuple(a...),
            length(lsizes), first(lsizes), last(lsizes))
    end
end
Network(lsizes::Int...) = Network(lsizes);
Network(lsizes::Vector{Int}) = Network(lsizes...);
Base.size(net::Network) = net.lsizes;
Base.size(net::Network, l::Int) = net.lsizes[l];
output(net::Network) = last(net.a);
input(net::Network) = first(net.a);
outlength(net::Network) = length(output(net));
inlength(net::Network) = length(input(net));

function feedforward!(net::Network, input::Vector{Float64})::Vector{Float64}
    @assert net.inlength == length(input)

    #first layer
    net.z[1] .= input;
    net.a[1] .= input;

    #Other layers
    for l in 2:net.depth
        net.z[l] .= net.w[l]*net.a[l-1] + net.b[l];
        net.a[l] .= act_function.(net.z[l]);
    end

    #The output
    return output(net);
end

feedforward!(net::Network, testdata) = feedforward!(net, input(testdata));


text:
module NMRSpecifyRegions

# Write your package code here.

# dependencies to NMRSpectraSimulator. Need to add explicitly since these packages are not on the Julia public registry.
#import GISSMOReader # https://github.com/AI4DBiological-Systems/GISSMOReader.jl

#import NMRSpectraSimulator # https://github.com/AI4DBiological-Systems/NMRSpectraSimulator.jl
import LightGraphs

include("../src/types.jl")
include("../src/utils.jl")

include("../src/IO/string_helpers.jl")
include("../src/IO/config_helpers.jl")

include("../src/regions/setup.jl")
include("../src/regions/overlap.jl")

end


text:
# ---
# title: 564. Find the Closest Palindrome
# id: problem564
# author: zhwang
# date: 2022-03-02
# difficulty: Hard
# categories: String
# link: <https://leetcode.com/problems/find-the-closest-palindrome/description/>
# hidden: true
# ---
#
# Given an integer n, find the closest integer (not including itself), which is
# a palindrome.
#
# The 'closest' is defined as absolute difference minimized between two
# integers.
#
# **Example 1:**
#
#
#
#     Input: "123"
#     Output: "121"
#
#
# **Note:**
#
#   1. The input **n** is a positive integer represented by string, whose length will not exceed 18.
#   2. If there is a tie, return the smaller one as answer.
#
#
## @lc code=start
using LeetCode

function nearest_palindromic(num::String)::String
    function remove!(data::Vector{T}, target::T) where {T}
        return deleteat!(data, findall(==(target), data))
    end
    target_num, n = parse(Int, num), length(num)
    half_num = parse(Int, num[1:cld(n, 2)])
    nums = [half_num + i for i in -1:1] ## integers
    nums = [s * reverse(s[1:(n ÷ 2)]) for s in string.(nums)] ## strings
    nums = [parse(Int, i) for i in nums] ## integers
    if n > 1 ## special cases
        push!(nums, foldl((i, j) -> 10 * i + j, fill(9, n - 1))) ## 99...9
        push!(nums, foldl((i, j) -> 10 * i + j, (1, fill(0, n - 1)..., 1))) ## 10...01
    end
    sort!(remove!(nums, target_num))
    return string(nums[argmin([abs(i - target_num) for i in nums])])
end
## @lc code=end


text:
# This is a minimal working example for running LMlovo on several
# workers in Julia and saving them into a shared matrix

using Random
using Revise
using Distributed
using SharedArrays

@everywhere using RAFF
@everywhere using Base.CoreLogging

@everywhere CoreLogging._min_enabled_level[] = CoreLogging.Warn + 1

@everywhere gmodel!(x, t_, g) = begin
    g[1] = exp(t_ * x[2])
    g[2] = t_ * x[1] * exp(t_ * x[2]);
end

@everywhere model(x, t) = x[1] * exp(t * x[2])

n = 2

np = 1000

p = 700

# Set the seed for generating the same data
Random.seed!(123456789)

data, xSol = RAFF.generateNoisyData(model, n, np, p)

result = SharedArray{Float64, 2}(n, np)

f = @sync @distributed for i = 1:np

    # Starting point
    x = zeros(Float64, 2)

    # Call function and store results
    s, x, iter, p, f = LMlovo(model, gmodel!, x, data, 2, i, MAXITER=1000)

    result[:, i] .= x

    println("Finished. p = $(p) and f = $(f).")

end