text:
"""
Utility for performing LU factorizations with UMFPACK, Pardiso, or MUMPS
"""
module LU_Factorizations

export LUfact
export MSolver
export PSolver
export USolver
export DEFAULT_LUPACK
export AbstractLUPACK

using ArnoldiMethodTransformations
using ArnoldiMethodTransformations: AbstractSolver
using LinearAlgebra
# using MPI
# using MUMPS
using Pardiso
using SparseArrays

import ..LUPACK
@eval const DEFAULT_LUPACK = $(LUPACK)()
"""
    DEFAULT_LUPACK = $DEFAULT_LUPACK
"""
DEFAULT_LUPACK

"""
    MSolver
"""
MSolver

"""
    PSolver
"""
PSolver

"""
    USolver
"""
USolver

"""
    AbstractLUPACK
"""
AbstractLUPACK = AbstractSolver

# colors for pretty printing
import ..PRINTED_COLOR_NUMBER
import ..PRINTED_COLOR_LIGHT

# main structure used internally
struct LUfact{M,TLU,T}
    LU::TLU
    Z::SparseMatrixCSC{T,Int}
end

LinearAlgebra.lu(A,::T,args...;kwargs...) where T<:AbstractSolver  = lu(A,T,args...;kwargs...)
function LinearAlgebra.lu(A::AbstractArray{T},::Type{USolver},args...;kwargs...) where T
    alu = lu(A,args...;kwargs...)
    return LUfact{USolver,typeof(alu),T}(alu,spzeros(T,size(A)...))
end
function LinearAlgebra.lu(A::AbstractArray{T},::Type{MSolver},args...;kwargs...) where T
    throw("MUMPS still under construction")
    try MPI catch; throw(ErrorException("MUMPS not loaded. Try again after `using MUMPS`")) end
    MPI.Initialized() ? nothing : MPI.Init()
    alu = mumps_factorize(α)
    return LUfact{MSolver,typeof(alu),T}(alu,spzeros(T,size(A)...))
end
function LinearAlgebra.lu(A::AbstractArray{T},::Type{PSolver},args...;kwargs...) where T
    try PardisoSolver catch; throw(ErrorException("Pardiso not loaded. Try again after `using Pardiso`")) end
    issym = issymmetric(A)
    type = eltype(A)
    alu = PardisoSolver()
    set_iparm!(alu,1,1) # don't revert to defaults
    set_iparm!(alu,12,1) # transpose b/c of CSR vs SCS
    x = Vector{type}(undef,size(A,1))
    y = Vector{type}(undef,size(A,1))
    if issym & (type<:Real)
        set_matrixtype!(alu,Pardiso.REAL_SYM)
        pardiso(alu,x,sparse(tril(A)),y)
    elseif issym & (type<:Complex)
        set_matrixtype!(alu,Pardiso.COMPLEX_SYM)
        pardiso(alu,x,sparse(tril(A)),y)
    elseif !issym & (type<:Real)
        set_matrixtype!(alu,Pardiso.REAL_NONSYM)
        pardiso(alu,x,sparse(A),y)
    else # if !issym & type<:Complex
        set_matrixtype!(alu,Pardiso.COMPLEX_NONSYM)
        pardiso(alu,x,sparse(A),y)
    end
    set_phase!(alu,12) # analyze and factorize
    set_phase!(alu,33) # set to solve for future calls
    return LUfact{PSolver,typeof(alu),T}(alu,spzeros(T,size(A)...))
end

# define action of :LUfact
LinearAlgebra.ldiv!(Y,A::LUfact{USolver},B) = ldiv!(Y, A.LU, B)
LinearAlgebra.ldiv!(A::LUfact{USolver},B) = ldiv!(A.LU, B)
LinearAlgebra.ldiv!(Y,A::LUfact{MSolver},B) = ldiv!(Y, A.LU, B)
LinearAlgebra.ldiv!(Y,A::LUfact{PSolver},B) = pardiso(A.LU,Y,A.Z,B)
function Base.:\(A::LUfact,B)
    Y = similar(B)
    ldiv!(Y,A,B)
    return Y
end

# pretty printing
function Base.show(io::IO,si::LUfact{M}) where M
    printstyled(io,"LUfact ",color=PRINTED_COLOR_LIGHT)
    print(io,"(using ")
    if M<:MSolver
        printstyled(io,"MUMPS",color=PRINTED_COLOR_NUMBER)
    elseif M<:PSolver
        printstyled(io,"Pardiso",color=PRINTED_COLOR_NUMBER)
    elseif M<:USolver
        printstyled(io,"UMFPACK",color=PRINTED_COLOR_NUMBER)
    end
    print(io,")")
end

end # module


text:
function get_model(model_type)
    s = Symbol(model_type)
    return getfield(_PM, s)
end

function get_solver(optimizer::String, nl::String="ipopt", mip::String="cbc",
    log_level::Int=0, time_limit::Float64=Inf, nl_time_limit::Float64=Inf,
    mip_time_limit::Float64=Inf, ipopt_tol::Float64=1e-8)

    if optimizer == "gurobi"
            solver = JuMP.optimizer_with_attributes(Gurobi.Optimizer, "TimeLimit" => time_limit, "OutputFlag" => log_level)
    end

    if optimizer == "ipopt"
            ipopt_tol = 1e-8
            solver = JuMP.optimizer_with_attributes(Ipopt.Optimizer, "print_level" => log_level, "max_cpu_time" => time_limit,
                "tol" => ipopt_tol)
    end

    if optimizer == "juniper" && nl == "ipopt" && mip == "cbc"
        mip_solver = JuMP.optimizer_with_attributes(Cbc.Optimizer, "logLevel" => log_level, "seconds" => mip_time_limit)
        nl_solver = JuMP.optimizer_with_attributes(Ipopt.Optimizer, "print_level" => log_level, "max_cpu_time" => nl_time_limit)
        solver = JuMP.optimizer_with_attributes(Juniper.Optimizer,
                     "nl_solver" => nl_solver,
                     "mip_solver" => mip_solver,
                     "log_levels" => [],
                     "time_limit" => time_limit)
    end

    if optimizer == "juniper" && nl == "gurobi" && mip == "cbc"
        mip_solver = JuMP.optimizer_with_attributes(Cbc.Optimizer, "logLevel" => log_level, "seconds" => mip_time_limit)
        nl_solver = JuMP.optimizer_with_attributes(Gurobi.Optimizer, "TimeLimit" => nl_time_limit)
        solver = JuMP.optimizer_with_attributes(Juniper.Optimizer,
                     "nl_solver" => nl_solver,
                     "mip_solver" => mip_solver,
                     "log_levels" => [],
                     "time_limit" => time_limit)
    end

    if optimizer == "juniper" && nl == "gurobi" && mip == "gurobi"
        mip_solver = JuMP.optimizer_with_attributes(Gurobi.Optimizer, "TimeLimit" => mip_time_limit)
        nl_solver = JuMP.optimizer_with_attributes(Gurobi.Optimizer, "TimeLimit" => nl_time_limit)
        solver = JuMP.optimizer_with_attributes(Juniper.Optimizer,
                     "nl_solver" => nl_solver,
                     "mip_solver" => mip_solver,
                     "log_levels" => [],
                     "time_limit" => time_limit)
    end

    if optimizer == "knitro"
        solver = JuMP.optimizer_with_attributes(KNITRO.Optimizer)
    end

    if optimizer == "cbc"
        solver = JuMP.optimizer_with_attributes(Cbc.Optimizer, "seconds" => time_limit)
    end

    if optimizer == "scip"
        solver = JuMP.optimizer_with_attributes(SCIP.Optimizer)
    end

    return solver

end

function load_pm_from_json(json_path)
    pm = Dict()
    open(json_path, "r") do f
        pm = JSON.parse(f)
    end
    for (idx, gen) in pm["gen"]
        if gen["model"] == 1
            pm["gen"][idx]["cost"] = convert(Array{Float64,1}, gen["cost"])
        end
    end
    if pm["correct_pm_network_data"]
        _PM.correct_network_data!(pm)
    end
    return pm
end

ext: jl


text:
# EXCLUDE FROM TESTING
using VectorEngine

function kernel()
    @veprintf("Hello world from VE!\n")
    return nothing
end

# Compile the kernel
vefunc = VectorEngine.vefunction(kernel)

# Launch kernel
vefunc()

synchronize()


text:
"""
    CapsuleCapsuleCollision

    collision between two capsules

    origin_parent: radius of contact for parent body
    origin_child: radius of contact for child body
    orientation_parent: radius of contact for parent body
    orientation_child: radius of contact for child body
    radius_parent: radius of contact for parent body
    radius_child: radius of contact for child body
    height_parent: radius of contact for parent body
    height_child: radius of contact for child body
"""
mutable struct CapsuleCapsuleCollision{T,O,I,OI} <: Collision{T,O,I,OI}
    origin_parent::SVector{I,T}
    origin_child::SVector{I,T}
    orientation_parent::Quaternion{T}
    orientation_child::Quaternion{T}
    radius_parent::T
    radius_child::T
    height_parent::T
    height_child::T
end

# distance
function distance(collision::CapsuleCapsuleCollision, xp, qp, xc, qc)
    # contact points
    pp = contact_point(:parent, collision, xp, qp, xc, qc)
    pc = contact_point(:child,  collision, xp, qp, xc, qc)

    # difference
    d = pp - pc

    return norm(pp - pc, 2)
end

function ∂distance∂x(gradient::Symbol, collision::CapsuleCapsuleCollision, xp, qp, xc, qc)
    # contact points
    pp = contact_point(:parent, collision, xp, qp, xc, qc)
    pc = contact_point(:child,  collision, xp, qp, xc, qc)

    # difference
    d = pp - pc

    X = ∂norm∂x(d) * (∂contact_point∂x(:parent, jacobian, collision, xp, qp, xc, qc) - ∂contact_point∂x(:child, jacobian, collision, xp, qp, xc, qc))

    return X
end

function ∂distance∂q(gradient::Symbol, collision::CapsuleCapsuleCollision, xp, qp, xc, qc)
    # contact points
    pp = contact_point(:parent, collision, xp, qp, xc, qc)
    pc = contact_point(:child,  collision, xp, qp, xc, qc)

    # difference
    d = pp - pc

    X = ∂norm∂x(d) * (∂contact_point∂q(:parent, jacobian, collision, xp, qp, xc, qc) - ∂contact_point∂q(:child, jacobian, collision, xp, qp, xc, qc))

    return X
end

# contact point in world frame
function contact_point(relative::Symbol, collision::CapsuleCapsuleCollision, xp, qp, xc, qc)
    # call mehrotra

    if relative == :parent
        return collision.ip.z[0 .+ (1:3)]
    elseif relative == :child
        return collision.ip.z[3 .+ (1:3)]
    end
end

function ∂contact_point∂x(relative::Symbol, jacobian::Symbol, collision::CapsuleCapsuleCollision, xp, qp, xc, qc)
    # call mehrotra

    if relative == :parent
        if jacobian == :parent
            X = collision.ip.δz[]
        elseif jacobian == :child
            X = collision.ip.δz[]
        end
    elseif relative == :child
        if jacobian == :parent
            X = collision.ip.δz[]
        elseif jacobian == :child
            X = collision.ip.δz[]
        end
    end

    return X
end

function ∂contact_point∂q(relative::Symbol, jacobian::Symbol, collision::CapsuleCapsuleCollision, xp, qp, xc, qc)
# call mehrotra

    if relative == :parent
        if jacobian == :parent
            X = collision.ip.δz[]
        elseif jacobian == :child
            X = collision.ip.δz[]
        end
    elseif relative == :child
        if jacobian == :parent
            X = collision.ip.δz[]
        elseif jacobian == :child
            X = collision.ip.δz[]
        end
    end

    return X
end


text:
using Sockets

pipename_for_wrapper = ARGS[2]
pipename_for_debugger = ARGS[1]
cwd = ARGS[3]
julia_env = ARGS[4]
pipename_for_crashreporting = ARGS[5]

conn = Sockets.connect(pipename_for_wrapper)

jl_cmd = joinpath(Sys.BINDIR, Base.julia_exename())

debugger_script = joinpath(@__DIR__, "run_debugger.jl")

cmd = Cmd(`$jl_cmd --color=yes --history-file=no --startup-file=no --project=$julia_env $debugger_script $pipename_for_debugger $pipename_for_crashreporting`, dir = cwd)

p = run(pipeline(cmd, stdin = stdin, stdout = stdout, stderr = stderr), wait = false)

@async begin
    l = readline(conn)

    if l == "TERMINATE"
        kill(p)
    else
        error("Invalid state.")
    end
end

wait(p)

println()
printstyled("Julia debuggee finished. Press ENTER to close this terminal.\n", bold = true)

readline()

max_stars_repo_name: cloudgiant/julia-vscode

max_issues_repo_name: cloudgiant/julia-vscode

max_forks_repo_name: cloudgiant/julia-vscode


text:
module EMLflash_cost
	using DanaTypes
	using DotPlusInheritance
	using Reexport
	@reexport using ...stage_separators.EMLflash
	import EMLtypes.length
	include("flash_cost/flash_cost.jl")
end


text:
# # Optimization of the default MARGO configuration

# ## Using `ClimateMARGO.jl`
using ClimateMARGO
using ClimateMARGO.Models
using ClimateMARGO.Optimization

# ## Loading preset configurations

# Load the pre-defined default MARGO parameters, which are described by the ClimateModelParameters struct
params = deepcopy(ClimateMARGO.IO.included_configurations["default"])

# Create a MARGO instance based on these parameters
m = ClimateModel(params);

# By default, the climate control timeseries are all set to zero.
m.controls

# ## Real-time climate optimization

# Let's optimize the controls with the default settings, which finds the cheapest combination of the control timeseries that keeps adapted temperatures below 2°C.
@time optimize_controls!(m);

# The optimization can be slow the first time since it has to compile.
# Let's re-initialize the model and see how fast it runs now that the the optimization routine has been precompiled.
m = ClimateModel(params);
@time optimize_controls!(m);

