text:
@testset "Change bounds" begin
    cp = test_LP()

    change_bound!(cp, 1, lower = -10, upper = 10)
    @test cp.xl[1] == -10
    @test cp.xu[1] == 10
    change_bounds!(cp, [1, 2]; lower = [-11, -12.2], upper = [11, 23.0])
    @test cp.xl[2] == -12.2
    @test cp.xu[1] == 11

    change_bound!(cp, "r1", lower = -101, upper = 101)
    @test cp.xl[1] == -101
    @test cp.xu[1] == 101
    change_bounds!(cp, ["r1", "r2"]; lower = [-113, -12.23], upper = [114, 233.0])
    @test cp.xl[2] == -12.23
    @test cp.xu[1] == 114
    change_bounds!(cp, ["r1", "r2"]; lower = [-114, nothing], upper = [nothing, 2333.0])
    @test cp.xl[1] == -114
    @test cp.xl[2] == -12.23
    @test cp.xu[1] == 114
    @test cp.xu[2] == 2333

    new_model = change_bound(cp, 1, lower = -10, upper = 10)
    @test new_model.xl[1] == -10
    @test new_model.xu[1] == 10
    new_model = change_bounds(cp, [1, 2]; lower = [-11, -12.2], upper = [11, 23.0])
    @test new_model.xl[2] == -12.2
    @test new_model.xu[1] == 11

    new_model = change_bound(cp, "r1", lower = -101, upper = 101)
    @test new_model.xl[1] == -101
    @test new_model.xu[1] == 101
    new_model = change_bounds(cp, ["r1", "r2"]; lower = [-113, -12.23], upper = [113, 1000])
    @test new_model.xl[2] == -12.23
    @test new_model.xu[1] == 113
    new_model =
        change_bounds(cp, ["r1", "r2"]; lower = [nothing, -10], upper = [110, nothing])
    @test new_model.xl[1] == -114
    @test new_model.xl[2] == -10
    @test new_model.xu[1] == 110
    @test new_model.xu[2] == 2333
end

@testset "Verify consistency" begin
    cp = test_LP()
    (new_reactions, new_mets) = verify_consistency(
        cp,
        reshape(cp.S[:, end], :, 1),
        [1.0, 2.0, 3.0, 4.0],
        [2.0],
        [-1.0],
        [1.0],
        ["r4"],
        ["m1", "m2", "m3", "m6"],
        [1],
        [4],
    )
    @test new_reactions == [1]
    @test new_mets == [4]

    (new_reactions, new_mets) = verify_consistency(
        cp,
        reshape(cp.S[:, end], :, 1),
        [1.0, 2.0, 3.0, 4.0],
        [2.0],
        [-1.0],
        [1.0],
        ["r1"],
        ["m1", "m2", "m3", "m6"],
        [],
        [4],
    )
    @test new_reactions == []
    @test new_mets == [4]
end

@testset "Add reactions (checking existence and consistency)" begin
    cp = test_LP()
    (new_cp, new_reactions, new_mets) = add_reactions(
        cp,
        cp.S[:, end],
        [1.0, 2.0, 3.0, 4.0],
        2.0,
        -1.0,
        1.0,
        check_consistency = true,
    )
    @test n_reactions(cp) + 1 == n_reactions(new_cp)

    (new_cp, new_reactions, new_mets) = add_reactions(
        cp,
        cp.S[:, end],
        [1.0, 2.0, 3.0, 4.0],
        2.0,
        -1.0,
        1.0,
        "r1",
        ["m1", "m2", "m3", "m6"],
        check_consistency = true,
    )
    @test n_reactions(cp) == n_reactions(new_cp)
    @test n_metabolites(cp) + 1 == n_metabolites(new_cp)
end

@testset "Add reactions" begin
    cp = test_LP()
    cp = add_reactions(cp, 2.0 * ones(4), 3 .* ones(4), 2.0, -1.0, 1.0)
    @test size(cp.S) == (8, 4)
    cp = add_reactions(cp, 2.0 * ones(4, 1), 3 .* ones(4), 2 .* ones(1), -ones(1), ones(1))
    @test size(cp.S) == (12, 5)
    cp = add_reactions(
        cp,
        2.0 * ones(4, 10),
        3 .* ones(4),
        2 .* ones(10),
        -ones(10),
        ones(10),
    )
    @test size(cp.S) == (16, 15)

    cp = test_sparseLP()
    @test size(cp.S) == (4000, 3000)
    cp = add_reactions(cp, 2.0 * sprand(4000, 0.5), 3 .* sprand(4000, 0.5), 2.0, -1.0, 1.0)
    @test size(cp.S) == (8000, 3001)
    cp = add_reactions(
        cp,
        2.0 * sprand(4000, 1, 0.5),
        3 .* sprand(4000, 0.5),
        2 .* sprand(1, 0.5),
        -sprand(1, 0.5),
        sprand(1, 0.5),
    )
    @test size(cp.S) == (12000, 3002)
    cp = add_reactions(
        cp,
        2.0 * sprand(4000, 1000, 0.5),
        3 .* sprand(4000, 0.5),
        2 .* sprand(1000, 0.5),
        -sprand(1000, 0.5),
        sprand(1000, 0.5),
    )
    @test size(cp.S) == (16000, 4002)

    cp = test_sparseLP()
    @test size(cp.S) == (4000, 3000)
    cp = add_reactions(cp, 2.0 * ones(4000), 3 .* ones(4000), 2.0, -1.0, 1.0)
    @test size(cp.S) == (8000, 3001)
    cp = add_reactions(
        cp,
        2.0 * ones(4000, 1),
        3 .* ones(4000),
        2 .* ones(1),
        -ones(1),
        ones(1),
    )
    @test size(cp.S) == (12000, 3002)
    cp = add_reactions(
        cp,
        2.0 * ones(4000, 1000),
        3 .* ones(4000),
        2 .* ones(1000),
        -ones(1000),
        ones(1000),
    )
    @test size(cp.S) == (16000, 4002)

    # proper subset of existing metabolites
    cp = test_LP()
    new_cp = add_reactions(cp, [-1.0], zeros(1), 1.0, 0.0, 1.0, "r4", ["m1"])
    @test n_reactions(cp) + 1 == n_reactions(new_cp)

    @test_throws DimensionMismatch add_reactions(
        cp,
        2.0 * ones(4000, 1),
        3 .* ones(4000),
        2 .* ones(2),
        -ones(1),
        ones(1),
    )
end

@testset "Remove reactions" begin
    cp = test_LP()
    cp = remove_reaction(cp, 2)
    @test size(cp.S) == (4, 2)
    cp = remove_reactions(cp, [2, 1])
    @test size(cp.S) == (4, 0)

    cp = test_LP()
    cp = remove_reaction(cp, "r1")
    @test size(cp.S) == (4, 2)
    cp = remove_reactions(cp, ["r2"])
    @test size(cp.S) == (4, 1)

    lp = CoreModel(
        [1.0 1 1 0; 1 1 1 0; 1 1 1 0; 0 0 0 1],
        collect(1.0:4),
        collect(1.0:4),
        collect(1.0:4),
        collect(1.0:4),
        ["r1"; "r2"; "r3"; "r4"],
        ["m1"; "m2"; "m3"; "m4"],
    )

    modLp = remove_reactions(lp, [4; 1])
    @test stoichiometry(modLp) == stoichiometry(lp)[:, 2:3]
    @test balance(modLp) == balance(lp)
    @test objective(modLp) == objective(lp)[2:3]
    @test bounds(modLp)[1] == bounds(lp)[1][2:3]
    @test bounds(modLp)[2] == bounds(lp)[2][2:3]
    @test reactions(modLp) == reactions(lp)[2:3]
    @test metabolites(modLp) == metabolites(lp)
end

@testset "Remove metabolites" begin
    model = load_model(CoreModel, model_paths["e_coli_core.json"])

    m1 = remove_metabolites(model, ["glc__D_e", "for_c"])
    m2 = remove_metabolite(model, "glc__D_e")
    m3 = remove_metabolites(model, Int.(indexin(["glc__D_e", "for_c"], metabolites(model))))
    m4 = remove_metabolite(model, first(indexin(["glc__D_e"], metabolites(model))))

    @test size(stoichiometry(m1)) == (70, 90)
    @test size(stoichiometry(m2)) == (71, 93)
    @test size(stoichiometry(m3)) == (70, 90)
    @test size(stoichiometry(m4)) == (71, 93)
    @test all((!in(metabolites(m1))).(["glc__D_e", "for_c"]))
    @test !(["glc__D_e"] in metabolites(m2))
    @test all((!in(metabolites(m3))).(["glc__D_e", "for_c"]))
    @test !(["glc__D_e"] in metabolites(m4))
end

@testset "Core in place modifications" begin
    toymodel = test_toyModel()

    rxn1 = Reaction("nr1"; metabolites = Dict("m1[c]" => -1, "m3[c]" => 1))
    rxn2 = Reaction("nr2"; metabolites = Dict("m1[c]" => -1, "m2[c]" => 1))
    rxn3 = Reaction("nr3"; metabolites = Dict("m2[c]" => -1, "m3[c]" => 1))
    rxn3.lb = 10

    add_reaction!(toymodel, rxn1)
    @test toymodel.S[1, 8] == -1
    @test toymodel.S[2, 8] == 1
    @test all(toymodel.S[3:end, 8] .== 0)
    @test size(toymodel.S) == (6, 8)
    @test toymodel.rxns[end] == "nr1"

    add_reactions!(toymodel, [rxn2, rxn3])
    @test size(toymodel.S) == (6, 10)
    @test toymodel.xl[end] == 10

    change_objective!(toymodel, "nr1")
    @test objective(toymodel)[8] == 1.0
    @test objective(toymodel)[7] == 0.0
end


text:
"""
LTESOM - Long Term Energy System Optimization Model
Copyright (C) 2021, Peking University, College of Engineering (COE),
Department of Industrial Engineering and Management
"""

### Working directory
cd(dirname(@__FILE__))

## Julia Enironment
using Pkg
env_path = "../env/julenv.jl"
# Run this line to activate the Julia virtual environment for GenX;
# skip it, if the appropriate package versions are installed
#include(env_path)
Pkg.activate("LTESOMEnv")
### Set relevant directory paths
src_path = "../src/"
settings_path = joinpath(pwd(), "Settings")
root_path = pwd()

### Load GenX and other modules
println("Loading packages")
push!(LOAD_PATH, src_path)

using GenX
using YAML

### Settings initialization
setup = Dict()
genx_settings = joinpath(settings_path, "genx_settings.yml")
setup_genx = YAML.load(open(genx_settings))
global_settings = joinpath(settings_path, "global_model_settings.yml")
setup_global = YAML.load(open(global_settings))
setup = merge(setup_genx, setup_global)

### Cluster time series inputs if necessary or specified by the user
TDRpath = joinpath(root_path, setup["TimeDomainReductionFolder"])
if setup["TimeDomainReduction"] == 1
    if (!isfile(TDRpath*"/Load_data.csv")) || (!isfile(TDRpath*"/Generators_variability.csv")) || (!isfile(TDRpath*"/Fuels_data.csv"))
        println("Clustering Time Series Data...")
        cluster_inputs(root_path, settings_path, setup)
    else
        println("Time Series Data Already Clustered.")
    end
end

### Configure solver
solver_settings_path = joinpath(settings_path, "solvers")
OPTIMIZER = configure_solver(setup["Solver"], solver_settings_path)

### Running a case
## Load inputs
inputs = load_inputs(setup, root_path)

## Generate model
EP = generate_model(setup, inputs, OPTIMIZER)

## Solve model
EP, solve_time = solve_model(EP, setup)
inputs["solve_time"] = solve_time # Store the model solve time in inputs

## Write power system output
out_path = "$root_path/Results"
out_path = write_outputs(EP, out_path, setup, inputs)

# Run MGA if the MGA flag is set to 1 else only save the least cost solution
# Only valid for power system analysis at this point
if setup["ModelingToGenerateAlternatives"] == 1
    mga(EP,root_path,setup,inputs,out_path)
end

ext: jl

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
@generated function _vrange(::Val{W}, ::Type{T}, ::Val{O}, ::Val{F}) where {W,T,O,F}
  t = Expr(:tuple)
  foreach(w -> push!(t.args, Expr(:call, :(Core.VecElement), T(F * w + O))), 0:W-1)
  Expr(:block, Expr(:meta, :inline), Expr(:call, :Vec, t))
end
@inline function vrange(::Val{W}, ::Type{T}, ::Val{O}, ::Val{F}) where {W,T,O,F}
  _vrange(Val{W}(), pick_integer(Val{W}(), T), Val{O}(), Val{F}())
end

function pick_integer_bytes(
  W::Int,
  preferred::Int,
  sirs::Int,
  minbytes::Int = min(preferred, 4),
)
  # SIMD quadword integer support requires AVX512DQ
  # preferred = AVX512DQ ? preferred :  min(4, preferred)
  max(minbytes, min(preferred, prevpow2(sirs ÷ W)))
end
"""
  vrange(::Val{W}, i::I, ::Val{O}, ::Val{F})

W - Vector width
i::I - dynamic offset
O - static offset
F - static multiplicative factor
"""
@generated function _vrangeincr(
  ::Val{W},
  i::I,
  ::Val{O},
  ::Val{F},
  ::StaticInt{SIRS},
) where {W,I<:Union{Integer,StaticInt},O,F,SIRS}
  isone(W) && return Expr(:block, Expr(:meta, :inline), :(Base.add_int(i, $(O % I))))
  bytes = pick_integer_bytes(W, sizeof(I), SIRS)
  bits = 8bytes
  jtypesym = Symbol(I <: Signed ? :Int : :UInt, bits)
  iexpr = bytes == sizeof(I) ? :i : Expr(:call, :%, :i, jtypesym)
  typ = "i$(bits)"
  vtyp = vtype(W, typ)
  rangevec = join(("$typ $(F*w + O)" for w ∈ 0:W-1), ", ")
  instrs = """
      %ie = insertelement $vtyp undef, $typ %0, i32 0
      %v = shufflevector $vtyp %ie, $vtyp undef, <$W x i32> zeroinitializer
      %res = add nsw $vtyp %v, <$rangevec>
      ret $vtyp %res
  """
  quote
    $(Expr(:meta, :inline))
    Vec($LLVMCALL($instrs, _Vec{$W,$jtypesym}, Tuple{$jtypesym}, $iexpr))
  end
end
@inline function vrangeincr(::Val{W}, i::I, ::Val{O}, ::Val{F}) where {W,I<:Union{Integer,StaticInt},O,F}
  _vrangeincr(Val{W}(), i, Val{O}(), Val{F}(), simd_integer_register_size())
end
@generated function vrangeincr(
  ::Val{W},
  i::T,
  ::Val{O},
  ::Val{F},
) where {W,T<:FloatingTypes,O,F}
  isone(W) && return Expr(:block, Expr(:meta, :inline), :(Base.add_float_fast(i, $(T(O)))))
  typ = LLVM_TYPES[T]
  vtyp = vtype(W, typ)
  rangevec = join(("$typ $(F*w+O).0" for w ∈ 0:W-1), ", ")
  instrs = """
      %ie = insertelement $vtyp undef, $typ %0, i32 0
      %v = shufflevector $vtyp %ie, $vtyp undef, <$W x i32> zeroinitializer
      %res = fadd fast $vtyp %v, <$rangevec>
      ret $vtyp %res
  """
  quote
    $(Expr(:meta, :inline))
    Vec($LLVMCALL($instrs, _Vec{$W,$T}, Tuple{$T}, i))
  end
