text:
#  Copyright 2016, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.
#############################################################################
# JuMP
# An algebraic modelling langauge for Julia
# See http://github.com/JuliaOpt/JuMP.jl
#############################################################################
# test/model.jl
# Testing Model printing, basic solving
# Must be run as part of runtests.jl, as it needs a list of solvers.
#############################################################################
using JuMP, FactCheck

# If solvers not loaded, load them (i.e running just these tests)
!isdefined(:lp_solvers) && include("solvers.jl")

# To ensure the tests work on Windows and Linux/OSX, we need
# to use the correct comparison operators
const leq = JuMP.repl[:leq]
const geq = JuMP.repl[:geq]
const  eq = JuMP.repl[:eq]

const TOL = 1e-4

modPath = joinpath(dirname(@__FILE__), "mod")

facts("[model] Check error cases") do
    @fact_throws Model(solver=:Foo)
    modErr = Model()
    @fact_throws setobjectivesense(modErr, :Maximum)
    @variable(modErr, errVar)
    @fact getvalue(errVar) --> isnan
    @fact_throws getdual(errVar)

    modErr = Model()
    @variable(modErr, x, Bin)
    @objective(modErr, Max, x)
    con = @constraint(modErr, x <= 0.5)
    solve(modErr)
    @fact_throws getdual(con)

    modErr = Model()
    @variable(modErr, 0 <= x <= 1)
    @objective(modErr, Max, x)
    @constraint(modErr, x <= -1)
    solve(modErr, suppress_warnings=true)
    @fact getvalue(x) --> isnan
end

facts("[model] Performance warnings") do
    m = Model()
    @variable(m, x[1:2], start=0)
    for i in 1:500
       getvalue(x)
    end
    q = 0
    for i in 1:30000
       q += 3x[1]
    end
end

facts("[model] Test printing a model") do
    modA = Model()
    x = Variable(modA, 0, Inf, :Cont)
    @variable(modA, y <= 5, Int)
    @variable(modA, 2 <= z <= 4)
    @variable(modA, 0 <= r[i=3:6] <= i)
    @objective(modA, Max, ((x + y)/2.0 + 3.0)/3.0 + z + r[3])
    @constraintref constraints[1:3]
    constraints[1] = @constraint(modA, 2 <= x+y <= 4)
    constraints[2] = @constraint(modA, sum{r[i],i=3:5} <= (2 - x)/2.0)
    constraints[3] = @constraint(modA, 6y + y <= z + r[6]/1.9)
    #####################################################################
    # Test LP writer
    writeLP(modA, modPath * "A.lp")
    modALP = if VERSION >= v"0.5.0-dev+1866" # leading zero, base Julia PR #14377
        String[
        "Maximize",
        "obj: 0.16666666666666666 VAR1 + 0.16666666666666666 VAR2 + 1 VAR3 + 1 VAR4",
        "Subject To",
        "c1: 1 VAR1 + 1 VAR2 >= 2",
        "c2: 1 VAR1 + 1 VAR2 <= 4",
        "c3: 1 VAR4 + 1 VAR5 + 1 VAR6 + 0.5 VAR1 <= 1",
        "c4: 6 VAR2 + 1 VAR2 - 1 VAR3 - 0.5263157894736842 VAR7 <= 0",
        "Bounds",
        "0 <= VAR1 <= +inf",
        "-inf <= VAR2 <= 5",
        "2 <= VAR3 <= 4",
        "0 <= VAR4 <= 3",
        "0 <= VAR5 <= 4",
        "0 <= VAR6 <= 5",
        "0 <= VAR7 <= 6",
        "General",
        "VAR2",
        "End"]
    else
        String[
        "Maximize",
        "obj: .16666666666666666 VAR1 + .16666666666666666 VAR2 + 1 VAR3 + 1 VAR4",
        "Subject To",
        "c1: 1 VAR1 + 1 VAR2 >= 2",
        "c2: 1 VAR1 + 1 VAR2 <= 4",
        "c3: 1 VAR4 + 1 VAR5 + 1 VAR6 + .5 VAR1 <= 1",
        "c4: 6 VAR2 + 1 VAR2 - 1 VAR3 - .5263157894736842 VAR7 <= 0",
        "Bounds",
        "0 <= VAR1 <= +inf",
        "-inf <= VAR2 <= 5",
        "2 <= VAR3 <= 4",
        "0 <= VAR4 <= 3",
        "0 <= VAR5 <= 4",
        "0 <= VAR6 <= 5",
        "0 <= VAR7 <= 6",
        "General",
        "VAR2",
        "End"]
    end
    modAfp = open(modPath * "A.lp")
    lineInd = 1
    while !eof(modAfp)
        line = readline(modAfp)
        @fact strip(line) --> strip(modALP[lineInd])
        lineInd += 1
    end
    close(modAfp)
    #####################################################################
    # Test MPS writer
    writeMPS(modA, modPath * "A.mps")
    if VERSION >= v"0.5.0-dev+1866" # leading zero, base Julia PR #14377
        modAMPS = String[
        "NAME   JuMPModel",
        "ROWS",
        " N  OBJ",
        " E  CON1",
        " L  CON2",
        " L  CON3",
        "COLUMNS",
        "    VAR1  CON1  1",
        "    VAR1  CON2  0.5",
        "    VAR1  OBJ  -0.16666666666666666",
        "    MARKER    'MARKER'                 'INTORG'",
        "    VAR2  CON1  1",
        "    VAR2  CON3  7",
        "    VAR2  OBJ  -0.16666666666666666",
        "    MARKER    'MARKER'                 'INTEND'",
        "    VAR3  CON3  -1",
        "    VAR3  OBJ  -1",
        "    VAR4  CON2  1",
        "    VAR4  OBJ  -1",
        "    VAR5  CON2  1",
        "    VAR5  OBJ  0",
        "    VAR6  CON2  1",
        "    VAR6  OBJ  0",
        "    VAR7  CON3  -0.5263157894736842",
        "    VAR7  OBJ  0",
        "RHS",
        "    rhs    CON1    2",
        "    rhs    CON2    1",
        "    rhs    CON3    0",
        "RANGES",
        "    rhs    CON1    2",
        "BOUNDS",
        "  PL BOUND VAR1",
        "  MI BOUND VAR2",
        "  UP BOUND VAR2 5",
        "  LO BOUND VAR3 2",
        "  UP BOUND VAR3 4",
        "  UP BOUND VAR4 3",
        "  UP BOUND VAR5 4",
        "  UP BOUND VAR6 5",
        "  UP BOUND VAR7 6",
        "ENDATA"]
    else
        modAMPS = String[
        "NAME   JuMPModel",
        "ROWS",
        " N  OBJ",
        " E  CON1",
        " L  CON2",
        " L  CON3",
        "COLUMNS",
        "    VAR1  CON1  1",
        "    VAR1  CON2  .5",
        "    VAR1  OBJ  -.16666666666666666",
        "    MARKER    'MARKER'                 'INTORG'",
        "    VAR2  CON1  1",
        "    VAR2  CON3  7",
        "    VAR2  OBJ  -.16666666666666666",
        "    MARKER    'MARKER'                 'INTEND'",
        "    VAR3  CON3  -1",
        "    VAR3  OBJ  -1",
        "    VAR4  CON2  1",
        "    VAR4  OBJ  -1",
        "    VAR5  CON2  1",
        "    VAR5  OBJ  0",
        "    VAR6  CON2  1",
        "    VAR6  OBJ  0",
        "    VAR7  CON3  -.5263157894736842",
        "    VAR7  OBJ  0",
        "RHS",
        "    rhs    CON1    2",
        "    rhs    CON2    1",
        "    rhs    CON3    0",
        "RANGES",
        "    rhs    CON1    2",
        "BOUNDS",
        "  PL BOUND VAR1",
        "  MI BOUND VAR2",
        "  UP BOUND VAR2 5",
        "  LO BOUND VAR3 2",
        "  UP BOUND VAR3 4",
        "  UP BOUND VAR4 3",
        "  UP BOUND VAR5 4",
        "  UP BOUND VAR6 5",
        "  UP BOUND VAR7 6",
        "ENDATA"]
    end
    modAfp = open(modPath * "A.mps")
    lineInd = 1
    while !eof(modAfp)
        line = readline(modAfp)
        @fact chomp(line) --> modAMPS[lineInd]
        lineInd += 1
    end
    close(modAfp)

    # Getter/setters
    @fact MathProgBase.numvar(modA) --> 7
    @fact MathProgBase.numlinconstr(modA) --> 3
    @fact MathProgBase.numquadconstr(modA) --> 0
    @fact MathProgBase.numconstr(modA) --> 3
    @fact getobjectivesense(modA) --> :Max
    setobjectivesense(modA, :Min)
    @fact getobjectivesense(modA) --> :Min
end

facts("[model] Quadratic MPS writer") do
    modQ = Model()
    @variable(modQ, x == 1)
    @variable(modQ, y ≥ 0)
    @objective(modQ, Min, x^2 - 2*x*y + y^2 + x)

    #####################################################################
    # Test MPS writer
    writeMPS(modQ, modPath * "Q.mps")
    modQMPS = String[
    "NAME   JuMPModel",
    "ROWS",
    " N  OBJ",
    "COLUMNS",
    "    VAR1  OBJ  1",
    "    VAR2  OBJ  0",
    "RHS",
    "BOUNDS",
    "  LO BOUND VAR1 1",
    "  UP BOUND VAR1 1",
    "  PL BOUND VAR2",
    "QMATRIX",
    "  VAR1 VAR1  2",
    "  VAR1 VAR2 -2",
    "  VAR2 VAR2  2",
    "ENDATA"]
    modQfp = open(modPath * "Q.mps")
    lineInd = 1
    while !eof(modQfp)
        line = readline(modQfp)
        @fact chomp(line) --> modQMPS[lineInd]
        lineInd += 1
    end
    close(modQfp)

    # Getter/setters
    @fact MathProgBase.numvar(modQ) --> 2
    @fact MathProgBase.numlinconstr(modQ) --> 0
    @fact MathProgBase.numquadconstr(modQ) --> 0
    @fact MathProgBase.numconstr(modQ) --> 0
    @fact getobjectivesense(modQ) --> :Min
end

facts("[model] Test solving a MILP") do
for solver in ip_solvers
context("With solver $(typeof(solver))") do
    modA = Model(solver=solver)
    @variable(modA, x >= 0)
    @variable(modA, y <= 5, Int)
    @variable(modA, 2 <= z <= 4)
    @variable(modA, 0 <= r[i=3:6] <= i)
    @objective(modA, Max, ((x + y)/2.0 + 3.0)/3.0 + z + r[3])
    @constraint(modA, 2 <= x+y)
    @constraint(modA, sum{r[i],i=3:5} <= (2 - x)/2.0)
    @constraint(modA, 7.0*y <= z + r[6]/1.9)

    @fact solve(modA)       --> :Optimal
    @fact modA.objVal       --> roughly(1.0+4.833334, TOL)
    @fact getvalue(x)       --> roughly(1.0, TOL)
    @fact getvalue(y)       --> roughly(1.0, TOL)
    @fact getvalue(z)       --> roughly(4.0, TOL)
    @fact getvalue(r)[3]    --> roughly(0.5, TOL)
    @fact getvalue(r)[4]    --> roughly(0.0, TOL)
    @fact getvalue(r)[5]    --> roughly(0.0, TOL)
    @fact getvalue(r)[6]    --> greater_than_or_equal(5.7-TOL)
    @fact getvalue(r)[6]    --> less_than_or_equal(6.0+TOL)
    @fact getobjective(modA).aff --> ((x + y)/2.0 + 3.0)/3.0 + z + r[3]
end # solver context
end # loop over solvers
end # facts block

facts("[model] Test solving an LP (Min)") do
for solver in lp_solvers
context("With solver $(typeof(solver))") do
    modA = Model(solver=solver)
    @variable(modA, x >= 0)
    @variable(modA, y <= 5)
    @variable(modA, 2 <= z <= 4)
    @variable(modA, 0 <= r[i=3:6] <= i)
    @objective(modA, Min, -((x + y)/2.0 + 3.0)/3.0 - z - r[3])
    @constraintref cons[1:3]
    cons[1] = @constraint(modA, x+y >= 2)
    cons[2] = @constraint(modA, sum{r[i],i=3:5} <= (2 - x)/2.0)
    cons[3] = @constraint(modA, 7.0*y <= z + r[6]/1.9)

    # Solution
    @fact solve(modA) --> :Optimal
    @fact getobjectivevalue(modA) --> roughly(-5.8446115, TOL)
    @fact getvalue(x)       --> roughly(0.9774436, TOL)
    @fact getvalue(y)       --> roughly(1.0225563, TOL)
    @fact getvalue(z)       --> roughly(4.0, TOL)
    @fact getvalue(r)[3]    --> roughly(0.5112781, TOL)
    @fact getvalue(r)[4]    --> roughly(0.0, TOL)
    @fact getvalue(r)[5]    --> roughly(0.0, TOL)
    @fact getvalue(r)[6]    --> roughly(6.0, TOL)

    # Reduced costs
    @fact getdual(x)    --> roughly( 0.0, TOL)
    @fact getdual(y)    --> roughly( 0.0, TOL)
    @fact getdual(z)    --> roughly(-1.0714286, TOL)
    @fact getdual(r)[3] --> roughly( 0.0, TOL)
    @fact getdual(r)[4] --> roughly(1.0, TOL)
    @fact getdual(r)[5] --> roughly(1.0, TOL)
    @fact getdual(r)[6] --> roughly(-0.03759398, TOL)

    # Row duals
    @fact getdual(cons)[1] --> roughly( 0.333333, TOL)
    @fact getdual(cons)[2] --> roughly(-1.0, TOL)
    @fact getdual(cons)[3] --> roughly(-0.0714286, TOL)
end # solver context
end # loop over solvers
end # facts block

