text:
# Using StructArrays.jl package

using StructArrays

sa = StructArray(records);

# We can index into the records as usual

sa[1:3]
#=
julia> sa[1:3]
3-element StructArray(::Array{String,1}, ::Array{String,1}, ::Array{String,1}, ::Array{Int64,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}) with eltype TripPayment:
 TripPayment("1", "2018-12-01 00:28:22", "2018-12-01 00:44:07", 2, 2.5, 12.0, 0.5, 0.5, 3.95, 0.0, 0.3, 17.25)
 TripPayment("1", "2018-12-01 00:52:29", "2018-12-01 01:11:37", 3, 2.3, 13.0, 0.5, 0.5, 2.85, 0.0, 0.3, 17.15)
 TripPayment("2", "2018-12-01 00:12:52", "2018-12-01 00:36:23", 1, 0.0, 2.5, 0.5, 0.5, 0.0, 0.0, 0.3, 3.8)
=#

# type of an indexed record

sa[1]
typeof(sa[1])
#=
julia> sa[1]
TripPayment(1, "2018-12-01 00:28:22", "2018-12-01 00:44:07", 2, 2.5, 12.0, 0.5, 0.5, 3.95, 0.0, 0.3, 17.25)

julia> typeof(sa[1])
TripPayment
=#

# extracting one field

sa.fare_amount
#=
julia> sa.fare_amount
100000-element Array{Float64,1}:
 12.0
 13.0
  2.5
  ⋮
=#

# How does it perform?

@btime mean(sa.fare_amount);
#=
julia> @btime mean(sa.fare_amount);
  27.175 μs (1 allocation: 16 bytes)
=#

# Note that it makes a copy of the data

Base.summarysize(records) / 1024 / 1024
Base.summarysize(sa) / 1024 / 1024
#=
julia> Base.summarysize(records) / 1024 / 1024
15.068092346191406

julia> Base.summarysize(sa) / 1024 / 1024
14.305671691894531
=#

# Observation: Columnar format is more memory-efficient.

records = nothing
GC.gc()
#= remove original
julia> records = nothing

julia> GC.gc()
=#


text:
struct WordleGuess
    guess::String
    result::Vector{Symbol}

    WordleGuess(guess::String, result::Vector{Symbol}) = let
        invalid_results = setdiff(result, [CORRECT, PRESENT, INCORRECT])

        if !isempty(invalid_results)
            message = strip("""
            The following result symbols are invalid: $(join(invalid_results .|> String, ", "))
            """)

            error(message)
        end

        if length(result) != 5
            error("Results must have length 5")
        end

        new(guess, result)
    end
end

show(io::IO, guess::WordleGuess) = print(io, guess.result |> join)

ext: jl


text:
################################################################################
# MIT License
#
# Copyright (c) 2021 Hajime Nakagami<nakagami@gmail.com>
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
################################################################################

const SQL_TYPE_TEXT = 452
const SQL_TYPE_VARYING = 448
const SQL_TYPE_SHORT = 500
const SQL_TYPE_LONG = 496
const SQL_TYPE_FLOAT = 482
const SQL_TYPE_DOUBLE = 480
const SQL_TYPE_D_FLOAT = 530
const SQL_TYPE_TIMESTAMP = 510
const SQL_TYPE_BLOB = 520
const SQL_TYPE_ARRAY = 540
const SQL_TYPE_QUAD = 550
const SQL_TYPE_TIME = 560
const SQL_TYPE_DATE = 570
const SQL_TYPE_INT64 = 580
const SQL_TYPE_INT128 = 32752
const SQL_TYPE_TIMESTAMP_TZ = 32754
const SQL_TYPE_TIME_TZ = 32756
const SQL_TYPE_DEC_FIXED = 32758
const SQL_TYPE_DEC64 = 32760
const SQL_TYPE_DEC128 = 32762
const SQL_TYPE_BOOLEAN = 32764
const SQL_TYPE_NULL = 32766

xsqlvar_type_length = Dict(
    SQL_TYPE_TEXT=>-1,
    SQL_TYPE_VARYING=>-1,
    SQL_TYPE_SHORT=>4,
    SQL_TYPE_LONG=>4,
    SQL_TYPE_FLOAT=>4,
    SQL_TYPE_TIME=>4,
    SQL_TYPE_DATE=>4,
    SQL_TYPE_DOUBLE=>8,
    SQL_TYPE_TIMESTAMP=>8,
    SQL_TYPE_BLOB=>8,
    SQL_TYPE_ARRAY=>8,
    SQL_TYPE_QUAD=>8,
    SQL_TYPE_INT64=>8,
    SQL_TYPE_INT128=>16,
    SQL_TYPE_TIMESTAMP_TZ=>12,
    SQL_TYPE_TIME_TZ=>8,
    SQL_TYPE_DEC64=>8,
    SQL_TYPE_DEC128=>16,
    SQL_TYPE_DEC_FIXED=>16,
    SQL_TYPE_BOOLEAN=>1,
)

mutable struct XSQLVAR
    sqltype::Int
    sqlscale::Int
    sqlsubtype::Int
    sqllen::Int
    null_ok::Bool
    fieldname::String
    relname::String
    ownname::String
    aliasname::String
end

function io_length(x::XSQLVAR)::Int
    x.sqltype == SQL_TYPE_TEXT ? x.sqllen : xsqlvar_type_length[x.sqltype]
end

function has_precision_scale(x::XSQLVAR)::Bool
    (x.sqltype == SQL_TYPE_SHORT ||
        x.sqltype == SQL_TYPE_LONG ||
        x.sqltype == SQL_TYPE_QUAD ||
        x.sqltype == SQL_TYPE_INT64 ||
        x.sqltype == SQL_TYPE_INT128 ||
        x.sqltype == SQL_TYPE_DEC64 ||
        x.sqltype == SQL_TYPE_DEC128 ||
        x.sqltype == SQL_TYPE_DEC_FIXED) && x.sqlscale != 0
end

function _parse_date(raw_value::Vector{UInt8})::Tuple{Int, Int, Int}
    nday = bytes_to_buint32(raw_value) + 678882

    century = div(4 * nday - 1,  146097)
    nday = 4 * nday - 1 - 146097 * century
    day = nday / 4

    nday = div(4 * day + 3, 1461)
    day = 4 * day + 3 - 1461 * nday
    day = div(day + 4, 4)

    month = div(5 * day - 3, 153)
    day = 5 * day - 3 - 153 * month
    day = div(day + 5,  5)
    year = 100 * century + nday
    if month < 10
        month += 3
    else
        month -= 9
        year += 1
    end

    (year, month, day)
end

function _parse_time(raw_value::Vector{UInt8})
    n = bytes_to_buint32(raw_value)
    s = div(n, 10000)
    m = div(s, 60)
    h = div(m, 60)
    m = m % 60
    s = s % 60

    (h, m, s, (n % 10000) * 100000)
end

function parse_date(raw_value::Vector{UInt8})::Date
    year, month, day = _parse_date(raw_value)
    Date(year, month, day)
end

function parse_time(raw_value::Vector{UInt8})::Time
    h, m, s, n = _parse_time(raw_value)
    Time(h, m, s, div(n, 1000000))
end

function parse_timestamp(raw_value::Vector{UInt8})::DateTime
    year, month, day = _parse_date(raw_value[1:4])
    h, m, s, n = _parse_time(raw_value[5:8])
    DateTime(year, month, day, h, m, s, div(n, 1000000))
end

function parse_time_tz(raw_value::Vector{UInt8})::ZonedDateTime
    h, m, s, n = _parse_time(raw_value[1:4])
    timezone = TimeZones.TimeZone(get_timezone_name_by_id_dict()[bytes_to_buint16(raw_value[5:6])])
    offset = TimeZones.TimeZone(get_timezone_name_by_id_dict()[bytes_to_buint16(raw_value[7:8])])
    zdt = ZonedDateTime(0, 1, 1, h, m, s, div(n, 1000000), timezone)
    astimezone(zdt, offset)
end

function parse_timestamp_tz(raw_value::Vector{UInt8})::ZonedDateTime
    year, month, day = _parse_date(raw_value[1:4])
    h, m, s, n = _parse_time(raw_value[5:8])
    timezone = TimeZones.TimeZone(get_timezone_name_by_id_dict()[bytes_to_buint16(raw_value[9:10])])
    offset = TimeZones.TimeZone(get_timezone_name_by_id_dict()[bytes_to_buint16(raw_value[11:12])])
    zdt = ZonedDateTime(year, month, day, h, m, s, div(n, 1000000), timezone)
    astimezone(zdt, offset)
end

function value(x::XSQLVAR, raw_value::Vector{UInt8})
    if x.sqltype == SQL_TYPE_TEXT
        if x.sqlsubtype == 1
            raw_value
        else
            String(raw_value)
        end
    elseif x.sqltype == SQL_TYPE_VARYING
        if x.sqlsubtype == 1
            raw_value
        else
            String(raw_value)
        end
    elseif x.sqltype == SQL_TYPE_SHORT
        i16::Int16 = bytes_to_bint16(raw_value)
        if x.sqlscale == 0
            i16
        else
            if i16 < 0
                Decimal(1, i16*-1, x.sqlscale)
            else
                Decimal(0, i16, x.sqlscale)
            end
        end
    elseif x.sqltype == SQL_TYPE_LONG
        i32::Int32 = bytes_to_bint32(raw_value)
        if x.sqlscale == 0
            i32
        else
            if i32 < 0
                Decimal(1, i32*-1, x.sqlscale)
            else
                Decimal(0, i32, x.sqlscale)
            end
        end
    elseif x.sqltype == SQL_TYPE_FLOAT
        reinterpret(Float32, raw_value)[1]
    elseif x.sqltype == SQL_TYPE_DOUBLE
        reinterpret(Float64, raw_value)[1]
    elseif x.sqltype == SQL_TYPE_TIMESTAMP
        parse_timestamp(raw_value)
    elseif x.sqltype == SQL_TYPE_BLOB
        raw_value
    elseif x.sqltype == SQL_TYPE_TIME
        parse_time(raw_value)
    elseif x.sqltype == SQL_TYPE_DATE
        parse_date(raw_value)
    elseif x.sqltype == SQL_TYPE_INT64
        i64::Int64 = bytes_to_bint64(raw_value)
        if x.sqlscale == 0
            i64
        else
            if i64 < 0
                Decimal(1, i64*-1, x.sqlscale)
            else
                Decimal(0, i64, x.sqlscale)
            end
        end
    elseif x.sqltype == SQL_TYPE_INT128
        i128::Int128 = bytes_to_bint128(raw_value)
        if x.sqlscale == 0
            i128
        else
            if i128 < 0
                Decimal(1, i128*-1, x.sqlscale)
            else
                Decimal(0, i128, x.sqlscale)
            end
        end
    elseif x.sqltype == SQL_TYPE_TIMESTAMP_TZ
        parse_timestamp_tz(raw_value)
    elseif x.sqltype == SQL_TYPE_TIME_TZ
        parse_time_tz(raw_value)
    elseif x.sqltype == SQL_TYPE_DEC_FIXED
        decimal_fiexed_to_decimal(value_value)
    elseif x.sqltype == SQL_TYPE_DEC64
        decimal64_to_decimal(value_value)
    elseif x.sqltype == SQL_TYPE_DEC128
        decimal128_to_decimal(value_value)
    elseif x.sqltype == SQL_TYPE_BOOLEAN
        raw_value[0] != 0
    elseif x.sqltype == SQL_TYPE_NULL
        missing
    end
end

function juliatype(x::XSQLVAR)
    if x.sqltype == SQL_TYPE_TEXT
        if x.sqlsubtype == 1
            T = Vector{UInt8}
        else
            T = String
        end
    elseif x.sqltype == SQL_TYPE_VARYING
        if x.sqlsubtype == 1
            T = Vector{UInt8}
        else
            T = String
        end
    elseif x.sqltype == SQL_TYPE_SHORT
        if x.sqlscale != 0
            T = Decimal
        else
            T = Int16
        end
    elseif x.sqltype == SQL_TYPE_LONG
        T = Int64
    elseif x.sqltype == SQL_TYPE_FLOAT
        T = Float32
    elseif x.sqltype == SQL_TYPE_DOUBLE
        T = Float64
    elseif x.sqltype == SQL_TYPE_TIMESTAMP
        T = DateTime
    elseif x.sqltype == SQL_TYPE_BLOB
        if x.sqlsubtype == 1    # TEXT
            T = String
        else
            T = Vector{UInt8}
        end
    elseif x.sqltype == SQL_TYPE_TIME
        T = Time
    elseif x.sqltype == SQL_TYPE_DATE
        T = Date
    elseif x.sqltype == SQL_TYPE_INT64
        if x.sqlscale != 0
            T = Decimal
        else
            T = Int64
        end
    elseif x.sqltype == SQL_TYPE_INT128
        if x.sqlscale != 0
            T = Decimal
        else
            T = Int128
        end
    elseif x.sqltype == SQL_TYPE_TIMESTAMP_TZ
        T = ZonedDateTime
    elseif x.sqltype == SQL_TYPE_TIME_TZ
        T = ZonedDateTime
    elseif x.sqltype == SQL_TYPE_DEC_FIXED
        T = Decimal
    elseif x.sqltype == SQL_TYPE_DEC64
        T = Decimal
    elseif x.sqltype == SQL_TYPE_DEC128
        T = Decimal
    elseif x.sqltype == SQL_TYPE_BOOLEAN
        T = Bool
    elseif x.sqltype == SQL_TYPE_NULL
        T = Missing
    end

    x.null_ok ? Union{Missing, T} : T
end


text:
immutable Assignment
    lhs
    op
    args
end

immutable FunctionAssignment
    lhs
    func
    args
end

immutable GeneratedFunction
    input_arguments::Vector{Symbol}
    output_arguments::Vector{Symbol}
    code::Expr
end