end
# @generated function vrangemul(::Val{W}, i::I, ::Val{O}, ::Val{F}) where {W,I<:Integer,O,F}
#     isone(W) && return Expr(:block, Expr(:meta,:inline), :(vmul(i, $(O % I))))
#     bytes = pick_integer_bytes(W, sizeof(T))
#     bits = 8bytes
#     jtypesym = Symbol(I <: Signed  ? :Int : :UInt, bits)
#     iexpr = bytes == sizeof(I) ? :i : Expr(:call, :%, :i, jtypesym)
#     typ = "i$(bits)"
#     vtyp = vtype(W, typ)
#     rangevec = join(("$typ $(F*w+O)" for w ∈ 0:W-1), ", ")
#     instrs = """
#         %ie = insertelement $vtyp undef, $typ %0, i32 0
#         %v = shufflevector $vtyp %ie, $vtyp undef, <$W x i32> zeroinitializer
#         %res = mul nsw $vtyp %v, <$rangevec>
#         ret $vtyp %res
#     """
#     quote
#         $(Expr(:meta,:inline))
#         Vec($LLVMCALL(instrs, _Vec{$W,$jtypesym}, Tuple{$jtypesym}, $iexpr))
#     end
# end
# @generated function vrangemul(::Val{W}, i::T, ::Val{O}, ::Val{F}) where {W,T<:FloatingTypes,O,F}
#     isone(W) && return Expr(:block, Expr(:meta,:inline), :(Base.FastMath.mul_fast(i, $(T(O)))))
#     typ = LLVM_TYPES[T]
#     vtyp = vtype(W, typ)
#     rangevec = join(("$typ $(F*w+O).0" for w ∈ 0:W-1), ", ")
#     instrs = """
#         %ie = insertelement $vtyp undef, $typ %0, i32 0
#         %v = shufflevector $vtyp %ie, $vtyp undef, <$W x i32> zeroinitializer
#         %res = fmul fast $vtyp %v, <$rangevec>
#         ret $vtyp %res
#     """
#     quote
#         $(Expr(:meta,:inline))
#         Vec($LLVMCALL(instrs, _Vec{$W,$T}, Tuple{$T}, i))
#     end
# end

@inline Vec(i::MM{W,X}) where {W,X} = vrangeincr(Val{W}(), data(i), Val{0}(), Val{X}())
@inline Vec(i::MM{W,X,StaticInt{N}}) where {W,X,N} =
  vrange(Val{W}(), Int, Val{N}(), Val{X}())
@inline Vec(i::MM{1}) = data(i)
@inline Vec(i::MM{1,<:Any,StaticInt{N}}) where {N} = N
@inline vconvert(::Type{Vec{W,T}}, i::MM{W,X}) where {W,X,T} =
  vrangeincr(Val{W}(), convert(T, data(i)), Val{0}(), Val{X}())
@inline vconvert(::Type{Vec{W,T}}, i::MM{W,X}) where {W,X,T<:IntegerTypesHW} =
  vrangeincr(Val{W}(), data(i) % T, Val{0}(), Val{X}())
@inline vconvert(::Type{T}, i::MM{W,X}) where {W,X,T<:NativeTypes} =
  vrangeincr(Val{W}(), convert(T, data(i)), Val{0}(), Val{X}())

# Addition
@inline vadd_fast(i::MM{W,X}, j::MM{W,Y}) where {W,X,Y} =
  MM{W}(vadd_fast(data(i), data(j)), StaticInt{X}() + StaticInt{Y}())
@inline vadd_fast(i::MM{W}, j::AbstractSIMDVector{W}) where {W} = vadd_fast(Vec(i), j)
@inline vadd_fast(i::AbstractSIMDVector{W}, j::MM{W}) where {W} = vadd_fast(i, Vec(j))
@inline vadd_nsw(i::MM{W,X}, j::MM{W,Y}) where {W,X,Y} =
  MM{W}(vadd_nsw(data(i), data(j)), StaticInt{X}() + StaticInt{Y}())
@inline vadd_nsw(i::MM{W}, j::AbstractSIMDVector{W}) where {W} = vadd_nsw(Vec(i), j)
@inline vadd_nsw(i::AbstractSIMDVector{W}, j::MM{W}) where {W} = vadd_nsw(i, Vec(j))

# @inline vadd(i::MM{W,X}, j::MM{W,Y}) where {W,X,Y} = vadd_fast(i, j)
# @inline vadd(i::MM{W}, j::AbstractSIMDVector{W}) where {W} = vadd_fast(i, j)
# @inline vadd(i::AbstractSIMDVector{W}, j::MM{W}) where {W} = vadd_fast(i, j)

# Subtraction
@inline vsub_fast(i::MM{W,X}, j::MM{W,Y}) where {W,X,Y} =
  MM{W}(vsub_fast(data(i), data(j)), StaticInt{X}() - StaticInt{Y}())
@inline vsub_fast(i::MM{W}, j::AbstractSIMDVector{W}) where {W} = vsub_fast(Vec(i), j)
@inline vsub_fast(i::AbstractSIMDVector{W}, j::MM{W}) where {W} = vsub_fast(i, Vec(j))

@inline vsub_nsw(i::MM{W,X}, j::MM{W,Y}) where {W,X,Y} =
  MM{W}(vsub_nsw(data(i), data(j)), StaticInt{X}() - StaticInt{Y}())
@inline vsub_nsw(i::MM{W}, j::AbstractSIMDVector{W}) where {W} = vsub_nsw(Vec(i), j)
@inline vsub_nsw(i::AbstractSIMDVector{W}, j::MM{W}) where {W} = vsub_nsw(i, Vec(j))
# Multiplication
@inline vmul_fast(i::MM{W}, j::AbstractSIMDVector{W}) where {W} = vmul_fast(Vec(i), j)
@inline vmul_fast(i::AbstractSIMDVector{W}, j::MM{W}) where {W} = vmul_fast(i, Vec(j))
@inline vmul_fast(i::MM{W}, j::MM{W}) where {W} = vmul_fast(Vec(i), Vec(j))
@inline vmul_fast(i::MM, j::IntegerTypesHW) = vmul_fast(Vec(i), j)
@inline vmul_fast(j::IntegerTypesHW, i::MM) = vmul_fast(j, Vec(i))

@inline vmul_nsw(i::MM{W}, j::AbstractSIMDVector{W}) where {W} = vmul_nsw(Vec(i), j)
@inline vmul_nsw(i::AbstractSIMDVector{W}, j::MM{W}) where {W} = vmul_nsw(i, Vec(j))
@inline vmul_nsw(i::MM{W}, j::MM{W}) where {W} = vmul_nsw(Vec(i), Vec(j))
@inline vmul_nsw(i::MM, j::IntegerTypesHW) = vmul_nsw(Vec(i), j)
@inline vmul_nsw(j::IntegerTypesHW, i::MM) = vmul_nsw(j, Vec(i))

# Division
@generated _floattype(::Union{StaticInt{R},Val{R}}) where {R} = R ≥ 8 ? :Float64 : :Float32
@inline floattype(::Val{W}) where {W} = _floattype(register_size() ÷ StaticInt{W}())

@inline vfloat(i::MM{W,X,I}) where {W,X,I} =
  Vec(MM{W,X}(floattype(Val{W}())(getfield(i, :i) % pick_integer(Val{W}(), I))))
@inline vfdiv(i::MM, j::T) where {T<:Real} = float(i) / j
@inline vfdiv(j::T, i::MM) where {T<:Real} = j / float(i)
@inline vfdiv_fast(i::MM, j::MM) = vfdiv_fast(float(i), float(j))
@inline vfdiv_fast(i::MM, j::T) where {T<:Real} = vfdiv_fast(float(i), j)
@inline vfdiv_fast(j::T, i::MM) where {T<:Real} = vfdiv_fast(j, float(i))

@inline vfdiv(i::MM, j::VecUnroll{N,W,T,V}) where {N,W,T,V} = float(i) / j
@inline vfdiv(j::VecUnroll{N,W,T,V}, i::MM) where {N,W,T,V} = j / float(i)

@inline vfdiv(i::MM, j::MM) = float(i) / float(j)
@inline vfdiv(vu::VecUnroll, m::MM) = vu * inv(m)
@inline vfdiv(m::MM, vu::VecUnroll) = Vec(m) / vu

@inline Base.:(<<)(i::MM{W,X,T}, j::StaticInt) where {W,X,T<:IntegerTypes} =
  MM{W}(getfield(i, :i) << j, StaticInt{X}() << j)
@inline Base.:(>>)(i::MM{W,X,T}, j::StaticInt) where {W,X,T<:IntegerTypes} =
  MM{W}(getfield(i, :i) >> j, StaticInt{X}() >> j)
@inline Base.:(>>>)(i::MM{W,X,T}, j::StaticInt) where {W,X,T<:IntegerTypes} =
  MM{W}(getfield(i, :i) >>> j, StaticInt{X}() >>> j)

@inline Base.:(<<)(i::MM{W,X,T}, j::StaticInt) where {W,X,T<:StaticInt} =
  MM{W}(getfield(i, :i) << j, StaticInt{X}() << j)
@inline Base.:(>>)(i::MM{W,X,T}, j::StaticInt) where {W,X,T<:StaticInt} =
  MM{W}(getfield(i, :i) >> j, StaticInt{X}() >> j)
@inline Base.:(>>>)(i::MM{W,X,T}, j::StaticInt) where {W,X,T<:StaticInt} =
  MM{W}(getfield(i, :i) >>> j, StaticInt{X}() >>> j)

# for (f,op) ∈ [
#     (:scalar_less, :(<)), (:scalar_greater,:(>)), (:scalar_greaterequal,:(≥)), (:scalar_lessequal,:(≤)), (:scalar_equal,:(==)), (:scalar_notequal,:(!=))
# ]
#     @eval @inline $f(i::MM, j::Real) = $op(data(i), j)
#     @eval @inline $f(i::Real, j::MM) = $op(i, data(j))
#     @eval @inline $f(i::MM, ::StaticInt{j}) where {j} = $op(data(i), j)
#     @eval @inline $f(::StaticInt{i}, j::MM) where {i} = $op(i, data(j))
#     @eval @inline $f(i::MM, j::MM) = $op(data(i), data(j))
#     @eval @inline $f(i, j) = $op(i, j)
# end

for f ∈ [:vshl, :vashr, :vlshr]
  @eval begin
    @inline $f(i::MM{W,X,T}, v::SignedHW) where {W,X,T<:SignedHW} = $f(Vec(i), v)
    @inline $f(i::MM{W,X,T}, v::SignedHW) where {W,X,T<:UnsignedHW} = $f(Vec(i), v)
    @inline $f(i::MM{W,X,T}, v::UnsignedHW) where {W,X,T<:SignedHW} = $f(Vec(i), v)
    @inline $f(i::MM{W,X,T}, v::UnsignedHW) where {W,X,T<:UnsignedHW} = $f(Vec(i), v)
    @inline $f(i::MM{W,X,T}, v::IntegerTypesHW) where {W,X,T<:StaticInt} = $f(Vec(i), v)

    @inline $f(v::SignedHW, i::MM{W,X,T}) where {W,X,T<:SignedHW} = $f(v, Vec(i))
    @inline $f(v::UnsignedHW, i::MM{W,X,T}) where {W,X,T<:SignedHW} = $f(v, Vec(i))
    @inline $f(v::SignedHW, i::MM{W,X,T}) where {W,X,T<:UnsignedHW} = $f(v, Vec(i))
    @inline $f(v::UnsignedHW, i::MM{W,X,T}) where {W,X,T<:UnsignedHW} = $f(v, Vec(i))
    @inline $f(v::IntegerTypesHW, i::MM{W,X,T}) where {W,X,T<:StaticInt} = $f(v, Vec(i))

    @inline $f(i::MM{W,X1,T1}, j::MM{W,X2,T2}) where {W,X1,X2,T1<:SignedHW,T2<:SignedHW} =
      $f(Vec(i), Vec(j))
    @inline $f(i::MM{W,X1,T1}, j::MM{W,X2,T2}) where {W,X1,X2,T1<:UnsignedHW,T2<:SignedHW} =
      $f(Vec(i), Vec(j))
    @inline $f(i::MM{W,X1,T1}, j::MM{W,X2,T2}) where {W,X1,X2,T1<:SignedHW,T2<:UnsignedHW} =
      $f(Vec(i), Vec(j))
    @inline $f(
      i::MM{W,X1,T1},
      j::MM{W,X2,T2},
    ) where {W,X1,X2,T1<:UnsignedHW,T2<:UnsignedHW} = $f(Vec(i), Vec(j))
    @inline $f(
      i::MM{W,X1,T1},
      j::MM{W,X2,T2},
    ) where {W,X1,X2,T1<:StaticInt,T2<:IntegerTypes} = $f(Vec(i), Vec(j))
    @inline $f(
      i::MM{W,X1,T1},
      j::MM{W,X2,T2},
    ) where {W,X1,X2,T1<:IntegerTypes,T2<:StaticInt} = $f(Vec(i), Vec(j))
    @inline $f(i::MM{W,X1,T1}, j::MM{W,X2,T2}) where {W,X1,X2,T1<:StaticInt,T2<:StaticInt} =
      $f(Vec(i), Vec(j))

    @inline $f(
      i::MM{W,X,T1},
      v::AbstractSIMDVector{W,T2},
    ) where {W,X,T1<:SignedHW,T2<:SignedHW} = $f(Vec(i), v)
    @inline $f(
      i::MM{W,X,T1},
      v::AbstractSIMDVector{W,T2},
    ) where {W,X,T1<:UnsignedHW,T2<:SignedHW} = $f(Vec(i), v)
    @inline $f(
      i::MM{W,X,T1},
      v::AbstractSIMDVector{W,T2},
    ) where {W,X,T1<:SignedHW,T2<:UnsignedHW} = $f(Vec(i), v)
    @inline $f(
      i::MM{W,X,T1},
      v::AbstractSIMDVector{W,T2},
    ) where {W,X,T1<:UnsignedHW,T2<:UnsignedHW} = $f(Vec(i), v)
    @inline $f(
      i::MM{W,X,T1},
      v::AbstractSIMDVector{W,T2},
    ) where {W,X,T1<:StaticInt,T2<:IntegerTypes} = $f(Vec(i), v)

    @inline $f(
      v::AbstractSIMDVector{W,T1},
      i::MM{W,X,T2},
    ) where {W,X,T1<:SignedHW,T2<:SignedHW} = $f(v, Vec(i))
    @inline $f(
      v::AbstractSIMDVector{W,T1},
      i::MM{W,X,T2},
    ) where {W,X,T1<:UnsignedHW,T2<:SignedHW} = $f(v, Vec(i))
    @inline $f(
      v::AbstractSIMDVector{W,T1},
      i::MM{W,X,T2},
    ) where {W,X,T1<:SignedHW,T2<:UnsignedHW} = $f(v, Vec(i))
    @inline $f(
      v::AbstractSIMDVector{W,T1},
      i::MM{W,X,T2},
    ) where {W,X,T1<:UnsignedHW,T2<:UnsignedHW} = $f(v, Vec(i))
    @inline $f(
      v::AbstractSIMDVector{W,T1},
      i::MM{W,X,T2},
    ) where {W,X,T1<:IntegerTypes,T2<:StaticInt} = $f(v, Vec(i))
  end