facts("[model] Test solving an LP (Max)") do
for solver in lp_solvers
context("With solver $(typeof(solver))") do
    modA = Model(solver=solver)
    @variable(modA, x >= 0)
    @variable(modA, y <= 5)
    @variable(modA, 2 <= z <= 4)
    @variable(modA, 0 <= r[i=3:6] <= i)
    @objective(modA, Max, ((x + y)/2.0 + 3.0)/3.0 + z + r[3])
    @constraintref cons[1:3]
    cons[1] = @constraint(modA, x+y >= 2)
    cons[2] = @constraint(modA, sum{r[i],i=3:5} <= (2 - x)/2.0)
    cons[3] = @constraint(modA, 7.0*y <= z + r[6]/1.9)

    # Solution
    @fact solve(modA) --> :Optimal
    @fact getobjectivevalue(modA) --> roughly(5.8446115, TOL)
    @fact getvalue(x)       --> roughly(0.9774436, TOL)
    @fact getvalue(y)       --> roughly(1.0225563, TOL)
    @fact getvalue(z)       --> roughly(4.0, TOL)
    @fact getvalue(r)[3]    --> roughly(0.5112781, TOL)
    @fact getvalue(r)[4]    --> roughly(0.0, TOL)
    @fact getvalue(r)[5]    --> roughly(0.0, TOL)
    @fact getvalue(r)[6]    --> roughly(6.0, TOL)

    # Reduced costs
    @fact getdual(x)    --> roughly( 0.0, TOL)
    @fact getdual(y)    --> roughly( 0.0, TOL)
    @fact getdual(z)    --> roughly( 1.0714286, TOL)
    @fact getdual(r)[3] --> roughly( 0.0, TOL)
    @fact getdual(r)[4] --> roughly(-1.0, TOL)
    @fact getdual(r)[5] --> roughly(-1.0, TOL)
    @fact getdual(r)[6] --> roughly( 0.03759398, TOL)

    # Row duals
    @fact getdual(cons)[1] --> roughly(-0.333333, TOL)
    @fact getdual(cons)[2] --> roughly( 1.0, TOL)
    @fact getdual(cons)[3] --> roughly( 0.0714286, TOL)
end # solver context
end # loop over solvers
end # facts block

facts("[model] Test binary variable handling") do
for solver in ip_solvers
context("With solver $(typeof(solver))") do
    modB = Model(solver=solver)
    @variable(modB, x, Bin)
    @objective(modB, Max, x)
    @constraint(modB, x <= 10)
    status = solve(modB)
    @fact status --> :Optimal
    @fact getvalue(x) --> roughly(1.0, TOL)
end
end
end

facts("[model] Test model copying") do
    source = Model()
    @variable(source, 2 <= x <= 5)
    @variable(source, 0 <= y <= 1, Int)
    @objective(source, Max, 3*x + 1*y)
    @constraint(source, x + 2.0*y <= 6)
    @constraint(source, x*x <= 1)
    addSOS2(source, [x, 2y])
    @SDconstraint(source, x*ones(3,3) >= eye(3,3))
    @SDconstraint(source, ones(3,3) <= 0)
    @variable(source, z[1:3])
    @variable(source, w[2:4]) # JuMPArray
    @variable(source, v[[:red],i=1:3;isodd(i)]) # JuMPDict
    JuMP.setsolvehook(source, m -> :Optimal)

    # uncomment when NLP copying is implemented
    # @NLconstraint(source, c[k=1:3], x^2 + y^3 * sin(x+k*y) >= 1)

    dest = copy(source)

    # uncomment when NLP copying is implemented
    # for name in source.nlpdata
    #     @fact source.name == dest.name --> true
    # end

    # Obj
    @objective(source, Max, 1x)
    @fact length(source.obj.aff.coeffs) --> 1
    @fact length(dest.obj.aff.coeffs) --> 2
    @objective(dest, Max, 1x)
    @fact length(source.obj.aff.coeffs) --> 1
    @fact length(dest.obj.aff.coeffs) --> 1
    @objective(dest, Max, 3*x + 1*y)
    @fact length(source.obj.aff.coeffs) --> 1
    @fact length(dest.obj.aff.coeffs) --> 2

    # Constraints
    source.linconstr[1].ub = 5.0
    @fact dest.linconstr[1].ub --> 6.0
    source.quadconstr[1].sense = :(>=)
    @fact dest.quadconstr[1].sense --> :(<=)
    source.sosconstr[1].terms[1] = Variable(source, 2)
    source.sosconstr[1].terms[2] = Variable(source, 1)
    source.sosconstr[1].weights = [2.0,1.0]
    source.sosconstr[1].sostype = :SOS1
    @fact dest.sosconstr[1].terms[1].col --> 1
    @fact dest.sosconstr[1].terms[2].col --> 2
    @fact dest.sosconstr[1].weights --> [1.0,2.0]
    @fact dest.sosconstr[1].sostype --> :SOS2
    @fact length(dest.sdpconstr) --> 2
    xx = copy(x, dest)
    @fact all(t -> isequal(t[1],t[2]), zip(dest.sdpconstr[1].terms, xx*ones(3,3) - eye(3,3))) --> true
    @fact all(t -> isequal(t[1],t[2]), zip(dest.sdpconstr[2].terms, convert(Matrix{AffExpr}, -ones(3,3)))) --> true

    @fact dest.solvehook(dest) --> :Optimal

    @fact Set(collect(keys(dest.varDict))) --> Set([:x,:y,:z,:w,:v])
    @fact isequal(dest.varDict[:x], Variable(dest, 1)) --> true
    @fact isequal(dest.varDict[:y], Variable(dest, 2)) --> true
    @fact all(t -> isequal(t[1], t[2]), zip(dest.varDict[:z], [Variable(dest, 3), Variable(dest, 4), Variable(dest, 5)])) --> true
    @fact all(t -> isequal(t[1], t[2]), zip(dest.varDict[:w].innerArray, [Variable(dest, 6), Variable(dest, 7), Variable(dest, 8)])) --> true
    td = dest.varDict[:v].tupledict
    @fact length(td) --> 2
    @fact isequal(td[:red,1], Variable(dest, 9))  --> true
    @fact isequal(td[:red,3], Variable(dest, 10)) --> true

    # Issue #358
    @fact typeof(dest.linconstr)  --> Array{JuMP.GenericRangeConstraint{JuMP.GenericAffExpr{Float64,JuMP.Variable}},1}
    @fact typeof(dest.quadconstr) --> Array{JuMP.GenericQuadConstraint{JuMP.GenericQuadExpr{Float64,JuMP.Variable}},1}

    JuMP.setprinthook(source, m -> 2)
    dest2 = copy(source)
    @fact dest2.printhook(dest2) --> 2

    # Test copying model with multiple variables of same name
    @variable(source, x)
    dest3 = copy(source)
    @fact_throws getvariable(dest3, :x)
    @fact dest3.varDict[:x] --> nothing

    addlazycallback(source, cb -> 3)
    @fact_throws copy(source)
end

type TemporaryExtensionTestType
    x::Int
end
facts("[model] Test extension copy") do
    source = Model()
    source.ext[:extensiontype] = TemporaryExtensionTestType(1)
    @fact_throws copy(source)
    source.ext[:extensiontype] = 1
    dest = copy(source)
    source.ext[:extensiontype] = 2
    @fact haskey(dest.ext, :extensiontype) --> true
    @fact dest.ext[:extensiontype] --> 1
end

facts("[model] Test variable/model 'hygiene'") do
context("Linear constraint") do
    modA = Model(); @variable(modA, x)
    modB = Model(); @variable(modB, y)
    @constraint(modA, x+y == 1)
    @fact_throws solve(modA)
end
context("Linear objective") do
    modA = Model(); @variable(modA, 0 <= x <= 1)
    modB = Model(); @variable(modB, 0 <= y <= 1)
    @objective(modA, Min, x + y)
    @fact_throws solve(modA)
end
context("Quadratic constraint") do
    modA = Model(); @variable(modA, x)
    modB = Model(); @variable(modB, y)
    @constraint(modB, x*y >= 1)
    @fact_throws solve(modB)
end
context("Affine in quadratic constraint") do
    modA = Model(); @variable(modA, x)
    modB = Model(); @variable(modB, y)
    @constraint(modB, y*y + x + y <= 1)
    @fact_throws solve(modB)
end
context("Quadratic objective") do
    modA = Model(); @variable(modA, x)
    modB = Model(); @variable(modB, y)
    @objective(modB, Min, x*y)
    @fact_throws solve(modB)
end
end

facts("[model] Test NaN checking") do
    mod = Model()
    @variable(mod, x)
    @objective(mod, Min, NaN*x)
    @fact_throws solve(mod)
    @objective(mod, Min, NaN*x^2)
    @fact_throws solve(mod)
    @objective(mod, Min, x)
    @constraint(mod, Min, NaN*x == 0)
    @fact_throws solve(mod)
end