# ## Visualizing the results

# Finally, let's plot the resulting temperature timeseries
using PyPlot
ClimateMARGO.Plotting.plot_state(m);
gcf()


text:
#Adopted from RDatasets

module Datasets

using Reexport
@reexport using DataFrames
using RLESUtils
@reexport using DataFrameSets

export dataset, list_datasets, list_dataset, datadir, load_meta, write_dataset

const DATAPATH = joinpath(dirname(@__FILE__), "..", "data")
const METAFILE = "_META.csv.gz"

#Load a particular dataframe
function dataset(package_name::AbstractString, dataset_name::AbstractString)
    dirpath = joinpath(DATAPATH, package_name)

    filename = joinpath(dirpath, string(dataset_name, ".csv.gz"))
    if !isfile(filename)
        error(@sprintf "Unable to locate file %s\n" filename)
    else
        return readtable(filename)
    end
end

#Load all dataframes in the package, returns a DFSet
function dataset(data_name::AbstractString)
    dirpath  = joinpath(DATAPATH, data_name)
    if !isdir(dirpath)
        error(@sprintf "No such directory %s\n" dirpath)
    end
    Ds = load_csvs(dirpath)
    Ds
end

#Load all dataframes in the package, returns a DFSet
function dataset(data_name::AbstractString, label::Symbol;
    transform::Function=identity)
    Ds = dataset(data_name)
    Dl = DFSetLabeled(Ds, label; transform=transform)
    Dl
end

function list_datasets()
    dsets = readdir(DATAPATH)
    filter!(x -> isdir(joinpath(DATAPATH, x)), dsets)
    dsets
end

function list_dataset(package_name::AbstractString)
   data = readdir(joinpath(DATAPATH, package_name))
   filter!(x -> endswith(x, ".csv.gz"), data)
   data
end

datadir(package_name::AbstractString="") = abspath(joinpath(DATAPATH, package_name))

function load_meta(data_name::AbstractString)
    dirpath  = joinpath(DATAPATH, data_name)
    if !isdir(dirpath)
        error(@sprintf "No such directory %s\n" dirpath)
    end
    M = readtable(joinpath(dirpath, METAFILE))
    M
end

function write_dataset(data_name::AbstractString, Ds::DFSet)
    dirpath  = joinpath(DATAPATH, data_name)
    m = getmeta(Ds)
    m[:id] = 1:nrow(m)
    save_csvs(dirpath, Ds)
end

end # module


text:
# Test naive implementation of entropy calculations towards efficient calculation of entropy on a manifold

using DocStringExtensions
using KernelDensityEstimate

include(joinpath(dirname(@__FILE__), "circularEntropyUtils.jl"))

using Gadfly, Colors
using Distributions
using Random

using TransformUtils
using Optim

const TU = TransformUtils

using ApproxManifoldProducts

# logmap_SO2(Rl::Matrix{Float64}) = sign(Rl[2,1])*acos(Rl[1,1])
# difftheta(wth1::Float64, wth2::Float64)::Float64 = logmap_SO2(TU.R(wth1)'*TU.R(wth2))
# addtheta(wth1::Float64, wth2::Float64) = TU.wrapRad( wth1+wth2 )

# some test points to work with
pts = 0.3*randn(30)
pts = [pts; TU.wrapRad.(0.6*randn(70) .- pi)]
shuffle!(pts)

# p = kde!(pts)

# N = length(pts)
CV  = zeros(100)
global i_global = 0
BW= range(0.05, 1.0, length=100)
for bw in BW
    global i_global
    i_global += 1
    CV[i_global] = manifoldLooCrossValidation(pts, bw, own=true, diffop=difftheta)
end

plot(x=BW, y=CV, Geom.line)

BW[findfirst(CV .== maximum(CV))]

# getBW(kde!(pts))[1,1]

# what we want
# pc = kde!(pts, bw, mani=S1)

# evaluate(pc)

## do the same with

PD = zeros(100)

global i_global = 0
for bw in BW
    global i_global
    i_global += 1
    pd = kde!(pts, [bw;])
    PD[i_global] = KernelDensityEstimate.entropy(pd) # TODO addop, diffop
end

## optimize the KDE bandwidth

# initial testing values
lower = 0.001
upper = 10.0

# pts = 10.0.+randn(100)

minEntropyLOOCV = (bw) -> -manifoldLooCrossValidation(pts, bw, own=true, diffop=difftheta)

# TODO Compare Optim.GoldenSection vs KDE.golden
@time res = optimize(minEntropyLOOCV, lower, upper, GoldenSection(), x_tol=0.001)

bw = res.minimizer

pc = kde!(pts, [bw])

plotKDECircular(pc)

# gg = (x)->pc([x;])[1]
# arr = [gg;]
# pl = plotCircBeliefs(arr)

## Construct circular KDE

pc2 = kde!_CircularNaiveCV(pts)

getBW(pc2)[:,1]

plotKDECircular(pc2)

#


text:
#=

  Set covering deployment in Julia ConstraintSolver.jl

  From http://mathworld.wolfram.com/SetCoveringDeployment.html
  """
  Set covering deployment (sometimes written "set-covering deployment"
  and abbreviated SCDP for "set covering deployment problem") seeks
  an optimal stationing of troops in a set of regions so that a
  relatively small number of troop units can control a large
  geographic region. ReVelle and Rosing (2000) first described
  this in a study of Emperor Constantine the Great's mobile field
  army placements to secure the Roman Empire.
  """

  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#
using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
const CS = ConstraintSolver
include("constraints_utils.jl")

function set_covering_deployment(problem,print_solutions=true,all_solutions=true)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   # "all_solutions"=> all_solutions,
                                                            "all_optimal_solutions"=>true,
                                                            "logging"=>[],

                                                            "traverse_strategy"=>:BFS,
                                                            # "traverse_strategy"=>:DFS,
                                                            # "traverse_strategy"=>:DBFS,

                                                            # "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            "branch_split"=>:InHalf,

                                                            # https://wikunia.github.io/ConstraintSolver.jl/stable/options/#branch_strategy-(:Auto)
                                                            # "branch_strategy" => :IMPS, # default
                                                            "branch_strategy" => :ABS, # Activity Based Search
                                                            "activity.decay" => 0.999, # default 0.999
                                                            "activity.max_probes" => 100, # default, 10
                                                            "activity.max_confidence_deviation" => 20, # default 20

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>6,

                                                            # "backtrack" => false, # default true
                                                            # "backtrack_sorting" => false, # default true

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            # "lp_optimizer" => glpk_optimizer,
                                                            # "lp_optimizer" => ipopt_optimizer,
                                        ))

    n, _ = size(problem)
    armies = ["Alexandria","Asia Minor","Britain","Byzantium","Gaul","Iberia","Rome","Tunis"]

    # First army
    @variable(model, xs[1:n], Bin)
    # Second army
    @variable(model, ys[1:n], Bin)

    # Objective
    @variable(model, 0 <= z <= 2*n, Int)

    # Constraint 1: There is always an army in a city (+ maybe a backup)
    #          Or rather: Is there a backup, there must be an
    #          an army
    for (x,y) in zip(xs,ys)
        @constraint(model,x>=y)
    end

    # Constraint 2: There should always be an backup army near
    # every city
    for i in 1:n
        @constraint(model, xs[i] + sum([y*m for (y,m) in zip(ys,problem[i,:]) ] ) >= 1)
    end

    # objective: minimize the number of armies
    @constraint(model, z == sum(xs + ys))

    @objective(model,Min,z)

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                println("solution #$sol")
                xs_val = convert.(Integer,JuMP.value.(xs; result=sol))
                ys_val = convert.(Integer,JuMP.value.(ys; result=sol))
                z_val = convert.(Integer,JuMP.value.(z; result=sol))
                println("z:$z_val")
                println("xs: $xs_val")
                println("ys: $ys_val")
                println("1st army: $([armies[i] for i in 1:n if xs_val[i] > 0])")
                println("2st army: $([armies[i] for i in 1:n if ys_val[i] > 0])")

                println()
            end
        end
    else
        println("status:$status")
    end

    return status
end

problem = resize_matrix(
   [[0, 1, 0, 1, 0, 0, 1, 1],
    [1, 0, 0, 1, 0, 0, 0, 0],
    [0, 0, 0, 0, 1, 1, 0, 0],
    [1, 1, 0, 0, 0, 0, 1, 0],
    [0, 0, 1, 0, 0, 1, 1, 0],
    [0, 0, 1, 0, 1, 0, 1, 1],
    [1, 0, 0, 1, 1, 1, 0, 1],
    [1, 0, 0, 0, 0, 1, 1, 0]]
)

@time set_covering_deployment(problem)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
"
    expands the given mps using the algorithm given in the vumps paper
"
@with_kw struct OptimalExpand<:Algorithm
    trscheme::TruncationScheme = truncdim(1)
end