end
# for f ∈ [:vand, :vor, :vxor, :vlt, :vle, :vgt, :vge, :veq, :vne, :vmin, :vmax, :vcopysign]
for f ∈ [:vand, :vor, :vxor, :veq, :vne, :vmin, :vmin_fast, :vmax, :vmax_fast, :vcopysign]
  @eval begin
    @inline $f(i::MM{W,X,T}, v::IntegerTypes) where {W,X,T<:IntegerTypes} = $f(Vec(i), v)
    @inline $f(v::IntegerTypes, i::MM{W,X,T}) where {W,X,T<:IntegerTypes} = $f(v, Vec(i))
    @inline $f(
      i::MM{W,X,T1},
      v::AbstractSIMDVector{W,T2},
    ) where {W,X,T1<:IntegerTypes,T2<:IntegerTypes} = $f(Vec(i), v)
    @inline $f(
      v::AbstractSIMDVector{W,T1},
      i::MM{W,X,T2},
    ) where {W,X,T1<:IntegerTypes,T2<:IntegerTypes} = $f(v, Vec(i))
    @inline $f(
      i::MM{W,X1,T1},
      j::MM{W,X2,T2},
    ) where {W,X1,X2,T1<:IntegerTypes,T2<:IntegerTypes} = $f(Vec(i), Vec(j))
  end
end
for f ∈ [:vdiv, :vrem]
  @eval begin
    @inline $f(
      i::MM{W,X,T1},
      v::AbstractSIMDVector{W,T2},
    ) where {W,X,T1<:SignedHW,T2<:IntegerTypes} = $f(Vec(i), v)
    @inline $f(
      v::AbstractSIMDVector{W,T1},
      i::MM{W,X,T2},
    ) where {W,X,T1<:SignedHW,T2<:IntegerTypes} = $f(v, Vec(i))
    @inline $f(
      i::MM{W,X1,T1},
      j::MM{W,X2,T2},
    ) where {W,X1,X2,T1<:SignedHW,T2<:IntegerTypes} = $f(Vec(i), Vec(j))
    @inline $f(
      i::MM{W,X,T1},
      v::AbstractSIMDVector{W,T2},
    ) where {W,X,T1<:UnsignedHW,T2<:IntegerTypes} = $f(Vec(i), v)
    @inline $f(
      v::AbstractSIMDVector{W,T1},
      i::MM{W,X,T2},
    ) where {W,X,T1<:UnsignedHW,T2<:IntegerTypes} = $f(v, Vec(i))
    @inline $f(
      i::MM{W,X1,T1},
      j::MM{W,X2,T2},
    ) where {W,X1,X2,T1<:UnsignedHW,T2<:IntegerTypes} = $f(Vec(i), Vec(j))
  end
end
for f ∈
    [:vlt, :vle, :vgt, :vge, :veq, :vne, :vmin, :vmax, :vmin_fast, :vmax_fast, :vcopysign]
  @eval begin
    # left floating
    @inline $f(i::MM{W,X,T}, v::IntegerTypes) where {W,X,T<:FloatingTypes} = $f(Vec(i), v)
    @inline $f(
      i::MM{W,X,T1},
      v::AbstractSIMDVector{W,T2},
    ) where {W,X,T1<:FloatingTypes,T2<:IntegerTypes} = $f(Vec(i), v)
    @inline $f(
      v::AbstractSIMDVector{W,T1},
      i::MM{W,X,T2},
    ) where {W,X,T1<:FloatingTypes,T2<:IntegerTypes} = $f(v, Vec(i))
    @inline $f(
      i::MM{W,X1,T1},
      j::MM{W,X2,T2},
    ) where {W,X1,X2,T1<:FloatingTypes,T2<:IntegerTypes} = $f(Vec(i), Vec(j))
    # right floating
    @inline $f(i::MM{W,X,T}, v::FloatingTypes) where {W,X,T<:IntegerTypes} = $f(Vec(i), v)
    @inline $f(v::IntegerTypes, i::MM{W,X,T}) where {W,X,T<:FloatingTypes} = $f(v, Vec(i))
    @inline $f(
      i::MM{W,X,T1},
      v::AbstractSIMDVector{W,T2},
    ) where {W,X,T1<:IntegerTypes,T2<:FloatingTypes} = $f(Vec(i), v)
    @inline $f(
      v::AbstractSIMDVector{W,T1},
      i::MM{W,X,T2},
    ) where {W,X,T1<:IntegerTypes,T2<:FloatingTypes} = $f(v, Vec(i))
    @inline $f(
      i::MM{W,X1,T1},
      j::MM{W,X2,T2},
    ) where {W,X1,X2,T1<:IntegerTypes,T2<:FloatingTypes} = $f(Vec(i), Vec(j))
    # both floating
    @inline $f(i::MM{W,X,T}, v::FloatingTypes) where {W,X,T<:FloatingTypes} = $f(Vec(i), v)
    @inline $f(
      i::MM{W,X,T1},
      v::AbstractSIMDVector{W,T2},
    ) where {W,X,T1<:FloatingTypes,T2<:FloatingTypes} = $f(Vec(i), v)
    @inline $f(
      v::AbstractSIMDVector{W,T1},
      i::MM{W,X,T2},
    ) where {W,X,T1<:FloatingTypes,T2<:FloatingTypes} = $f(v, Vec(i))
    @inline $f(
      i::MM{W,X1,T1},
      j::MM{W,X2,T2},
    ) where {W,X1,X2,T1<:FloatingTypes,T2<:FloatingTypes} = $f(Vec(i), Vec(j))
  end
  if f === :copysign
    @eval begin
      @inline $f(v::Float32, i::MM{W,X,T}) where {W,X,T<:IntegerTypes} = $f(v, Vec(i))
      @inline $f(v::Float32, i::MM{W,X,T}) where {W,X,T<:FloatingTypes} = $f(v, Vec(i))
      @inline $f(v::Float64, i::MM{W,X,T}) where {W,X,T<:IntegerTypes} = $f(v, Vec(i))
      @inline $f(v::Float64, i::MM{W,X,T}) where {W,X,T<:FloatingTypes} = $f(v, Vec(i))
    end
  else
    @eval begin
      @inline $f(v::FloatingTypes, i::MM{W,X,T}) where {W,X,T<:IntegerTypes} = $f(v, Vec(i))
      @inline $f(v::FloatingTypes, i::MM{W,X,T}) where {W,X,T<:FloatingTypes} =
        $f(v, Vec(i))
    end
  end
end

@inline vadd_fast(i::MM{W,Zero}, j::MM{W,Zero}) where {W} =
  vrange(Val{W}(), Int, Val{0}(), Val{2}())
@inline vadd_nsw(i::MM{W,Zero}, j::MM{W,Zero}) where {W} =
  vrange(Val{W}(), Int, Val{0}(), Val{2}())


text:
__precompile__()

module Laziness

import Base: convert, show

export @delay, force

type Promise
    value
    evaluated::Bool
end

Promise(a::Any) = Promise(a, true)

Promise(f::Function) = Promise(f, false)

"""
    @delay(ex)

Introduce laziness to julia.

# Examples

```jldoctest
julia> function foo(x)
           @delay x+1
       end
foo (generic function with 1 method)

julia> p = foo(2)
Promise
value → #[...]

julia> force(p)
3

julia> p
Evaluated promise
value → 3
```

"""

macro delay(ex)
    quote Promise(()->$(esc(ex))) end
end

"""
    force(p::Promise) -> Any

To 'force' a promise means to run its code.
The results will be saved, which means forcing twice will run only once.

See also: `@delay`
"""

force(p::Promise) = begin
    if p.evaluated
        p.value
    else
        p.value = p.value()
        p.evaluated = true
        p.value
    end
end

convert(::Type{T}, p::Promise) where{T} = convert(T, force(p))

Base.show(io::IO, p::Promise) = begin
    println(io, p.evaluated ? "Evaluated promise" : "Promise")
    print(io, "value → $(p.value)")
end

end #module

max_stars_repo_name: sclereid/collection

max_issues_repo_name: sclereid/collection

max_forks_repo_name: sclereid/collection


text:
#!/usr/bin/env julia

module TestSwitch

using FactCheck
using Switch

# basic switch constructs
function switch_with_fallthrough(x::Int)
    xs = Int[]
    @switch x begin
        @case 1
            push!(xs, 1)

        @case 2
            push!(xs, 2)

        @case 3
            push!(xs, 3)

        @default
            push!(xs, 4)
    end
    return xs
end

# fallthrough with weird default case placement
function switch_with_leading_default(x::Int)
    xs = Int[]
    @switch x begin
        @default
            push!(xs, 4)

        @case 1
            push!(xs, 1)

        @case 2
            push!(xs, 2)

        @case 3
            push!(xs, 3)
    end
    return xs
end

# switching with expressions
function switch_with_expressions(x::Int)
    xs = Int[]
    @switch x + 1 begin
        @case 1 + 1
            push!(xs, 1)

        @case 2 + 1
            push!(xs, 2)

        @case 3 + 1
            push!(xs, 3)

        @default
            push!(xs, 4)
    end
    return xs
end

# switching with breaks
function switch_with_break(x::Int)
    xs = Int[]
    @switch x begin
        @case 1
            push!(xs, 1)
            break

        @case 2
            push!(xs, 2)
            break

        @case 3
            push!(xs, 3)
            break

        @default
            push!(xs, 4)
    end
    return xs
end

facts("switching") do
    @fact switch_with_fallthrough(1) => Int[1, 2, 3, 4]
    @fact switch_with_fallthrough(3) => Int[3, 4]
    @fact switch_with_fallthrough(100) => Int[4]

    @fact switch_with_leading_default(1) => Int[1, 2, 3]
    @fact switch_with_leading_default(3) => Int[3]
    @fact switch_with_leading_default(100) => Int[4, 1, 2, 3]

    @fact switch_with_expressions(1) => Int[1, 2, 3, 4]
    @fact switch_with_expressions(3) => Int[3, 4]
    @fact switch_with_expressions(100) => Int[4]

    @fact switch_with_break(1) => Int[1]
    @fact switch_with_break(3) => Int[3]
    @fact switch_with_break(100) => Int[4]
end

exitstatus()

end


text:
module test_pkgs_lazy

using Base.Test
using Lazy: @>

@test '가'         == @> 0xac00 Char
@test ['가', '가'] == @> 0xac00 Char fill(2)

end # module test_pkgs_lazy


text:
function IntervalsSet(arg0::BSPTree, arg1::jdouble)
    return IntervalsSet((BSPTree, jdouble), arg0, arg1)
end

function IntervalsSet(arg0::Collection, arg1::jdouble)
    return IntervalsSet((Collection, jdouble), arg0, arg1)
end

function IntervalsSet(arg0::jdouble)
    return IntervalsSet((jdouble,), arg0)
end

function IntervalsSet(arg0::jdouble, arg1::jdouble, arg2::jdouble)
    return IntervalsSet((jdouble, jdouble, jdouble), arg0, arg1, arg2)
end

function apply_transform(obj::AbstractRegion, arg0::Transform)
    return jcall(obj, "applyTransform", AbstractRegion, (Transform,), arg0)
end

function as_list(obj::IntervalsSet)
    return jcall(obj, "asList", List, ())
end

function build_new(obj::IntervalsSet, arg0::BSPTree)
    return jcall(obj, "buildNew", IntervalsSet, (BSPTree,), arg0)
end

function check_point(obj::AbstractRegion, arg0::JVector)
    return jcall(obj, "checkPoint", Region_Location, (JVector,), arg0)
end

function check_point(obj::AbstractRegion, arg0::Point)
    return jcall(obj, "checkPoint", Region_Location, (Point,), arg0)
end

function contains(obj::AbstractRegion, arg0::Region)
    return jcall(obj, "contains", jboolean, (Region,), arg0)
end

function copy_self(obj::AbstractRegion)
    return jcall(obj, "copySelf", AbstractRegion, ())
end

function for_each(obj::Iterable, arg0::Consumer)
    return jcall(obj, "forEach", void, (Consumer,), arg0)
end

function get_barycenter(obj::AbstractRegion)
    return jcall(obj, "getBarycenter", Point, ())
end

function get_boundary_size(obj::AbstractRegion)
    return jcall(obj, "getBoundarySize", jdouble, ())
end

function get_inf(obj::IntervalsSet)
    return jcall(obj, "getInf", jdouble, ())
end

function get_size(obj::AbstractRegion)
    return jcall(obj, "getSize", jdouble, ())
end

function get_sup(obj::IntervalsSet)
    return jcall(obj, "getSup", jdouble, ())
end

function get_tolerance(obj::AbstractRegion)
    return jcall(obj, "getTolerance", jdouble, ())
end

function get_tree(obj::AbstractRegion, arg0::jboolean)
    return jcall(obj, "getTree", BSPTree, (jboolean,), arg0)
end

function intersection(obj::AbstractRegion, arg0::SubHyperplane)
    return jcall(obj, "intersection", SubHyperplane, (SubHyperplane,), arg0)
end

function is_empty(obj::AbstractRegion)
    return jcall(obj, "isEmpty", jboolean, ())
end

function is_empty(obj::AbstractRegion, arg0::BSPTree)
    return jcall(obj, "isEmpty", jboolean, (BSPTree,), arg0)
end

function is_full(obj::AbstractRegion)
    return jcall(obj, "isFull", jboolean, ())
end

function is_full(obj::AbstractRegion, arg0::BSPTree)
    return jcall(obj, "isFull", jboolean, (BSPTree,), arg0)
end

function iterator(obj::IntervalsSet)
    return jcall(obj, "iterator", JIterator, ())
end

function project_to_boundary(obj::IntervalsSet, arg0::Point)
    return jcall(obj, "projectToBoundary", BoundaryProjection, (Point,), arg0)
end

function spliterator(obj::Iterable)
    return jcall(obj, "spliterator", Spliterator, ())
end


text:
function dfs(graph, v)
    n = nv(graph)
    visited = zeros(Bool, n)
    list = Int64[]
    cycleclosures = Vector{Int64}[]
    dfs!(graph, v, 0, list, cycleclosures, visited)

    # removes double entries of cycle connections and keeps the first found pair
    cycleclosures = cycleclosures[sortperm(sort.(cycleclosures))[1:2:length(cycleclosures)]]

    return list, cycleclosures
end

function dfs!(graph, v, p, list, cycleclosures, visited)
    if !visited[v]
        visited[v] = true
        for node in all_neighbors(graph, v)
            if node != p && visited[node]
                push!(cycleclosures, [v; node])
            end
            dfs!(graph, node, v, list, cycleclosures, visited)
        end
        push!(list, v)
    end

    return
end