facts("[model] Test column-wise modeling") do
    mod = Model()
    @variable(mod, 0 <= x <= 1)
    @variable(mod, 0 <= y <= 1)
    @objective(mod, Max, 5x + 1y)
    @constraint(mod, con[i=1:2], i*x + y <= i+5)
    @variable(mod, 0 <= z1 <= 1, objective=10.0, inconstraints=con, coefficients=[1.0,-2.0])
    @variable(mod, 0 <= z2 <= 1, objective=10.0, inconstraints=Any[con[i] for i in 1:2], coefficients=[1.0,-2.0])
    @fact solve(mod) --> :Optimal
    @fact getvalue(z1) --> roughly(1.0, TOL)
    @fact getvalue(z2) --> roughly(1.0, TOL)

    # do a vectorized version as well
    mod = Model()
    @variable(mod, 0 <= x <= 1)
    @variable(mod, 0 <= y <= 1)
    obj = [5,1]'*[x,y]
    @objective(mod, Max, obj[1])
    A = [1 1
         2 1]
    @constraint(mod, A*[x,y] .<= [6,7])
    @variable(mod, 0 <= z1 <= 1, objective=10.0, inconstraints=con, coefficients=[1.0,-2.0])
    @variable(mod, 0 <= z2 <= 1, objective=10.0, inconstraints=Any[con[i] for i in 1:2], coefficients=[1.0,-2.0])
    @fact solve(mod) --> :Optimal
    @fact getvalue(z1) --> roughly(1.0, TOL)
    @fact getvalue(z2) --> roughly(1.0, TOL)

    # vectorized with sparse matrices
    mod = Model()
    @variable(mod, 0 <= x <= 1)
    @variable(mod, 0 <= y <= 1)
    # TODO: get this to work by adding method for Ac_mul_B!
    # obj = sparse([5,1])'*[x,y]
    obj = sparse([5,1]')*[x,y]
    @objective(mod, Max, obj[1])
    A = sparse([1 1
                2 1])
    @constraint(mod, A*[x,y] .<= [6,7])
    @variable(mod, 0 <= z1 <= 1, objective=10.0, inconstraints=con, coefficients=[1.0,-2.0])
    @variable(mod, 0 <= z2 <= 1, objective=10.0, inconstraints=Any[con[i] for i in 1:2], coefficients=[1.0,-2.0])
    @fact solve(mod) --> :Optimal
    @fact getvalue(z1) --> roughly(1.0, TOL)
    @fact getvalue(z2) --> roughly(1.0, TOL)
end

facts("[model] Test all MPS paths") do
    mod = Model()
    @variable(mod, free_var)
    @variable(mod, int_var, Bin)
    @variable(mod, low_var >= 5)
    @constraint(mod, free_var == int_var)
    @constraint(mod, free_var - int_var >= 0)
    @objective(mod, Max, free_var*int_var + low_var)
    writeMPS(mod,"test.mps")
end

facts("[model] Test all LP paths") do
    mod = Model()
    @variable(mod, free_var)
    @objective(mod, Max, free_var*free_var)
    @fact_throws writeLP(mod,"test.lp")
    @objective(mod, Max, free_var)
    @constraint(mod, free_var - 2*free_var == 0)
    @constraint(mod, free_var + 2*free_var >= 1)
    writeLP(mod,"test.lp")
end

facts("[model] Test semi-continuous variables") do
for solver in semi_solvers
context("With solver $(typeof(solver))") do
    mod = Model(solver=solver)
    @variable(mod, x >= 3, SemiCont)
    @variable(mod, y >= 2, SemiCont)
    @constraint(mod, x + y >= 1)
    @objective(mod, Min, x+y)
    solve(mod)
    @fact getvalue(x) --> roughly(0.0, TOL)
    @fact getvalue(y) --> roughly(2.0, TOL)
end; end; end

facts("[model] Test semi-integer variables") do
for solver in semi_solvers
context("With solver $(typeof(solver))") do
    mod = Model(solver=solver)
    @variable(mod, x >= 3, SemiInt)
    @variable(mod, y >= 2, SemiInt)
    @constraint(mod, x + y >= 2.5)
    @objective(mod, Min, x+1.1y)
    solve(mod)
    @fact getvalue(x) --> roughly(3.0, TOL)
    @fact getvalue(y) --> 0.0
end; end; end

facts("[model] Test fixed variables don't leak through MPB") do
for solver in lp_solvers
context("With solver $(typeof(solver))") do
    mod = Model(solver=solver)
    @variable(mod, 0 <= x[1:3] <= 2)
    @variable(mod, y[k=1:2] == k)
    @objective(mod, Min, x[1] + x[2] + x[3] + y[1] + y[2])
    solve(mod)
    for i in 1:3
        @fact getvalue(x[i]) --> roughly(0, TOL)
    end
    for k in 1:2
        @fact getvalue(y[k]) --> roughly(k, TOL)
    end
end; end
for solver in ip_solvers
context("With solver $(typeof(solver))") do
    mod = Model(solver=solver)
    @variable(mod, x[1:3], Bin)
    @variable(mod, y[k=1:2] == k)
    JuMP.build(mod)
    @fact MathProgBase.getvartype(internalmodel(mod)) --> [:Bin,:Bin,:Bin,:Cont,:Cont]
end; end; end

facts("[model] Test SOS constraints") do
for solver in sos_solvers
context("With solver $(typeof(solver))") do
    modS = Model(solver=solver)

    @variable(modS, x[1:3], Bin)
    @variable(modS, y[1:5], Bin)
    @variable(modS, z)
    @variable(modS, w)

    @objective(modS, Max, z+w)

    a = [1,2,3]
    b = [5,4,7,2,1]

    @constraint(modS, z == sum{a[i]*x[i], i=1:3})
    @constraint(modS, w == sum{b[i]*y[i], i=1:5})

    @fact_throws constructSOS([x[1]+y[1]])
    @fact_throws constructSOS([1z])

    addSOS1(modS, [a[i]x[i] for i in 1:3])
    addSOS2(modS, [b[i]y[i] for i in 1:5])

    @fact_throws addSOS1(modS, [x[1], x[1]+x[2]])

    @fact solve(modS) --> :Optimal
    @fact modS.objVal --> roughly(15.0, TOL)
    @fact getvalue(z) --> roughly( 3.0, TOL)
    @fact getvalue(w) --> roughly(12.0, TOL)

    m = Model(solver=solver)
    ub = [1,1,2]
    @variable(m, x[i=1:3] <= ub[i])
    @objective(m, Max, 2x[1]+x[2]+x[3])
    addSOS1(m, [x[1],2x[2]])
    addSOS1(m, [x[1],2x[3]])

    @fact solve(m) --> :Optimal
    @fact getvalue(x)[:] --> roughly([0.0,1.0,2.0], TOL)
    @fact getobjectivevalue(m) --> roughly(3.0, TOL)
end; end; end

facts("[model] Test vectorized model creation") do

    A = sprand(50,10,0.15)
    B = sprand(50, 7,0.2)
    modV = Model()
    @variable(modV, x[1:10])
    @variable(modV, y[1:7])
    @constraint(modV, A*x + B*y .<= 1)
    obj = (x'*2A')*(2A*x) + (B*2y)'*(B*(2y))
    @objective(modV, Max, obj[1])

    modS = Model()
    @variable(modS, x[1:10])
    @variable(modS, y[1:7])
    for i in 1:50
        @constraint(modS, sum{A[i,j]*x[j], j=1:10} + sum{B[i,k]*y[k], k=1:7} <= 1)
    end
    AA, BB = 4A'*A, 4B'*A
    @objective(modS, Max, sum{AA[i,j]*x[i]*x[j], i=1:10,j=1:10} + sum{BB[i,j]*y[i]*y[j], i=1:7, j=1:7})

    @fact JuMP.prepConstrMatrix(modV) --> JuMP.prepConstrMatrix(modS)
    @fact JuMP.prepProblemBounds(modV) --> JuMP.prepProblemBounds(modS)
end

facts("[model] Test MIQP vectorization") do
    n = 1000
    p = 4
    function bestsubset(solver,X,y,K,M,integer)
        mod = Model(solver=solver)
        @variable(mod, β[1:p])
        if integer
            @variable(mod, z[1:p], Bin)
        else
            @variable(mod, 0 <= z[1:p] <= 1)
        end
        obj = (y-X*β)'*(y-X*β)
        @objective(mod, Min, obj[1])
        @constraint(mod, eye(p)*β .<=  M*eye(p)*z)
        @constraint(mod, eye(p)*β .>= -M*eye(p)*z)
        @constraint(mod, sum(z) == K)
        solve(mod)
        return getvalue(β)[:]
    end
    include(joinpath("data","miqp_vector.jl")) # loads X and q
    y = X * [100, 50, 10, 1] + 20*q
    for solver in quad_solvers
        @fact bestsubset(solver,X,y,2,500,false) --> roughly([101.789,49.414,8.63904,1.72663], 10TOL)
        @fact bestsubset(solver,sparse(X),y,2,500,false) --> roughly([101.789,49.414,8.63904,1.72663], 10TOL)
    end
    for solver in quad_mip_solvers
        y = X * [100, 50, 10, 1] + 20*q
        @fact bestsubset(solver,X,y,2,500,true) --> roughly([106.25,53.7799,0.0,0.0], 1.0)
        @fact bestsubset(solver,sparse(X),y,2,500,true) --> roughly([106.25,53.7799,0.0,0.0], 1.0)
    end
end

facts("[model] Test setsolver") do
    m = Model()
    @variable(m, x[1:5])
    @constraint(m, con[i=1:5], x[6-i] == i)
    @fact solve(m) --> :Optimal

    for solver in lp_solvers
        setsolver(m, solver)
        @fact m.solver --> solver
        @fact (m.internalModel == nothing) --> true
        @fact solve(m) --> :Optimal
        @fact m.solver --> solver
        @fact (m.internalModel == nothing) --> false
    end
end

facts("[model] Setting solve hook") do
    m = Model()
    @variable(m, x ≥ 0)
    dummy = [1]
    kwarglist = Any[]
    function solvehook(m::Model; kwargs...)
        dummy[1] += 1
        append!(kwarglist, kwargs)
        :DidntDoAnything
    end
    JuMP.setsolvehook(m, solvehook)
    solve(m)
    @fact dummy --> [2]
    @fact kwarglist --> Any[(:suppress_warnings,false)]
end

facts("[model] Setting print hook") do
    m = Model()
    @variable(m, x ≥ 0)
    dummy = [1]
    function printhook(io::IO, m::Model)
        dummy[1] += 1
    end
    JuMP.setprinthook(m, printhook)
    print(m)
    @fact dummy --> [2]
end

facts("[model] Test linearindex") do
    m = Model()
    @variable(m, x[1:5])
    for i in 1:5
        @fact isequal(Variable(m,linearindex(x[i])),x[i]) --> true
    end
end

facts("[model] Test LinearConstraint from ConstraintRef") do
    m = Model()
    @variable(m, x)
    @constraint(m, constr, x == 1)
    real_constr = LinearConstraint(constr)
    @fact real_constr.terms --> @LinearConstraint(x == 1).terms
end

facts("[model] Test getvalue on OneIndexedArrays") do
    m = Model()
    @variable(m, x[i=1:5], start=i)
    @fact typeof(x) --> Vector{Variable}
    @fact typeof(getvalue(x)) --> Vector{Float64}
end

facts("[model] Relaxation keyword argument to solve") do
    m = Model()
    @variable(m, 1.5 <= y <= 2, Int)
    @variable(m, z, Bin)
    @variable(m, 0.5 <= w <= 1.5, Int)
    @variable(m, 1 <= v <= 2)

    @objective(m, Min, y + z + w + v)

    # Force LP solver since not all MIP solvers
    # return duals (i.e. Cbc)
    setsolver(m, lp_solvers[1])
    @fact solve(m, relaxation=true) --> :Optimal
    @fact getvalue(y) --> 1.5
    @fact getvalue(z) --> 0
    @fact getvalue(w) --> 0.5
    @fact getvalue(v) --> 1
    @fact getdual(y) --> 1
    @fact getdual(z) --> 1
    @fact getdual(w) --> 1
    @fact getdual(v) --> 1
    @fact getobjectivevalue(m) --> 1.5 + 0 + 0.5 + 1

    # Let JuMP choose solver again
    setsolver(m, JuMP.UnsetSolver())
    @fact solve(m) --> :Optimal
    @fact getvalue(y) --> 2
    @fact getvalue(z) --> 0
    @fact getvalue(w) --> 1
    @fact getvalue(v) --> 1
    @fact getobjectivevalue(m) --> 2 + 0 + 1 + 1

    @variable(m, 1 <= x <= 2, SemiCont)
    @variable(m, -2 <= t <= -1, SemiInt)

    addSOS1(m, [x, 2y, 3z, 4w, 5v, 6t])
    @objective(m, Min, x + y + z + w + v - t)

    @fact solve(m, relaxation=true) --> :Optimal

    @fact getvalue(x) --> 0
    @fact getvalue(y) --> 1.5
    @fact getvalue(z) --> 0
    @fact getvalue(w) --> 0.5
    @fact getvalue(v) --> 1
    @fact getvalue(t) --> 0
    @fact getobjectivevalue(m) --> 0 + 1.5 + 0 + 0.5 + 1 + 0
end

facts("[model] Unrecognized keyword argument to solve") do
    m = Model()
    @fact_throws solve(m, this_should_throw=true)
end

facts("[model] Solve MIP relaxation with continuous solvers") do
for solver in lp_solvers
context("With solver $(typeof(solver))") do
    m = Model(solver=solver)
    @variable(m, x, Bin)
    @constraint(m, x >= 0.5)
    @objective(m, Min, x)
    @fact solve(m, relaxation=true) --> :Optimal
    @fact getvalue(x) --> roughly(0.5, TOL)
end; end
for solver in nlp_solvers
context("With solver $(typeof(solver))") do
    m = Model(solver=solver)
    @variable(m, x, Bin)
    @objective(m, Min, x)
    @NLconstraint(m, x >= 0.5)
    @fact solve(m, relaxation=true) --> :Optimal
    @fact getvalue(x) --> roughly(0.5, TOL)
end; end
end

facts("[model] Nonliteral exponents in @constraint") do
    m = Model()
    @variable(m, x)
    foo() = 2
    @constraint(m, x^(foo()) + x^(foo()-1) + x^(foo()-2) == 1)
    @constraint(m, (x-1)^(foo()) + (x-1)^2 + (x-1)^1 + (x-1)^0 == 1)
    @constraint(m, sum{x, i in 1:3}^(foo()) == 1)
    @constraint(m, sum{x, i in 1:3}^(foo()-1) == 1)
    @fact m.quadconstr[1].terms --> x^2 + x
    @fact m.quadconstr[2].terms --> x^2 + x^2 - x - x - x - x + x + 1
    @fact m.quadconstr[3].terms --> x^2 + x^2 + x^2 + x^2 + x^2 + x^2 + x^2 + x^2 + x^2 - 1
    @fact m.quadconstr[4].terms --> QuadExpr(x + x + x - 1)
end


text:
x_shift = deepcopy(x_sol)
for x in x_shift
    x[3] += 0.01
end
z = [minimal_to_maximal(env.mechanism, x) for x in x_shift]

t = 1 #10, 20, 30, 41
set_robot(env.vis, env.mechanism, z[41])

# ## visualize
x_view = [[x_shift[1] for t = 1:15]..., x_shift..., [x_shift[end] for t = 1:15]...]
visualize(env, x_view)

set_camera!(vis, cam_pos=[0,-50,0], zoom=30)
set_floor!(env.vis, z=0.01)

ext: jl


text:
using LinearAlgebra, SolverCore
mutable struct problem_summary # for backward compatibility
    status::Symbol
    it_count::Int64
    total_time::Float64
    fval::Float64
    con_vio::Float64
    dual_feas::Float64
    comp::Float64

    function problem_summary()
        return new(:Blank,-2,NaN,NaN,NaN)
    end
end
mutable struct problem_summary2
    status::Symbol
    it_count::Int64
    total_time::Float64
    fval::Float64
    con_vio::Float64
    dual_feas::Float64
    comp::Float64
    dual_max::Float64

    #size
    number_variables::Int64
    number_constraints::Int64
    #toal evaluations
    total_fval_evaluation::Int64
    total_grad_evaluation::Int64
    total_jac_evaluation::Int64
    total_cons_evaluation::Int64
    total_hess_evaluation::Int64

    function problem_summary2()
        return new(:Blank,-2,NaN,NaN,NaN,NaN,NaN,0,0,0,0,0,0,0)
    end
end

function cps(ps::problem_summary)
    ps_new = problem_summary2()
    ps_new.status = ps.status
    ps_new.it_count = ps.it_count
    ps_new.total_time = ps.total_time
    ps_new.fval = ps.fval
    ps_new.con_vio = ps.con_vio
    ps_new.dual_feas = ps.dual_feas
    ps_new.comp = ps.comp
    ps_new.dual_max = NaN

    ps_new.number_variables = ps.number_variables
    ps_new.number_variables = ps.number_variables
    #toal evaluations
    ps_new.total_fval_evaluation = ps.total_fval_evaluation
    ps_new.total_grad_evaluation = ps.total_grad_evaluation
    ps_new.total_jac_evaluation = ps.total_jac_evaluation
    ps_new.total_cons_evaluation = ps.total_cons_evaluation
    ps_new.total_hess_evaluation = ps.total_hess_evaluation

    return ps_new
end

function cps(ls::Dict{String,problem_summary})
    new_ls = Dict{String,problem_summary2}()
    for (problem_name,info) in ls
      new_ls[problem_name] = cps(info)
    end

    return new_ls
end

function cps(ls::Dict{String,problem_summary2})
    return ls
end

#= 31st May
type problem_summary2
    status::Symbol
    it_count::Int64
    total_time::Float64

    function problem_summary2()
        return new()
    end
end
=#

function write_summary(stream::IOStream, summary::Dict{String, problem_summary2})
    write(stream, "problem_name \t status \t it count \t total_time \n");
    for (problem_name, info) in summary
      # write(stream, "$(OnePhase.pd(problem_name)) \t $(OnePhase.pd(info.status)) \t $(OnePhase.pd(info.it_count)) \t $(OnePhase.rd(info.total_time))\n");
       write(stream, "$(OnePhase.pd(problem_name)) \t $(OnePhase.pd(info.status)) \t $(OnePhase.pd(info.it_count)) \t $(OnePhase.rd(info.total_time))\n");
    end
end

function write_summary(file_name::String, summary::Dict{String, problem_summary2})
    summary_file = open(file_name, "w")
    write_summary(summary_file, summary)
    close(summary_file)
end

function set_info_me!(info::problem_summary2, hist, status::Symbol, iter)
    set_info_me!(info, hist, status)
    info.total_fval_evaluation = iter.nlp.nlp.counters.neval_obj
    info.total_grad_evaluation = iter.nlp.nlp.counters.neval_grad
    info.total_jac_evaluation = iter.nlp.nlp.counters.neval_jac
    info.total_cons_evaluation =  iter.nlp.nlp.counters.neval_cons
    info.total_hess_evaluation = iter.nlp.nlp.counters.neval_hess
end

function set_info_me!(info::problem_summary2, hist, status::Symbol)
    info.it_count = hist[end].t
    info.fval = hist[end].fval;
    info.con_vio = hist[end].con_vio;
    info.dual_feas = hist[end].norm_grad_lag;
    info.comp = hist[end].comp;
    info.dual_max = hist[end].y_norm
    info.status = status
end

function set_cutest_info_ipopt!(info::problem_summary2, stats::SolverCore.GenericExecutionStats{Float64, Vector{Float64}}, nlp_raw::AbstractNLPModel, x::Array{Float64,1})
  num_vars = length(nlp_raw.meta.lvar)
  x_true = x[1:num_vars]

  info.fval = stats.objective
  info.total_fval_evaluation = stats.counters.neval_obj
  info.total_grad_evaluation = stats.counters.neval_grad
  info.total_jac_evaluation = stats.counters.neval_jac
  info.total_cons_evaluation =  stats.counters.neval_cons
  info.total_hess_evaluation = stats.counters.neval_hess

  a = cons(nlp_raw, x_true);
  info.con_vio = max(0.0, maximum(nlp_raw.meta.lvar - x_true), maximum(x_true - nlp_raw.meta.uvar), maximum(nlp_raw.meta.lcon - a), maximum(a - nlp_raw.meta.ucon))

  info.dual_feas = norm(grad(nlp_raw, x_true) + jac(nlp_raw, x_true)' * stats.multipliers + stats.multipliers_U - stats.multipliers_L, Inf);
  info.comp = maximum(stats.multipliers_U .* min.(1e16, abs.(nlp_raw.meta.uvar - x_true) ) ) + maximum( stats.multipliers_L .* min.(1e16, abs.(x_true - nlp_raw.meta.lvar)) )
  info.dual_max = max(maximum(stats.multipliers_U),maximum(stats.multipliers_L), maximum(abs.(stats.multipliers)))
end

max_stars_repo_name: ohinder/OnePhase

max_issues_repo_name: ohinder/OnePhase

max_forks_repo_name: ohinder/OnePhase


text:
include("RbfModel.jl")
include("ExactModel.jl")
include("TaylorModel.jl")
include("LagrangeModel.jl")

struct SurrogateWrapper{
        O <: AbstractObjective,
        M <: SurrogateModel,
        I <: SurrogateMeta} <: AbstractSurrogateWrapper
    objf :: O
    model :: M
    meta :: I
end

function init_wrapper( :: Type{<:SurrogateWrapper}, objf :: AbstractObjective,
    model :: SurrogateModel, meta :: SurrogateMeta)
    return SurrogateWrapper( objf, model, meta )
end

num_outputs( sw :: SurrogateWrapper ) :: Int = num_outputs(get_objf(sw))
get_objf( sw :: SurrogateWrapper ) = sw.objf
get_model( sw :: SurrogateWrapper ) = sw.model
get_meta( sw :: SurrogateWrapper ) = sw.meta

struct SurrogateContainer{
        T <: AbstractVector{<:AbstractSurrogateWrapper} } <: AbstractSurrogateContainer
    surrogates :: T

    "A dict of tuples where the key is a MOP output index and the
    value is a tuple where the first entry is the index of the
    corresponding surrogate in `surrogates` and the second entry is the
    index of the surrogate model output (or 1 if it's a scalar models)."
    output_model_mapping :: Base.ImmutableDict{Int,Tuple{Int,Int}}

    "A dict of tuples where the key is an output index of the
    SurrogateContainer (i.e., corresponding to the internal sorting) and the
    value is a tuple where the first entry is the index of the
    corresponding surrogate in `surrogates` and the second entry is the
    index of the surrogate model output (or 1 if it's a scalar models)."
    sc_output_model_mapping :: Base.ImmutableDict{Int,Tuple{Int,Int}}
end

list_of_wrappers( sc :: SurrogateContainer ) = sc.surrogates
get_wrapper( sc, i ) = sc.surrogates[i]

_output_model_mapping( sc :: SurrogateContainer, l ) = sc.output_model_mapping[l]
_sc_output_model_mapping( sc :: SurrogateContainer, l ) = sc.sc_output_model_mapping[l]

@doc "Return a SurrogateContainer initialized from the information provided in `mop`."
function init_surrogates( ::Type{<:SurrogateContainer}, mop :: AbstractMOP,
    id :: AbstractIterData, db :: AbstractDB, ac :: AbstractConfig )

    @logmsg loglevel2 "Initializing surrogate models."

    # init fields for SurrogateContainer
    sw_array = SurrogateWrapper[]   # TODO make AbstractSurrogateWrapper type configurable?
    output_model_mapping = Dict{Int,Tuple{Int,Int}}()
    sc_output_model_mapping = Dict{Int,Tuple{Int,Int}}()

    meta_array = []

    for objf ∈ list_of_objectives(mop)
        push!( meta_array, prepare_init_model( model_cfg( objf ), objf, mop, id, db, ac; meta_array ))
    end

    @logmsg loglevel2 "Evaluation of unevaluated results."
    eval_missing!(db, mop)

    L = 1   # counter of scalarized outputs of `mop`
    for (i, objf) ∈ enumerate(list_of_objectives(mop))
        meta_dat = meta_array[i]
        model, meta = _init_model( model_cfg(objf), objf, mop, id, db, ac, meta_dat )
        push!( sw_array, init_wrapper(SurrogateWrapper, objf, model, meta) )

        for (j, ℓ) in enumerate( output_indices( objf, mop ) )
            output_model_mapping[ℓ] = (i,j)
            sc_output_model_mapping[L] = (i,j)
            L += 1
        end
    end

    return SurrogateContainer(
        sw_array,
        Base.ImmutableDict(output_model_mapping...),
        Base.ImmutableDict(sc_output_model_mapping...)
    )
end

function update_surrogates!( sc :: SurrogateContainer, mop :: AbstractMOP,
    id :: AbstractIterData, db :: AbstractDB, ac :: AbstractConfig; ensure_fully_linear :: Bool = false ) :: Nothing
    @logmsg loglevel2 "Updating surrogate models."

    meta_array = SurrogateMeta[]
    for sw ∈ sc.surrogates
        push!(meta_array, prepare_update_model(get_model(sw), get_objf(sw), get_meta(sw), mop, id, db, ac; ensure_fully_linear, meta_array ) )
    end
    eval_missing!(db, mop)

    for (si,sw) ∈ enumerate(sc.surrogates)
        meta_dat = meta_array[si]
        new_model, new_meta = update_model( get_model(sw), get_objf(sw), meta_dat, mop, id, db, ac; ensure_fully_linear )
        new_sw = init_wrapper(
            SurrogateWrapper,
            get_objf(sw),
            new_model,
            new_meta
        )
        sc.surrogates[si] = new_sw
    end
    nothing
end

function improve_surrogates!( sc :: SurrogateContainer, mop :: AbstractMOP,
    id :: AbstractIterData, db :: AbstractDB, ac :: AbstractConfig;
    ensure_fully_linear :: Bool = false, kwargs... ) :: Nothing
    @logmsg loglevel2 "Improving surrogate models."

    meta_array = SurrogateMeta[]
    for sw ∈ sc.surrogates
        push!(meta_array, prepare_improve_model(get_model(sw), get_objf(sw), get_meta(sw), mop, id, db, ac; ensure_fully_linear ) )
    end
    eval_missing!(db, mop)

    for (si,sw) ∈ enumerate(sc.surrogates)
        meta_dat = meta_array[si]
        new_model, new_meta = improve_model( get_model(sw), get_objf(sw), meta_dat, mop, id, db, ac; ensure_fully_linear, kwargs... )
        new_sw = SurrogateWrapper(
            get_objf(sw),
            new_model,
            new_meta
        )
        sc.surrogates[si] = new_sw
    end
    nothing
end


text:
using IndexedTables, SparseArrays, BenchmarkTools, Random

N = 10_000

Random.seed!(666)
s = table((a = rand(1:3, N), b = rand(1:3, N), c1 = rand(N), d1 = rand(N)), pkey = (:a, :b))
t = table((a = rand(1:3, N), b = rand(1:3, N), c2 = rand(N), d2 = rand(N)), pkey = (:a, :b))

join(s, t);
join(merge, s, t);
groupjoin(s, t);
@time join(s, t);
@time join(merge, s, t);
@time groupjoin(s, t);

Random.seed!(666)
idx = Columns(p=rand(1:100, N), q=rand(1:100, N))
t = NDSparse(idx, rand(N))
t2 = NDSparse(Columns(q=rand(1:100, N)), rand(N))
@btime broadcast(+, $t, $t2)

using SparseArrays
Random.seed!(666)
S = sprand(1000, 1000,.1)
v = rand(1000)
nd = convert(NDSparse, S)
ndv = convert(NDSparse,v)
@btime broadcast(*, $nd, $ndv)


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

name = "FastJet_Julia_Wrapper"
version = v"0.7.0"

# Collection of sources required to build FastJet_Julia_Wrapper
sources = [
	GitSource("https://github.com/jstrube/FastJet_Julia_Wrapper.git", "6d4181ca351c7a40348745097c87afd97bf9ce62"; unpack_target="FastJet_Julia_Wrapper"),
    ArchiveSource("https://julialang-s3.julialang.org/bin/linux/x64/1.3/julia-1.3.1-linux-x86_64.tar.gz", "faa707c8343780a6fe5eaf13490355e8190acf8e2c189b9e7ecbddb0fa2643ad"; unpack_target="julia-x86_64-linux-gnu"),
    ArchiveSource("https://github.com/Gnimuc/JuliaBuilder/releases/download/v1.3.0/julia-1.3.0-x86_64-apple-darwin14.tar.gz", "f2e5359f03314656c06e2a0a28a497f62e78f027dbe7f5155a5710b4914439b1"; unpack_target="julia-x86_64-apple-darwin14"),
]

# Bash recipe for building across all platforms
script = raw"""
case "$target" in
	arm-linux-gnueabihf|x86_64-linux-gnu)
        Julia_PREFIX=${WORKSPACE}/srcdir/julia-$target/julia-1.3.1
        ;;
    x86_64-apple-darwin14|x86_64-w64-mingw32)
        Julia_PREFIX=${WORKSPACE}/srcdir/julia-$target/juliabin
        ;;
esac

cd ${WORKSPACE}/srcdir/FastJet_Julia_Wrapper/FastJet_Julia_Wrapper*
mkdir build && cd build
cmake -DJulia_PREFIX=${Julia_PREFIX} -DCMAKE_INSTALL_PREFIX=${prefix} -DCMAKE_FIND_ROOT_PATH=${prefix} -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DCMAKE_BUILD_TYPE=Release ..
VERBOSE=ON cmake --build . --config Release --target install
install_license $WORKSPACE/srcdir/FastJet_Julia_Wrapper/FastJet_Julia_Wrapper/LICENSE.md
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = Platform[
    Linux(:x86_64, libc=:glibc),
    MacOS(:x86_64)
]
platforms = expand_cxxstring_abis(platforms)

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

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("libcxxwrap_julia_jll"),
    Dependency("FastJet_jll")
]

build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version=v"7")

max_stars_repo_name: tbeason/Yggdrasil

max_issues_repo_name: tbeason/Yggdrasil

max_forks_repo_name: tbeason/Yggdrasil


text:
### A Pluto.jl notebook ###
# v0.12.6

using Markdown
using InteractiveUtils

# ╔═╡ acf67f82-1a85-11eb-2bdf-abc05bfb68eb
begin
	using DataFrames
	using Statistics
	using CSV
	using FreqTables
	using Pipe
end

# ╔═╡ 8061b6ba-1a85-11eb-1a73-130f152eb8d8
VERSION

# ╔═╡ ad304f82-1a85-11eb-1b1e-bba3c3d49c45
isfile.(["Project.toml", "Manifest.toml"])

# ╔═╡ 4356bc2c-1b0c-11eb-069b-155c8cf6512e
ENV["LINES"], ENV["COLUMNS"] = 15, 200

# ╔═╡ 8849a66c-1adc-11eb-3c30-371f47ba053e
md"""
# Working with groups of rows
"""

# ╔═╡ 23c9e8ac-1b0c-11eb-0834-1bdf38850a22
df = CSV.File("auto2.csv") |> DataFrame

# ╔═╡ 55f5d282-1b0c-11eb-2664-f9739e646807
md"""
We want to group our data frame by :brand
"""

# ╔═╡ 23aeac54-1b0c-11eb-3dbb-4566edc7b0e8
gdf = groupby(df, :brand)

# ╔═╡ 238fbd4c-1b0c-11eb-0d12-2b9942b58740
gdf[("ford",)]

# ╔═╡ 2334882c-1b0c-11eb-39d6-51ac9fa55c4f
md"""
Using the combine function we can easily calculate some aggregates by group:
"""

# ╔═╡ fce8af84-1b0b-11eb-0c3f-f17c2350aa0d
brand_mpg = combine(gdf, :mpg => mean)

# ╔═╡ fccdd632-1b0b-11eb-0ca3-19c005a0bc47
combine(gdf, :mpg => mean => :mean_mpg)

# ╔═╡ fcb3d11a-1b0b-11eb-0e42-8fdc9591d38e
sort!(brand_mpg, :mpg_mean, rev=true)

# ╔═╡ fc976412-1b0b-11eb-3718-39db026e228e
md"""
A typical data cleaning task is to check its consistency. In this case it would be making sure that each brand has a unique origin. We will go it several ways.
"""

# ╔═╡ fc7d9802-1b0b-11eb-20ce-15226fd3a6f5
freqtable(df, :brand, :origin)

# ╔═╡ 1d627fae-1b0e-11eb-36ce-a1d67dd1e3ad
orig_brand = @pipe df |>
				   groupby(_, :brand) |>
                   combine(_, :origin => x -> length(unique(x)))

# ╔═╡ 1d7e8a3c-1b0e-11eb-3dfc-7fdf2376b161
extrema(orig_brand.origin_function)

# ╔═╡ 1d9a9f74-1b0e-11eb-3c64-e79a33c15318
origin_brand2 = @pipe df |>
                      groupby(_, [:origin, :brand]) |>
                      combine(_, nrow)

# ╔═╡ 1db60408-1b0e-11eb-0098-7dab96566a24
origin_vs_brand = unstack(origin_brand2, :brand, :origin, :nrow)

# ╔═╡ 1dd04eb2-1b0e-11eb-1636-3ba9f2ee7752
coalesce.(origin_vs_brand, 0)

# ╔═╡ 1debc2bc-1b0e-11eb-3be3-d981c8e2e6a5
names(origin_vs_brand)

# ╔═╡ 1e05f80a-1b0e-11eb-10b7-cfb0c4e0cf8b
propertynames(origin_vs_brand)

# ╔═╡ fc62b794-1b0b-11eb-0a82-f3cadd0ffb1c
origin_vs_brand3 = @pipe df |>
                         groupby(_, :origin) |>
                         combine(_, :brand => x -> Ref(unique(x)))

# ╔═╡ fc450492-1b0b-11eb-1257-1f35afd6ba2f
@pipe df |> groupby(_, :origin) |> combine(_, :brand => unique)

# ╔═╡ fc26545c-1b0b-11eb-025a-93aca39ca5fb
flatten(origin_vs_brand3, :brand_function)

# ╔═╡ fc091054-1b0b-11eb-00f2-4b26f065152b

# ╔═╡ fbea2a7c-1b0b-11eb-3dfd-d1edcb81816f

# ╔═╡ c5110f36-1aff-11eb-2712-e56fcfd44535

# ╔═╡ c4f62900-1aff-11eb-319c-6fbed55179ff

# ╔═╡ a79492ea-1ac8-11eb-098c-f12b9d38f313

# ╔═╡ a76f13ea-1ac8-11eb-3800-b1c27e74a466

# ╔═╡ a76f037a-1ac8-11eb-199c-cb332db6cca1

# ╔═╡ a7282752-1ac8-11eb-0f05-83cace66d3b0

# ╔═╡ ab7c9592-1a85-11eb-342f-3dab4c1d54a6

# ╔═╡ Cell order:
# ╠═8061b6ba-1a85-11eb-1a73-130f152eb8d8
# ╠═ad304f82-1a85-11eb-1b1e-bba3c3d49c45
# ╠═acf67f82-1a85-11eb-2bdf-abc05bfb68eb
# ╠═4356bc2c-1b0c-11eb-069b-155c8cf6512e
# ╟─8849a66c-1adc-11eb-3c30-371f47ba053e
# ╠═23c9e8ac-1b0c-11eb-0834-1bdf38850a22
# ╟─55f5d282-1b0c-11eb-2664-f9739e646807
# ╠═23aeac54-1b0c-11eb-3dbb-4566edc7b0e8
# ╠═238fbd4c-1b0c-11eb-0d12-2b9942b58740
# ╟─2334882c-1b0c-11eb-39d6-51ac9fa55c4f
# ╠═fce8af84-1b0b-11eb-0c3f-f17c2350aa0d
# ╠═fccdd632-1b0b-11eb-0ca3-19c005a0bc47
# ╠═fcb3d11a-1b0b-11eb-0e42-8fdc9591d38e
# ╟─fc976412-1b0b-11eb-3718-39db026e228e
# ╠═fc7d9802-1b0b-11eb-20ce-15226fd3a6f5
# ╠═1d627fae-1b0e-11eb-36ce-a1d67dd1e3ad
# ╠═1d7e8a3c-1b0e-11eb-3dfc-7fdf2376b161
# ╠═1d9a9f74-1b0e-11eb-3c64-e79a33c15318
# ╠═1db60408-1b0e-11eb-0098-7dab96566a24
# ╠═1dd04eb2-1b0e-11eb-1636-3ba9f2ee7752
# ╠═1debc2bc-1b0e-11eb-3be3-d981c8e2e6a5
# ╠═1e05f80a-1b0e-11eb-10b7-cfb0c4e0cf8b
# ╠═fc62b794-1b0b-11eb-0a82-f3cadd0ffb1c
# ╠═fc450492-1b0b-11eb-1257-1f35afd6ba2f
# ╠═fc26545c-1b0b-11eb-025a-93aca39ca5fb
# ╠═fc091054-1b0b-11eb-00f2-4b26f065152b
# ╠═fbea2a7c-1b0b-11eb-3dfd-d1edcb81816f
# ╠═c5110f36-1aff-11eb-2712-e56fcfd44535
# ╠═c4f62900-1aff-11eb-319c-6fbed55179ff
# ╠═a79492ea-1ac8-11eb-098c-f12b9d38f313
# ╠═a76f13ea-1ac8-11eb-3800-b1c27e74a466
# ╠═a76f037a-1ac8-11eb-199c-cb332db6cca1
# ╠═a7282752-1ac8-11eb-0f05-83cace66d3b0
# ╠═ab7c9592-1a85-11eb-342f-3dab4c1d54a6

max_stars_repo_name: martin-fabbri/colab-notebooks

max_issues_repo_name: martin-fabbri/colab-notebooks

max_forks_repo_name: martin-fabbri/colab-notebooks


text:
using GnuTLS
using Base.Test

h = HMACState(MD5,"")
update(h,"")
@test hexdigest!(h) == "74e6f7298a9c2d168935f58c001bad88"
h = HMACState(SHA1,"")
update(h,"")
@test hexdigest!(h) == "fbdb1d1b18aa6c08324b7d64b71fb76370690e1d"
h = HMACState(SHA256,"")
update(h,"")
@test hexdigest!(h) == "b613679a0814d9ec772f95d778c35fc5ff1697c493715653c6c712144292c5ad"
h = HMACState(MD5,"key")
update(h,"The quick brown fox jumps over the lazy dog")
@test hexdigest!(h) == "80070713463e7749b90c2dc24911e275"
h = HMACState(SHA1,"key")
update(h,"The quick brown fox jumps over the lazy dog")
@test hexdigest!(h) == "de7c9b85b8b78aa6bc8a7a36f70a90701c9db4d9"
h = HMACState(SHA256,"key")
update(h,"The quick brown fox jumps over the lazy dog")
@test hexdigest!(h) == "f7bc83f430538424b13298e6aa6fb143ef4d59a14946175997479dbc2d1a3cd8"


text:
using ..Atmos
using ..ConfigTypes
using ..DiagnosticsMachine

@diagnostics_group(
    "AtmosGCMDefault",          # name
    AtmosGCMConfigType,         # configuration type
    Nothing,                    # params type
    (_...) -> nothing,          # initialization function
    InterpolateAfterCollection, # if/when to interpolate
    # various pointwise variables
    u,
    v,
    w,
    rho,
    temp,
    pres,
    thd,
    et,
    ei,
    ht,
    hi,
    #vort, TODO
    # moisture related
    qt,
    ql,
    qv,
    qi,
    thv,
    thl,
)


text:
module KNITRO

const _DEPS_FILE = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
if isfile(_DEPS_FILE)
    include(_DEPS_FILE)
else
    error("KNITRO.jl not properly installed. Please run `] build KNITRO`")
end

using Libdl, SparseArrays
import Base: show

const IS_KNITRO_LOADED = endswith(libknitro, Libdl.dlext)

const KNITRO_VERSION = if IS_KNITRO_LOADED
    len = 15
    out = zeros(Cchar, len)
    ccall((:KTR_get_release, libknitro), Any, (Cint, Ptr{Cchar}), len, out)
    res = String(strip(String(convert(Vector{UInt8}, out)), '\0'))
    VersionNumber(split(res, " ")[2])
else
    VersionNumber(0, 0, 0) # Fake a version for AutoMerge
end

if KNITRO_VERSION != v"0.0.0" && KNITRO_VERSION < v"11.0"
    error(
        "You have installed version $KNITRO_VERSION of Artelys Knitro, which is not supported
  by KNITRO.jl. We require a Knitro version greater than 11.0.
  ",
    )
end

has_knitro() = IS_KNITRO_LOADED
knitro_version() = KNITRO_VERSION

include("libknitro.jl")
include("C_wrapper.jl")
include("callbacks.jl")
include("MOI_wrapper.jl")

end


text:
"""
    bothat(img; dims=coords_spatial(img), r=1)
    bothat(img, se)

Performs morphological bottom-hat transform for given image, i.e., `closing(img, se) - img`.

$(_docstring_se)

This bottom-hat transform, also known as _black_ top-hat transform, can be used to extract
small black elements and details from an image. To extract white details, the _white_
top-hat transform [`tophat`](@ref) can be used.

# Examples

```jldoctest; setup = :(using ImageMorphology)
julia> img = falses(7, 7); img[3:5, 3:5] .= true; img[4, 6] = true; img[4, 4] = false; img
7×7 BitMatrix:
 0  0  0  0  0  0  0
 0  0  0  0  0  0  0
 0  0  1  1  1  0  0
 0  0  1  0  1  1  0
 0  0  1  1  1  0  0
 0  0  0  0  0  0  0
 0  0  0  0  0  0  0

julia> Int.(bothat(img))
7×7 Matrix{$Int}:
 0  0  0  0  0  0  0
 0  0  0  0  0  0  0
 0  0  0  0  0  0  0
 0  0  0  1  0  0  1
 0  0  0  0  0  0  0
 0  0  0  0  0  0  0
 0  0  0  0  0  0  0

julia> Int.(bothat(img, strel_diamond(img))) # use diamond shape SE
7×7 Matrix{$Int}:
 0  0  0  0  0  0  0
 0  0  0  0  0  0  0
 0  0  0  0  0  0  0
 0  0  0  1  0  0  0
 0  0  0  0  0  0  0
 0  0  0  0  0  0  0
 0  0  0  0  0  0  0
```

See also [`bothat!`](@ref) for the in-place version.
"""
bothat(img; kwargs...) = bothat!(similar(img, maybe_floattype(eltype(img))), img, similar(img); kwargs...)
bothat(img, se) = bothat!(similar(img, maybe_floattype(eltype(img))), img, se, similar(img))

"""
    bothat!(out, img, buffer; [dims], [r])
    bothat!(out, img, se, buffer)

The in-place version of [`bothat`](@ref) with input image `img` and output image `out`. The
intermediate dilation result is stored in `buffer`.
"""
function bothat!(out, img, buffer; dims=coords_spatial(img), r=nothing)
    return bothat!(out, img, strel_box(img, dims; r), buffer)
end
function bothat!(out, img, se, buffer)
    closing!(out, img, se, buffer)
    @. out = out - img
    return out
end
bothat!(out::AbstractArray{<:Color3}, img, se, buffer) = throw(ArgumentError("color image is not supported"))


text:
#!/usr/bin/env julia
# Geof Hannigan

# Purpose of this script is to find all SAM file seconary alignments
# whose quality match their associated primary alignment.

# Need to bring in the relevant packages before use
# using Pkg
# Pkg.add("Bio")
# Pkg.add("ArgParse")

# Use the package BioAlignments
# https://biojulia.net/BioAlignments.jl/latest/hts-files.html#SAM-and-BAM-Records-1
using BioAlignments
using DelimitedFiles
using ArgParse

# Setup the command line arguments
s = ArgParseSettings(
	version = "0.1",
	add_version = true
	)
@add_arg_table s begin
	"-i", "--input"
		help = "Input file name."
		required = true
	"-o", "--output"
		help = "Output file name."
		required = true
	"-s", "--samout"
		help = "Output file name for sam format."
		required = true
	"-v", "--verbose"
        help = "Show verbose output (includes dictionary values)."
        action = :store_true
end

parsed_args = parse_args(ARGS, s)

println("Parsed args:")
for (arg,val) in parsed_args
    println("  $arg  =>  $val")
end

# Open the SAM file
reader = open(SAM.Reader, parsed_args["input"])

# Testing can look like this:
# reader = open(SAM.Reader, "test1.sam")

# Make a count hash to capture multi-mappers
# Here we do it only at the read level
counts = Dict{String,Int64}()
for record in reader
	if SAM.ismapped(record) && SAM.isnextmapped(record)
		# Get the name
		tname = SAM.tempname(record)
		tflag = SAM.flag(record)
		flagresult = collect(last(bitstring(tflag), 12))
		if Int64(flagresult[6]) == 49
			rpair = "First"
			# Make sure to include position for better ID
			tpos = SAM.position(record)
			trpos = SAM.nextposition(record)
		elseif Int64(flagresult[5]) == 49
			rpair = "Second"
			# Reverse for the other paired read
			trpos = SAM.position(record)
			tpos = SAM.nextposition(record)
		else
			# println("Looks like its not paired")
			rpair = "Fail"
		end

		# For the multimappers, we only want to look at the sequence read name
		# and whether it if first or second. Location doesn't matter.

		tcombo = string(tname, "-", rpair)

		# First or second in pair
		# println(tcombo)
		if tcombo in keys(counts)
			counts[tcombo] += 1
		else
			counts[tcombo] = 1
		end
	end
end

if parsed_args["verbose"]
	println("Printing Counts:")
	for (key,val) in counts
    	println("  $key  =>  $val")
	end
end

# Get a dictionary of the primary alignments
# Should only be one primary alignment per read
# Key is the read, value is the reference

reader = open(SAM.Reader, parsed_args["input"])
# reader = open(SAM.Reader, "test1.sam")

primaryd = Dict{String,Array}()
for record in reader
	if SAM.ismapped(record) && SAM.isnextmapped(record)
		# This needs to be defined so that the scope carries
		# through to the rest of the loop
		mapq = 255

		tname = SAM.tempname(record)
		# println(tname)
		refn = SAM.refname(record)
		# println(refn)
		# Exception handling in case mapping quality is 255,
		# which means it does not exist
		try
			mapq = Int64(SAM.mappingquality(record))
		catch
			# 255 means it does not exist
			mapq = 255
		end
		# println(mapq)
		cig = SAM.cigar(record)

		tflag = SAM.flag(record)
		flagresult = collect(last(bitstring(tflag), 12))
		if Int64(flagresult[6]) == 49
			rpair = "First"
			# Make sure to include position for better ID
			tpos = SAM.position(record)
			trpos = SAM.nextposition(record)
		elseif Int64(flagresult[5]) == 49
			rpair = "Second"
			# Reverse for the other paired read
			trpos = SAM.position(record)
			tpos = SAM.nextposition(record)
		else
			# println("Looks like its not paired")
			rpair = "Fail"
		end

		tcombo = string(tname, "-", rpair)

		# println(tcombo)
		if tcombo in keys(primaryd)
			# println("Warning: Potentially multiple primary alignments!")
		else
			if SAM.isprimary(record)
				# println("Adding")
				# println(mapq)
				primaryd[tcombo] = [refn, tpos, trpos, mapq, cig]
			else
				# println("Not a primary alignment.")
			end
		end
	end
end

if parsed_args["verbose"]
	println("Printing Primary Aligment Stats:")
	for (key,val) in primaryd
	    println("  $key  =>  $val")
	end
end

# Find matching best hits among each read
# Report the read and alignment pair
# reader = open(SAM.Reader, "test1.sam")
reader = open(SAM.Reader, parsed_args["input"])

# Count primary pairs
primarypairs = Dict{String,Int64}()

for record in reader
	# Only look at what mapped, always
	if SAM.ismapped(record) && SAM.isnextmapped(record)
		mapq = 255

		#Get the ID information
		tname = SAM.tempname(record)
		refn = SAM.refname(record)
		try
			mapq = Int64(SAM.mappingquality(record))
		catch
			# 255 means it does not exist
			mapq = 255
		end
		cig = SAM.cigar(record)

		tflag = SAM.flag(record)
		flagresult = collect(last(bitstring(tflag), 12))
		if Int64(flagresult[6]) == 49
			rpair = "First"
			# Make sure to include position for better ID
			tpos = SAM.position(record)
			trpos = SAM.nextposition(record)
		elseif Int64(flagresult[5]) == 49
			rpair = "Second"
			# Reverse for the other paired read
			trpos = SAM.position(record)
			tpos = SAM.nextposition(record)
		else
			# println("Looks like its not paired")
			rpair = "Fail"
		end

		tcombo = string(tname, "-", rpair, "-", tpos, "-", trpos)

		# We want to count the read-pair combos because we only want
		# to keep those with two (a full paired match)

		combotr = string(tname, "-", refn, "-", tpos, "-", trpos)
		multimapperid = string(tname, "-", rpair)

		# And establish the counts starting at zero
		if combotr in keys(primarypairs)
			#nothing
		else
			primarypairs[combotr] = 0
		end

		# println("Searching ", tcombo, "-", refn)
		# println("Multimapper ID is ", multimapperid)
		# Check if the read is a multimapper
		if counts[multimapperid] > 1
			# Check if it is a secondary hit
			if SAM.isprimary(record)
				# println(tcombo, "---", refn)
				# Keep track of primary alignment counts
				primarypairs[combotr] += 1
			else
				# Determine if this hit is the same as the
				# primary alignment of the read
				aref = primaryd[multimapperid][1]
				refstart = primaryd[multimapperid][2]
				refend = primaryd[multimapperid][3]
				aqual = primaryd[multimapperid][4]
				acigar = primaryd[multimapperid][5]
				if isequal(tpos, refstart) && isequal(trpos, refend)
						# nothing
				else
					if isequal(mapq, aqual)
						if isequal(cig, acigar)
							# println(tcombo, "---", refn)
							primarypairs[combotr] += 1
						end
					end
				end
			end
		elseif counts[multimapperid] == 1
			# println(tcombo, "---", refn)
			primarypairs[combotr] += 1
		end
	end
end

if parsed_args["verbose"]
	println("Printing Primary Pairs:")
	for (key,val) in primarypairs
	    println("  $key  =>  $val")
	end
end

# Print output to a 2 dim array for now
outarray = Array{Char}(undef, 0, 4)
filterarray = Array{String}(undef, 0, 1)

for (key, value) in primarypairs
	if isequal(value, 2)
		karray = split(key, "-")
		tread = karray[1]
		tref = karray[2]
		spos = karray[3]
		srpos = karray[4]
		# println(karray)
		global outarray = [outarray; [tread tref spos srpos]]
		# push!(outarray, karray)
		# Also make an array of values that pass the two test
		global filterarray = [filterarray; [key]]
	end
end

# Write the final array to tsv file for reference
# writedlm( "testout.txt",  outarray, '\t')
writedlm(parsed_args["output"],  outarray, '\t')

# We also want to write the sam file
# First get in the sam file and pull the header
# reader = open(SAM.Reader, "test1.sam")
reader = open(SAM.Reader, parsed_args["input"])
samhead = header(reader)

samw = SAM.Writer(open(parsed_args["samout"], "w"), samhead)

for record in reader
	# Only look at what mapped, always
	if SAM.ismapped(record) && SAM.isnextmapped(record)
		mapq = 255

		#Get the ID information
		tname = SAM.tempname(record)
		refn = SAM.refname(record)
		try
			mapq = Int64(SAM.mappingquality(record))
		catch
			# 255 means it does not exist
			mapq = 255
		end

		tflag = SAM.flag(record)
		flagresult = collect(last(bitstring(tflag), 12))
		if Int64(flagresult[6]) == 49
			rpair = "First"
			# Make sure to include position for better ID
			tpos = SAM.position(record)
			trpos = SAM.nextposition(record)
		elseif Int64(flagresult[5]) == 49
			rpair = "Second"
			# Reverse for the other paired read
			trpos = SAM.position(record)
			tpos = SAM.nextposition(record)
		else
			# println("Looks like its not paired")
			rpair = "Fail"
		end

		combotr = string(tname, "-", refn, "-", tpos, "-", trpos)

		# Does the combo string match the array

		if combotr in filterarray
			# println(filterarray)
			# println(combotr)
			write(samw, record)
		end

	end
end

close(samw)


text:
if get(ARGS,1,"Array") == "CuArray"
    using CUDA
    const ArrayType = CUDA.CuArray
else
    const ArrayType = Array
end

include("testhelper.jl")
include("problems.jl")

include("integrator.jl")
include("convergence.jl")
include("callbacks.jl")

#=
@testset "ODE Tests: Basic" begin
    runmpi(joinpath(@__DIR__, "basic.jl"))
end
=#
# FIXME: Should consolodate all convergence tests into single
# testset --- this test is slightly redundant
# @testset "ODE Tests: Convergence" begin
#     runmpi(joinpath(@__DIR__, "ode_tests_convergence.jl"))
# end


text:
"""
    isinsertable(dict::AbstractDictionary)

Return `true` if `dict` supports the insertable interface, or `false` otherwise. The primary
functions `dict` needs to implement for the insertable interface are:

 * `insert!(dict, i, value)` - add a new `value` at index `i` (will error if index exists)
 * `delete!(dict, i)` - remove element at index `i` (will error if index does not exist)

Functions `get!`, `set!` and `unset!` are also provided for common operations where you are
not sure if an index exists or not. New insertable dictionaries are primarily generated via
the `empty` function.

See also `issettable` and `istokenizable`.
"""
isinsertable(::AbstractDictionary) = false

"""
    isinsertable(indices::AbstractIndices)

Return `true` if `indices` supports the insertable interface, or `false` otherwise. The
primary functions a map needs to implement for the insertable interface are:

 * `insert!(indices, i)` - add new index `i` to `indices` (will error if index exists)
 * `delete!(indices, i)` - remove an existing index `i` from `indices` (will error if index does not exist).

Functions `set!` and `unset!` are also provided for common operations where you are not sure
if an index exists or not. New insertable indices are primarily generated via the `empty`
function.
"""
isinsertable(::AbstractIndices) = false

### Underlying token interface functions

"""
    gettoken!(dict, i)

Return the tuple `(hadindex, token)`, where `hadindex` is `true` if `i` previously existed
in `dict`, or `false` otherwise (in which case `dict` was mutated to insert a slot for the
new index `i`). The `token` may be used to retrieve a value using the `gettokenvalue` or set
a corresponding value via the `settokenvalue!`.

See also `gettoken` and `deletetoken!`.
"""
@propagate_inbounds function gettoken!(d::AbstractDictionary{I}, i) where {I}
    return gettoken!(d, convert(I, i))
end

function gettoken!(d::AbstractDictionary{I}, i::I) where {I}
    if isinsertable(d)
        error("gettoken! needs to be defined for insertable dictionary: $(typeof(d))")
    else
        error("dictionary not insertable: $(typeof(d))")
    end
end

function gettoken!(d::AbstractIndices{I}, i::I) where {I}
    if isinsertable(d)
        error("gettoken! needs to be defined for insertable indices: $(typeof(d))")
    else
        error("indices not insertable: $(typeof(d))")
    end
end

"""
    deletetoken!(dict, token)

Remove the slot of the dictionary at `token`.

See also `gettoken` and `gettoken!`.
"""
function deletetoken!(d::AbstractDictionary, token)
    if isinsertable(d)
        error("deletetoken! needs to be defined for insertable dictionary: $(typeof(d))")
    else
        error("dictionary not insertable: $(typeof(d))")
    end
end

function deletetoken!(d::AbstractIndices, token)
    if isinsertable(d)
        error("deletetoken! needs to be defined for insertable indices: $(typeof(d))")
    else
        error("indices not insertable: $(typeof(d))")
    end
end

### User-facing scalar insertion/deletion

"""
    insert!(indices::AbstractIndices, i)

Insert the new index `i` into `indices`. An error is thrown if `i` already exists.
"""
@propagate_inbounds function Base.insert!(indices::AbstractIndices{I}, i) where {I}
    i2 = convert(I, i)
    if !isequal(i, i2)
        throw(ArgumentError("$i is not a valid key for type $I"))
    end
    return insert!(indices, i2)
end

function Base.insert!(indices::AbstractIndices{I}, i::I) where {I}
    (hadindex, token) = gettoken!(indices, i)
    if hadindex
        throw(IndexError("Indices already contains index: $i"))
    end
    return indices
end

"""
    insert!(dict::AbstractDictionary, i, value)

Insert the `value` at new index `i` into `dict`. An error is thrown if index `i` already
exists.

Hint: Use `setindex!` to update an existing value, and `set!` to perform an "upsert"
(update-or-insert) operation.
"""
function Base.insert!(d::AbstractDictionary{I}, i, value) where {I}
    i2 = convert(I, i)
    if !isequal(i, i2)
        throw(ArgumentError("$i is not a valid key for type $I"))
    end
    return insert!(d, i2, value)
end

function Base.insert!(d::AbstractDictionary{I, T}, i::I, value) where {I, T}
    return insert!(d, i, convert(T, value))
end

function Base.insert!(d::AbstractDictionary{I, T}, i::I, value::T) where {I, T}
    (hadtoken, token) = gettoken!(d, i)
    if hadtoken
        throw(IndexError("Dictionary already contains index: $i"))
    end
    @inbounds settokenvalue!(d, token, value)
    return d
end

# Since `AbstractIndices <: AbstractDictionary` we should still obey the supertype's interface where possible...
function Base.insert!(indices::AbstractIndices{I}, i1::I, i2::I) where {I}
    if isinsertable(indices)
        if isequal(i1, i2)
            return insert!(indices, i1)
        else
            error("Attempted to set distinct key-value pair ($i1, $i2) to indices: $(typeof(indices))")
        end
    else
        error("indices not insertable: $(typeof(indices))")
    end
end

"""
    set!(dict::AbstractDictionary, i, value)

Update or insert the `value` at index `i` into `dict`. Sometimes referred to as an "upsert"
operation.

Hint: Use `setindex!` to exclusively update an existing value, and `insert!` to exclusively
insert a new value. See also `get!`.
"""
@propagate_inbounds function set!(d::AbstractDictionary{I}, i, value) where {I}
    i2 = convert(I, i)
    if !isequal(i, i2)
        throw(ArgumentError("$i is not a valid key for type $I"))
    end
    return set!(d, i2, value)
end

function set!(d::AbstractDictionary{I, T}, i::I, value) where {I, T}
    return set!(d, i, convert(T, value))
end

function set!(d::AbstractDictionary{I, T}, i::I, value::T) where {I, T}
    (hadtoken, token) = gettoken!(d, i)
    @inbounds settokenvalue!(d, token, value)
    return d
end

# Since `AbstractIndices <: AbstractDictionary` we should still obey the supertype's interface where possible...
function set!(indices::AbstractIndices{I}, i1::I, i2::I) where {I}
    if isinsertable(indices)
        if isequal(i1, i2)
            return set!(indices, i1)
        else
            error("Attempted to set distinct key-value pair ($i1, $i2) to indices: $(typeof(indices))")
        end
    else
        error("indices not insertable: $(typeof(indices))")
    end
end

"""
    set!(indices::AbstractIndices, i)

Insert a new value `i` into `indices` if it doesn't exist, or do nothing otherwise.
"""
@propagate_inbounds function set!(indices::AbstractIndices{I}, i) where {I}
    i2 = convert(I, i)
    if !isequal(i, i2)
        throw(ArgumentError("$i is not a valid key for type $I"))
    end
    return set!(indices, i2)
end

function set!(indices::AbstractIndices{I}, i::I) where {I}
    gettoken!(indices, i)
    return indices
end

"""
    get!(dict::AbstractDictionary, i, default)

Return the value `dict[i]` if index `i` exists. Otherwise, a new index `i` is inserted and
set to `default`, which is returned.

See also `get`, `set!`.
"""
@propagate_inbounds function Base.get!(d::AbstractDictionary{I}, i, default) where {I}
    i2 = convert(I, i)
    if !isequal(i, i2)
        throw(ArgumentError("$i is not a valid key for type $I"))
    end
    return get!(d, i2, default)
end

function Base.get!(d::AbstractDictionary{I, T}, i::I, default) where {I, T}
    return get!(d, i, convert(T, default))
end

function Base.get!(d::AbstractDictionary{I, T}, i::I, default::T) where {I, T}
    (hadindex, token) = gettoken!(d, i)
    if hadindex
        return @inbounds gettokenvalue(d, token)
    else
        @inbounds settokenvalue!(d, token, default)
        return default
    end
end

"""
    get!(f::Union{Function, Type}, dict::AbstractDictionary, i)

Return the value `dict[i]` if index `i` exists. Otherwise, a new index `i` is inserted and
set to the value `f()`, which is returned.
"""
function Base.get!(f::Callable, d::AbstractDictionary{I}, i) where {I}
    i2 = convert(I, i)
    if !isequal(i, i2)
        throw(ArgumentError("$i is not a valid key for type $I"))
    end
    get!(f, d, i2)
end

function Base.get!(f::Callable, d::AbstractDictionary{I, T}, i::I) where {I, T}
    (hadindex, token) = gettoken!(d, i)
    if hadindex
        return @inbounds gettokenvalue(d, token)
    else
        default = convert(T, f())
        @inbounds settokenvalue!(d, token, default)
        return default
    end
end

"""
    delete!(indices::AbstractIndices, i)

Delete the index `i` from `indices`. An error is thrown if `i` does not exist.

    delete!(dict::AbstractDictionary, i)

Delete the index `i` from `dict`. An error is thrown if `i` does not exist.

See also `unset!`, `insert!`.
"""
@propagate_inbounds function Base.delete!(d::AbstractDictionary{I}, i) where {I}
    i2 = convert(I, i)
    if !isequal(i, i2)
        throw(ArgumentError("$i is not a valid key for type $I"))
    end
    return delete!(d, i2)
end

function Base.delete!(d::AbstractDictionary{I}, i::I) where {I}
    (hasindex, token) = gettoken(d, i)
    if !hasindex
        throw(IndexError("Index doesn't exist: $i"))
    end
    @inbounds deletetoken!(d, token)
    return d
end

"""
    unset!(indices::AbstractIndices, i)

Delete the index `i` from `indices` if it exists, or do nothing otherwise.

    unset!(dict::AbstractDictionary, i)

Delete the index `i` from `dict` if it exists, or do nothing otherwise.

See also `delete!`, `set!`.
"""
@propagate_inbounds function unset!(indices::AbstractDictionary{I}, i) where {I}
    i2 = convert(I, i)
    if !isequal(i, i2)
        throw(ArgumentError("$i is not a valid key for type $I"))
    end
    return unset!(indices, i2)
end

function unset!(d::AbstractDictionary{I}, i::I) where {I}
    (hasindex, token) = gettoken(d, i)
    if hasindex
        @inbounds deletetoken!(d, token)
    end
    return d
end

### Non-scalar insertion/deletion
function Base.merge!(d::AbstractDictionary, d2::AbstractDictionary)
    for (i, v) in pairs(d2)
        set!(d, i, v)
    end
    return d
end

if isdefined(Base, :mergewith) # Julia 1.5+
    function Base.mergewith!(combiner::Callable, d::AbstractDictionary, d2::AbstractDictionary)
        for (i, v) in pairs(d2)
            (hasindex, token) = gettoken!(d, i)
            if hasindex
                @inbounds settokenvalue!(d, token, combiner(gettokenvalue(d, token), v))
            else
                @inbounds settokenvalue!(d, token, v)
            end
        end
        return d
    end
end

# TODO some kind of exclusive merge (throw on key clash like `insert!`)

# TODO: deletemany!, unsetmany! (which is setdiff! on indices)

### Indices ("sets") versions of above

function Base.union!(s1::AbstractIndices, itr)
    # Optimized to handle repeated values in `itr` - e.g. if `itr` is already sorted
    x = iterate(itr)
    if x === nothing
        return s1
    end
    (i, s) = x
    set!(s1, i)
    i_old = i
    while x !== nothing
        (i, s) = x
        !isequal(i, i_old) && set!(s1, i)
        i_old = i
        x = iterate(itr, s)
    end
    return s1
end

function Base.intersect!(s1::AbstractIndices, s2::AbstractIndices)
    return filter!(in(s2), s1)
end

function Base.setdiff!(s1::AbstractIndices, s2::AbstractIndices)
    for i in s2
        unset!(s1, i)
    end
    return s1
end

function Base.symdiff!(s1::AbstractIndices, s2::AbstractIndices)
    for i in s2
        (hastoken, token) = gettoken!(s1, i)
        if hastoken
            @inbounds deletetoken!(s1, token)
        end
    end
    return s1
end

## Filtering

# These generic implementations are gimped.

# # `filter!` is basically a programmatic version of `intersect!`.
# function Base.filter!(pred, indices::AbstractIndices)
#     for i in copy(indices)
#         if !pred(i)
#             delete!(indices, i)
#         end
#     end
#     return indices
# end

# # Dictionary version is similar
# function Base.filter!(pred, dict::AbstractDictionary)
#     for (i, v) in pairs(copy(dict))
#         if !pred(v)
#             delete!(dict, i)
#         end
#     end
#     return dict
# end

# This implementation is faster when deleting indices does not invalidate tokens/iteration,
# and is opt-in only. Works for both dictionaries and indices
function Base.unsafe_filter!(pred, dict::AbstractDictionary)
    for token in tokens(dict)
        @inbounds if !pred(gettokenvalue(dict, token))
            deletetoken!(dict, token)
        end
    end
end

# Factories for insertable Dictionaries/Indices

"""
    empty(inds::AbstractIndices, I::Type)
    empty(dict::AbstractDictionary, I::Type)

Return an empty, insertable `AbstractIndices` of element type `I` (even when the first
argument is a dictionary). The default container is `Indices{I}`, but the output may
depend on the first argument.

    empty(inds::AbstractIndices)

Return an empty, insertable `AbstractIndices` of element type `eltype(inds)`.
"""
Base.empty(inds::AbstractIndices) = empty(inds, eltype(inds))

"""
    empty(inds::AbstractIndices, I::Type, T::Type)
    empty(dict::AbstractDictionary, I::Type, T::Type)

Return an empty, insertable `AbstractDictionary` of with indices of type `I` and elements
of type `T` (even when the first argument is are indices). The default container is
`Dictionary{I}`, but the output may depend on the first argument.

    empty(dict::AbstractDictionary)

Return an empty, insertable `AbstractDictionary` with indices of type `keytype(dict)` and
elements of type `eltype(inds)`.
"""
Base.empty(d::AbstractDictionary) = empty(keys(d), keytype(d), eltype(d))

Base.empty(d::AbstractDictionary, ::Type{I}) where {I} = empty(keys(d), I)

Base.empty(::AbstractDictionary, ::Type{I}, ::Type{T}) where {I, T} = Dictionary{I, T}()


text:
module TestStoreApi

using ..MyPetStore
using Swagger
using Test
using Dates

function test(uri)
    println("testing StoreApi...")
    client = Swagger.Client(uri)
    api = StoreApi(client)

    println("   - getInventory")
    inventory = getInventory(api)
    @test isa(inventory, Dict{String,Int32})
    @test !isempty(inventory)

    println("   - placeOrder")
    @test_throws Swagger.ValidationException Order(; id=10, petId=10, quantity=2, shipDate=DateTime(2017, 03, 12), status="invalid_status", complete=false)
    order = Order(; id=10, petId=10, quantity=2, shipDate=DateTime(2017, 03, 12), status="placed", complete=false)
    neworder = placeOrder(api, order)
    @test neworder.id == 10

    println("   - getOrderById")
    @test_throws Swagger.ValidationException getOrderById(api, 0)
    order = getOrderById(api, 10)
    @test isa(order, Order)
    @test order.id == 10

    println("   - deleteOrder")
    @test deleteOrder(api, 10) == nothing

    nothing
end

end # module TestStoreApi


text:
using IterativelyHypergraphEmbedding
using SimpleHypergraphs
using LinearAlgebra
using Test

@testset "IterativelyHypergraphEmbedding.jl" begin
  h = Hypergraph{Int, Int}(4, 0)
  add_hyperedge!(h, vertices=Dict(1=>1, 2=>1, 3=>1))
  add_hyperedge!(h, vertices=Dict(3=>1, 4=>1))
  dims = 10
  hn_vec, he_vec = iteratively_hypergraph_embedding(h, dims)
  @test size(hn_vec) == (dims, nhv(h))
  @test size(he_vec) == (dims, nhe(h))

  # Test that the hypernodes and hyperedges are embedded in the hypersphere.
  for vec in eachcol(hn_vec)
    @test abs(norm(vec) - 1) < 1e-8
  end
  for vec in eachcol(he_vec)
    @test abs(norm(vec) - 1) < 1e-8
  end

end


text:
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using ViewerGL
GL = ViewerGL

function randlines(n=300, t=0.4)
	#n = 100 #1000 #1000 #20000
	#t = 0.4 #0.15 #0.4 #0.15
	V = zeros(Float64,2,2*n)
	EV = [zeros(Int64,2) for k=1:n]
	for k=1:n
		v1 = rand(Float64,2)
		v2 = rand(Float64,2)
		vm = (v1+v2)/2
		transl = rand(Float64,2)
		V[:,k] = (v1-vm)*t + transl
		V[:,n+k] = (v2-vm)*t + transl
		EV[k] = [k,n+k]
	end

	V = GL.normalize2(V)
	model = (V,EV)
	Sigma = Lar.spaceindex(model)

	model = V,EV;
	W,EW = Lar.fragmentlines(model);
	U,EVs = Lar.biconnectedComponent((W,EW::Lar.Cells));
	EV = convert(Lar.Cells, cat(EVs))
	V,FVs,EVs = Lar.arrange2D(U,EV)
end

# ////////////////////////////////////////////////////////////
# generation of 2D arrangement
V,FVs,EVs = randlines()

# native OpenGL visualization
GL.VIEW(GL.GLExplode(V,FVs,1.2,1.2,1.2,1));
GL.VIEW(GL.GLExplode(V,FVs,1.2,1.2,1.2,3,1));
GL.VIEW(GL.GLExplode(V,FVs,1.2,1.2,1.2,99,1));
GL.VIEW(GL.GLExplode(V,FVs,1.,1.,1.,99,1));
GL.VIEW(GL.GLExplode(V,EVs,1.2,1.2,1.2,1,1));


text:
"""
```
    ContrastStretching <: AbstractHistogramAdjustmentAlgorithm
    ContrastStretching(; t = 0.5,  slope = 1.0)

    adjust_histogram([T,] img, f::ContrastStretching)
    adjust_histogram!([out,] img, f::ContrastStretching)
```

Returns an image where intensities below `t` are compressed into a narrower
range of dark intensities, and values above `t` are compressed into a narrower
band of light intensities.

# Details

Contrast stretching is a transformation that  enhances or reduces (for `slope` >
1 or < 1, respectively) the contrast near saturation (0 and 1).
It is given by the relation
```math
f(x) = \\frac{1}{1 + \\left(\\frac{t}{x} \\right)^s}, \\; s \\in \\mathbb{R},
```
where ``s`` represents the `slope` argument.

# Options

Various options for the parameters of the `adjust_histogram` and
`ContrastStretching` type are described in more detail below.

## Choices for `img`

The function can handle a variety of input types. The returned
image depends on the input type.

For colored images, the input is converted to the
[YIQ](https://en.wikipedia.org/wiki/YIQ)  type and the intensities of the Y
channel are stretched to the specified range. The modified Y channel is then
combined with the I and Q channels and the resulting image converted to the same
type as the input.

## Choice for `t`

The value of `t` needs to be in the unit interval. If left unspecified a
default value of 0.5 is utilized.

## Choice for `slope`

The value of `slope` can be any real number. If left unspecified a
default value of 1.0 is utilized.

# Example

```julia
using ImageContrastAdjustment, ImageView, TestImages

img = testimage("mandril_gray")
ret = adjust_histogram(img, ContrastStretching(t = 0.6, slope = 3))

```

# References
1. Gonzalez, R. C., Woods, R. E., & Eddins, S. L. (2004). *Digital image processing using MATLAB* (Vol. 624). Upper Saddle River, New Jersey: Pearson-Prentice-Hall.

"""
Base.@kwdef struct ContrastStretching{T₁ <: Union{Real,AbstractGray},
                                      T₂ <: Union{Real,AbstractGray}}  <: AbstractHistogramAdjustmentAlgorithm
     t::T₁ = 0.5
     slope::T₂ = 1.0
end

function (f::ContrastStretching)(out::GenericGrayImage, img::GenericGrayImage)
    T = eltype(out)
    ϵ = eps(T)
    out .= img
    map!(out,out) do val
        if isnan(val)
            return val
        else
            newval = contrast_stretch(val, f.t, f.slope, ϵ)
            return  T <: Integer ? round(Int, newval ) : newval
        end
    end
end

function (f::ContrastStretching)(out::AbstractArray{<:Color3}, img::AbstractArray{<:Color3})
    T = eltype(img)
    yiq = convert.(YIQ, img)
    yiq_view = channelview(yiq)
    adjust_histogram!(view(yiq_view,1,:,:), f)
    out .= convert.(T, yiq)
end

(f::ContrastStretching)(out::GenericGrayImage, img::AbstractArray{<:Color3}) =
    f(out, of_eltype(Gray, img))

function contrast_stretch(x, t, s, ϵ)
    return 1 / (1 + (t / (x+ϵ))^s)
end


text:
module BangBangImpl

# TODO: Let BangBang depend on MicroCollections.
using BangBang: BangBang, append!!, push!!, setindex!!, union!!
using BangBang.Experimental: mergewith!!

using ..MicroCollections:
    AbstractMicroDict,
    AbstractMicroSet,
    AbstractMicroVector,
    SingletonVector,
    UndefArray,
    upcast

BangBang.append!!(A::AbstractMicroVector, B::AbstractVector) = append!!(upcast(A), B)
# TODO: Avoid allocating array twice if `upcast(A)` is going to be
# widen anyway.

BangBang.push!!(A::AbstractMicroVector, x) = append!!(A, SingletonVector((x,)))

struct _NoValue end

function BangBang.Experimental.mergewith!!(combine, dict::AbstractMicroDict, other)
    udict = foldl(push!!, pairs(other), init = upcast(dict))
    if length(dict) == 0
        return udict
    else
        k, v = first(dict)
        vo = get(other, k, _NoValue())
        if vo isa _NoValue
            return udict
        else
            return setindex!!(udict, combine(v, vo), k)
        end
    end
end

BangBang.setindex!!(dict::AbstractMicroDict, v, k) = setindex!!(upcast(dict), v, k)

BangBang.union!!(A::AbstractMicroSet, B) = union!!(upcast(A), B)

BangBang.implements(::typeof(resize!), ::Type{<:UndefArray}) = false

end  # module


text:
module PIPG
using StaticArrays
using SparseArrays
using LinearAlgebra
using LazyArrays

include("solver.jl")
include("preconditioning.jl")
include("MOI_wrapper/MOI_wrapper.jl")

# Write your package code here.

end


text:
#########################
# Spin Operator Methods #
#########################
immutable HalfSpin{T<:Integer}
    val::T
    HalfSpin(val) = val >= 0 ? new(val) : error("Spin must be positive")
end

HalfSpin{T}(val::T) = HalfSpin{T}(val)

spin(s::HalfSpin) = s
spin(j::Integer) = HalfSpin(2*j)
function spin(j)
    if isinteger(j)
        return spin(@compat(Int(j))) # j = 0.0, 1.0, 2.0...
    elseif isinteger(2*j)
        return HalfSpin(@compat(Int(2*j))) # j = .5, 1.5, 2.5...
    else
        error("Spin must be a multiple of 1/2.")
    end
end

spin_value(s::HalfSpin) = s.val/2

#####################
# Auxiliary Methods #
#####################
spin_coeffs(j, scalar) = [scalar*sqrt(j*(j+1)-(x+1)*x) for x in j-1:-1:-j]

function spinjyx_mat(s::HalfSpin, a, b)
    N = s.val + 1
    j = spin_value(s)
    up_coeffs = spin_coeffs(j, a)
    down_coeffs = b == 1 ? up_coeffs : b * up_coeffs
    return spdiagm((up_coeffs, down_coeffs), (1, -1), N, N)
end

spinjx_mat(j) = spinjyx_mat(spin(j), .5, 1)
spinjy_mat(j) = spinjyx_mat(spin(j), -.5im, -1)

function spinjpm_mat(s::HalfSpin, k)
    N = s.val + 1
    return spdiagm(spin_coeffs(spin_value(s), 1), k, N, N)
end

spinjp_mat(j) = spinjpm_mat(spin(j), 1)
spinjm_mat(j) = spinjpm_mat(spin(j), -1)

function spinjz_mat(s::HalfSpin)
    j = spin_value(s)
    return spdiagm(j:-1:-j)
end

spinjz_mat(j) = spinjz_mat(spin(j))

################################
# Jx, Jy, Jz, Jp, Jm operators #
################################
# These Wikipedia pages give a nice overview of spin/ladder operators:
# http://en.wikipedia.org/wiki/Spin_%28physics%29#Operator
# http://en.wikipedia.org/wiki/Ladder_operator#Angular_momentum
spinjx(j) = QuArray(spinjx_mat(j))
spinjy(j) = QuArray(spinjy_mat(j))
spinjz(j) = QuArray(spinjz_mat(j))
spinjp(j) = QuArray(spinjp_mat(j))
spinjm(j) = QuArray(spinjm_mat(j))

############################
#  Pauli spin 1/2 matrices #
############################
const sigmax = 2.0 * full(spinjx(1/2))
const sigmay = 2.0 * full(spinjy(1/2))
const sigmaz = 2.0 * full(spinjz(1/2))
const sigmap = full(sigmax) + im * full(sigmay)
const sigmam = full(sigmax) - im * full(sigmay)

# TODO : unicode sigmay, sigmaz
const σₓ = sigmax

export spinjx,
    spinjy,
    spinjz,
    spinjp,
    spinjm,
    sigmax,
    sigmay,
    sigmaz,
    sigmap,
    sigmam,
    σₓ


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

name = "Xorg_inputproto"
version = v"2.3.2"

# Collection of sources required to build inputproto
sources = [
    ArchiveSource("https://www.x.org/archive/individual/proto/inputproto-$(version).tar.bz2",
               "893a6af55733262058a27b38eeb1edc733669f01d404e8581b167f03c03ef31d")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/inputproto-*/
CPPFLAGS="-I${prefix}/include"
# When compiling for things like ppc64le, we need newer `config.sub` files
update_configure_scripts
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} --enable-malloc0returnsnull=no
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [p for p in supported_platforms() if Sys.islinux(p) || Sys.isfreebsd(p)]

products = Product[
]

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

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

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
l(s)=println(s)

function main()
  l("Julia")
end

main()

max_stars_repo_name: justletterh/el

max_issues_repo_name: justletterh/el

max_forks_repo_name: justletterh/el


text:
"""
    Polysegment <: AbstractImageBinarizationAlgorithm
    Polysegment()

    binarize([T,] img, f::Polysegment)
    binarize!([out,] img, f::Polysegment)

Uses the *polynomial segmentation* technique to group the image pixels
into two categories (foreground and background).

# Output

Return the binarized image as an `Array{Gray{T}}` of size `size(img)`. If
`T` is not specified, it is inferred from `out` and `img`.

# Details

The approach involves constructing a univariate second-degree polynomial
such that the two roots of the polynomial represent the graylevels
of two cluster centers (i.e the foreground and background). Pixels are then
assigned to the foreground or background depending on which cluster center
is closest.

# Arguments

The function argument is described in more detail below.

##  `img::AbstractArray`

The image that needs to be binarized. The image is automatically converted
to `Gray`.

# Example

Binarize the "cameraman" image in the `TestImages` package.

```julia
using TestImages, ImageBinarization

img = testimage("cameraman")
img_binary = binarize(img, Polysegment())
```

## Reference

1. R. E. Vidal, "Generalized Principal Component Analysis (GPCA): An Algebraic Geometric Approach to Subspace Clustering and Motion Segmentation." Order No. 3121739, University of California, Berkeley, Ann Arbor, 2003.
"""
struct Polysegment <: AbstractImageBinarizationAlgorithm end

function (f::Polysegment)(out::GenericGrayImage, img::GenericGrayImage)
  # Construct data matrix for second-degree polynomial (Equation 2.3) in [1].
  x = vec(img)
  Lₙ = hcat(ones(length(x)), x, x.^2)
  F = svd(Lₙ)
  c = F.Vt[end,:]
  p = Polynomial(vec(c))
  μ₁, μ₂ = sort(roots(p))
  # Binarize the image.
  @simd for i in CartesianIndices(img)
    val = img[i]
    out[i] = (val-μ₁)^2 < (val-μ₂)^2 ? 0 : 1
  end
  out
end

function (f::Polysegment)(out::GenericGrayImage, img::AbstractArray{<:Color3})
    # map `img` to grayspace while keeping the storage type
    # TODO: this mapping operation can be done lazily to reduce memory allocation
    GT = base_color_type(eltype(out)){eltype(eltype(img))}
    f(out, GT.(img))
end


text:
# FWI on Overthrust model using minConf library
# Author: Philipp Witte, pwitte@eoas.ubc.ca
# Date: December 2017
#

using Statistics, Random, LinearAlgebra
using JUDI, HDF5, SegyIO, PyPlot, IterativeSolvers

# Load starting model
n,d,o,m0 = read(h5open("../../data/overthrust_model.h5","r"), "n", "d", "o", "m0")
model0 = Model((n[1],n[2]), (d[1],d[2]), (o[1],o[2]), m0)

# Bound constraints
v0 = sqrt.(1 ./ model0.m)
vmin = ones(Float32,model0.n) .* 1.3f0
vmax = ones(Float32,model0.n) .* 6.5f0
vmin[:,1:21] .= v0[:,1:21]   # keep water column fixed
vmax[:,1:21] .= v0[:,1:21]

# Slowness squared [s^2/km^2]
mmin = vec((1f0 ./ vmax).^2)
mmax = vec((1f0 ./ vmin).^2)

# Load data
block = segy_read("../../data/overthrust_shot_records.segy")
d_obs = judiVector(block)

# Set up wavelet
src_geometry = Geometry(block; key="source")
wavelet = ricker_wavelet(src_geometry.t[1],src_geometry.dt[1],0.008f0)  # 8 Hz wavelet
q = judiVector(src_geometry,wavelet)

############################### FWI ###########################################

# Set up operators
ntComp = get_computational_nt(q.geometry,d_obs.geometry,model0) # no. of computational time steps
info = Info(prod(model0.n),d_obs.nsrc,ntComp)
Pr = judiProjection(info,d_obs.geometry)
Ps = judiProjection(info,q.geometry)
F = judiModeling(info,model0)
J = judiJacobian(Pr*F*Ps',q)

# Optimization parameters
maxiter = 10
maxiter_GN = 5
fhistory_GN = zeros(Float32,maxiter)
proj(x) = reshape(median([vec(mmin) vec(x) vec(mmax)]; dims=2),model0.n)

# Gauss-Newton method
for j=1:maxiter
    println("Iteration: ",j)

    # # Model predicted data for subset of sources
    d_pred = Pr*F*Ps'*q
    fhistory_GN[j] = .5f0*norm(d_pred - d_obs)^2

    # GN update direction
    p = lsqr!(similar(model0.m), J, d_pred - d_obs; maxiter=maxiter_GN, verbose=true)

    # update model and bound constraints
    model0.m .= proj(model0.m .- reshape(p, model0.n))  # alpha=1
end

figure(); imshow(sqrt.(1f0./model0.m)'); title("FWI with Gauss-Newton")

max_stars_repo_name: nogueirapeterson/JUDI

max_issues_repo_name: nogueirapeterson/JUDI

max_forks_repo_name: nogueirapeterson/JUDI


text:
import HTTP

abstract type Authentication end

abstract type DICOMType end

struct BaseURL
    uri::HTTP.URI
end

struct Client{A <: Authentication}
    base_url::BaseURL
    auth::A
end

# struct Credential
#     secret_buffer::Base.SecretBuffer
# end

struct AnonymousAuth <: Authentication end


text:
# core methods

###
### neighbors
###
@propagate_inbounds function outneighbors(g::AdjacencyMap, v::Integer)
    cptr = getcolptr(g)
    rvals = rowvals(g)
    out = unsafe_inneighbors(nvertices(g), cptr, rvals, v)
    @inbounds for i in cptr[v]:cptr[v + 1] - 1
        pushfirst!(out, i)
    end
    return out
end

inneighbors(g::AdjMap, v::Integer) = outneighbors(g, v)
function inneighbors(g::AdjDiMap, v::Integer)
    unsafe_inneighbors(nvertices(g), getcolptr(g), rowvals(g), v)
end
function unsafe_inneighbors(N, cptr, rval, v)
    nzinds = Int[]
    ptrI = 1
    @inbounds for j in OneTo(N)
        rowI = v
        ptrA = Int(cptr[j])
        stopA = Int(cptr[j+1]-1)
        if ptrA <= stopA
            if rval[ptrA] <= rowI
                ptrA = searchsortedfirst(rval, rowI, ptrA, stopA, Base.Order.Forward)
                if ptrA <= stopA && rval[ptrA] == rowI
                    push!(nzinds, j)
                end
            end
            ptrI += 1
        end
    end
    return nzinds
end

neighbors(g::AdjMap, v::Integer) = outneighbors(g, v)

## all_neighbors
all_neighbors(g::AdjMap, v::Integer) = outneighbors(g, v)
all_neighbors(g::AdjDiMap, v::Integer) = union(outneighbors(g, v), inneighbors(g, v))

## common_neighbors
function common_neighbors(g::AdjacencyMap, u::Integer, v::Integer)
    intersect(neighbors(g, u), neighbors(g, v))
end

###
### density
###
function density(g::AdjDiMap)
    N = nvertices(g)
    return nedges(g) / (N * (N - 1))
end
function density(g::AdjMap)
    N = nvertices(g)
    return (2 * nedges(g)) / (N * (N - 1))
end

## num_self_loops
function num_self_loops(g::AdjacencyMap)
    if nvertices(g) == 0
        return 0
    else
        return sum(v -> has_edge(g, v, v), vertices(g))
    end
end

## degree_histogram
function degree_histogram(g::AdjacencyMap, degfn=degree)
    hist = Dict{eltype(g),Int}()
    for v in vertices(g)        # minimize allocations by
        for d in degfn(g, v)    # iterating over vertices
            hist[d] = get(hist, d, 0) + 1
        end
    end
    return hist
end

###
### induced_subgraph
###
#=
function induced_subgraph(g::AdjacencyMap, elist::AbstractVector{U}) where {U <: AbstractEdge}
    h = zero(g)
    T = eltype(T)
    newvid = Dict{T,T}()
    vmap = Vector{T}()

    for e in elist
        u, v = Tuple(e)
        for i in (u, v)
            if !haskey(newvid, i)
                add_vertex!(h)
                newvid[i] = nvertices(h)
                push!(vmap, i)
            end
        end
        add_edge!(h, newvid[u], newvid[v])
    end
    return h, vmap
end
=#


text:
module Flux3D

using Flux,
    Zygote,
    NearestNeighbors,
    LinearAlgebra,
    Statistics,
    CUDA,
    FileIO,
    SparseArrays,
    Requires,
    Meshing
using Base: tail
import Zygote: @nograd, @ignore
import GeometryBasics
import Flux: @functor, functor, gpu, cpu, Chain

export gpu, cpu, Chain, visualize

# representation
include("rep/utils.jl")
include("rep/pcloud.jl")
include("rep/mesh.jl")
include("rep/voxels.jl")

# conversions
include("conversions.jl")

# transforms
include("transforms/utils.jl")
include("transforms/pcloud_func.jl")
include("transforms/mesh_func.jl")
include("transforms/transforms.jl")

# metrics
include("metrics/pcloud.jl")
include("metrics/mesh.jl")

# Dataset module
include("datasets/Dataset.jl")
using .Dataset
export CustomDataset, ModelNet10, ModelNet40

# visualization
@init @require AbstractPlotting = "537997a7-5e4e-5d89-9595-2241ea00577e" include(
    "visualize.jl",
)

# models
include("models/utils.jl")
include("models/dgcnn.jl")
include("models/pointnet.jl")

# borowed from Flux.jl
const use_cuda = Ref(false)
@init begin
    use_cuda[] = CUDA.functional() # Can be overridden after load with `Flux.use_cuda[] = false`
    if CUDA.functional()
        if !CUDA.has_cudnn()
            @warn "CUDA.jl found cuda, but did not find libcudnn. Some functionality will not be available."
        end
    end
end

end # module


text:
include("6.5.10 (c+d x)^m (a+b sech)^n.jl")
include("6.5.11 (e x)^m (a+b sech(c+d x^n))^p.jl")

max_stars_repo_path: src/rules/6 Hyperbolic functions/6.5 Hyperbolic secant/.jl


text:
using SafeTestsets
@safetestset "Potato test" begin include("testsets.jl") end


text:
module TestArrays

using Modia
using Modia.StaticArrays

@usingModiaPlot

LinearODE = Model(
    A = parameter | [-1.0  0.0;
                      0.0 -2.0],
    x = Var(init=[1.0, 2.0]),
    equations = :[der(x) = A*x]
)

linearODE = @instantiateModel(LinearODE, logCode=true)
simulate!(linearODE, stopTime = 2, log=false, logStates=false, requiredFinalStates=[0.13533533463680386, 0.036632273646086545])
plot(linearODE, ["x", "der(x)"], figure=1, heading="LinearODE with length(x)=2")

simulate!(linearODE, stopTime = 2, log=false, logStates=false, requiredFinalStates=[0.14886882881582736, 0.038462894626776434, 0.00768439894426358],
            merge = Map(A = [-1.0  0.0  0.0;
                              0.0 -2.0   0.0;
                              0.0  0.0 -3.0], x=[1.1, 2.1, 3.1]))
plot(linearODE, ["x", "der(x)"], figure=2, heading="LinearODE with length(x)=3")

LinearODE2 = Model(
    A = parameter | SMatrix{2,2}([-1.0  0.0;
                                   0.0 -2.0]),
    x = Var(init = SVector{2}(1.0, 2.0)),
    equations = :[der(x) = A*x]
)

linearODE2 = @instantiateModel(LinearODE2, logCode=true)
simulate!(linearODE2, stopTime = 2, log=false, logStates=false, requiredFinalStates=[0.13533533463680386, 0.036632273646086545])
@show typeof(linearODE2.evaluatedParameters[:A])
plot(linearODE2, ["x", "der(x)"], figure=3, heading="LinearODE2 with static arrays.")

end


text:
@deprecate read(t::AITask, precision::DataType, num_samples_per_chan::Integer = -1) read(t::AITask, num_samples_per_chan, precision)
@deprecate read(t::DITask, precision::DataType, num_samples_per_chan::Integer = -1) read(t::DITask, num_samples_per_chan, precision)

# see also `analog_input` and `analog_input_configs`


text:
using Mill, Flux
using Distributions

function sample_ds1(n = 128)
    d1 = Binomial(40, 0.5)
    d2 = Binomial(40, 0.1)
    y = rand(1:2, n)
    bs = [rand(y[i] == 1 ? d1 : d2) for i in 1:n]
    BagNode(rand(2, sum(bs)) |> ArrayNode, length2bags(bs)), Flux.onehotbatch(y, 1:2)
end

function sample_ds2(n = 128)
    d1 = Binomial(80, 0.5)
    d2 = MixtureModel(Binomial[Binomial(80, 0.95), Binomial(80, 0.1)])
    y = rand(1:2, n)
    bs = [rand(y[i] == 1 ? d1 : d2) for i in 1:n]
    BagNode(rand(2, sum(bs)) |> ArrayNode, length2bags(bs)), Flux.onehotbatch(y, 1:2)
end

ds = sample_ds2
m = reflectinmodel(ds(1)[1]; b=Dict("" => k -> Chain(Dense(k, 10, relu), Dense(10, 2, relu))))
ps = Flux.params(m)
opt = ADAM()
loss(x, y) = Flux.logitcrossentropy(m(x).data, y)
acc(x, y) = mean(Flux.onecold(m(x).data) .== Flux.onecold(y))

Flux.@epochs 10 begin
    Flux.train!(loss, ps, [ds() for _ in 1:1000], opt)
    val = ds()
    println(acc(val...))
end