function changebonds(state::InfiniteMPS, H::Hamiltonian,alg::OptimalExpand,envs=environments(state,H))
    #determine optimal expansion spaces around bond i
    pexp = PeriodicArray(map(1:length(state)) do i
        ACAR = _permute_front(state.AC[i])*_permute_tail(state.AR[i+1])
        AC2 = ac2_prime(ACAR,i,state,envs)

        #Calculate nullspaces for AL and AR
        NL = leftnull(state.AL[i])
        NR = rightnull(_permute_tail(state.AR[i+1]))

        #Use this nullspaces and SVD decomposition to determine the optimal expansion space
        intermediate = adjoint(NL)*AC2*adjoint(NR)
        (U,S,V) = tsvd(intermediate,trunc=alg.trscheme,alg=SVD())

        (NL*U,V*NR)
    end)

    newstate = copy(state);

    #do the actual expansion
    for i in 1:length(state)
        al = _permute_tail(catdomain(newstate.AL[i],pexp[i][1]))
        lz = TensorMap(zeros,_lastspace(pexp[i-1][1])',domain(al))
        newstate.AL[i] = _permute_front(catcodomain(al,lz))

        ar = _permute_front(catcodomain(_permute_tail(newstate.AR[i+1]),pexp[i][2]))
        rz = TensorMap(zeros,codomain(ar),space(pexp[i+1][2],1))
        newstate.AR[i+1] = catdomain(ar,rz)

        l = TensorMap(zeros,codomain(newstate.CR[i]),space(pexp[i][2],1))
        newstate.CR[i] = catdomain(newstate.CR[i],l)
        r = TensorMap(zeros,_lastspace(pexp[i][1])',domain(newstate.CR[i]))
        newstate.CR[i] = catcodomain(newstate.CR[i],r)

        newstate.AC[i] = newstate.AL[i]*newstate.CR[i]
    end

    return newstate,envs
end

function changebonds(state::InfiniteMPS,H::InfiniteMPO,alg,envs=environments(state,H))
    (nmstate,envs) = changebonds(convert(MPSMultiline,state),convert(MPOMultiline,H),alg,envs);
    return (convert(InfiniteMPS,nmstate),envs)
end

function changebonds(state::MPSMultiline, H,alg::OptimalExpand,envs=environments(state,H))
    #=
        todo : merge this with the MPSCentergauged implementation
    =#
    #determine optimal expansion spaces around bond i
    pexp = PeriodicArray(map(Iterators.product(1:size(state,1),1:size(state,2))) do (i,j)

        ACAR = _permute_front(state.AC[i-1,j])*_permute_tail(state.AR[i-1,j+1])
        AC2 = ac2_prime(ACAR,i-1,j,state,envs)

        #Calculate nullspaces for AL and AR
        NL = leftnull(state.AL[i,j])
        NR = rightnull(_permute_tail(state.AR[i,j+1]))

        #Use this nullspaces and SVD decomposition to determine the optimal expansion space
        intermediate = adjoint(NL)*AC2*adjoint(NR)
        (U,S,V) = tsvd(intermediate,trunc=alg.trscheme,alg=SVD())

        (NL*U,V*NR)
    end)

    newstate = copy(state);

    #do the actual expansion
    for (i,j) in Iterators.product(1:size(state,1),1:size(state,2))
        al = _permute_tail(catdomain(newstate.AL[i,j],pexp[i,j][1]))
        lz = TensorMap(zeros,_lastspace(pexp[i,j-1][1])',domain(al))

        newstate.AL[i,j] = _permute_front(catcodomain(al,lz))

        ar = _permute_front(catcodomain(_permute_tail(newstate.AR[i,j+1]),pexp[i,j][2]))
        rz = TensorMap(zeros,codomain(ar),space(pexp[i,j+1][2],1))
        newstate.AR[i,j+1] = catdomain(ar,rz)

        l = TensorMap(zeros,codomain(newstate.CR[i,j]),space(pexp[i,j][2],1))
        newstate.CR[i,j] = catdomain(newstate.CR[i,j],l)
        r = TensorMap(zeros,_lastspace(pexp[i,j][1])',domain(newstate.CR[i,j]))
        newstate.CR[i,j] = catcodomain(newstate.CR[i,j],r)

        newstate.AC[i,j] = newstate.AL[i,j]*newstate.CR[i,j]
    end

    return newstate,envs
end

changebonds(state::Union{FiniteMPS,MPSComoving}, H::Hamiltonian, alg::OptimalExpand,envs=environments(state,H)) = changebonds!(copy(state),H,alg,envs)
function changebonds!(state::Union{FiniteMPS,MPSComoving}, H::Hamiltonian,alg::OptimalExpand,envs=environments(state,H))
    #inspired by the infinite mps algorithm, alternative is to use https://arxiv.org/pdf/1501.05504.pdf

    #the idea is that we always want to expand the state in such a way that there are zeros at site i
    #but "optimal vectors" at site i+1
    #so during optimization of site i, you have access to these optimal vectors :)

    for i in 1:(length(state)-1)
        ACAR = _permute_front(state.AC[i])*_permute_tail(state.AR[i+1])
        AC2 = ac2_prime(ACAR,i,state,envs)

        #Calculate nullspaces for left and right
        NL = leftnull(state.AC[i])
        NR = rightnull(_permute_tail(state.AR[i+1]))

        #Use this nullspaces and SVD decomposition to determine the optimal expansion space
        intermediate = adjoint(NL) * AC2 * adjoint(NR);
        (U,S,V) = tsvd(intermediate,trunc=alg.trscheme,alg=SVD())

        ar_re = V*NR;
        ar_le = TensorMap(zeros,codomain(state.AC[i]),space(V,1))

        (nal,nc) = leftorth(catdomain(state.AC[i],ar_le),alg=QRpos())
        nar = _permute_front(catcodomain(_permute_tail(state.AR[i+1]),ar_re));

        state.AC[i] = (nal,nc)
        state.AC[i+1] = (nc,nar)
    end

    (state,envs)
end


text:
module RTS

# IDEAS:
# 1. An RTS + Clicker combo:
#   Exponential-scaling game, but shown visually, not just via text. You start
#    out really tiny, growing your web/nest/city/etc, eventually you see some
#    enemies and you kill them, all the while expanding and zooming out.
#   Then you somehow increase the scale of the units you're placing. Maybe you
#    can combine your on-screen units into a structure, or maybe you literally
#    group some units and *define* the structure yourself. That would be cool!
#    It's like emergent scaling or something.

using SimpleDirectMediaLayer
SDL2 = SimpleDirectMediaLayer

using ApplicationBuilderAppUtils

# True if this file is being run through the interpreter, and false if being
# compiled.
debug = true

# Override SDL libs locations if this script is being compiled for mac .app builds
if get(ENV, "COMPILING_APPLE_BUNDLE", "false") == "true"
    #  (note that you can still change these values b/c no functions have
    #  actually been called yet, and so the original constants haven't been
    #  "compiled in".)
    eval(SDL2, :(libSDL2 = "libSDL2.dylib"))
    eval(SDL2, :(libSDL2_ttf = "libSDL2_ttf.dylib"))
    eval(SDL2, :(libSDL2_mixer = "libSDL2_mixer.dylib"))
    debug = false
end

const assets = "assets" # directory path for game assets relative to pwd().

include("timing.jl")
include("config.jl")
include("../assets/configs.jl")
include("player.jl")
include("game.jl")
include("display.jl")
include("keyboard.jl")
include("menu.jl")

const kGAME_NAME = "Paddle Battle"
const kSAFE_GAME_NAME = "PaddleBattle"
const kBUNDLE_ORGANIZATION = "nhdalyMadeThis"

# -------- Opening a window ---------------
# Forward reference for @cfunction
function windowEventWatcher end

# Note: These are all Atomics, since they can be modified by the
# windowEventWatcher callback, which can run in another thread!
winWidth, winHeight = Threads.Atomic{Int32}(800), Threads.Atomic{Int32}(600)
winWidth_highDPI, winHeight_highDPI = Threads.Atomic{Int32}(800), Threads.Atomic{Int32}(600)
function makeWinRenderer()
    global winWidth, winHeight, winWidth_highDPI, winHeight_highDPI

    win = SDL2.CreateWindow(kGAME_NAME,
        Int32(SDL2.WINDOWPOS_CENTERED()), Int32(SDL2.WINDOWPOS_CENTERED()), winWidth[], winHeight[],
        UInt32(SDL2.WINDOW_ALLOW_HIGHDPI|SDL2.WINDOW_OPENGL|SDL2.WINDOW_FULLSCREEN_DESKTOP|SDL2.WINDOW_SHOWN));
    SDL2.SetWindowMinimumSize(win, minWinWidth, minWinHeight)
    SDL2.AddEventWatch(@cfunction(windowEventWatcher, Cint, (Ptr{Nothing}, Ptr{SDL2.Event})), win);

    # Find out how big the created window actually was (depends on the system):
    winWidth[], winHeight[], winWidth_highDPI[], winHeight_highDPI[] = getWindowSize(win)
    #cam.w[], cam.h[] = winWidth_highDPI, winHeight_highDPI

    renderer = SDL2.CreateRenderer(win, Int32(-1), UInt32(SDL2.RENDERER_ACCELERATED | SDL2.RENDERER_PRESENTVSYNC))
    SDL2.SetRenderDrawBlendMode(renderer, UInt32(SDL2.BLENDMODE_BLEND))
    return win,renderer
end

# This huge function handles all window events. I believe it needs to be a
# callback instead of just the regular pollEvent because the main thread is
# paused while resizing, whereas this callback continues to trigger.
function windowEventWatcher(data_ptr::Ptr{Cvoid}, event_ptr::Ptr{SDL2.Event})::Cint
    global winWidth, winHeight, cam, window_paused, renderer, win
    ev = unsafe_load(event_ptr, 1)
    ee = ev._Event
    t = UInt32(ee[4]) << 24 | UInt32(ee[3]) << 16 | UInt32(ee[2]) << 8 | ee[1]
    t = SDL2.Event(t)
    if (t == SDL2.WindowEvent)
        event = unsafe_load( Ptr{SDL2.WindowEvent}(pointer_from_objref(ev)) )
        winevent = event.event;  # confusing, but that's what the field is called.
        if (winevent == SDL2.WINDOWEVENT_RESIZED || winevent == SDL2.WINDOWEVENT_SIZE_CHANGED)
            curPaused = window_paused[]
            window_paused[] = 1  # Stop game playing so resizing doesn't cause problems.
            winID = event.windowID
            eventWin = SDL2.GetWindowFromID(winID);
            if (eventWin == data_ptr)
                w,h,w_highDPI,h_highDPI = getWindowSize(eventWin)
                winWidth[], winHeight[] = w, h
                winWidth_highDPI[], winHeight_highDPI[] = w_highDPI, h_highDPI
                cam.w[], cam.h[] = winWidth[], winHeight[]
                recenterButtons!()
            end
            # Note: render after every resize event. I tried limiting it with a
            # timer, but it's hard to tune (too infrequent and the screen
            # blinks) & it didn't seem to reduce cpu significantly.
            render(sceneStack[end], renderer, eventWin)
            SDL2.GL_SwapWindow(eventWin);
            window_paused[] = curPaused  # Allow game to resume now that resizing is done.
        elseif (winevent == SDL2.WINDOWEVENT_FOCUS_LOST || winevent == SDL2.WINDOWEVENT_HIDDEN || winevent == SDL2.WINDOWEVENT_MINIMIZED)
            # Stop game playing so resizing doesn't cause problems.
            if !debug  # For debug builds, allow editing while playing
                window_paused[] = 1
            end
        elseif (winevent == SDL2.WINDOWEVENT_FOCUS_GAINED || winevent == SDL2.WINDOWEVENT_SHOWN)
            window_paused[] = 0
        end
        # Note that window events pause the game, so at the end of any window
        # event, restart the timer so it doesn't have a HUGE frame.
        start!(timer)
    end
    return 0
end

function getWindowSize(win)
    w,h,w_highDPI,h_highDPI = Int32[0],Int32[0],Int32[0],Int32[0]
    SDL2.GetWindowSize(win, w, h)
    SDL2.GL_GetDrawableSize(win, w_highDPI, h_highDPI)
    return w[],h[],w_highDPI[],h_highDPI[]
end

# Having a QuitException is useful for testing, since an exception will simply
# pause the interpreter. For release builds, the catch() block will call quitSDL().
struct QuitException <: Exception end

function quitSDL(win)
    # Need to close the callback before quitting SDL to prevent it from hanging
    # https://github.com/n0name/2D_Engine/issues/3
    SDL2.DelEventWatch(cfunction(windowEventWatcher, Cint, Tuple{Ptr{Cvoid}, Ptr{SDL2.Event}}), win);
    SDL2.Mix_CloseAudio()
    SDL2.TTF_Quit()
    SDL2.Quit()
end

# -------------- Game ------------------------------

# Game State Globals
renderer = win = nothing
p1 = Player()
p2 = Player()
cam = nothing
scoreA = 0
scoreB = 0
paused_ = true # start paused to show the initial menu.
paused = Ref(paused_)
window_paused = Threads.Atomic{UInt8}(0) # Whether or not the game should be running (if lost focus)
game_started_ = true # start paused to show the initial menu.
game_started = Ref(game_started_)
playing_ = true
playing = Ref(playing_)
debugText = false
audioEnabled = true
last_10_frame_times = [1.]
timer = WallTimer()
i = 1

sceneStack = []  # Used to keep track of the current scene

include("../assets/game_configs.jl")

e = nothing  # FOR DEBUGGING ONLY!

struct UserError
    msg
end

"""
    runSceneGameLoop(scene, renderer, win, inSceneVar::Ref{Bool})
The main game loop. Implements the poll, render, update loop, delegating to the
current scene (pushes it to the top of sceneStack).
 - Polls SDL events and passes them to `handleEvents!(scene, e, t)`.
 - Calls `render(scene, renderer, win)`.
 - Calls `performUpdates!(scene, dt)`.
This loop continues until the provided `inSceneVar` is false, then pops the
scene off the sceneStack.
"""
function runSceneGameLoop(scene, renderer, win, inSceneVar::Ref{Bool})
    global last_10_frame_times, i, e
    push!(sceneStack, scene)
    start!(timer)
    while (inSceneVar[])
        # Don't run if game is paused by system (resizing, lost focus, etc)
        while window_paused[] != 0  # Note that this will be fixed by windowEventWatcher
            _ = pollEvent!()
            sleep(0.5)  # maybe increase this?
        end
        # Reload config for debug
        if (debug) reloadConfigsFiles() end

        # Handle Events
        errorMsg = ""
        try
            hadEvents = true
            while hadEvents
                e,hadEvents = pollEvent!()
                t = getEventType(e)
                handleEvents!(scene,e,t)
            end
        catch e
            if isa(e, UserError)
                errorMsg = e.msg
            else
                rethrow()
            end
        end

        # Render
        render(scene, renderer, win)
        if (debug && debugText) renderFPS(renderer,last_10_frame_times) end
        if (errorMsg != "") renderText(renderer, cam, errorMsg, UIPixelPos(winWidth[]*1/2, 200)) end
        SDL2.RenderPresent(renderer)

        # Update
        dt = elapsed(timer)
        # Don't let the game proceed at fewer than this frames per second. If an
        # update takes too long, allow the game to actually slow, rather than
        # having too big of frames.
        min_fps = 20.0
        dt = min(dt, 1.0/min_fps)
        start!(timer)
        if debug
            last_10_frame_times = push!(last_10_frame_times, dt)
            if length(last_10_frame_times) > 10; pop!(last_10_frame_times) ; end
        end

        performUpdates!(scene, dt)
        #sleep(0.01)

        if (playing[] == false)
            throw(QuitException())
        end

        i += 1
    end
    pop!(sceneStack)
end
# Scenes must overload these functions:
#  render(scene, renderer, win)
#  handleEvents!(scene, e,t)
# Scenes can optionally overload this function:
function performUpdates!(scene, dt) end  # default

# ------------------------------------------

function pollEvent!()
    #SDL2.Event() = [SDL2.Event(NTuple{56, Uint8}(zeros(56,1)))]
    SDL_Event() = Array{UInt8}(zeros(56))
    e = SDL_Event()
    success = (SDL2.PollEvent(e) != 0)
    return e,success
end
function getEventType(e::Array{UInt8})
    # HAHA This is still pretty janky, but I guess that's all you can do w/ unions.
    bitcat(UInt32, e[4:-1:1])
end
function getEventType(e::SDL2.Event)
    e._Event[1]
end

function bitcat(outType::Type{T}, arr)::T where T<:Number
    out = zero(outType)
    for x in arr
        out = out << sizeof(x)*8
        out |= convert(T, x)  # the `convert` prevents signed T from promoting to Int64.
    end
    out
end

function renderFPS(renderer,last_10_frame_times)
    fps = Int(floor(1.0/mean(last_10_frame_times)))
    txt = "FPS: $fps"
    renderText(renderer, cam, txt, UIPixelPos(winWidth[]*1/5, 200))
end

# -------------- Game Scene ---------------------
struct GameScene end

curFoodParticles = []
function handleEvents!(scene::GameScene, e,t)
    global playing,paused
    # Handle Events
    if (t == SDL2.KEYDOWN || t == SDL2.KEYUP);  handleGameKeyPress(e,t);
    elseif (t == SDL2.MOUSEWHEEL); handleMouseScroll(e)
    #elseif (t == SDL2.MOUSEBUTTONUP || t == SDL2.MOUSEBUTTONDOWN)
    #elseif (t == SDL2.MOUSEMOTION); handleMousePan(e)
    elseif (t == SDL2.QUIT);  playing[] = false;
    end

    if (paused[])
         pause!(timer)
         enterPauseGameLoop(renderer,win)
         unpause!(timer)
    end
end
function handleMouseScroll(e)
    my = bitcat(SDL2.Sint32, e[24:-1:21])

    # TODO: move cam based on mouse pos during scroll
    #  want to keep mouse pos over same WorldPos while scrolling.
    #  ie: toScreenPos(toWorldPos(ScreenPixelPos(mx,my),cam),cam) ==
    #         toScreenPos(toWorldPos(ScreenPixelPos(mx,my),cam2),cam2)
    # zoom
    aspectRatio = cam.w[] / cam.h[]
    cam.w[] -= my * kCamZoomRate  # If scrolling up, zoom in (shrink camera).
    cam.h[] -= my * kCamZoomRate  # (positive `my` is scrolling up)
    if cam.w[] <= kCamMinSize
        cam.w[] = kCamMinSize
        cam.h[] = kCamMinSize / aspectRatio
    elseif cam.h[] <= kCamMinSize
        cam.w[] = kCamMinSize * aspectRatio
        cam.h[] = kCamMinSize
    end
end

#function handleGameMouseClick!(e, clickType)
#    mx = bitcat(UInt32, e[24:-1:21])
#    my = bitcat(UInt32, e[28:-1:25])
#    mButton = e[17]
#    if mButton == SDL2.BUTTON_RIGHT
#        p = toWorldPos(UIPixelPos(mx,my), cam)
#    end
#end

function handleGameKeyPress(e,t)
    global paused,debugText
    keySym = getKeySym(e)
    keyDown = (t == SDL2.KEYDOWN)
    keyRepeat = (getKeyRepeat(e) != 0)
    println("repeat: $(getKeyRepeat(e))")

    if (keySym == keySettings[:keyP1Collector])
        buyCollector(p1)
    elseif (keySym == keySettings[:keyP2Collector])
        buyCollector(p2)
    elseif (keySym == keySettings[:keyP1Fighter])
        buyFighter(p1)
    elseif (keySym == keySettings[:keyP2Fighter])
        buyFighter(p2)
    elseif (keySym == keySettings[:keyP1Attack] && keyDown && !keyRepeat)
        attackToggle(p2, p1)
    elseif (keySym == keySettings[:keyP2Attack] && keyDown && !keyRepeat)
        attackToggle(p1, p2)
    elseif (keySym == SDL2.SDLK_RIGHT)
        cam.pos += Vector2D(kCamPanRate,0)
    elseif (keySym == SDL2.SDLK_LEFT)
        cam.pos += Vector2D(-kCamPanRate,0)
    elseif (keySym == SDL2.SDLK_UP)
        cam.pos += Vector2D(0,kCamPanRate)
    elseif (keySym == SDL2.SDLK_DOWN)
        cam.pos += Vector2D(0,-kCamPanRate)
    else
        # Fallback
        handlePauseKeyPress(e,t)
    end
end

unitRenderColor(::Type{Fighter}) = kFighterColor
unitRenderColor(::Type{Collector}) = kCollectorColor
playerRenderColor(p::Player) = (if (p === p1) kP1Color elseif (p === p2) kP2Color end)

function render(scene::GameScene, renderer, win)
    global scoreA,scoreB,last_10_frame_times,paused,playing

    SetRenderDrawColor(renderer, kBackgroundColor)
    SDL2.RenderClear(renderer)

    renderScore(renderer)

    #renderFoodParticles(cam,renderer, curFoodParticles)

    for u in p1.units.units
        render(u, kP1Color, cam, renderer)
    end
    for u in p2.units.units
        render(u, kP2Color, cam, renderer)
    end

    # UI text at bottom
    renderText(renderer, cam, "P1: $(display_key_setting(:keyP1Collector)): collector (\$$(build_cost(Collector)))   $(display_key_setting(:keyP1Fighter)): scout (\$$(build_cost(Fighter))) $(display_key_setting(:keyP1Attack)): Attack",
               UIPixelPos(5, winHeight[] - kUIFontSize)
               ; fontSize=kUIFontSize, align=leftJustified)
    renderText(renderer, cam, "P2: $(display_key_setting(:keyP2Collector)): collector (\$$(build_cost(Collector)))   $(display_key_setting(:keyP2Fighter)): scout (\$$(build_cost(Fighter))) $(display_key_setting(:keyP2Attack)): Attack",
               UIPixelPos(winWidth[]-5, winHeight[] - kUIFontSize)
               ; fontSize=kUIFontSize, align=rightJustified)

    # BuildOp UI
    function drawBuildOps(p, xPos)
        buildOpsHeight = kBuildOpsRenderHeight
        for b in p.build_ops
            pos = UIPixelPos(xPos, winHeight[] - buildOpsHeight)
            percent = time_remaining(b) / b.buildLength
            renderProgressBar(percent, cam, renderer, pos,
                    UIPixelDims(200, 10), blendAlphaColors(playerRenderColor(p), unitRenderColor(b.unitType)),
                    kBuildOpsBgColor, healthBarOutlineColor)
            buildOpsHeight += 10
        end
    end
    drawBuildOps(p1, winWidth[]/4)
    drawBuildOps(p2, winWidth[]*3/4)

end

function renderScore(renderer)
    # Size the text with a single-digit score so it doesn't move when score hits double-digits.
    txtDims = toUIPixelDims(sizeText("Player 1: 0", defaultFontName, defaultFontSize), cam)
    hcat_render_text(["Player 1: $(floor(p1.money))","Player 2: $(floor(p2.money))"], renderer, cam,
         100, UIPixelPos(screenCenterX(), 20)
         ; fixedWidth = txtDims.w)
end

function performUpdates!(scene::GameScene, dt)
    global p1, p2
    reloadConfigsFiles(["../assets/configs.jl", "../assets/game_configs.jl"])
    update!(p1, dt)
    update!(p2, dt)

    performAttackUpdate!(dt)

    moveCamIfMouseOnEdge!(dt)
    #updateFoodParticles(dt)
end

foodSheetDriftVel = Vector2D(0,0)
function updateFoodParticles(dt)
    global curFoodParticles,worldFoodSheetOffset,foodSheetDriftVel
    curFoodParticles = foodDistribution(cam.pos, WorldDims(cam.w[],cam.h[]))

    worldFoodSheetOffset += foodSheetDriftVel * dt
    #range = (kFoodSheetDriftMaxSpeed_sqrd - magSqrd(foodSheetDriftVel))
    driftAccel = Vector2D(rand(-kFoodSheetDriftAccel:0.01:kFoodSheetDriftAccel, 2)...) * dt
    foodSheetDriftVel += driftAccel
    if magSqrd(foodSheetDriftVel) > kFoodSheetDriftMaxSpeed_sqrd
        foodSheetDriftVel -= foodSheetDriftVel * (0.1 * dt)
    end

    collectorDims = WorldDims(collectorRenderWidth,collectorRenderWidth)
    for p in curFoodParticles
        for c in p1.units.collectors
            if overlapping(toWorldPos(p.pos), p.size, c.pos, collectorDims)
                destroyFoodParticle(p)
                p1.money += p.amountFood
            end
        end
    end
end

function moveCamIfMouseOnEdge!(dt)
    x,y = Int[1], Int[1]
    #SDL2.PumpEvents()
    SDL2.GetMouseState(pointer(x), pointer(y))
    if (winWidth[] - x[]) < kMouseEdgeDetectionWidth
        cam.pos += Vector2D(kCamPanRate,0) * dt
    elseif x[] < kMouseEdgeDetectionWidth
        cam.pos += Vector2D(-kCamPanRate,0) * dt
    end
    if y[] < kMouseEdgeDetectionWidth
        cam.pos += Vector2D(0,kCamPanRate) * dt
    elseif (winHeight[] - y[]) < kMouseEdgeDetectionWidth
        cam.pos += Vector2D(0,-kCamPanRate) * dt
    end
end

function enterWinnerGameLoop(renderer,win, winnerName)
    # Reset the buttons to the beginning of the game.
    buttons[:bRestart].enabled = false # Nothing to restart
    buttons[:bNewContinue].text = "New Game"
    global paused,game_started; paused[] = true; game_started[] = false;

    scene = PauseScene("$winnerName wins!!", "")
    runSceneGameLoop(scene, renderer, win, paused)

    # When the pause scene returns, reset the game before starting.
    resetGame()
end
function resetGame()
    global scoreA,scoreB, p1, p2, cam, attackingMap
    cam = Camera(WorldPos(0,0),
             Threads.Atomic{Float32}(winWidth[]),
             Threads.Atomic{Float32}(winHeight[]))

    scoreB = scoreA = 0
    p1 = Player()
    add_unit!(p1.units, Collector(p1.units, WorldPos(-250,50)))
    add_unit!(p1.units, Fighter(p1.units, WorldPos(-200,200)))
    add_unit!(p1.units, Fighter(p1.units, WorldPos(-200,-200)))

    p2 = Player()
    add_unit!(p2.units, Collector(p2.units, WorldPos(250,50)))
    add_unit!(p2.units, Fighter(p2.units, WorldPos(200,200)))
    add_unit!(p2.units, Fighter(p2.units, WorldPos(200,-200)))

    attackingMap = Dict(p1=>false, p2=>false)
end

mutable struct KeyControls
    rightDown::Bool
    leftDown::Bool
    KeyControls() = new(false,false)
end
const paddleAKeys = KeyControls()
const paddleBKeys = KeyControls()
mutable struct GameControls
    escapeDown::Bool
    GameControls() = new(false)
end
const gameControls = GameControls()

randScreenPos(cam) = UIPixelPos(rand(0:winWidth[]), rand(0:winHeight[]))
randWorldPosOnScreen(cam) = toWorldPos(randScreenPos(cam), cam)
getKeySym(e) = bitcat(UInt32, e[24:-1:21])
getKeyRepeat(e) = bitcat(UInt8, e[14:-1:14])
function handlePauseKeyPress(e,t)
    global paused,debugText
    keySym = getKeySym(e)
    keyDown = (t == SDL2.KEYDOWN)

    if (keySym == SDL2.SDLK_ESCAPE)
        if (!gameControls.escapeDown && keyDown)
            if game_started[]  # Escape shouldn't start the game.
                paused[] = !paused[]
            end
        end
        gameControls.escapeDown = keyDown
    elseif (keySym == SDL2.SDLK_BACKQUOTE)
        keyDown && (debugText = !debugText)
    end
end

function randPosAroundCollectors(p, randRange)
    if isempty(p.units.collectors)
        randPos = randWorldPosOnScreen(cam)
    else
        randCollectorPos = rand(p.units.collectors).pos
        randPos = WorldPos(rand(-randRange:randRange) + randCollectorPos.x,
                        rand(-randRange:randRange) + randCollectorPos.y)
    end
end
function randPosAroundFighters(p, randRange)
    if isempty(p.units.fighters)
        throw(UserError("Must have a scout to create Collector."))
    else
        randFighterPos = rand(p.units.fighters).pos
        randPos = WorldPos(rand(-randRange:randRange) + randFighterPos.x,
                        rand(-randRange:randRange) + randFighterPos.y)
    end
end
function buyCollector(p)
    purchase_collector!(p, randPosAroundFighters(p, kRandPurchaseCollectorPosRange))
end
function buyFighter(p)
    randRange = 300
    if isempty(p.units.collectors)
        audioEnabled && SDL2.Mix_PlayChannel( Int32(-1), badKeySound, Int32(0) )
    else
        purchase_fighter!(p, randPosAroundCollectors(p, kRandPurchaseFighterPosRange))
    end
end

attackingMap = Dict(p1=>false, p2=>false)
function attackToggle(pTarget, p)
    attackingMap[p] = !attackingMap[p]
    if attackingMap[p]
        # sick units on pTarget's units
        for f in p.units.fighters
            if isempty(pTarget.units.units)
                break
            end
            f.attackTargetUnit = rand(pTarget.units.units)
        end
    else
        for f in p.units.fighters
            f.attackTargetUnit = nothing
        end
    end
end
function performAttackUpdate!(dt)
    # Move units towards pTarget's units
    units_to_destroy = []
    for f in [p1.units.fighters..., p2.units.fighters...]
        # Stop attacking if target already dead or not attacking anything
        if f.attackTargetUnit == nothing || f.attackTargetUnit.health <= 0
            f.attackTargetUnit = nothing
            continue
        end
        targetVec = f.attackTargetUnit.pos - f.pos
        if (magSqrd(targetVec) > 1) targetVec = unitVec(targetVec) end
        f.pos += targetVec * kFighterSpeed * dt

        if overlapping(f.pos,kFighterSize, f.attackTargetUnit.pos, kFighterSize)
            attack!(f.attackTargetUnit, f)
            push!(units_to_destroy, f)
        end
    end
    for f in units_to_destroy
        set_health!(f, 0)
        destroy_unit!(f)
    end
end

# -------------- Pause Scene ---------------------
struct PauseScene
    titleText::String
    subtitleText::String
end
function enterPauseGameLoop(renderer,win)
    global paused
    scene = PauseScene("$kGAME_NAME", "Main Menu")
    runSceneGameLoop(scene, renderer, win, paused)
end
function handleEvents!(scene::PauseScene, e,t)
    global playing,paused
    # Handle Events
    if (t == SDL2.KEYDOWN || t == SDL2.KEYUP);  handlePauseKeyPress(e,t);
    elseif (t == SDL2.MOUSEBUTTONUP || t == SDL2.MOUSEBUTTONDOWN)
        b = handleMouseClickButton!(e,t);
        if (b != nothing); run(b); end
    elseif (t == SDL2.QUIT);
        playing[]=false; paused[]=false;
    end
end

heartIcon = nothing
jlLogoIcon = nothing
function render(scene::PauseScene, renderer, win)
    global heartIcon, jlLogoIcon
    if heartIcon == nothing || jlLogoIcon == nothing
        heart_surface = SDL2.LoadBMP("assets/heart.bmp")
        heartIcon = SDL2.CreateTextureFromSurface(renderer, heart_surface) # Will be C_NULL on failure.
        SDL2.FreeSurface(heart_surface)
        jlLogo_surface = SDL2.LoadBMP("assets/jllogo.bmp")
        jlLogoIcon = SDL2.CreateTextureFromSurface(renderer, jlLogo_surface) # Will be C_NULL on failure.
        SDL2.FreeSurface(jlLogo_surface)
    end
    screenRect = SDL2.Rect(0,0, winWidth_highDPI[], winWidth_highDPI[])
    # First render the scene under the pause menu so it looks like the pause is over it.
    if (length(sceneStack) > 1) render(sceneStack[end-1], renderer, win) end
    color = kBackgroundColor
    SDL2.SetRenderDrawColor(renderer, Int64(color.r), Int64(color.g), Int64(color.b), 200) # transparent
    SDL2.RenderFillRect(renderer, Ref(screenRect))
    renderText(renderer, cam, scene.titleText, screenOffsetFromCenter(0,-149)
               ; fontSize=kPauseSceneTitleFontSize)
    renderText(renderer, cam, scene.subtitleText, screenOffsetFromCenter(0,-109); fontSize = kPauseSceneSubtitleFontSize)
    for b in values(buttons)
        render(b, cam, renderer)
    end
    renderText(renderer, cam, "Player 1 Controls", UIPixelPos(paddleAControlsX(),winHeight[]-169); fontSize = kControlsHeaderFontSize)
    renderText(renderer, cam, "Player 2 Controls", UIPixelPos(paddleBControlsX(),winHeight[]-169); fontSize = kControlsHeaderFontSize)
    renderText(renderer, cam, kCopyrightNotices[1], UIPixelPos(screenCenterX(), winHeight[] - 20);
            fontName="assets/fonts/FiraCode/ttf/FiraCode-Regular.ttf",
            fontSize=10)
    renderText(renderer, cam, kCopyrightNotices[2], UIPixelPos(screenCenterX(), winHeight[] - 8);
            fontName="assets/fonts/FiraCode/ttf/FiraCode-Regular.ttf",
            fontSize=10)

    _, heartPos, _, jlLogoPos =
      hcat_render_text(kProgrammedWithJuliaText, renderer, cam,
         0, UIPixelPos(screenCenterX(), winHeight[] - 35);
          fontName="assets/fonts/FiraCode/ttf/FiraCode-Regular.ttf",
          fontSize=16)
    render(heartIcon, heartPos, cam, renderer; size=UIPixelDims(16,16))
    render(jlLogoIcon, jlLogoPos, cam, renderer; size=UIPixelDims(16,16))
end

clickedButton = nothing
function handleMouseClickButton!(e, clickType)
    global clickedButton
    mx = bitcat(UInt32, e[24:-1:21])
    my = bitcat(UInt32, e[28:-1:25])
    mButton = e[17]
    if mButton != SDL2.BUTTON_LEFT
        return
    end
    didClickButton = false
    for b in values(buttons)
        if mouseOnButton(UIPixelPos(mx,my),b,cam)
            if (clickType == SDL2.MOUSEBUTTONDOWN)
                clickedButton = b
                didClickButton = true
                break
            elseif clickedButton == b && clickType == SDL2.MOUSEBUTTONUP
                clickedButton = nothing
                didClickButton = true
                return b
            end
        end
    end
    if clickedButton != nothing && clickType == SDL2.MOUSEBUTTONUP && didClickButton == false
        clickedButton = nothing
    end
    return nothing
end

function mouseOnButton(m::UIPixelPos, b::CheckboxButton, cam)
    return mouseOnButton(m, b.button, cam)
end
function mouseOnButton(m::UIPixelPos, b::AbstractButton, cam)
    if (!b.enabled) return false end
    topLeft = topLeftPos(b.pos, b.dims)
    if m.x > topLeft.x && m.x <= topLeft.x + b.dims.w &&
        m.y > topLeft.y && m.y <= topLeft.y + b.dims.h
        return true
    end
    return false
end

function load_audio_files()
    global pingSound, scoreSound, badKeySound
    pingSound = SDL2.Mix_LoadWAV( "$assets/ping.wav" );
    scoreSound = SDL2.Mix_LoadWAV( "$assets/score.wav" );
    badKeySound = SDL2.Mix_LoadWAV( "$assets/ping.wav" );
end

function game_main(ARGS)
    global renderer, win, paused,game_started, cam
    win = nothing
    try
        SDL2.init()
        ApplicationBuilderAppUtils.cd_to_bundle_resources()
        #init_prefspath()
        #load_prefs_backup()
        load_audio_files()
        music = SDL2.Mix_LoadMUS( "$assets/music.wav" );
        win,renderer = makeWinRenderer()
        global paused,game_started; paused[] = true; game_started[] = false;
        # Warm up
        for i in 1:3
            pollEvent!()
            SDL2.SetRenderDrawColor(renderer, 200, 200, 200, 255)
            SDL2.RenderClear(renderer)
            SDL2.RenderPresent(renderer)
            #sleep(0.01)
        end
        audioEnabled && SDL2.Mix_PlayMusic( music, Int32(-1) )
        recenterButtons!()
        resetGame();  # Initialize game stuff.
        #println("Preferences file: $(repr(prefsfile))")
        playing[] = paused[] = true
        scene = GameScene()
        runSceneGameLoop(scene, renderer, win, playing)
    catch e
        if isa(e, QuitException)
            quitSDL(win)
        else
            rethrow()  # Every other kind of exception
        end
    end
        return 0
end
Base.@ccallable function julia_main(ARGS::Vector{String})::Cint
    game_main(ARGS)
end

#julia_main([""])  # no julia_main if currently compiling.

end


text:
"""
    `individual_summary(df,xvar,yvar; summary = mean, err = :MouseID)`
    computes a summary statistic (default = mean) for each individual
    (default = :MouseID) used to later calculate the error across individual
"""

function individual_summary(df,xvar,yvar; summary = mean, err = :MouseID)
    #evaluate if it's going over multiple grouping columns
    multiple = typeof(xvar) <: AbstractVector && sizeof(xvar) > 1
    if multiple
        gdc = groupby(df,vcat(xvar,err))
    else
        gdc = groupby(df,[xvar,err])
    end
    df1 = combine(gdc,yvar => summary => yvar)
    sort!(df1,xvar)
    if multiple
        label_df = unique(df1[:,xvar])
    else
        label_df = unique(df1[:,[xvar]])
    end
    label_df[!,:xpos] = 1:nrow(label_df)
    leftjoin(df1,label_df; on = xvar)
    try
        firstval = union(df1[:,xvar])[1]
        df1[!,:xpos] = [v == firstval ? 1 : 2  for v in df1[:,xvar]]
    catch
        println("Can't define x position")
    end
    return df1
end

"""
    `group_summary(df1,xvar,yvar; normality = true)`
    using bootstrap computes central measure plus 95% CI
    if normality is true it uses mean as central measure
    else it uses median
"""

function group_summary(df1,xvar,yvar; normality = true)
    if normality
        central = mean
    else
        central = median
    end
    df2 = combine(groupby(df1,xvar)) do dd
        # using bootstrap to calculate the central point and the 95% CI
        b = bootstrap(central, dd[:,yvar], BasicSampling(100000))
        #bootstrap confint returns a tuple containing a tuple with the statistic
        #estimates and the CI lower and uppper bounds
        m, lower, upper = confint(b,PercentileConfInt(0.95))[1]
        ci1 = m - lower
        ci2 = upper - m
        (Central = m, ERR = (ci1,ci2))
    end
    return df2
end

function MWU_test(df1,group,yvar)
    cases = union(df1[:,group])
    case1 = df1[df1[:,group] .== cases[1], yvar]
    case2 = df1[df1[:,group] .== cases[2], yvar]
    test = MannWhitneyUTest(case1,case2)
end

function survivalrate_algorythm(var; step = 0.05, xaxis = nothing)
    isnothing(xaxis) && (xaxis = range(extrema(var)..., step = step))
    survival = 1 .- ecdf(var).(xaxis)
    return (Xaxis = collect(xaxis), fy = survival)
end

function cumulative_algorythm(var; step = 0.05, xaxis = nothing)
    isnothing(xaxis) && (xaxis = range(extrema(var)..., step = step))
    cum = ecdf(var).(xaxis)
    return (Xaxis = collect(xaxis), fy = cum)
end

function hazardrate_algorythm(var; step = 0.05, xaxis = nothing)
    isnothing(xaxis) && (xaxis = range(extrema(var)..., step = step))
    survival = 1 .- ecdf(var).(xaxis)
    hazard = -pushfirst!(diff(survival),0)./survival
    return (Xaxis = collect(xaxis), fy = hazard)
end


text:
using Oceananigans.Fields: validate_field_tuple_grid

#####
##### Forms for IncompressibleModel constructor
#####

DiffusivityFields(diffusivities::NamedTuple, arch, grid, tracer_names, bcs, closure) =
    validate_field_tuple_grid("diffusivities", diffusivities, grid)

DiffusivityFields(::Nothing, arch, grid, tracer_names, bcs, closure) =
    DiffusivityFields(arch, grid, tracer_names, bcs, closure)

#####
##### Closures without precomputed diffusivities
#####

DiffusivityFields(arch, grid, tracer_names, bcs, closure) = nothing

#####
##### Closure tuples
#####

DiffusivityFields(arch, grid, tracer_names, bcs, closure_tuple::Tuple) =
    Tuple(DiffusivityFields(arch, grid, tracer_names, bcs, closure) for closure in closure_tuple)


text:
#-----------------------------------------------------------------------
#   TOLERANCES
#-----------------------------------------------------------------------

δtol = 1e-8
tol  = 1e-13

#-----------------------------------------------------------------------
#   GEODESIC SOLVER TEST
#-----------------------------------------------------------------------

Z0 = [0.;0.;0.;0.;1.;1.;0.;0.]      # Index raised in p's here
ZF = [1.;1.;0.;0.;-1.;1.;0.;0.]     # Index lowered in p's here
@test squirrel.gsolve(  Z0[1:4] , Z0[5:8] , η , δtol
                        , AutoVern7(Rodas5()) ) ≈ ZF

Z0 = [0.;0.;0.;0.;1.;0.;1.;0.]      # Index raised in p's here
ZF = [1.;0.;1.;0.;-1.;0.;1.;0.]     # Index lowered in p's here
@test squirrel.gsolve(  Z0[1:4] , Z0[5:8] , η , δtol
                        , AutoVern7(Rodas5()) ) ≈ ZF

Z0 = [0.;0.;0.;0.;1.;0.;0.;1.]      # Index raised in p's here
ZF = [1.;0.;0.;1.;-1.;0.;0.;1.]     # Index lowered in p's here
@test squirrel.gsolve(  Z0[1:4] , Z0[5:8] , η , δtol
                        , AutoVern7(Rodas5()) ) ≈ ZF

#-----------------------------------------------------------------------
#   V34 TEST
#-----------------------------------------------------------------------

v3 = rand(3)
@test length(squirrel.V34(v3)) ≈ 4
@test squirrel.V34(v3)[2:4] == v3

#-----------------------------------------------------------------------
#   VidF TEST
#-----------------------------------------------------------------------

Zi = zeros(8,4)
Zi[6:8,1] = ones(3)
Zi[6:8,2] = 2*ones(3)
Zi[6:8,3] = 3*ones(3)
Zi[6:8,4] = 4*ones(3)

Vid = squirrel.VidF( Zi )

@test Vid[1:3] ≈ ones(3)
@test Vid[4:6] ≈ 2*ones(3)
@test Vid[7:9] ≈ 3*ones(3)
@test Vid[10:12] ≈ 4*ones(3)

#-----------------------------------------------------------------------
#   INITIAL DATA GENERATION
#-----------------------------------------------------------------------

X   = [-0.9938120402728282 -0.9938120402728282 -0.9938120402728282 -0.9938120402728282;
        0.43329524598249924 0.6468681369385507 0.615089937176196 0.5122363212193108;
        0.738110796898415 0.7597257376287113 0.2839892724192345 0.0952119165738411;
        0.5369302954880543 0.15569297625663295 0.7506786499110514 0.8654592211384085]

Xc  = [0.006187959727171788, 0.008549046454540871, 0.004581697263911051, 0.006340316336591607]

Zi  = [-0.9938120402728282 -0.9938120402728282 -0.9938120402728282 -0.9938120402728282;
        0.43329524598249924 0.6468681369385507 0.615089937176196 0.5122363212193108;
        0.738110796898415 0.7597257376287113 0.2839892724192345 0.0952119165738411;
        0.5369302954880543 0.15569297625663295 0.7506786499110514 0.8654592211384085;
        1.0 1.0 1.0 1.0;
       -0.4247461995279584 -0.6383190904840098 -0.6065408907216551 -0.5036872747647699;
       -0.733529099634504 -0.7551440403648003 -0.27940757515532344 -0.09063021930993004;
       -0.5305899791514627 -0.14935265992004135 -0.7443383335744598 -0.8591189048018169]

V   = [-1.0 -1.0 -1.0 -1.0;
        0.4247461995279584 0.6383190904840098 0.6065408907216551 0.5036872747647699;
        0.733529099634504 0.7551440403648003 0.27940757515532344 0.09063021930993004;
        0.5305899791514627 0.14935265992004135 0.7443383335744598 0.8591189048018169]

Vid  = vcat(-V[2:4,1],-V[2:4,2],-V[2:4,3],-V[2:4,4])

#-----------------------------------------------------------------------
#   ZERO FUNCTION TEST
#-----------------------------------------------------------------------

@test squirrel.zF( Vid , Zi , η , δtol ) ≈ zeros(12) atol=tol

#-----------------------------------------------------------------------
#   GEJAC TEST
#-----------------------------------------------------------------------

for i=1:4
    local res = squirrel.gejac( Zi[1:4,i] , Zi[5:8,i] , η , δtol )
    @test res[1][1:4] ≈ Xc atol=tol
    @test res[1][6:8] ≈ -V[2:4,i] atol=tol
    @test res[2][1,:] ≈ -V[2:4,i] atol=tol
end

#-----------------------------------------------------------------------
#   JACOBIAN CALCULATOR TEST
#-----------------------------------------------------------------------
global resx = squirrel.geocJ( Zi , η , δtol )

@test resx[1][6:8,:]       ≈ -V[2:4,:] atol=tol

for i=1:4
    @test resx[1][1:4,i]   ≈ Xc atol=tol
end

@test resx[2][1,1:3]       ≈ -V[2:4,1] atol=tol
@test resx[2][5,1:3]       ≈ -V[2:4,1] atol=tol
@test resx[2][9,1:3]       ≈ -V[2:4,1] atol=tol

@test resx[2][1,4:6]       ≈ V[2:4,2] atol=tol
@test resx[2][5,7:9]       ≈ V[2:4,3] atol=tol
@test resx[2][9,10:12]     ≈ V[2:4,4] atol=tol

@test resx[2][5:12,  4:6]  ≈ zeros(8,3)
@test resx[2][1:4,   7:9]  ≈ zeros(4,3)
@test resx[2][9:12,  7:9]  ≈ zeros(4,3)
@test resx[2][1:8, 10:12]  ≈ zeros(8,3)

@test resx[2][2:4,1:3]     ≈ 1.0*I(3)
@test resx[2][6:8,1:3]     ≈ 1.0*I(3)
@test resx[2][10:12,1:3]   ≈ 1.0*I(3)
@test resx[2][2:4,4:6]     ≈ -1.0*I(3)
@test resx[2][6:8,7:9]     ≈ -1.0*I(3)

#-----------------------------------------------------------------------
#   zFc TEST
#-----------------------------------------------------------------------

Zf = ones(8,4)
Zf[:,2] = -ones(8)
Zf[:,3] = -2*ones(8)
Zf[:,4] = -3*ones(8)

@test squirrel.zFc( Zf )[1:4]  ≈ 2*ones(4)
@test squirrel.zFc( Zf )[5:8]  ≈ 3*ones(4)
@test squirrel.zFc( Zf )[9:12] ≈ 4*ones(4)

#-----------------------------------------------------------------------
#   INITIAL DATA CORRECTOR TEST
#-----------------------------------------------------------------------

Zf = squirrel.idf( Zi , η , δtol , 24 )
@test Zf[1:4] ≈ Zi[1:4] atol=δtol*mean(Zi)/10
@test Zf[6:8] ≈ Zi[6:8] atol=δtol*mean(Zi)/10

#-----------------------------------------------------------------------
#   INITIAL DATA CONSTRUCTOR
#-----------------------------------------------------------------------

gk   = x->gks(x,0.5) ;

XE   = [4.8552418244084965 4.845669922167221 4.856689375629473 4.8206419799059494 4.823208150136851 4.856202508967677  ;
        7.299483910135271 6.8054911167005425 7.1810230696314425 6.153377751048404 7.746460952899268 7.227281564130913  ;
        10.349124252528748 9.710191586980727 10.434190218681847 10.14678773174398 9.282597162251408 10.217530411064208 ;
        5.265994805337563 6.471968677435302 5.6150675825992185 5.951777443574898 5.391088065951989 6.0625298198886775  ] ;

Xtar = [5.855953201888279 ; 6.78444938070708 ; 9.525600927112148 ; 5.492794659774428] ;

#-----------------------------------------------------------------------
#   SINGLE LOCATOR TEST
#-----------------------------------------------------------------------

Xs = squirrel.locator4( XE[:,1:4] ,
                        squirrel.locator4FHC22(XE[:,1:4])[1] , gk ,
                        1e-10 , 24 , false )

@test Xs ≈ Xtar  atol=5e-14

#-----------------------------------------------------------------------
#   MULTI LOCATOR TEST
#-----------------------------------------------------------------------

ne  = 6

Xsc  = squirrel.locator(  XE , gk , 1e-10 , ne , 24 , 2e1 , Double64  )

@test Xsc ≈ Xtar  atol=5e-14


text:
@testset "48.rotate-image.jl" begin
    @test rotate_48([
        1 2 3
        4 5 6
        7 8 9
    ]) == [7 4 1
        8 5 2
        9 6 3]
    @test rotate_48([
        5 1 9 11
        2 4 8 10
        13 3 6 7
        15 14 12 16]) == [
        15 13 2 5
        14 3 4 1
        12 6 8 9
        16 7 10 11]
    @test rotate_48(fill(1, 1, 1)) == fill(1, 1, 1)
    @test rotate_48([
        1 2
        3 4]) == [
        3 1
        4 2]
end


text:
# This file is a part of EmpiricalDistributions.jl, licensed under the MIT License (MIT).

import Test
Test.@testset "Package EmpiricalDistributions" begin

include("test_utils.jl")
include("test_uv_binned_dist.jl")
include("test_mv_binned_dist.jl")

end # testset


text:
function applyperipheryics!(system::LegSystem,old=Dict("a"=>0),restart="no")
    numupper = 0;
    if restart == "yes"
        branches = old["branches"];
    end
    for i = 1:length(system.branches.ID)
        if isempty(system.branches.children[i])
            if restart == "no"
                system.branches.term[i].P[1,1] = (
                    system.branches.P[i][system.solverparams.JL,1]*mmHgToPa);
            elseif restart == "yes"
                term = branches["term"][i];
                system.branches.term[i].P[1,1] = term["P"][end,1];
                system.branches.term[i].Q[1,1] = term["Q"][end,1];
            end
        end
    end
end

max_stars_repo_name: dcanuto/leg-hemodynamics

max_issues_repo_name: dcanuto/leg-hemodynamics

max_forks_repo_name: dcanuto/leg-hemodynamics


text:
# Julia wrapper for header: /usr/include/libavutil/des.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_des_alloc,
    av_des_init,
    av_des_crypt,
    av_des_mac

function av_des_alloc()
    ccall((:av_des_alloc, libavutil), Ptr{AVDES}, ())
end

function av_des_init(d, key, key_bits::Integer, decrypt::Integer)
    ccall((:av_des_init, libavutil), Cint, (Ptr{AVDES}, Ptr{UInt8}, Cint, Cint), d, key, key_bits, decrypt)
end

function av_des_crypt(d, dst, src, count::Integer, iv, decrypt::Integer)
    ccall((:av_des_crypt, libavutil), Cvoid, (Ptr{AVDES}, Ptr{UInt8}, Ptr{UInt8}, Cint, Ptr{UInt8}, Cint), d, dst, src, count, iv, decrypt)
end

function av_des_mac(d, dst, src, count::Integer)
    ccall((:av_des_mac, libavutil), Cvoid, (Ptr{AVDES}, Ptr{UInt8}, Ptr{UInt8}, Cint), d, dst, src, count)
end


text:
import Downloads
using CSV, DataFrames
using Test

function learn_basic()
    url = "https://raw.githubusercontent.com/Arkoniak/UrlDownload.jl/master/data/ext.csv"

    # Download to a file
    tmp = Downloads.download(url)
    csv = CSV.File(tmp)
    rm(tmp)
    df = DataFrame(csv)

    @test [1,3] == df[!,:x]
    @test [2,4] == df[!,:y]

    # Download to memory
    io = IOBuffer()
    Downloads.download(url, io)
    csv = CSV.File(take!(io))
    df = DataFrame(csv)

    @test [1,3] == df[!,:x]
    @test [2,4] == df[!,:y]
end

@testset "Base" begin
    learn_basic()
end

nothing


text:
# code for Video 4.3:
fname = "iris.csv"
using DataFrames
data = readtable(fname, separator = ',')

### scatter plots:
using Gadfly
set_default_plot_size(15cm, 15cm/golden)
p = plot(data, x="sepal_width", y="sepal_length", Geom.point,
         Guide.xlabel("Sepal Width"),
         Guide.ylabel("Sepal Length"),
         Guide.title("Scatterplot of Sepal Length versus Width") )
draw(PDF("Length_versus_Width.pdf", 14cm, 10cm), p)

cor(data[:sepal_width], data[:sepal_length]) # -0.10936924995064931

# add colour dimension:
p = plot(data, x="sepal_width", y="sepal_length", color="species",
         Guide.xlabel("Sepal Width"),
         Guide.ylabel("Sepal Length"),
         Guide.title("Scatterplot of Sepal Length versus Width") )
draw(PDF("Length_versus_Width_Color.pdf", 14cm, 10cm), p)

cor(data[1:50, :sepal_width], data[1:50, :sepal_length]) # 0.7467803732639265
cor(data[51:100, :sepal_width], data[51:100, :sepal_length]) # 0.5259107172828241
cor(data[101:150, :sepal_width], data[101:150, :sepal_length]) # 0.4572278163941129

### linear regression:
x = convert(Array, data[1:50, :sepal_width])
y = convert(Array, data[1:50, :sepal_length])
a, b = linreg(x, y)
# 2-element Array{Float64,1}:
#  2.64466
#  0.690854
# sepal_length = 2.64466 + 0.690854 * sepal_width

# adding a linear regression fit to the plot:
using Winston
pl = plot(data[1:50, :sepal_width], data[1:50, :sepal_length], "g^") # use the symbol for green triangles
s = Slope(0.690854, (0, 2.64466), colorant="blue", linekind="dashed") # this is the linear fit
add(pl, s)
display(pl) # not necessary in the REPL
savefig(pl, "Winston_Linreg.png")

using PyPlot
x = convert(Array, data[1:50, :petal_width])
y = convert(Array, data[1:50, :petal_length])
a, b = linreg(x, y)
# 2-element Array{Float64,1}:
#  1.34304
#  0.495739
cor(x, y) # correlation = 0.30630821115803575
plot(x, y, "o")
plot(x, [a+b*i for i in x])

### histogram:
Gadfly.plot(data, x="petal_width", Geom.histogram(bincount=30))

### boxplot with Gadfly:
Gadfly.plot(data, x = "species", y = "sepal_length", Geom.boxplot)

max_stars_repo_path: Section 4/4_3.jl

max_stars_repo_name: lytemar/Julia-for-Data-Science-Video

max_issues_repo_path: Section 4/4_3.jl

max_issues_repo_name: lytemar/Julia-for-Data-Science-Video

max_forks_repo_path: Section 4/4_3.jl

max_forks_repo_name: lytemar/Julia-for-Data-Science-Video


text:
module IOUtils

"""
    WriteBlockingIO(buf=0; spawn=false) do io ... end
    WriteBlockingIO(buf=0)

Create a blocking IO object which blocks character entry (when the buffer is full). This is
useful if you have a function that may take a really long time to potentially write
thousands or millions of characters, but you only want the first 10.

If you create the WriteBlockingIO object with an associated function, that function will be
spawned concurrently. If you close the io object, an exception will be thrown on subsequent
writes, killing the spawned Task if unhandled.

# Example
```julia-repl
julia> io = WriteBlockingIO() do io
           println(io, collect(1:1_000_000))
       end
WriteBlockingIO(Channel{UInt8}(sz_max:0,sz_curr:1))

julia> String(take_up_to_n!(io, 10))
"[1, 2, 3, "

julia> close(io)

julia> String(take_up_to_n!(io, 10))
ERROR: ArgumentError: read failed, WriteBlockingIO is already closed.
```
"""
struct WriteBlockingIO <: IO
    ch :: Channel{UInt8}
    WriteBlockingIO(buf=0) = new(Channel{UInt8}(buf))
    function WriteBlockingIO(f::Function, buf=0; spawn=false)
        new(Channel{UInt8}(ch -> f(new(ch)), buf; spawn=spawn))
    end
end

Base.write(io::WriteBlockingIO, c::UInt8) = put!(io.ch, c)
Base.close(io::WriteBlockingIO) = close(io.ch)
Base.isopen(io::WriteBlockingIO) = isopen(io.ch)

Base.take!(io::WriteBlockingIO) = collect(io.ch)
function take_up_to_n!(io::WriteBlockingIO, n)
    if !isopen(io.ch) && !isready(io.ch)
        throw(ArgumentError("read failed, WriteBlockingIO is already closed."))
    end
    out = UInt8[]
    for i in 1:n
        if !isopen(io.ch) && !isready(io.ch) break end
        try
            push!(out, take!(io.ch))
        catch e
            if e isa InvalidStateException
                break
            else
                rethrow()
            end
        end
    end
    out
end

end # module IOUtils


text:
using Kannaadi
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

@test எழுத்துவகை_பிரி('ஆ') == :உயிர்
@test எழுத்துவகை_பிரி('ம') == :உயிர்மெய்
@test எழுத்துவகை_பிரி("மை") == :உயிர்மெய்
@test எழுத்துவகை_பிரி("ப்") == :மெய்
@test எழுத்துவகை_பிரி('ௐ') == :எழுத்தல்ல
@test_throws ArgumentError எழுத்துவகை_பிரி('x')
@test_throws ArgumentError எழுத்துவகை_பிரி("blah")
@test_throws ArgumentError எழுத்துவகை_பிரி("நிரை")

@test மாத்திரை_அளவெடு('ஆ') == 2
@test மாத்திரை_அளவெடு('ம') == 1
@test மாத்திரை_அளவெடு("மை") == 2
@test மாத்திரை_அளவெடு("ப்") == 0
@test_throws ArgumentError மாத்திரை_அளவெடு('ௐ')
@test_throws ArgumentError மாத்திரை_அளவெடு('x')
@test_throws ArgumentError மாத்திரை_அளவெடு("blah")
@test_throws ArgumentError மாத்திரை_அளவெடு("நிரை")

@test துப்புரவு_செய்("தொகுப்பு (Package)") == "தொகுப்பு "
@test துப்புரவு_செய்("தமிழ் and English (2) உண்டு") == "தமிழ்    உண்டு"

@test தமிழுரையா("வணக்கம்") == true
@test தமிழுரையா("வணக்கம் தமிழகம்") == false
@test தமிழுரையா("vanakkam") == false


text:
using Parameters

abstract type AbstractRateUncertainty end
export AbstractRateUncertainty

struct EmptyRateUncertainty <: AbstractRateUncertainty end
export EmptyRateUncertainty

@with_kw struct ConstantRateUncertainty <: AbstractRateUncertainty
    u::AbstractFloat = nothing
end
(c::ConstantRateUncertainty)(T::AbstractFloat=nothing) = c.u
export ConstantRateUncertainty

@with_kw struct EnergyRateUncertainty <: AbstractRateUncertainty
    dE::AbstractFloat=nothing
end
(er::EnergyRateUncertainty)(T::AbstractFloat=nothing) = exp(dE/(8.314*T))
export EnergyRateUncertainty


text:
function ϕ(A, v, p)
    m = LinearAlgebra.checksquare(A)
    length(v) == m || throw(DimensionMismatch("second dimension of A, $m, does not match length of x, $(length(v))"))
    p == 0 && return exp(A)*v
    A′ = fill!(similar(A, m+p, m+p), 0)
    copyto!(view(A′, 1:m, 1:m), A)
    copyto!(view(A′, 1:m, m+1), v)
    for k = 1:p-1
        A′[m+k, m+k+1] = 1
    end
    return exp(A′)[1:m, end]
end

@testset "Lanczos - expintegrator full" begin
    @testset for T in (Float32, Float64, ComplexF32, ComplexF64)
        @testset for orth in (cgs2, mgs2, cgsr, mgsr)
            A = rand(T,(n,n)) .- one(T)/2
            A = (A+A')/2
            V = one(A)
            W = zero(A)
            alg = Lanczos(orth = orth, krylovdim = n, maxiter = 2, tol = 10*n*eps(real(T)))
            for k = 1:n
                W[:,k], =  @inferred exponentiate(A, 1, view(V,:,k), alg)
            end
            @test W ≈ exp(A)

            pmax = 5
            for t in (rand(real(T)), -rand(real(T)), im*randn(real(T)), randn(real(T))+im*randn(real(T)))
                for p = 1:pmax
                    u = ntuple(i->rand(T, n), p+1)
                    w, info = @inferred expintegrator(A, t, u, alg)
                    w2 = exp(t*A)*u[1]
                    for j = 1:p
                        w2 .+= t^j*ϕ(t*A, u[j+1], j)
                    end
                    @test info.converged > 0
                    @test w2 ≈ w
                end
            end
        end
    end
end

@testset "Arnoldi - expintegrator full" begin
    @testset for T in (Float32, Float64, ComplexF32, ComplexF64)
        @testset for orth in (cgs2, mgs2, cgsr, mgsr)
            A = rand(T,(n,n)) .- one(T)/2
            V = one(A)
            W = zero(A)
            alg = Arnoldi(orth = orth, krylovdim = n, maxiter = 2, tol = 10*n*eps(real(T)))
            for k = 1:n
                W[:,k], =  @inferred exponentiate(A, 1, view(V,:,k), alg)
            end
            @test W ≈ exp(A)

            pmax = 5
            for t in (rand(real(T)), -rand(real(T)), im*randn(real(T)), randn(real(T))+im*randn(real(T)))
                for p = 1:pmax
                    u = ntuple(i->rand(T, n), p+1)
                    w, info = @inferred expintegrator(A, t, u, alg)
                    w2 = exp(t*A)*u[1]
                    for j = 1:p
                        w2 .+= t^j*ϕ(t*A, u[j+1], j)
                    end
                    @test info.converged > 0
                    @test w2 ≈ w
                end
            end
        end
    end
end

@testset "Lanczos - expintegrator iteratively" begin
    @testset for T in (Float32, Float64, ComplexF32, ComplexF64)
        @testset for orth in (cgs2, mgs2, cgsr, mgsr)
            A = rand(T,(N,N)) .- one(T)/2
            A = (A+A')/2
            s = norm(eigvals(A), 1)
            rmul!(A, 1/(10*s))
            pmax = 5
            for t in (rand(real(T)), -rand(real(T)), im*randn(real(T)), randn(real(T))+im*randn(real(T)))
                for p = 1:pmax
                    u = ntuple(i->rand(T, N), p+1)
                    w1, info = @inferred expintegrator(A, t, u...; maxiter = 100, krylovdim = n)
                    @assert info.converged > 0
                    w2 = exp(t*A)*u[1]
                    for j = 1:p
                        w2 .+= t^j*ϕ(t*A, u[j+1], j)
                    end
                    @test w2 ≈ w1
                    w1, info = @inferred expintegrator(A, t, u...; maxiter = 100, krylovdim = n, tol = 1e-3)
                    @test norm(w1 - w2) < 1e-2*abs(t)
                end
            end
        end
    end
end

@testset "Arnoldi - expintegrator iteratively" begin
    @testset for T in (Float32, Float64, ComplexF32, ComplexF64)
        @testset for orth in (cgs2, mgs2, cgsr, mgsr)
            A = rand(T,(N,N)) .- one(T)/2
            s = norm(eigvals(A), 1)
            rmul!(A, 1/(10*s))
            pmax = 5
            for t in (rand(real(T)), -rand(real(T)), im*randn(real(T)), randn(real(T))+im*randn(real(T)))
                for p = 1:pmax
                    u = ntuple(i->rand(T, N), p+1)
                    w1, info = @inferred expintegrator(A, t, u...; maxiter = 100, krylovdim = n)
                    @test info.converged > 0
                    w2 = exp(t*A)*u[1]
                    for j = 1:p
                        w2 .+= t^j*ϕ(t*A, u[j+1], j)
                    end
                    @test w2 ≈ w1
                    w1, info = @inferred expintegrator(A, t, u...; maxiter = 100, krylovdim = n, tol = 1e-3)
                    @test norm(w1 - w2) < 1e-2*abs(t)
                end
            end
        end
    end
end

@testset "Arnoldi - expintegrator fixed point branch" begin
    @testset for T in (ComplexF32, ComplexF64) # less probable that :LR eig is degenerate
        A = rand(T,(N,N))
        v₀ = rand(T, N)
        λs, vs, infoR = eigsolve(A, v₀, 1, :LR)
        @test infoR.converged > 0
        r = vs[1]
        A = A - λs[1]*I
        λs, vs, infoL = eigsolve(A', v₀, 1, :LR)
        @test infoL.converged > 0
        l = vs[1]
        w1, info1 = expintegrator(A, 1000., v₀)
        @test info1.converged > 0
        @test abs(dot(r, w1))/norm(r)/norm(w1) ≈ 1 atol=1e-4
        v₁ = rand(T, N)
        v₁ -= r*dot(l,v₁)/dot(l,r)
        w2, info2 = expintegrator(A, 1000., v₀, v₁)
        @test info2.converged > 0
        @test A*w2 ≈ -v₁
    end
end


text:
using CairoMakie

@testset "series_sol" begin
    @named m = samplesystem()
    @named me = embed_system(m)
    prob = ODEProblem(me, [m.x => 1.1], (0.0,1.0), [])
    sol = solve(prob,Tsit5())
    fig = Figure(); ax = CairoMakie.Axis(fig[1,1]; xlabel = "time (yr)");
    series_sol!(ax, sol, [m.x, m.RHS])
    #axislegend(ax, unique=true, position=:lb)
    fig[1,2] = Legend(fig, ax, unique=true)
    #display(fig)
    @test true
end;


text:
import Base: ∘
export ∘, select, stretch, periodic, shift

"""
    ∘(f::GP, g)

Constructs the GP f′ given by f′(x) := f(g(x))
"""
∘(f::AbstractGP, g) = CompositeGP((∘, f, g), f.gpc)

const comp_args = Tuple{typeof(∘), AbstractGP, Any}

mean_vector((_, f, g)::comp_args, x::AV) = mean_vector(f, g.(x))

cov((_, f, g)::comp_args, x::AV) = cov(f, g.(x))
cov_diag((_, f, g)::comp_args, x::AV) = cov_diag(f, g.(x))

cov((_, f, g)::comp_args, x::AV, x′::AV) = cov(f, g.(x), g.(x′))
cov_diag((_, f, g)::comp_args, x::AV, x′::AV) = cov_diag(f, g.(x), g.(x′))

cov((_, f, g)::comp_args, f′::AbstractGP, x::AV, x′::AV) = cov(f, f′, g.(x), x′)
cov(f::AbstractGP, (_, f′, g)::comp_args, x::AV, x′::AV) = cov(f, f′, x, g.(x′))

cov_diag((_, f, g)::comp_args, f′::AbstractGP, x::AV, x′::AV) = cov_diag(f, f′, g.(x), x′)
cov_diag(f::AbstractGP, (_, f′, g)::comp_args, x::AV, x′::AV) = cov_diag(f, f′, x, g.(x′))

"""
    Stretch{T<:Union{Real, AbstractMatrix{<:Real}}}

Stretch all elements of the inputs by `l`.
"""
struct Stretch{T<:Union{Real, AbstractMatrix{<:Real}}}
    l::T
end
(s::Stretch)(x) = s.l * x
broadcasted(s::Stretch{<:Real}, x::StepRangeLen) = s.l .* x
broadcasted(s::Stretch{<:Real}, x::ColVecs) = ColVecs(s.l .* x.X)
broadcasted(s::Stretch{<:AbstractMatrix}, x::ColVecs) = ColVecs(s.l * x.X)

"""
    stretch(f::AbstractGP, l::Union{AbstractVecOrMat, Real})

Equivalent to `f ∘ Stretch(l)`
"""
stretch(f::AbstractGP, l::Real) = f ∘ Stretch(l)
stretch(f::AbstractGP, a::AbstractVector) = stretch(f, Diagonal(a))
stretch(f::AbstractGP, A::AbstractMatrix) = f ∘ Stretch(A)

"""
    Select{Tidx}

Use inside an input transformation meanfunction / crosskernel to improve peformance.
"""
struct Select{Tidx}
    idx::Tidx
end
(f::Select)(x) = x[f.idx]
broadcasted(f::Select, x::ColVecs) = ColVecs(x.X[f.idx, :])
broadcasted(f::Select{<:Integer}, x::ColVecs) = x.X[f.idx, :]

function broadcasted(f::Select, x::AbstractVector{<:CartesianIndex})
    out = Matrix{Int}(undef, length(f.idx), length(x))
    for i in f.idx, n in eachindex(x)
        out[i, n] = x[n][i]
    end
    return ColVecs(out)
end
@adjoint function broadcasted(f::Select, x::AV{<:CartesianIndex})
    return broadcasted(f, x), Δ->(nothing, nothing)
end
function broadcasted(f::Select{<:Integer}, x::AV{<:CartesianIndex})
    out = Matrix{Int}(undef, length(x))
    for n in eachindex(x)
        out[n] = x[n][f.idx]
    end
    return out
end

"""
    select(f::AbstractGP, idx)

Select the dimensions of the input to `f` given by `idx`.
"""
select(f::AbstractGP, idx) = f ∘ Select(idx)

"""
    Periodic{Tf<:Real}

Make a kernel or mean function periodic by projecting into two dimensions.
"""
struct Periodic{Tf<:Real}
    f::Tf
end
(p::Periodic)(t::Real) = [cos((2π * p.f) * t), sin((2π * p.f) * t)]
function broadcasted(p::Periodic, x::AbstractVector{<:Real})
    return ColVecs(vcat(cos.((2π * p.f) .* x)', sin.((2π * p.f) .* x)'))
end

"""
    periodic(g::AbstractGP, f::Real)

Produce an AbstractGP with period `f`.
"""
periodic(g::AbstractGP, f::Real) = g ∘ Periodic(f)

#
# Translations of GPs through their input spaces.
#

struct Shift{Ta<:Union{Real, AV{<:Real}}}
    a::Ta
end
(f::Shift{<:Real})(x::Real) = x - f.a

broadcasted(f::Shift, x::ColVecs) = ColVecs(x.X .- f.a)

"""
    shift(f::AbstractGP, a::Real)
    shift(f::AbstractGP, a::AbstractVector{<:Real})

Returns the GP `g` given by `g(x) = f(x - a)`
"""
shift(f::AbstractGP, a) = f ∘ Shift(a)


text:
using mPulseAPI
using Base.Test

# Check environment
if !haskey(ENV, "mPulseAPIToken")
    println("Set the `mPulseAPIToken' environment variable to run tests.  See https://soasta.github.io/mPulseAPI.jl/apiToken/index.html for details on how to get a token.")
    exit(-2)
end

if !haskey(ENV, "mPulseAPITenant")
    println("Set the `mPulseAPITenant' environment variable to run tests against your specific mPulse tenant.")
    exit(-3)
end

mPulseAPIToken  = ENV["mPulseAPIToken"]
mPulseAPITenant = ENV["mPulseAPITenant"]
mPulseAPIAlert  = get(ENV, "mPulseAPIAlert", "")

DA_mPulseAPIAlert  = get(ENV, "DA_mPulseAPIAlert", "")

verbosity = (get(ENV, "mPulseAPIVerbose", "false") == "true")

endpoint = get(ENV, "mPulseAPIEndpoint", "")

if !isempty(endpoint)
    mPulseAPI.setEndpoints(endpoint)
end

mPulseAPI.setVerbose(verbosity)

include("repository-tests.jl")

include("query-tests.jl")

include("zzz_change-url-tests.jl")


text:
"""
FinEtools (C) 2017-2020, Petr Krysl

Finite Element tools.  Julia implementation  of the finite element method
for continuum mechanics. Package for problems of vibration in fluids.
"""
module FinEtoolsVibInFluids

__precompile__(true)

include("LaplBEM.jl")
include("AlgoVibInFluidsModule.jl")

end # module

max_stars_repo_name: PetrKryslUCSD/FinEtoolsVibInFluids

max_issues_repo_name: PetrKryslUCSD/FinEtoolsVibInFluids

max_forks_repo_name: PetrKryslUCSD/FinEtoolsVibInFluids


text:
#############################################
# Create some temporary files & directories #
#############################################
dir = mktempdir()
file = joinpath(dir, "afile.txt")
close(open(file,"w")) # like touch, but lets the operating system update the timestamp for greater precision on some platforms (windows)

#######################################################################
# This section tests some of the features of the stat-based file info #
#######################################################################
@test isdir(dir)
@test !isfile(dir)
@test !islink(dir)
@test !isdir(file)
@test isfile(file)
@test !islink(file)
@test isreadable(file)
@test iswritable(file)
# Here's something else that might be UNIX-specific?
run(`chmod -w $file`)
@test !iswritable(file)
run(`chmod +w $file`)
@test !isexecutable(file)
@test filesize(file) == 0
# On windows the filesize of a folder is the accumulation of all the contained
# files and is thus zero in this case.
@windows_only begin
    @test filesize(dir) == 0
end
@unix_only begin
    @test filesize(dir) > 0
end
@test int(time()) >= int(mtime(file)) >= int(mtime(dir)) >= 0 # 1 second accuracy should be sufficient

# rename file
newfile = joinpath(dir, "bfile.txt")
mv(file, newfile)
@test !ispath(file)
@test isfile(newfile)
file = newfile

# Test renaming directories
a_tmpdir = mktempdir()
b_tmpdir = joinpath(dir, "b_tmpdir")

# grab a_tmpdir's file info before renaming
a_stat = stat(a_tmpdir)

# rename, then make sure b_tmpdir does exist and a_tmpdir doesn't
mv(a_tmpdir, b_tmpdir)
@test isdir(b_tmpdir)
@test !ispath(a_tmpdir)

# get b_tmpdir's file info and compare with a_tmpdir
b_stat = stat(b_tmpdir)
@test Base.samefile(a_stat, b_stat)

rmdir(b_tmpdir)

#######################################################################
# This section tests file watchers.                                   #
#######################################################################
function test_file_poll(channel,timeout_s)
    rc = poll_file(file, iround(timeout_s/10), timeout_s)
    put!(channel,rc)
end

function test_timeout(tval)
    tic()
    channel = RemoteRef()
    @async test_file_poll(channel,tval)
    tr = take!(channel)
    t_elapsed = toq()
    @test !tr
    @test tval <= t_elapsed
end

function test_touch(slval)
    tval = slval*1.1
    channel = RemoteRef()
    @async test_file_poll(channel, tval)
    sleep(tval/10)  # ~ one poll period
    f = open(file,"a")
    write(f,"Hello World\n")
    close(f)
    tr = take!(channel)
    @test tr
end

function test_monitor(slval)
    FsMonitorPassed = false
    fm = FileMonitor(file) do args...
        FsMonitorPassed = true
    end
    sleep(slval/2)
    f = open(file,"a")
    write(f,"Hello World\n")
    close(f)
    sleep(slval)
    @test FsMonitorPassed
    close(fm)
end

function test_monitor_wait(tval)
    fm = watch_file(file)
    @async begin
        sleep(tval)
        f = open(file,"a")
        write(f,"Hello World\n")
        close(f)
    end
    fname, events = wait(fm)
    @test fname == basename(file)
    @test events.changed
end

# Commented out the tests below due to issues 3015, 3016 and 3020
test_timeout(0.1)
test_timeout(1)
# the 0.1 second tests are too optimistic
#test_touch(0.1)
test_touch(2)
#test_monitor(0.1)
test_monitor(2)
test_monitor_wait(0.1)

##########
#  mmap  #
##########

s = open(file, "w")
write(s, "Hello World\n")
close(s)
s = open(file, "r")
@test isreadonly(s) == true
c = mmap_array(Uint8, (11,), s)
@test c == "Hello World".data
close(s)
s = open(file, "r+")
@test isreadonly(s) == false
c = mmap_array(Uint8, (11,), s)
c[5] = uint8('x')
msync(c)
close(s)
s = open(file, "r")
str = readline(s)
close(s)
@test beginswith(str, "Hellx World")
c=nothing; gc(); gc(); # cause munmap finalizer to run & free resources

s = open(file, "w")
write(s, [0xffffffffffffffff,
          0xffffffffffffffff,
          0xffffffffffffffff,
          0x000000001fffffff])
close(s)
s = open(file, "r")
@test isreadonly(s)
b = mmap_bitarray((17,13), s)
@test b == trues(17,13)
@test_throws ErrorException mmap_bitarray((7,3), s)
close(s)
s = open(file, "r+")
b = mmap_bitarray((17,19), s)
rand!(b)
msync(b)
b0 = copy(b)
close(s)
s = open(file, "r")
@test isreadonly(s)
b = mmap_bitarray((17,19), s)
@test b == b0
close(s)
b=nothing; b0=nothing; gc(); gc(); # cause munmap finalizer to run & free resources

#######################################################################
# This section tests temporary file and directory creation.           #
#######################################################################

# my_tempdir = tempdir()
# @test isdir(my_tempdir) == true

# path = tempname()
# @test ispath(path) == false

# (file, f) = mktemp()
# print(f, "Here is some text")
# close(f)
# @test isfile(file) == true
# @test readall(file) == "Here is some text"

emptyfile = joinpath(dir, "empty")
touch(emptyfile)
emptyf = open(emptyfile)
@test isempty(readlines(emptyf))
close(emptyf)
rm(emptyfile)

# Test copy file
afile = joinpath(dir, "a.txt")
touch(afile)
af = open(afile, "r+")
write(af, "This is indeed a test")

bfile = joinpath(dir, "b.txt")
cp(afile, bfile)

a_stat = stat(afile)
b_stat = stat(bfile)
@test a_stat.mode == b_stat.mode
@test a_stat.size == b_stat.size

close(af)
rm(afile)
rm(bfile)

###################
# FILE* interface #
###################

f = open(file, "w")
write(f, "Hello, world!")
close(f)
f = open(file, "r")
FILEp = convert(CFILE, f)
buf = Array(Uint8, 8)
str = ccall(:fread, Csize_t, (Ptr{Void}, Csize_t, Csize_t, Ptr{Void}), buf, 1, 8, FILEp.ptr)
@test bytestring(buf) == "Hello, w"
@test position(FILEp) == 8
seek(FILEp, 5)
@test position(FILEp) == 5
close(f)

############
# Clean up #
############
rm(file)
rmdir(dir)
@test !ispath(file)
@test !ispath(dir)

max_stars_repo_name: tobig/julia

max_issues_repo_name: tobig/julia

max_forks_repo_name: tobig/julia


text:
# -----------------------------------------------------------------------------
# Diagonally-dominant matrices [JuMP]
# -----------------------------------------------------------------------------

# Ahmadi, A., & Majumdar, A. (2017). Dsos and sdsos optimization: More
# tractable alternatives to sum of squares and semidefinite optimization.
# SIAM Journal on Applied Algebra and Geometry, 3. https://doi.org/
# 10.1137/18M118935X

# JuMP constraints for diagonal-dominance -------------------------------------

"""
Constraint for weak column diagonal dominance with nonnegative diagonal entries.
A matrix M is diagonally-dominant, or M ∈ dd, if:
        |M[i,i]| ≥ ∑_{i≂̸j} |M[i,j]|, ∀i.
However, for our problem, dd also imposes nonnegative diagonal entries.
"""
struct dd end
struct ddConstraint <: AbstractConstraint
        mat::AbstractArray
end

JuMP.build_constraint(::Function, m::AbstractArray, ::dd) = ddConstraint(m)

function JuMP.add_constraint(model::Model, c::ddConstraint, ::String = "")
        X = c.mat
        n = size(X, 1)
        W = X - Diagonal(X)
        Z = @variable(model, [1:n, 1:n], Symmetric)
        @constraint(model, W .>= -Z)
        @constraint(model, W .<= Z)
        @constraint(model, diag(X) .>= sum(Z, dims = 2))
end

# JuMP constraints for the DD family of cones ---------------------------------
# DD(U) := {M ∈ Symmetric | M = U'QU for some dd matrix Q}.
# One can think of DD(U) as the set of matrices that are dd after a change of
# coordinates via the matrix U.
struct DD
        u::AbstractArray
end
struct DDConstraint <: AbstractConstraint
        mat::AbstractArray
        u::AbstractArray
end

JuMP.build_constraint(::Function, m::AbstractArray, d::DD) =
        DDConstraint(m, d.u)

function JuMP.add_constraint(model::Model, c::DDConstraint, ::String = "")
        n = size(c.mat, 1)
        Q = @variable(model, [1:n, 1:n], Symmetric)
        @constraint(model, Q in dd())
        xs = JuMP.vectorize(c.mat - c.u' * Q * c.u, SymmetricMatrixShape(n))
	@constraint(model, xs in MOI.Zeros(length(xs)))
end


text:
"""
Read in inference configuration.
"""
function read_inference_configuration(data::Dict{String, Any})
    ic = InferenceConfiguration()
	log_message("created inference config object")
    if haskey(data, "inference_configuration")
        ic_data = data["inference_configuration"]
        if haskey(ic_data, "model") ic.model_config = Meta.eval(Symbol(ic_data["model"])) end
        if haskey(ic_data, "convergence_criterion") ic.convergence_criterion = ic_data["convergence_criterion"] end
        if haskey(ic_data, "min_iter") ic.min_iter = ic_data["min_iter"] end
        if haskey(ic_data, "max_iter") ic.max_iter = ic_data["max_iter"] end
        if haskey(ic_data, "do_gridsearch") ic.do_gridsearch = ic_data["do_gridsearch"] end
        if haskey(ic_data, "do_cross_validation") ic.do_cross_validation = ic_data["do_cross_validation"] end
        if haskey(ic_data, "compute_wpc_index") ic.compute_wpc_index = ic_data["compute_wpc_index"] end
        if haskey(ic_data, "do_variance_filtering") ic.do_variance_filtering = ic_data["do_variance_filtering"] end
        if haskey(ic_data, "subsume_pathways") ic.subsume_pathways = ic_data["subsume_pathways"] end
    end
    ic
end

function read_model_type(data::Dict{String, Any})
	log_message("in read_model_type")
	if haskey(data, "model")
		try
			model_type = eval(Meta.parse(data["model"]))
			log_message("model type is $model_type")
			return model_type
		catch
			log_message("unknown model type: $(mc_data["model"])")
			throw(ArgumentError("Unknown model: $(mc_data["model"])"))
		end
    end
	log_message("no model specified in $data")
    throw(ArgumentError("You must specify a model."))
end

function read_model_configuration(pm::PredictionModel, data::Dict{String, Any})
	log_message("read_model_configuration not implemented for $(typeof(pm)), skipping...")
end

function init!(pm::PredictionModel, ic::InferenceConfiguration)
	log_message("init! method not implemented for $(typeof(pm)), skipping...")
end

function update_model_config!(pm::PredictionModel, mc::ModelConfiguration)
	log_message("update_model_config! not implemented for $(typeof(pm)), skipping...")
end

function post_init!(pm::PredictionModel)
	log_message("post_init! method not implemented, skipping...")
end

function train(experiment_id::String, data::Dict{String, Any})
	try
	    (pm, inference_config, model_config) = create_prediction_model(experiment_id, data)
	    log_progress(experiment_id, "created prediction model")

	    # actual training procedure

	    # create a temp dir for storing the results
        temp_dir = mktempdir()
	    inference_config.target_dir = joinpath(temp_dir, "results")

	    log_progress(experiment_id, "running model computations")
	    # do cross validation or single model training
	    if inference_config.do_cross_validation
	        run_cross_validation(pm, inference_config, model_config)
	    else
	        run_model(pm, inference_config, model_config)
	    end

	    log_progress(experiment_id, "preparing results for download")
	    # compile the results into a tar ball
	    tar_dir = mktempdir()
	    tar_file = joinpath(tar_dir, "results.tgz")
	    log_message("archiving results at $tar_file")
	    tar_results_cmd = `tar czf $(tar_file) -C $temp_dir results`
	    run(tar_results_cmd)
	    # remove the output files
	    rm(temp_dir, force = true, recursive = true)
	    Utils.store_result(experiment_id, tar_file)
	    log_progress(experiment_id, "finished, results ready for download")
	catch ex
		exception_name = string(ex)
		Core.println(exception_name)
        st = map(string, stacktrace(catch_backtrace()))
		log_progress(experiment_id, exception_name)
		for ste in st
			Core.println(ste)
			log_progress(experiment_id, ste)
		end
	end
end

function predict(experiment_id::String, data::Dict{String, Any})
    (pm, inference_config, model_config) = create_prediction_model(experiment_id, data)

    if model_config == nothing
        throw(ArgumentError("Prediction not implemented for gridsearch"))
    end

    log_progress(experiment_id, "training model")
    # do one inference run on the model with provided parameters
    # TODO: do the training via additional api call, then store the most recent paramters object
    (lls, errs, test_errs, params, convergence) = inference(pm, model_config, inference_config = inference_config)

    # read cell lines from data for which we want a prediction
    log_progress(experiment_id, "predicting efficacy on provided cell lines")
    prediction_cell_lines = create_cell_lines(experiment, data, for_prediction = true)
    (predictions, ranks) = predict_outcomes(pm, params, prediction_cell_lines, held_out = true)

    # create a result object
    # TODO: do this asynchronously and download results?
    result_dictionary = Dict{String, Any}()
    for (t, drug) in enumerate(keys(pm.data.results))
        result_dictionary[drug.id] = Dict{String, Any}()
        for (cl_idx, cl) in enumerate(prediction_cell_lines)
            result_dictionary[drug.id][cl.id] = Dict{String, Any}()
            result_dictionary[drug.id][cl.id]["prediction"] = predictions[drug][cl_idx]
            result_dictionary[drug.id][cl.id]["rank"] = ranks[drug][cl_idx]
        end
    end
    log_progress(experiment, "finished")
    result_dictionary
end


text:
@testset "ancestors" begin
    mtg = read_mtg("files/simple_plant.mtg")
    width_all = [nothing, nothing, nothing, 1.0, 6.0, nothing, 7.0]

    # Using a leaf node from the mtg:
    leaf_node = mtg.children[2].children[3].children[4].children[5]

    @test ancestors(leaf_node, :Width; type = Union{Nothing,Float64}) == reverse(width_all[1:4])
    @test ancestors(leaf_node, :Width) == reverse(width_all[1:4])

    d = ancestors(leaf_node, :Width, scale = 3)
    @test typeof(d) == Vector{Any}
    @test length(d) == 1
    @test d[1] == width_all[4]
    d_typed = ancestors(leaf_node, :Width, type = Union{Nothing,Float64})
    @test typeof(d_typed) == Vector{Union{Nothing,Float64}}
    @test ancestors(leaf_node, :Width, symbol = ("Leaf", "Internode")) == width_all[[4]]

    @test ancestors(leaf_node, :Width, symbol = ("Leaf", "Internode"), self = true) ==
          width_all[5:-1:4]
end


text:
using Documenter
using VortexFilaments

ENV["GKSwstype"] = "nul"

makedocs(
    modules = [VortexFilaments],
    sitename = "VortexFilaments.jl",
    doctest = true,
    clean = true,
    pages = [
        "Home" => "index.md",
        "Library" => ["public.md",
                      "private.md"]
    ],
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true",
        mathengine = MathJax2(Dict(
            :TeX => Dict(
                :equationNumbers => Dict(:autoNumber => "AMS"),
                :Macros => Dict()
            )
        ))
    )
)

deploydocs(
    repo = "github.com/diederikb/VortexFilaments.jl.git",
    target = "build"
)