text:
# The functions in this file are meant to be used for testing the vortex
# evolution routines written using Chebyshev series.  They are not designed to
# fit with other rest of code.

struct JoukowskyMap
    c::ComplexF64
    α::Float64
end

function (m::JoukowskyMap)(ζ, derivative=0)
    if derivative == 0
        return joukowsky(ζ, m.c, m.α)
    elseif derivative == 1
        return joukowsky_d(ζ, m.α)
    elseif derivative == 2
        return joukowsky_d2(ζ, m.α)
    else
        throw(ArgumentError("JoukowskyMap is only implemented up to 2nd derivative"))
    end
end

joukowsky(ζ, c, α) = c + 0.5exp(im*α)*(ζ + 1/ζ)
joukowsky_d(ζ, α) = 0.5exp(im*α)*(1 - 1/ζ^2)
joukowsky_d2(ζ, α) = exp(im*α)/ζ^3

W_motion(ζ, ċ, α, Ω) = im*(imag(ċ*exp(-im*α))/ζ^2 + 0.5Ω/ζ^3)

function W_vortex(ζ, ζs, Γs)
    wᵥ = zero(Complex)
    for (ζᵥ, Γᵥ) in zip(ζs, Γs)
        wᵥ += Γᵥ*(1/(ζ - ζᵥ) - 1/(ζ - 1.0/conj(ζᵥ)))
    end
    return wᵥ/(2π*im)
end


text:
incx = 1
incy = 1
incz = 1


text:
# functions for creating a new mesh object using the same underlying
# apf::Mesh as an existing mesh object.

"""
  This function returns a new mesh object that uses the same underlying
  apf::Mesh as an existing mesh object

  **Inputs**

   * mesh: the existing mesh
   * sbp: an SBP operator
   * opts: options dictonary
   * sbpface: an `AbstractFace` that is compatible with `sbp`. If `sbp` is
              the same `sbp` used to create `mesh`, then this argument is not
              required

  **Outputs**

   * newmesh: the new mesh object

  Note that `sbp` can be a different SBP operator than was used to create `mesh`,
  either a different type of a different degree.

"""
function copy_mesh(mesh::PumiMesh, sbp::AbstractOperator, opts,
                   sbpface::AbstractFace=mesh.sbpface)

  error("generic fallback reached")
end

function copy_mesh(mesh::PumiMeshDG2, sbp::AbstractOperator, opts,
                   sbpface::AbstractFace=mesh.sbpface)

  return PumiMeshDG2(mesh, sbp, opts, sbpface)
end

function copy_mesh(mesh::PumiMeshDG3, sbp::AbstractOperator, opts,
                   sbpface::AbstractFace=mesh.sbpface)

  return PumiMeshDG3(mesh, sbp, opts, sbpface)

end


text:
@api_default function raw_query(api::API, company::Company, q; options...)
    res = hget(api, "/company/$(company.Id)/query", params=Dict("query"=>q); options...)
    Tree(JSON.parse(String(res.data); floattype=Decimal))
end

@api_default function query(api::API, company::Company, q; options...)
    r = raw_query(api, company, q; options...).x
    [map(x->from_json(key_to_type(kv[1]), x), kv[2]) for kv in filter(
        (k,_)->k ∉ ("startPosition", "maxResults", "totalCount"), r["QueryResponse"])]
end


text:
module TestSort

using DataFrames, Random, Test, CategoricalArrays

@testset "standard tests" begin
    dv1 = [9, 1, 8, missing, 3, 3, 7, missing]
    dv2 = [9, 1, 8, missing, 3, 3, 7, missing]
    dv3 = Vector{Union{Int, Missing}}(1:8)
    cv1 = CategoricalArray(dv1, ordered=true)

    d = DataFrame(dv1 = dv1, dv2 = dv2, dv3 = dv3, cv1 = cv1)

    @test sortperm(d) == sortperm(dv1)
    @test sortperm(d[:, [:dv3, :dv1]]) == sortperm(dv3)
    @test sort(d, :dv1)[!, :dv3] == sort(d, "dv1")[!, "dv3"] == sortperm(dv1)
    @test sort(d, :dv2)[!, :dv3] == sortperm(dv1)
    @test sort(d, :cv1)[!, :dv3] == sortperm(dv1)
    @test sort(d, [:dv1, :cv1])[!, :dv3] == sortperm(dv1)
    @test sort(d, [:dv1, :dv3])[!, :dv3] == sortperm(dv1)

    df = DataFrame(rank=rand(1:12, 1000),
                   chrom=rand(1:24, 1000),
                   pos=rand(1:100000, 1000))

    @test issorted(sort(df))
    @test issorted(sort(df, rev=true), rev=true)
    @test issorted(sort(df, [:chrom, :pos])[:, [:chrom, :pos]])
    @test issorted(sort(df, ["chrom", "pos"])[:, ["chrom", "pos"]])

    ds = sort(df, [order(:rank, rev=true), :chrom, :pos])
    @test issorted(ds, [order(:rank, rev=true), :chrom, :pos])
    @test issorted(ds, rev=(true, false, false))

    ds = sort(df, [order("rank", rev=true), "chrom", "pos"])
    @test issorted(ds, [order("rank", rev=true), "chrom", "pos"])
    @test issorted(ds, rev=(true, false, false))

    ds2 = sort(df, [:rank, :chrom, :pos], rev=(true, false, false))
    @test issorted(ds2, [order(:rank, rev=true), :chrom, :pos])
    @test issorted(ds2, rev=(true, false, false))

    @test ds2 == ds

    ds2 = sort(df, ["rank", "chrom", "pos"], rev=(true, false, false))
    @test issorted(ds2, [order("rank", rev=true), "chrom", "pos"])
    @test issorted(ds2, rev=(true, false, false))

    @test ds2 == ds

    sort!(df, [:rank, :chrom, :pos], rev=(true, false, false))
    @test issorted(df, [order(:rank, rev=true), :chrom, :pos])
    @test issorted(df, rev=(true, false, false))

    @test df == ds

    sort!(df, ["rank", "chrom", "pos"], rev=(true, false, false))
    @test issorted(df, [order("rank", rev=true), "chrom", "pos"])
    @test issorted(df, rev=(true, false, false))

    @test df == ds

    @test_throws ArgumentError sort(df, (:rank, :chrom, :pos))

    df = DataFrame(x = [3, 1, 2, 1], y = ["b", "c", "a", "b"])
    @test !issorted(df, :x)
    @test issorted(sort(df, :x), :x)

    df = DataFrame(x = [3, 1, 2, 1], y = ["b", "c", "a", "b"])
    @test !issorted(df, "x")
    @test issorted(sort(df, "x"), "x")

    x = DataFrame(a=1:3, b=3:-1:1, c=3:-1:1)
    @test issorted(x)
    @test !issorted(x, [:b, :c])
    @test !issorted(x[:, 2:3], [:b, :c])
    @test issorted(sort(x, [2, 3]), [:b, :c])
    @test issorted(sort(x[:, 2:3]), [:b, :c])

    x = DataFrame(a=1:3, b=3:-1:1, c=3:-1:1)
    @test issorted(x)
    @test !issorted(x, ["b", "c"])
    @test !issorted(x[:, 2:3], ["b", "c"])
    @test issorted(sort(x, [2, 3]), ["b", "c"])
    @test issorted(sort(x[:, 2:3]), ["b", "c"])

    # Check that columns that shares the same underlying array are only permuted once PR#1072
    df = DataFrame(a=[2, 1])
    df.b = df.a
    sort!(df, :a)
    @test df == DataFrame(a=[1, 2], b=[1, 2])

    x = DataFrame(x=[1, 2, 3, 4], y=[1, 3, 2, 4])
    sort!(x, :y)
    @test x.y == [1, 2, 3, 4]
    @test x.x == [1, 3, 2, 4]

    @test_throws ArgumentError sort(x, by=:x)

    Random.seed!(1)
    # here there will be probably no ties
    df_rand1 = DataFrame(rand(100, 4), :auto)
    # but here we know we will have ties
    df_rand2 = copy(df_rand1)
    df_rand2.x1 = shuffle([fill(1, 50); fill(2, 50)])
    df_rand2.x4 = shuffle([fill(1, 50); fill(2, 50)])

    # test sorting by 1 column
    for df_rand in [df_rand1, df_rand2]
        # testing sort
        for n1 in names(df_rand)
            # passing column name
            @test sort(df_rand, n1) == df_rand[sortperm(df_rand[:, n1]), :]
            # passing vector with one column name
            @test sort(df_rand, [n1]) == df_rand[sortperm(df_rand[:, n1]), :]
            # passing vector with two column names
            for n2 in setdiff(names(df_rand), [n1])
                @test sort(df_rand, [n1, n2]) ==
                      df_rand[sortperm(collect(zip(df_rand[:, n1],
                                                   df_rand[:, n2]))), :]
            end
        end
        # testing if sort! is consistent with issorted and sort
        ref_df = df_rand
        for n1 in names(df_rand)
            df_rand = copy(ref_df)
            @test sort!(df_rand, n1) == sort(ref_df, n1)
            @test issorted(df_rand, n1)
            df_rand = copy(ref_df)
            @test sort!(df_rand, [n1]) == sort(ref_df, [n1])
            @test issorted(df_rand, [n1])
            for n2 in setdiff(names(df_rand), [n1])
                df_rand = copy(ref_df)
                @test sort!(df_rand, [n1, n2]) == sort(ref_df, [n1, n2])
                @test issorted(df_rand, n1)
                @test issorted(df_rand, [n1, n2])
            end
        end
    end
end

@testset "non standard selectors" begin
    Random.seed!(1234)
    df = DataFrame(rand(1:2, 1000, 4), :auto)
    for f in [sort, sort!, sortperm, issorted]
        @test f(df) == f(df, :) == f(df, All()) == f(df, Cols(:)) == f(df, r"x") ==
              f(df, Between(1, 4)) == f(df, Not([]))
    end
end

@testset "view kwarg test" begin
    df = DataFrame(rand(3, 4), :auto)
    @test sort(df) isa DataFrame
    @inferred sort(df)
    @test sort(view(df, 1:2, 1:2)) isa DataFrame
    @test sort(df, view=false) isa DataFrame
    @test sort(view(df, 1:2, 1:2), view=false) isa DataFrame
    @test sort(df, view=true) isa SubDataFrame
    @test sort(df, view=true) == sort(df)
    @test sort(view(df, 1:2, 1:2), view=true) isa SubDataFrame
    @test sort(view(df, 1:2, 1:2), view=true) == sort(view(df, 1:2, 1:2))
end

end # module


text:
# testing json export by
# exporting, then importing
@testset "all in one" begin
    # toDict

    ## FottedShape

    # plane
    p1 = SVector(15.6, 0, -13.7)
    n1 = normalize(SVector(34,45,7))
    s_plane = FittedPlane(p1, n1)
    d_plane = RANSAC.toDict(s_plane)

    @test d_plane == Dict("type"=>"plane", "point"=>p1, "normal"=>n1)

    # sphere
    s_sphere = FittedSphere(p1, 13.23444, true)
    d_sphere = RANSAC.toDict(s_sphere)

    @test d_sphere == Dict("type"=>"sphere", "radius"=>13.23444, "center"=>p1, "outwards"=>true)

    # cylinder
    a1 = SVector(-17.1, 8, 2.42)
    s_cylinder = FittedCylinder(a1, p1, 0.13, false)
    d_cylinder = RANSAC.toDict(s_cylinder)

    @test d_cylinder == Dict("type"=>"cylinder", "axis"=>a1, "center"=>p1, "radius"=>0.13, "outwards"=>false)

    # cone
    ap1 = SVector(0.0,0,0)
    ax1 = normalize(SVector(-1.5, 7, 2))
    s_cone = FittedCone(ap1, ax1, 0.785, true)
    d_cone = RANSAC.toDict(s_cone)

    @test d_cone == Dict("type"=>"cone", "apex"=>ap1, "axis"=>ax1, "opang"=>0.785, "outwards"=>true)

    ## ExtractedShape
    sc1 = ExtractedShape(s_plane, [1])
    @test RANSAC.toDict(sc1) == d_plane

    ss1 = ExtractedShape(s_cone, [1,2,3])
    @test RANSAC.toDict(ss1) == d_cone

    ## Array of primitives
    sa = [s_plane, s_sphere, s_cylinder, s_cone]
    sa_dict = [d_plane, d_sphere, d_cylinder, d_cone]

    d_sa = RANSAC.toDict(sa)

    @test d_sa == Dict("primitives"=>sa_dict)

    ## Array of ExtractedShape
    sc_a = [sc1, ss1]
    sc_a_dict = [d_plane, d_cone]

    @test RANSAC.toDict(sc_a) == Dict("primitives"=>sc_a_dict)

    # exportJSON

    ## with indentation
    io1 = IOBuffer()
    exportJSON(io1, s_cone, 2)
    string_io1 = String(take!(io1))
    #@test string_io1 == "{\n  \"opang\": 0.785,\n  \"outwards\": true,\n  \"axis\": [\n    -0.20180183819889375,\n    0.9417419115948374,\n    0.269069117598525\n  ],\n  \"apex\": [\n    0.0,\n    0.0,\n    0.0\n  ],\n  \"type\": \"cone\"\n}\n"
    pd1 = JSON.parse(string_io1)
    @test pd1 == d_cone

    io2 = IOBuffer()
    exportJSON(io2, [s_plane, s_cone], 1)
    string_io2 = String(take!(io2))
    #@test string_io2 == "{\n \"primitives\": [\n  {\n   \"point\": [\n    15.6,\n    0.0,\n    -13.7\n   ],\n   \"normal\": [\n    0.5982430416161189,\n    0.7917922609625102,\n    0.1231676850386127\n   ],\n   \"type\": \"plane\"\n  },\n  {\n   \"opang\": 0.785,\n   \"outwards\": true,\n   \"axis\": [\n    -0.20180183819889375,\n    0.9417419115948374,\n    0.269069117598525\n   ],\n   \"apex\": [\n    0.0,\n    0.0,\n    0.0\n   ],\n   \"type\": \"cone\"\n  }\n ]\n}\n"
    pd2 = JSON.parse(string_io2)
    @test pd2 == Dict("primitives"=>[d_plane, d_cone])

    ## Without indentation
    io3 = IOBuffer()
    exportJSON(io3, s_cylinder)
    string_io3 = String(take!(io3))
    #@test string_io3 == "{\"outwards\":false,\"axis\":[-17.1,8.0,2.42],\"radius\":0.13,\"center\":[15.6,0.0,-13.7],\"type\":\"cylinder\"}"
    pd3 = JSON.parse(string_io3)
    @test pd3 == d_cylinder

    io4 = IOBuffer()
    exportJSON(io4, [s_plane, s_sphere, s_cone, s_plane])
    string_io4 = String(take!(io4))
    #@test string_io4 == "{\"primitives\":[{\"point\":[15.6,0.0,-13.7],\"normal\":[0.5982430416161189,0.7917922609625102,0.1231676850386127],\"type\":\"plane\"},{\"outwards\":true,\"radius\":13.23444,\"center\":[15.6,0.0,-13.7],\"type\":\"sphere\"},{\"opang\":0.785,\"outwards\":true,\"axis\":[-0.20180183819889375,0.9417419115948374,0.269069117598525],\"apex\":[0.0,0.0,0.0],\"type\":\"cone\"},{\"point\":[15.6,0.0,-13.7],\"normal\":[0.5982430416161189,0.7917922609625102,0.1231676850386127],\"type\":\"plane\"}]}"
    pd4 = JSON.parse(string_io4)
    @test pd4 == Dict("primitives"=>[d_plane, d_sphere, d_cone, d_plane])

    ## Write to file w&wo indentation
    f1, fio1 = mktemp(pwd())
    close(fio1)
    open(f1, "w") do io
        exportJSON(io, [s_plane, s_sphere, s_cone, s_plane], 4)
    end
    str_f1 = read(f1, String)
    pd_f1 = JSON.parse(str_f1)
    @test pd_f1 == Dict("primitives"=>[d_plane, d_sphere, d_cone, d_plane])
    Base.Filesystem.rm(f1)

    f2, fio2 = mktemp(pwd())
    close(fio2)
    open(f2, "w") do io
        exportJSON(io, [s_plane, s_sphere, s_cone, s_plane, s_cylinder])
    end
    str_f2 = read(f2, String)
    pd_f2 = JSON.parse(str_f2)
    @test pd_f2 == Dict("primitives"=>[d_plane, d_sphere, d_cone, d_plane, d_cylinder])
    Base.Filesystem.rm(f2)
