text:
export setupDDParallel
export getWorkerForSubDomainMultiColor,solveDDParallel

function setupDDParallel(AT::Union{SparseMatrixCSC,DomainDecompositionOperatorConstructor},DDparam::DomainDecompositionParam,workerList::Array{Int64})
if workerList==[]
	return setupDDSerial(AT,DDparam);
else
	ActualWorkers = intersect(workerList,workers());
	if length(ActualWorkers)<length(workerList)
		warn("setupDDParallel: workerList included indices of non-existing workers.")
	end
end
numWorkers = length(ActualWorkers);

M 			= DDparam.Mesh;
numDomains 	= DDparam.numDomains;
overlap 	= DDparam.overlap;
n = M.n;

DDPreconditioners 	= Array{RemoteChannel}(undef,prod(numDomains));
DDparam.GlobalIndices = Array{Array{DDIndType}}(undef,prod(numDomains));

# ix = 1; nextidx() = (idx=ix; ix+=1; idx)

@sync begin
	for p = ActualWorkers
		@async begin
			# while true
				# idx = nextidx()
				# if idx > numWorkers
					# break
				# end
				for ii=1:length(DDPreconditioners)
					i = cs2loc(ii,numDomains);
					worker_index = getWorkerForSubDomainMultiColor(i,numDomains,numWorkers);

					# println("For cell ",i," we need worker ",ActualWorkers[worker_index],", and we have ",p)
					if ActualWorkers[worker_index] == p
						IIp = DDparam.getIndicesOfCell(numDomains,overlap, i,n);
						if isa(AT,SparseMatrixCSC)==true
							AI = sparse(AT[IIp,IIp]');
							DDPreconditioners[ii] = initRemoteChannel(performSetup,p,i, AI,  DDparam);
						else
							subparams = AT.getSubParams(AT.problem_param, M,i,numDomains,overlap);
							AI = DomainDecompositionOperatorConstructor(subparams,AT.getSubParams,AT.getOperator,AT.getDirichletMass);
							DDPreconditioners[ii] = initRemoteChannel(performSetup,p,i, AI,  DDparam);
						end
						IIp = convert(Array{DDIndType},IIp);
						DDparam.GlobalIndices[ii] = IIp;
						wait(DDPreconditioners[ii]);
					end
				end
			# end
		end
	end
end
DDparam.PrecParams = DDPreconditioners;
return DDparam;
end

function solveDDParallel(AT::SparseMatrixCSC,b::Array,x::Array,DDparam::DomainDecompositionParam,workerList,niter=1,doTranspose::Int64=0)
if workerList==[]
	return solveDD(AT,b,x,DDparam,niter,doTranspose);
else
	ActualWorkers = intersect(workerList,workers());
	if length(ActualWorkers)<length(workerList)
		error("setupDDParallel: workerList included indices of non-existing workers:",workerList,",",ActualWorkers)
	end
end
numWorkers = length(ActualWorkers);

numDomains 	= DDparam.numDomains;
overlap 	= DDparam.overlap;
dim 		= length(DDparam.numDomains);

for k=1:niter
	for color = 1:2^dim
	@sync begin
			# ix = 1; nextidx() = (idx=ix; ix+=1; idx)
			for p = ActualWorkers
				@async begin
					# while true
						# idx = nextidx()
						# if idx > numWorkers
							# break
						# end
						for ii=1:prod(numDomains)
							i = cs2loc(ii,numDomains);
							if ActualWorkers[getWorkerForSubDomainMultiColor(i,numDomains,numWorkers)] == p && cellColor(i)==color
								IIp = DDparam.GlobalIndices[ii];
								r = computeResidualAtIdx(AT,b,x,IIp);
								if DDparam.PrecParams[ii].where != p
									error("Sending to wrong worker!!!");
								end
								# println("before sending worker ",p," to solve subdomain",i,", for color ",cellColor(i))
								(t,DDparam.PrecParams[ii]) = remotecall_fetch(solveSubDomain, p, r, DDparam.PrecParams[ii],doTranspose);
								x[IIp] = x[IIp] + t;
								# println("after sending worker ",p," to solve subdomain",i,", for color ",cellColor(i))
							end
						end
					# end
				end
			end
		end
	end
	if niter > 1
		println(k,": ",norm(b - AT'*x)/norm(b))
	end
end
return x,DDparam;
end

function solveSubDomain(r,precRef::RemoteChannel,doTranspose::Int64)
	prec = take!(precRef);
	t = zeros(eltype(r),size(r));
	(t,prec.Ainv) = solveLinearSystem!(prec.A_i,r,t,prec.Ainv,doTranspose);
	precRef = put!(precRef,prec)
	return (t,precRef)
end

function getWorkerForSubDomainMultiColor(iiloc::Array{Int64},NumCells,NumWorkers)
	NumCellsWorkers = div.(NumCells.+1,2);
	iiploc = div.(iiloc.+1,2);
	iip = loc2cs(iiploc,NumCellsWorkers);
	p = mod(iip-1,NumWorkers)+1;
	return p;
end


text:
#hello world
println("hello world")

print("hello")
println(" world")

ext: jl

max_stars_repo_name: KorfLab/julie

max_issues_repo_name: KorfLab/julie

max_forks_repo_name: KorfLab/julie


text:
export
    Factor,
    FactorObject,
    FactorGraphStruct,
    FactorGraph,
    assocvariables,
    assocfactors,
    linkedfactors,
    chainstruct,
    chaingraph

"""
    Factor

Encapsulation of a factor, made out of the function that computes the next
bounce, the gradient of the log likelihood associated with the factor and
the index of the factor in the list of factors of a factor graph.
See also: `FactorGraph`, `FactorGraphStruct`.
"""
struct Factor
    nextevent::Function # compute the corresponding first arrival time
    gll::Function       # gradient of log likelihood
    index::Int          # factor index
end

"""
    FactorGraphStruct

Encapsulation of the structure of a factor graph. It is made of a list of
factors each corresponding to a list of variables attached to the factor.
See also: `Factor`, `FactorGraph`.
"""
struct FactorGraphStruct
    flist::Vector{Vector{Int}} # each entry is a list of vars
    # -- implicit
    vlist::Vector{Vector{Int}} # each entry is a list of factors
    nfactors::Int
    nvars::Int
    # -- constructor
    function FactorGraphStruct(flist)
        # finding how many variables there are
        maxvar = 0
        for f in flist
            maxvar = max(maxvar, maximum(f))
        end
        # populating the lists of factors for each variable
        vars = [Vector{Int}(undef, 0) for i in 1:maxvar]
        for (fi, f) in enumerate(flist)
            for e in f
                push!(vars[e], fi)
            end
        end
        new(flist, vars, length(flist), maxvar)
    end
end

"""
    FactorGraph

Encapsulation of a factor graph. It is made out of a structure
(FactorGraphStruct) and an array of factors corresponding to the given
structure.
See also: `Factor`, `FactorGraphStruct`.
"""
struct FactorGraph
    structure::FactorGraphStruct
    factors::Vector{Factor}
end
FactorGraph(fgs::Vector{Vector{Int}}, factors::Vector{Factor}) =
    FactorGraph(FactorGraphStruct(fgs), factors)

# access the variables associated to a factor or vice versa
assocvariables(fgs::FactorGraphStruct, fidx::Int) = fgs.flist[fidx]
assocfactors(fgs::FactorGraphStruct, vidx::Int) = fgs.vlist[vidx]
assocvariables(fg::FactorGraph, fidx::Int) = assocvariables(fg.structure, fidx)
assocfactors(fg::FactorGraph, vidx::Int) = assocfactors(fg.structure, vidx)

"""
    linkedfactors(fgs, fidx)

Return the list of factors that have intersecting variable set with the factor
indexed by `fidx`. It does *not* contain the current factor fidx for efficiency
reason (the current factor can be updated quicker directly from the general
simulation loop).
"""
function linkedfactors(fgs::FactorGraphStruct, fidx::Int)
    lf = Vector{Int}(undef, 0)
    for v ∈ fgs.flist[fidx]
        push!(lf, fgs.vlist[v]...)
    end
    unq = unique(lf)                    # this still contains fidx
    deleteat!(unq, findall((in)(fidx), unq))   # this does not
end
linkedfactors(fg::FactorGraph, fidx::Int) = linkedfactors(fg.structure, fidx)

# ==============================================================================
# Definition of standard FG structures

chainstruct(nvars::Int) = FactorGraphStruct([[i,i+1] for i ∈ 1:nvars-1])
chaingraph(lf::Vector{Factor}) = FactorGraph(chainstruct(length(lf) + 1), lf)

# row based reading (left right then top down)
# gridstruct(nvars::Int) = error("Not implemented yet")


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "TtH"
version = v"0.1.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("http://hutchinson.belmont.ma.us/tth/tth-noncom/tth_C.tgz", "83c1f39fbf3377fb43e3d01d042302fa91f8758aa9acc10e22fe8af140f0126c")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/tth_C
mkdir -p ${bindir}
gcc -o ${bindir}/tth${exeext} tth.c
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    ExecutableProduct("tth", :tth),
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "LibVPX"
version = v"1.8.0"

# Collection of sources required to build LibVPX
sources = [
    "https://github.com/webmproject/libvpx/archive/v1.8.0.tar.gz" =>
    "86df18c694e1c06cc8f83d2d816e9270747a0ce6abe316e93a4f4095689373f6",

    "./patches"
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/libvpx-*/
sed -i 's/cp -p/cp/' build/make/Makefile
patch -p1 < $WORKSPACE/srcdir/macos.patch
mkdir libvpx-build
cd libvpx-build
apk add diffutils yasm

export CONFIG_OPTS="--enable-shared --disable-static"
export TARGET=generic-gnu
if [[ "${target}" == i686-linux-* ]]; then
    export CROSS=x86-linux-gcc
elif [[ "${target}" == x86_64-linux-* ]]; then
    export CROSS=x86_64-linux-gcc
elif [[ "${target}" == arm-linux-* ]]; then
    export CROSS=armv7-linux-gcc
elif [[ "${target}" == powerpc64le-linux-* ]]; then
    export CROSS=ppc64le-linux-gcc
elif [[ "${target}" == x86_64-apple-* ]]; then
    export CROSS=x86_64-darwin14-gcc
    export CC=/opt/${target}/bin/${target}-gcc
    export CXX=/opt/${target}/bin/${target}-g++
    export LD=/opt/${target}/bin/${target}-ld
    export TARGET=$CROSS
elif [[ "${target}" == i686-w64-mingw32 ]]; then
    export CROSS=x86-win32-gcc
    export CONFIG_OPTS=""
    export TARGET=$CROSS
elif [[ "${target}" == x86_64-w64-mingw32 ]]; then
    export CROSS=x86_64-win64-gcc
    export CONFIG_OPTS=""
    export TARGET=$CROSS
elif [[ "${target}" == x86_64-unknown-* ]]; then
    export CROSS=x86_64-unknown-gcc
    export CC=/opt/${target}/bin/${target}-gcc
    export CXX=/opt/${target}/bin/${target}-g++
    export LD=/opt/${target}/bin/${target}-ld
    export CONFIG_OPTS="--enable-shared --disable-static --disable-multithread"
elif [[ "${target}" == aarch64-linux-* ]]; then
    export CROSS=aarch64-linux-gcc
else
    export CROSS=$target
fi

../configure --prefix=$prefix --target=${TARGET} --as=yasm ${CONFIG_OPTS}
echo "SRC_PATH_BARE=.." >> config.mk
echo "target=libs" >> config.mk
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libvpx", :libvpx)
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
#=>
     Julia can be yet more snappy when register-focused use
        (bitshifts, offsets) keeps counts in Int32 registers.
     Advantages are entirely system dependant, often found.
     The functions in this file allow us to choose the Signed
        type used for indexing most parsimoneously -- and here
        we use the smallest system integer type that suffices.
<=#

bitsof(::Type{T}) where T = sizeof(T) << 3

const TYPES_FOR_INDEXING = [ Int8, Int16, Int32, Int64, Int128 ];

function type_for_indexing(posint::T) where T<:Signed
    posint <= 0 && throw(ErrorException("posint ($(posint)) must be > 0"))

    bits_required = bitsof(T) - leading_zeros(posint)
    type_index = nextpow2(bits_required + 1) >>> 2
    type_index = trailing_zeros(type_index)
    type_index = max(1, min(5, type_index))

    return TYPES_FOR_INDEXING[ type_index ]
end


text:
# Conradi, C., Shiu, A.,
# Dynamics of post-translational modification systems: recent progress and future directions
# Eq. 3.4
using Logging

using StructuralIdentifiability

logger = Logging.SimpleLogger(stdout, Logging.Info)
global_logger(logger)

ode = @ODEmodel(
    x1'(t) = -k1 * x1(t) * x2(t) + k2 * x4(t) + k4 * x6(t),
    x2'(t) = -k1 * x1(t) * x2(t) + k2 * x4(t) + k3 * x4(t),
    x3'(t) = k3 * x4(t) + k5 * x6(t) - k6 * x3(t) * x5(t),
    x4'(t) = k1 * x1(t) * x2(t) - k2 * x4(t) - k3 * x4(t),
    x5'(t) = k4 * x6(t) + k5 * x6(t) - k6 * x3(t) * x5(t),
    x6'(t) = -k4 * x6(t) - k5 * x6(t) + k6 * x3(t) * x5(t),
    y1(t) = x3(t),
    y2(t) = x2(t)
)

@time println(assess_identifiability(ode))


text:
__precompile__()
module ForSyDeIO

include("ForSyDeIO/Models.jl")
include("ForSyDeIO/Drivers.jl")

using ForSyDeIO.Drivers

export load_model, write_model

# Write your package code here.

end

max_stars_repo_name: forsyde/forsyde-io

max_issues_repo_name: forsyde/forsyde-io

max_forks_repo_name: forsyde/forsyde-io


text:
"""
    setup_growth(Ωm, ΩΛ[, a1=1e-2])

Return the solutions of coupled differential equations given in `growth!`,
with the initial condition such that ``δ(a) = a`` at `a = a1`. The default value is `a1 = 1e-2`.

# Arguments
- `Ωm::Real`: present-day total matter density parameter.
- `ΩΛ::Real`: present-day dark energy density parameter (for cosmological constant).

# Optional arguments
- `a1::Real=1e-2`: initial scale factor, at which the initial condition is set as δ(a1) = a1.

# Examples
```julia-repl
julia> Ωm, ΩΛ = 0.3, 0.7
(0.3, 0.7)

julia> sol = MatterPower.setup_growth(Ωm, ΩΛ)
retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 14-element Array{Float64,1}:
 0.01
 0.022952118470240847
 0.03152976635130906
 0.04919933323749941
 ⋮
 0.7292247303730331
 0.994509796947469
 1.0
u: 14-element Array{Array{Float64,1},1}:
 [0.01, -0.05477231965144438]
 [0.022953661875878814, -0.08298398087457971]
 [0.0315315317884963, -0.09726149943169236]
 [0.049200231832575124, -0.12149028570497]
 ⋮
 [0.6428440417117532, -0.39916939875846624]
 [0.7768316486186253, -0.3997506971222691]
 [0.7790350080966789, -0.3994861471523371]

julia> redshift = 1.2
1.2

julia> a = 1/(1+redshift)
0.45454545454545453

julia> δ = sol(a)[1]
0.43811688788577136

julia> θ = sol(a)[2]
-0.3526347043144974
```
"""
function setup_growth(Ωm::Real, ΩΛ::Real, a1::Real=1e-2)
    Ωk = 1 - Ωm - ΩΛ
    E1 = √(Ωm / a1^3 + Ωk / a1^2 + ΩΛ)
    u0 = [a1; -a1^2 * E1]
    tspan = (a1, 1.0)
    prob = ODEProblem(growth!, u0, tspan, [Ωm, ΩΛ])
    sol = solve(prob, Tsit5())
end

"""
    growth!(du, u, p, a)

Coupled differential equations to obtain the growth factor of linear density fluctuations.

1. Continuity Equation, ``δ' + θ = 0``

2. Euler equation, ``θ' + (a'/a)θ - k^2 ϕ = 0``
3. Poisson equation, ``k^2ϕ = -4πG a^2 ρm δ``

where θ = div(velocity), ϕ is the Newtonian gravitational potential, and the primes denote conformal time derivatives.

To solve these equations, we rewrite them using derivatives with respect to the scale factor, `a`.
1. `du[1]` = ``dδ/da = -θ/a^2/E(a)``
2. `du[2]` = ``dθ/da = -θ/a + q^2ϕ/a^2/E(a)``
3. ``q^2ϕ = -(3/2)Ωm δ/a``

where

- `u[1]` = δ, `u[2]` = θ = div(velocity)/H0
- `p[1]` = Ωm, `p[2]` = ΩΛ
- ``E(a) = H(a)/H0 = √(Ωm/a^3 + Ωk/a^2 + ΩΛ)``
- ``q = k/H0``

These equations will be solved by Julia's ODE solver as

`ODEProblem(growth!, u0, tspan, p)`

See this [documentation](https://github.com/SciML/OrdinaryDiffEq.jl) for how to use this ODE solver.
"""
function growth!(du, u, p, a)
    δ, θ = u
    Ωm, ΩΛ = p
    Ωk = 1 - Ωm - ΩΛ
    E = √(Ωm / a^3 + Ωk / a^2 + ΩΛ) # E(a) = H(a)/H0
    q2ϕ = -(3 / 2) * Ωm * δ / a # Poisson equation
    du[1] = dδda = -θ / a^2 / E # Continuity equation
    du[2] = dθda = -θ / a + q2ϕ / a^2 / E # Euler equation
end


text:
title = "LyceumBase"
weight = 5


text:
using ConditionalJuMP
using ConditionalJuMP: switch!,  _getvalue, UnhandledComplementException, hascomplement, getindicator!, disjunction!, isjump, Conditional
using IntervalArithmetic: Interval
using JuMP
using Cbc
using Test
using Random
using Pkg
using Nullables

# Macro hygiene
module MacroHygieneTest
    using ConditionalJuMP: @?
    using Test

    function f(x)
        x + 1
    end

    @testset "macro hygiene" begin
        x = 5
        @test @?(f(x) <= 0) == false
    end
end

@testset "ConditionalJuMP" begin
    @testset "simplification" begin
        test_cases = [
            ("small amount of variables", 100, 20, 1000),
            ("large amount of variables", 10000, 2000, 100),
            ("huge amount of variables", 100000, 20000, 1),
        ]
        for (test_name, num_vars_total, num_vars_in_expr, num_trials) in test_cases
            @testset "$(test_name)" begin
                m = Model()
                Random.seed!(42)
                @variable m q[1:num_vars_total]
                for i in 1:num_trials
                    N = rand(1:num_vars_in_expr)
                    x = randn(N)' * rand(q, N)
                    s1 = copy(x)
                    s2 = copy(x)
                    ConditionalJuMP.simplify_dict!(s1)
                    ConditionalJuMP.simplify_inplace!(s2)
                    @test s1.constant == s2.constant
                    @test length(s1.vars) == length(s2.vars)
                    @test length(s1.coeffs) == length(s2.coeffs)
                    @test length(s1.vars) == length(unique(s1.vars))

                    s1_sorted = sort(collect(zip(s1.vars, s1.coeffs)), by=p -> (p[1].col, p[2]))
                    s2_sorted = sort(collect(zip(s2.vars, s2.coeffs)), by=p -> (p[1].col, p[2]))
                    for (v1, v2) in zip(s1_sorted, s2_sorted)
                        @test v1[1] == v2[1]
                        @test v1[2] ≈ v2[2]
                        # directly testing equality of the two sorted lists is too strict and leads to spurious errors.
                    end
                    for v in s1.vars
                        @test count(x -> x == v, s1.vars) == 1
                    end
                end
            end
        end
    end

    @testset "bounds" begin
        m = Model()
        @variable m 1 <= x <= 3
        e = 2 * x + 1
        @test upperbound(e) == 7
        @test lowerbound(e) == 3
        @test upperbound(e - 2x) == 1
        @test lowerbound(e - 2x) == 1
        @test lowerbound(AffExpr(2)) == 2
        @test upperbound(AffExpr(2)) == 2

        i = interval(e, false)
        @test i == Interval(3, 7)
        @test lowerbound(i) == lowerbound(e)
        @test upperbound(i) == upperbound(e)

        i = interval(5.0)
        @test lowerbound(i) == 5.0
        @test upperbound(i) == 5.0

        i = interval(x)
        @test lowerbound(i) == 1
        @test upperbound(i) == 3
    end

    @testset "simple model" begin
        m = Model()
        @variable m -1 <= x <= 10
        c1 = @?(x <= 0)
        c2 = @?(x >= 0)
        c3 = @?(x + 5 <= 5)
        c4 = @?(x + 1 == 10)

        print(IOBuffer(), c1)
        print(IOBuffer(), c2)
        print(IOBuffer(), c1 & c2)

        @testset "equality" begin
            @test c1 != c2
            @test c2 != c3
            @test c1 == c3
            @test c1 !== c3
        end

        @testset "intersect" begin
            @test (c1 & c3) == c1
            @test length(c1 & c2) == 2
        end

        @testset "complement" begin
            @test !c1 == c2
            @test !c2 == c1
            @test !hascomplement(c4)
            @test_throws UnhandledComplementException !c4
        end

        @testset "hashing" begin
            @test hash(c1) == hash(c3)
            @test hash(c1) != hash(c2)
        end

        @testset "value" begin
            @test isnull(_getvalue(c1))
            @test isnull(_getvalue(c2))
            setvalue(x, -0.1)
            @test get(_getvalue(c1)) == -0.1
            @test get(_getvalue(c2)) == 0.1
        end

        @testset "bounds" begin
            @test lowerbound(x) == -1
            @test upperbound(x) == 10
            @test lowerbound(x + 1) == 0
            @test upperbound(x + 10) == 20
        end

        @testset "indicators" begin
            z1 = getindicator!(m, c1)
            z2 = getindicator!(m, c2)
            z3 = getindicator!(m, c3)
            @test z1 === z3
            @test z1 == !z2
        end

        @testset "isjump" begin
            @test isjump(x)
            @test !isjump(1)
            @test isjump(2x)
            @test isjump(c1)
            @test isjump(c1 & c2)
            @test isjump(c1 => c2)
            @test isjump(c1 => Conditional())
        end
    end

    @testset "disjunction" begin
        m = Model(solver=CbcSolver(logLevel=0))
        @variable m -5 <= x <= 5
        c1 = @?(x >= -4) & @?(x <= -2)
        c2 = @?(x >= 1) & @?(x <= 3)
        @objective m Min x
        disjunction!(m, (c1, c2))
        solve(m)
        @test getvalue(x) ≈ -4
        @objective m Max x
        solve(m)
        @test getvalue(x) ≈ 3
        @constraint m x <= 0
        solve(m)
        @test getvalue(x) ≈ -2
    end

    @testset "implies" begin
        @testset "min" begin
            m = Model(solver=CbcSolver(logLevel=0))
            @variable m -5 <= x <= 5
            @variable m -5 <= y <= 5
            @implies m x ≤ 0 => y == 3
            @implies m x ≥ 0 => y == 1
            @objective m Min x

            # test that the two implications are being represented
            # by a single binary variable
            @test sum(m.colCat .== :Bin) == 1
            solve(m)
            @test getvalue(x) ≈ lowerbound(x)
            @test getvalue(y) ≈ 3
        end

        @testset "max" begin
            m = Model(solver=CbcSolver(logLevel=0))
            @variable m -5 <= x <= 5
            @variable m -5 <= y <= 5
            @implies m x ≤ 0 => y == 3
            @implies m x ≥ 0 => y == 1
            @objective m Max x

            # test that the two implications are being represented
            # by a single binary variable
            @test sum(m.colCat .== :Bin) == 1
            solve(m)
            @test getvalue(x) ≈ upperbound(x)
            @test getvalue(y) ≈ 1
        end

        @testset "non-disjoint" begin
            m = Model(solver=CbcSolver(logLevel=0))
            @variable m -5 <= x <= 5
            @variable m -5 <= y <= 5
            @implies m x ≤ -1 => y == 3
            @implies m x ≥ 1 => y == 1
            @objective m Max x
            @test sum(m.colCat .== :Bin) == 2
            solve(m)
            @test getvalue(x) ≈ upperbound(x)
            @test getvalue(y) ≈ 1
        end
    end

    @testset "vector" begin
        @testset "min" begin
            m = Model(solver=CbcSolver(logLevel=0))
            @variable m -5 <= x <= 5
            @variable m -5 <= y[1:2] <= 5
            @implies m x ≤ 0 => y .== [3, 3.5]
            @objective m Min x
            @test sum(m.colCat .== :Bin) == 1
            solve(m)
            @test getvalue(x) ≈ lowerbound(x)
            @test getvalue.(y) ≈ [3, 3.5]
        end

        @testset "max" begin
            m = Model(solver=CbcSolver(logLevel=0))
            @variable m -5 <= x <= 5
            @variable m -5 <= y[1:2] <= 5
            @implies m x ≥ 0 => y .== [1, -1]
            @objective m Max x
            @test sum(m.colCat .== :Bin) == 1
            solve(m)
            @test getvalue(x) ≈ upperbound(x)
            @test getvalue.(y) ≈ [1, -1]
        end
    end

    @testset "dynamics" begin
        function update(x)
            @ifelse(x <= 0, 1, -1)
        end

        @testset "mixed integer" begin
            m = Model(solver=CbcSolver(logLevel=0))
            @variable m -0.5 <= x <= 0.5

            ys = [x]
            for i in 1:3
                push!(ys, update(ys[end]))
            end

            @objective m Max sum(ys)
            @test sum(m.colCat .== :Bin) == 3
            solve(m)
            @test getvalue.(ys) ≈ [0, 1, -1, 1]
        end

        @testset "seeded positive" begin
            m = Model(solver=CbcSolver(logLevel=0))
            @variable m -0.5 <= x <= 0.5

            ys = [x]
            for i in 1:3
                push!(ys, update(ys[end]))
            end

            @objective m Max sum(ys)
            setvalue(x, 0.5)
            warmstart!(m, true)
            @test sum(m.colCat .== :Bin) == 0
            solve(m)
            @test getvalue.(ys) ≈ [0.5, -1, 1, -1]
        end

        @testset "seeded negative" begin
            m = Model(solver=CbcSolver(logLevel=0))
            @variable m -0.5 <= x <= 0.5

            ys = [x]
            for i in 1:3
                push!(ys, update(ys[end]))
            end

            @objective m Max sum(ys)
            setvalue(x, -0.5)
            warmstart!(m, true)
            @test sum(m.colCat .== :Bin) == 0
            solve(m)
            @test getvalue.(ys) ≈ [0, 1, -1, 1]
        end

        @testset "vector" begin
            function update2(x)
                @ifelse(x[1] <= 0,
                    [1, 2],
                    [3, 4])
            end

            m = Model(solver=CbcSolver(logLevel=0))
            @variable(m, -10 <= x[1] <= 10)
            y = update2(x)
            @constraint(m, x[1] == -5)
            solve(m)
            @test getvalue(y) ≈ [1, 2]
        end
    end

    function model_latex(m::Model)
        io = IOBuffer()
        show(io, MIME"text/latex"(), m)
        seekstart(io)
        read(io, String)
    end

    @testset "fixing and unfixing" begin
        m = Model(solver=CbcSolver(logLevel=0))
        @variable m -1 <= x <= 1
        @variable m -1 <= y <= 1
        @implies m x ≤ 0 => y == 0.5
        setvalue(x, 0.5)
        @objective(m, Min, x + 0.1y)
        out1 = model_latex(m)

        # Use the value of x as a hint
        warmstart!(m)
        @test sum(m.colCat .== :Bin) == 1
        @test model_latex(m) == out1
        solve(m)
        @test getvalue(x) ≈ -1
        @test getvalue(y) ≈ 0.5

        # Now use the current value of x to fix the indicators
        setvalue(x, 0.5)
        warmstart!(m, true)
        @test sum(m.colCat .== :Bin) == 0
        solve(m)
        @test getvalue(x) ≈ 0
        @test getvalue(y) ≈ -1

        # Now un-fix the indicators and make sure we get the original
        # model back
        setvalue(x, 0.5)
        warmstart!(m)
        @test sum(m.colCat .== :Bin) == 1
        @test model_latex(m) == out1
        solve(m)
        @test getvalue(x) ≈ -1
        @test getvalue(y) ≈ 0.5
    end

    @testset "disjunctions" begin
        @testset "simple three-way disjunction" begin
            m = Model(solver=CbcSolver(logLevel=0))
            @variable m -1 <= x <= 1
            @variable m -1 <= y <= 1
            c1 = @?(x <= -0.5)
            c3 = @?(x >= 0.5)
            c2 = !c1 & !c3
            @implies(m,
                c1 => nothing,
                c2 => y == -0.5,
                c3 => nothing)
            @constraint(m, x == 0)
            warmstart!(m)
            solve(m)
            @test getvalue(x) ≈ 0
            @test getvalue(y) ≈ -0.5
        end

        @testset "two-way disjunction" begin
            m = Model(solver=CbcSolver(logLevel=0))
            @variable m -1 <= x <= 1
            @variable m -1 <= y <= 1
            c1 = @?(x <= -0.5)
            c2 = @?(x >= 0.5)
            @implies(m,
                c1 => y == 0.1,
                c2 => nothing
                )
            # Test that the two-way disjunction is represented with just one
            # binary variable
            @test sum(m.colCat .== :Bin) == 1
        end

        @testset "missing disjunction" begin
            m = Model(solver=CbcSolver(logLevel=0))
            @variable m -1 <= x <= 1
            @variable m -1 <= y <= 1
            c1 = @?(x <= -0.5)
            c3 = @?(x >= 0.5)
            c2 = !c1 & !c3
            @test_throws(
                ConditionalJuMP.UnhandledComplementException,
                @implies(m, c2 => y == -0.5)
                )
        end

        @testset "disjunction with switch" begin
            @testset "case 1" begin
                m = Model(solver=CbcSolver(logLevel=0))
                @variable m -1 <= x <= 1
                c1 = @?(x <= -0.5)
                c3 = @?(x >= 0.5)
                c2 = !c1 & !c3
                y = switch!(m, c1=>0.1, c2=>0.2, c3=>0.3)
                @constraint(m, x == -1)
                warmstart!(m)
                solve(m)
                @test getvalue(x) ≈ -1
                @test getvalue(y) ≈ 0.1
            end
            @testset "case 2" begin
                m = Model(solver=CbcSolver(logLevel=0))
                @variable m -1 <= x <= 1
                c1 = @?(x <= -0.5)
                c3 = @?(x >= 0.5)
                c2 = !c1 & !c3
                y = switch!(m, c1=>0.1, c2=>0.2, c3=>0.3)
                @constraint(m, x == -0)
                warmstart!(m)
                solve(m)
                @test getvalue(x) ≈ 0
                @test getvalue(y) ≈ 0.2
            end
            @testset "case 3" begin
                m = Model(solver=CbcSolver(logLevel=0))
                @variable m -1 <= x <= 1
                c1 = @?(x <= -0.5)
                c3 = @?(x >= 0.5)
                c2 = !c1 & !c3
                y = switch!(m, c1=>0.1, c2=>0.2, c3=>0.3)
                @constraint(m, x == 1)
                warmstart!(m)
                solve(m)
                @test getvalue(x) ≈ 1
                @test getvalue(y) ≈ 0.3
            end
        end

        @testset "disjunction with switch macro" begin
            function f(x)
                @switch(
                    (x <= -0.5) => 0.1,
                    ((x >= -0.5) & (x <= 0.5)) => 0.2,
                    (x >= 0.5) => 0.3
                )
            end

            @test @inferred(f(-1)) == 0.1
            @test @inferred(f(0)) == 0.2
            @test @inferred(f(1)) == 0.3

            @testset "case 1" begin
                m = Model(solver=CbcSolver(logLevel=0))
                @variable m -1 <= x <= 1
                y = f(x)
                @constraint(m, x == -1)
                warmstart!(m)
                solve(m)
                @test getvalue(x) ≈ -1
                @test getvalue(y) ≈ 0.1
            end
            @testset "case 2" begin
                m = Model(solver=CbcSolver(logLevel=0))
                @variable m -1 <= x <= 1
                y = f(x)
                @constraint(m, x == -0)
                warmstart!(m)
                solve(m)
                @test getvalue(x) ≈ 0
                @test getvalue(y) ≈ 0.2
            end
            @testset "case 3" begin
                m = Model(solver=CbcSolver(logLevel=0))
                @variable m -1 <= x <= 1
                y = f(x)
                @constraint(m, x == 1)
                warmstart!(m)
                solve(m)
                @test getvalue(x) ≈ 1
                @test getvalue(y) ≈ 0.3
            end
        end

        @testset "disjunctions with arrays" begin
            function g(x)
                @switch(
                    (x <= -0.5) => [0.1, 1.1],
                    ((x >= -0.5) & (x <= 0.5)) => [0.2, 2.2],
                    (x >= 0.5) => [0.3, 3.3]
                )
            end

            @test @inferred(g(-1)) == [0.1, 1.1]
            @test @inferred(g(0)) == [0.2, 2.2]
            @test @inferred(g(1)) == [0.3, 3.3]

            @testset "case 1" begin
                m = Model(solver=CbcSolver(logLevel=0))
                @variable m -1 <= x <= 1
                y = g(x)
                @constraint(m, x == -1)
                warmstart!(m)
                solve(m)
                @test getvalue(x) ≈ -1
                @test getvalue(y) ≈ [0.1, 1.1]
            end
            @testset "case 2" begin
                m = Model(solver=CbcSolver(logLevel=0))
                @variable m -1 <= x <= 1
                y = g(x)
                @constraint(m, x == -0)
                warmstart!(m)
                solve(m)
                @test getvalue(x) ≈ 0
                @test getvalue(y) ≈ [0.2, 2.2]
            end
            @testset "case 3" begin
                m = Model(solver=CbcSolver(logLevel=0))
                @variable m -1 <= x <= 1
                y = g(x)
                @constraint(m, x == 1)
                warmstart!(m)
                solve(m)
                @test getvalue(x) ≈ 1
                @test getvalue(y) ≈ [0.3, 3.3]
            end
        end

        @testset "simple 1-D example" begin
            m = Model(solver=CbcSolver(logLevel=0))
            @variable(m, -5 <= x <= 5)
            @disjunction(m, ((x >= -2) & (x <= -1)), ((x >= 1.5) & (x <= 2)))
            # `s` behaves like a poor-man's quadratic objective
            @variable(m, s)
            @constraints m begin
                s >= x
                s >= -x
                s >= 2x - 1
                s >= -2x - 1
            end
            @objective(m, Min, s)
            solve(m)
            @test getvalue(x) ≈ -1

            setvalue(x, 2)
            warmstart!(m, true)
            solve(m)
            @test getvalue(x) ≈ 1.5
        end

        @testset "complementarity" begin
            m = Model(solver=CbcSolver(logLevel=0))
            @variable(m, -1 <= x <= 1)
            @variable(m, -1 <= y <= 1)
            @disjunction(m, x == 0, y == 0)
            @objective(m, Min, x + 2y)
            solve(m)
            @test getvalue(x) ≈ 0
            @test getvalue(y) ≈ -1

            # Check that we can warmstart even when neither case is precisely satisfied
            setvalue(y, 0.01)
            setvalue(x, -0.5)
            warmstart!(m, true)
            solve(m)
            @test getvalue(x) ≈ -1
            @test getvalue(y) ≈ 0
        end

        @testset "disjunction with one entry" begin
            m = Model(solver=CbcSolver(logLevel=0))
            @variable m -1 <= x <= 1
            @disjunction m x >= 0.25
            @objective m Min x
            solve(m)
            @test getvalue(x) ≈ 0.25
        end

        @testset "disjunction with overlapping regions" begin
            m = Model(solver=CbcSolver(logLevel=0))
            @variable m -2 <= x[1:2] <= 2
            @disjunction(m, x[2] >= 0.5, x[2] <= -0.5, x[1] <= 1)
            @constraint m x[1] == 0.8
            @constraint m x[2] == 0.6
            @test solve(m) == :Optimal
            @test getvalue(x) ≈ [0.8, 0.6]
        end

        @testset "test that warmstart does not introduce new variables" begin
            m = Model(solver=CbcSolver(logLevel=0))
            @variable m -1 <= x <= 1
            @disjunction m (x == 0.5) (x == -0.5)
            @objective m Min x
            @test length(m.colCat) == 2
            warmstart!(m, false)
            @test length(m.colCat) == 2
            solve(m)
            @test getvalue(x) ≈ -0.5
            warmstart!(m, true)
            @test length(m.colCat) == 2
            solve(m)
            @test getvalue(x) ≈ -0.5
        end
    end

    @testset "examples" begin
        @testset "friction with @switch" begin
            include("../examples/friction.jl")
        end

        @testset "block with wall" begin
            if haskey(Pkg.installed(), "Gurobi")
                include("../examples/block_with_wall.jl")
                run_mpc(State(1.0, -1.0), 10)
            end
        end

        @testset "complementarity" begin
            include("../examples/complementarity.jl")
        end

        @testset "polyhedra" begin
            include("../examples/polyhedra.jl")
        end
    end

    @testset "constant objective" begin
        m = Model(solver=CbcSolver())
        @variable m x >= 0
        @objective m Min x - 1
        ConditionalJuMP.handle_constant_objective!(m)
        @test getobjective(m).aff.constant == 0
        @test sum(m.colCat .== :Cont) == 1
        @test sum(m.colCat .== :Fixed) == 1
        solve(m)
        @test getvalue(x) ≈ 0
    end
end


text:
using BudAutomata
using Test

@testset "BudAutomata.jl" begin
    # Write your own tests here.
end


text:
using TikzPictures
using Base.Test

using TikzPictures
tp = TikzPicture("\\draw (0,0) -- (10,10);\n\\draw (10,0) -- (0,10);\n\\node at (5,5) {tikz \$\\sqrt{\\pi}\$};", options="scale=0.25", preamble="")
save(TEX("test"), tp)
if success(`lualatex -v`)
    save(PDF("test"), tp)
    save(SVG("test"), tp)
end


text:
using SimJulia

function test_callback(ev::Event)
  println("I am a callback function running in $(typeof(environment(ev)))")
end

sim = Simulation()
ev = Event(sim)
println(typeof(ev))
cb = @callback test_callback(ev)
remove_callback(cb, ev)
println(state(ev))
show(value(ev))
println()
succeed(ev, value="Hi")
@callback test_callback(ev)
println(state(ev))
println(value(ev))
run(sim)
try
  @callback test_callback(ev)
catch exc
  println("$exc has been thrown!")
end
println(state(ev))
println(value(ev))


text:
doc"""
    Exponential(θ)

The *Exponential distribution* with scale parameter `θ` has probability density function

$f(x; \theta) = \frac{1}{\theta} e^{-\frac{x}{\theta}}, \quad x > 0$

```julia
Exponential()      # Exponential distribution with unit scale, i.e. Exponential(1)
Exponential(b)     # Exponential distribution with scale b

params(d)          # Get the parameters, i.e. (b,)
scale(d)           # Get the scale parameter, i.e. b
rate(d)            # Get the rate parameter, i.e. 1 / b
```

External links

* [Exponential distribution on Wikipedia](http://en.wikipedia.org/wiki/Exponential_distribution)

"""
immutable Exponential{T<:Real} <: ContinuousUnivariateDistribution
    θ::T		# note: scale not rate

    Exponential(θ::Real) = (@check_args(Exponential, θ > zero(θ)); new(θ))
end

Exponential{T<:Real}(θ::T) = Exponential{T}(θ)
Exponential(θ::Integer) = Exponential(Float64(θ))
Exponential() = Exponential(1.0)

@distr_support Exponential 0.0 Inf

### Conversions
convert{T <: Real, S <: Real}(::Type{Exponential{T}}, θ::S) = Exponential(T(θ))
convert{T <: Real, S <: Real}(::Type{Exponential{T}}, d::Exponential{S}) = Exponential(T(d.θ))

#### Parameters

scale(d::Exponential) = d.θ
rate(d::Exponential) = 1 / d.θ

params(d::Exponential) = (d.θ,)
@inline partype{T<:Real}(d::Exponential{T}) = T

#### Statistics

mean(d::Exponential) = d.θ
median(d::Exponential) = logtwo * d.θ
mode{T<:Real}(d::Exponential{T}) = zero(T)

var(d::Exponential) = d.θ^2
skewness{T<:Real}(d::Exponential{T}) = T(2)
kurtosis{T<:Real}(d::Exponential{T}) = T(6)

entropy(d::Exponential) = 1 + log(d.θ)

#### Evaluation

zval(d::Exponential, x::Real) = x / d.θ
xval(d::Exponential, z::Real) = z * d.θ

pdf(d::Exponential, x::Real) = (λ = rate(d); x < 0 ? zero(λ) : λ * exp(-λ * x))
function logpdf{T<:Real}(d::Exponential{T}, x::Real)
    (λ = rate(d); x < 0 ? -T(Inf) : log(λ) - λ * x)
end

cdf{T<:Real}(d::Exponential{T}, x::Real) = x > 0 ? -expm1(-zval(d, x)) : zero(T)
ccdf{T<:Real}(d::Exponential{T}, x::Real) = x > 0 ? exp(-zval(d, x)) : zero(T)
logcdf{T<:Real}(d::Exponential{T}, x::Real) = x > 0 ? log1mexp(-zval(d, x)) : -T(Inf)
logccdf{T<:Real}(d::Exponential{T}, x::Real) = x > 0 ? -zval(d, x) : zero(T)

quantile(d::Exponential, p::Real) = -xval(d, log1p(-p))
cquantile(d::Exponential, p::Real) = -xval(d, log(p))
invlogcdf(d::Exponential, lp::Real) = -xval(d, log1mexp(lp))
invlogccdf(d::Exponential, lp::Real) = -xval(d, lp)

gradlogpdf{T<:Real}(d::Exponential{T}, x::Real) = x > 0 ? -rate(d) : zero(T)

mgf(d::Exponential, t::Real) = 1/(1 - t * scale(d))
cf(d::Exponential, t::Real) = 1/(1 - t * im * scale(d))

#### Sampling

rand(d::Exponential) = xval(d, randexp())

#### Fit model

immutable ExponentialStats <: SufficientStats
    sx::Float64   # (weighted) sum of x
    sw::Float64   # sum of sample weights

    ExponentialStats(sx::Real, sw::Real) = new(sx, sw)
end

suffstats{T<:Real}(::Type{Exponential}, x::AbstractArray{T}) = ExponentialStats(sum(x), length(x))
suffstats{T<:Real}(::Type{Exponential}, x::AbstractArray{T}, w::AbstractArray{Float64}) = ExponentialStats(dot(x, w), sum(w))

fit_mle(::Type{Exponential}, ss::ExponentialStats) = Exponential(ss.sx / ss.sw)


text:
# test for moveDownLevel
# Claudia April 2015

# starting topology: Case G
include("../examples/case_g_example.jl");
currT = deepcopy(net);
printEdges(currT)

# real network: Case H
df = readtable("CaseH_output.csv")
d = readTableCF(df)

optBL!(currT,d)

include("../examples/case_f_example.jl");
currT = deepcopy(net);
printEdges(currT)

# real network: Case H
df = readtable("CaseH_output.csv")
d = readTableCF(df)

optBL!(currT,d)


text:
include("DataDictionary.jl")
include("Kinetics.jl")
include("Balances.jl")

import PyPlot
using ODE
using DelimitedFiles

const plt = PyPlot


text:
# See estimate/posterior.jl for a test of the wrapper for dsgevecm_likelihood that
# takes a DSGEVECM as an input.
using DSGE, FileIO, Test, ModelConstructors
fp = dirname(@__FILE__)
@testset "Likelihood function of DSGE-VECM p(Y | θ, λ)" begin
    matdata = load(joinpath(fp, "../../reference/test_dsgevecm_likelihood.jld2"))

    # note nobs in the Matlab code is actually T in the Julia code (n_observations)
    # and lambdaT in the Matlab code is λ * T in the Julia code.
    @test matdata["lnpy"] ≈ DSGE.dsgevecm_likelihood(matdata["YYYY"],
                                                     matdata["XXYY"], matdata["XXXX"],
                                                     matdata["yyyyd"], matdata["xxyyd"],
                                                     matdata["xxxxd"], Int(matdata["nobs"]),
                                                     matdata["lambdaT"] / matdata["nobs"])
    @test matdata["lnpyinf"] ≈ DSGE.dsgevecm_likelihood(matdata["YYYY"],
                                                        matdata["XXYY"], matdata["XXXX"],
                                                        matdata["yyyyd"], matdata["xxyyd"],
                                                        matdata["xxxxd"], Int(matdata["nobs"]),
                                                        Inf)
    @test_throws DomainError matdata["lnpyinf"] ≈ DSGE.dsgevecm_likelihood(matdata["YYYY"],
                                                                           matdata["XXYY"], matdata["XXXX"],
                                                                           matdata["yyyyd"], matdata["xxyyd"],
                                                                           matdata["xxxxd"], Int(matdata["nobs"]),
                                                                           -1.)
end


text:
"""
    ThermInt([rng::AbstractRNG]; n_steps=30, n_samples=2000, n_warmup=500)
    ThermInt([rng::AbstractRNG], schedule; n_samples=2000, n_warmup=500)

- `schedule` can be any iterable object
- `n_steps` is the number of steps for the schedule using the formula
`(1:n_steps) ./ n_steps).^5`

A `ThermInt` object can then be used as a function:
```julia
alg = ThermInt(30)
alg(loglikelihood, logprior, x_init)
```
"""
struct ThermInt{AD,TRNG,V}
    schedule::V
    n_samples::Int
    n_warmup::Int
    rng::TRNG
end

function ThermInt(rng::AbstractRNG, schedule; n_samples::Int=2000, n_warmup::Int=500)
    return ThermInt{ADBACKEND[],typeof(rng),typeof(schedule)}(
        schedule, n_samples, n_warmup, rng
    )
end

function ThermInt(schedule; n_samples::Int=2000, n_warmup::Int=500)
    return ThermInt(GLOBAL_RNG, schedule; n_samples=n_samples, n_warmup=n_warmup)
end

function ThermInt(rng::AbstractRNG; n_steps::Int, n_samples::Int=2000, n_warmup::Int=500)
    return ThermInt(rng, range(0, 1, length=n_steps) .^ 5; n_samples=n_samples, n_warmup=n_warmup)
end

function ThermInt(; n_steps::Int=30, n_samples::Int=2000, n_warmup::Int=500)
    return ThermInt(
        GLOBAL_RNG, range(0, 1, length=n_steps) .^ 5; n_samples=n_samples, n_warmup=n_warmup
    )
end

struct TIParallelThreads end

function (alg::ThermInt)(
    loglikelihood, logprior, x_init::AbstractVector; progress=true, kwargs...
)
    p = ProgressMeter.Progress(length(alg.schedule); enabled=progress, desc="TI Sampling :")
    ΔlogZ = [
        begin
            ProgressMeter.next!(p)
            evaluate_loglikelihood(loglikelihood, logprior, alg, x_init, β; kwargs...)
        end for β in alg.schedule
    ]
    return trapz(alg.schedule, ΔlogZ)
end

function (alg::ThermInt)(
    loglikelihood,
    logprior,
    x_init::AbstractVector,
    ::TIParallelThreads;
    progress=true,
    kwargs...,
)
    Threads.nthreads() > 1 || @warn "Only one thread available, parallelization will not happen. Start Julia with `--threads n`"
    nsteps = length(alg.schedule)
    nthreads = min(Threads.nthreads(), nsteps)
    ΔlogZ = zeros(Float64, nsteps)
    algs = [deepcopy(alg) for _ in 1:nthreads]
    p = ProgressMeter.Progress(length(alg.schedule); enabled=progress, desc="TI Sampling :")
    Threads.@threads for i in 1:nsteps
        id = Threads.threadid()
        ΔlogZ[i] = evaluate_loglikelihood(
            loglikelihood, logprior, algs[id], x_init, alg.schedule[i]; kwargs...
        )
        ProgressMeter.next!(p)
    end
    return trapz(alg.schedule, ΔlogZ)
end

function evaluate_loglikelihood(loglikelihood, logprior, alg::ThermInt, x_init, β::Real)
    powerlogπ(θ) = β * loglikelihood(θ) + logprior(θ)
    samples = sample_powerlogπ(powerlogπ, alg, x_init)
    x_init .= samples[end] # Update the initial sample to be the last one of the chain
    return mean(loglikelihood, samples)
end

function sample_powerlogπ(powerlogπ, alg::ThermInt, x_init)
    D = length(x_init)
    metric = DiagEuclideanMetric(D)
    hamiltonian = get_hamiltonian(metric, powerlogπ, alg)

    initial_ϵ = find_good_stepsize(hamiltonian, x_init)
    integrator = Leapfrog(initial_ϵ)

    proposal = AdvancedHMC.NUTS{MultinomialTS,GeneralisedNoUTurn}(integrator)
    adaptor = StanHMCAdaptor(MassMatrixAdaptor(metric), StepSizeAdaptor(0.8, integrator))

    samples, stats = sample(
        alg.rng,
        hamiltonian,
        proposal,
        x_init,
        alg.n_samples,
        adaptor,
        alg.n_warmup;
        verbose=false,
        progress=false,
    )
    return samples
end

function get_hamiltonian(metric, powerlogπ, ::ThermInt{:ForwardDiff})
    return Hamiltonian(metric, powerlogπ, ForwardDiff)
end
function get_hamiltonian(metric, powerlogπ, ::ThermInt{:Zygote})
    return Hamiltonian(metric, powerlogπ, Zygote)
end
function get_hamiltonian(metric, powerlogπ, ::ThermInt{:ReverseDiff})
    return Hamiltonian(metric, powerlogπ, ReverseDiff)
end


text:
module GP

include("covariance.jl")
include("predict.jl")

end


text:
using SFML

function thread1()
    while true
        println("In thread1")
        sleep(1)
    end
end

function thread2()
    while true
        println("In thread 2")
        sleep(1)
    end
end

# const callback_c = Base.SingleAsyncWork(data -> thread1())
# async_send_pthread(func::Ptr{Void}) = (ccall(:uv_async_send, Cint, (Ptr{Void},), func); C_NULL)
# const c_async_send_pthread = cfunction(async_send_pthread, Ptr{Void}, (Ptr{Void},))

# t = Thread(c_async_send_pthread, callback_c.handle)
t = Thread(thread1)
launch(t)

# SFML.wait(t)

t2 = Thread(thread2)
launch(t2)

while true
    println("Main thread")
    sleep(1)
end


text:
##  COMPOSITE AND SURRUGOTE MODEL TYPES

# For example, we want to define

# abstract type ProbabilisticComposite <: Probabilistic end
# struct ProbabilisticSurrogate <: Probabilistic end
# Probabilistic() = ProbablisiticSurrogate()

# but also want this for all the abstract `Model` subtypes:

const COMPOSITE_TYPES = Symbol[]
const SURROGATE_TYPES = Symbol[]
const composite_types = Any[]
const surrogate_types = Any[]

for T in MLJModelInterface.ABSTRACT_MODEL_SUBTYPES
    composite_type_name = string(T, "Composite") |> Symbol
    surrogate_type_name = string(T, "Surrogate") |> Symbol

    @eval(abstract type $composite_type_name <: $T end)
    @eval(struct $surrogate_type_name <: $T end)

    push!(COMPOSITE_TYPES, composite_type_name)
    push!(SURROGATE_TYPES, surrogate_type_name)
    push!(composite_types, @eval($composite_type_name))
    push!(surrogate_types, @eval($surrogate_type_name))

    # shorthand surrogate constructor:
    @eval($T() = $surrogate_type_name())
end

const Surrogate = Union{surrogate_types...}
const Composite = Union{composite_types...}

MLJModelInterface.is_wrapper(::Type{<:Union{Composite,Surrogate}}) = true
MLJModelInterface.package_name(::Type{<:Union{Composite,Surrogate}}) = "MLJBase"
for T in surrogate_types
    MLJModelInterface.load_path(::Type{T}) = string("MLJBase.", T)
end

export MLJType, Model, Surrogate, Composite
for T in vcat(MMI.ABSTRACT_MODEL_SUBTYPES, COMPOSITE_TYPES, SURROGATE_TYPES)
    @eval(export $T)
end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export segmentedcontrol

"""
    segmentedcontrol(;kwargs...)

A SegmentedControl component.
Horizontal control made of multiple segments, alternative to RadioGroup. For more information, see: https://mantine.dev/core/segmented-control/
Keyword arguments:
- `id` (String; optional): The ID of this component, used to identify dash components in callbacks
- `className` (String; optional): Often used with CSS to style elements with common properties
- `color` (a value equal to: "dark", "gray", "red", "pink", "grape", "violet", "indigo", "blue", "cyan", "teal", "green", "lime", "yellow", "orange"; optional): Active control color from theme.colors, defaults to white in light color scheme and theme.colors.dark[9] in dark
- `data` (required): Data based on which controls are rendered. data has the following type: Array of lists containing elements 'label', 'value'.
Those elements have the following types:
  - `label` (String; required): The option's label
  - `value` (String; required): Option's values
- `fullWidth` (Bool; optional): True if component should have 100% width
- `radius` (a value equal to: "xs", "sm", "md", "lg", "xl" | Real; optional): Border-radius from theme or number to set border-radius in px
- `size` (a value equal to: "xs", "sm", "md", "lg", "xl"; optional): Controls font-size, paddings and height
- `style` (Dict; optional): Inline style override
- `value` (String; optional): Current selected value
"""
function segmentedcontrol(; kwargs...)
        available_props = Symbol[:id, :className, :color, :data, :fullWidth, :radius, :size, :style, :value]
        wild_props = Symbol[]
        return Component("segmentedcontrol", "SegmentedControl", "dash_mantine_components", available_props, wild_props; kwargs...)
end


text:
using libcg
using Test

@testset "libcg.jl" begin
    # Write your tests here.
end

max_stars_repo_name: simonbyrne/libcg

max_issues_repo_name: simonbyrne/libcg

max_forks_repo_name: simonbyrne/libcg


text:
"""
NEOSServer(;email)

Construct a `NEOSServer` object. The `email` argument should take the users valid
email address (required for solvers like CPLEX).
"""
mutable struct NEOSServer
    useragent::String
    host::String
    contenttype::String
    email::String
    NEOSServer(;email="") = new("NEOS.jl", "https://neos-server.org:3333", "text/xml", email)
end

function addemail!(s::NEOSServer, email::String)
    s.email = email
end

struct NEOSJob
    number::Int
    password::String
end

_add_text(value, arg::String) = add_text(value, arg)
_add_text(value, arg) = add_text(new_child(value, "int"), string(arg))

function buildxml(s::NEOSServer, name::String, args...)
    xml = XMLDocument()
    mthd = create_root(xml, "methodCall")
    mname = new_child(mthd, "methodName")
    add_text(mname, name)
    if length(args) > 0
	params = new_child(mthd, "params")
	for a in args
	    param = new_child(params, "param")
	    value = new_child(param, "value")
	    _add_text(value, a)
	end
    end
    return string(xml)
end

function apimethod(s::NEOSServer, name::String, args...)
    xml = buildxml(s, name, args...)
    res = send(s, xml)
    parseresponse(res)
end

function send(s::NEOSServer, xml::String)
    hdrs = Dict{String, String}("user-agent" => s.useragent, "host" => s.host,
	                        "content-type" => s.contenttype, "content-length" => string(length(xml)))
    res = HTTP.request("POST", s.host; headers=hdrs, body=xml)
end

function getvalues!(values, node)
    types = ["int", "i4", "string", "double", "base64", "dateTime.iso8601"]
    for child in child_nodes(node)
	if name(child) in types
	    push!(values, content(child))
	else
	    getvalues!(values, child)
	end
    end
end

function parseresponse(res)
    if res.status != 200
	error("XML-RPC failed with code: $(res.status)")
    end
    parameters = Any[]
    xml = parse_string(String(res.body))

    xroot = root(xml)
    getvalues!(parameters, xroot)
    return parameters
end

# ===========================================
#
# 	NEOS API Methods
#
function neosHelp(s::NEOSServer)
    return apimethod(s, "help")[1]
end

messages = [:emailHelp, :welcome, :version, :ping, :printQueue]
for m in messages
    m_str = string(m)
    @eval function ($m)(s::NEOSServer)
	return apimethod(s, $(m_str))[1]
    end
end

lists = [:listAllSolvers, :listCategories]
for m in lists
    m_str = string(m)
    @eval 	function ($m)(s::NEOSServer)
	return apimethod(s, $(m_str))
    end
end

function getSolverTemplate(s::NEOSServer, category::String, solvername::String, inputMethod::String)
    apimethod(s, "getSolverTemplate", category, solvername, inputMethod)[1]
end

function listSolversInCategory(s::NEOSServer, category::String)
    apimethod(s, "listSolversInCategory", category)
end

function submitJob(s::NEOSServer, xmlstring::String, print_level::Int = 1)
    res = apimethod(s, "submitJob", xmlstring)
    if print_level >= 1
	println("==================")
	println("NEOS Job submitted")
	println("number:\t$(res[1])")
	println("pwd:\t$(res[2])")
	println("==================")
    end
    return NEOSJob(parse(Int, res[1]), res[2])
end

jobapimethods = [:getJobStatus, :killJob]
for m in jobapimethods
    m_str = string(m)
    @eval 	function ($m)(s::NEOSServer, j::NEOSJob)
	return apimethod(s, $(m_str), j.number, j.password)[1]
    end
end

function getJobInfo(s::NEOSServer, j::NEOSJob)
    apimethod(s, "getJobInfo", j.number, j.password)
end

for s in ["", "NonBlocking"]
    _final_str = "getFinalResults$s"
    _final = Symbol(_final_str)
    _intermediate_str = "getIntermediateResults$s"
    _intermediate = Symbol(_intermediate_str)
    @eval function ($_final)(s::NEOSServer, j::NEOSJob)
	decode_to_string(apimethod(s, $(_final_str), j.number, j.password)[1])
    end

    @eval function ($(_intermediate))(s::NEOSServer, j::NEOSJob; offset=0)
	(results, offset) = apimethod(s, $(_intermediate_str), j.number, j.password, offset)
	decode_to_string(results), parse(Int, offset)
    end
end

function decode_to_string(s)
    String(base64decode(replace(s, "\n" => "")))
end

struct NEOSData
    input::String
    category::String
    model::String
    data::String
    commands::String
    comments::String
end

function createNEOSJob(input::String, category::String, model::String, data::String, commands::String, comments::String)
    return NEOSData(input, category, model, data, commands, comments)
end

function getNEOSJobXML(solver::String, m :: NEOSData)
    category = m.category
    input = m.input
    model = m.model
    data = m.data
    commands = m.commands
    comments = m.comments

    return "<document>\n<category>$category</category>\n<solver>$solver</solver>\n<inputMethod>$input</inputMethod>\n\n<model><![CDATA[\n$model\n]]></model>\n\n<data><![CDATA[\n$data\n]]></data>\n\n<commands><![CDATA[\n$commands\n]]></commands>\n\n<comments><![CDATA[\n$comments\n]]></comments>\n\n</document>\n"
end


text:
import LowRankModels: ConvergenceHistory, update_ch!
export ConvergenceHistory, update_ch!

# type ConvergenceHistory
#     name::AbstractString
#     objective::Array
#     dual_objective::Array
#     primal_residual::Array
#     dual_residual::Array
#     times::Array
#     stepsizes::Array
#     optval
# end
# ConvergenceHistory(name::AbstractString,optval=0) = ConvergenceHistory(name,Float64[],Float64[],Float64[],Float64[],Float64[],Float64[],optval)

# function update!(ch::ConvergenceHistory, dt; obj=0, dual_obj=0, stepsize=0, pr=0, dr=0)
#     push!(ch.objective,obj)
#     push!(ch.dual_objective,dual_obj)
#     push!(ch.primal_residual,pr)
#     push!(ch.dual_residual,dr)
#     push!(ch.stepsizes,stepsize)
#     if isempty(ch.times)
#         push!(ch.times,dt)
#     else
#         push!(ch.times,ch.times[end]+dt)
#     end
# end

# function show(ch::ConvergenceHistory)
#     print("Convergence History for $(ch.name)\n\n")
#     @printf "%16s%16s\n" "time (s)" "objective"
#     for i=1:length(ch.objective)
#         @printf "%16.2e%16.4e\n" ch.times[i] ch.objective[i]
#     end
# end

max_stars_repo_name: madeleineudell/FirstOrderOptimization


text:
# TODO:
# - parse for discrete valves and add labels for those links (valves as part of pipes are
#   done)
# x parse for demand at nodes and add demand labels
# x add fill color for each node that scales with elevation
# x add colorbar for elevation colors (how to do it?)
# x add an indication of edge length
# x add coordinates for the nodes
# - allow user to choose an output style that uses coordinates

"""
Build networkx graph object from a WaterModels network dictionary parsed from an EPANET file.
"""
function build_graph(data::Dict{String,Any})
    # presumes input is data dict from WaterModels.parse_file()

    # TODO:
    # - check that input dictionary has the expected fields

    G = pgv.AGraph(strict=false, directed=true)

    add_nodes!(G, data["node"])
    demand!(G, data["junction"])

    add_links!(G, data["pipe"])
    add_links!(G, data["pump"])
    add_links!(G, data["valve"])

    return G
end # funtion build_graph

function add_nodes!(G::PyCall.PyObject, nodes::Dict{String,Any})
    # This needs to be redone. Nodes are now disctinct from junctions, tanks, and reservoirs. JJS 7/22/20

    # determine max and min elevations and coordinates
    elmin = 1e16
    elmax = 1e-16
    xmin = 1e16
    xmax = 1e-16
    ymin = 1e16
    ymax = 1e-16
    for (key,node) in nodes
        elmin = min(elmin, node["elevation"])
        elmax = max(elmax, node["elevation"])
        if haskey(node, "coordinates")
            coord = node["coordinates"]
            xmin = min(xmin, coord[1])
            xmax = max(xmax, coord[1])
            ymin = min(ymin, coord[2])
            ymax = max(ymax, coord[2])
        end
    end
    xspan = xmax-xmin
    yspan = ymax-ymin
    # scale to use to get "good" position results with graphviz output -- should
    # calculate this systematically from the number of nodes?
    scale = 20

    # add elmin and elmax as graph attributes for use in creating a colorbar in
    # colorbar (or via write_visualization)
    PyCall.set!(G."graph_attr", "elmin", elmin)
    PyCall.set!(G."graph_attr", "elmax", elmax)

    # create and populate networkx nodes
    for (key,node) in nodes
        node_type = node["source_id"][1]
        if node_type == "reservoir"
            name = "Rsvr\n"*node["name"]
        elseif node_type == "tank"
            name = "Tank\n"*node["name"]
        else
            name = node["name"]
        end

        # color by elevation
        elev = node["elevation"]
        clr = HSV(get(viridis, elev, (elmin, elmax)))
        # convert clr to a string that is in a form suitable for graphviz
        clrstr = string(clr.h/360)*" "*string(clr.s)*" "*string(clr.v)

        # change font color depending on the background color
        if clr.v < 0.6
            fntclr = "white"
        else
            fntclr = "black"
        end

        if haskey(node, "coordinates")
            coord = node["coordinates"]
            xval = scale*(coord[1] - xmin)/xspan
            yval = scale*(coord[2] - ymin)/yspan
            posstr = string(xval)*","*string(yval)*"!" # exclamation point forces exact
                                                     # positioning when using Neato
            G.add_node(node["index"], label=name, elevation=elev, pos=posstr,
                       style="filled", fillcolor=clrstr, fontcolor=fntclr)
        else
            G.add_node(node["index"], label=name, elevation=elev, style="filled",
                       fillcolor=clrstr, fontcolor=fntclr)
        end
    end
end

""" Add demand label for junctions."""
function demand!(G::PyCall.PyObject, junctions::Dict{String,Any})
    for (key,junction) in junctions
        #### are the keys in data["node"] alwyas equal to string(index) ??? will
        #### presume so for now, JJS 7/7/20
        ## I don't think this bit is needed with pygraphviz, JJS 7/17/20
        # if typeof(key) == String
        #     index = parse(Int, key)
        # else
        #     index = key
        # end

        nodeobj = @pycall G."get_node"(key)::PyObject
        name = get(nodeobj.attr, "label")
        dem = @sprintf("%2.2g", junction["demand"])
        PyCall.set!(nodeobj.attr, "label", name*"\nd = "*dem)
    end
end

function add_links!(G::PyCall.PyObject, links::Dict{String,Any})
    for (key,link) in links
        link_type = link["source_id"][1]
        if link_type == "pump"
            G.add_edge(link["node_fr"], link["node_to"], link["index"],
                       label="P "*link["name"], color="red", style="bold")
        elseif link_type == "valve"
            println("valves not yet implemented")
        else # it is a pipe
            length = @sprintf("%2.2g", link["length"])
            type = link["status"]
            if type != "Open"
                G.add_edge(link["node_fr"], link["node_to"], link["index"],
                           label=type*"\n"*length)
            else
                G.add_edge(link["node_fr"], link["node_to"], link["index"], label=length)
            end
        end
    end
end

"""
Write out to a file a visualization for a WaterModels network dictionary parsed from an
EPANET file. `basefilename` should not include an extension and will be appended with
`_w_cb.pdf` in the final output file, which is a multi-page PDF. Use `del_files=false` to
keep the intermediate files.
"""
function write_visualization(data::Dict{String,Any}, basefilename::String;
                             del_files::Bool=true)
    # TODO:
    # - pass through arguments to `write_graph`

    #gvfile = basefilename*".gv"
    pdffile = basefilename*".pdf"
    cbfile = basefilename*"_cbar.pdf"
    outfile = basefilename*"_w_cb.pdf"

    G = build_graph(data)
    write_graph(G, pdffile)
    colorbar(G, cbfile)

    # add option -dAutoRotatePages=/None ?
    run(`gs -sDEVICE=pdfwrite -dNOPAUSE -dQUIET -dBATCH -sOutputFile=$outfile $cbfile $pdffile`)
    # delete the intermediate files
    if del_files
        run(`rm $pdffile $cbfile`)
    end
end

"""
Use graphviz (via pygraphviz) to output a visualization to a file for a graph
"""
function write_graph(G::PyCall.PyObject, filename::String)
    # TODO:
    # - allow different graphviz programs to be used, e.g., neato
    # - allow other arguments to be passed through to graphviz
    G.draw(filename, prog="dot")
end

"""
create and save a colorbar that represents the node elevations
"""
function colorbar(G::PyCall.PyObject, filename::String)
    elmin = parse(Float64, get(G.graph_attr, "elmin"))
    elmax = parse(Float64, get(G.graph_attr, "elmax"))
    elmid = elmin + (elmax-elmin)/2

    x = reshape(collect(range(0.0, stop=1.0, length=100)), (1,:))
    Plots.heatmap(x, c=:viridis, size=(500,100), legend=:none, yaxis=false)
    Plots.plot!(xticks=(0:50:100, [elmin, elmid, elmax]))
    Plots.title!("Elevation")
    Plots.savefig(filename)
end


text:
module SIMDIndexing

using SIMD
using Requires

const VecOfIntegers{N, T <: SIMD.IntegerTypes} = Vec{N, T}

abstract type AbstractIndex end
const Index = Union{Int, VecOfIntegers, VecRange, AbstractIndex}

function __init__()
    @require(MappedArrays="dbb5928d-eab1-5f90-85c2-b9b0edb7c900",
             include("mappedarrays.jl"))
    @require(FillArrays="1a297f60-69ca-5386-bcde-b61e274b549b",
             include("fillarrays.jl"))
end

end # module


text:
export residual, residual!, jac_residual, jac_structure_residual, jac_structure_residual!
export jac_coord_residual!, jac_coord_residual, jprod_residual, jprod_residual!
export jtprod_residual, jtprod_residual!, jac_op_residual, jac_op_residual!
export hess_residual, hess_structure_residual, hess_structure_residual!
export hess_coord_residual!, hess_coord_residual, jth_hess_residual
export hprod_residual, hprod_residual!, hess_op_residual, hess_op_residual!
export project!, project, get_meta

"""
    Fx = residual(model, x)
Computes ``F(x)``, the residual at x.
"""
function NLPModels.residual(model::AbstractVIModel{T, S}, x::AbstractVector{T}) where {T, S}
  @lencheck model.meta.nvar x
  Fx = S(undef, model.meta.nvar)
  residual!(model, x, Fx)
end

"""
    Fx = residual!(model, x, Fx)
Computes ``F(x)``, the residual at x.
"""
function residual! end

"""
    Jx = jac_residual(model, x)
Computes ``J(x)``, the Jacobian of the residual at x.
"""
function NLPModels.jac_residual(model::AbstractVIModel, x::AbstractVector)
  @lencheck model.meta.nvar x
  rows, cols = jac_structure_residual(model)
  vals = jac_coord_residual(model, x)
  sparse(rows, cols, vals, model.meta.nvar, model.meta.nvar)
end

"""
    (rows,cols) = jac_structure_residual!(model, rows, cols)
Returns the structure of the constraint's Jacobian in sparse coordinate format in place.
"""
function jac_structure_residual! end

"""
    (rows,cols) = jac_structure_residual(model)
Returns the structure of the constraint's Jacobian in sparse coordinate format.
"""
function NLPModels.jac_structure_residual(model::AbstractVIModel)
  rows = Vector{Int}(undef, model.meta.nnzj)
  cols = Vector{Int}(undef, model.meta.nnzj)
  jac_structure_residual!(model, rows, cols)
end

"""
    vals = jac_coord_residual!(model, x, vals)
Computes the Jacobian of the residual at `x` in sparse coordinate format, rewriting
`vals`. `rows` and `cols` are not rewritten.
"""
function jac_coord_residual! end

"""
    (rows,cols,vals) = jac_coord_residual(model, x)
Computes the Jacobian of the residual at `x` in sparse coordinate format.
"""
function NLPModels.jac_coord_residual(model::AbstractVIModel, x::AbstractVector)
  @lencheck model.meta.nvar x
  vals = Vector{eltype(x)}(undef, model.meta.nnzj)
  jac_coord_residual!(model, x, vals)
end

"""
    Jv = jprod_residual(model, x, v)
Computes the product of the Jacobian of the residual at x and a vector, i.e.,  ``J(x)v``.
"""
function NLPModels.jprod_residual(
  model::AbstractVIModel{T, S},
  x::AbstractVector{T},
  v::AbstractVector,
) where {T, S}
  @lencheck model.meta.nvar x v
  Jv = S(undef, model.meta.nvar)
  jprod_residual!(model, x, v, Jv)
end

"""
    Jv = jprod_residual!(model, x, v, Jv)
Computes the product of the Jacobian of the residual at x and a vector, i.e.,  ``J(x)v``, storing it in `Jv`.
"""
function jprod_residual! end

"""
    Jv = jprod_residual!(model, rows, cols, vals, v, Jv)
Computes the product of the Jacobian of the residual given by `(rows, cols, vals)`
and a vector, i.e.,  ``J(x)v``, storing it in `Jv`.
"""
function NLPModels.jprod_residual!(
  model::AbstractVIModel,
  rows::AbstractVector{<:Integer},
  cols::AbstractVector{<:Integer},
  vals::AbstractVector,
  v::AbstractVector,
  Jv::AbstractVector,
)
  @lencheck model.meta.nnzj rows cols vals
  @lencheck model.meta.nvar v Jv
  increment!(model, :neval_jprod_residual)
  coo_prod!(rows, cols, vals, v, Jv)
end

"""
    Jtv = jtprod_residual(model, x, v)
Computes the product of the transpose of the Jacobian of the residual at x and a vector, i.e.,  ``J(x)^Tv``.
"""
function NLPModels.jtprod_residual(
  model::AbstractVIModel{T, S},
  x::AbstractVector{T},
  v::AbstractVector,
) where {T, S}
  @lencheck model.meta.nvar x v
  Jtv = S(undef, model.meta.nvar)
  jtprod_residual!(model, x, v, Jtv)
end

"""
    Jtv = jtprod_residual!(model, x, v, Jtv)
Computes the product of the transpose of the Jacobian of the residual at x and a vector, i.e.,  ``J(x)^Tv``, storing it in `Jtv`.
"""
function jtprod_residual! end

"""
    Jtv = jtprod_residual!(model, rows, cols, vals, v, Jtv)
Computes the product of the transpose of the Jacobian of the residual given by `(rows, cols, vals)`
and a vector, i.e.,  ``J(x)^Tv``, storing it in `Jv`.
"""
function NLPModels.jtprod_residual!(
  model::AbstractVIModel,
  rows::AbstractVector{<:Integer},
  cols::AbstractVector{<:Integer},
  vals::AbstractVector,
  v::AbstractVector,
  Jtv::AbstractVector,
)
  @lencheck model.meta.nnzj rows cols vals
  @lencheck model.meta.nvar v Jtv
  increment!(model, :neval_jtprod_residual)
  coo_prod!(cols, rows, vals, v, Jtv)
end

"""
    Jx = jac_op_residual(model, x)
Computes ``J(x)``, the Jacobian of the residual at x, in linear operator form.
"""
function NLPModels.jac_op_residual(model::AbstractVIModel{T, S}, x::AbstractVector{T}) where {T, S}
  @lencheck model.meta.nvar x
  Jv = S(undef, model.meta.nvar)
  Jtv = S(undef, model.meta.nvar)
  return jac_op_residual!(model, x, Jv, Jtv)
end

"""
    Jx = jac_op_residual!(model, rows, cols, vals, Jv, Jtv)
Computes ``J(x)``, the Jacobian of the residual given by `(rows, cols, vals)`, in linear operator form. The
vectors `Jv` and `Jtv` are used as preallocated storage for the operations.
"""
function NLPModels.jac_op_residual!(
  model::AbstractVIModel,
  rows::AbstractVector{<:Integer},
  cols::AbstractVector{<:Integer},
  vals::AbstractVector,
  Jv::AbstractVector,
  Jtv::AbstractVector,
)
  @lencheck model.meta.nnzj rows cols vals
  @lencheck model.meta.nvar Jv Jtv
  prod! = @closure (res, v, α, β) -> begin
    jprod_residual!(model, rows, cols, vals, v, Jv)
    if β == 0
      @. res = α * Jv
    else
      @. res = α * Jv + β * res
    end
    return res
  end
  ctprod! = @closure (res, v, α, β) -> begin
    jtprod_residual!(model, rows, cols, vals, v, Jtv)
    if β == 0
      @. res = α * Jtv
    else
      @. res = α * Jtv + β * res
    end
    return res
  end
  return LinearOperator{eltype(vals)}(
    model.meta.nvar,
    model.meta.nvar,
    false,
    false,
    prod!,
    ctprod!,
    ctprod!,
  )
end

"""
    H = hess_residual(model, x, v)
Computes the linear combination of the Hessians of the residuals at `x` with coefficients
`v`.
A `Symmetric` object wrapping the lower triangle is returned.
"""
function NLPModels.hess_residual(model::AbstractVIModel, x::AbstractVector, v::AbstractVector)
  @lencheck model.meta.nvar x v
  rows, cols = hess_structure_residual(model)
  vals = hess_coord_residual(model, x, v)
  Symmetric(sparse(rows, cols, vals, model.meta.nvar, model.meta.nvar), :L)
end

"""
    (rows,cols) = hess_structure_residual(model)
Returns the structure of the residual Hessian.
"""
function NLPModels.hess_structure_residual(model::AbstractVIModel)
  rows = Vector{Int}(undef, model.meta.nnzh)
  cols = Vector{Int}(undef, model.meta.nnzh)
  hess_structure_residual!(model, rows, cols)
end

"""
    hess_structure_residual!(model, rows, cols)
Returns the structure of the residual Hessian in place.
"""
function hess_structure_residual! end

"""
    vals = hess_coord_residual!(model, x, v, vals)
Computes the linear combination of the Hessians of the residuals at `x` with coefficients
`v` in sparse coordinate format, rewriting `vals`.
"""
function hess_coord_residual! end

"""
    vals = hess_coord_residual(model, x, v)
Computes the linear combination of the Hessians of the residuals at `x` with coefficients
`v` in sparse coordinate format.
"""
function NLPModels.hess_coord_residual(model::AbstractVIModel, x::AbstractVector, v::AbstractVector)
  @lencheck model.meta.nvar x v
  vals = Vector{eltype(x)}(undef, model.meta.nnzh)
  hess_coord_residual!(model, x, v, vals)
end

"""
    Hj = jth_hess_residual(model, x, j)
Computes the Hessian of the j-th residual at x.
"""
function NLPModels.jth_hess_residual(model::AbstractVIModel, x::AbstractVector, j::Int)
  @lencheck model.meta.nvar x
  increment!(model, :neval_jhess_residual)
  decrement!(model, :neval_hess_residual)
  v = [i == j ? one(eltype(x)) : zero(eltype(x)) for i = 1:(model.meta.nvar)]
  return hess_residual(model, x, v)
end

"""
    Hiv = hprod_residual(model, x, i, v)
Computes the product of the Hessian of the i-th residual at x, times the vector v.
"""
function NLPModels.hprod_residual(
  model::AbstractVIModel{T, S},
  x::AbstractVector{T},
  i::Int,
  v::AbstractVector,
) where {T, S}
  @lencheck model.meta.nvar x
  Hv = S(undef, model.meta.nvar)
  hprod_residual!(model, x, i, v, Hv)
end

"""
    Hiv = hprod_residual!(model, x, i, v, Hiv)
Computes the product of the Hessian of the i-th residual at x, times the vector v, and stores it in vector Hiv.
"""
function hprod_residual! end

"""
    Hop = hess_op_residual(model, x, i)
Computes the Hessian of the i-th residual at x, in linear operator form.
"""
function NLPModels.hess_op_residual(
  model::AbstractVIModel{T, S},
  x::AbstractVector{T},
  i::Int,
) where {T, S}
  @lencheck model.meta.nvar x
  Hiv = S(undef, model.meta.nvar)
  return hess_op_residual!(model, x, i, Hiv)
end

get_meta(model::AbstractVIModel) = model.meta

"""
    project!(Px::AbstractVector{T}, model::VIModel, x::AbstractVector{T}) where {T}

Compute the projection of d over X in-place, i.e.,
```math
    min_x 0.5 | d - x |²₂ s.t. x ∈ X
```
"""
function project! end

"""
    project(model::VIModel, x::AbstractVector{T}) where {T}

Compute the projection of d over X, i.e.,
```math
    min_x 0.5 | d - x |²₂ s.t. x ∈ X
```
```
"""
function project(model::AbstractVIModel{T, S}, x::AbstractVector{T}) where {T, S}
  @lencheck model.meta.nvar x
  Px = S(undef, model.meta.nvar)
  project!(model, x, Px)
end


text:
#=
This page defines the functions:
fwd_prop! - for feeding forward through the network
bwd_prop! - for feeding backward through the network to calculate gradients

Different methods are defined for different layer types and the network as
a whole.
=#

################################################################################
#
# Forward Propagate
#
################################################################################

########################################
# Layer -> OutputLayer
########################################

# LinearOutputLayer does not have an activation function
# Write weighted sum directly to A
function fwd_prop!(ℓ::LinOutputLayer, prev_ℓ::HDL)

    # A⁽²⁾ = b⁽²⁾ .+ W⁽²⁾ A⁽¹⁾
    A_mul_B!(ℓ.A, ℓ.W, prev_ℓ.Ar)
    broadcast!(+, ℓ.A, ℓ.b, ℓ.A)

    # No activation
end

function fwd_prop!(ℓ::LinOutputLayer, prev_ℓ::LDL)

    # Z⁽²⁾ = b⁽²⁾ .+ W⁽²⁾ A⁽¹⁾
    A_mul_B!(ℓ.A, ℓ.W, prev_ℓ.A)
    broadcast!(+, ℓ.A, ℓ.b, ℓ.A)

    # No activation
end

# LogisticOutputLayer
# fwd_prop! method covered by the general LD method defined below

# Softmax Layer activation function different to standard
# Vector `u` used for the max values of each row
function fwd_prop!(ℓ::SoftmaxOutputLayer, prev_ℓ::HDL)

    # Z⁽²⁾ = b⁽²⁾ .+ W⁽²⁾ A⁽¹⁾
    A_mul_B!(ℓ.Z, ℓ.W, prev_ℓ.Ar)
    broadcast!(+, ℓ.Z, ℓ.b, ℓ.Z)

    # A⁽²⁾ = f(Z⁽²⁾)
    ℓ.f!(ℓ.A, ℓ.Z, ℓ.u)
end

function fwd_prop!(ℓ::SoftmaxOutputLayer, prev_ℓ::LDL)

    # Z⁽²⁾ = b⁽²⁾ .+ W⁽²⁾ A⁽¹⁾
    A_mul_B!(ℓ.Z, ℓ.W, prev_ℓ.A)
    broadcast!(+, ℓ.Z, ℓ.b, ℓ.Z)

    # A⁽²⁾ = f(Z⁽²⁾)
    ℓ.f!(ℓ.A, ℓ.Z, ℓ.u)
end

########################################
# Layer -> Low Dimension Layer
########################################

# LD -> LD
function fwd_prop!(ℓ::LDL, prev_ℓ::LDL)

    # Z⁽²⁾ = b⁽²⁾ .+ W⁽²⁾ A⁽¹⁾
    A_mul_B!(ℓ.Z, ℓ.W, prev_ℓ.A)
    broadcast!(+, ℓ.Z, ℓ.b, ℓ.Z)

    # A⁽²⁾ = f(Z⁽²⁾)
    ℓ.f!(ℓ.A, ℓ.Z)
end

# HD -> LD
function fwd_prop!(ℓ::LDL, prev_ℓ::HDL)

    # Z⁽²⁾ = b⁽²⁾ .+ W⁽²⁾ A⁽¹⁾
    A_mul_B!(ℓ.Z, ℓ.W, prev_ℓ.Ar)
    broadcast!(+, ℓ.Z, ℓ.b, ℓ.Z)

    # A⁽²⁾ = f(Z⁽²⁾)
    ℓ.f!(ℓ.A, ℓ.Z)
end

########################################
# Layer -> Word Embedding Layer
########################################

# LD -> Word Embedding
function fwd_prop!(ℓ::WordEmbeddingLayer, prev_ℓ::LDL)

    # Embed word features
    embed!(ℓ.A, ℓ.W, prev_ℓ.A)
end

# HD -> Word Embedding
function fwd_prop!(ℓ::WordEmbeddingLayer, prev_ℓ::HDL)

    # Embed word features
    embed!(ℓ.A, ℓ.W, prev_ℓ.Ar)
end

########################################
# Layer -> Sparse Encoder Layer
########################################

# LD -> Sparse Encoder
function fwd_prop!(ℓ::SparseEncoderLayer, prev_ℓ::LDL)

    # Z⁽²⁾ = b⁽²⁾ .+ W⁽²⁾ A⁽¹⁾
    A_mul_B!(ℓ.Z, ℓ.W, prev_ℓ.A)
    broadcast!(+, ℓ.Z, ℓ.b, ℓ.Z)

    # A⁽²⁾ = f(Z⁽²⁾)
    ℓ.f!(ℓ.A, ℓ.Z)

    # Average of activations - "probability" of each neuron firing
    mean!(ℓ.p̂, ℓ.A)
end

# HD -> Sparse Encoder
function fwd_prop!(ℓ::SparseEncoderLayer, prev_ℓ::HDL)

    # Z⁽²⁾ = b⁽²⁾ .+ W⁽²⁾ A⁽¹⁾
    A_mul_B!(ℓ.Z, ℓ.W, prev_ℓ.A)
    broadcast!(+, ℓ.Z, ℓ.b, ℓ.Z)

    # A⁽²⁾ = f(Z⁽²⁾)
    ℓ.f!(ℓ.A, ℓ.Z)

    # Average of activations - "probability" of each neuron firing
    mean!(ℓ.p̂, ℓ.A)
end

########################################
# High Dimension Layer -> High Dimension Layer
########################################

# Convolution
function fwd_prop!(ℓ::ConvolutionLayer, prev_ℓ::HDL)

    # Z⁽²⁾ = b⁽²⁾ .+ W⁽²⁾ ∗ A⁽¹⁾
    conv!(ℓ.Z, ℓ.W, prev_ℓ.A)
    bcastbias!(ℓ.Z, ℓ.b)

    # A⁽²⁾ = f(Z⁽²⁾)
    ℓ.f!(ℓ.A, ℓ.Z)
end

# Mean Pooling
function fwd_prop!(ℓ::MeanPoolLayer, prev_ℓ::HDL)
    meanpool!(ℓ.A, prev_ℓ.A)
end

# Max Pooling
function fwd_prop!(ℓ::MaxPoolLayer, prev_ℓ::HDL)
    maxpool!(ℓ.A, prev_ℓ.A, ℓ.maxinds)
end

########################################
# Whole Network
########################################

function fwd_prop!(net::NeuralNet)

    # Feed forward each hidden layer and the output layer
    for (ℓ, prev_ℓ) in zip(net.layers[2:end], net.layers[1:end-1])
        fwd_prop!(ℓ, prev_ℓ)
    end
end

################################################################################
#
# Backward Propagate
#
################################################################################

########################################
# Layer -> OutputLayer
########################################

function bwd_prop!(ℓ::OutputLayer, prev_ℓ::HDL)

    # Delta of final layer
    # Δ⁽ᴸ⁾ = 1/m (A⁽ᴸ⁾ - Y)               (recall: A⁽ᴸ⁾ = H)
    ℓ.δ!(ℓ.Δ, ℓ.A, ℓ.y, ℓ.m)

    # Gradient of bias vector
    # ∇b⁽ᴸ⁾ = sum(Δ⁽ᴸ⁾, 2)
    sum!(ℓ.∇b, ℓ.Δ)

    # Gradient of weight matrix
    # ∇W⁽ᴸ⁾ = Δ⁽ᴸ⁾ A⁽ᴸ⁻¹⁾ᵀ
    A_mul_Bt!(ℓ.∇W, ℓ.Δ, prev_ℓ.Ar)
end

function bwd_prop!(ℓ::OutputLayer, prev_ℓ::LDL)

    # Delta of final layer
    # Δ⁽ᴸ⁾ = 1/m (A⁽ᴸ⁾ - Y)               (recall: A⁽ᴸ⁾ = H)
    ℓ.δ!(ℓ.Δ, ℓ.A, ℓ.y, ℓ.m)

    # Gradient of bias vector
    # ∇b⁽ᴸ⁾ = sum(Δ⁽ᴸ⁾, 2)
    sum!(ℓ.∇b, ℓ.Δ)

    # Gradient of weight matrix
    # ∇W⁽ᴸ⁾ = Δ⁽ᴸ⁾ A⁽ᴸ⁻¹⁾ᵀ
    A_mul_Bt!(ℓ.∇W, ℓ.Δ, prev_ℓ.A)
end

########################################
# Low -> Low -> Low
########################################

# Low -> Low -> Low
function bwd_prop!(ℓ::LDL, prev_ℓ::LDL, next_ℓ::LDL)

    # Δ⁽²⁾ = W⁽³⁾ᵀ Δ⁽³⁾ ⊙ f'(Z⁽²⁾)
    # Broken into 2 steps:
    # (a) - Δ⁽²⁾ := W⁽³⁾ᵀ Δ⁽³⁾
    # (b) - Δ⁽²⁾ := Δ⁽²⁾ ⊙ f'(Z⁽²⁾)
    At_mul_B!(ℓ.Δ, next_ℓ.W, next_ℓ.Δ)
    ℓ.∇f!(ℓ.Δ, ℓ.Z)

    # Gradient of bias vector
    # ∇b⁽²⁾ = sum(Δ⁽²⁾, 2)
    sum!(ℓ.∇b, ℓ.Δ)

    # Gradient of weight matrix
    # ∇W⁽²⁾ = Δ⁽²⁾ A⁽¹⁾
    A_mul_Bt!(ℓ.∇W, ℓ.Δ, prev_ℓ.A)
end

# Low -> Word embedding -> Low
function bwd_prop!(ℓ::WordEmbeddingLayer, prev_ℓ::LDL, next_ℓ::LDL)

    # Δ⁽²⁾ = W⁽³⁾ᵀ Δ⁽³⁾
    At_mul_B!(ℓ.Δ, next_ℓ.W, next_ℓ.Δ)

    # Gradient of weight matrix
    embed_grad!(ℓ.∇W, ℓ.Δ, prev_ℓ.A)
end

# Low -> Sparse -> Low
function bwd_prop!(ℓ::SparseEncoderLayer, prev_ℓ::LDL, next_ℓ::LDL)

    # Δ⁽²⁾ = (W⁽³⁾ᵀ Δ⁽³⁾ .+ spgrad⁽²⁾) ⊙ f'(Z⁽²⁾)
    # Broken into 3 steps:
    # (a) - Δ⁽²⁾ := W⁽³⁾ᵀ Δ⁽³⁾
    # (b) - Δ⁽²⁾ := Δ⁽²⁾ .+ spgrad⁽²⁾
    # (c) - Δ⁽²⁾ := Δ⁽²⁾ ⊙ f'(Z⁽²⁾)
    At_mul_B!(ℓ.Δ, next_ℓ.W, next_ℓ.Δ)
    ℓ.δ!(ℓ.Δ, ℓ.spgrad, ℓ.p̂, ℓ.p, ℓ.β, ℓ.m)
    ℓ.∇f!(ℓ.Δ, ℓ.Z)

    # Gradient of bias vector
    # ∇b⁽²⁾ = sum(Δ⁽²⁾, 2)
    sum!(ℓ.∇b, ℓ.Δ)

    # Gradient of weight matrix
    # ∇W⁽²⁾ = Δ⁽²⁾ A⁽¹⁾
    A_mul_Bt!(ℓ.∇W, ℓ.Δ, prev_ℓ.A)
end

########################################
# High -> Low -> Low
########################################

# High -> Low -> Low
function bwd_prop!(ℓ::LDL, prev_ℓ::HDL, next_ℓ::LDL)

    # Δ⁽²⁾ = W⁽³⁾ᵀ Δ⁽³⁾ ⊙ f'(Z⁽²⁾)
    # Broken into 2 steps:
    # (a) - Δ⁽²⁾ := W⁽³⁾ᵀ Δ⁽³⁾
    # (b) - Δ⁽²⁾ := Δ⁽²⁾ ⊙ f'(Z⁽²⁾)
    At_mul_B!(ℓ.Δ, next_ℓ.W, next_ℓ.Δ)
    ℓ.∇f!(ℓ.Δ, ℓ.Z)

    # Gradient of bias vector
    # ∇b⁽²⁾ = sum(Δ⁽²⁾, 2)
    sum!(ℓ.∇b, ℓ.Δ)

    # Gradient of weight matrix
    # ∇W⁽²⁾ = Δ⁽²⁾ A⁽¹⁾
    A_mul_Bt!(ℓ.∇W, ℓ.Δ, prev_ℓ.Ar)
end

# High -> Word embedding -> Low
function bwd_prop!(ℓ::WordEmbeddingLayer, prev_ℓ::HDL, next_ℓ::LDL)

    # Δ⁽²⁾ = W⁽³⁾ᵀ Δ⁽³⁾
    At_mul_B!(ℓ.Δ, next_ℓ.W, next_ℓ.Δ)

    # Gradient of weight matrix
    embed_grad!(ℓ.∇W, ℓ.Δ, prev_ℓ.Ar)
end

# High -> Sparse -> Low
function bwd_prop!(ℓ::SparseEncoderLayer, prev_ℓ::HDL, next_ℓ::LDL)

    # Δ⁽²⁾ = (W⁽³⁾ᵀ Δ⁽³⁾ .+ spgrad⁽²⁾) ⊙ f'(Z⁽²⁾)
    # Broken into 3 steps:
    # (a) - Δ⁽²⁾ := W⁽³⁾ᵀ Δ⁽³⁾
    # (b) - Δ⁽²⁾ := Δ⁽²⁾ .+ spgrad⁽²⁾
    # (c) - Δ⁽²⁾ := Δ⁽²⁾ ⊙ f'(Z⁽²⁾)
    At_mul_B!(ℓ.Δ, next_ℓ.W, next_ℓ.Δ)
    ℓ.δ!(ℓ.Δ, ℓ.spgrad, ℓ.p̂, ℓ.p, ℓ.β, ℓ.m)
    ℓ.∇f!(ℓ.Δ, ℓ.Z)

    # Gradient of bias vector
    # ∇b⁽²⁾ = sum(Δ⁽²⁾, 2)
    sum!(ℓ.∇b, ℓ.Δ)

    # Gradient of weight matrix
    # ∇W⁽²⁾ = Δ⁽²⁾ A⁽¹⁾
    A_mul_Bt!(ℓ.∇W, ℓ.Δ, prev_ℓ.Ar)
end

########################################
# High -> High -> Low
########################################

# High -> Convolution -> Low
function bwd_prop!(ℓ::ConvolutionLayer, prev_ℓ::HDL, next_ℓ::LDL)

    # Δ⁽²⁾ = W⁽³⁾ᵀ Δ⁽³⁾ ⊙ f'(Z⁽²⁾)
    # Broken into 2 steps:
    # (a) - Δ⁽²⁾ := W⁽³⁾ᵀ Δ⁽³⁾
    # (b) - Δ⁽²⁾ := Δ⁽²⁾ ⊙ f'(Z⁽²⁾)
    At_mul_B!(ℓ.Δr, next_ℓ.W, next_ℓ.Δ)
    ℓ.∇f!(ℓ.Δr, ℓ.Zr)

    # Gradient of bias vector
    # ∇bₖ⁽²⁾ = ∑ Δₖ⁽²⁾
    biasgrad!(ℓ.∇b, ℓ.Δ)

    # Gradient of weight matrices
    # ∇Wₖ⁽²⁾ = Δₖ⁽²⁾ ∗ A⁽¹⁾
    convgrad!(ℓ.∇W, ℓ.Δ, prev_ℓ.A)

    return
end

# High -> Pool -> Low
function bwd_prop!(ℓ::PoolLayer, prev_ℓ::HDL, next_ℓ::LDL)

    # Δ⁽²⁾ = W⁽³⁾ᵀ Δ⁽³⁾
    At_mul_B!(ℓ.Δr, next_ℓ.W, next_ℓ.Δ)
end

########################################
# High -> High -> High
########################################

# High -> Convolution -> MeanPool
function bwd_prop!(ℓ::ConvolutionLayer, prev_ℓ::HDL, next_ℓ::MeanPoolLayer)

    # Δ⁽²⁾ = upsample( Δ⁽³⁾ ) ⊙ f'(Z⁽²⁾)
    # Broken into 2 steps:
    # (a) - Δ⁽²⁾ := upsample( Δ⁽³⁾ )
    # (b) - Δ⁽²⁾ := Δ⁽²⁾ ⊙ f'(Z⁽²⁾)
    upsample_mean!(ℓ.Δ, next_ℓ.Δ)
    ℓ.∇f!(ℓ.Δr, ℓ.Zr)

    # Gradient of bias vector
    # ∇bₖ⁽²⁾ = ∑ Δₖ⁽²⁾
    biasgrad!(ℓ.∇b, ℓ.Δ)

    # Gradient of weight matrices
    # ∇Wₖ⁽²⁾ = Δₖ⁽²⁾ ∗ A⁽¹⁾
    convgrad!(ℓ.∇W, ℓ.Δ, prev_ℓ.A)
end

# High -> Convolution -> MaxPool
function bwd_prop!(ℓ::ConvolutionLayer, prev_ℓ::HDL, next_ℓ::MaxPoolLayer)

    # Δ⁽²⁾ = upsample( Δ⁽³⁾ ) ⊙ f'(Z⁽²⁾)
    # Broken into 2 steps:
    # (a) - Δ⁽²⁾ := upsample( Δ⁽³⁾ )
    # (b) - Δ⁽²⁾ := Δ⁽²⁾ ⊙ f'(Z⁽²⁾)
    upsample_max!(ℓ.Δ, next_ℓ.Δ, next_ℓ.maxinds)
    ℓ.∇f!(ℓ.Δr, ℓ.Zr)

    # Gradient of bias vector
    # ∇bₖ⁽²⁾ = ∑ Δₖ⁽²⁾
    biasgrad!(ℓ.∇b, ℓ.Δ)

    # Gradient of weight matrices
    # ∇Wₖ⁽²⁾ = Δₖ⁽²⁾ ∗ A⁽¹⁾
    convgrad!(ℓ.∇W, ℓ.Δ, prev_ℓ.A)
end

# High -> Pool -> High
function bwd_prop!(ℓ::PoolLayer, prev_ℓ::HDL, next_ℓ::ConvolutionLayer)

    # Delta from convolution layer
    # Δₖ⁽²⁾ = W⁽³⁾ ∗ Δₖ⁽³⁾
    conv!(ℓ.Δ, next_ℓ.W, next_ℓ.Δ)
end
# TODO - CHECK THIS !!!

########################################
# Whole Network
########################################

function bwd_prop!(net::NeuralNet)

    # Final layer
    ℓ = net.layers[end]
    prev_ℓ = net.layers[end-1]
    bwd_prop!(ℓ, prev_ℓ)

    # Hidden layers
    for (prev_ℓ, ℓ, next_ℓ) in zip(net.layers[end-2:-1:1], net.layers[end-1:-1:2], net.layers[end:-1:3])
        bwd_prop!(ℓ, prev_ℓ, next_ℓ)
    end

    # Regularisation
    net.reg_grad!(net)
end


text:
type ModelSaver
    file
    bestscore
end

ModelSaver(file) = ModelSaver(file, 0.0)

function (s::ModelSaver)(model, score)
    if s.bestscore < score
        println("\nSAVER: BEST SCORE $score")
        print("SAVER: SAVING MODEL TO $(s.file)... ")
        open(io -> serialize(io, model), s.file, "w")
        s.bestscore = score
        println("DONE")
    end
end

max_stars_repo_name: hshindo/JukaiNLP

max_issues_repo_name: hshindo/JukaiNLP

max_forks_repo_name: hshindo/JukaiNLP


text:
module Templates

using CrypticCrosswords: DerivedSolution, explain
using HTTP: HTTP

sanitize(s::AbstractString) = replace(s, r"[<>]" => "")
sanitize(s) = sanitize(string(s)::AbstractString)

function explain_to_string(derivation::DerivedSolution)
    io = IOBuffer()
    explain(io, derivation)
    String(take!(io))
end

function Answer(s::DerivedSolution)
    """

    $(sanitize(s.output)) $(round(Int, s.similarity * 100))%

 """
end

function Explanation(s::DerivedSolution)
    sanitize(explain_to_string(s))
end

function Solution(s::DerivedSolution)
    """

    $(Answer(s))
    $(Explanation(s))

"""
end

function Results(derived_solutions)
    """

Results

    $(join(Solution(s) for s in derived_solutions))

    """
end

function ClueInput(clue="", length="", pattern="")
    if length === nothing
        length = ""
    end
    """
    <form action="solve" method="GET">

                Clue
                <input id="clue" type="text" name="clue" value="$(sanitize(clue))" style="width: 40em; max-width: 100%"/>

                Length (optional)
                <input id="length" type="number" name="length" value="$(string(length))" style="width: 9em"/>

                Regex (optional)
                <input id="pattern" type="text" name="pattern" value="$(sanitize(pattern))" style="width: 10em"/>

            <input type="submit" value="Solve" style="margin-top: 1em"/>

    </form>
    """
end

function HomeLink()
    """

CrypticCrosswords.jl

    """
end

function Intro()
    """

    This is a general tool for solving cryptic (or "British-style") crossword clues, written entirely in the Julia programming language. You can find the source code for the solver on Github at rdeits/CrypticCrosswords.jl.

    """
end

function Example(clue, length, pattern, answer)
    """

        <form action="solve" method="GET">
            <input type="submit" id="try" value="Try"/>
            <input type="text" name="clue" value="$(sanitize(clue))" hidden />
            <input type="number" name="length" value="$(sanitize(length))" hidden />
            <input type="text" name="pattern" value="$(sanitize(pattern))" hidden />
            $(clue) ($(length)) $pattern → $(uppercase(answer))
        </form>

    """
end

function Examples()
    EXAMPLES = [
        ("Couch is unfinished until now", 4, "", "sofa"),
        ("Spin broken shingle", 7, "", "english"),
        ("Initially babies are naked", 4, "", "bare"),
        ("At first, congoers like us eschew solving hints", 5, "", "clues"),
        ("Initial meetings disappoint rosemary internally", 6, "", "intros"),
        ("M's Rob Titon pitching slider?", 10, "", "trombonist"),
        ("Aerial worker Anne on the way up", 7, "", "antenna"),
        ("In glee over unusual color", 10, "^o", "olive green")
    ]
    """

Examples

            $(join(Example(c...) for c in EXAMPLES))

    """
end

function Style()
    """

    """
end

function Index(body)
    """
    <!DOCTYPE html>

        <meta charset="utf-8" />
        <link rel="shortcut icon" href="favicon.ico" />
        <meta name="viewport" content="width=device-width, initial-scale=1" />
        <meta name="theme-color" content="#000000" />
        <title>CrypticCrosswords.jl</title>
        $(Style())

                $(body)

    """
end

end


text:
struct RepetitionBlueprint
    origin::Int
    destination::Int
    size::Int
end

struct HeterozygosityBlueprint
    position::Int
    alleles::Vector{DNA}
end

position(h::HeterozygosityBlueprint) = h.position
alleles(h::HeterozygosityBlueprint) = h.alleles
Base.isless(a::HeterozygosityBlueprint, b::HeterozygosityBlueprint) = a.position < b.position

struct ChromosomeBlueprint
    seqlen::Int
    available::Set{UnitRange{Int}}
    ncopies::Int
    hets::Vector{HeterozygosityBlueprint}
    intra_repetitions::Vector{RepetitionBlueprint}
end

function ChromosomeBlueprint(len::Int, ncopies::Int)
    return ChromosomeBlueprint(len, Set([1:len]), ncopies, Vector{HeterozygosityBlueprint}(), Vector{RepetitionBlueprint}())
end

function ChromosomeBlueprint(cp::ChromosomeBlueprint, v::Vector{RepetitionBlueprint})
    return ChromosomeBlueprint(cp.seqlen, cp.available, cp.ncopies, cp.hets, v)
end

function ChromosomeBlueprint(cp::ChromosomeBlueprint, v::Vector{HeterozygosityBlueprint})
    return ChromosomeBlueprint(cp.seqlen, cp.available, cp.ncopies, v, cp.intra_repetitions)
end

function ChromosomeBlueprint(cp::ChromosomeBlueprint, s::Set{UnitRange{Int}})
    return ChromosomeBlueprint(cp.seqlen, s, cp.ncopies, cp.hets, cp.intra_repetitions)
end

"Get the number of copies of a chromosome planned in a chromosome blueprint."
ncopies(cp::ChromosomeBlueprint) = cp.ncopies

"Get a copy of the available blocks of sequence in a chromosome blueprint."
available_blocks(cp::ChromosomeBlueprint) = copy(cp.available)

"Get a reference to the available blocks of sequence in a chromosome blueprint."
available_blocks_unsafe(cp::ChromosomeBlueprint) = cp.available

"Get a copy of the repetitions planned in a chromosome blueprint."
repetitions(cp::ChromosomeBlueprint) = copy(cp.intra_repetitions)

"Get a refernece to the repetitions planned in a chromosome blueprint."
repetitions_unsafe(cp::ChromosomeBlueprint) = cp.intra_repetitions

"Get a copy of the heterozygosity planned in a chromosome blueprint."
hets(cp::ChromosomeBlueprint) = copy(cp.hets)

"Get a reference to the heterozygosity planned in a chromosome blueprint."
hets_unsafe(cp::ChromosomeBlueprint) = cp.hets

"Get the length (in bp) of the chromosome planned in a chromosome blueprint."
Base.length(cp::ChromosomeBlueprint) = cp.seqlen

function Base.show(io::IO, cp::ChromosomeBlueprint)
    println(io, "Specification for a chromosome:")
    println(io, " Number of copies: ", ncopies(cp))
    println(io, " Length of chromosomes: ", length(cp))
    ablocks = sort!(collect(available_blocks_unsafe(cp)))
    asum = sum(length(x) for x in ablocks)
    println(io, ' ', asum, "bp available in the following stretches:")
    if length(ablocks) <= 20
        for b in ablocks
            println("  ", first(b), "bp .. ", last(b), "bp")
        end
    else
        for i in 1:10
            b = ablocks[i]
            println("  ", first(b), "bp .. ", last(b), "bp")
        end
        println("    ⋮")
        li = lastindex(ablocks)
        for i in (li - 9):li
            b = ablocks[i]
            println("  ", first(b), "bp .. ", last(b), "bp")
        end
    end
    hlen = length(hets_unsafe(cp))
    if hlen > 0
        println(' ', length(hets_unsafe(cp)), " heterozygous positions at:")
        hs = sort!(hets(cp))
        if hlen <= 20
            for h in hs
                println("  location: ", position(h), "bp  bases: ", string(alleles(h)))
            end
        else
            for i in 1:10
                h = hs[i]
                println("  location: ", position(h), "bp  bases: ", string(alleles(h)))
            end
            println("    ⋮")
            li = lastindex(hs)
            for i in (li - 9):li
                h = hs[i]
                println("  location: ", position(h), "bp  bases: ", string(alleles(h)))
            end
        end
    end
end

function count_het(cb::ChromosomeBlueprint)
    counts = zeros(Int, ncopies(cb), ncopies(cb))
    for het in hets_unsafe(cb)
        for i in 1:ncopies(cb)
            for j in i:ncopies(cb)
                inc = ifelse(het.alleles[i] != het.alleles[j], 1, 0)
                counts[i, j] = counts[j, i] = (counts[i, j] + inc)
            end
        end
    end
    return counts
end

function Base.copy(cp::ChromosomeBlueprint)
    a = available_blocks(cp)
    s = copy(cp.snps)
    r = repetitions(cp)
    return ChromosomeBlueprint(cp.seqlen, a, cp.ncopies, s, r)
end

# Choosing and consuming available regions of the genome

## Internals

function find_big_block(blocks::Set{UnitRange{Int}}, region::UnitRange{Int})
    for block in blocks
        if contains(block, region)
            return block
        end
    end
    throw(ArgumentError(string("Region ", region, " is not contained by an available block.")))
end

function consume_region!(blocks::Set{UnitRange{Int}}, region::UnitRange{Int})
    bigblock = find_big_block(blocks, region)
    pop!(blocks, bigblock)
    l, r = complimentary(bigblock, region)
    if !isempty(l)
        push!(blocks, l)
    end
    if !isempty(r)
        push!(blocks, r)
    end
    return blocks
end

## Public API

"""
    suggest_regions(cp::ChromosomeBlueprint, size::Int, n::Int)

A useful utility function to assist planning features in a chromosome blueprint.

This method returns a vector of non-overlapping, regions of the chromosome planned
represented by the ChromosomeBlueprint `cb`.

These regions are free regions: they are untouched by any other planned features,
and so may be used when planning other features (See also: [`plan_repetition`](@ref), [`plan_het`](@ref)).

The regions will be `size`bp in length.

!!! warning
    This function was designed for use interactively in a julia session.

    If this method cannot find `n` free regions of the size you've asked for, it
    will still give an output vector containing the regions it did manage to find,
    but it will issue a warning to the terminal. This will get increasingly
    likely as you fill the chromosome blueprint up with features.

    If you are using this method in a script or program where you depend on a
    reliable number of regions, either add a check to make sure you got the
    number of regions you need, or use this method interactively, and hard code
    an appropriate output into your script.
"""
function suggest_regions(cp::ChromosomeBlueprint, size::Int, n::Int)
    suitable_blocks = filter(x -> length(x) ≥ size, available_blocks_unsafe(cp))
    suggested_regions = Vector{UnitRange{Int}}(undef, n)
    nfound = 0
    while nfound < n
        if isempty(suitable_blocks) # We've run out of big enough available blocks...
            @warn string("Could only suggest ", nfound, " available regions of ", size, "bp or longer")
            break
        end
        bigblock = rand(suitable_blocks)
        pop!(suitable_blocks, bigblock)
        smallblock = pick_available_subregion(bigblock, size)
        l, r = complimentary(bigblock, smallblock)
        if length(l) ≥ size
            push!(suitable_blocks, l)
        end
        if length(r) ≥ size
            push!(suitable_blocks, r)
        end
        suggested_regions[nfound += 1] = smallblock
    end
    return resize!(suggested_regions, nfound)
end

function consume_regions(cp::ChromosomeBlueprint, blocks::Vararg{UnitRange{Int}})
    avb = available_blocks(cp)
    for b in blocks
        consume_region!(avb, b)
    end
    return ChromosomeBlueprint(cp, avb)
end

"Returns true if the `region` of the chromosome in the blueprint is unoccupied by a planned feature."
function isavailable(cp::ChromosomeBlueprint, region::UnitRange{Int})
    for block in available_blocks_unsafe(cp)
        contains(block, region) && return true
    end
    return false
end

function checkavailability(cp::ChromosomeBlueprint, from::UnitRange{Int}, to::UnitRange{Int})
    if !isavailable(cp, from)
        throw(ArgumentError("The first region is not available"))
    end
    if !isavailable(cp, to)
        throw(ArgumentError("The second region is not available"))
    end
end

function checkavailability(cp::ChromosomeBlueprint, interval::UnitRange{Int})
    if !isavailable(cp, interval)
        throw(ArgumentError("requested region not available: $interval"))
    end
end


text:
ROOT = pwd()
include(joinpath(ROOT,"..","src_PowSysMod", "PowSysMod_body.jl"))

"""Test
Test to launch in complex-modeler/src_test folder for example D:\repo\complex-modeler\src_test>julia global_test.jl

# Arguments
- type of instance (GOCInput, MatpowerInput, IIDMInput)
if GOCInput : julia global_test.jl GOCInput <instances folder> <instance name> <scenario name>
if MatpowerInput or IIDMInput : julia global_test.jl MatpowerInput <instances folder> <instance name>

"""
function read_arguments(ARGS)
    if length(ARGS)==0
        error("Arguments must be :
        Argument 1 is the type of instance (GOCInput, MatpowerInput or IIDMInput)\n
         julia global_test.jl GOCInput <instances folder> <instance name> <scenario name>
         julia global_test.jl MatpowerInput <instances folder> <instance name>
         example : julia global_test.jl MatpowerInput ..\data_Matpower\matpower case9.m
         julia global_test.jl IIDMInput <instances folder> <instance name>

        ")
    end
    inputtype = ARGS[1]
    if inputtype=="GOCInput"
        inputtype = GOCInput
        if length(ARGS)!=5
            error("Arguments must be : julia global_test.jl GOCInput <instances folder> <instance name> <scenario name> <otuput folder>")
        else
            instance_path = joinpath(pwd(), ARGS[2], ARGS[3],ARGS[4])
            output_path = ARGS[3]
        end

    elseif inputtype=="MatpowerInput"
        inputtype = MatpowerInput
        if length(ARGS)!=4
            error("Arguments must be : julia global_test.jl MatpowerInput <instances folder> <instance name> <otuput folder>")
        else
            instance_path = joinpath(pwd(), ARGS[2], ARGS[3])
            output_path = ARGS[4]
        end

    elseif inputtype=="IIDMInput"
        inputtype = IIDMInput
        if length(ARGS)!=4
            error("Arguments must be : julia global_test.jl IIDMInput <instances folder> <instance name> <otuput folder>")
        else
            instance_path = joinpath(pwd(), ARGS[2], ARGS[3])
            output_path = ARGS[4]
        end
    else
        error("Unexpected input type.")
    end
    return inputtype, instance_path, output_path
end

function main(ARGS)
    inputtype, instance_path, output_path = read_arguments(ARGS)

    ## Building network description
    OPFpbs = load_OPFproblems(inputtype, instance_path)

    ## Setting Generators labels for coupling constraints if GOC instance
    inputtype == GOCInput && introduce_Sgenvariables!(OPFpbs)

    ## Building polynomila optim problem
    pb_global = build_globalpb!(OPFpbs)

    ## Exporting
    export_to_dat(pb_global, output_path)
end

main(ARGS)

max_stars_repo_name: Gimaju/GOC

max_issues_repo_name: Gimaju/GOC

max_forks_repo_name: Gimaju/GOC