# Close to single assignment form
type FlattenedAST
    input_variables::Set{Symbol}
    variables::Vector{Symbol}  # cleaned version
    intermediate::Vector{Symbol}  # generated vars
    code # ::Vector{Assignment}
end

import Base.show
function show(io::IO, flatAST::FlattenedAST)
    println(io, "input vars: ", flatAST.input_variables)
    println(io, "intermediate vars: ", flatAST.intermediate)
    println(io, "code: ")
    println(io, flatAST.code)
end

FlattenedAST() = FlattenedAST(Set{Symbol}(), [], [], [])

export FlattenedAST

export flatten!

doc"""`flatten!` adds information about any variables
generated to the `FlattenedAST` object. It returns the object
at the top of the current piece of tree."""
# process numbers
function flatten!(flatAST::FlattenedAST, ex)
    return ex  # nothing to do to the AST; return the number
end

function flatten!(flatAST::FlattenedAST, ex::Symbol)  # symbols are leaves
    push!(flatAST.input_variables, ex)  # add the discovered symbol as an input variable
    return ex
end

function flatten!(flatAST::FlattenedAST, ex::Expr)

    if ex.head == :$   # process constants of form $a
        process_constant!(flatAST, ex)

    elseif ex.head == :call
        process_call!(flatAST, ex)

    elseif ex.head == :(=)
        process_assignment!(flatAST, ex)

    elseif ex.head == :block
        process_block!(flatAST, ex)

    elseif ex.head == :tuple
        process_tuple!(flatAST, ex)

    elseif ex.head == :return
        process_return!(flatAST, ex)
    end
end

function process_constant!(flatAST::FlattenedAST, ex)
    return :(esc($(ex.args[1])))  # return the value of the constant
end

function process_block!(flatAST::FlattenedAST, ex)

    local top

    for arg in ex.args

        isa(arg, LineNumberNode) && continue

        top = flatten!(flatAST, arg)

    end

    return top  # last variable assigned
end

function process_tuple!(flatAST::FlattenedAST, ex)

    top_args = [flatten!(flatAST, arg) for arg in ex.args]

    @show "Tuple arguments", top_args

    return top_args

end

function process_assignment!(flatAST::FlattenedAST, ex)
    process_call!(flatAST, ex.args[2], ex.args[1])
end

function process_call!(flatAST::FlattenedAST, ex, new_var=nothing)
    # new_var is an optional variable name to assign the result of the call to
    # if none is given, then a new, unique variable name is created

    op = ex.args[1]

    if isa(op, Expr) && op.head == :line
        return
    end

    # rewrite +(a,b,c) as +(a,+(b,c)):
    # TODO: Use @match here!

    # if op in (:+, :*) && length(ex.args) > 3
    #     return insert_variables( :( ($op)($(ex.args[2]), ($op)($(ex.args[3:end]...) )) ))
    # end

    # new_code = quote end
    # current_args = []  # the arguments in the current expression that will be added
    # all_vars = Set{Symbol}()  # all variables contained in the sub-expressions
    # generated_variables = Symbol[]

    top_args = []
    for arg in ex.args[2:end]

        isa(arg, LineNumberNode) && continue

        top = flatten!(flatAST, arg)
        push!(top_args, top)
    end

    top_level_code = quote end

    #@show op

    if op ∈ keys(rev_ops)  # standard operator
        if new_var == nothing
            new_var = make_symbol()
        end

        push!(flatAST.intermediate, new_var)

        top_level_code = Assignment(new_var, op, top_args)

    else
        if haskey(registered_functions, op)
            println("Processing function $op")

            # make enough new variables for all the returned arguments:
            new_vars = Symbol[]

            for var in registered_functions[op].generated
                push!(new_vars, make_symbol())
            end

            append!(flatAST.intermediate, new_vars)

            top_level_code = FunctionAssignment(new_vars, op, top_args)

            new_var = new_vars[1:length(registered_functions[op].return_arguments)]

        else

            throw(ArgumentError("Function $op not supported"))
        end
    end

    push!(flatAST.code, top_level_code)

    return new_var

end

function flatten!(ex)
    flatAST = FlattenedAST()
    top_var = flatten!(flatAST, ex)

    return top_var, flatAST
end

function emit_forward_code(a::Assignment)
    :($(a.lhs) = $(a.op)($(a.args...) ) )
end

make_tuple(args) = Expr(:tuple, args...)

function emit_forward_code(a::FunctionAssignment)
    f = a.func
    args = make_tuple(a.lhs)
    :($args = $(f).forward($(a.args...) ) )
end

function emit_forward_code(code) #code::Vector{Assignment})
    new_code = quote end
    new_code.args = vcat([emit_forward_code(line) for line in code])
    return new_code
end

function emit_backward_code(a::Assignment)
    return_args = [a.lhs, a.args...]
    rev_op = rev_ops[a.op]  # find reverse operation

    rev_code = :($(rev_op)($(return_args...)))

    # delete non-symbols in return args:
    for (i, arg) in enumerate(return_args)
        if !(isa(arg, Symbol))
            return_args[i] = :_
        end
    end

    return_tuple = Expr(:tuple, return_args...)  # make tuple out of array
# or: :($(return_args...),)

    return :($(return_tuple) = $(rev_code))

end

function emit_backward_code(a::FunctionAssignment)
    f = a.func
    args = make_tuple(a.lhs)
    :($args = $(f).backward($(a.args...) ) )
end

function emit_backward_code(code) #::Vector{Assignment})
    new_code = quote end
    new_code.args = vcat([emit_backward_code(line) for line in reverse(code)])
    return new_code
end

function forward_pass(flatAST::FlattenedAST)

    @show flatAST.input_variables
    @show flatAST.intermediate

    input_variables = sort(collect(flatAST.input_variables))
    input_variables = setdiff(input_variables, flatAST.intermediate)  # remove local variables
    flatAST.variables = input_variables

    generated_code = emit_forward_code(flatAST.code)
    #make_function(input_variables, flatAST.intermediate, generated_code)
    return GeneratedFunction(input_variables, flatAST.intermediate, generated_code)
end

function backward_pass(flatAST::FlattenedAST)

    generated_code = emit_backward_code(flatAST.code)
    # make_function([flatAST.variables; flatAST.intermediate],
    #                 flatAST.variables,
    #                 generated_code)
    # # reverse input_variables and intermediate?

    all_variables = [flatAST.variables; flatAST.intermediate]
    return GeneratedFunction(all_variables,
                            flatAST.variables,
                            generated_code)

end


text:
using MPI, Test
include(joinpath("..", "..", "testhelpers.jl"))

@testset "SingleStackUtils" begin
    runmpi(joinpath(@__DIR__, "ssu_tests.jl"))
end

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
using OrdinaryDiffEq, Test, LinearAlgebra, Statistics

# create mass matrix problems
function make_mm_probs(mm_A, ::Val{iip}) where iip
  # iip
  function mm_f(du,u,p,t)
    update_coefficients!(mm_A,OrdinaryDiffEq.constvalue.(u),p,t)
    mm_b = vec(sum(mm_A; dims=2))
    mul!(du,mm_A,u)
    du .+= t * mm_b
    nothing
  end
  mm_g(du,u,p,t) = (@. du = u + t; nothing)

  # oop
  mm_f(u,p,t) = (update_coefficients!(mm_A,OrdinaryDiffEq.constvalue.(u),p,t);
                 mm_A * (u .+ t))
  mm_g(u,p,t) = u .+ t

  mm_analytic(u0, p, t) = @. 2 * u0 * exp(t) - t - 1

  u0 = ones(3)
  tspan = (0.0, 1.0)

  prob = ODEProblem(ODEFunction{iip,true}(mm_f, analytic=mm_analytic, mass_matrix=mm_A), u0, tspan)
  prob2 = ODEProblem(ODEFunction{iip,true}(mm_g, analytic=mm_analytic), u0, tspan)

  prob, prob2
end

function _norm_dsol(alg,prob,prob2,dt=1/10)
  sol = solve(prob,  alg,dt=dt,adaptive=false)
  sol2 = solve(prob2,alg,dt=dt,adaptive=false)
  norm(sol .- sol2)
end

function update_func1(A,u,p,t)
    A[1,1] = cos(t)
    A[2,1] = sin(t)*u[1]
    A[3,1] = t^2
    A[1,2] = cos(t)*sin(t)
    A[2,2] = cos(t)^2 + u[3]
    A[3,2] = sin(t)*u[2]
    A[1,3] = sin(t)
    A[2,3] = t^2
    A[3,3] = t*cos(t) + 1
end

function update_func2(A,u,p,t)
    A[1,1] = cos(t)
    A[2,1] = sin(t)
    A[3,1] = t^2
    A[1,2] = cos(t)*sin(t)
    A[2,2] = cos(t)^2
    A[3,2] = sin(t)
    A[1,3] = sin(t)
    A[2,3] = t^2
    A[3,3] = t*cos(t) + 1
end

almost_I = Matrix{Float64}(1.01I, 3, 3)
mm_A = Float64[-2 1 4; 4 -2 1; 2 1 3]
dependent_M1 = DiffEqArrayOperator(ones(3,3),update_func=update_func1)
dependent_M2 = DiffEqArrayOperator(ones(3,3),update_func=update_func2)

@testset "Mass Matrix Accuracy Tests" for mm in (almost_I, mm_A)
  # test each method for exactness
  for iip in (false, true)
    prob, prob2 = make_mm_probs(mm, Val(iip))

    println("SDIRKs")
    @test _norm_dsol(ImplicitEuler(),prob,prob2) ≈ 0 atol=1e-12
    @test _norm_dsol(Trapezoid(),prob,prob2) ≈ 0 atol=1e-12
    @test _norm_dsol(RadauIIA5(),prob,prob2) ≈ 0 atol=1e-12
    @test _norm_dsol(ImplicitMidpoint(extrapolant = :constant),prob,prob2) ≈ 0 atol=1e-10

    println("BDFs")
    @test _norm_dsol(ABDF2(),prob,prob2) ≈ 0 atol=1e-12

    @test _norm_dsol(QBDF1(),prob,prob2) ≈ 0 atol=1e-12
    @test _norm_dsol(QBDF2(),prob,prob2) ≈ 0 atol=1e-12
    @test _norm_dsol(QBDF(),prob,prob2) ≈ 0 atol=1e-12

    @test _norm_dsol(QNDF1(),prob,prob2) ≈ 0 atol=1e-12
    @test _norm_dsol(QNDF2(),prob,prob2) ≈ 0 atol=1e-12
    @test _norm_dsol(QNDF(),prob,prob2) ≈ 0 atol=1e-12

    println("Rosenbrocks")
    @test _norm_dsol(Rosenbrock23(),prob,prob2) ≈ 0 atol=1e-11
    @test _norm_dsol(Rosenbrock32(),prob,prob2) ≈ 0 atol=1e-11
    @test _norm_dsol(ROS3P(),prob,prob2) ≈ 0 atol=1e-11
    @test _norm_dsol(Rodas3(),prob,prob2) ≈ 0 atol=1e-11
    @test _norm_dsol(RosShamp4(),prob,prob2) ≈ 0 atol=1e-10
    @test _norm_dsol(Veldd4(),prob,prob2) ≈ 0 atol=1e-10
    @test _norm_dsol(Velds4(),prob,prob2) ≈ 0 atol=1e-10
    @test _norm_dsol(GRK4T(),prob,prob2) ≈ 0 atol=1e-10
    @test _norm_dsol(GRK4A(),prob,prob2) ≈ 0 atol=1e-10
    @test _norm_dsol(Ros4LStab(),prob,prob2) ≈ 0 atol=1e-10
    @test _norm_dsol(RosenbrockW6S4OS(),prob,prob2) ≈ 0 atol=1e-10
    @test _norm_dsol(ROS34PW1a(),prob,prob2) ≈ 0 atol=1e-10
    @test _norm_dsol(ROS34PW1b(),prob,prob2) ≈ 0 atol=1e-10
    @test _norm_dsol(ROS34PW2(),prob,prob2) ≈ 0 atol=1e-10
    @test _norm_dsol(ROS34PW3(),prob,prob2) ≈ 0 atol=1e-10
    @test _norm_dsol(Rodas4(),prob,prob2) ≈ 0 atol=1e-9
    @test _norm_dsol(Rodas42(),prob,prob2) ≈ 0 atol=1e-9
    @test _norm_dsol(Rodas4P(),prob,prob2) ≈ 0 atol=1e-9
    @test _norm_dsol(Rodas5(),prob,prob2) ≈ 0 atol=1e-7
  end
end

@testset "Mass Matrix Functional Iteration Tests" for mm in (almost_I,) # Requires contraction map
  for iip in (false, true)
    @show iip
    prob, prob2 = make_mm_probs(mm, Val(iip))

    sol = solve(prob,ImplicitEuler(
                          nlsolve=NLFunctional()),dt=1/10,adaptive=false,reltol=1e-7,abstol=1e-10)
    sol2 = solve(prob2,ImplicitEuler(nlsolve=NLFunctional()),dt=1/10,adaptive=false,reltol=1e-7,abstol=1e-10)
    @test norm(sol .- sol2) ≈ 0 atol=1e-7

    sol = solve(prob, ImplicitMidpoint(extrapolant = :constant,
                          nlsolve=NLFunctional()),dt=1/10,reltol=1e-7,abstol=1e-10)
    sol2 = solve(prob2,ImplicitMidpoint(extrapolant = :constant, nlsolve=NLFunctional()),dt=1/10,reltol=1e-7,abstol=1e-10)
    @test norm(sol .- sol2) ≈ 0 atol=1e-7

    sol = solve(prob,ImplicitEuler(nlsolve=NLAnderson()),dt=1/10,adaptive=false)
    sol2 = solve(prob2,ImplicitEuler(nlsolve=NLAnderson()),dt=1/10,adaptive=false)
    @test norm(sol .- sol2) ≈ 0 atol=1e-7
    @test norm(sol[end] .- sol2[end]) ≈ 0 atol=1e-7

    sol = solve(prob, ImplicitMidpoint(extrapolant = :constant, nlsolve=NLAnderson()),dt=1/10,reltol=1e-7,abstol=1e-10)
    sol2 = solve(prob2,ImplicitMidpoint(extrapolant = :constant, nlsolve=NLAnderson()),dt=1/10,reltol=1e-7,abstol=1e-10)
    @test norm(sol .- sol2) ≈ 0 atol=1e-7
    @test norm(sol[end] .- sol2[end]) ≈ 0 atol=1e-7
  end