end


text:
# This file contains functions used in the high level service definitions
_merge(a::AbstractDict, b::AbstractDict) = merge(a, b)
_merge(a::AbstractString, b::AbstractString) = b

function _flatten_query(service::String, query::AbstractDict{String,<:Any})
    return _flatten_query!(Pair{String,String}[], service, query)
end

function _flatten_query!(
    result::Vector{Pair{String,String}},
    service::String,
    query::AbstractDict{String,<:Any},
    prefix::String="",
)
    for (k, v) in query
        if v isa AbstractDict
            _flatten_query!(result, service, v, string(prefix, k, "."))
        elseif v isa AbstractArray
            for (i, j) in enumerate(v)
                suffix = service in ("ec2", "sqs") ? "" : ".member"
                prefix_key = string(prefix, k, suffix, ".", i)

                if j isa AbstractDict
                    _flatten_query!(result, service, j, string(prefix_key, "."))
                else
                    push!(result, Pair(prefix_key, string(j)))
                end
            end
        else
            push!(result, Pair(string(prefix, k), string(v)))
        end
    end

    return result
end

"""
Escape special AWS S3 characters in the path of the uri properly.

AWS S3 allows for various special characters in file names, these characters are not being
properly escaped before we make the requests.

We cannot call `HTTP.escapeuri(request.uri)` because this will escape `/` characters which
are used in the filepathing for sub-directories.
"""
function _clean_s3_uri(uri::AbstractString)
    chars_to_clean = (
        ' ' => "%20",
        '!' => "%21",
        ''' => "%27",
        '(' => "%28",
        ')' => "%29",
        '*' => "%2A",
        '+' => "%2B",
        ',' => "%2C",
        ':' => "%3A",
        '=' => "%3D",
        '@' => "%40",
    )
    parsed_uri = URIs.URI(uri)
    cleaned_path = reduce(replace, chars_to_clean; init=parsed_uri.path)
    return string(URIs.URI(parsed_uri; path=cleaned_path))
end

function _extract_common_kw_args(service, args)
    return (
        service=service.signing_name,
        api_version=service.api_version,
        response_stream=_pop!(args, "response_stream", nothing),
        headers=LittleDict{String,String}(_pop!(args, "headers", [])),
        http_options=_pop!(args, "http_options", LittleDict{Symbol,String}()),
        backend=_pop!(args, "backend", DEFAULT_BACKEND[]),

        # Deprecated keywords
        return_stream=_pop!(args, "return_stream", nothing),
        return_raw=_pop!(args, "return_raw", nothing),
        response_dict_type=_pop!(args, "response_dict_type", nothing),
    )
end

function _delete_legacy_response_kw_args!(args)
    used_kw_args = intersect(
        keys(args), ("return_headers", "return_stream", "return_raw", "response_dict_type")
    )

    if !isempty(used_kw_args)
        Base.depwarn(
            "The parameter(s) \"$(join(used_kw_args, "\", \""))\" are no longer " *
            "supported when the feature `use_response_type` is enabled.",
            :_delete_legacy_response_kw_args!,
        )
    end

    delete!(args, "return_headers")
    delete!(args, "return_stream")
    delete!(args, "return_raw")
    delete!(args, "response_dict_type")
    return args
end

# Use this until the three arg pop! is available for LittleDict
# https://github.com/JuliaCollections/OrderedCollections.jl/pull/59
function _pop!(dict::AbstractDict{String,<:Any}, kw, default)
    if haskey(dict, kw)
        val = dict[kw]
        delete!(dict, kw)

        return val
    else
        return default
    end
end

function _generate_rest_resource(request_uri::String, args::AbstractDict{String,<:Any})
    # There maybe a time where both $k and $k+ are in the request_uri, in which case this needs to be updated
    # From looking around, I have not seen an example yet, however that doesn't mean it doesn't exist

    for (k, v) in args
        if occursin("{$k}", request_uri)
            request_uri = replace(request_uri, "{$k}" => v)
        elseif occursin("{$k+}", request_uri)
            request_uri = replace(request_uri, "{$k+}" => HTTP.escapepath(v))
        end
    end

    return request_uri
end

function _assignment_to_kw!(expr::Expr)
    if expr.head === :(=)
        expr.head = :kw
    else
        throw(ArgumentError("Expected assignment expression, instead found: `$expr`"))
    end

    return expr
end

function _assignment_to_kw!(x)
    return throw(ArgumentError("Expected assignment expression, instead found: `$x`"))
end


text:
include("./src/functions/Functions.jl")
include("./src/evaluate/Evaluate.jl")
include("./src/linesearch/Linesearch.jl")

using MAT

function logregGD(X, y, lambda, epsilon=1e-4, maxiter=1000,maxLSiter=25, verbose=false)
	m = size(X,1) # i have m datapoints
	n = size(X,2) # n = dimension of the problem
	theta = zeros(784,1)
	multiplier = 0.5
	f_likelihood = f_logreg(X,y,theta,lambda)
	g_likelihood = g_logreg(X,y,theta,lambda)
	optimal = false

	likelihood = Inf
	for i in 1:maxiter
        likelihood = f_likelihood(theta)
        gradient = g_likelihood(theta)
		alpha = 1

		# check optimality with inf norm
		if norm(gradient,Inf) < epsilon
			return (theta,"optimal")
		end
		# find step size
		lineSearchIter = 0
        # Linesearch
        for j in 1:maxLSiter
            f_xk_alpha_k = f_likelihood(theta + ( alpha * -gradient))
            f_xk_plus_grad = likelihood + (epsilon*alpha*gradient' * -gradient)[1]
            #@printf("%f \t < \t %f\n",f_xk_alpha_k,f_xk_plus_grad)

            if  ! (f_xk_alpha_k < f_xk_plus_grad)
                alpha = alpha * multiplier
                lineSearchIter += 1
            else
                break
            end
        end

		if verbose
            @printf("%i \t\t %f \t %f \t %f \t\t %i \n",i,likelihood,norm(gradient,Inf),alpha,lineSearchIter)
		end
		# descent
		theta = theta + alpha * -gradient
	end
	return (theta,"maxIter")
end

dataLocation = "data/mnist67.scale.1k.mat"

# 784 dimensions | 1000 data points
file = matopen(dataLocation)
tmp = read(file, "X")

X = full(tmp)
y = read(file,"y")

theta, state = logregGD(X,y,0,1e-4,1000,20,true)
@printf("Optimization finished with state %s.\n",state)
accuracy = evaluate(X,y,theta)
@printf("Accuracy of solution : %f %%.\n",accuracy*100)

theta, state = logregGD(X,y,0.1,1e-4,1000,20,true)
@printf("Optimization finished with state %s.\n",state)
accuracy = evaluate(X,y,theta)
@printf("Accuracy of solution : %f %%.\n",accuracy*100)

theta, state = logregGD(X,y,100,1e-4,1000,20,true)
@printf("Optimization finished with state %s.\n",state)
accuracy = evaluate(X,y,theta)
@printf("Accuracy of solution : %f %%.\n",accuracy*100)

theta, state = logregGD(X,y,1000,1e-4,1000,20,true)
@printf("Optimization finished with state %s.\n",state)
accuracy = evaluate(X,y,theta)
@printf("Accuracy of solution : %f %%.\n",accuracy*100)

max_stars_repo_name: zagazao/num-opt

max_issues_repo_name: zagazao/num-opt

max_forks_repo_name: zagazao/num-opt


text:
#= NOTE createReservedWords2.jl also has this function
F: get rid of the trailing newlines, comments, blank lines in a file;
I: file path, read line-by-line;
O: array of strings, each element maps to a line in the file;
B: load the whole file at the same time, may require large memory for large file;
A: read line-by-line, ingnore blank lines and pure comments, throw away in-line comments;
=#
function getRidOfNewlineHashcommentBlankline(inputFilePath::AbstractString)
  f = open(inputFilePath, "r")
  lines = readlines(f)
  newlines = Array{Any, 1}()  # container for results
  for (id, val) in enumerate(lines)
    val = chomp(val)  # trailing newline
    hashId = findfirst(isequal('#'), val)
    if something(hashId,-1) != -1  # get rid of comments
      val = val[1:hashId-1]
    end
    if length(split(val)) != 0  # non-empty
      # throw away ugly trailing spaces, although can be done when splitting
      while(val[end] == ' ')
        val = val[1:end-1]
      end
      push!(newlines, (val, id))
      println(val)
    end
  end
  # println(">"^20)
  # println(newlines)
  return newlines
end

#= NOTE tokenization -> assume can be done just by spaces
FIOBA: using "split" directly
=#
function sentenceTokenization1(inputSenArray::Array)
  tokenizedInputSentencesArray = Array{Array{String,1}, 1}()
  for sen in inputSentencesArray
    tokens = split(sen)
    push!(tokenizedInputSentencesArray, tokens)
  end
  println(typeof(tokenizedInputSentencesArray), size(tokenizedInputSentencesArray))
  return tokenizedInputSentencesArray
end
#= NOTE
F: tokenizing a sentence by processing it char-by-char.
I: output of getRidOfNewlineHashcommentBlankline.
O: Array of Array of tokens.
B: .
A: tokenizing a sentence by processing it char-by-char, punctuation and left/right
   parenthesis are treated as independent tokens, other tokens are seperated by
   space or punctuation set.
=#
function sentenceTokenization2(inputSenArray::Array)
  puncSet = Set(['|', '&', ',', '(', ')'])
  tokenizedInputSentencesArray = Array{Any, 1}()
  for (sen, id) in inputSenArray
    sen = sen*" "  # add a space to end the sentence
    tokens = Array{String,1}()
    currentToken = ""
    for char in sen
      if char == ' '  # encounter " ", push! if have someting
        if !isempty(currentToken)
          push!(tokens, currentToken)
          currentToken = ""
        end
      elseif in(char, puncSet)
        if !isempty(currentToken)
          push!(tokens, currentToken)
          currentToken = ""
        end
        push!(tokens, string(char))
      else
        currentToken *= string(char)
        # println(currentToken)
      end
    end
    push!(tokenizedInputSentencesArray, (tokens, id))
    println(tokens)
  end
  println(typeof(tokenizedInputSentencesArray), size(tokenizedInputSentencesArray))
  return tokenizedInputSentencesArray
end

#= Token classifier 1:
F: get rid of negligible reserved words, tagging or normalizing
I: array{String, 1}, usually from getRidOfNewlineHashcommentBlankline
O: Array of Arrays of Tuples --> Tuple(token, tag)
B: tag all unknown tokens as BioSym, @FIXME how to format for grammar searching?
A: use Dict[tag] = Set{words}, e.g. createReservedWords1.jl
=#
function tokenClassification1(sentencesArray::Array, reservedWords::Dict)
  taggedSentencesArray = Array{Array, 1}()  # container for tagged array
  TAGs = ["induce", "repress", "BioSym"]  # can be incorporated into createReservedWords.jl
  for tokens in sentencesArray  # tag each sentence
    taggedTokens = Array{Tuple,1}()  # Tuple(token, tag)
    for token in tokens  # tag each token in a sentence
      if !(token in reservedWords["negligible"])  # significant words
        breakEndsForLoop = false
        for tag in TAGs  # find a tag
          if token in reservedWords[tag]
            push!(taggedTokens, (token, tag))
            breakEndsForLoop = true
            break
          end
        end
        if !breakEndsForLoop  # token of unknown tagging type FIXME: add a boolean value to indicate status
          # push!(taggedTokens, (token, "UNK"))
          push!(taggedTokens, (token, "BioSym"))  # TODO tag general biological symbols
        end
      end
    end
    push!(taggedSentencesArray, taggedTokens)
  end
  return taggedSentencesArray
end
#= Token classifier 2:
F: get rid of negligible reserved words, tagging or normalizing
I: array{String, 1}, usually from getRidOfNewlineHashcommentBlankline
O: Array of Arrays of Tuples --> Tuple(token, tag)
B: tag all unknown tokens as BioSym, @FIXME how to format for grammar searching?
A: use Dict[word] = tag, e.g. createReservedWords2.jl
=#
function tokenClassification2(sentencesArray::Array, reservedWordsPath::AbstractString)
  include(reservedWordsPath)
  taggedSentencesArray = Array{Any, 1}()  # container for tagged array
  for (tokens, id) in sentencesArray  # tag each sentence
    taggedTokens = Array{Tuple,1}()  # Tuple(token, tag)
    for token in tokens  # tag each token in a sentence
      if haskey(reservedWords, token)  # a reserved word
        tag = reservedWords[token]
        if tag != "negligible"  # important reserved word --> keep it
          push!(taggedTokens, (token, tag))
        end
      else  # non-reserved word --> BioSym
        push!(taggedTokens, (token, "BioSym"))
      end
    end
    push!(taggedSentencesArray, (taggedTokens, id))
  end
  return taggedSentencesArray
end

# # NOTE test
# println("\n---------loading sentences------------")
# inputSentencesArray = getRidOfNewlineHashcommentBlankline("../test/fbacase.txt")
# tokenizedInputSentencesArray = sentenceTokenization2(inputSentencesArray)
#
# println("\n------------normalizing and tagging------------")
# # include("reservedWords.jl")
# taggedSentencesArray = tokenClassification2(tokenizedInputSentencesArray,
#                                             "reservedWords.jl")
# print(taggedSentencesArray)
# println("-----------")
# # reshape for output observation
# printTagSen = [["$(y[1])/$(y[2])" for y in x[1]] for x in taggedSentencesArray]
# println(typeof(printTagSen))
# for ts in printTagSen
#   println(join(ts, "  "))
# end


text:
#greedy algorithm

function greedy{T<:Real}(distmat::Matrix{T})
  dmsz = size(distmat)
  numCities = dmsz[1]
  ordDist = sortperm(vec(distmat))
  G = Graph(numCities)
  k = 1
  #main loop
  while ne(G) < numCities-1
    dnode,anode = ind2sub(dmsz,ordDist[k])
    if dnode != anode
      d°,a° = degree(G,dnode),degree(G,anode)
      if d° < 2 && a° < 2
        if d° == 0 && a° == 0
          add_edge!(G,dnode,anode)
        elseif independent(G,dnode,anode)
          add_edge!(G,dnode,anode)
        end
      end
    end
    k+=1
  end
  #final edge addition
  add_edge!(G,findin(degree(G),1)...)
  return saw(G,1,numCities)
end


text:
function test_cross_validation_accuracy()
    println("-- Testing cross validation with accuracy metrics")

    m = [NaN 3 NaN 1 2 1 NaN 4
         1 2 NaN NaN 3 2 NaN 3
         NaN NaN NaN NaN NaN NaN NaN NaN
         NaN 2 3 3 NaN 5 NaN 1]
    data = DataAccessor(m)

    fold = 5

    # MF(data, 2)
    @test cross_validation(fold, MAE, MF, data, 2) <= 2.5

    # UserMean(data)
    @test cross_validation(fold, MAE, UserMean, data) <= 2.5

    # ItemMean(data)
    @test cross_validation(fold, MAE, ItemMean, data) <= 2.5

    # UserKNN(data, 2, true)
    @test cross_validation(fold, MAE, UserKNN, data, 2, true) <= 2.5
end

function test_cross_validation_ranking()
    println("-- Testing cross validation with ranking metrics")

    m = [NaN 3 NaN 1 2 1 NaN 4
         1 2 NaN NaN 3 2 NaN 3
         NaN 2 3 3 NaN 5 NaN 1]
    data = DataAccessor(m)

    # 5-fold, top-4 recommendation
    fold = 5
    k = 4

    # MF(data, 2)
    @test cross_validation(fold, Recall, k, MF, data, 2) <= 0.5

    # MostPopular(data)
    @test cross_validation(fold, Recall, k, MostPopular, data) <= 0.5

    # UserKNN(data, 2, true)
    @test cross_validation(fold, Recall, k, UserKNN, data, 2, true) <= 0.5
end

test_cross_validation_accuracy()
test_cross_validation_ranking()


text:
module MaxHelpingHandRotatingBumper

using ..Ahorn, Maple

@pardef RotatingBumper(x1::Integer, y1::Integer, x2::Integer=x1 + 16, y2::Integer=y1, speed::Number=360.0, attachToCenter::Bool=false, notCoreMode::Bool=false, wobble::Bool=false) =
    Entity("MaxHelpingHand/RotatingBumper", x=x1, y=y1, nodes=Tuple{Int, Int}[(x2, y2)], speed=speed, attachToCenter=attachToCenter, notCoreMode=notCoreMode, wobble=wobble)

const placements = Ahorn.PlacementDict(
    "Bumper (Rotating) (max480's Helping Hand)" => Ahorn.EntityPlacement(
        RotatingBumper,
        "point",
        Dict{String, Any}(),
        function(entity::RotatingBumper)
            x, y = Int(entity.data["x"]), Int(entity.data["y"])
            entity.data["nodes"] = Tuple{Int, Int}[(x + 32, y)]
        end
    )
)

Ahorn.nodeLimits(entity::RotatingBumper) = 1, 1

sprite = "objects/Bumper/Idle22.png"

function Ahorn.selection(entity::RotatingBumper)
    nx, ny = Int.(entity.data["nodes"][1])
    x, y = Ahorn.position(entity)

    return [Ahorn.getSpriteRectangle(sprite, x, y), Ahorn.getSpriteRectangle(sprite, nx, ny)]
end

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::RotatingBumper, room::Maple.Room)
    x, y = Ahorn.position(entity)

    Ahorn.drawSprite(ctx, sprite, x, y)
end

function Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::RotatingBumper, room::Maple.Room)
    speed = get(entity.data, "speed", 360.0)
    dir = speed > 0 ? 1 : -1

    centerX, centerY = Int.(entity.data["nodes"][1])
    x, y = Ahorn.position(entity)

    radius = sqrt((centerX - x)^2 + (centerY - y)^2)

    Ahorn.drawCircle(ctx, centerX, centerY, radius, Ahorn.colors.selection_selected_fc)
    Ahorn.drawArrow(ctx, centerX + radius, centerY, centerX + radius, centerY + 0.001 * dir, Ahorn.colors.selection_selected_fc, headLength=6)
    Ahorn.drawArrow(ctx, centerX - radius, centerY, centerX - radius, centerY + 0.001 * -dir, Ahorn.colors.selection_selected_fc, headLength=6)
end

end

max_stars_repo_name: Kitty-Cats/MaxHelpingHand

max_issues_repo_name: Kitty-Cats/MaxHelpingHand

max_forks_repo_name: Kitty-Cats/MaxHelpingHand


text:
export expMod, boundMod, fMod,identityMod

"""
	sigma,dsigma = expMod(model)

	maps model parameter to conductivity via

	sigma(m) = exp(m)

"""
function expMod(m)
	return exp.(m), sdiag(exp.(m))
end

"""
	sigma,dsigma = fMod(model;f::Function=identity,df::Function=m->sparse(1.0I,length(m),length(m)))

	maps model parameter to conductivity via

	sigma(m) = f(m) and dsigma(m) = sdiag(df(m))

"""
function fMod(m;f::Function=identity,df::Function=m->ones(eltype(m),length(m)))
	return f(m),sdiag(df(m))
end

function identityMod(m)
	return m, UniformScaling(one(eltype(m)))
end

export boundMod
"""
	sigma,dsigma = boundMod(m;boundLow=0.0,boundHigh=1.0)

	maps model parameter to conductivity via

		sigma = 0.5*(boundHigh-boundLow) * (tanh(m)+1.0 ) + boundLow

"""
function boundMod(m; boundLow = 0.0, boundHigh = 1.0)
	u        = tanh.(m)
	d        = 0.5 * (boundHigh - boundLow)
	sigma    = d * (u .+ 1.0) .+ boundLow
	dsigmadm = sdiag(d * (1.0 .- u .* u))
	return sigma, dsigmadm
end


text:
using ChaosTools
using Test
using StatsBase

test_value = (val, vmin, vmax) -> @test vmin <= val <= vmax

println("\nTesting generalized entropy (genentropy) & linear scaling...")
@testset "Generalized Dimensions" begin
  @testset "Henon Map" begin
    ds = Systems.henon()
    ts = trajectory(ds, 200000)
    mat = Matrix(ts)
    # Test call with dataset
    genentropy(1, 0.001, ts)
    es = 10 .^ range(-0, stop = -3, length = 7)
    dd = zero(es)
    for q in [0,2,1, 2.56]
      for (i, ee) in enumerate(es)
        dd[i] = genentropy(q, ee, mat)
      end
      linr, dim = linear_region(-log.(es), dd)
      test_value(dim, 1.1, 1.3)
    end
  end
  @testset "Lorenz System" begin
    ds = Systems.lorenz()
    ts = trajectory(ds, 5000)
    es = 10 .^ range(1, stop = -3, length = 11)
    dd = zero(es)
    for q in [0,1,2]
      for (i, ee) in enumerate(es)
        dd[i] = genentropy(q, ee, ts)
      end
      linr, dim = linear_region(-log.(es), dd)
      if q == 0
        test_value(dim, 1.7, 1.9)
      else
        test_value(dim, 1.85, 2.0)
      end
    end
  end
end

println("\nTesting dimension calls (all names)...")
@testset "Dimension calls" begin
  ds = Systems.henon()
  ts = trajectory(ds, 20000)
  # Test call with dataset
  test_value(generalized_dim(1.32, ts), 1.1, 1.3)
  test_value(capacity_dim(ts), 1.1, 1.3)
  test_value(information_dim(ts), 1.1, 1.3)
end

println("\nTesting permutation entropy...")
@testset "Permutation Entropy" begin
    @testset "Trivial entropies" begin
        @test permentropy(zeros(10), 3) == 0
        @test permentropy(ones(10), 3) == 0
        @test permentropy(collect(1:10), 3) == 0
        @test permentropy(collect(10:-1:1), 3) == 0
    end
    @testset "Examples" begin
        # Examples from Bandt & Pompe (2002):
        for (xs, order, desired) in [
                ([4, 7, 9, 10, 6, 11, 3], 2,
                 -(4/6)log2(4/6) - (2/6)log2(2/6)),
                ([4, 7, 9, 10, 6, 11, 3], 3,
                 -2(2/5)log2(2/5) - (1/5)log2(1/5)),
            ]
            @test permentropy(xs, order; base=2) ≈ desired
        end
    end
    @testset "User Interface" begin
        order = Int(typemax(UInt8)) + 1
        @test_throws ArgumentError permentropy([], order)
    end
end


text:
using JuliaDB
using IndexedTables
using Base.Test

@testset "join" begin
    t1 = NDSparse(Columns([1,1,2,2], [1,2,1,2]), [1,2,3,4])
    t2 = NDSparse(Columns([0,2,2,3], [1,1,2,2]), [1,2,3,4])

    j1 = innerjoin(t1,t2)
    j2 = innerjoin(+, t1,t2)

    lj1 = leftjoin(t1,t2)
    lj2 = leftjoin(+, t1,t2)

    mj1 = merge(t1,t2)
    mj2 = merge(t2,t1)

    for n=1:5
        d1 = distribute(t1, n)
        for n2 = 1:5
            d2 = distribute(t2, n2)
            @test collect(innerjoin(d1, d2)) == j1
            @test collect(innerjoin(+, d1, d2)) == j2

            @test collect(leftjoin(d1, d2)) == lj1
            @test collect(leftjoin(+, d1, d2)) == lj2

            @test collect(merge(d1, d2)) == mj1
            @test collect(merge(d2, d1)) == mj2
        end
    end

    t1 = NDSparse([:msft,:ibm,:ge], [1,3,4], [100,200,150])
    t2 = NDSparse([:ibm,:msft,:aapl,:ibm], [0,0,0,2], [100,99,101,98])

    aj = asofjoin(t1,t2)

    for n=1:5
        d1 = distribute(t1, n)
        for n2 = 1:5
            d2 = distribute(t2, n2)
            @test collect(asofjoin(d1, d2)) == aj
        end
    end
end

@testset "broadcast" begin
    t1 = NDSparse(Columns([1,2,3,4],[1,1,2,2]), [5,6,7,8])
    t2 = NDSparse(Columns([0,1,2,4]),[0,10,100,10000])

    for n=1:4
        for m=1:4
            d1 = distribute(t1, n)
            d2 = distribute(t2, m)

            @test collect(d1 .+ d2) == t1 .+ t2
            @test collect(broadcast(+, d1, d2, dimmap=(0,1))) == broadcast(+, t1, t2, dimmap=(0,1))
        end
    end
end


text:
@testset "text_file" begin

# test of text_file
sc = SparkContext(master="local")

txt = text_file(sc, @__FILE__)
nums  = map(txt, it -> length(it))

@test reduce(nums, +) == 182

close(sc)

end


text:
using LightXML

macro logmsg(s)
end

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

name2sym(name) = join(map(uppercasefirst, split(name, '-')))

const GEN_TOP = """
# This file is automatically generated.
# Do not edit this file by hand.
# Make changes to gen.jl or the source specification instead.
"""

const GEN_BOTTOM = "# end generated code"

const CLS_TOP = """# Classes
const CLASS_MAP = Dict{TAMQPClassId,ClassSpec}("""

const CLS_BOTTOM = """) # CLASS_MAP")

function make_classmethod_map()
    cmmap = Dict{Tuple{Symbol,Symbol},MethodSpec}()
    for v in values(CLASS_MAP)
        for m in values(v.method_map)
            cmmap[(v.name,m.name)] = m
        end
    end
    cmmap
end

const CLASSNAME_MAP = Dict{Symbol,ClassSpec}(v.name => v for v in values(CLASS_MAP))
const CLASSMETHODNAME_MAP = make_classmethod_map()
# end Classes
"""

const DOMAIN_TOP = "# Domains"
const DOMAIN_BOTTOM = "# end Domains\n"

const domainmap = Dict{String,Type}(
    "bit" => TAMQPBit,
    "octet" => TAMQPOctet,
    "short" => TAMQPShortInt,
    "long" => TAMQPLongInt,
    "longlong" => TAMQPLongLongInt,
    "shortstr" => TAMQPShortStr,
    "longstr" => TAMQPLongStr,
    "timestamp" => TAMQPTimeStamp,
    "table" => TAMQPFieldTable,
    "class-id" => TAMQPClassId,
    "method-id" => TAMQPMethodId
)

const precreated_consts = ["FrameEnd"]
const clsindent = " "^4
const methindent = " "^8
const argsindent = " "^12

function gen_spec(specfile)
    xdoc = parse_file(specfile)
    amqp = root(xdoc)

    println("# Source: ", specfile)
    println(GEN_TOP)

    println("const AMQP_VERSION = v", '"', attribute(amqp, "major"), '.', attribute(amqp, "minor"), '.', attribute(amqp, "revision"), '"')
    println("const AMQP_DEFAULT_PORT = ", attribute(amqp, "port"))
    println("")

    println("# Constants")
    for constant in get_elements_by_tagname(amqp, "constant")
        has_attribute(constant, "class") && continue
        constantname = name2sym(attribute(constant, "name"))
        (constantname in precreated_consts) && continue
        constantvalue = attribute(constant, "value")
        println("const ", constantname, " = ", constantvalue)
    end
    println("")

    println("# Error Codes")
    for constant in get_elements_by_tagname(amqp, "constant")
        !has_attribute(constant, "class") && continue
        cls = attribute(constant, "class")
        name = attribute(constant, "name")
        constantname = name2sym("$cls-$name")
        constantvalue = attribute(constant, "value")
        println("const ", constantname, " = ", constantvalue)
    end
    println("")

    # domains
    println(DOMAIN_TOP)
    for domain in get_elements_by_tagname(amqp, "domain")
        name = attribute(domain, "name")
        (name in keys(domainmap)) && continue
        name = "TAMQP" * name2sym(name)
        typ = domainmap[attribute(domain, "type")].name
        println("const $name = $typ")
    end
    println(DOMAIN_BOTTOM)

    # classes and methods
    println(CLS_TOP)
    clssep = ""
    for cls in get_elements_by_tagname(amqp, "class")
        clsname = Symbol(name2sym(attribute(cls, "name")))
        clsidx = parse(Int, attribute(cls, "index"))

        println(clsindent, clssep, "$clsidx => ClassSpec($clsidx, :$clsname, Dict{TAMQPMethodId, MethodSpec}(")
        isempty(clssep) && (clssep = ", ")

        methsep = ""
        for meth in get_elements_by_tagname(cls, "method")
            methname = Symbol(name2sym(attribute(meth, "name")))
            methidx = parse(Int, attribute(meth, "index"))
            methargs = Pair{Symbol,Type}[]
            methrespelem = find_element(meth, "response")
            methresp = (methrespelem === nothing) ? :Nothing : Symbol(name2sym(attribute(methrespelem, "name")))
            println(methindent, methsep, "$methidx => MethodSpec($methidx, :$methname, :$methresp, Pair{Symbol,DataType}[")
            isempty(methsep) && (methsep = ", ")

            argssep = ""
            for arg in get_elements_by_tagname(meth, "field")
                fieldname = Symbol(name2sym(attribute(arg, "name")))
                fielddomain = attribute(arg, "domain")
                if fielddomain === nothing
                    fielddomain = attribute(arg, "type")
                end
                fieldtype = (fielddomain in keys(domainmap)) ? domainmap[fielddomain].name : ("TAMQP" * name2sym(fielddomain))
                println(argsindent, argssep, ":$fieldname => $fieldtype")
                isempty(argssep) && (argssep = ", ")
            end

            println(methindent, "]) # method $methname")
        end

        println(clsindent, ")) # class $clsname")
    end
    println(CLS_BOTTOM)

    println(GEN_BOTTOM)
end

gen_spec("amqp0-9-1.extended.xml")


text:
struct TrajectoryOptimizationData{T}
    states::Vector{Vector{T}}
    actions::Vector{Vector{T}}
    parameters::Vector{Vector{T}}
    objective::Objective{T}
    dynamics::Vector{Dynamics{T}}
    constraints::Constraints{T}
    bounds::Bounds{T}
    duals_dynamics::Vector{Vector{T}}
    duals_constraints::Vector{Vector{T}}
    state_dimensions::Vector{Int}
    action_dimensions::Vector{Int}
    parameter_dimensions::Vector{Int}
end

function TrajectoryOptimizationData(objective::Objective{T}, dynamics::Vector{Dynamics{T}}, constraints::Constraints{T}, bounds::Bounds{T};
    parameters=[zeros(num_parameter) for num_parameter in dimensions(dynamics)[3]]) where T

    state_dimensions, action_dimensions, parameter_dimensions = dimensions(dynamics)

    states = [zeros(num_state) for num_state in state_dimensions]
    actions = [zeros(nu) for nu in action_dimensions]

    duals_dynamics = [zeros(num_state) for num_state in state_dimensions[2:end]]
    duals_constraints = [zeros(con.num_constraint) for con in constraints]

    TrajectoryOptimizationData(
        states,
        actions,
        parameters,
        objective,
        dynamics,
        constraints,
        bounds,
        duals_dynamics,
        duals_constraints,
        state_dimensions,
        action_dimensions,
        parameter_dimensions)
end

TrajectoryOptimizationData(objective::Objective, dynamics::Vector{Dynamics}) = TrajectoryOptimizationData(objective, dynamics, [Constraint() for t = 1:length(dynamics)], [Bound() for t = 1:length(dynamics)])

struct TrajectoryOptimizationIndices
    objective_hessians::Vector{Vector{Int}}
    dynamics_constraints::Vector{Vector{Int}}
    dynamics_jacobians::Vector{Vector{Int}}
    dynamics_hessians::Vector{Vector{Int}}
    stage_constraints::Vector{Vector{Int}}
    stage_jacobians::Vector{Vector{Int}}
    stage_hessians::Vector{Vector{Int}}
    general_constraint::Vector{Int}
    general_jacobian::Vector{Int}
    general_hessian::Vector{Int}
    states::Vector{Vector{Int}}
    actions::Vector{Vector{Int}}
    state_action::Vector{Vector{Int}}
    state_action_next_state::Vector{Vector{Int}}
end

function indices(objective::Objective{T}, dynamics::Vector{Dynamics{T}}, constraints::Constraints{T}, general::GeneralConstraint{T},
    key::Vector{Tuple{Int,Int}}, num_state::Vector{Int}, num_action::Vector{Int}, num_trajectory::Int) where T
    # Jacobians
    dynamics_constraints = constraint_indices(dynamics,
        shift=0)
    dynamics_jacobians = jacobian_indices(dynamics,
        shift=0)
    stage_constraints = constraint_indices(constraints,
        shift=num_constraint(dynamics))
    stage_jacobians = jacobian_indices(constraints,
        shift=num_jacobian(dynamics))
    general_constraint = constraint_indices(general,
        shift=(num_constraint(dynamics) + num_constraint(constraints)))
    general_jacobian = jacobian_indices(general,
        shift=(num_jacobian(dynamics) + num_jacobian(constraints)))

    # Hessian of Lagrangian
    objective_hessians = hessian_indices(objective, key, num_state, num_action)
    dynamics_hessians = hessian_indices(dynamics, key, num_state, num_action)
    stage_hessians = hessian_indices(constraints, key, num_state, num_action)
    general_hessian = hessian_indices(general, key, num_trajectory)

    # indices
    x_idx = state_indices(dynamics)
    u_idx = action_indices(dynamics)
    xu_idx = state_action_indices(dynamics)
    xuy_idx = state_action_next_state_indices(dynamics)

    return TrajectoryOptimizationIndices(
        objective_hessians,
        dynamics_constraints,
        dynamics_jacobians,
        dynamics_hessians,
        stage_constraints,
        stage_jacobians,
        stage_hessians,
        general_constraint,
        general_jacobian,
        general_hessian,
        x_idx,
        u_idx,
        xu_idx,
        xuy_idx)
end

struct NLPData{T} <: MOI.AbstractNLPEvaluator
    trajopt::TrajectoryOptimizationData{T}
    num_variables::Int
    num_constraint::Int
    num_jacobian::Int
    num_hessian_lagrangian::Int
    variable_bounds::Vector{Vector{T}}
    constraint_bounds::Vector{Vector{T}}
    indices::TrajectoryOptimizationIndices
    jacobian_sparsity
    hessian_lagrangian_sparsity
    hessian_lagrangian::Bool
    general_constraint::GeneralConstraint{T}
    parameters::Vector{T}
    duals_general::Vector{T}
end

function primal_bounds(bounds::Bounds{T}, num_variables::Int, state_indices::Vector{Vector{Int}}, action_indices::Vector{Vector{Int}}) where T
    lower, upper = -Inf * ones(num_variables), Inf * ones(num_variables)

    for (t, bnd) in enumerate(bounds)
        length(bnd.state_lower)  > 0 && (lower[state_indices[t]]  = bnd.state_lower)
        length(bnd.state_upper)  > 0 && (upper[state_indices[t]]  = bnd.state_upper)
        length(bnd.action_lower) > 0 && (lower[action_indices[t]] = bnd.action_lower)
        length(bnd.action_upper) > 0 && (upper[action_indices[t]] = bnd.action_upper)
    end
    return lower, upper
end

function constraint_bounds(constraints::Constraints{T}, general::GeneralConstraint{T},
    num_dynamics::Int, num_stage::Int, indices::TrajectoryOptimizationIndices) where T
    # total constraints
    total_constraint = num_dynamics + num_stage + general.num_constraint
    # bounds
    lower, upper = zeros(total_constraint), zeros(total_constraint)
    # stage
    for (t, con) in enumerate(constraints)
        lower[indices.stage_constraints[t][con.indices_inequality]] .= -Inf
    end
    # general
    lower[collect(num_dynamics + num_stage .+ general.indices_inequality)] .= -Inf
    return lower, upper
end

function NLPData(trajopt::TrajectoryOptimizationData;
    evaluate_hessian=false,
    general_constraint=GeneralConstraint())

    # number of variables
    total_variables = sum(trajopt.state_dimensions) + sum(trajopt.action_dimensions)

    # number of constraints
    num_dynamics = num_constraint(trajopt.dynamics)
    num_stage = num_constraint(trajopt.constraints)
    num_general = num_constraint(general_constraint)
    total_constraints = num_dynamics + num_stage + num_general

    # number of nonzeros in constraint Jacobian
    num_dynamics_jacobian = num_jacobian(trajopt.dynamics)
    num_constraint_jacobian = num_jacobian(trajopt.constraints)
    num_general_jacobian = num_jacobian(general_constraint)
    total_jacobians = num_dynamics_jacobian + num_constraint_jacobian + num_general_jacobian

    # constraint Jacobian sparsity
    sparsity_dynamics = sparsity_jacobian(trajopt.dynamics, trajopt.state_dimensions, trajopt.action_dimensions,
        row_shift=0)
    sparsity_constraint = sparsity_jacobian(trajopt.constraints, trajopt.state_dimensions, trajopt.action_dimensions,
        row_shift=num_dynamics)
    sparsity_general = sparsity_jacobian(general_constraint, total_variables, row_shift=(num_dynamics + num_stage))
    jacobian_sparsity = collect([sparsity_dynamics..., sparsity_constraint..., sparsity_general...])

    # Hessian of Lagrangian sparsity
    sparsity_objective_hessians = sparsity_hessian(trajopt.objective, trajopt.state_dimensions, trajopt.action_dimensions)
    sparsity_dynamics_hessians = sparsity_hessian(trajopt.dynamics, trajopt.state_dimensions, trajopt.action_dimensions)
    sparsity_constraint_hessian = sparsity_hessian(trajopt.constraints, trajopt.state_dimensions, trajopt.action_dimensions)
    sparsity_general_hessian = sparsity_hessian(general_constraint, total_variables)
    hessian_lagrangian_sparsity = [sparsity_objective_hessians..., sparsity_dynamics_hessians..., sparsity_constraint_hessian..., sparsity_general_hessian...]
    hessian_lagrangian_sparsity = !isempty(hessian_lagrangian_sparsity) ? hessian_lagrangian_sparsity : Tuple{Int,Int}[]
    hessian_sparsity_key = sort(unique(hessian_lagrangian_sparsity))

    # number of nonzeros in Hessian of Lagrangian
    num_hessian_lagrangian = length(hessian_lagrangian_sparsity)

    # indices
    idx = indices(
        trajopt.objective,
        trajopt.dynamics,
        trajopt.constraints,
        general_constraint,
        hessian_sparsity_key,
        trajopt.state_dimensions,
        trajopt.action_dimensions,
        total_variables)

    # primal variable bounds
    primal_lower, primal_upper = primal_bounds(trajopt.bounds, total_variables, idx.states, idx.actions)

    # nonlinear constraint bounds
    constraint_lower, constraint_upper = constraint_bounds(trajopt.constraints, general_constraint, num_dynamics, num_stage, idx)

    NLPData(trajopt,
        total_variables,
        total_constraints,
        total_jacobians,
        num_hessian_lagrangian,
        [primal_lower, primal_upper],
        [constraint_lower, constraint_upper],
        idx,
        jacobian_sparsity,
        hessian_sparsity_key,
        evaluate_hessian,
        general_constraint,
        vcat(trajopt.parameters...),
        zeros(general_constraint.num_constraint))
end

struct SolverData
    nlp_bounds::Vector{MOI.NLPBoundsPair}
    block_data::MOI.NLPBlockData
    optimizer::Ipopt.Optimizer
    variables::Vector{MOI.VariableIndex}
end

function SolverData(nlp::NLPData;
    options=Options())

    # solver
    nlp_bounds = MOI.NLPBoundsPair.(nlp.constraint_bounds...)
    block_data = MOI.NLPBlockData(nlp_bounds, nlp, true)

    # instantiate NLP solver
    optimizer = Ipopt.Optimizer()

    # set NLP optimizer options
    for name in fieldnames(typeof(options))
        optimizer.options[String(name)] = getfield(options, name)
    end

    z = MOI.add_variables(optimizer, nlp.num_variables)

    for i = 1:nlp.num_variables
        MOI.add_constraint(optimizer, z[i], MOI.LessThan(nlp.variable_bounds[2][i]))
        MOI.add_constraint(optimizer, z[i], MOI.GreaterThan(nlp.variable_bounds[1][i]))
    end

    MOI.set(optimizer, MOI.NLPBlock(), block_data)
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    return SolverData(nlp_bounds, block_data, optimizer, z)
end

function trajectory!(states::Vector{Vector{T}}, actions::Vector{Vector{T}},
    trajectory::Vector{T},
    state_indices::Vector{Vector{Int}}, action_indices::Vector{Vector{Int}}) where T
    for (t, idx) in enumerate(state_indices)
        states[t] .= @views trajectory[idx]
    end
    for (t, idx) in enumerate(action_indices)
        actions[t] .= @views trajectory[idx]
    end
end

function duals!(duals_dynamics::Vector{Vector{T}}, duals_constraints::Vector{Vector{T}}, duals_general::Vector{T}, duals,
    dynamics_indices::Vector{Vector{Int}}, constraint_indices::Vector{Vector{Int}}, general_indices::Vector{Int}) where T
    for (t, idx) in enumerate(dynamics_indices)
        duals_dynamics[t] .= @views duals[idx]
    end
    for (t, idx) in enumerate(constraint_indices)
        duals_constraints[t] .= @views duals[idx]
    end
    duals_general .= duals[general_indices]
end


text:
@testset "blockinv" begin
    M = [1 2 3;
         3 5 6;
         7 8 9]
    @test blockinv(inv(M), [3]) ≈ [1/9]

    M = [1 2 12 1;
         5 6 7 3;
         9 10 11 12;
         13 14 15 16];
    @test blockinv(inv(M), [3,4]) ≈ inv(M[[3,4], [3,4]])
    @test blockinv(inv(M), [4,3]) ≈ inv(M[[4,3], [4,3]])

    n = 6
    i = Int(n/2) + 1 : n
    for rep = 1 : 10
        M = rand(n, n)
        if abs(det(M)) > 1e-6 && abs(det(M[i,i])) > 1e-6
            if !(blockinv(inv(M), i) ≈ inv(M[i,i]))
                @show M
            end
            @test blockinv(inv(M), i) ≈ inv(M[i,i])
        end
    end

    n = 6
    i = Int(n/2) + 1 : n
    for rep = 1 : 10
        M = rand(n, n)
        M = Hermitian(0.5M + 0.5M')
        if abs(det(M)) > 1e-6 && abs(det(M[i,i])) > 1e-6
            if !(blockinv(inv(M), i) ≈ inv(M[i,i]))
                @show M
            end
            @test blockinv(inv(M), i) ≈ inv(M[i,i])
            @test blockinv(inv(M), i) isa Hermitian
        end
    end
end


text:
using Test
using LinearAlgebra
using QuantumControlBase: MinimalObjective
using QuantumControlBase.Functionals
using QuantumControlBase.WeylChamber
using QuantumControlBase.TestUtils

⊗ = kron
const 𝕚 = 1im

function ket(i::Int64; N)
    Ψ = zeros(ComplexF64, N)
    Ψ[i+1] = 1
    return Ψ
end

function ket(i::Int64; N)
    Ψ = zeros(ComplexF64, N)
    Ψ[i+1] = 1
    return Ψ
end

function ket(indices::Int64...; N)
    Ψ = ket(indices[1]; N=N)
    for i in indices[2:end]
        Ψ = Ψ ⊗ ket(i; N=N)
    end
    return Ψ
end

function ket(label::AbstractString; N)
    indices = [parse(Int64, digit) for digit in label]
    return ket(indices...; N=N)
end

const CNOT = [
    1  0  0  0
    0  1  0  0
    0  0  0  1
    0  0  1  0
]

const CPHASE8 = [
    1  0  0  0
    0  1  0  0
    0  0  1  0
    0  0  0  exp(𝕚 * π / 8)
]

const sqrt_SWAP = [
    1      0        0      0
    0   1/2-𝕚/2  1/2+𝕚/2   0
    0   1/2+𝕚/2  1/2-𝕚/2   0
    0      0        0      1
]

const SWAP = sqrt_SWAP^2

@testset "find-gate" begin

    # This was originally going to test a find_gate routine, but that routine
    # turned out to be a very intuitive one-liner, so it was better not to
    # define it. Still, it's a good idea to test out the pattern of obtaining
    # transformed basis states and reconstructing the gate from the transformed
    # basis states.
    #
    # In particular, getting the transformed basis states is actually a little
    # bit un-intuitive: it is a matrix-vector product with the transpose gate!

    N = 10
    basis = [ket("00"; N), ket("01"; N), ket("10"; N), ket("11"; N)]

    #! format: off
    U = [
        #  00 01 10 11 :j
        1  0  0  0  # i: 00
        0  0  𝕚  0  # i: 01   # =  𝕚 |01⟩⟨10|  # U_01_10 = ⟨01|U|10⟩
        0 -𝕚  0  0  # i: 10   # = -𝕚 |10⟩⟨01|
        0  0  0  1  # i: 11
        # columns are image vectors
    ]
    #! format: on
    @test norm(U * U' - I) < 1e-15
    U_basis = [ket("00"; N), -𝕚 * ket("10"; N), 𝕚 * ket("01"; N), ket("11"; N)]
    U_basis2 = transpose(U) * basis
    @test maximum(norm.(U_basis - U_basis2)) < 1e-15
    gate = [basis[i] ⋅ U_basis[j] for i = 1:4, j = 1:4]
    @test norm(gate - U) < 1e-14

    U = exp(1im * random_hermitian_matrix(4, 1.0))
    U_basis = transpose(U) * basis
    gate = [basis[i] ⋅ U_basis[j] for i = 1:4, j = 1:4]
    @test norm(gate - U) < 1e-14

end

@testset "Local Invariants" begin

    g1, g2, g3 = local_invariants(CNOT)
    @test g1 ≈ 0 atol = 1e-15
    @test g2 ≈ 0 atol = 1e-15
    @test g3 ≈ 1 atol = 1e-15

    g1, g2, g3 = local_invariants(CPHASE8)
    @test g1 ≈ cos(0.5 * π / 8)^2 atol = 1e-15
    @test g2 ≈ 0 atol = 1e-15
    @test g3 ≈ (1 + 2 * g1) atol = 1e-15

    g1, g2, g3 = local_invariants(sqrt_SWAP)
    @test g1 ≈ 0 atol = 1e-15
    @test g2 ≈ (1 / 4) atol = 1e-15
    @test g3 ≈ 0 atol = 1e-15

    g1, g2, g3 = local_invariants(sqrt_SWAP')
    @test g1 ≈ 0 atol = 1e-15
    @test g2 ≈ (-1 / 4) atol = 1e-15
    @test g3 ≈ 0 atol = 1e-15

    g1, g2, g3 = local_invariants(SWAP)
    @test g1 ≈ -1 atol = 1e-15
    @test g2 ≈ 0 atol = 1e-15
    @test g3 ≈ -3 atol = 1e-15

end

@testset "Weyl Chamber Coordinates" begin

    c1, c2, c3 = weyl_chamber_coordinates(CNOT)
    @test c1 ≈ 0.5 atol = 1e-15
    @test c2 ≈ 0 atol = 1e-15
    @test c3 ≈ 0 atol = 1e-15
    Uc = canonical_gate(c1, c2, c3)
    @test norm(collect(weyl_chamber_coordinates(Uc)) - [c1, c2, c3]) < 1e-14

    c1, c2, c3 = weyl_chamber_coordinates(CPHASE8)
    @test c1 ≈ (1 / 16) atol = 1e-15  #  = γ / 2   for   γ = π/8
    @test c2 ≈ 0 atol = 1e-15
    @test c3 ≈ 0 atol = 1e-15
    # the Weyl chambers coordinates of the canonical gate are unstable: The
    # bottom surfaces of the W0 and W0* regions are equivalent, and floating
    # point errors can push π/8 to the equivalent 1-π/8

    c1, c2, c3 = weyl_chamber_coordinates(sqrt_SWAP)
    @test c1 ≈ (1 / 4) atol = 1e-15
    @test c2 ≈ (1 / 4) atol = 1e-15
    @test c3 ≈ (1 / 4) atol = 1e-15
    Uc = canonical_gate(c1, c2, c3)
    @test norm(collect(weyl_chamber_coordinates(Uc)) - [c1, c2, c3]) < 1e-14

    c1, c2, c3 = weyl_chamber_coordinates(sqrt_SWAP')
    @test c1 ≈ (3 / 4) atol = 1e-15
    @test c2 ≈ (1 / 4) atol = 1e-15
    @test c3 ≈ (1 / 4) atol = 1e-15
    Uc = canonical_gate(c1, c2, c3)
    @test norm(collect(weyl_chamber_coordinates(Uc)) - [c1, c2, c3]) < 1e-14

    c1, c2, c3 = weyl_chamber_coordinates(SWAP)
    @test c1 ≈ (1 / 2) atol = 1e-15
    @test c2 ≈ (1 / 2) atol = 1e-15
    @test c3 ≈ (1 / 2) atol = 1e-15
    Uc = canonical_gate(c1, c2, c3)
    @test norm(collect(weyl_chamber_coordinates(Uc)) - [c1, c2, c3]) < 1e-14

end

@testset "Gate Concurrence" begin
    @test gate_concurrence(CNOT) ≈ 1.0
    @test gate_concurrence(sqrt_SWAP) ≈ 1.0
    @test gate_concurrence(sqrt_SWAP') ≈ 1.0
    @test gate_concurrence(CPHASE8) ≈ abs(sin(0.5 * π / 8))
    @test gate_concurrence(SWAP) ≈ 0.0 atol = 1e-15

    @test gate_concurrence(0.2, 0.1, 0.1) < 1.0
    @test gate_concurrence(0.8, 0.1, 0.1) < 1.0
    @test gate_concurrence(0.5, 0.4, 0.1) ≈ 1.0
    @test gate_concurrence(0.5, 0.4, 0.2) < 1.0

end

@testset "PE Evaluation" begin

    @test D_PE(CNOT) ≈ 0.0
    @test D_PE(SWAP) ≈ -2.0
    @test D_PE(SWAP; absolute_square=true) ≈ 4.0
    @test D_PE(I(4)) ≈ 2.0
    @test D_PE(I(4); absolute_square=true) ≈ 4.0

end

@testset "Weyl Chamber Region" begin

    @test in_weyl_chamber(CNOT; region="W")
    @test in_weyl_chamber(CNOT; region="PE")
    @test weyl_chamber_region(CNOT) == "PE"
    @test !in_weyl_chamber(CNOT; region="W0")
    @test !in_weyl_chamber(CNOT; region="W0*")
    @test !in_weyl_chamber(CNOT; region="W1")

    @test in_weyl_chamber(I(4); region="W0")

    @test in_weyl_chamber(SWAP; region="W1")
    @test weyl_chamber_region(SWAP) == "W1"

    @test in_weyl_chamber(CPHASE8; region="W0")

    @test in_weyl_chamber(0.2, 0.1, 0.1)
    @test in_weyl_chamber(0.2, 0.1, 0.1; region="W0")
    @test !in_weyl_chamber(0.2, 0.1, 0.1; region="W1")

    @test in_weyl_chamber(0.8, 0.1, 0.1; region="W0*")

    @test in_weyl_chamber(0.5, 0.4, 0.4; region="W1")

    @test !in_weyl_chamber(1.1, 0.0, 0.0)
    @test !in_weyl_chamber(0.0, 0.6, 0.0)
    @test !in_weyl_chamber(0.0, 0.0, 0.6)
    @test !in_weyl_chamber(-0.1, 0.0, 0.0)
    @test !in_weyl_chamber(0.0, -0.1, 0.0)
    @test !in_weyl_chamber(0.0, 0.0, -0.1)
    @test !in_weyl_chamber(1.0, 0.1, 0.0)
    @test !in_weyl_chamber(0.6, 0.45, 0.0)
    @test !in_weyl_chamber(0.5, 0.0, 0.1)

    @test weyl_chamber_region(0.0, 0.0, -0.1) == ""

end

@testset "PE chi states" begin

    Ψ1 = random_state_vector(4)
    Ψ2 = random_state_vector(4)
    Ψ3 = random_state_vector(4)
    Ψ4 = random_state_vector(4)

    objectives = [
        MinimalObjective(
            initial_state=ket(lbl; N=2),
            generator=random_hermitian_real_matrix(4, 1.0)
        ) for lbl in ("00", "01", "10", "11")
    ]

    J_T = gate_functional(D_PE; unitarity_weight=0.0)
    chi_pe! = make_chi(J_T, objectives)
    ϕ = [Ψ1, Ψ2, Ψ3, Ψ4]
    χ_zygote = [similar(ϕₖ) for ϕₖ in ϕ]
    chi_pe!(χ_zygote, ϕ, objectives)

    gate_chi_pe! = make_gate_chi(D_PE, objectives; unitarity_weight=0.0)
    χ_gate_zygote = [similar(ϕₖ) for ϕₖ in ϕ]
    gate_chi_pe!(χ_gate_zygote, ϕ, objectives)

    gate_chi_pe_fd! =
        make_gate_chi(D_PE, objectives; use_finite_differences=true, unitarity_weight=0.0)
    χ_gate_fd = [similar(ϕₖ) for ϕₖ in ϕ]
    gate_chi_pe_fd!(χ_gate_fd, ϕ, objectives)

    vec_angle(v⃗, w⃗) = acos((v⃗ ⋅ w⃗) / (norm(v⃗) * norm(w⃗)))

    # Does Zygote gate gradient and FD gate gradient match?
    @test norm(χ_gate_zygote[1] - χ_gate_fd[1]) < 1e-8
    @test norm(χ_gate_zygote[2] - χ_gate_fd[2]) < 1e-8
    @test norm(χ_gate_zygote[3] - χ_gate_fd[3]) < 1e-8
    @test norm(χ_gate_zygote[4] - χ_gate_fd[4]) < 1e-8

    # Does the Zygote gate gradient and Zygote full gradient match?
    @test norm(χ_gate_zygote[1] - χ_zygote[1]) < 1e-8
    @test norm(χ_gate_zygote[2] - χ_zygote[2]) < 1e-8
    @test norm(χ_gate_zygote[3] - χ_zygote[3]) < 1e-8
    @test norm(χ_gate_zygote[4] - χ_zygote[4]) < 1e-8

    # Does the Zygote full gradient and the FD gate gradient match?
    @test norm(χ_zygote[1] - χ_gate_fd[1]) < 1e-8
    @test norm(χ_zygote[2] - χ_gate_fd[2]) < 1e-8
    @test norm(χ_zygote[3] - χ_gate_fd[3]) < 1e-8
    @test norm(χ_zygote[4] - χ_gate_fd[4]) < 1e-8

end

@testset "Concurrence-chi ($region)" for region in ("W0", "W0*", "W1")

    # Note: the χ states are known to fail for the PE region (but they'd be
    # zero anyway)

    J_T_val = 0.0
    local J_T_U
    local J_T
    local ϕ

    basis = [ket(lbl; N=2) for lbl in ("00", "01", "10", "11")]
    H = random_hermitian_real_matrix(4, 1.0)

    objectives = [MinimalObjective(initial_state=Ψ, generator=H) for Ψ in basis]

    c1, c2, c3 = rand(3) .* [1, 0.5, 0.5]
    while !in_weyl_chamber(c1, c2, c3; region)
        c1, c2, c3 = rand(3) .* [1, 0.5, 0.5]
    end

    U = canonical_gate(c1, c2, c3)
    ϕ = U' * basis

    J_T_U = U -> 1 - gate_concurrence(U)
    J_T = gate_functional(J_T_U)

    gate_chi_pe! = make_gate_chi(J_T_U, objectives)
    χ_gate_zygote = [similar(ϕₖ) for ϕₖ in ϕ]
    gate_chi_pe!(χ_gate_zygote, ϕ, objectives)

    gate_chi_pe_fd! = make_gate_chi(J_T_U, objectives; use_finite_differences=true)
    χ_gate_fd = [similar(ϕₖ) for ϕₖ in ϕ]
    gate_chi_pe_fd!(χ_gate_fd, ϕ, objectives)

    # Does Zygote gate gradient and FD gate gradient match?
    @test norm(χ_gate_zygote[1] - χ_gate_fd[1]) < 1e-8
    @test norm(χ_gate_zygote[2] - χ_gate_fd[2]) < 1e-8
    @test norm(χ_gate_zygote[3] - χ_gate_fd[3]) < 1e-8
    @test norm(χ_gate_zygote[4] - χ_gate_fd[4]) < 1e-8

end


text:
import GLFW  # we focus exclusively on the GLFW approach
using ModernGL

# The GLFW package calls GLFW.Init() automatically when loaded, so we
# don't need to call it explicitly.

# Specify minimum versions
GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, 3)  # minimum OpenGL v. 3
GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, 2)  # minimum OpenGL v. 3.2
GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE)
GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE)

# Create the window
GLFW.WindowHint(GLFW.RESIZABLE, GL_FALSE)
window = GLFW.CreateWindow(800, 600, "Context creation")  # windowed
# window = GLFW.CreateWindow(800, 600, "Context creation", GLFW.GetPrimaryMonitor())  # fullscreen
GLFW.MakeContextCurrent(window)

# This is a touch added from
#   http://www.opengl-tutorial.org/beginners-tutorials/tutorial-1-opening-a-window/:
# Retain keypress events until the next call to GLFW.GetKey, even if
# the key has been released in the meantime
GLFW.SetInputMode(window, GLFW.STICKY_KEYS, GL_TRUE)

# Draw nothing, while waiting for a close event
while !GLFW.WindowShouldClose(window)
    GLFW.SwapBuffers(window)
    GLFW.PollEvents()
    if GLFW.GetKey(window, GLFW.KEY_ESCAPE) == GLFW.PRESS
        GLFW.SetWindowShouldClose(window, true)
    end
end

GLFW.DestroyWindow(window)  # needed if you're running this from the REPL

# GLFW.Terminate() is called automatically upon exit


text:
#=
There are exactly ten ways of selecting three from five, 12345:

123, 124, 125, 134, 135, 145, 234, 235, 245, and 345

In combinatorics, we use the notation, 5C3 = 10.

In general,

nCr =
n!
r!(n−r)!
,where r ≤ n, n! = n×(n−1)×...×3×2×1, and 0! = 1.
It is not until n = 23, that a value exceeds one-million: 23C10 = 1144066.

How many, not necessarily distinct, values of  nCr, for 1 ≤ n ≤ 100, are greater than one-million?
=#
function combine(n, r)
  factorial(n)/(factorial(r)factorial(n-r))
end

function calc()
  sum = 0
  for n in 1:100, r in reverse(1:n)
    if combine(BigInt(n), BigInt(r)) > 1_000_000
      sum = sum + 1
    end
  end
  sum
end
@time println(calc())

max_stars_repo_name: daniel-beard/JuliaProjectEuler

max_issues_repo_name: daniel-beard/JuliaProjectEuler

max_forks_repo_name: daniel-beard/JuliaProjectEuler