text:
import ..isdir_nothrow, ..Registry.RegistrySpec, ..isurl

###############
# PackageSpec #
###############
"""
Parser for PackageSpec objects.
"""
function parse_package(args::Vector{QString}, options; add_or_dev=false)::Vector{PackageSpec}
    args = PackageToken[PackageToken(pkgword) for pkgword in package_lex(args)]
    return parse_package_args(args; add_or_dev=add_or_dev)
end

struct VersionToken
    version::String
end

struct Rev
    rev::String
end

struct Subdir
    dir::String
end

const PackageIdentifier = String
const PackageToken = Union{PackageIdentifier, VersionToken, Rev, Subdir}

    # Match a git repository URL. This includes uses of `@` and `:` but
    # requires that it has `.git` at the end.
let url = raw"((git|ssh|http(s)?)|(git@[\w\-\.]+))(:(//)?)([\w\.@\:/\-~]+)(\.git$)(/)?",

    # Match a `NAME=UUID` package specifier.
    name_uuid = raw"[^@\#\s:]+\s*=\s*[^@\#\s:]+",

    # Match a `#BRANCH` branch or tag specifier.
    branch = raw"\#\s*[^@\#\s]*",

    # Match an `@VERSION` version specifier.
    version = raw"@\s*[^@\#\s]*",

    # Match a `:SUBDIR` subdir specifier.
    subdir = raw":[^@\#\s]+",

    # Match any other way to specify a package. This includes package
    # names, local paths, and URLs that don't match the `url` part. In
    # order not to clash with the branch, version, and subdir
    # specifiers, these cannot include `@` or `#`, and `:` is only
    # allowed if followed by `/` or `\`. For URLs matching this part
    # of the regex, that means that `@` (e.g. user names) and `:`
    # (e.g. port) cannot be used but it doesn't have to end with
    # `.git`.
    other = raw"([^@\#\s:] | :(/|\\))+"

    # Combine all of the above.
    global const package_id_re = Regex(
        "$url | $name_uuid | $branch | $version | $subdir | $other", "x")
end

function package_lex(qwords::Vector{QString})::Vector{String}
    words = String[]
    for qword in qwords
        qword.isquoted ?
            push!(words, qword.raw) :
            append!(words, map(m->m.match, eachmatch(package_id_re, qword.raw)))
    end
    return words
end

PackageToken(word::String)::PackageToken =
    first(word) == '@' ? VersionToken(word[2:end]) :
    first(word) == '#' ? Rev(word[2:end]) :
    first(word) == ':' ? Subdir(word[2:end]) :
    String(word)

function parse_package_args(args::Vector{PackageToken}; add_or_dev=false)::Vector{PackageSpec}
    # check for and apply PackageSpec modifier (e.g. `#foo` or `@v1.0.2`)
    function apply_modifier!(pkg::PackageSpec, args::Vector{PackageToken})
        (isempty(args) || args[1] isa PackageIdentifier) && return
        modifier = popfirst!(args)
        if modifier isa Subdir
            pkg.subdir = modifier.dir
            (isempty(args) || args[1] isa PackageIdentifier) && return
            modifier = popfirst!(args)
        end

        if modifier isa VersionToken
            pkg.version = modifier.version
        elseif modifier isa Rev
            pkg.rev = modifier.rev
        else
            pkgerror("Package name/uuid must precede subdir specifier `[$arg]`.")
        end
    end

    pkgs = PackageSpec[]
    while !isempty(args)
        arg = popfirst!(args)
        if arg isa PackageIdentifier
            pkg = parse_package_identifier(arg; add_or_develop=add_or_dev)
            apply_modifier!(pkg, args)
            push!(pkgs, pkg)
        # Modifiers without a corresponding package identifier -- this is a user error
        else
            arg isa VersionToken ?
                pkgerror("Package name/uuid must precede version specifier `@$arg`.") :
            arg isa Rev ?
                pkgerror("Package name/uuid must precede revision specifier `#$(arg.rev)`.") :
                pkgerror("Package name/uuid must precede subdir specifier `[$arg]`.")
        end
    end
    return pkgs
end

let uuid = raw"(?i)[0-9a-z]{8}-[0-9a-z]{4}-[0-9a-z]{4}-[0-9a-z]{4}-[0-9a-z]{12}(?-i)",
    name = raw"(\w+)(?:\.jl)?"
    global const name_re = Regex("^$name\$")
    global const uuid_re = Regex("^$uuid\$")
    global const name_uuid_re = Regex("^$name\\s*=\\s*($uuid)\$")
end
# packages can be identified through: uuid, name, or name+uuid
# additionally valid for add/develop are: local path, url
function parse_package_identifier(word::AbstractString; add_or_develop=false)::PackageSpec
    if add_or_develop
        if isurl(word)
            return PackageSpec(; url=word)
        elseif any(occursin.(['\\','/'], word)) || word == "." || word == ".."
            if casesensitive_isdir(expanduser(word))
                return PackageSpec(; path=normpath(expanduser(word)))
            else
                pkgerror("`$word` appears to be a local path, but directory does not exist")
            end
        end
        if occursin(name_re, word) && casesensitive_isdir(expanduser(word))
            @info "Use `./$word` to add or develop the local directory at `$(Base.contractuser(abspath(word)))`."
        end
    end
    if occursin(uuid_re, word)
        return PackageSpec(;uuid=UUID(word))
    elseif occursin(name_re, word)
        return PackageSpec(String(match(name_re, word).captures[1]))
    elseif occursin(name_uuid_re, word)
        m = match(name_uuid_re, word)
        return PackageSpec(String(m.captures[1]), UUID(m.captures[2]))
    else
        pkgerror("Unable to parse `$word` as a package.")
    end
end

################
# RegistrySpec #
################
function parse_registry(raw_args::Vector{QString}, options; add=false)
    regs = RegistrySpec[]
    foreach(x -> push!(regs, parse_registry(x; add=add)), unwrap(raw_args))
    return regs
end

# Registries can be identified through: uuid, name, or name+uuid
# when updating/removing. When adding we can accept a local path or url.
function parse_registry(word::AbstractString; add=false)::RegistrySpec
    word = expanduser(word)
    registry = RegistrySpec()
    if add && isdir_nothrow(word) # TODO: Should be casesensitive_isdir
        if isdir(joinpath(word, ".git")) # add path as url and clone it from there
            registry.url = abspath(word)
        else # put the path
            registry.path = abspath(word)
        end
    elseif occursin(uuid_re, word)
        registry.uuid = UUID(word)
    elseif occursin(name_re, word)
        registry.name = String(match(name_re, word).captures[1])
    elseif occursin(name_uuid_re, word)
        m = match(name_uuid_re, word)
        registry.name = String(m.captures[1])
        registry.uuid = UUID(m.captures[2])
    elseif add
        # Guess it is a url then
        registry.url = String(word)
    else
        pkgerror("`$word` cannot be parsed as a registry")
    end
    return registry
end

#
# # Other
#
function parse_activate(args::Vector{QString}, options)
    isempty(args) && return [] # nothing to do
    if length(args) == 1
        x = first(args)
        if x.isquoted
            return [x.raw]
        end
        x = x.raw
        if first(x) == '@'
            options[:shared] = true
            return [x[2:end]]
        else
            return [expanduser(x)]
        end
    end
    return args # this is currently invalid input for "activate"
end

#
# # Option Maps
#
function do_preserve(x::String)
    x == "all"    && return Types.PRESERVE_ALL
    x == "direct" && return Types.PRESERVE_DIRECT
    x == "semver" && return Types.PRESERVE_SEMVER
    x == "none"   && return Types.PRESERVE_NONE
    x == "tiered" && return Types.PRESERVE_TIERED
    pkgerror("`$x` is not a valid argument for `--preserve`.")
end


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons, Josef Kircher
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

###############
# Prepare FMU #
###############

import FMIImport.FMICore: fmi3FMUState

myFMU = fmi3Load("BouncingBall", "ModelicaReferenceFMUs", "0.0.14")
inst = fmi3InstantiateCoSimulation!(myFMU; loggingOn=true)
@test inst != 0

@test fmi3EnterInitializationMode(inst) == 0
@test fmi3ExitInitializationMode(inst) == 0

###########################
# Testing state functions #
###########################

if fmi3CanGetSetState(myFMU) && fmi3CanSerializeFMUState(myFMU)
    @test fmi3GetFloat64(inst, "h") == 1.0
    FMUState = fmi3GetFMUState(inst)
    @test typeof(FMUState) == fmi3FMUState
    len = fmi3SerializedFMUStateSize(inst, FMUState)
    @test len > 0
    serial = fmi3SerializeFMUState(inst, FMUState)
    @test length(serial) == len
    @test typeof(serial) == Array{Cuchar,1}

    fmi3SetFloat64(inst, "h", 10.0)
    FMUState = fmi3GetFMUState(inst)
    @test fmi3GetFloat64(inst, "h") == 10.0

    FMUState2 = fmi3DeSerializeFMUState(inst, serial)
    @test typeof(FMUState2) == fmi3FMUState
    fmi3SetFMUState(inst, FMUState2)
    @test fmi3GetFloat64(inst, "h") == 1.0
    fmi3SetFMUState(inst, FMUState)
    @test fmi3GetFloat64(inst, "h") == 10.0
    fmi3FreeFMUState!(inst, FMUState)
    fmi3FreeFMUState!(inst, FMUState2)
else
    @info "The FMU provided from the tool `$(ENV["EXPORTINGTOOL"])` does not support state get, set, serialization and deserialization. Skipping related tests."
end

############
# Clean up #
############

@test fmi3Terminate(inst) == 0
fmi3Unload(myFMU)

ext: jl


text:
"""
    solve_rachford_rice(K, z, [V]; <keyword arguments>)

Compute vapor mole fraction `V` for given equilibrium constants `K` and mole fractions `z`.

# Arguments
`K` - Equal length to `z`, containing the equilibrium constants for each component.
`z` - Mole fractions. Should sum up to unity.

# Keyword arguments
`tol = 1e-12` - Tolerance for solve.
`maxiter` - Maximum number of iterations
`ad` - Use automatic differentiation (ForwardDiff) instead of analytical gradient.

# Examples
```julia-repl
julia> solve_rachford_rice([0.5, 1.5], [0.3, 0.7])
0.8000000000000002
```
"""
function solve_rachford_rice(K, z, V = NaN; tol=1e-12, maxiter=1000, ad = false)
    # Lowest bound for solution
    V_min = 1/(1 - maximum(K))
    # Largest bound for solution
    V_max = 1/(1 - minimum(K))
    if isnan(V)
        V = (V_min + V_max)/2
    end
    i = 1
    n = length(z)
    while i < maxiter
        if ad
            f(V) = objectiveRR(V, K, z)
            r = f(V)
            Jo = ForwardDiff.derivative(f, V)
            V = V - Jo\r
        else
            denum = 0.0
            r = 0.0
            @inbounds for i = 1:n
                ΔK = K[i]-1
                a = z[i]*ΔK
                b = (1 + V*ΔK)
                r += a/b
                denum += z[i]*ΔK^2/b^2
            end
            V = V + r/denum
        end
        if abs(r) < tol
            break
        end
        if V < V_min || V > V_max
            obj = (v) -> objectiveRR(v, K, z)
            V = find_zero(obj, (V_min, V_max), Bisection())
            break
        end
        i += 1
    end
    return V
end

function objectiveRR(V, K, z)
    eq = 0.0
    for (i, k) in enumerate(K)
        @inbounds eq = eq + ((k - 1.0)*z[i])/(1.0 + V*(k - 1.0))
    end
    return eq
end


text: (b"[encoded data, 65536 characters]", (1024, 4))

max_stars_repo_name: joshua-maros/julia-blue-noise

max_issues_repo_name: joshua-maros/julia-blue-noise

max_forks_repo_name: joshua-maros/julia-blue-noise


text:
module DGHeatEquationTests

using Gridap
using ForwardDiff
using LinearAlgebra
using Test
using GridapODEs.ODETools
using GridapODEs.TransientFETools
using Gridap.FESpaces: get_algebraic_operator

θ = 0.2

u(x,t) = (1.0-x[1])*x[1]*(1.0-x[2])*x[2]*t
u(t::Real) = x -> u(x,t)
f(t) = x -> ∂t(u)(x,t)-Δ(u(t))(x)

L= 1.0
n = 2
domain = (0,L,0,L)
partition = (n,n)
model = CartesianDiscreteModel(domain,partition)