end

# Singular mass matrices

@testset "Mass Matrix Tests with Singular Mass Matrices" begin
  function f!(du, u, p, t)
      du[1] = u[2]
      du[2] = u[2] - 1.
      return
  end

  u0 = [0.,1.]
  tspan = (0.0, 1.0)

  M = fill(0., 2,2)
  M[1,1] = 1.

  m_ode_prob = ODEProblem(ODEFunction(f!;mass_matrix=M), u0, tspan)
  @test_nowarn sol = solve(m_ode_prob, Rosenbrock23())

  M = [0.637947  0.637947
       0.637947  0.637947]

  inv(M) # not caught as singular

  function f2!(du, u, p, t)
      du[1] = u[2]
      du[2] = u[1]
      return
  end
  u0 = zeros(2)

  m_ode_prob = ODEProblem(ODEFunction(f2!;mass_matrix=M), u0, tspan)
  println("Rosenbrocks")
  sol1  = @test_nowarn solve(m_ode_prob, Rosenbrock23(), reltol=1e-10, abstol=1e-10)
  sol2  = @test_nowarn solve(m_ode_prob, RadauIIA5(), reltol=1e-10, abstol=1e-10)
  sol3  = @test_nowarn solve(m_ode_prob, Cash4(), reltol=1e-10, abstol=1e-10)
  println("SDIRKs")
  sol4  = @test_nowarn solve(m_ode_prob, ImplicitEuler(), reltol=1e-10, abstol=1e-10)
  sol5  = @test_nowarn solve(m_ode_prob, ImplicitMidpoint(), dt = 1/1000, reltol=1e-10, abstol=1e-10)
  sol6  = @test_nowarn solve(m_ode_prob, Trapezoid(), reltol=1e-10, abstol=1e-10)
  println("BDFs")
  sol7  = @test_nowarn solve(m_ode_prob, QNDF1(), reltol=1e-10, abstol=1e-10)
  sol8  = @test_nowarn solve(m_ode_prob, QBDF1(), reltol=1e-10, abstol=1e-10)
  sol9  = @test_nowarn solve(m_ode_prob, QNDF2(), reltol=1e-10, abstol=1e-10)
  sol10 = @test_nowarn solve(m_ode_prob, QBDF2(), reltol=1e-10, abstol=1e-10)
  sol11 = @test_nowarn solve(m_ode_prob, ABDF2(), reltol=1e-10, abstol=1e-10)
  sol12 = @test_nowarn solve(m_ode_prob, QBDF(), reltol=1e-10, abstol=1e-10)
  sol13 = @test_nowarn solve(m_ode_prob, QNDF(), reltol=1e-10, abstol=1e-10)

  @test sol1[end] ≈ sol2[end] atol=1e-9
  @test sol1[end] ≈ sol3[end] atol=1e-9
  @test sol1[end] ≈ sol4[end] atol=1e-9
  @test sol1[end] ≈ sol5[end] atol=1e-9
  @test sol1[end] ≈ sol6[end] atol=1e-9
  @test sol1[end] ≈ sol7[end] atol=1e-9
  @test sol1[end] ≈ sol8[end] atol=1e-9
  @test sol1[end] ≈ sol9[end] atol=1e-9
  @test sol1[end] ≈ sol10[end] atol=1e-9
  @test sol1[end] ≈ sol11[end] atol=1e-9
  @test sol1[end] ≈ sol12[end] atol=1e-9
  @test sol1[end] ≈ sol13[end] atol=1e-9
end

function _norm_dsol2(alg,prob,prob2; kwargs...)
  sol = solve(prob,  alg; kwargs...)
  sol2 = solve(prob2,alg; kwargs...)
  norm(sol[end] .- sol2[end])
end
@testset "Dependent Mass Matrix Tests" for mm in (dependent_M1, dependent_M2)
  # test each method for exactness
  for iip in (false, true)
    @show iip
    prob, prob2 = make_mm_probs(mm, Val(iip))
    eulersol = solve(prob, ImplicitEuler(nlsolve=NLNewton(κ=1e-10)), reltol=1e-10)
    @test _norm_dsol2(ImplicitEuler(nlsolve=NLNewton(κ=1e-10)),prob,prob2,reltol=1e-10) ≈ 0 atol=5e-4
    @test _norm_dsol2(ImplicitMidpoint(nlsolve=NLNewton(κ=1e-10)),prob,prob2,tstops=eulersol.t) ≈ 0 atol=1e-6
    @test_skip _norm_dsol(RadauIIA5(),prob,prob2) ≈ 0 atol=1e-12

    @test_skip _norm_dsol(QNDF1(),prob,prob2) ≈ 0 atol=1e-7
    @test_skip _norm_dsol(QBDF1(),prob,prob2) ≈ 0 atol=1e-7
    @test_skip _norm_dsol(QNDF2(),prob,prob2) ≈ 0 atol=1e-7
    @test_skip _norm_dsol(QBDF2(),prob,prob2) ≈ 0 atol=1e-7
    @test_skip _norm_dsol(ABDF2(),prob,prob2) ≈ 0 atol=1e-7
    @test_skip _norm_dsol(QBDF(),prob,prob2) ≈ 0 atol=1e-7
    @test_skip _norm_dsol(QNDF(),prob,prob2) ≈ 0 atol=1e-7
  end

  # test functional iteration
  @test_skip for iip in (false, true)
    @show iip
    prob, prob2 = make_mm_probs(mm, Val(iip))

    sol = solve(prob,ImplicitEuler(
                          nlsolve=NLFunctional()),dt=1/10,adaptive=false,reltol=1e-7,abstol=1e-10)
    sol2 = solve(prob2,ImplicitEuler(nlsolve=NLFunctional()),dt=1/10,adaptive=false,reltol=1e-7,abstol=1e-10)
    @test_broken norm(sol .- sol2) ≈ 0 atol=1e-7

    sol = solve(prob, ImplicitMidpoint(extrapolant = :constant,
                          nlsolve=NLFunctional()),dt=1/10,reltol=1e-7,abstol=1e-10)
    sol2 = solve(prob2,ImplicitMidpoint(extrapolant = :constant, nlsolve=NLFunctional()),dt=1/10,reltol=1e-7,abstol=1e-10)
    @test_broken norm(sol .- sol2) ≈ 0 atol=1e-7

    sol = solve(prob,ImplicitEuler(nlsolve=NLAnderson()),dt=1/10,adaptive=false)
    sol2 = solve(prob2,ImplicitEuler(nlsolve=NLAnderson()),dt=1/10,adaptive=false)
    @test norm(sol .- sol2) ≈ 0 atol=1e-7
    @test norm(sol[end] .- sol2[end]) ≈ 0 atol=1e-7

    sol = solve(prob, ImplicitMidpoint(extrapolant = :constant, nlsolve=NLAnderson()),dt=1/10,reltol=1e-7,abstol=1e-10)
    sol2 = solve(prob2,ImplicitMidpoint(extrapolant = :constant, nlsolve=NLAnderson()),dt=1/10,reltol=1e-7,abstol=1e-10)
    @test norm(sol .- sol2) ≈ 0 atol=1e-7
    @test norm(sol[end] .- sol2[end]) ≈ 0 atol=1e-7
  end
end


text:
module CovariantMatrix

using ProgressMeter
using JuLIP
using JuLIP: sites
using ACE
using ACE: BondEnvelope, cutoff_env, cutoff_radialbasis, ACEBasis
using LinearAlgebra: norm, dot
using StaticArrays
using NeighbourLists
using LinearAlgebra
using ACEatoms

#import ChainRules
#import ChainRules: rrule, NoTangent
using Zygote
using ACEds.Utils: toMatrix

export CovSpeciesMatrixBasis, CovMatrixBasis, MatrixModel, evaluate, evaluate!, Sigma, Gamma

abstract type MatrixModel end

function evaluate(m::MatrixModel, at::AbstractAtoms; nlist=nothing)
    B = allocate_B(m, length(at))
    evaluate!(B, m, at; nlist=nlist)
    return B
end

allocate_B(basis::MatrixModel, n_atoms::Int) = [zeros(SVector{3,Float64},n_atoms,n_atoms) for n=1:length(basis)]

struct CovMatrixBasis <: MatrixModel
    onsite_basis
    offsite_basis
    inds::Dict{AtomicNumber, UnitRange{Int}}
    r_cut::Real
    offsite_env::BondEnvelope
end

cutoff(basis::CovMatrixBasis) =  maximum([cutoff_onsite(basis), cutoff_offsite(basis)])
cutoff_onsite(basis::CovMatrixBasis) = basis.r_cut
cutoff_offsite(basis::CovMatrixBasis) = cutoff_env(basis.offsite_env)
Base.length(basis::CovMatrixBasis) = length(basis.onsite_basis) + length(basis.offsite_basis)

CovMatrixBasis(onsite, offsite, r_cut::Real, offsite_env::BondEnvelope) = CovMatrixBasis(onsite, offsite, _get_basisinds(onsite,offsite), r_cut, offsite_env)

_get_basisinds(m::CovMatrixBasis) = _get_basisinds(m.onsite_basis,m.offsite_basis)

function _get_basisinds(onsite_basis,offsite_basis)
    inds = Dict{Bool, UnitRange{Int}}()
    inds[true] = 1:length(onsite_basis)
    inds[false] = (length(onsite_basis)+1):(length(onsite_basis)+length(offsite_basis))
    return inds
end

function evaluate!(B::AbstractVector{Matrix{SVector{3, Float64}}}, m::CovMatrixBasis, at::AbstractAtoms; nlist=nothing)
    if nlist === nothing
        nlist = neighbourlist(at, cutoff(m))
    end
    for k=1:length(at)
        evaluate_onsite!(B[m.inds[true]], m, at, k, nlist)
        evaluate_offsite!(B[m.inds[false]],m, at, k, nlist)
    end
end

function evaluate_onsite!(B::AbstractVector{Matrix{SVector{3, Float64}}}, m::CovMatrixBasis, at::AbstractAtoms, k::Int, nlist::PairList)
    Js, Rs = NeighbourLists.neigs(nlist, k)
    Zs = at.Z[Js]
    onsite_cfg = [ ACE.State(rr = r, mu = z)  for (r,z) in zip( Rs,Zs) if norm(r) <= m.r_cut] |> ACEConfig
    B_vals = ACE.evaluate(m.onsite_basis, onsite_cfg) # can be improved by pre-allocating memory
    for (b,b_vals) in zip(B,B_vals)
        b[k,k] += real(b_vals.val)
    end
end

function evaluate_offsite!(B::AbstractVector{Matrix{SVector{3, Float64}}}, m::CovMatrixBasis, at::AbstractAtoms, k::Int, nlist::PairList)
    Js, Rs = NeighbourLists.neigs(nlist, k)
    Zs = at.Z[Js]
    bondatoms =  [(j = j,r=r,z=z) for (j,r,z) in zip(Js,Rs,Zs ) if norm(r)<= m.offsite_env.r0cut] # atoms within max bond length
    for ba in bondatoms
        config = [ ACE.State(rr = r, rr0 = ba.r, be = (j==ba.j ? :bond : :env ), mu = z)  for (j,r,z) in zip(Js, Rs,Zs)]
        bond_config = [c for c in config if filter(m.offsite_env, c)] |> ACEConfig
        B_vals = ACE.evaluate(m.offsite_basis, bond_config) # can be improved by pre-allocating memory
        for (b,b_vals) in zip(B,B_vals)
            b[ba.j,k]+= real(b_vals.val)
        end
    end
end

struct CovSpeciesMatrixBasis <: MatrixModel
    models::Dict{AtomicNumber, CovMatrixBasis}  # model = basis
    inds::Dict{Tuple{AtomicNumber,Bool}, UnitRange{Int}}
end

CovSpeciesMatrixBasis(models::Dict{AtomicNumber, CovMatrixBasis}) = CovSpeciesMatrixBasis(models, _get_basisinds(models))

cutoff(basis::CovSpeciesMatrixBasis) = maximum(cutoff, values(basis.models))
cutoff_onsite(basis::CovSpeciesMatrixBasis) = maximum(cutoff_onsite, values(basis.models))
cutoff_offsite(basis::CovSpeciesMatrixBasis) = maximum(cutoff_offsite, values(basis.models))

Base.length(basis::CovSpeciesMatrixBasis) = sum(length, values(basis.models))

function _get_basisinds(models::Dict{AtomicNumber, CovMatrixBasis})
    inds = Dict{Tuple{AtomicNumber,Bool}, UnitRange{Int}}()
    i0 = 0
    onsite = true
    for (z, mo) in models
        len = length(mo)
        for onsite in [true,false]
            inds[(z,onsite)] = i0 .+ _get_basisinds(mo)[onsite]
        end
        i0 += len
    end
    return inds
 end

function evaluate!(B::AbstractVector{Matrix{SVector{3, Float64}}}, m::CovSpeciesMatrixBasis, at::AbstractAtoms; nlist=nothing)
    if nlist === nothing
        nlist = neighbourlist(at, cutoff(m))
    end
    for k=1:length(at)
        z0 = at.Z[k]
        evaluate_onsite!(B[m.inds[(z0,true)]], m.models[z0], at, k, nlist)
        evaluate_offsite!(B[m.inds[(z0,false)]], m.models[z0], at, k, nlist)
    end
end

Sigma(params::Vector{T}, basis) where {T <: Real} = sum( p * b for (p, b) in zip(params, basis) )

function Sigma(params_vec::Vector{Vector{T}}, basis) where {T <: Real}
    return hcat( [Sigma(params, basis) for params in params_vec]... )
end

function Sigma(params, basis, n_rep)
    n_basis = length(basis)
    @assert length(params) == n_rep * n_basis
    return hcat( [Sigma(params[((i-1)*n_basis+1):(i*n_basis)], basis) for i=1:n_rep]... )
end

MSMatrix{T} = Union{SMatrix{N,M,T}, Matrix{T}} where {N, M, T}
function Gamma(params, B::AbstractVector{MSMatrix{T}}; n_rep = 1) where {T<: Real}
    S = Sigma(params, B, n_rep)
    return outer(S,S)
end

#

function outer( A::Matrix{SVector{3, Float64}}, B::Matrix{SVector{3, Float64}})
    return sum([ kron(a, b') for (i,a) in enumerate(ac), (j,b) in enumerate(bc)]  for  (ac,bc) in zip(eachcol(A),eachcol(B)))
end
function outer( A::MSMatrix{T}, B::MSMatrix{T}) where {T}
    return A * B'
end

function get_dataset(model, raw_data; inds = nothing)
    return @showprogress [
        begin
            B = evaluate(model,at)
            if inds ===nothing
                inds = 1:length(at)
            end
            B_dense = toMatrix.(map(b -> b[inds,inds],B))
            (B = B_dense,Γ=toMatrix(Γ))
        end
        for (at,Γ) in raw_data ]
end

#=
function rrule(::typeof(Sigma), params, basis)
    val = Sigma(params,basis)

    function pb(dW)
       @assert dW isa Matrix
       # grad_params( <dW, val> )
       #grad_params = [dot(dW, b) for b in m.basis]
       grad_params = Zygote.gradient(p -> dot(dW, Sigma(p, basis)), params)[1]
       return NoTangent(), grad_params
    end

    return val, pb
 end
=#

#=
Sigma(params, basis) = sum( p * b for (p, b) in zip(params, basis) )

function rrule(::typeof(Sigma), params, basis)
    val = Sigma(params,basis)

    function pb(dW)
       @assert dW isa Matrix
       # grad_params( <dW, val> )
       #grad_params = [dot(dW, b) for b in m.basis]
       grad_params = Zygote.gradient(p -> dot(dW, Sigma(p, basis)), m.params)[1]
       return NoTangent(), grad_params
    end

    return val, pb
 end

function Gamma(params::AbstractVector{Float64}, B::AbstractVector{Matrix{T}}) where {T}
    return gamma_inner(params,B)
end

function gamma_inner(params,B)
    S = Sigma(params, B)
    return outer(S,S)
end

function outer( A::Matrix{SVector{3, Float64}}, B::Matrix{SVector{3, Float64}})
    return sum([ kron(a, b') for (i,a) in enumerate(ac), (j,b) in enumerate(bc)]  for  (ac,bc) in zip(eachcol(A),eachcol(B)))
end
function outer( A::Matrix{Float64}, B::Matrix{Float64})
    return A * B'
end
=#

#=
function ChainRules.rrule(::typeof(Sigma), params::AbstractVector{Float64}, B::AbstractVector{Matrix{Float64}})
   val = Sigma(params,B)

   function pb(dW)
      @assert dW isa Matrix
      #grad_params = [dot(dW, b) for b in B]
      #grad_params = [sum(dW.* b) for b in m.basis]
      grad_params = Zygote.gradient(p -> dot(dW, sigma_inner(p, B)), params)[1]
      return NoTangent(), grad_params
   end

   return val, pb
end

function ChainRules.rrule(::typeof(Gamma), params::AbstractVector{Float64}, B::AbstractVector{Matrix{Float64}})
    val = Gamma(params,B)

    function pb(dW)
       @assert dW isa Matrix
       #grad_params = [dot(dW, b) for b in B]
       #grad_params = [sum(dW.* b) for b in m.basis]
       grad_params = Zygote.gradient(p -> dot(dW, gamma_inner(p, B)), params)[1]
       return NoTangent(), grad_params
    end

    return val, pb
 end
=#

end

max_stars_repo_name: ACEsuit/ACEds

max_issues_repo_name: ACEsuit/ACEds

max_forks_repo_name: ACEsuit/ACEds


text:
using IterTools, SparseArrays, Test

import Base: IteratorSize, IteratorEltype
import Base: IsInfinite, SizeUnknown, HasLength, HasShape, HasEltype, EltypeUnknown

import Base.Iterators: take, countfrom, drop, peek

include("testing_macros.jl")

@testset "IterTools" begin
@testset "iterators" begin
    @testset "firstrest" begin
        # Ranges/generators have different rest states vs array/tuples
        test_base_cases = [
            (1:1, 1),
            (1:3, 1),
            ([1], 2),
            ([1, 2, 3], 2),
            ((1,), 2),
            ((1, 2, 3), 2),
            ((i for i in 1:1), 1),
            ((i for i in 1:3), 1),
        ]
        @testset "$xs" for (xs, s) in test_base_cases
            f, r = firstrest(xs)
            @test f == first(xs)
            @test collect(r) == collect(Iterators.rest(xs, s))
        end

        test_empty_cases = [
            (1:0, 1),
            (Int[], 2),
            ((), 2),
            ((i for i in 1:0), 1),
        ]

        @testset "$xs" for (xs, s) in test_empty_cases
            @test_throws ArgumentError firstrest(xs)
        end
    end

    @testset "takestrict" begin
        itr = 1:10
        take_itr = takestrict(itr, 5)
        @test eltype(take_itr) == Int
        @test IteratorEltype(take_itr) isa HasEltype
        @test IteratorSize(take_itr) isa HasLength
        @test length(take_itr) == 5
        @test collect(take_itr) == collect(1:5)

    end

    @testset "ncycle" begin
        ncy1 = ncycle(0:3,3)

        i = 0
        for j in ncy1
            @test j == i % 4
            i += 1
        end

        @test eltype(ncy1) == Int
        i = 0
        for j in collect(ncy1)
            @test j == i % 4
            i += 1
        end
    end

    @testset "repeatedly" begin
        i = 0
        for j = repeatedly(() -> 1, 10)
            @test j == 1
            i += 1
        end
        @test i == 10
        for j = repeatedly(() -> 1)
            @test j == 1
            i += 1
            i <= 10 || break
        end
    end

    @testset "distinct" begin
        x = [5, 2, 2, 1, 2, 1, 1, 2, 4, 2]
        unique_x = unique(x)
        di0 = distinct(x)
        @test eltype(di0) == Int
        @test collect(di0) == unique_x

        # repeated operations on the same Distinct iterator should function identically
        @test collect(di0) == unique_x
    end

    @testset "partition" begin
        pa0 = partition(take(countfrom(1), 6), 2)
        @test eltype(pa0) == Tuple{Int, Int}
        @test collect(pa0) == [(1,2), (3,4), (5,6)]

        pa1 = partition(take(countfrom(1), 4), 2, 1)
        @test eltype(pa1) == Tuple{Int, Int}
        @test collect(pa1) == [(1,2), (2,3), (3,4)]

        pa2 = partition(take(countfrom(1), 8), 2, 3)
        @test eltype(pa2) == Tuple{Int, Int}
        @test collect(pa2) == [(1,2), (4,5), (7,8)]

        pa3 = partition(take(countfrom(1), 0), 2, 1)
        @test eltype(pa3) == Tuple{Int, Int}
        @test collect(pa3) == []

        pa4 = partition(1:8, 1)
        @test eltype(pa4) == Tuple{Int}
        @test collect(pa4) == [(1,), (2,), (3,), (4,), (5,), (6,), (7,), (8,)]

        @test_throws ArgumentError partition(take(countfrom(1), 8), 2, 0)
    end

    @testset "imap" begin
        function test_imap(expected, input...)
            result = collect(imap(+, input...))
            @test result == expected
        end

        @testset "empty arrays" begin
            test_imap(
                Any[],
                []
            )

            test_imap(
                Any[],
                Union{}[]
            )
        end

        @testset "simple operation" begin
            test_imap(
                Any[1,2,3],
                [1,2,3]
            )
        end

        @testset "multiple arguments" begin
            test_imap(
                Any[5,7,9],
                [1,2,3],
                [4,5,6]
            )
        end

        @testset "different-length arguments" begin
            test_imap(
                Any[2,4,6],
                [1,2,3],
                countfrom(1)
            )
        end
    end

    @testset "groupby" begin
        function test_groupby(input, expected)
            result = collect(groupby(x -> x[1], input))
            @test result == expected
        end

        @testset "empty arrays" begin
            test_groupby(
                [],
                Any[]
            )

            test_groupby(
                Union{}[],
                Any[]
            )
        end

        @testset "singletons" begin
            test_groupby(
                ["xxx"],
                Any[["xxx"]]
            )
        end

        @testset "typical operation" begin
            test_groupby(
                ["face", "foo", "bar", "book", "baz"],
                Any[["face", "foo"], ["bar", "book", "baz"]]
            )
        end

        @testset "trailing singletons" begin
            test_groupby(
                ["face", "foo", "bar", "book", "baz", "xxx"],
                Any[["face", "foo"], ["bar", "book", "baz"], ["xxx"]]
            )
        end

        @testset "leading singletons" begin
            test_groupby(
                ["xxx", "face", "foo", "bar", "book", "baz"],
                Any[["xxx"], ["face", "foo"], ["bar", "book", "baz"]]
            )
        end

        @testset "middle singletons" begin
            test_groupby(
                ["face", "foo", "xxx", "bar", "book", "baz"],
                Any[["face", "foo"], ["xxx"], ["bar", "book", "baz"]]
            )
        end
    end

    @testset "subsets" begin
        @testset "all lengths" begin
            s0 = subsets(Any[])
            @test eltype(eltype(s0)) == Any
            @test collect(s0) == Vector{Any}[Any[]]

            s1 = subsets([:a])
            @test eltype(eltype(s1)) == Symbol
            @test collect(s1) == Vector{Symbol}[Symbol[], Symbol[:a]]

            s2 = subsets([:a, :b, :c])
            @test eltype(eltype(s2)) == Symbol
            @test collect(s2) == Vector{Symbol}[
                Symbol[], Symbol[:a], Symbol[:b], Symbol[:a, :b], Symbol[:c],
                Symbol[:a, :c], Symbol[:b, :c], Symbol[:a, :b, :c],
            ]
        end

        @testset "specific length" begin
            sk0 = subsets(Any[],0)
            @test eltype(eltype(sk0)) == Any
            @test collect(sk0) == Vector{Any}[Any[]]

            sk1 = subsets([:a, :b, :c], 1)
            @test eltype(eltype(sk1)) == Symbol
            @test collect(sk1) == Vector{Symbol}[Symbol[:a], Symbol[:b], Symbol[:c]]

            sk2 = subsets([:a, :b, :c], 2)
            @test eltype(eltype(sk2)) == Symbol
            @test collect(sk2) == Vector{Symbol}[Symbol[:a,:b], Symbol[:a,:c], Symbol[:b,:c]]

            sk3 = subsets([:a, :b, :c], 3)
            @test eltype(eltype(sk3)) == Symbol
            @test collect(sk3) == Vector{Symbol}[Symbol[:a,:b,:c]]

            sk4 = subsets([:a, :b, :c], 4)
            @test eltype(eltype(sk4)) == Symbol
            @test collect(sk4) == Vector{Symbol}[]

            @testset for i in 1:5
                sk5 = subsets(collect(1:4), i)
                @test eltype(eltype(sk5)) == Int
                @test length(collect(sk5)) == binomial(4, i)
            end

            @testset "implicit conversions" begin
                sk10 = subsets(1:4, 3)
                @test eltype(eltype(sk10)) == Int
                @test length(collect(sk10)) == binomial(4, 3)

                sk11 = subsets(1:3, Int32(2))
                @test eltype(eltype(sk11)) == Int
                @test length(collect(sk11)) == binomial(3, 2)
            end
        end

        @testset "specific static length" begin
            sk0 = subsets([:a, :b, :c], Val{0}())
            @test collect(sk0) == [()]

            sk1 = subsets([:a, :b, :c], Val{1}())
            @test eltype(eltype(sk1)) == Symbol
            @test collect(sk1) == [(:a,), (:b,), (:c,)]

            sk2 = subsets([:a, :b, :c], Val{2}())
            @test eltype(eltype(sk2)) == Symbol
            @test collect(sk2) == [(:a, :b), (:a, :c), (:b, :c)]

            sk3 = subsets([:a, :b, :c], Val{3}())
            @test eltype(eltype(sk3)) == Symbol
            @test collect(sk3) == [(:a, :b, :c)]

            sk4 = subsets([:a, :b, :c], Val{4}())
            @test eltype(eltype(sk4)) == Symbol
            @test collect(sk4) == []

            sk5 = subsets([:a, :b, :c], Val{5}())
            @test eltype(eltype(sk5)) == Symbol
            @test collect(sk5) == []

            @testset for i in 1:6
                sk5 = subsets(collect(1:4), Val{i}())
                @test eltype(eltype(sk5)) == Int
                @test length(collect(sk5)) == binomial(4, i)
            end

            function collect_pairs(x)
                p = Vector{NTuple{2, eltype(x)}}(undef, binomial(length(x), 2))
                idx = 1
                for i = 1:length(x)
                    for j = i+1:length(x)
                        p[idx] = (x[i], x[j])
                        idx += 1
                    end
                end
                return p
            end
            @testset for n = 1:10
                @test collect(subsets(1:n, Val{2}())) == collect_pairs(1:n)
            end
        end
    end

    @testset "nth" begin
        @testset for xs in Any[
            [1, 2, 3],
            1:3,
            reshape(1:3, 3, 1),
            sparse(reshape(1:3, 1, 3)),
        ]
            @test nth(xs, 3) == 3
            @test_throws BoundsError nth(xs, 0)
            @test_throws BoundsError nth(xs, 4)
        end

        @testset for xs in Any[take(1:3, 3), drop(-1:3, 2)]
            @test nth(xs, 3) == 3
            @test_throws BoundsError nth(xs, 0)
        end

        s = subsets([1, 2, 3])
        @test_throws BoundsError nth(s, 0)
        @test_throws BoundsError nth(s, length(s) + 1)

        # #100
        @test nth(drop(repeatedly(() -> 1), 1), 1) == 1
    end

    @testset "takenth" begin
        tn0 = takenth(Any[], 10)
        @test eltype(tn0) == Any
        @test collect(tn0) == Any[]

        tn1 = takenth(Int[], 10)
        @test eltype(tn1) == Int
        @test collect(tn1) == Int[]

        @test_throws ArgumentError takenth([], 0)

        tn2 = takenth(10:20, 3)
        @test eltype(tn2) == Int
        @test collect(tn2) == [12,15,18]

        tn3 = takenth(10:20, 1)
        @test eltype(tn3) == Int
        @test collect(tn3) == collect(10:20)
    end

    @testset "iterated" begin
        times_called = Ref(0)

        function iter_func(x)
            times_called[] += 1
            x + 1
        end

        itd = iterated(iter_func, 3)
        @test IteratorSize(itd) isa IsInfinite
        @test collect(take(itd, 4)) == 3:6
        # the first item is just the seed so iter_func isn't called
        @test times_called[] == 3
    end

    @testset "peekiter" begin
        pi0 = peekiter(1:10)
        @test IteratorEltype(pi0) isa HasEltype
        @test eltype(pi0) == Int
        @test collect(pi0) == collect(1:10)

        pi1 = peekiter([])
        @test IteratorEltype(pi1) isa HasEltype
        @test eltype(pi1) == eltype([])
        @test collect(pi1) == collect([])

        it = peekiter([:a, :b, :c])
        @test IteratorEltype(it) isa HasEltype
        @test eltype(it) == Symbol
        x, s = iterate(it)
        @test x == :a
        @test peek(it, s) == Some(:b)

        @test iterate(peekiter([])) === nothing

        it = peekiter(1:10)
        x, s = iterate(it)
        @test x == 1
        @test peek(it, s) == Some(2)

        it = peekiter(1:1)
        x, s = iterate(it)
        @test x == 1
        @test peek(it, s) === nothing
        @test iterate(it, s) === nothing
    end

    @testset "ivec" begin
        irange = 1:12
        vector = collect(irange)
        @test collect(ivec(irange)) == vector
        @test collect(ivec(vector)) == vector

        matrix = reshape(vector, 3, 4)
        @test collect(ivec(matrix)) == vector

        ndarray = reshape(vector, 2, 2, 3)
        @test collect(ivec(ndarray)) == vector
    end

    @testset "flagfirst" begin
        v = rand(1:10, 20)
        Tv = typeof(v)
        ff = flagfirst(v)
        Tff = typeof(ff)
        @test IteratorEltype(Tff) ≡ IteratorEltype(v)
        @test eltype(Tff) ≡ Tuple{Bool, eltype(v)}
        @test collect(flagfirst(v)) ==
            collect(zip(vcat([true], fill(false, length(v) - 1)), v))

        @test collect(flagfirst(Int[])) == Tuple{Bool,Int}[]
    end

    @testset "takewhile" begin
        @test collect(takewhile(x -> x^2 < 10, 1:10)) == Any[1, 2, 3]
        @test collect(takewhile(x -> x^2 < 10, Iterators.countfrom(1))) == Any[1, 2, 3]
        @test collect(takewhile(x -> x^2 < 10, 5:10)) == Any[]
        @test collect(takewhile(x -> true, 5:10)) == collect(5:10)
    end

    @testset "properties" begin
        p1 = properties(1 + 2im)
        @test IteratorEltype(p1) isa HasEltype
        @test eltype(p1) == Any
        @test IteratorSize(p1) isa HasLength
        @test length(p1) == 2
        @test collect(p1) == Any[(:re, 1), (:im, 2)]

        ntp = (a = "", b = 1, c = 2.0)
        p2 = properties(ntp)
        @test collect(p2) == Tuple.(collect(pairs(ntp)))

         # HasLength used as an example no-field struct
        p3 = properties(HasLength())
        @test collect(p3) == Any[]
    end

    @testset "propertyvalues" begin
        pv1 = propertyvalues(1 + 2im)
        @test IteratorEltype(pv1) isa HasEltype
        @test eltype(pv1) == Any
        @test IteratorSize(pv1) isa HasLength
        @test length(pv1) == 2
        @test collect(pv1) == Any[1, 2]

        tp = ("", 1, 2.0)
        pv2 = propertyvalues(tp)

        # getproperty for tuples wasn't introduced until 1.2
        # https://github.com/JuliaLang/julia/pull/31324
        @static if VERSION < v"1.2.0-DEV.460"
            @test_broken collect(pv2) == collect(tp)
        else
            @test collect(pv2) == collect(tp)
        end

        # HasLength used as an example no-field struct
        pv3 = propertyvalues(HasLength())
        @test collect(pv3) == Any[]
    end

    @testset "fieldvalues" begin
        fv1 = fieldvalues(1 + 2im)
        @test IteratorEltype(fv1) isa HasEltype
        @test eltype(fv1) == Any
        @test IteratorSize(fv1) isa HasLength
        @test length(fv1) == 2
        @test collect(fv1) == Any[1, 2]

        tp = ("", 1, 2.0)
        fv2 = fieldvalues(tp)
        @test collect(fv2) == collect(tp)

        # HasLength used as an example no-field struct
        fv3 = fieldvalues(HasLength())
        @test collect(fv3) == Any[]
    end

    @testset "traits overriding defaults" begin
        iters = [
            firstrest(1:10),
            takestrict(1:10, 5),
            repeatedly(() -> 1, 10),
            distinct([5, 2, 2, 1, 2, 1, 1, 2, 4, 2]),
            partition(take(countfrom(1), 6), 2),
            groupby(x -> x[1], [(1,2),(3,4),(1,4)]),
            imap(+, [1,2,3], [3, 2, 1]),
            subsets([:a, :b, :c]),
            iterated(sin, 3),
            nth(1:3, 2),
            takenth(1:10, 2),
            peekiter(1:10),
            ncycle(0:3,3),
            ivec(ones(3,3)),
            flagfirst(1:10),
            takewhile(x -> x^2 < 10, 1:10),
            properties(1 + 2im),
            propertyvalues(1 + 2im),
            fieldvalues(1 + 2im)
        ]
        for iter in iters
            @test IteratorSize(iter) == IteratorSize(typeof(iter))
            # indirect way to test the same thing, through generators
            g = (x for x in iter)
            @test IteratorSize(g) == IteratorSize(iter)
            @test IteratorSize(typeof(g)) == IteratorSize(typeof(iter))
        end
    end
end
end


text:
module AccuracyAtTopDual

using LIBSVM
using ProgressMeter
using Random
using Roots
using UnPack
using Statistics

import MLKernels
import KernelFunctions

using Base.Iterators: partition

abstract type Surrogate end
abstract type Model end
abstract type AbstractTopPush{S<:Surrogate} <: Model end
abstract type AbstractPatMat{S<:Surrogate} <: Model end

function model_name(model::M) where {M<:Model}
    params = values(parameters(model))
    return string(M.name.name, "(", join(params, ", "), ")")
end
Base.show(io::IO, model::Model) = print(io, model_name(model))

export Hinge, Quadratic, KernelType, Linear, Gaussian
export Model, TopPush, TopPushK, τFPL, SVM, PatMat, PatMatNP
export solve!, predict

include("projections.jl")
include("kernels.jl")
include("utilities.jl")
include("toppush.jl")
include("patmat.jl")
include("svm.jl")
include("solve.jl")

end # module


text:
using JuLIPMaterials
using Test

tests = [
    "test_greensfunctions.jl",
    "test_gfcorrectors3D.jl",
    "test_dislocations.jl",
    "test_forceconstants.jl"
    # "test_cauchyborn1.jl",
]

println("≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡")
println("   JuLIPMaterials: Start Tests   ")
println("≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡")
@testset "JuLIPMaterials" begin
   for t in tests
      ts = @testset "$(t[6:end-3])" begin
         include(t)
      end
      println(ts)
   end
end


text:
Model_NCS = HA_solver.NCS_growth_params()

Ev, grid = HA_solver.Solve_Smolyak(Model_NCS)
@testset "Solve_Smolyak" begin

  #test for NCS_growth_params, check if params make sense
  Model_NCS.β = 2.0
  @test_throws AssertionError HA_solver.Solve_Smolyak(Model_NCS)
  Model_NCS.β = 0.984

  Model_NCS.σ = 0.0
  @test_throws AssertionError HA_solver.Solve_Smolyak(Model_NCS)
  Model_NCS.σ = 2.0

  Model_NCS.rho = -1.0
  @test_throws AssertionError HA_solver.Solve_Smolyak(Model_NCS)
  Model_NCS.rho = 0.979

  Model_NCS.α = -0.5
  @test_throws AssertionError HA_solver.Solve_Smolyak(Model_NCS)
  Model_NCS.α = 0.33

  Model_NCS.δ = -0.08
  @test_throws AssertionError HA_solver.Solve_Smolyak(Model_NCS)
  Model_NCS.δ = 0.025

  #check if the model solves, if it does not solve (solution not converge), that will throw error

end

@testset "simulate_smolyak_NGS" begin
  #test for NCS_growth_params, check if params make sense

  Model_NCS.sigma_e = -10.0
  @test_throws AssertionError HA_solver.simulate_smolyak_NGS(Model_NCS, Ev)
  Model_NCS.sigma_e = 0.0072
  Model_NCS.rho = -1.0
  @test_throws AssertionError HA_solver.simulate_smolyak_NGS(Model_NCS, Ev)
  Model_NCS.rho = 0.979
  Model_NCS.α = -0.5
  @test_throws AssertionError HA_solver.simulate_smolyak_NGS(Model_NCS, Ev)
  Model_NCS.α = 0.33

  Model_NCS.δ = -0.08
  @test_throws AssertionError HA_solver.simulate_smolyak_NGS(Model_NCS, Ev)
  Model_NCS.δ = 0.025
  #check if the model solves, if it does not solve (solution not converge), that will throw error
  sim_result = HA_solver.simulate_smolyak_NGS(Model_NCS, Ev)
  #check if the simulation results have sense
  @test abs(sim_result[1] - Model_NCS.rho)<=0.06
  @test abs(sim_result[2] - Model_NCS.sigma_e)<=0.002
  #interest rate for this model should be small
  @test sim_result[3]-1.0<=0.05
  #if no error if thrown than finish tests
  @test 1>0

end


text:
using Test, SafeTestsets

# @time @safetestset "Data simulation test" begin include("simulation_test.jl") end
# @time @safetestset "RTrees Test" begin include("RTrees_test.jl") end
# @time @safetestset "XMLConstructor Test" begin include("xml_test.jl") end
# @time @safetestset "General workflow test" begin include("workflow_test.jl") end


text:
type AFArray
  ptr::af_array

  function AFArray(ptr)
    a = new(ptr)
    finalizable() && finalizer(a, release)
    a
  end
end

function AFArray{T,N}(::Type{T}, dims::NTuple{N,Int})
  out = af_array[0]
  dims = dim_t[dims...]
  af_create_handle(out, N, dims, aftype(T))
  AFArray(out[1])
end
AFArray{T}(::Type{T}, dims...) = AFArray(T, dims)

function AFArray{T,N}(data::Array{T,N})
  out = af_array[0]
  dims = dim_t[size(data)...]
  af_create_array(out, data, N, dims, aftype(T))
  AFArray(out[1])
end

Base.unsafe_convert(::Type{af_array}, a::AFArray) = a.ptr

Base.show(io::IO, in::AFArray) = show(io, to_host(in))

function jl_size(in::AFArray)
  dims = dim_t[0, 0, 0, 0]
  af_get_dims(pointer(dims,1), pointer(dims,2), pointer(dims,3), pointer(dims,4), in)
  dims[4] > 1 && return Int(dims[1]), Int(dims[2]), Int(dims[3]), Int(dims[4])
  dims[3] > 1 && return Int(dims[1]), Int(dims[2]), Int(dims[3])
  dims[2] > 1 && return Int(dims[1]), Int(dims[2])
  (Int(dims[1]), )
end
jl_size(in::AFArray, dim::Int) = jl_size(in)[dim]

function Base.size(in::AFArray)
  dims = dim_t[0, 0, 0, 0]
  af_get_dims(pointer(dims,1), pointer(dims,2), pointer(dims,3), pointer(dims,4), in)
  Int(dims[1]), Int(dims[2]), Int(dims[3]), Int(dims[4])
end
Base.size(in::AFArray, dim::Int) = size(in)[dim]

function Base.length(in::AFArray)
  out = dim_t[0]
  af_get_elements(out, in)
  Int(out[1])
end

function Base.eltype(in::AFArray)
  out = af_dtype[0]
  af_get_type(out, in)
  dtype = out[1]
  dtype == f32 && return Float32
  dtype == c32 && return Complex{Float32}
  dtype == f64 && return Float64
  dtype == c64 && return Complex{Float64}
  dtype == b8 && return Bool
  dtype == s32 && return Int32
  dtype == u32 && return UInt32
  dtype == u8 && return UInt8
  dtype == s64 && return Int64
  dtype == u64 && return UInt64
  throw("error: eltype")
end

function Base.ndims(in::AFArray)
  out = Cuint[0]
  af_get_numdims(out, in)
  Int(out[1])
end

Base.similar(in::AFArray) = AFArray(eltype(in), size(in))

function copy(in::AFArray)
  out = af_array[0]
  af_copy_array(out, in)
  AFArray(out[1])
end

eval(in::AFArray) = af_eval(in)

function to_host(in::AFArray)
  out = Array(eltype(in), jl_size(in))
  af_get_data_ptr(out, in)
  out
end

function refcount(in::AFArray)
  out = Int32[0]
  af_get_data_ref_count(out, in)
  Int(out[1])
end

release(in::AFArray) = af_release_array(in)

function retain(in::AFArray)
  out = af_array[0]
  af_retain_array(out, in)
  AFArray(out[1])
end

function fill{T}(value::T, dims)
  out = af_array[0]
  dims = dim_t[dims...]
  af_constant(out, value, length(dims), dims, aftype(T))
  AFArray(out[1])
end

zeros{T}(::Type{T}, dims) = fill(T(0), dims)
zeros{T}(::Type{T}, dims...) = fill(T(0), dims)
zeros(in::AFArray) = zeros(eltype(in), size(in))

ones{T}(::Type{T}, dims) = fill(T(1), dims)
ones{T}(::Type{T}, dims...) = fill(T(1), dims)
ones(in::AFArray) = ones(eltype(in), size(in))

function range{T,N}(::Type{T}, dims::NTuple{N,Int}, seq_dim::Int=0)
  out = af_array[0]
  dims = dim_t[dims...]
  af_range(out, N, dims, seq_dim-1, aftype(T))
  AFArray(out[1])
end
range{T}(::Type{T}, dims::Int...) = range(T, dims)

function rand{T}(::Type{T}, dims)
  out = af_array[0]
  dims = dim_t[dims...]
  af_randu(out, length(dims), dims, aftype(T))
  AFArray(out[1])
end
rand{T}(::Type{T}, dims...) = rand(T, dims)

function randn{T}(::Type{T}, dims)
  out = af_array[0]
  dims = dim_t[dims...]
  af_randn(out, length(dims), dims, aftype(T))
  AFArray(out[1])
end
randn{T}(::Type{T}, dims...) = randn(T, dims)

#af_write_array

#function eye{T,N}(::Type{AFArray{T}}, dims::NTuple{N,Int})
#  out = af_array[0]
#  dims = dim_t[dims...]
#  af_identity(out, N, dims, aftype(T))
#  AFArray(out[1])
#end
#eye{T}(::Type{AFArray{T}}, dims...) = eye(T, dims)

#function iota{T,N}(::Type{AFArray{T}}, dims::NTuple{N,Int}, tdims::NTuple{N,Int})
#  out = af_array[0]
#  dims = dim_t[dims...]
#  tdims = dim_t[tdims...]
#  af_iota(out, N, dims, N, tdims, aftype(T))
#  AFArray(out[1])
#end

function cast{T}(in::AFArray, ::Type{T})
  out = af_array[0]
  af_cast(out, in, aftype(T))
  AFArray(out[1])
end

function Base.vec(in::AFArray)
  out = af_array[0]
  af_flat(out, in)
  AFArray(out[1])
end

function flip(in::AFArray, dim::Int)
  dim > 0 || throw("Invalid dimension: $(dim).")
  out = af_array[0]
  af_flip(out, in, dim-1)
  AFArray(out[1])
end

function Base.cat(ins::Vector{AFArray}, dim::Int)
  dim > 0 || throw("Invalid dimension: $(dim).")
  length(ins) == 0 && throw("Invalid input length: $(length(ins)).")
  out = af_array[0]
  af_join_many(out, dim-1, length(ins), ins)
  AFArray(out[1])
end

function Base.reshape(in::AFArray, dims)
  out = af_array[0]
  dims = dim_t[dims...]
  af_moddims(out, in, length(dims), dims)
  AFArray(out[1])
end
Base.reshape(in::AFArray, dims...) = reshape(in, dims)

function moddims(in::AFArray, dims)
  out = af_array[0]
  dims = dim_t[dims...]
  af_moddims(out, in, length(dims), dims)
  AFArray(out[1])
end

function reorder(in::AFArray, x::Int, y::Int, z::Int, w::Int)
  out = af_array[0]
  af_reorder(out, in, x, y, z, w)
  AFArray(out[1])
end

function Base.transpose(in::AFArray; conjugate::Bool=false)
  out = af_array[0]
  af_transpose(out, in, conjugate)
  AFArray(out[1])
end

function lookup(in::AFArray, indices::AFArray, dim::Int)
  out = af_array[0]
  af_lookup(out, in, indices, dim-1)
  AFArray(out[1])
end


text:
using BinaryBuilder, Pkg

name = "MKL"
version = v"2021.1.1"

# Bash recipes for building across all platforms
script = read(joinpath(@__DIR__, "script.sh"), String)
script_macos = read(joinpath(@__DIR__, "script_macos.sh"), String)

non_reg_ARGS = filter(arg -> arg != "--register", ARGS)

platform_sources = [
    (
        platform = Platform("x86_64", "linux"; libc="glibc"),
        source = ArchiveSource(
            "https://anaconda.org/intel/mkl/2021.1.1/download/linux-64/mkl-2021.1.1-intel_52.tar.bz2",
            "bfb0fd056576cad99ae1d9c69ada2745420da9f9cf052551d5b91f797538bda2";
            unpack_target = "mkl-x86_64-linux-gnu"
        ),
        # We need to run autofix on Linux, because here libmkl_rt doesn't
        # have a soname, so we can't ccall it without specifying the path:
        # https://github.com/JuliaSparse/Pardiso.jl/issues/69
        autofix = true,
        script = script,
    ),
    (
        platform = Platform("i686", "linux"; libc="glibc"),
        source = ArchiveSource(
            "https://anaconda.org/intel/mkl/2021.1.1/download/linux-32/mkl-2021.1.1-intel_52.tar.bz2",
            "7b6f55a30886154bd96d4b4c6b7428494a59397b87779b58e5b3de00250343f9";
            unpack_target = "mkl-i686-linux-gnu"
        ),
        autofix = true,
        script = script,
    ),
    (
        platform = Platform("x86_64", "macos"),
        source = ArchiveSource(
            "https://anaconda.org/intel/mkl/2021.1.1/download/osx-64/mkl-2021.1.1-intel_50.tar.bz2",
            "819fb8875909d4d024e2a936c54b561aebd1e3aebe58fc605c70aa1ad9a66b70";
            unpack_target = "mkl-x86_64-apple-darwin14"
        ),
        # Need to disable autofix: updating linkage of libmkl_intel_thread.dylib on
        # macOS causes runtime issues:
        # https://github.com/JuliaPackaging/Yggdrasil/issues/915.
        autofix = false,
        script = script_macos,
    ),
    (
        platform = Platform("i686", "windows"),
        source = ArchiveSource(
            "https://anaconda.org/intel/mkl/2021.1.1/download/win-32/mkl-2021.1.1-intel_52.tar.bz2",
            "dba6a12a481407ec55fba9895b68afacb15f044905dcb5e185db341b688e6177";
            unpack_target = "mkl-i686-w64-mingw32"
        ),
        autofix = false,
        script = script,
    ),
    (
        platform = Platform("x86_64", "windows"),
        source = ArchiveSource(
            "https://anaconda.org/intel/mkl/2021.1.1/download/win-64/mkl-2021.1.1-intel_52.tar.bz2",
            "4024391b8a45836d5a7ee92405b7767874b3c3bbf2f490349fda042db3b60dfd";
            unpack_target = "mkl-x86_64-w64-mingw32"
        ),
        autofix = false,
        script = script,
    ),
]

# The products that we will ensure are always built
products = [
    LibraryProduct(["libmkl_core", "mkl_core", "mkl_core.1"], :libmkl_core),
    LibraryProduct(["libmkl_rt", "mkl_rt", "mkl_rt.1"], :libmkl_rt),
]

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

non_reg_ARGS = filter(arg -> arg != "--register", ARGS)
include("../../fancy_toys.jl")
filter!(p -> should_build_platform(triplet(first(p))), platform_sources)

for (idx, (platform, source, autofix, script)) in enumerate(platform_sources)
    # Use "--register" only on the last invocation of build_tarballs
    if idx < length(platform_sources)
        args = non_reg_ARGS
    else
        args = ARGS
    end
    build_tarballs(args, name, version, [source], script, [platform], products, dependencies; lazy_artifacts = true, autofix = autofix)
end

max_stars_repo_name: plut/Yggdrasil

max_issues_repo_name: plut/Yggdrasil

max_forks_repo_name: plut/Yggdrasil


text:
"""
    @option [alias::String] <struct def>

Define an option struct type. This will auto-generate methods that parse a given `Dict{String}`
object (the keys must be of type `String`) into an instance of the struct type you defined. One
can use `alias` string to distinguish multiple possible option type for the same field.

# Special Types

- `Maybe{T}`: this type is equivalent to `Union{Nothing, T}` and
    is treated specially in `@option`, it will always have a default
    value of `nothing` if not specified.
- [`Reflect`](@ref): this type is treated specially to allow one to use a field
    to store the corresponding type information.

!!! compat "Configurations 0.16"
    from v0.16.0 Configurations stops overloading the `Base.show` method for you,
    if you need pretty printing of your option types, consider overloading
    the `Base.show(io::IO, mime::MIME, x)` method to `pprint_struct(io, mime, x)` provided by
    [GarishPrint](https://github.com/Roger-luo/GarishPrint.jl)

!!! compat "Configurations 0.12"
    from v0.12.0 the field alias feature is removed due to the syntax conflict with
    field docstring. Please refer to [#17](https://github.com/Roger-luo/Configurations.jl/issues/17).

# Example

One can define option type via `@option` macro with or without an alias.

```julia-repl
julia> "Option A"
       @option "option_a" struct OptionA
           name::String
           int::Int = 1
       end

julia> "Option B"
       @option "option_b" struct OptionB
           opt::OptionA = OptionA(;name = "Sam")
           float::Float64 = 0.3
       end
julia> option = from_dict(OptionB, d)
OptionB(;
    opt = OptionA(;
        name = "Roger",
        int = 2,
    ),
    float = 0.33,
)
```

when there are multiple possible option type for one field,
one can use the alias to distinguish them

```julia-repl
julia> @option struct OptionD
           opt::Union{OptionA, OptionB}
       end

julia> d1 = Dict{String, Any}(
               "opt" => Dict{String, Any}(
                   "option_b" => d
               )
           );

julia> from_dict(OptionD, d1)
OptionD(;
    opt = OptionB(;
        opt = OptionA(;
            name = "Roger",
            int = 2,
        ),
        float = 0.33,
    ),
)
```
"""
macro option(ex)
    return esc(option_m(__module__, ex))
end

macro option(alias::String, ex)
    return esc(option_m(__module__, ex, alias))
end

function option_m(mod::Module, ex, type_alias=nothing)
    ex = macroexpand(mod, ex)
    def = JLKwStruct(ex, type_alias)
    return codegen_option_type(mod, def)
end

"""
    validate_option_def(def::JLKwStruct)

Validate the option definition.
"""
function validate_option_def(mod::Module, def::JLKwStruct)
    if def.typealias !== nothing
        isempty(def.typevars) ||
            throw(ArgumentError(
                "only concrete type definition can have type alias"
        ))
    end
    has_duplicated_reflect_type(mod, def) &&
        throw(ArgumentError("struct fields contain duplicated `Reflect` type"))
    return
end

"""
    codegen_option_type(mod::Module, def::JLKwStruct)

Generate the `Configurations` option type definition from
a given `JLKwStruct` created by [`Expronicon`](https://github.com/Roger-luo/Expronicon.jl).
"""
function codegen_option_type(mod::Module, def::JLKwStruct)
    # preprocess
    validate_option_def(mod, def)
    add_field_defaults!(mod, def)

    quote
        $(codegen_ast(def))
        Core.@__doc__ $(def.name)
        $(codegen_create(def))
        $(codegen_is_option(def))
        $(codegen_convert(def))
        $(codegen_field_default(def))
        $(codegen_type_alias(def))
        $(codegen_isequal(def))
        $(codegen_from_dict_specialize(def))
        nothing
    end
end

"""
    add_field_defaults!(m::Module, def::JLKwStruct)

Add default value for `Maybe` and `Reflect` type.
"""
function add_field_defaults!(m::Module, def::JLKwStruct)
    for field in def.fields
        if is_reflect_type_expr(m, field.type)
            field.default = Reflect()
        elseif is_maybe_type_expr(m, field.type) && field.default === no_default
            field.default = nothing
        end
    end
    return def
end

"""
    has_duplicated_reflect_type(m::Module, def::JLKwStruct)

Check if the definition has duplicated reflect type.
"""
function has_duplicated_reflect_type(m::Module, def::JLKwStruct)
    has_reflect_type = false
    for field in def.fields
        if is_reflect_type_expr(m, field.type)
            has_reflect_type && return true
            has_reflect_type = true
        end
    end
    return false
end

"""
    is_reflect_type_expr(m::Module, @nospecialize(ex))

Check if the expression `ex` evaluates to a [`Reflect`](@ref).
"""
function is_reflect_type_expr(m::Module, @nospecialize(ex))
    if isdefined(m, :Reflect) && (getfield(m, :Reflect) === Reflect)
        ex === :Reflect && return true
    end
    # no need to check definition
    ex == Reflect && return true
    ex == GlobalRef(Configurations, :Reflect) && return true
    ex == :($Configurations.Reflect) && return true
    ex == :($Configurations.$Reflect) && return true
    ex == :(Configurations.$Reflect) && return true
    if isdefined(m, :Configurations)
        ex == :(Configurations.Reflect) && return true
    end
    return false
end

"""
    is_maybe_type_expr(m::Module, @nospecialize(ex))

Check if the expression `ex` evaluates to a `Maybe{T}`.
"""
function is_maybe_type_expr(m::Module, @nospecialize(ex))
    if isdefined(m, :Maybe) && (getfield(m, :Maybe) === Maybe)
        _is_maybe_type_expr(ex) && return true
    end

    if ex isa GlobalRef && ex.mod === Configurations
        ex.name === :Maybe && return true
    end

    if ex isa Type && ex isa Union && Nothing <: ex
        return true
    end

    ex isa Expr || return false
    if ex.head === :.
        if ex.args[1] === Configurations || ex.args[1] === :Configurations
            return _is_maybe_type_expr(ex.args[2])
        end
    elseif ex.head === :curly
        return is_maybe_type_expr(m, ex.args[1])
    end
    return false
end

function _is_maybe_type_expr(@nospecialize(ex))
    ex === Maybe && return true
    ex === :Maybe && return true
    if ex isa QuoteNode
        ex.value === :Maybe && return true
        ex.value isa Type && ex.value <: Maybe && return true
    end

    ex isa Expr || return false
    if ex.head === :curly
        ex.args[1] === :Maybe && return true
        ex.args[1] isa Type && ex.args[1] <: Maybe && return true
    end
    return false
end

"""
    codegen_is_option(x::JLKwStruct)

Generate the [`is_option`](@ref) method.
"""
function codegen_is_option(def::JLKwStruct)
    quote
        $Configurations.is_option(::$(def.name)) = true
        $Configurations.is_option(::Type{<:$(def.name)}) = true
    end
end

"""
    codegen_convert(x::JLKwStruct)

Generate `Base.convert` from `AbstractDict{String}` to the given option type.
"""
function codegen_convert(def::JLKwStruct)
    quote
        function $Base.convert(::Type{<:$(def.name)}, d::AbstractDict{String})
            return $Configurations.from_dict($(def.name), d)
        end
    end
end

"""
    codegen_field_default(def::JLKwStruct)

Generate [`field_default`](@ref) overload to support the default value
reflection.
"""
function codegen_field_default(def::JLKwStruct)
    obj = gensym(:x)
    msg = Expr(:string, "type $(def.name) does not have field ", obj)
    err = :(error($msg))
    ret = JLIfElse()
    ret.otherwise = err

    prev_field_names = Symbol[]

    for (k, field) in enumerate(def.fields)
        vars = Symbol[]

        for name in prev_field_names
            if has_symbol(field.default, name)
                push!(vars, name)
            end
        end

        cond = :($obj == $(QuoteNode(field.name)))

        if isempty(vars) # const default
            ret[cond] = field.default
        else
            jlfn = JLFunction(; args=vars, body=field.default)
            fn = gensym(:fn)
            ret[cond] = quote
                $fn = $(codegen_ast(jlfn))
                $PartialDefault($fn, $vars, $(Expr(:quote, field.default)))
            end
        end

        push!(prev_field_names, field.name)
    end

    type = gensym(:type)
    typevars = name_only.(def.typevars)
    ub = isempty(def.typevars) ? def.name : Expr(:curly, def.name, typevars...)

    return codegen_ast(
        JLFunction(;
            name=:($Configurations.field_default),
            args=[:(::Type{$type}), :($obj::Symbol)],
            body=codegen_ast(ret),
            whereparams=[typevars..., :($type <: $ub)],
        ),
    )
end

# function codegen_field_alias(def::JLKwStruct)
# end

"""
    codegen_isequal(x::JLKwStruct)

Generate `Base.:(==)` to overload comparison operator to [`compare_options`](@ref)
for given option type.
"""
function codegen_isequal(x::JLKwStruct)
    return :(Base.:(==)(a::$(x.name), b::$(x.name)) = $compare_options(a, b))
end

"""
    codegen_type_alias(def::JLKwStruct)

Generate type alias method [`type_alias`](@ref).
"""
function codegen_type_alias(def::JLKwStruct)
    @gensym TYPE_ALIAS_MAP
    ret = quote
        const $TYPE_ALIAS_MAP = Dict{String, Any}()
        $Configurations.get_type_alias_map(::Type{<:$(def.name)}) = $TYPE_ALIAS_MAP

        # the type can only be a concrete type if it has alias
        # if the type is not concrete we will return nothing
        $Configurations.type_alias(::Type{<:$(def.name)}) = $(def.typealias)
    end

    if !(def.typealias === nothing)
        push!(ret.args, :($TYPE_ALIAS_MAP[$(def.typealias)] = $(def.name)))
    end
    return ret
end

"""
    codegen_create(def::JLKwStruct)

Generate [`Configurations.create`](@ref) overload.
"""
function codegen_create(def::JLKwStruct)
    return codegen_ast_kwfn(def, :($Configurations.create))
end

"""
    codegen_from_dict_specialize(def::JLKwStruct)

Generate the specialized `from_dict` for the given definition.
"""
function codegen_from_dict_specialize(def::JLKwStruct)
    quote
        @generated function $Configurations.from_dict_specialize(
            ::Type{T}, d::AbstractDict{String}
        ) where {T<:$(def.name)}
            return $Configurations.from_dict_generated(T, :d)
        end
    end
end

"""
    @type_alias <type> <name::String>

Define a type alias for option type `type`. The corresponding
`type` must be a concrete type in order to map this Julia type
to a human readable markup language (e.g TOML, YAML, etc.).
"""
macro type_alias(type, name::String)
    esc(type_alias_m(type, name))
end

function type_alias_m(type, name::String)
    quote
        $Configurations.type_alias(::Type{$type}) = $name
        $Configurations.set_type_alias($type, $name)
    end
end


text:
"""
    function get_PTDF_load(sys::System, PTDF_matrix ::Array)

Takes the full PTDF matrix of the system and only selects the columns that correspond to the
load nodes (PQ). It returns the parsed PTDF Matrix sorted by bus.

# Arguments
- `sys::System`:                            Power system in p.u. (from PowerSystems.jl)
- `PTDF_matrix::Any`:                       PTDF matrix as a direct output from
                                            PTDF_matrix = PTDF(sys), see PowerSystems.jl

# Throws
- `ERROR`:                                  PTDF_matrix has no field data. The argument must
                                            be a direct output from PTDF_matrix = PTDF(sys)
                                            see PowerSystems.jl

"""
function get_PTDF_load(sys::System, PTDF_matrix ::Array)

    # Get components of the Bus
    buses = get_components(Bus, sys)
    numberof_PQnodes = 0
    (nl,nb)=size(PTDF_matrix)
    PQ_buses = Array{Int64}(undef,0)
    lines = get_components(Line, sys)
    PTDF_PQLoad = Array{Float64}(undef,(nl,0))
    for bus in buses
        if get_bustype(bus) == BusTypes.PQ
            numberof_PQnodes = numberof_PQnodes + 1
            resize!(PQ_buses, numberof_PQnodes)
            PQ_buses[numberof_PQnodes] = bus.number
        end
    end
    if numberof_PQnodes == 0
        PTDF_PQLoad = Array{Float64}(undef,(nl,0))
    else
        sort!(PQ_buses)
        PTDF_PQLoad = PTDF_matrix[:,PQ_buses]
    end

    return PQ_buses, PTDF_PQLoad
end

max_stars_repo_name: invenia/MaximumGeneratorProfitTCRDD

max_issues_repo_name: invenia/MaximumGeneratorProfitTCRDD

max_forks_repo_name: invenia/MaximumGeneratorProfitTCRDD


text:
# This file is for different constraints

abstract type AbstractLinearProgram end
struct StandardLP            <: AbstractLinearProgram end
struct LinearRelaxedLP       <: AbstractLinearProgram end
struct TriangularRelaxedLP   <: AbstractLinearProgram end
struct BoundedMixedIntegerLP <: AbstractLinearProgram end
struct SlackLP               <: AbstractLinearProgram end
struct MixedIntegerLP        <: AbstractLinearProgram end

Base.Broadcast.broadcastable(LP::AbstractLinearProgram) = Ref(LP)

function init_model_vars(model::Model,
                         problem::Problem,
                         encoding::AbstractLinearProgram
                         ;
                         kwargs...)
    model[:network] = problem.network
    model[:input]   = problem.input
    model[:output]  = problem.output

    init_vars(model, problem.network, :z, with_input=true)
    for (k, v) in kwargs
        model[k] = v
    end
    _init_unique(model, problem, encoding)
    model
end

_init_unique(m::Model, prob::Problem, encoding::AbstractLinearProgram) = nothing
_init_unique(m::Model, prob::Problem, encoding::SlackLP) = init_vars(m, prob.network, :slack)
_init_unique(m::Model, prob::Problem, encoding::MixedIntegerLP) = init_vars(m, prob.network, :δ, binary=true)
_init_unique(m::Model, prob::Problem, encoding::TriangularRelaxedLP) = _insert_bounds(m, prob, encoding)
function _init_unique(m::Model, prob::Problem, encoding::BoundedMixedIntegerLP)
    init_vars(m, prob.network, :δ, binary=true)
    _insert_bounds(m, prob, encoding)
end

function _insert_bounds(m::Model, prob::Problem, encoding::Union{TriangularRelaxedLP, BoundedMixedIntegerLP})
    if !haskey(object_dictionary(m), :bounds)
        before_act = get!(object_dictionary(m), :before_act, true)
        m[:bounds] = get_bounds(prob, !before_act)
    end
end

model_params(LP::BoundedMixedIntegerLP, m::Model, i::Integer) = (_ẑᵢ₊₁(m, i), m[:z][i+1], m[:δ][i], _ẑᵢ₊₁_bound(m, i)...)
model_params(LP::TriangularRelaxedLP,   m::Model, i::Integer) = (_ẑᵢ₊₁(m, i), m[:z][i+1], _ẑᵢ₊₁_bound(m, i)...)
model_params(LP::StandardLP,            m::Model, i::Integer) = (_ẑᵢ₊₁(m, i), m[:z][i+1], m[:δ][i])
model_params(LP::LinearRelaxedLP,       m::Model, i::Integer) = (_ẑᵢ₊₁(m, i), m[:z][i+1], m[:δ][i])
model_params(LP::MixedIntegerLP,        m::Model, i::Integer) = (_ẑᵢ₊₁(m, i), m[:z][i+1], m[:δ][i], -m[:M], m[:M])
model_params(LP::SlackLP,               m::Model, i::Integer) = (_ẑᵢ₊₁(m, i), m[:z][i+1], m[:δ][i], m[:slack][i])

# helper function to get the bounds relevant for the i-th layer constraint.
# If before_act=true, then the model is storing the pre-activation
# bounds, and we return the i+1th bound (goes with ẑᵢ₊₁). Otherwise
# NOTE we forward propagate (affine map) the ith bound. We assume the model
# stores post-activation bounds since this is the current default for get_bounds
function _ẑᵢ₊₁_bound(m, i)
    layer = m[:network].layers[i]
    # let's assume we're post-activation if the parameter isn't set,
    # since that's the default for get_bounds
    if get(object_dictionary(m), :before_act, false)
        ẑ_bound = m[:bounds][i+1]
    else
        ẑ_bound = approximate_affine_map(layer, m[:bounds][i])
    end
    low(ẑ_bound), high(ẑ_bound)
end
_ẑᵢ₊₁(m, i) = affine_map(m[:network].layers[i], m[:z][i])

# Any encoding passes through here first:
function encode_network!(model::Model, network::Network, encoding::AbstractLinearProgram)
    get!(object_dictionary(model), :network, network)
    for (i, layer) in enumerate(network.layers)
        encode_layer!(encoding, model, layer, model_params(encoding, model, i)...)
    end
end

function encode_layer!(LP::AbstractLinearProgram, model::Model, i::Integer)
    L = model[:network].layers[i]
    encode_layer!(LP, model, L, model_params(LP, model, i)...)
    nothing
end

# For an Id Layer, any encoding type defaults to this:
function encode_layer!(::AbstractLinearProgram, model::Model, layer::Layer{Id}, ẑᵢ, zᵢ, args...)
    @constraint(model, zᵢ .== ẑᵢ)
    nothing
end

function encode_layer!(LP::AbstractLinearProgram, model::Model, layer::Layer{ReLU}, args...)
    encode_relu.(LP, model, args...)
    nothing
end

# SlackLP is slightly different, because we need to keep track of the slack variables
function encode_layer!(SLP::SlackLP, model::Model, layer::Layer{Id}, ẑᵢ, zᵢ, δᵢⱼ, sᵢ)
    @constraint(model, zᵢ .== ẑᵢ)
    # We need identity layer slack variables so that the algorithm doesn't
    # "get confused", but they are set to 0 because they're not relevant
    @constraint(model, sᵢ .== 0.0)
    return nothing
end

function encode_ij(LP, model, i, j)
    L = model[:network].layers[i]
    params = model_params(LP, model, i)
    if L.activation isa Id
        ẑᵢⱼ = _ẑᵢ₊₁(model, i)[j]
        zᵢⱼ = model[:z][i+1][j]
        @constraint(model, ẑᵢⱼ == zᵢⱼ)
    elseif L.activation isa ReLU
        # ridiculous hack so it works with scalar params (like :M)
        _getindex_hack(x, j) = x isa Number ? x : x[j]
        encode_relu(LP, model, _getindex_hack.(params, j)...)
    end
end

function encode_relu(::SlackLP, model, ẑᵢⱼ, zᵢⱼ, δᵢⱼ, sᵢⱼ)
    if δᵢⱼ
        @constraint(model, zᵢⱼ == ẑᵢⱼ + sᵢⱼ)
        @constraint(model, ẑᵢⱼ + sᵢⱼ >= 0.0)
    else
        @constraint(model, zᵢⱼ == sᵢⱼ)
        @constraint(model, ẑᵢⱼ <= sᵢⱼ)
    end
end

function encode_relu(::BoundedMixedIntegerLP, model, ẑᵢⱼ, zᵢⱼ, δᵢⱼ, l̂ᵢⱼ, ûᵢⱼ)
    if l̂ᵢⱼ >= 0.0
        @constraint(model, zᵢⱼ == ẑᵢⱼ)
    elseif ûᵢⱼ <= 0.0
        @constraint(model, zᵢⱼ == 0.0)
    else
        @constraints(model, begin
                                zᵢⱼ >= 0.0
                                zᵢⱼ >= ẑᵢⱼ
                                zᵢⱼ <= ûᵢⱼ * δᵢⱼ
                                zᵢⱼ <= ẑᵢⱼ - l̂ᵢⱼ * (1 - δᵢⱼ)
                            end)
    end
end

function encode_relu(::MixedIntegerLP, args...)
    encode_relu(BoundedMixedIntegerLP(), args...)
end

function encode_relu(::TriangularRelaxedLP, model, ẑᵢⱼ, zᵢⱼ, l̂ᵢⱼ, ûᵢⱼ)
    if l̂ᵢⱼ > 0.0
        @constraint(model, zᵢⱼ == ẑᵢⱼ)
    elseif ûᵢⱼ < 0.0
        @constraint(model, zᵢⱼ == 0.0)
    else
        @constraints(model, begin
                                zᵢⱼ >= 0.0
                                zᵢⱼ >= ẑᵢⱼ
                                zᵢⱼ <= (ẑᵢⱼ - l̂ᵢⱼ) * ûᵢⱼ / (ûᵢⱼ - l̂ᵢⱼ)
                            end)
    end
end

function encode_relu(::LinearRelaxedLP, model, ẑᵢⱼ, zᵢⱼ, δᵢⱼ)
    @constraint(model, zᵢⱼ == (δᵢⱼ ? ẑᵢⱼ : 0.0))
end

function encode_relu(::StandardLP, model, ẑᵢⱼ, zᵢⱼ, δᵢⱼ)
    if δᵢⱼ
        @constraint(model, ẑᵢⱼ >= 0.0)
        @constraint(model, zᵢⱼ == ẑᵢⱼ)
    else
        @constraint(model, ẑᵢⱼ <= 0.0)
        @constraint(model, zᵢⱼ == 0.0)
    end
end

#=
Add input/output constraints to model
=#
function add_complementary_set_constraint!(model::Model, output::HPolytope, z::Vector{VariableRef})
    out_A, out_b = tosimplehrep(output)
    # Needs to take the complementary of output constraint
    n = length(constraints_list(output))
    if n == 1
        # Here the output constraint is a half space
        halfspace = first(constraints_list(output))
        add_complementary_set_constraint!(model, halfspace, z)
    else
        error("Non-convex constraints are not supported. Please make sure that the
            output set is a HalfSpace (or an HPolytope with a single constraint) so that the
            complement of the output is convex. Got $n constraints.")
    end
    return nothing
end

function add_complementary_set_constraint!(m::Model, H::HalfSpace, z::Vector{VariableRef})
    a, b = tosimplehrep(H)
    @constraint(m, a * z .>= b)
    return nothing
end
function add_complementary_set_constraint!(m::Model, PC::PolytopeComplement, z::Vector{VariableRef})
    add_set_constraint!(m, PC.P, z)
    return nothing
end

function add_set_constraint!(m::Model, set::Union{HPolytope, HalfSpace}, z::Vector{VariableRef})
    A, b = tosimplehrep(set)
    @constraint(m, A * z .<= b)
    return nothing
end

function add_set_constraint!(m::Model, set::Hyperrectangle, z::Vector{VariableRef})
    @constraint(m, z .<= high(set))
    @constraint(m, z .>= low(set))
    return nothing
end

function add_set_constraint!(m::Model, PC::PolytopeComplement, z::Vector{VariableRef})
    add_complementary_set_constraint!(m, PC.P, z)
    return nothing
end


text:
# We use [dcraw](http://www.cybercom.net/~dcoffin/dcraw/) to extract from RAW images
# the blue channel pixels and convert to 16bit pgm without exposure manipulation.

# The website notes: "Unless otherwise noted in the source code, these programs
# are free for all uses, although I would like to receive credit for them."
# The program dcraw is written by Dave Coffin.

# On mac, compile dcraw.c with
# 	`llvm-gcc -o dcraw dcraw.c -lm -DNO_JPEG -DNO_LCMS -DNO_JASPER`
# source [http://vkphotoblog.blogspot.be/2014/05/dcraw-921-for-os-x-mavericks-users.html]

# ### dcraw options
# * -W: avoids stretching ("No matter how dark an image is, dcraw's auto-exposure stretches it so
#       that one percent of its pixels appear white. The "-W" option avoids this behavior.")
# * -4: is for linear 16-bit, same as -6 -W -g 1 1 (with -g for gamma correction)
# * -j: don't stretch or rotate raw pixels
# * -t [0-7]: flip image (0=none)
# * -D: document mode without scaling (totally raw), while -d scales to 16bit (eg from 14bit).
#       Either use -D and then reinterpret in julia or use -d. We now use -d otherwise we need to
#       extract image property for bit depth (12bit, 14bit, ...).
# * -r 0 0 0 1: select only the blue channel (see Brusa and Bunker, 2014). This option selects from a bayer RGGB layout.
# * -v: for verbose output

const OVEREXP = 0.005

const RAW_EXT = String[".3fr", ".ari", ".arw", ".bay", ".crw", ".cr2",
".cap", ".dcs", ".dcr", ".dng",
".drf", ".eip", ".erf", ".fff", ".iiq", ".k25", ".kdc", ".mdc", ".mef", ".mos", ".mrw",
".nef", ".nrw", ".obm", ".orf", ".pef", ".ptx", ".pxn", ".r3d", ".raf", ".raw", ".rwl",
".rw2", ".rwz", ".sr2", ".srf", ".srw", ".tif", ".x3f"]

const DCRAW_DIR = joinpath(dirname(pathof(LeafAreaIndex)), "dcraw")

@static if Sys.iswindows()
    const DCRAW_EXE = joinpath(DCRAW_DIR, "dcraw-9.26-ms-64-bit.exe")
end

@static if Sys.isunix()
    const DCRAW_EXE = joinpath(DCRAW_DIR, "dcraw")
    is_user_executable(file) = isodd(uperm(file))
    is_user_executable(DCRAW_EXE) || chmod(DCRAW_EXE, 0o755)
end

@assert isfile(DCRAW_EXE)

"""Read in the raw blue channel from a raw image. The image is copied to
a temporary directory and run through dcraw with options `-d -4 -j -t 0 -r 0 0 0 1`."""
rawblueread(filepath::String; kwargs...) = readraw(filepath; kwargs...)

rawcolourread(filepath::String; kwargs...) = readraw(filepath; colour=true, kwargs...)

function readraw(filepath::String; overwrite=false, rmcopy=true, rmpxm=true,
    destdir=joinpath(tempdir(), "pxm"), colour=false, kwargs...)

    @assert isfile(filepath)
    ext = lowercase(last(splitext(filepath)))
    @assert ext ∈ RAW_EXT
    isdir(destdir) || mkdir(destdir)

    if colour
        raw2pxm = f -> run(`$DCRAW_EXE -4 -j -t 0 $f`)
        pxm_ext = ".ppm"
    else
        #raw2pxm = f -> run(`$DCRAW_EXE -d -4 -j -t 0 -r 0 0 0 1 $f`)
        #pxm_ext = ".pgm"
		# add De-Bayering
		raw2pxm = f -> run(`$DCRAW_EXE -4 -j -t 0 -q 3 -o 0 $f`)
		pxm_ext = ".ppm"
    end

    copyfile = joinpath(destdir, last(splitdir(filepath)))
    pxm = first(splitext(copyfile)) * pxm_ext

    if !overwrite && isfile(pxm)
        return FileIO.load(pxm)
    end

    cp(filepath, copyfile; force=true)
    raw2pxm(copyfile)
    img = FileIO.load(pxm)
    rmcopy && rm(copyfile)
    rmpxm && rm(pxm)
    img
end

isoverexposed(img::AbstractArray) = sum(img .== 1) > OVEREXP * length(img)
isoverexposed(polim::PolarImage) = isoverexposed(pixels(polim))


text:
function mkelconcg(dgnodes)

# remove duplicate nodes in mesh.p1
ns,dim,nt = size(dgnodes);
A = reshape(permutedims(dgnodes,[1,3,2]),ns*nt,dim);

# B = flip(A,1);
# [~,I]=unique(B,'rows');
# B = B(sort(I),:);
# B = flip(B,1);
# B = unique(A,dims=1);
# b = xiny(A,B,1);

B,~,b = uniquerows(A);
# display(maximum(abs.(b-b1)))
# ii = 100:200;
# display(b[ii]')
# display(b1[ii]')
# #display(b1[31:60]')
# error("here")

# CG mesh
cgnodes = B';
cgelcon = reshape(b,ns, nt);

return cgnodes,cgelcon

end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
"""To read a simple edge list, each line being an (i, j) pair"""
function readIJ(filename::AbstractString, sep=',')
  edlist = readdlm(filename, sep)
  n = maximum(edlist)
  m = size(edlist)
  edlist = convert(Array{Int64,2}, edlist)
  a = sparse(edlist[:,1],edlist[:,2],ones(m[1]),n,n)
  a = a + a'

  return a
end # readIJ

"""To read a simple edge list, each line being an (i, j, v) pair.
The parens should not be there in the format, just commas separating.
To generate this format in Matlab, you just need to be careful to write the
vertex indices with sufficient precision.  For example, you can do this

```
>> [ai,aj,av] = find(triu(a));
>> dlmwrite('graph.txt',[ai,aj,av],'precision',9);
```
"""
function readIJV(filename::AbstractString)
  data = readdlm(filename,',')
  n = maximum(data[:,1:2])
  m = size(data)
  edlist = convert(Array{Int64,2}, data[:,1:2])
  wts = convert(Array{Float64,1}, data[:,3])

  a = sparse(edlist[:,1],edlist[:,2],wts,n,n)
  a = a + a'

  return a
end # readIJV

"""Writes the upper portion of a matrix in ijv format, one row for each edge,
separated by commas.  Only writes the upper triangular portion.
The result can be read from Matlab like this:

```
>> dl = dlmread('graph.txt');
>> a = sparse(dl(:,1),dl(:,2),dl(:,3));
>> n = max(size(a))
>> a(n,n) = 0;
>> a = a + a';
```
"""
function writeIJV(filename::AbstractString, mat)

  (ai,aj,av) = findnz(triu(mat))
  fh = open(filename,"w")
  for i in 1:length(ai)
    write(fh, "$(ai[i]),$(aj[i]),$(av[i])\n")
  end
  close(fh)

end #write IJV

" Writes a sparse matrix to a given text file "
function writeToFile{Tv,Ti}(filename::ASCIIString, a::SparseMatrixCSC{Tv,Ti})
    f = open(filename, "w")

    println(f, a.n, " ", a.m, " ", length(a.nzval))

    mat = tril(a);

    pos = 1
    for i in 1:length(mat.nzval)
        while mat.colptr[pos + 1] <= i
            pos = pos + 1
        end

        println(f, mat.rowval[i], " ", pos, " ", mat.nzval[i])
    end

    close(f)
end

" Reads a spare matrix from a given text file "
function readFromFile(filename::ASCIIString)
    r = readdlm(filename, ' ')

    m::Int64 = r[1,1]
    n::Int64 = r[1,2]
    nz::Int64 = r[1,3]
    U = Int64[]
    V = Int64[]
    W = Float64[]

    for i in 2:(nz+1)
        push!(U, r[i,1])
        push!(V, r[i,2])
        push!(W, r[i,3])

        push!(U, r[i,2])
        push!(V, r[i,1])
        push!(W, r[i,3])
    end

    return sparse(U,V,W,m,n)
end


text:
@theme TracTheme Dict(
    :name => "Trac",
    :description => "Based on the default trac highlighter design.",
    :comments => "Based on the theme from Pygments.",
    :tokens => Dict(
        COMMENT             => S"italic; fg: 999988",
        COMMENT_PREPROC     => S"bold; fg: 999999",
        COMMENT_SPECIAL     => S"bold; fg: 999999",

        ERROR               => S"bg: e3d2d2; fg: a61717",

        GENERIC_DELETED     => S"bg: ffdddd; fg: 000000",
        GENERIC_EMPH        => S"italic",
        GENERIC_ERROR       => S"fg: aa0000",
        GENERIC_HEADING     => S"fg: 999999",
        GENERIC_INSERTED    => S"bg: ddffdd; fg: 000000",
        GENERIC_OUTPUT      => S"fg: 888888",
        GENERIC_PROMPT      => S"fg: 555555",
        GENERIC_STRONG      => S"bold",
        GENERIC_SUBHEADING  => S"fg: aaaaaa",
        GENERIC_TRACEBACK   => S"fg: aa0000",

        KEYWORD             => S"bold",
        KEYWORD_TYPE        => S"fg: 445588",

        NAME_ATTRIBUTE      => S"fg: 008080",
        NAME_BUILTIN        => S"fg: 999999",
        NAME_CLASS          => S"bold; fg: 445588",
        NAME_CONSTANT       => S"fg: 008080",
        NAME_ENTITY         => S"fg: 800080",
        NAME_EXCEPTION      => S"bold; fg: 990000",
        NAME_FUNCTION       => S"bold; fg: 990000",
        NAME_NAMESPACE      => S"fg: 555555",
        NAME_TAG            => S"fg: 000080",
        NAME_VARIABLE       => S"fg: 008080",

        NUMBER              => S"fg: 009999",

        OPERATOR            => S"bold",

        STRING              => S"fg: bb8844",
        STRING_REGEX        => S"fg: 808000",

        TEXT                => S"",

        WHITESPACE          => S"fg: bbbbbb",
    )
)


text:
include("features.jl")
include("feature_extractors.jl")
include("lidar_sensors.jl")


text:
using BelgianDataset
using Test

@testset "BelgianDataset.jl" begin
    # join_partecipant_common_and_extra_datasets("resources/2010_Willem_BELGIUM_participant_common.csv","resources/2010_Willem_BELGIUM_participant_extra.csv")
    analyze_contact_data("resources/2010_Willem_BELGIUM_participant_common.csv","resources/2010_Willem_BELGIUM_participant_extra.csv","resources/2010_Willem_BELGIUM_contact_common.csv","resources/2010_Willem_BELGIUM_contact_extra.csv")
end

max_stars_repo_name: Luigi-Crisci/Hypergraphs-Epidemics

max_issues_repo_name: Luigi-Crisci/Hypergraphs-Epidemics

max_forks_repo_name: Luigi-Crisci/Hypergraphs-Epidemics