order = 2

reffe = ReferenceFE(lagrangian,Float64,order)
V0 = FESpace(
  model,
  reffe,
  conformity=:L2
)
U = TransientTrialFESpace(V0)

Ω = Triangulation(model)
degree = 2*order
dΩ = Measure(Ω,degree)

Γ = BoundaryTriangulation(model)
dΓ = Measure(Γ,degree)
nb = get_normal_vector(Γ)

Λ = SkeletonTriangulation(model)
dΛ = Measure(Λ,degree)
ns = get_normal_vector(Λ)

a(u,v) = ∫(∇(v)⋅∇(u))dΩ
b(v,t) = ∫(v*f(t))dΩ
m(u,v) = ∫(v*u)dΩ

h = 1.0 / n
γ = order*(order+1)
a_Γ(u,v) = ∫( (γ/h)*v*u - v*(∇(u)⋅nb) - (∇(v)⋅nb)*u )dΓ
b_Γ(v,t) = ∫( (γ/h)*v*u(t) - (∇(v)⋅nb)*u(t) )dΓ

a_Λ(u,v) = ∫( (γ/h)*jump(v*ns)⊙jump(u*ns) - jump(v*ns)⊙mean(∇(u)) - mean(∇(v))⊙jump(u*ns) )dΛ

res(t,u,v) = a(u,v) + m(∂t(u),v) + a_Γ(u,v) + a_Λ(u,v) - b(v,t) - b_Γ(v,t)
jac(t,u,du,v) = a(du,v) + a_Γ(du,v) + a_Λ(du,v)
jac_t(t,u,dut,v) = m(dut,v)

op = TransientFEOperator(res,jac,jac_t,U,V0)

t0 = 0.0
tF = 1.0
dt = 0.1

U0 = U(0.0)
uh0 = interpolate_everywhere(u(0.0),U0)

ls = LUSolver()
using Gridap.Algebra: NewtonRaphsonSolver
ode_solver = ThetaMethod(ls,dt,θ)

sol_t = solve(ode_solver,op,uh0,t0,tF)

l2(w) = w*w

tol = 1.0e-6
_t_n = t0

for (uh_tn, tn) in sol_t
  global _t_n
  _t_n += dt
  e = u(tn) - uh_tn
  el2 = sqrt(sum( ∫(l2(e))dΩ ))
  @test el2 < tol
end

end #module


text:
@testset "---- Crank-Nicolson Convergence Tests, Complex Step Jacobian -----" begin
start_dir = pwd()

cd("./m1")
fname = "input_vals1.jl"
mesh, sbp, eqn, opts = solvePDE(fname)

cd("../m2")
fname = "input_vals1.jl"
mesh, sbp, eqn, opts = solvePDE(fname)

cd("../m3")
fname = "input_vals1.jl"
mesh, sbp, eqn, opts = solvePDE(fname)

cd("..")
include("calc_line.jl")

slope = calc_line()
# println("slope = ", slope)

data = readdlm("err_data.dat")
err_vals = data[:, 2]
#println("err_vals = ", err_vals)

slope_val = 2.00
slope_margin = 0.1

@test  slope  > slope_val - slope_margin
@test  slope  < slope_val + slope_margin

err_val = 0.09095728504176116
slope_fac = 1.25
# println("err_vals[1] = ", err_vals[1])
@test  err_vals[1]  > err_val/slope_fac
@test  err_vals[1]  < err_val*slope_fac

end


text:
module MRIBase

using Graphics: @mustimplement
using AbstractNFFTs
using NFFTTools # for density compensation weights in trajectory

include("Trajectories/Trajectories.jl")
include("Sequences/Sequence.jl")
include("Datatypes/Datatypes.jl")

end # module


text:
#############################################################################
# solve_dual.jl
# Calculate the dual problem of the given convex programming problems in their
# conic forms.
# Conic problems having the form
# minimize c'*x
# st       b - Ax ∈ K, where K is a cone,
# their dual problems should be
# maximize -b'*y
# st       c + Aᵀy ∈ K^*, where K^* is the dual cone of K.
#############################################################################

import MathProgBase
using Mosek
using ECOS

export solvedualMosek!

function solvedualMosek!(p::Problem)

  c, A, b, dual_var_cones, var_to_ranges, vartypes, conic_constraints = conic_problem(p)
  for i = 1:length(dual_var_cones)
    if dual_var_cones[i][1] == :Free
      dual_var_cones[i] = (:Zero, dual_var_cones[i][2])
    elseif dual_var_cones[i][1] == :Zero
      dual_var_cones[i] = (:Free, dual_var_cones[i][2])
    elseif dual_var_cones[i][1] == :ExpPrimal
      dual_var_cones[i] = (:ExpDual, dual_var_cones[i][2])
    elseif dual_var_cones[i][1] == :ExpDual
      dual_var_cones[i][1] = (:ExpPrimal, dual_var_cones[i][2])
    end
  end
  dual_constr_cones = fill((:Zero, 1:size(A, 2)),1)

  m = ConicModel(MosekSolver(MSK_IPAR_LOG = 0, MSK_DPAR_INTPNT_CO_TOL_REL_GAP = 1e-6))
  loadproblem!(m, vec(full(b)), -A', vec(full(c)), dual_constr_cones, dual_var_cones)
  optimize!(m)

  solution = try
    MathProgBase.getsolution(m)
  catch
    fill(NaN, numvar(m))
  end

  # return cones, vec(full(c)), A, vec(full(b))
  return dot(vec(full(b)), -solution)
end

# function solvedualECOS!(p::Problem)
#
#   c, A, b, dual_var_cones, var_to_ranges, vartypes, conic_constraints = conic_problem(p)
#   for i = 1:length(dual_var_cones)
#     if dual_var_cones[i][1] == :Free
#       dual_var_cones[i] = (:Zero, dual_var_cones[i][2])
#     elseif dual_var_cones[i][1] == :Zero
#       dual_var_cones[i] = (:Free, dual_var_cones[i][2])
#     elseif dual_var_cones[i][1] == :ExpPrimal
#       dual_var_cones[i] = (:ExpDual, dual_var_cones[i][2])
#     elseif dual_var_cones[i][1] == :ExpDual
#       dual_var_cones[i][1] = (:ExpPrimal, dual_var_cones[i][2])
#     end
#   end
#   dual_constr_cones = fill((:Zero, 1:size(A, 2)),1)
#   m = ConicModel(ECOSSolver(verbose=false, abstol = 1e-4))
#   loadproblem!(m, vec(full(b)), -A', vec(full(c)), dual_constr_cones, dual_var_cones)
#
#   optimize!(m)
#
#   solution = try
#     MathProgBase.getsolution(m)
#   catch
#     fill(NaN, numvar(m))
#   end
#
#   # return cones, vec(full(c)), A, vec(full(b))
#   return dot(vec(full(b)), -solution)
# end


text:
function print_parameters_to_screen(p::Param)

    #Grid information
    ti          = p.t_sim_init
    tf          = p.t_sim_final

    t_save      = p.out_every_t
    deltat      = p.deltat

    rmin        = p.rmin
    rmax        = p.rmax
    Nr          = p.rnodes

    Nt          = p.thnodes

    dr      = (rmax - rmin)/(Nr-1)
    dt      = 2*pi/Nt

    #Print information
    println("---------------------------------------------")
    println("Parameter Reader:")
    println("---------------------------------------------")
    println("Grid Information:")
    println("r      = [" , rmin , "," , rmax ,"]" )
    println("t      = [" , 0    , "," , 2*pi ,"]" )
    println("dr     = "  , dr)
    println("dtheta = "  , dt)
    println("---------------------------------------------")
    println("Time Information:")
    println("t      = [" , ti , "," , tf ,"]" )
    println("dt     = "  , deltat)
    println("---------------------------------------------")
    println("Initial Conditions:")
    println("A0     = " , p.A0)
    println("σ      = " , p.σ)
    println("r0     = " , p.r0)
    println("ω      = " , p.ω)
    println("m      = " , p.m)
    println("---------------------------------------------")
    println("Cylinder data:")
    println("Z           = " , p.Z0)
    println("epsilon     = " , p.ϵ )
    println("Omega       = " , p.Ω )
    println("---------------------------------------------")
    println("saving to file  = " ,  string(p.folder ,"/", p.fname , ".h5"))
    println("---------------------------------------------")

end


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# src/eago_optimizer/node_bb.jl
# Defines storage for a node in the B&B tree & utilities functions
#############################################################################

@enum(BranchDirection, BD_NONE, BD_NEG, BD_POS)

# Used internally to set & get variables in full problem space or just branch variables
struct FullVar end
struct BranchVar end
Base.Broadcast.broadcastable(d::FullVar) = Ref(d)
Base.Broadcast.broadcastable(d::BranchVar) = Ref(d)

"""
$(TYPEDEF)

Stores information associated with each node in Branch & Bound tree.

$(TYPEDFIELDS)
"""
struct NodeBB
    "Lower bounds of variable box."
    lower_variable_bounds::Vector{Float64}
    "Upper bounds of variable box."
    upper_variable_bounds::Vector{Float64}
    "Is dimension integer valued"
    is_integer::BitVector
    "Are all dimensions continuous or fixed"
    continuous::Bool
    "Lower bound of problem solution on nodeBB"
    lower_bound::Float64
    "Upper bound of problem solution on nodeBB"
    upper_bound::Float64
    "Depth of node in B&B tree."
    depth::Int
    "Depth of first parent in B&B tree that was continuously valued"
    cont_depth::Int
    "Unique id for each node."
    id::Int
    "Whether last branch was negative or positive in direction"
    branch_direction::BranchDirection
    "Dimension of last branch"
    last_branch::Int
    "Extent of last branch (using for psuedocost calculation)"
    branch_extent::Float64
end

# Constructors
function NodeBB(l::Vector{Float64}, u::Vector{Float64}, d::BitVector)
    NodeBB(l, u, d, any(d), -Inf, Inf, 1, -1, 1, BD_NONE, -1, 0.0)
end
NodeBB() = NodeBB(Float64[], Float64[], BitVector(), false, -Inf, Inf, 0, -1, 1, BD_NONE, -1, 0.0)
NodeBB(x::NodeBB) = NodeBB(copy(x.lower_variable_bounds), copy(x.upper_variable_bounds), copy(x.is_integer),
                           x.continuous, x.lower_bound, x.upper_bound, x.depth, x.cont_depth, x.id,
                           x.branch_direction, x.last_branch, x.branch_extent)

# Copy utilities
Base.copy(x::NodeBB) = NodeBB(copy(x.lower_variable_bounds), copy(x.upper_variable_bounds), copy(x.is_integer),
                              x.continuous, x.lower_bound, x.upper_bound, x.depth, x.cont_depth, x.id,
                              x.branch_direction, x.last_branch, x.branch_extent)

# using alternative name as to not interfere with ordering...
function uninitialized(x::NodeBB)
    flag = isempty(x.lower_variable_bounds)
    flag &= isempty(x.upper_variable_bounds)
    flag &= isempty(x.is_integer)
    flag &= x.lower_bound == -Inf
    flag &= x.upper_bound == Inf
    flag &= x.depth == 0
    flag &= x.cont_depth == -1
    flag &= x.id == 1
    return flag
end

# Access functions for broadcasting data easily
@inline lower_variable_bounds(x::NodeBB) = x.lower_variable_bounds
@inline upper_variable_bounds(x::NodeBB) = x.upper_variable_bounds
@inline lower_variable_bounds(x::NodeBB, i::Int) = x.lower_variable_bounds[i]
@inline upper_variable_bounds(x::NodeBB, i::Int) = x.upper_variable_bounds[i]
@inline lower_variable_bounds(x::NodeBB, id::Int, nid::Int) = x.lower_variable_bounds[id:nid]
@inline upper_variable_bounds(x::NodeBB, id::Int, nid::Int) = x.upper_variable_bounds[id:nid]
@inline is_integer(x::NodeBB) = x.is_integer
@inline is_integer(x::NodeBB, id::Int) = x.is_integer[id]
@inline continuous(x::NodeBB) = x.continuous
@inline lower_bound(x::NodeBB) = x.lower_bound
@inline upper_bound(x::NodeBB) = x.upper_bound
@inline depth(x::NodeBB) = x.depth
@inline cont_depth(x::NodeBB) = x.cont_depth

# Iterations Functions
Base.isless(x::NodeBB, y::NodeBB) = x.lower_bound < y.lower_bound
Base.length(x::NodeBB) = length(x.lower_variable_bounds)
function Base.isempty(x::NodeBB)
    for i = 1:length(x)
        lower = @inbounds x.lower_variable_bounds[i]
        upper = @inbounds x.upper_variable_bounds[i]
        (lower > upper) && (return true)
    end
    return false
end

"""
$(FUNCTIONNAME)

Checks that node `x` and `y` have equal domains withing a tolerance of `atol`.
"""
function same_box(x::NodeBB, y::NodeBB, r::Float64)
    (isempty(x.lower_variable_bounds) ⊻ isempty(y.lower_variable_bounds)) && (return false)
    (isempty(x.upper_variable_bounds) ⊻ isempty(y.upper_variable_bounds)) && (return false)
    for i = 1:length(x)
        ~isapprox(x.lower_variable_bounds[i], y.lower_variable_bounds[i], atol=r) && (return false)
        ~isapprox(x.upper_variable_bounds[i], y.upper_variable_bounds[i], atol=r) && (return false)
    end
    return true
end

# Compute middle & diameter
@inline diam(x::NodeBB) = upper_variable_bounds(x) - lower_variable_bounds(x)
@inline diam(x::NodeBB, i::Int) = upper_variable_bounds(x,i) - lower_variable_bounds(x,i)
@inline mid(x::NodeBB) = 0.5*(upper_variable_bounds(x) + lower_variable_bounds(x))
@inline mid(x::NodeBB, i::Int) = 0.5*(upper_variable_bounds(x,i) + lower_variable_bounds(x,i))

max_stars_repo_name: MatthewStuber/EAGO

max_issues_repo_name: MatthewStuber/EAGO

max_forks_repo_name: MatthewStuber/EAGO


text:
module TestMeasureUtils
using SimString
using Test

@testset "Test Similarity Scores" begin
    X = [1, 2, 3]
    Y = [1, 2, 4, 5]
    @test SimString.similarity_score(Dice(), X, Y)        ≈ 0.5714285714285714
    @test SimString.similarity_score(Jaccard(), X, Y)     ≈ 0.4
    @test SimString.similarity_score(Cosine(), X, Y)      ≈ 0.5773502691896258
    @test SimString.similarity_score(Overlap(), X, Y)     ≈ 0.6666666666666666
    @test SimString.similarity_score(ExactMatch(), Y, X)  == 0.0
end

@testset "Test Minimum Candidate Feature Size" begin
    @test SimString.minimum_feature_size(Dice(), 5, 1.)      == 5
    @test SimString.minimum_feature_size(Dice(), 5, 0.5)     == 2

    @test SimString.minimum_feature_size(Jaccard(), 5, 1.)   == 5
    @test SimString.minimum_feature_size(Jaccard(), 5, 0.5)  == 3

    @test SimString.minimum_feature_size(Cosine(), 5, 1.)    == 5
    @test SimString.minimum_feature_size(Cosine(), 5, 0.5)   == 2

    @test SimString.minimum_feature_size(Overlap(), 5, 1.)   == 1
    @test SimString.minimum_feature_size(Overlap(), 5, 0.5)  == 1

    @test SimString.minimum_feature_size(ExactMatch(), 5, 1.) == 5
    @test SimString.minimum_feature_size(ExactMatch(), 5, 0.5)== 5
end

@testset "Test Maximum Candidate Feature Size" begin
    db = DictDB(CharacterNGrams(3, " "))
    append!(db, ["foo", "bar", "fooo"])

    @test SimString.maximum_feature_size(Dice(), db, 5, 1.)      == 5
    @test SimString.maximum_feature_size(Dice(), db, 5, 0.5)     == 15

    @test SimString.maximum_feature_size(Jaccard(), db, 5, 1.)   == 5
    @test SimString.maximum_feature_size(Jaccard(), db, 5, 0.5)  == 10

    @test SimString.maximum_feature_size(Cosine(), db, 5, 1.)    == 5
    @test SimString.maximum_feature_size(Cosine(), db, 5, 0.5)   == 20

    @test SimString.maximum_feature_size(Overlap(), db, 5, 1.)   == 6
    @test SimString.maximum_feature_size(Overlap(), db, 5, 0.5)  == 6

    @test SimString.maximum_feature_size(ExactMatch(), db, 5, 1.)   == 5
    @test SimString.maximum_feature_size(ExactMatch(), db, 5, 0.5)  == 5
end

@testset "Test Minimum Feature Overlap" begin
    @test SimString.minimum_overlap(Dice(), 5, 5, 1.0)           == 13
    @test SimString.minimum_overlap(Dice(), 5, 20, 1.0)          == 50
    @test SimString.minimum_overlap(Dice(), 5, 5, 0.5)           == 7

    @test SimString.minimum_overlap(Jaccard(), 5, 5, 1.0)        == 5
    @test SimString.minimum_overlap(Jaccard(), 5, 20, 1.0)       == 13
    @test SimString.minimum_overlap(Jaccard(), 5, 5, 0.5)        == 4

    @test SimString.minimum_overlap(Cosine(), 5, 5, 1.0)         == 5
    @test SimString.minimum_overlap(Cosine(), 5, 20, 1.0)        == 10
    @test SimString.minimum_overlap(Cosine(), 5, 5, 0.5)         == 3

    @test SimString.minimum_overlap(Overlap(), 5, 5, 1.0)        == 5
    @test SimString.minimum_overlap(Overlap(), 5, 20, 1.0)       == 5
    @test SimString.minimum_overlap(Overlap(), 5, 5, 0.5)        == 3

    @test SimString.minimum_overlap(ExactMatch(), 5, 5, 1.0)     == 5
    @test SimString.minimum_overlap(ExactMatch(), 5, 20, 1.0)    == 5
end

end  # module


text:
using Base.Test

include("custom-set.jl")

# canonical data
@testset "empty" begin
    @test  isempty(CustomSet([]))
    @test !isempty(CustomSet([1]))
end

@testset "in (contains)" begin
    @test !(1 in CustomSet([]))
    @test   1 in CustomSet([1, 2, 3])
    @test !(4 in CustomSet([1, 2, 3]))
end

@testset "subset" begin
    @test  issubset(CustomSet([]), CustomSet([]))
    @test  issubset(CustomSet([]), CustomSet([1]))
    @test !issubset(CustomSet([1]), CustomSet([]))
    @test  issubset(CustomSet([1, 2, 3]), CustomSet([1, 2, 3]))
    @test  issubset(CustomSet([1, 2, 3]), CustomSet([4, 1, 2, 3]))
    @test !issubset(CustomSet([1, 2, 3]), CustomSet([4, 1, 3]))
end

@testset "disjoint" begin
    @test  disjoint(CustomSet([]), CustomSet([]))
    @test  disjoint(CustomSet([]), CustomSet([1]))
    @test  disjoint(CustomSet([1]), CustomSet([]))
    @test !disjoint(CustomSet([1, 2]), CustomSet([2, 3]))
    @test  disjoint(CustomSet([1, 2]), CustomSet([3, 4]))
end

@testset "equal" begin
    @test CustomSet([]) == CustomSet([])
    @test CustomSet([]) != CustomSet([1, 2, 3])
    @test CustomSet([1, 2, 3]) != CustomSet([])
    @test CustomSet([1, 2]) == CustomSet([2, 1])
    @test CustomSet([1, 2, 3]) != CustomSet([1, 2, 4])
end

@testset "add" begin
    @test begin
        custom_set = CustomSet([])
        push!(custom_set, 3)
        custom_set == CustomSet([3])
    end
    @test begin
        custom_set = CustomSet([1, 2, 4])
        push!(custom_set, 3)
        custom_set == CustomSet([1, 2, 3, 4])
    end
    @test begin
        custom_set = CustomSet([1, 2, 3])
        push!(custom_set, 3)
        custom_set == CustomSet([1, 2, 3])
    end
end

@testset "intersection" begin
    @testset "in-place" begin
        @test begin
            cs1 = CustomSet([])
            cs2 = CustomSet([])
            intersect!(cs1, cs2)
            isempty(cs1)
        end
        @test begin
            cs1 = CustomSet([])
            cs2 = CustomSet([3, 2, 5])
            intersect!(cs1, cs2)
            isempty(cs1)
        end
        @test begin
            cs1 = CustomSet([1, 2, 3, 4])
            cs2 = CustomSet([])
            intersect!(cs1, cs2)
            isempty(cs1)
        end
        @test begin
            cs1 = CustomSet([1, 2, 3])
            cs2 = CustomSet([4, 5, 6])
            intersect!(cs1, cs2)
            isempty(cs1)
        end
        @test begin
            cs1 = CustomSet([1, 2, 3, 4])
            cs2 = CustomSet([3, 2, 5])
            intersect!(cs1, cs2)
            cs1 == CustomSet([2, 3])
        end
    end
    @testset "not in-place" begin
        @test isempty(intersect(CustomSet([]), CustomSet([])))
        @test isempty(intersect(CustomSet([]), CustomSet([3, 2, 5])))
        @test isempty(intersect(CustomSet([1, 2, 3, 4]), CustomSet([])))
        @test isempty(intersect(CustomSet([1, 2, 3]), CustomSet([4, 5, 6])))
        @test intersect(CustomSet([1, 2, 3, 4]), CustomSet([3, 2, 5])) == CustomSet([2, 3])
    end
end

@testset "complement (difference)" begin
    @testset "in-place" begin
        @test begin
            cs1 = CustomSet([])
            cs2 = CustomSet([])
            complement!(cs1, cs2)
            isempty(cs1)
        end
        @test begin
            cs1 = CustomSet([])
            cs2 = CustomSet([3, 2, 5])
            complement!(cs1, cs2)
            isempty(cs1)
        end
        @test begin
            cs1 = CustomSet([1, 2, 3, 4])
            cs2 = CustomSet([])
            complement!(cs1, cs2)
            cs1 == CustomSet([1, 2, 3, 4])
        end
        @test begin
            cs1 = CustomSet([3, 2, 1])
            cs2 = CustomSet([2, 4])
            complement!(cs1, cs2)
            cs1 == CustomSet([1, 3])
        end
    end
    @testset "not in-place" begin
        @test isempty(complement(CustomSet([]), CustomSet([])))
        @test isempty(complement(CustomSet([]), CustomSet([3, 2, 5])))
        @test complement(CustomSet([1, 2, 3, 4]), CustomSet([])) == CustomSet([1, 2, 3, 4])
        @test complement(CustomSet([3, 2, 1]), CustomSet([2, 4])) == CustomSet([1, 3])
    end
end

@testset "union" begin
    @testset "in-place" begin
        @test begin
            cs1 = CustomSet([])
            cs2 = CustomSet([])
            union!(cs1, cs2)
            isempty(cs1)
        end
        @test begin
            cs1 = CustomSet([])
            cs2 = CustomSet([2])
            union!(cs1, cs2)
            cs1 == CustomSet([2])
        end
        @test begin
            cs1 = CustomSet([1, 3])
            cs2 = CustomSet([])
            union!(cs1, cs2)
            cs1 == CustomSet([1, 3])
        end
        @test begin
            cs1 = CustomSet([1, 3])
            cs2 = CustomSet([2, 3])
            union!(cs1, cs2)
            cs1 == CustomSet([3, 2, 1])
        end
    end
    @testset "not in-place" begin
        @test isempty(union(CustomSet([]), CustomSet([])))
        @test union(CustomSet([]), CustomSet([2])) == CustomSet([2])
        @test union(CustomSet([1, 3]), CustomSet([])) == CustomSet([1, 3])
        @test union(CustomSet([1, 3]), CustomSet([2, 3])) == CustomSet([3, 2, 1])
    end
end

# language specific tests
@testset "implements correct abstract type" begin
    @test CustomSet <: Base.AbstractSet
end

@testset "length" begin
    @test length(CustomSet([])) == 0
    @test length(CustomSet([1, 2, 3])) == 3
end

@testset "iterable" begin
    @test begin
        cs1 = CustomSet([1, 2, 3, 4])
        cs2 = CustomSet([])
        for element in cs1
            push!(cs2, element)
        end
        cs1 == cs2
    end
end

@testset "copy" begin
    @test begin
        cs1 = CustomSet([1, 2, 3])
        cs2 = copy(cs1)
        push!(cs1, 4)
        cs2 == CustomSet([1, 2, 3])
    end
end


text:
module OpenImageIO

import Base: show, ==

using CxxWrap

wrap_module(
    joinpath(Pkg.dir("OpenImageIO"),
             joinpath("deps", "usr", "lib", "libJlOpenImageIO"))
)

export read_image, read_image!

#-------------------------------------------------------------------------------
# TypeDesc
function show(io::IO, desc::TypeDesc)
    print(io, "TypeDesc(\"$(_print(desc))\")")
end

# OIIO has its own simple type system for describing pixel types and short
# arrays thereof.  The following defines the machinery to convert between the
# OIIO type system and julia types.
#
# Note that OIIO treats integer types more like normalized fixed point numbers
# for the purposes of color channel conversion inside the OIIO library itself.
# For example, UInt8 is treated like FixedPointNumbers.Normed{UInt8,8}.

_type_to_basetype = Dict(
    (UInt8   => TypeDesc_UINT8),
    (Int8    => TypeDesc_INT8),
    (UInt16  => TypeDesc_UINT16),
    (Int16   => TypeDesc_INT16),
    (UInt32  => TypeDesc_UINT32),
    (Int32   => TypeDesc_INT32),
    (UInt64  => TypeDesc_UINT64),
    (Int64   => TypeDesc_INT64),
    (Float16 => TypeDesc_HALF),
    (Float32 => TypeDesc_FLOAT),
    (Float64 => TypeDesc_DOUBLE),
    (String  => TypeDesc_STRING),
    (Void    => TypeDesc_PTR),
)

_basetype_to_type = Dict(b=>a for (a,b) in _type_to_basetype)

const BaseType = Union{[a for (a,b) in _type_to_basetype]...}
const ColorChannelType = Union{UInt8,Int8,UInt16,Int16,UInt32,Int32,UInt64,Int64,Float16,Float32,Float64}

function TypeDesc(basetype, aggregate=TypeDesc_SCALAR,
                  vecsem=TypeDesc_NOXFORM, arraylen=0)
    TypeDesc(basetype, aggregate, vecsem, arraylen)
end

function TypeDesc(basetype::TypeDesc_BASETYPE, arraylen::Integer)
    TypeDesc(basetype, TypeDesc_SCALAR, TypeDesc_NOXFORM, arraylen)
end

TypeDesc(::Type{T}) where {T<:BaseType} = TypeDesc(_type_to_basetype[T])
TypeDesc(::Type{NTuple{N,T}}) where {N,T<:BaseType} = TypeDesc(_type_to_basetype[T], N)

# Get image channel type -
function channeltype(desc::TypeDesc)
    if aggregate(desc) != TypeDesc_SCALAR || vecsemantics(desc) != TypeDesc_NOXFORM || arraylen(desc) > 0
        throw(ArgumentError("Cannot convert $desc to an image channel scalar"))
    end
    _basetype_to_type[basetype(desc)]
end

#channeltype(::Array{NTuple{<:Any,T}}) where {T} = T
#imageshape(::Array{NTuple{N,<:Any}, M}) where {N,M} = N
#nchannels(::Array{T,2}) where {N} = N

#-------------------------------------------------------------------------------
# ImageSpec
function show(io::IO, spec::ImageSpec)
    print(io, "ImageSpec$(serialize(spec, ImageSpec_SerialText, ImageSpec_SerialDetailedHuman))")
end

#-------------------------------------------------------------------------------
# Convenience functions for reading from ImageInput

# Compute native color tuple type for the image.  For now, we only allow
# homogenous color tuples.
function native_color_tuple_type(spec::ImageSpec)
    N = Int(nchannels(spec))
    NTuple{N, channeltype(format(spec))}
end

"""
    alloc_image(spec::ImageSpec)

Allocate an array to hold color data for an image described by `spec` with a
color channel type which is the same as the native image storage.  Elements of
the array are colors which are treated as simple tuples of the channel type.
The resulting array is either one, two, or three dimensional, depending on the
dimensionality of the underlying image data.
"""
function alloc_image(spec::ImageSpec)
    C = native_color_tuple_type(spec)
    if depth(spec) > 1
        Array{C,3}(height(spec), width(spec), depth(spec))
    elseif width(spec) > 1
        Array{C,2}(height(spec), width(spec))
    else
        # Is this sensible?  Would width=1 normally be used to
        # represent a one dimensional image, rather height=1?
        Array{C,1}(height(spec))
    end
end

"""
    read_image!(input::ImageInput, data)

Read pixel data for the current subimage of `input` into the array `data`
containing color tuples.  The size of `data` should match with the size of the
input image.  Supported color channel types are any subtypes of
`OIIO.ColorChannelType`.
"""
function read_image!(input::ImageInput,
                     data::Array{NTuple{Nchan,Tchan}}) where {Nchan,Tchan}
    s = spec(input)
    if nchannels(s) != Nchan
        throw(ArgumentError("Given number of color channels $Nchan is different number of image channels $(nchannels(s))"))
    end
    if !(Tchan <: ColorChannelType)
        throw(ArgumentError("Color channel type $Tchan is not supported by OpenImageIO"))
    end
    data_size = (size(data,1), size(data,2), size(data,3))
    expected_size = (height(s), width(s), depth(s))
    if data_size != expected_size
        throw(ArgumentError("Array size $data_size different from image size $expected_size"))
    end
    # Julia column-major array ordering is opposite from what would be computed
    # by OIIO's AutoStride mechanism.
    ystride = sizeof(NTuple{Nchan,Tchan})
    xstride = ystride*size(data,1)
    zstride = xstride*size(data,2)
    success = unsafe_read_image!(input, 0, Nchan, TypeDesc(Tchan),
                                 convert(Ptr{Void},pointer(data)),
                                 xstride, ystride, zstride)
    if !success
        throw(ErrorException(geterror(input)))
    end
    data
end

read_image(input::ImageInput) = read_image!(input, alloc_image(spec(input)))

end


text:
"""Bitstring function with a split-mode that splits the bitstring into
sign, exponent and significant bits."""
function Base.bitstring(x::Float16,mode::Symbol)
    if mode == :split
        bstr = bitstring(x)
        return "$(bstr[1]) $(bstr[2:6]) $(bstr[7:end])"
    else
        bitstring(x)
    end
end

"""Bitstring function with a split-mode that splits the bitstring into
sign, exponent and significant bits."""
function Base.bitstring(x::Float32,mode::Symbol)
    if mode == :split
        bstr = bitstring(x)
        return "$(bstr[1]) $(bstr[2:9]) $(bstr[10:end])"
    else
        bitstring(x)
    end
end

"""Bitstring function with a split-mode that splits the bitstring into
sign, exponent and significant bits."""
function Base.bitstring(x::Float64,mode::Symbol)
    if mode == :split
        bstr = bitstring(x)
        return "$(bstr[1]) $(bstr[2:12]) $(bstr[13:end])"
    else
        bitstring(x)
    end
end


text:
# Chirp Z-Transform

mutable struct CZTFilter{T}
    dataLen::Int #Expected data length
    nPhones::Int #Expected number of phones
    fftIter::Vector{Int} #preallocated vector for fft iteration
    phonesIter::Vector{Int} #preallocated vector for phones iteration
    outLen::Int #Output length
    wk2::Vector{Complex{T}}
    Awk2::Vector{Complex{T}}
    Fwk2::Vector{Complex{T}}
    tempData::Array{Complex{T}}
    nFFT::Int
    CZTFilter{T}() where T = new{T}()
end

function (czt::CZTFilter{T})(dataIn::Array, dataOut::Array{Complex{T}}) where {T <: Real}
    dataLen = size(dataIn,1)
    if dataLen!=czt.dataLen
        error("Data Signal length: $dataLen Filter expected: $(czt.dataLen)")
    end
    if czt.nPhones != size(dataIn,2)
        error("Number of phones: $(size(dataIn,2)). Expected: $(czt.nPhones)")
    end
    vdI = view(dataIn,1:dataLen,:)
    fill!(czt.tempData, 0.0 + 0.0im)
    czt.tempData[1:dataLen,:] = vdI .* czt.Awk2; #zero out temp data
    fft!(czt.tempData,1)
    # fwk = czt.Fwk2
    for j in czt.phonesIter
      czt.tempData[:,j] .*= czt.Fwk2
    end
    ifft!(czt.tempData,1)
    vtD = view(czt.tempData, czt.dataLen:(czt.dataLen+czt.outLen-1),:)
    copyto!(dataOut, vtD)
    dataOut .*= czt.wk2
    nothing
end

# function (czt::CZTFilter{T})(dataIn::Array, dataOut::Array{Complex{T}}) where {T <: Real}
#     dataLen = size(dataIn,1)
#     if dataLen!=czt.dataLen
#         error("Data Signal length: $dataLen Filter expected: $(czt.dataLen)")
#     end
#     if czt.nPhones != size(dataIn,2)
#         error("Number of phones: $(size(dataIn,2)). Expected: $(czt.nPhones)")
#     end
#     for i in 1:czt.dataLen
#         for j in czt.phonesIter
#             czt.tempData[i,j] = dataIn[i,j] * czt.Awk2[i]
#         end
#     end
#     fft!(czt.tempData,1)
#     for i in czt.fftIter
#         for j in czt.phonesIter
#           czt.tempData[i,j] = czt.tempData[i,j] * czt.Fwk2[i]
#         end
#     end
#     ifft!(czt.tempData,1)
#     for i in czt.dataLen:(czt.dataLen+czt.outLen-1)
#       for j in czt.phonesIter
#         dataOut[i-czt.dataLen+1,j] = czt.tempData[i,j]*czt.wk2[i-czt.dataLen+1]
#       end
#     end
#     nothing
# end

#orig - starting point in the plane
#theoretically allocating all memory here
function prepCZTFilter(n::Int, nPhones::Int, w::Complex{T}, m::Int=-1, a::Complex{T}=1.0+0im) where {T}
    # , T::Type=Float64)
    out = CZTFilter{T}()
    out.dataLen=n
    if m<0
        m=n
    end
    out.outLen=m
    out.nPhones = nPhones
    out.phonesIter = 1:out.nPhones #linspace(1,nPhones,nPhones)
    out.nFFT = nextpow(2,n+m)
    out.fftIter = 1:out.nFFT # linspace(1,out.nFFT,out.nFFT)
    out.tempData = zeros(Complex{T},out.nFFT,nPhones)
    k = LinRange(0,max(n,m)-1,max(n,m))
    wk2 = w.^(k.^2/2); #@show size(wk2)
    out.Awk2 = ((a.^-k) .* wk2)[1:n]
    out.Fwk2 = zeros(Complex{T},out.nFFT)
    out.Fwk2[1:n] = 1.0 ./wk2[n:-1:1]
    out.Fwk2[n:m+n-1] = 1.0 ./wk2[1:m]
    #out.Fwk2[1:n+1] = 1./append!(wk2[n:-1:2], wk2[1:m])    #alternative for two lines above
    fft!(out.Fwk2)
    out.wk2 = wk2[1:m]
    return out
end

#


text:
# In distributed memory mode, the model may be too large to load into memory on
# a single node. It is therefore split by one of its dimensions into blobs.
#
# All matrix operations are done on the node where required, just that not all
# blobs are loaded at once at any time. This is more efficient for the ALS
# algorithm, because:
# - ALS computation is already distributed over users and items
# - it is more efficient to transfer a full blob once to a node, instead of
#   communicating to a remote node once for each step
type DistModel <: Model
    UFile::FileSpec
    PFile::FileSpec
    nfactors::Int
    lambda::Float64
    U::Nullable{ModelFactor}
    P::Nullable{ModelFactor}
    lambdaI::Nullable{ModelFactor}
    Pinv::Nullable{ModelFactor}
end

nusers(model::DistModel) = size(get(model.U))[2]
nitems(model::DistModel) = size(get(model.P))[2]
nfactors(model::DistModel) = model.nfactors

share!(model::DistModel) = nothing
localize!(model::DistModel) = nothing

function clear(model::DistModel)
    model.U = nothing
    model.P = nothing
    model.lambdaI = nothing
    model.Pinv = nothing
end

function ensure_loaded(model::DistModel)
    model.U = read_input(model.UFile)
    model.P = read_input(model.PFile)
    model.lambdaI = model.lambda * eye(model.nfactors)
    nothing
end

# TODO
#function pinv(model::DistModel)
#end
#mul_pinv(model::DistModel, v) = v * pinv(model)

function prep{TI<:DistInputs}(inp::TI, nfacts::Int, lambda::Float64, model_dir::AbstractString, maxcache::Int=def_cache())
    ensure_loaded(inp)
    t1 = time()
    @logmsg("preparing inputs...")

    nu = nusers(inp)
    ni = nitems(inp)

    Udir = DenseBlobs(joinpath(model_dir, "U"); maxcache=maxcache)
    Pdir = DenseBlobs(joinpath(model_dir, "P"); maxcache=maxcache)
    isdir(model_dir) || mkdir(model_dir)
    isdir(Udir.name) || mkdir(Udir.name)
    isdir(Pdir.name) || mkdir(Pdir.name)

    Usz = (nfacts, nu)
    Psz = (nfacts, ni)
    U = create(Udir, Float64, Usz, zeros, min(_max_items(Float64,Usz), ceil(Int, nu/nworkers())))
    P = create(Pdir, Float64, Psz, rand, min(_max_items(Float64,Psz), ceil(Int, ni/nworkers())))

    for idx in 1:ni
        P[1,idx] = mean(all_user_ratings(inp, idx))
    end
    save(U)
    save(P)

    lambdaI = lambda * eye(nfacts)
    model = DistModel(Udir, Pdir, nfacts, lambda, U, P, lambdaI, nothing)

    t2 = time()
    @logmsg("prep time: $(t2-t1)")
    model
end


text:
#####
##### Halo filling for flux, periodic, and no-penetration boundary conditions.
#####

# For flux boundary conditions we fill halos as for a *no-flux* boundary condition, and add the
# flux divergence associated with the flux boundary condition in a separate step. Note that
# ranges are used to reference the data copied into halos, as this produces views of the correct
# dimension (eg size = (1, Ny, Nz) for the west halos).

  _fill_west_halo!(c, ::FBC, H, N) = @views @. c.parent[1:H, :, :] = c.parent[1+H:1+H,  :, :]
 _fill_south_halo!(c, ::FBC, H, N) = @views @. c.parent[:, 1:H, :] = c.parent[:, 1+H:1+H,  :]
_fill_bottom_halo!(c, ::FBC, H, N) = @views @. c.parent[:, :, 1:H] = c.parent[:, :,  1+H:1+H]

 _fill_east_halo!(c, ::FBC, H, N) = @views @. c.parent[N+H+1:N+2H, :, :] = c.parent[N+H:N+H, :, :]
_fill_north_halo!(c, ::FBC, H, N) = @views @. c.parent[:, N+H+1:N+2H, :] = c.parent[:, N+H:N+H, :]
  _fill_top_halo!(c, ::FBC, H, N) = @views @. c.parent[:, :, N+H+1:N+2H] = c.parent[:, :, N+H:N+H]

# Periodic boundary conditions
  _fill_west_halo!(c, ::PBC, H, N) = @views @. c.parent[1:H, :, :] = c.parent[N+1:N+H, :, :]
 _fill_south_halo!(c, ::PBC, H, N) = @views @. c.parent[:, 1:H, :] = c.parent[:, N+1:N+H, :]
_fill_bottom_halo!(c, ::PBC, H, N) = @views @. c.parent[:, :, 1:H] = c.parent[:, :, N+1:N+H]

 _fill_east_halo!(c, ::PBC, H, N) = @views @. c.parent[N+H+1:N+2H, :, :] = c.parent[1+H:2H, :, :]
_fill_north_halo!(c, ::PBC, H, N) = @views @. c.parent[:, N+H+1:N+2H, :] = c.parent[:, 1+H:2H, :]
  _fill_top_halo!(c, ::PBC, H, N) = @views @. c.parent[:, :, N+H+1:N+2H] = c.parent[:, :, 1+H:2H]

# Recall that, by convention, the first grid point (k=1) in an array with a no-penetration boundary
# condition lies on the boundary, where as the last grid point (k=Nz) lies in the domain.

  _fill_west_halo!(c, ::NPBC, H, N) = @views @. c.parent[1:1+H, :, :] = 0
 _fill_south_halo!(c, ::NPBC, H, N) = @views @. c.parent[:, 1:1+H, :] = 0
_fill_bottom_halo!(c, ::NPBC, H, N) = @views @. c.parent[:, :, 1:1+H] = 0

 _fill_east_halo!(c, ::NPBC, H, N) = @views @. c.parent[N+H+1:N+2H, :, :] = 0
_fill_north_halo!(c, ::NPBC, H, N) = @views @. c.parent[:, N+H+1:N+2H, :] = 0
  _fill_top_halo!(c, ::NPBC, H, N) = @views @. c.parent[:, :, N+H+1:N+2H] = 0

# Generate functions that implement flux, periodic, and no-penetration boundary conditions
sides = (:west, :east, :south, :north, :top, :bottom)
coords = (:x, :x, :y, :y, :z, :z)

for (x, side) in zip(coords, sides)
    outername = Symbol(:fill_, side, :_halo!)
    innername = Symbol(:_fill_, side, :_halo!)
    H = Symbol(:H, x)
    N = Symbol(:N, x)
    @eval begin
        $outername(c, bc::Union{FBC, PBC, NPBC}, arch::AbstractArchitecture, grid::AbstractGrid, args...) =
            $innername(c, bc, grid.$(H), grid.$(N))
    end
end

#####
##### Halo filling for value and gradient boundary conditions
#####

@inline linearly_extrapolate(c₀, ∇c, Δ) = c₀ + ∇c * Δ

@inline bottom_gradient(bc::GBC, c¹, Δ, i, j, args...) = getbc(bc, i, j, args...)
@inline top_gradient(bc::GBC, cᴺ, Δ, i, j, args...) = getbc(bc, i, j, args...)

@inline south_gradient(bc::GBC, c¹, Δ, i, k, args...) = getbc(bc, i, k, args...)
@inline north_gradient(bc::GBC, cᴺ, Δ, i, k, args...) = getbc(bc, i, k, args...)

@inline bottom_gradient(bc::VBC, c¹, Δ, i, j, args...) = ( c¹ - getbc(bc, i, j, args...) ) / (Δ/2)
@inline top_gradient(bc::VBC, cᴺ, Δ, i, j, args...) =    ( getbc(bc, i, j, args...) - cᴺ ) / (Δ/2)

@inline left_gradient(bc::VBC, c¹, Δ, i, k, args...) =  ( c¹ - getbc(bc, i, k, args...) ) / (Δ/2)
@inline right_gradient(bc::VBC, cᴺ, Δ, i, k, args...) = ( getbc(bc, i, k, args...) - cᴺ ) / (Δ/2)

function fill_bottom_halo!(c, bc::Union{VBC, GBC}, arch, grid, args...)
    @launch device(arch) config=launch_config(grid, :xy) _fill_bottom_halo!(c, bc, grid, args...)
    return nothing
end

function fill_top_halo!(c, bc::Union{VBC, GBC}, arch, grid, args...)
    @launch device(arch) config=launch_config(grid, :xy) _fill_top_halo!(c, bc, grid, args...)
    return nothing
end

function fill_south_halo!(c, bc::Union{VBC, GBC}, arch, grid, args...)
    @launch device(arch) config=launch_config(grid, :xz) _fill_south_halo!(c, bc, grid, args...)
    return nothing
end

function fill_north_halo!(c, bc::Union{VBC, GBC}, arch, grid, args...)
    @launch device(arch) config=launch_config(grid, :xz) _fill_north_halo!(c, bc, grid, args...)
    return nothing
end

function _fill_bottom_halo!(c, bc::Union{VBC, GBC}, grid, args...)
    @loop_xy i j grid begin
        @inbounds ∇c = bottom_gradient(bc, c[i, j, 1], grid.Δz, i, j, grid, args...)
        @unroll for k in (1 - grid.Hz):0
            Δ = (k - 1) * grid.Δz  # separation between bottom grid cell and halo is negative
            @inbounds c[i, j, k] = linearly_extrapolate(c[i, j, 1], ∇c, Δ)
        end
    end
    return nothing
end

function _fill_top_halo!(c, bc::Union{VBC, GBC}, grid, args...)
    @loop_xy i j grid begin
        @inbounds ∇c = top_gradient(bc, c[i, j, grid.Nz], grid.Δz, i, j, grid, args...)
        @unroll for k in (grid.Nz + 1) : (grid.Nz + grid.Hz)
            Δ = (k - grid.Nz) * grid.Δz
            @inbounds c[i, j, k] = linearly_extrapolate(c[i, j, grid.Nz], ∇c, Δ)
        end
    end
    return nothing
end

function _fill_south_halo!(c, bc::Union{VBC, GBC}, grid, args...)
    @loop_xz i k grid begin
        @inbounds ∇c = south_gradient(bc, c[i, 1, k], grid.Δy, i, k, grid, args...)
        @unroll for j in (1 - grid.Hy):0
            Δ = (j - 1) * grid.Δy  # separation between southern-most grid cell and halo is negative
            @inbounds c[i, j, k] = linearly_extrapolate(c[i, 1, k], ∇c, Δ)
        end
    end
    return nothing
end

function _fill_north_halo!(c, bc::Union{VBC, GBC}, grid, args...)
    @loop_xz i k grid begin
        @inbounds ∇c = north_gradient(bc, c[i, grid.Ny, k], grid.Δy, i, k, grid, args...)
        @unroll for j in (grid.Ny + 1) : (grid.Ny + grid.Hy)
            Δ = (k - grid.Ny) * grid.Δy
            @inbounds c[i, j, k] = linearly_extrapolate(c[i, grid.Ny, k], ∇c, Δ)
        end
    end
    return nothing
end

#####
##### General halo filling functions
#####

"Fill halo regions in x, y, and z for a given field."
function fill_halo_regions!(c::AbstractArray, fieldbcs, arch, grid, args...)
      fill_west_halo!(c, fieldbcs.x.left,  arch, grid, args...)
      fill_east_halo!(c, fieldbcs.x.right, arch, grid, args...)

     fill_south_halo!(c, fieldbcs.y.left,  arch, grid, args...)
     fill_north_halo!(c, fieldbcs.y.right, arch, grid, args...)

     fill_bottom_halo!(c, fieldbcs.z.bottom, arch, grid, args...)
        fill_top_halo!(c, fieldbcs.z.top,    arch, grid, args...)
    return nothing
end

"""
    fill_halo_regions!(fields, bcs, arch, grid)

Fill halo regions for all fields in the `NamedTuple` `fields` according
to the corresponding `NamedTuple` of `bcs`.
"""
function fill_halo_regions!(fields::NamedTuple{S}, bcs::NamedTuple{S}, arch, grid, args...) where S
    for (field, fieldbcs) in zip(fields, bcs)
        fill_halo_regions!(field, fieldbcs, arch, grid, args...)
    end
    return nothing
end

"""
    fill_halo_regions!(fields, bcs, arch, grid)

Fill halo regions for each field in the tuple `fields` according
to the single instance of `FieldBoundaryConditions` in `bcs`, possibly recursing into
`fields` if it is a nested tuple-of-tuples.
"""
function fill_halo_regions!(fields::Union{Tuple, NamedTuple}, bcs::FieldBoundaryConditions, arch, grid, args...)
    for field in fields
        fill_halo_regions!(field, bcs, arch, grid, args...)
    end
    return nothing
end

fill_halo_regions!(::Nothing, args...) = nothing

#####
##### Halo zeroing functions
#####

  zero_west_halo!(c, H, N) = @views @. c[1:H, :, :] = 0
 zero_south_halo!(c, H, N) = @views @. c[:, 1:H, :] = 0
zero_bottom_halo!(c, H, N) = @views @. c[:, :, 1:H] = 0

 zero_east_halo!(c, H, N) = @views @. c[N+H+1:N+2H, :, :] = 0
zero_north_halo!(c, H, N) = @views @. c[:, N+H+1:N+2H, :] = 0
  zero_top_halo!(c, H, N) = @views @. c[:, :, N+H+1:N+2H] = 0

function zero_halo_regions!(c::AbstractArray, grid)
      zero_west_halo!(c, grid.Hx, grid.Nx)
      zero_east_halo!(c, grid.Hx, grid.Nx)
     zero_south_halo!(c, grid.Hy, grid.Ny)
     zero_north_halo!(c, grid.Hy, grid.Ny)
       zero_top_halo!(c, grid.Hz, grid.Nz)
    zero_bottom_halo!(c, grid.Hz, grid.Nz)
    return nothing
end

function zero_halo_regions!(fields::Tuple, grid)
    for field in fields
        zero_halo_regions!(field, grid)
    end
    return nothing
end


text:
######## SET UNIT CLASSES #######
#=
Here, we set classes of unit names for the quantities to be defined later.
The second argument is the preferred units in which to display the quantity.
The third argument is the dimension of the quantity.
=#
@displayedunits TimeType           "s"      𝐓
@displayedunits TemperatureType    "K"      𝚯
@displayedunits MassType           "kg"     𝐌
@displayedunits LengthType         "m"      𝐋
@displayedunits PressureType       "Pa"     𝐌/𝐓^2/𝐋
@displayedunits DensityType        "kg/m^3" 𝐌/𝐋^3
@displayedunits SpecificEnergyType "J/kg"   𝐋^2/𝐓^2
@displayedunits SpecificHeatType   "J/kg/K" 𝐋^2/𝐓^2/𝚯
@displayedunits VelocityType       "m/s"    𝐋/𝐓
@displayedunits AccelerationType   "m/s^2"    𝐋/𝐓^2
@displayedunits AngularVelocityType "rad/s" 𝐓^-1
@displayedunits AreaType           "m^2"    𝐋^2
@displayedunits VolumeType         "m^3"    𝐋^3
@displayedunits MassFlowRateType   "kg/s"   𝐌/𝐓
@displayedunits FlowRateType       "m^3/s"  𝐋^3/𝐓
@displayedunits KinematicViscosityType "m^2/s" 𝐋^2/𝐓
@displayedunits ViscosityType      "kg/m/s" 𝐌/𝐓/𝐋
@displayedunits ForcePerLengthType "N/m"    𝐌/𝐓^2
@displayedunits ForceType          "N"      𝐌*𝐋/𝐓^2
@displayedunits SpecificWeightType "N/m^3"  𝐌/𝐋^2/𝐓^2
@displayedunits AreaMomentType      "m^4"    𝐋^4
@displayedunits MomentType          "N*m"    𝐌*𝐋^2/𝐓^2

#=
Quantities and properties

Note: to add a new quantity, you define it here, assigning the units from
the list above of QtyType
=#
@dimvar SpecificHeatPressure   SpecificHeatType
@dimvar SpecificHeatVolume     SpecificHeatType
@dimvar GasConstant            SpecificHeatType
@dimvar Viscosity              ViscosityType
@dimvar KinematicViscosity     KinematicViscosityType
@dimvar SurfaceTension         ForcePerLengthType
@dimvar SpecificWeight         SpecificWeightType
@dimvar VaporPressure          PressureType
@dimvar BulkModulus            PressureType

@nondimvar SpecificHeatRatio
@nondimvar SpecificGravity

# Dimensional variables
@dimvar Pressure               PressureType
@dimvar PressureDifference     PressureType
@dimvar GaugePressure          PressureType
@dimvar StagnationPressure     PressureType
@dimvar Density                DensityType
@dimvar StagnationDensity      DensityType
@dimvar Temperature            TemperatureType
@dimvar StagnationTemperature  TemperatureType
@dimvar InternalEnergy         SpecificEnergyType
@dimvar StagnationInternalEnergy SpecificEnergyType
@dimvar Enthalpy               SpecificEnergyType
@dimvar StagnationEnthalpy     SpecificEnergyType
@dimvar SoundSpeed             VelocityType
@dimvar StagnationSoundSpeed   VelocityType
@dimvar Velocity               VelocityType
@dimvar UVelocity              VelocityType
@dimvar VVelocity              VelocityType
@dimvar WVelocity              VelocityType
@dimvar AngularVelocity        AngularVelocityType
@dimvar Acceleration           AccelerationType
@dimvar Gravity                AccelerationType
@dimvar Entropy                SpecificHeatType
@dimvar MassFlowRate           MassFlowRateType
@dimvar HeatFlux               SpecificEnergyType
@dimvar VolumeFlowRate         FlowRateType
@dimvar Force                  ForceType
@dimvar ForcePerDepth          ForcePerLengthType
@dimvar Moment                 MomentType
@dimvar Time                   TimeType
@dimvar Head                   LengthType

@dimvar Area                   AreaType
@dimvar Diameter               LengthType
@dimvar Length                 LengthType
@dimvar Height                 LengthType
@dimvar Depth                  LengthType
@dimvar Volume                 VolumeType
@dimvar RoughnessHeight        LengthType

# Non-dimensional variables
@nondimvar MachNumber
@nondimvar ReynoldsNumber
@nondimvar DragCoefficient
@nondimvar LiftCoefficient
@nondimvar PressureCoefficient
@nondimvar FrictionFactor

# a generic catchall
@nondimvar DimensionlessParameter

@nondimvar FLOverD
@nondimvar PressureRatio
@nondimvar StagnationPressureRatio
@nondimvar TemperatureRatio
@nondimvar DensityRatio
@nondimvar VelocityRatio
@nondimvar AreaRatio
@nondimvar RoughnessRatio

max_stars_repo_name: UCLAMAEThreads/ThermofluidQuantities

max_issues_repo_name: UCLAMAEThreads/ThermofluidQuantities

max_forks_repo_name: UCLAMAEThreads/ThermofluidQuantities


text:
################################################################################
#
#               ZZModUInt.jl : Z/nZ modelled with UInts
#
# This file is part of Hecke.
#
# Copyright (c) 2015, 2016: Claus Fieker, Tommy Hofmann
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# * Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
#
#  Copyright (C) 2015, 2016 Tommy Hofmann, Claus Fieker
#
################################################################################

export ZZModUInt, UIntMod

################################################################################
#
#  Basics
#
################################################################################

elem_type(::ZZModUInt) = UIntMod

parent_type(::Type{UIntMod}) = ZZModUInt

zero(R::ZZModUInt) = R(0)

Base.promote_rule(::Type{UIntMod}, ::Type{Int}) = UIntMod

parent(x::UIntMod) = x.parent

@inline function iszero(x::UIntMod)
  return x.m == 0
end

################################################################################
#
#  String I/O
#
################################################################################

Base.show(io::IO, x::UIntMod) = Base.show(io, x.m)

function Base.show(io::IO, R::ZZModUInt)
  print(io, "Residue ring of ZZ of size $(R.mod.n) (UInt)\n")
end

################################################################################
#
#  Binary operations
#
################################################################################

for (fn, fc) in ((:+, "nmod_add"), (:*, "nmod_mul"), (:-, "nmod_sub"),
      (://, "nmod_div"))
  @eval begin
    function ($fn)(x::UIntMod, y::UIntMod)
      z = UIntMod(parent(x), (ccall(($fc, :libflint), UInt,
              (UInt, UInt, nmod_struct),
              x.m, y.m, x.parent.mod)))
      return z
    end
  end
end

divexact(x::UIntMod, y::UIntMod) = //(x, y)

################################################################################
#
#  Unary operations
#
################################################################################

function -(x::UIntMod)
  z = UIntMod(parent(x), ccall((:nmod_neg, :libflint), UInt,
              (UInt, nmod_struct), x.m, x.parent.mod))
  return z
end

################################################################################
#
#  Comparison
#
################################################################################

function ==(x::UIntMod, y::UIntMod)
  x.parent == y.parent && x.m == y.m
end

################################################################################
#
#  Parent object overloading
#
################################################################################

function call(R::ZZModUInt)
  return ZZModUInt(R)
end

function call(R::ZZModUInt, n::UInt)
  return UIntMod(R, n % R.mod.n)
end

function call(R::ZZModUInt, n::fmpz)
  u = ccall((:fmpz_tdiv_ui, :libflint), UInt, (Ptr{fmpz}, UInt), &n, R.mod.n)
  return UIntMod(R, u)
end

call(R::ZZModUInt, n::Int) = R(UInt(n))

function call(R::ZZModUInt, x::GenRes{fmpz})
  R.mod.n != parent(x).modulus && error("moduli must be equal")
  return R(UInt(x.data))
end


text:
include("misc.jl")

function rbfBasis(X1, X2, sigma)
  return exp.(-distancesSquared(X1,X2)/(2*sigma^2))
end

function leastSquaresRBFL2(X, y, l, sigma)
	# Make RBF basis
	n = size(X,1)
	Z = rbfBasis(X, X, sigma)

	# Find regression weights minimizing squared error with L2 regularization
	w = (Z'*Z + eye(size(Z,2))*l)\(Z'*y)

	display(size(w))

	# Make linear prediction function
	predict(Xtilde) = rbfBasis(Xtilde, X, sigma)*w

	# Return model
	return LinearModel(predict,w)
end


text:
module IsingModels

using Random

export IsingModel2D,randomdata!,allup!,alldown!
export mcsweeps,mcsweep
export glaubersweeps
export calc_magnetization,calc_avg_magnetization

include("ising.jl")
include("callbacks.jl")
include("mc.jl")
include("glauber.jl")
include("observables.jl")

end


text:
using DelayEmbeddings, OrdinaryDiffEq

ti = time()

const diffeq = (alg = Vern9(), atol = 1e-9, rtol = 1e-9, maxiters = typemax(Int))

include("dataset_tests.jl")
include("reconstruction_tests.jl")
include("R_delay.jl")
include("R_dimension.jl")

ti = time() - ti
println("\nTest took total time of:")
println(round(ti, digits=3), " seconds or ", round(ti/60, digits=3), " minutes")


text:
# Configuration for nearly neutral simuation
# Example run from command line:  julia run.jl examples/ia_example2
function dfemod(x)
  dfe_mod(x,fit_inc=1.2)
end
const nn_simtype = 1    # nearly neutral infinite alleles
@everywhere const N_list     = [40,80]     # sample size n list, popsize N may be larger
const N_mu_list= [2.0,4.0]   # Mutation rate as a multiple of 1.0/N
const L = 1           # number of loci
const ngens = 10      # Generations after burn-in
const burn_in= 2.0    # generations of burn_in as a multiple of 1/mu
const use_poplist=false   # save a list of populations as a check of correctness
dfe=dfemod
dfe_str = "dfe_mod fit_inc: 1.2"

max_stars_repo_name: um-tech-evolution/nearly-neutral

max_issues_repo_name: um-tech-evolution/nearly-neutral

max_forks_repo_name: um-tech-evolution/nearly-neutral


text:
# Elements of Chemical Reaction Engineering
# Fifth Edition
# H. SCOTT FOGLER
# Chapter 13: Unsteady-State Nonisothermal Reactor Design
# Section 13.5: Nonisothermal Multiple Reactions
# Example 13–5 Multiple Reactions in a Semibatch Reactor
# p. 658

function semibatch_reactor(; store_results=false::Bool)
    datadir = nothing
    if store_results
        datadir = generate_data_subdir(@__FILE__)
    end

    # initial conditions and timepoints
    t0 = 0.0f0
    tf = 0.4f0  # Bradfoard uses 0.4
    Δt = 0.01f0
    tspan = (t0, tf)

    # state: CA, CB, CC, T, Vol
    u0 = [1.0f0, 0.0f0, 0.0f0, 290.0f0, 100.0f0]

    # control constraints
    # F = volumetric flow rate
    # V = exchanger temperature
    # F = 240 & V = 298 in Fogler's book
    # F ∈ (0, 250) & V ∈ (200, 500) in Bradford 2017
    control_ranges = [(100.0f0, 700.0f0), (0.0f0, 400.0f0)]

    # state constraints
    # T ∈ (0, 420]
    # Vol ∈ (0, 200]
    T_up = 380.0f0
    V_up = 100.0f0

    system_params = (
        CpA=30.0f0,
        CpB=60.0f0,
        CpC=20.0f0,
        CpH2SO4=35.0f0,
        N0H2S04=100.0f0,
        T0=305.0f0,
        CA0=4.0f0,
        HRA=-6500.0f0,
        HRB=8000.0f0,
        E1A=9500.0f0 / 1.987f0,
        E2B=7000.0f0 / 1.987f0,
        A1=1.25f0,
        A2=0.08f0,
        UA=35000.0f0,  # 45000,  Bradford value
        Tr1=320.0f0,  # 420, Bradford value
        Tr2=290.0f0,  # 400, Bradford value
    )

    function system!(du, u, p, t, controller; input=:state)
        @argcheck input in (:state, :time)

        CpA, CpB, CpC, CpH2SO4, N0H2S04, T0, CA0, HRA, HRB, E1A, E2B, A1, A2, UA, Tr1, Tr2 = values(
            system_params
        )

        # neural network outputs controls taken by the system
        CA, CB, CC, T, Vol = u  # states
        # controls
        if input == :state
            c_F, c_T = controller(u, p)
        elseif input == :time
            c_F, c_T = controller(t, p)
        end
        k1A = A1 * exp(E1A * ((1.0f0 / Tr1) - (1.0f0 / T)))
        k2B = A2 * exp(E2B * ((1.0f0 / Tr2) - (1.0f0 / T)))

        k1CA = k1A * CA
        k2CB = k2B * CB
        F_Vol = c_F / Vol

        ra = -k1CA
        rb = 0.5f0 * k1CA - k2CB
        rc = 3.0f0 * k2CB

        num =
            UA * (c_T - T) - CA0 * c_F * CpA * (T - T0) +
            (HRA * (-k1CA) + HRB * (-k2CB)) * Vol
        den = (CA * CpA + CpB * CB + CpC * CC) * Vol + N0H2S04 * CpH2SO4

        # dynamics of the controlled system
        dCA = ra + (CA0 - CA) * F_Vol
        dCB = rb - CB * F_Vol
        dCC = rc - CC * F_Vol
        dT = num / den
        dVol = c_F

        # update in-place
        @inbounds begin
            du[1] = dCA
            du[2] = dCB
            du[3] = dCC
            du[4] = dT
            du[5] = dVol
        end
    end

    # simulate the system with constant controls as in Fogler's
    # to reproduce his results and verify correctness
    fogler_ref = [240.0f0, 298.0f0]  # reference values in Fogler
    fogler_timespan = (0.0f0, 1.5f0)
    fixed_controlODE = ControlODE((u, p) -> fogler_ref, system!, u0, fogler_timespan)
    @info "Fogler's case: final time state" solve(fixed_controlODE, nothing).u[end]
    plot_simulation(
        fixed_controlODE,
        nothing;
        only=:states,
        vars=[1, 2, 3],
    )
    plot_simulation(
        fixed_controlODE,
        nothing;
        only=:states,
        vars=[4, 5],
    )

    function collocation(
        u0;
        num_supports::Integer=length(controlODE.tsteps),
        nodes_per_element::Integer=3,
        constrain_states::Bool=false,
    )
        optimizer = optimizer_with_attributes(
            Ipopt.Optimizer, "print_level" => 0, "check_derivatives_for_naninf" => "yes"
        )
        model = InfiniteModel(optimizer)
        method = OrthogonalCollocation(nodes_per_element)
        @infinite_parameter(
            model, t in [t0, tf], num_supports = num_supports, derivative_method = method
        )

        @variables(
            model,
            begin
                # state variables
                x[1:5], Infinite(t)
                # control variables
                c[1:2], Infinite(t)
            end
        )

        # tricks to make IPOPT work...
        almost_zero = 1.0f-3
        close_to_infinity = 1.0f3

        set_start_value_function(x[1], t -> 0.5f0)
        set_start_value_function(x[2], t -> 0.5f0)
        set_start_value_function(x[3], t -> 0.5f0)
        set_start_value_function(x[4], t -> 2.0f2)
        set_start_value_function(x[5], t -> 1.0f2)

        set_start_value_function(c[1], t -> 2.0f2)
        set_start_value_function(c[2], t -> 1.0f2)

        @constraints(
            model,
            begin
                almost_zero <= x[1] <= close_to_infinity
                almost_zero <= x[2] <= close_to_infinity
                almost_zero <= x[3] <= close_to_infinity
                almost_zero <= x[4] <= close_to_infinity
                almost_zero <= x[5] <= close_to_infinity
            end
        )

        # fixed_parameters
        CpA, CpB, CpC, CpH2SO4, N0H2S04, T0, CA0, HRA, HRB, E1A, E2B, A1, A2, UA, Tr1, Tr2 = values(
            system_params
        )

        initial_conditions = @constraints(
            model,
            begin
                x[1](0) == u0[1]
                x[2](0) == u0[2] + almost_zero
                x[3](0) == u0[3] + almost_zero
            end
        )

        control_constraints = @constraints(
            model,
            begin
                control_ranges[1][1] <= c[1] <= control_ranges[1][2]
                control_ranges[2][1] <= c[2] <= control_ranges[2][2]
            end
        )

        if constrain_states
            state_constraints = @constraints(
                model,
                begin
                    x[4] <= T_up
                    x[5] <= V_up
                end
            )
        end

        dynamic_constraints = @constraints(
            model,
            begin
                ∂(x[1], t) ==
                -(A1 * exp(E1A * ((1.0f0 / Tr1) - (1.0f0 / x[4]))) * x[1]) +
                (CA0 - x[1]) * c[1] / x[5]
                ∂(x[2], t) ==
                0.5f0 * (A1 * exp(E1A * ((1.0f0 / Tr1) - (1.0f0 / x[4]))) * x[1]) -
                (A2 * exp(E2B * ((1.0f0 / Tr2) - (1.0f0 / x[4]))) * x[2]) -
                x[2] * c[1] / x[5]
                ∂(x[3], t) ==
                3.0f0 * (A2 * exp(E2B * ((1.0f0 / Tr2) - (1.0f0 / x[4]))) * x[2]) -
                x[3] * c[1] / x[5]
                ∂(x[4], t) ==
                (
                    UA * (c[2] - x[4]) - CA0 * c[1] * CpA * (x[4] - T0) +
                    (
                        HRA * (-(A1 * exp(E1A * ((1.0f0 / Tr1) - (1.0f0 / x[4]))) * x[1])) +
                        HRB * (-(A2 * exp(E2B * ((1.0f0 / Tr2) - (1.0f0 / x[4]))) * x[2]))
                    ) * x[5]
                ) / ((x[1] * CpA + CpB * x[2] + CpC * x[3]) * x[5] + N0H2S04 * CpH2SO4)
                ∂(x[5], t) == c[1]
            end
        )

        @objective(model, Max, x[2](tf))

        optimize!(model)
        jump_model = optimizer_model(model)

        # list possible termination status: model |> termination_status |> typeof
        @info solution_summary(jump_model; verbose=false)
        if Int(termination_status(jump_model)) ∉ (1, 4)  # OPTIMAL = 1, LOCALLY_SOLVED = 4
            @warn raw_status(jump_model) termination_status(jump_model)
        end

        states = hcat(value.(x)...) |> permutedims
        controls = hcat(value.(c)...) |> permutedims

        return model, supports(t), states, controls
    end

    # set arquitecture of neural network controller
    controller = FastChain(
        (x, p) -> [1f2x[1], 1f2x[2], 1f2x[3], x[4], x[5]],
        FastDense(5, 8, tanh_fast),
        FastDense(8, 8, tanh_fast),
        FastDense(8, 2),
        # (x, p) -> [240f0, 298f0],
        scaled_sigmoids(control_ranges),
    )

    control_profile, infopt_model, times_collocation, states_collocation, controls_collocation = collocation_preconditioner(
        u0,
        collocation;
        plot=false,
        nodes_per_element=2,
        constrain_states=false,
    )
    _, _, _, states_collocation_constrained, controls_collocation_constrained = collocation_preconditioner(
        u0,
        collocation;
        plot=false,
        nodes_per_element=2,
        constrain_states=true,
    )
    #=
    plt.figure()
    plt.plot(times_collocation, states_collocation[1, :]; label="s1")
    plt.plot(times_collocation, states_collocation[2, :]; label="s2")
    plt.plot(times_collocation, states_collocation[3, :]; label="s3")
    plt.plot(
        times_collocation, states_collocation_constrained[1, :];
        label="s1_constrained", color=plt.gca().lines[1].get_color(), ls="dashdot"
    )
    plt.plot(
        times_collocation, states_collocation_constrained[2, :];
        label="s2_constrained", color=plt.gca().lines[2].get_color(), ls="dashdot"
    )
    plt.plot(
        times_collocation, states_collocation_constrained[3, :];
        label="s3_constrained", color=plt.gca().lines[3].get_color(), ls="dashdot"
    )
    plt.legend()
    plt.show()

    plt.figure()
    plt.plot(times_collocation, states_collocation[4, :]; label="s4")
    plt.plot(times_collocation, states_collocation[5, :]; label="s5")
    plt.plot(
        times_collocation, states_collocation_constrained[4, :];
        label="s4_constrained", color=plt.gca().lines[1].get_color(), ls="dashdot"
    )
    plt.plot(
        times_collocation, states_collocation_constrained[5, :];
        label="s5_constrained", color=plt.gca().lines[2].get_color(), ls="dashdot"
    )
    plt.legend()
    plt.show()

    plt.figure()
    plt.plot(times_collocation, controls_collocation[1, :]; label="c1")
    plt.plot(times_collocation, controls_collocation[2, :]; label="c2")
    plt.plot(
        times_collocation, controls_collocation_constrained[1, :];
        label="c1_constrained", color=plt.gca().lines[1].get_color(), ls="dashdot"
    )
    plt.plot(
        times_collocation, controls_collocation_constrained[2, :];
        label="c2_constrained", color=plt.gca().lines[2].get_color(), ls="dashdot"
    )
    plt.legend()
    plt.show()
    =#

    controlODE = ControlODE(controller, system!, u0, tspan; tsteps)

    θ = preconditioner(
        controlODE,
        control_profile;
        progressbar=true,
        plot_progress=false,
        # control_range_scaling=[range[end] - range[1] for range in control_ranges],
    )

    plot_simulation(controlODE, θ; only=:states)
    plot_simulation(controlODE, θ; only=:controls)
    store_simulation("precondition", controlODE, θ; datadir)

    # objective function splitted componenets to optimize
    function losses(controlODE, params; α=1.0f-3, δ=1.0f1, ρ=1.0f-2, kwargs...)

        # integrate ODE system
        sol_array = Array(solve(controlODE, params; kwargs...))

        # running cost
        out_temp = map(x -> relaxed_log_barrier(T_up - x; δ), sol_array[4, 1:end])
        out_vols = map(x -> relaxed_log_barrier(V_up - x; δ), sol_array[5, 1:end])

        # terminal cost
        # L = - (100 x₁ - x₂) + penalty  # Bradford
        objective = -sol_array[2, end]

        # integral penalty
        state_penalty = α * Δt * (sum(out_temp) + sum(out_vols))
        control_penalty = 0.0f0
        regularization = ρ * mean(abs2, θ)
        return objective, state_penalty, control_penalty, regularization
    end

    # α: penalty coefficient
    # δ: barrier relaxation coefficient
    α0, δ0 = 1.0f-5, 1.0f1
    barrier_iterations = 0:8
    αs = [α0 for _ in barrier_iterations]
    δs = [δ0 * 0.8f0^i for i in barrier_iterations]
    θ = constrained_training(
        controlODE,
        losses;
        starting_params=θ,
        αs,
        δs,
        show_progressbar=true,
        # plots_callback,
        datadir,
    )

    @info "Final states"
    plot_simulation(controlODE, θ; only=:states, vars=[1, 2, 3])
    plot_simulation(
        controlODE, θ; only=:states, vars=[4, 5], yrefs=[T_up, V_up]
    )

    @info "Final controls"
    plot_simulation(controlODE, θ; only=:controls)#  only=:states, vars=[1,2,3])

    @info "Final loss" losses(controlODE, θ; δ=δs[end])

    final_objective, final_state_penalty, final_control_penalty, final_regularization = losses(
        controlODE, θ; δ=δs[end]
    )

    return store_simulation(
        "constrained",
        controlODE,
        θ;
        datadir,
        metadata=Dict(
            :objective => final_objective,
            :state_penalty => final_state_penalty,
            :cotrol_penalty => final_control_penalty,
            :regularization => final_regularization,
            :num_params => length(θ),
            :layers => controller_shape(controller),
            :deltas => δs,
            :t0 => t0,
            :tf => tf,
            :Δt => Δt,
        ),
    )
end  # function wrapper


text:
Base.print(io::IO, bp::Bipartition) = print(io, mapreduce(x -> x ? "●" : "○", *, bp.v))

function Base.:(==)(a::Bipartition, b::Bipartition)
	cv1, cv2 = .! a.v, .! b.v

	return (a.v == b.v) || (a.v == cv2) || (cv1 == b.v) || (cv1 == cv2)
end

Base.show(io::IO, bp::Bipartition) = print(io, string(bp))

function Base.show(io::IO, ::MIME"text/plain", bp::Bipartition)
	summary(io, bp)
	print(io, "  ", string(bp))

	return nothing
end

function update_bipartition!(p::Node, q::Node, n::Int)
	br = getbranch(p, q)
	v = falses(n)
	if q.species > 0
		v[q.species] = true
	end
	for link in q.links
		link.to == q && continue
		v .|= link.branch.bipart.v
	end
	br.bipart = Bipartition(v)

	return nothing
end

function compute_bipartitions!(tree::AbstractTree, dir::SpeciesDirectory)
	n = length(dir)
	for (p, q) in preorder(tree)
		p == q && continue
		br = getbranch(p, q)
		v = falses(n)
		v[filter(x -> x > 0, getspecies(p, q))] .= true
		v = v[1] ? .!(v) : v
		br.bipart = Bipartition(v)
	end

	return nothing
end

"""
	compute_bipartitions!(tree)

Compute the bipartitions of species (in internal or external nodes) in a tree. Use the `bipartitions` function to retreive the bipartitions after computing them.
"""
function compute_bipartitions!(tree::AbstractTree)

	compute_bipartitions!(tree, tree.dir) # this should work

	return nothing
end

function _clear_bipartitions!(p::AbstractNode, q::AbstractNode, nullbp::Bipartition)
	for link in q.links
		link.to == p && continue
		link.branch.bipart = nullbp
		_clear_bipartitions!(q, link.to, nullbp)
	end

	return nothing
end

"""
	clear_bipartitions!(tree)

Remove the `Bipartition` objects attached to the branches of a `tree`.
"""
function clear_bipartitions!(tree)
	nullbp = Bipartition([false])
	for link in tree.origin.links
		_clear_bipartitions!(tree.origin, link.to, nullbp)
	end

	return nothing
end

"""
	istrivial(`bipartition`) -> Bool

A `bipartition` is trivial if it separates just one species from the rest.
"""
function istrivial(bp::Bipartition)::Bool
	i = length(bp.v)
	s = sum(bp.v)
	return (s == i - 1 || s == 1)
end

"""
	isinformative(bipartition) -> Bool

A `bipartition` is informative if it separates at least one species from the rest.

Uninformative bipartitions are not even bipartitions proper, as they don't define disjoint sets of species.
"""
isinformative(bp::Bipartition)::Bool = sum(bp.v) > 0

"""
	are_compatible(bipartition1, bipartition2) -> Bool

Return true if two bipartitions can be present in the same tree.
"""
function are_compatible(bp1::Bipartition, bp2::Bipartition)
	cv1 = .! bp1.v
	cv2 = .! bp2.v

	(bp1.v == bp1.v .| bp2.v) && return true
	(bp1.v == bp1.v .| cv2) && return true
	(cv1 == cv1 .| bp2.v) && return true
	(cv1 == cv1 .| cv2) && return true

	return false
end

"""
	are_conflicting(bipartition1, bipartition2) -> Bool

Return true two bipartitions cannot be present in the same tree.
"""
are_conflicting(bp1::Bipartition, bp2::Bipartition) = ! are_compatible(bp1, bp2)


text:
module URLUtils

# from: https://github.com/JuliaWeb/HttpCommon.jl/blob/release-0.3/src/HttpCommon.jl

# All characters that remain unencoded in URI encoding
#                                   ( AKA URL encoding
#                                     AKA percent-encoding )
#
const URIwhitelist = Set(Any['A','B','C','D','E','F','G','H','I',
                         'J','K','L','M','N','O','P','Q','R',
                         'S','T','U','V','W','X','Y','Z',
                         'a','b','c','d','e','f','g','h','i',
                         'j','k','l','m','n','o','p','q','r',
                         's','t','u','v','w','x','y','z',
                         '0','1','2','3','4','5','6','7','8',
                         '9','-','_','.','~'])

# decodeURI
#
# Decode URI encoded strings
#
function decodeURI(encoded::String)
    enc = split(replace(encoded,"+"," "),"%")
    decoded = enc[1]
    for c in enc[2:end]
        decoded = @compat string(decoded, Char((parse(Int, c[1:2], 16))),
                                 c[3:end])
    end
    decoded
end

# encodeURI
#
# Convert strings to URI encoding
#
function encodeURI(decoded::String)
    encoded = ""
    for c in decoded
        encoded = @compat encoded * string(c in URIwhitelist ? c :
                                           "%" * uppercase(hex(Int(c))))
    end
    encoded
end

# parsequerystring
#
# Convert a valid querystring to a Dict:
#
#    q = "foo=bar&baz=%3Ca%20href%3D%27http%3A%2F%2Fwww.hackershool.com%27%3Ehello%20world%21%3C%2Fa%3E"
#    parsequerystring(q)
#    # => Dict("foo"=>"bar","baz"=>"hello world!")
#
function parsequerystring(query::String)
    q = Dict{String,String}()
    if !('=' in query)
        return throw("Not a valid query string: $query, must contain at least one key=value pair.")
    end
    for set in split(query, "&")
        key, val = split(set, "=")
        q[decodeURI(key)] = decodeURI(val)
    end
    q
end

end # module


text:
using PowerSystems

################## OMIB data #####################

nodes_OMIB= [Bus(1 , #number
                 "Bus 1", #Name
                 "REF" , #BusType (REF, PV, PQ)
                 0, #Angle in radians
                 1.06, #Voltage in pu
                 (min=0.94, max=1.06), #Voltage limits in pu
                 69), #Base voltage in kV
             Bus(2 , "Bus 2"  , "PV" , 0 , 1.045 , (min=0.94, max=1.06), 69)]

branch_OMIB = [Line("Line1", #name
                    true, #available
                    0.0, #active power flow initial condition (from-to)
                    0.0, #reactive power flow initial condition (from-to)
                    Arc(from=nodes_OMIB[1], to=nodes_OMIB[2]), #Connection between buses
                    0.01938, #resistance in pu
                    0.05917, #reactance in pu
                    (from=0.0264, to=0.0264), #susceptance in pu
                    18.046, #rate in MW
                    1.04)]  #angle limits (-min and max)

Dynbranch_OMIB = [DynLine("Line1", #name
                    true, #available
                    Arc(from=nodes_OMIB[1], to=nodes_OMIB[2]), #Connection between buses
                    0.01938, #resistance in pu
                    0.05917, #reactance in pu
                    (from=0.0264, to=0.0264))]  #susceptance in pu

loads = [PowerLoad("Bus2", #name
                   true, #available
                   nodes_OMIB[1], #bus
                   PowerSystems.ConstantPower, #model
                   0.217, #active power in pu
                   0.127, #reactive power in pu
                   0.217, #max active power in pu
                   0.127)] #max reactive power in pu

################## Kundur machine case data #####################

nodes_kundur= [Bus(1 , "Bus 1"  , "PQ" , 0 , 1.02  , (min=0.94, max=1.06), 24),
                    Bus(2 , "Bus 2"  , "REF" , 0 , 1.01 , (min=0.94, max=1.06), 24)]

branch_kundur = [Line("Line1", true, 0.0, 0.0, Arc(from=nodes_kundur[1], to=nodes_kundur[2]),
                      0.01938, 0.05917, (from=0.0, to=0.0), 18.046, 1.04)]

loads_kundur = [PowerLoad("Bus1", true, nodes_kundur[1], PowerSystems.ConstantPower, 0.0, 0.0, 0.0, 0.0)]

################ Benchmark 3 Bus case ########################

nodes_benchmark = [Bus(1 , "Bus 1"  , "REF" , 0 , 1.02  , (min=0.94, max=1.06), 138),
                    Bus(2 , "Bus 2"  , "PV" , 0 , 1.00 , (min=0.94, max=1.06), 138),
                    Bus(3 , "Bus 3"  , "PV" , 0 , 1.00 , (min=0.94, max=1.06), 138)]

branch_benchmark = [Line("Line1", true, 0.0, 0.0, Arc(from=nodes_benchmark[1], to=nodes_benchmark[3]), 0.01, 0.12, (from=0.0, to=0.0), 100, 1.04),
                    Line("Line2", true, 0.0, 0.0, Arc(from=nodes_benchmark[1], to=nodes_benchmark[2]), 0.01, 0.12, (from=0.0, to=0.0), 100, 1.04),
                    Line("Line3", true, 0.0, 0.0, Arc(from=nodes_benchmark[2], to=nodes_benchmark[3]), 0.01, 0.12, (from=0.0, to=0.0), 100, 1.04)]

branch_benchmark_noline = [Line("Line2", true, 0.0, 0.0, Arc(from=nodes_benchmark[1], to=nodes_benchmark[2]), 0.01, 0.12, (from=0.0, to=0.0), 100, 1.04),
                           Line("Line3", true, 0.0, 0.0, Arc(from=nodes_benchmark[2], to=nodes_benchmark[3]), 0.01, 0.12, (from=0.0, to=0.0), 100, 1.04)]

loads_benchmark = [PowerLoad("Bus1", true, nodes_benchmark[1], PowerSystems.ConstantPower, 1.5, 0.8, 1.5, 0.8),
                   PowerLoad("Bus2", true, nodes_benchmark[2], PowerSystems.ConstantPower, 1.5, 0.7, 1.5, 0.8),
                   PowerLoad("Bus3", true, nodes_benchmark[3], PowerSystems.ConstantPower, 0.5, 0.3, 0.5, 0.3)]

################ Benchmark 2 Bus case ########################

nodes_bench = [Bus(1 , "Bus 1"  , "REF" , 0 , 1.02  , (min=0.94, max=1.06), 138),
                    Bus(2 , "Bus 2"  , "PV" , 0 , 1.00 , (min=0.94, max=1.06), 138)]

branch_bench = [Line("Line1", true, 0.0, 0.0, Arc(from=nodes_bench[1], to=nodes_bench[2]), 0.01, 0.12, (from=0.0, to=0.0), 100, 1.04)]

loads_bench = [PowerLoad("Bus1", true, nodes_bench[2], PowerSystems.ConstantPower, 1.0, 0.3, 1.5, 0.8)]

############## Static Kundur Machines for PF #######################

machine_Kundur = ThermalStandard("Bus1", true, nodes_kundur[1], 2.22, 0,
                    TechThermal(5.55, PowerSystems.ST, PowerSystems.COAL, (min=0.0, max=5.55), (min=-1.0, max=1.0), nothing, nothing),
                    ThreePartCost((430.292599, 2000.0), 0.0, 0.0, 0.0)
                    )
infinite_bus = ThermalStandard("Bus2", true, nodes_kundur[2], -2.22, 0,
                    TechThermal(5.55, PowerSystems.ST, PowerSystems.COAL, (min=-500.55, max=500.55), (min=-5.55, max=5.55), nothing, nothing),
                    ThreePartCost((430.292599, 2000.0), 0.0, 0.0, 0.0)
                    )


text:
using ObjectDetectors
using GeometryBasics

@testset "PolygonDetection" begin
    vertices = [Point2(2.0, 2.0), Point2(4.0, 2.0), Point2(3.0, 3.0)]
    core_args = ("class", 0.7, 256, 256, MACHINE, "detector")
    detection = PolygonDetection(vertices, core_args...)

    @test centroid(detection) == [3, 7 / 3]
    bb = bounding_box(detection)
    @test bb.origin == [2.0, 2.0]
    @test bb.widths == [2.0, 1.0]

    @test_throws ArgumentError PolygonDetection([Point2(2.0, 2.0), Point2(4.0, 2.0)], core_args...)
end


text:
using Cairo
using LSystems

c = CairoRGBSurface(485, 720);
cr = CairoContext(c);

save(cr);

set_source_rgb(cr, 1.0, 1.0, 1.0);
rectangle(cr, 0.0, 0.0, 485.0, 720.0);
fill(cr);

restore(cr);
save(cr);

x, y, a, d = 310., 560., 90, 15
set_line_width(cr, 1);
move_to(cr, x, y);

dragon_start = "A"
dragon_trans = Dict([
('A', "A+BF"),
('B', "FA-B")
])

function determine_new_position(x, y, a, d)
    x -= sind(a) * d
    y -= cosd(a) * d
    return x, y
end

level = 10
drawing_instructions = @task dol_system(dragon_start, dragon_trans, level)

for v in drawing_instructions
    if v == 'F'
        x, y = determine_new_position(x, y, a, d)
        line_to(cr, x, y)
    elseif v == '+'
        a -= 90
    elseif v == '-'
        a += 90
    elseif v == '['
        push!(s, [x, y, a])
    elseif v == ']'
        x, y, a = pop!(s)
        move_to(cr, x, y)
    end
end

stroke(cr);
write_to_png(c, "dragon_$level.png");