text:
using .KernelAbstractions
using .KernelAbstractions.CUDA
CUDA.allowscalar(false)

export compute_ba_J_cuda, compute_ba_J_multithread

println("NiBundleAdjustment: You find CUDA!")

@kernel function ba_kernel(reproj_err_d, obs, cams, X, w, feats)
    i = @index(Global)
    idx = (2*(i-1))+1
    @inbounds j = obs[1,i]
    @inbounds l = obs[2,i]
    @inbounds cam = cams[j]
    out! = zero(cam.x0)
    diff = zero(cam.x0)
    @inbounds res = compute_reproj_err(out!, diff, cam, X[l], w[i], feats[i])
    _, _, gcam, gX, gw, _ = (~compute_reproj_err)(P2(GVar(res[1].x, 1.0), GVar(res[1].y, 0.0)),
                    GVar(res[2]), GVar(res[3]), GVar(res[4]), GVar(res[5]), GVar(res[6]))
    NiBundleAdjustment.packgrad!(reproj_err_d, idx, gcam, gX, gw)
    _, _, gcam, gX, gw, _ = (~compute_reproj_err)(P2(GVar(res[1].x, 0.0), GVar(res[1].y, 1.0)),
                    GVar(res[2]), GVar(res[3]), GVar(res[4]), GVar(res[5]), GVar(res[6]))
    NiBundleAdjustment.packgrad!(reproj_err_d, idx+1, gcam, gX, gw)
end

@kernel function ba_kernel_w(w_err_d, w)
    i = @index(Global)
    @inbounds w_err_d[i] = grad(Grad(NiBundleAdjustment.compute_w_err)(Val(1), 0.0, w[i])[3])
end

function compute_ba_J_multithread(::Val{:NiLang}, cams::AbstractArray{<:Camera{T}}, X::AbstractArray{<:P3}, w, obs, feats::AbstractArray{<:P2}; nthread=4) where T
    p = size(obs,2)
    reproj_err_d = zeros(2*p, 15)
    event1 = ba_kernel(CPU(), nthread)(reproj_err_d, obs, cams, X, w, feats; ndrange=p)
    wait(event1)
    w_err_d = zeros(1,p)
    event2 = ba_kernel_w(CPU(), nthread)(w_err_d, w; ndrange=p)
    wait(event2)
    (reproj_err_d, w_err_d)
end

function compute_ba_J_cuda(::Val{:NiLang}, cams::CuArray{<:Camera{T}}, X::CuArray{<:P3}, w::CuArray, obs::CuArray, feats::CuArray{<:P2}; blockdim=256) where T
    p = size(obs,2)
    reproj_err_d = CUDA.zeros(T, 2*p, 15)
    event1 = ba_kernel(CUDADevice(), blockdim)(reproj_err_d, obs, cams, X, w, feats; ndrange=p)
    wait(event1)
    w_err_d = CUDA.zeros(T, 1, p)
    event2 = ba_kernel_w(CUDADevice(), blockdim)(w_err_d, w; ndrange=p)
    wait(event2)
    (reproj_err_d, w_err_d)
end


text:
@testset "$(rpad("Jacobian",80))" begin

    n = 1
    T = Float64

    x = [T(π),]
    j = reshape(2x, 1, 1)

    function F(x::Vector, b::Vector)
        b[:] = x.^2
    end

    function J(x::Vector, A::Matrix)
        A[:,:] = 2x
    end

    set_config(:jacobian_autodiff, true)
    JPAD = getJacobianParameters(nothing, F, T, n)

    set_config(:jacobian_autodiff, false)
    JPFD = getJacobianParameters(nothing, F, T, n)

    JPUS = getJacobianParameters(J, F, T, n)

    @test typeof(JPAD) <: JacobianParametersAD
    @test typeof(JPFD) <: JacobianParametersFD
    @test typeof(JPUS) <: JacobianParametersUser

    function test_jac(j1, j2, atol)
        for i in eachindex(j1,j2)
            @test j1[i] ≈ j2[i] atol=atol
        end
    end

    jad = zero(j)
    jfd = zero(j)
    jus = zero(j)

    computeJacobian(x, jad, JPAD)
    computeJacobian(x, jfd, JPFD)
    computeJacobian(x, jus, JPUS)

    test_jac(jad, j, eps())
    test_jac(jfd, j, 1E-7)
    test_jac(jus, j, 0)

    jad2 = zero(j)
    jfd2 = zero(j)

    computeJacobianAD(x, jad2, F)
    computeJacobianFD(x, jfd2, F, get_config(:jacobian_fd_ϵ))

    test_jac(jad, jad2, 0)
    test_jac(jfd, jfd2, 0)

end

ext: jl


text:
#!/cvmfs/soft.computecanada.ca/easybuild/software/2020/avx2/Core/julia/1.5.2/bin/julia
#SBATCH --mail-user=mmcleod2@ualberta.ca
#SBATCH --mail-type=ALL
#SBATCH -o onedtmaze_gpi.out # Standard output
#SBATCH -e onedtmaze_gpi.err # Standard error
#SBATCH --mem-per-cpu=3000M # Memory request of 3GB
#SBATCH --time=03:30:00 #105 hours total
#SBATCH --ntasks=32
#SBATCH --account=def-amw8

using Pkg
Pkg.activate(".")

include(joinpath(ENV["SLURM_SUBMIT_DIR"], "parallel/parallel_config.jl"))
reproduce_config_experiment("configs/revised_onedtmaze/GPI_behaviour.toml", save_path="~/scratch/curiosity")

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
"""
Validation PSSE/MultiGen:
This case study defines a four bus system multimachine with several units in one bus.
The fault trips a generator
"""

##################################################
############### SOLVE PROBLEM ####################
##################################################

raw_file = joinpath(TEST_FILES_DIR, "benchmarks/psse/MultiGen/FourBusMulti.raw")
dyr_file = joinpath(TEST_FILES_DIR, "benchmarks/psse/MultiGen/FourBus_multigen.dyr")
gen_trip_csv_file =
    joinpath(TEST_FILES_DIR, "benchmarks/psse/MultiGen/gen_trip_results.csv")
line_trip_csv_file =
    joinpath(TEST_FILES_DIR, "benchmarks/psse/MultiGen/line_trip_results.csv")

@testset "Test 35 MultiGen ResidualModel" begin
    path = (joinpath(pwd(), "test-multigen"))
    !isdir(path) && mkdir(path)
    try
        sys = System(raw_file, dyr_file)
        for l in get_components(PSY.PowerLoad, sys)
            PSY.set_model!(l, PSY.LoadModels.ConstantImpedance)
        end

        dc = get_dynamic_injector(get_component(ThermalStandard, sys, "generator-102-SW"))
        # Define Simulation Problem
        sim = Simulation!(
            ResidualModel,
            sys, #system
            path,
            (0.0, 20.0), #time span
            GeneratorTrip(1.0, dc),
        )

        # Solve problem
        @test execute!(sim, IDA(), dtmax = 0.005, saveat = 0.005) ==
              PSID.SIMULATION_FINALIZED
        results = read_results(sim)

        # # Obtain results
        M = get_csv_data(gen_trip_csv_file)
        t_psse = M[:, 1]
        for (ix, b) in enumerate(101:104)
            t, Vt = get_voltage_magnitude_series(results, b)
            Vt_psse = M[:, 47 + 2 * (ix - 1)]
            @test LinearAlgebra.norm(t - round.(t_psse, digits = 3)) == 0.0
            @test LinearAlgebra.norm(Vt - Vt_psse, Inf) <= 1e-3
        end

        for (ix, g) in enumerate(get_components(ThermalStandard, sys))
            if isa(get_dynamic_injector(g), PSY.DynamicInverter)
                continue
            end
            gen_name = get_name(g)
            t, ω = get_state_series(results, (gen_name, :ω))
            ω_psse = M[:, 5 + 5 * (ix - 1)] .+ 1.0
            @test LinearAlgebra.norm(t - round.(t_psse, digits = 3)) == 0.0
            @test LinearAlgebra.norm(ω - ω_psse, Inf) <= 0.1 # relaxed due to inconsistent PSSE outputs
        end

    finally
        @info("removing test files")
        rm(path, force = true, recursive = true)
    end
end

@testset "Test 35 MultiGen MassMatrixModel" begin
    path = (joinpath(pwd(), "test-psse-gast"))
    !isdir(path) && mkdir(path)
    try
        sys = System(raw_file, dyr_file)
        for l in get_components(PSY.PowerLoad, sys)
            PSY.set_model!(l, PSY.LoadModels.ConstantImpedance)
        end

        dc = get_dynamic_injector(get_component(ThermalStandard, sys, "generator-102-SW"))
        # Define Simulation Problem
        sim = Simulation!(
            MassMatrixModel,
            sys, #system
            path,
            (0.0, 20.0), #time span
            GeneratorTrip(1.0, dc),
        )

        # Solve problem
        @test execute!(sim, Rodas4(), dtmax = 0.005, saveat = 0.005) ==
              PSID.SIMULATION_FINALIZED
        results = read_results(sim)

        M = get_csv_data(gen_trip_csv_file)
        t_psse = M[:, 1]
        for (ix, b) in enumerate(101:104)
            t, Vt = get_voltage_magnitude_series(results, b)
            Vt_psse = M[:, 47 + 2 * (ix - 1)]
            @test LinearAlgebra.norm(t - round.(t_psse, digits = 3)) == 0.0
            @test LinearAlgebra.norm(Vt - Vt_psse, Inf) <= 1e-3
        end

        for (ix, g) in enumerate(get_components(ThermalStandard, sys))
            if isa(get_dynamic_injector(g), PSY.DynamicInverter)
                continue
            end
            gen_name = get_name(g)
            t, ω = get_state_series(results, (gen_name, :ω))
            ω_psse = M[:, 5 + 5 * (ix - 1)] .+ 1.0
            @test LinearAlgebra.norm(t - round.(t_psse, digits = 3)) == 0.0
            @test LinearAlgebra.norm(ω - ω_psse, Inf) <= 0.1 # relaxed due to inconsistent PSSE outputs
        end
    finally
        @info("removing test files")
        rm(path, force = true, recursive = true)
    end
end


text:
export roots

function _roots(f, region, contractor, tol = default_tolerance, target_task_number = nothing)
    search = ThreadedRootSearch(f, region, contractor, tol, target_task_number)
    return last(search)
end

"""
    roots(f, region, method = Krawczyk(), tol = 1e-7, maxtasks = nothing)

Find all isolated roots of a function `f:R^n → R^n` in a region `region`, if the
number of roots is finite.

Inputs:
  - `f`: function whose roots will be found
  - `region`: region in which roots are searched given by an interval or a
    vector of intervals for multidimensional problems
  - `method`: function that, when applied to the function `f`, determines
    the status of a given box `region`. It returns the new box and a symbol
    indicating the status. Current possible values are `Bisection()`, `Newton()`
    and `Krawczyk()`
  - `tol`: Absolute tolerance. If a region smaller than `tol` cannot be proven
    to contain or not contain a root, it is returned as indeterminate.
  - `maxtasks`: Limit the number of tasks spawned concurrently. By default set
    to `10^4 * nthreads()`.
"""
function roots(
    f,
    region,
    method::Union{Krawczyk,Newton} = default_contractor,
    tol = default_tolerance,
    target_task_number = nothing
)
    contractor = GradientContractor(f, method, region)
    return _roots(f, region, contractor, tol, target_task_number)
end

function roots(
    f,
    region,
    method::Bisection,
    tol = default_tolerance,
    target_task_number = nothing
)
    return _roots(f, region, TrivialContractor(), tol, target_task_number)
end


text:
function _append(file, txt)
    open(file, "a") do io
        print(io, txt)
    end
end


text:
using BitInformation
using JLD2, NetCDF

path = "/network/aopp/chaos/pred/kloewer/esowc/cams/"
filelist = filter(x->endswith(x,"_v3.nc"),readdir(path))

n = length(filelist)
varnames = fill("",n)
nbits = 32

IC = fill(0.0,n,nbits)

for (i,file) in enumerate(filelist)
    varname = split(split(file,"cams_")[end],"_2019")[1]
    varnames[i] = varname
    print("$varname")
    ncfile = NetCDF.open(joinpath(path,file))

    # find the variable name by size
    var = [var for var in ncfile.vars if prod(size(var[2])) == 900*451*137][1][1]
    X = ncfile.vars[var][:,6:end-5,:]  # exclude north/southpole

    # convert biased exponent to signed exponent
    BitInformation.signed_exponent!(X)

    IC[i,:] = bitinformation(X,:all_dimensions)
    print("X, ")

    @save joinpath(path,"..","analysis/bitinformation_all_filt.jld2") varnames IC
end


text:
module IteratorInterfaceExtensions

export getiterator, isiterable, IteratorSize2

isiterable(x::T) where {T} = Base.isiterable(T)

function getiterator(x)
    if !isiterable(x)
        error("Can't get iterator for non iterable source.")
    end
    return x
end

struct HasLengthAfterStart <: Base.IteratorSize end

IteratorSize2(x) = IteratorSize2(typeof(x))
IteratorSize2(::Type{T}) where {T} = Base.IteratorSize(T)

end # module


text:
# LLVM operations and intrinsics
module Intrinsics

# Note, that in the functions below, some care needs to be taken when passing
# Julia Bools to LLVM. Julia passes Bools as LLVM i8 but expect them to only
# have the last bit as non-zero. Failure to comply with this can give weird errors
# like false !== false where the first false is the result of some computation.

# Note, no difference is made between Julia usigned integers and signed integers
# when passed to LLVM. It is up to the caller to make sure that the correct
# intrinsic is called (e.g uitofp vs sitofp).

import ..SIMD: SIMD, VE, LVec, FloatingTypes
# Inlcude Bool in IntegerTypes
const IntegerTypes = Union{SIMD.IntegerTypes, Bool}

const d = Dict{DataType, String}(
    Bool         => "i8",
    Int8         => "i8",
    Int16        => "i16",
    Int32        => "i32",
    Int64        => "i64",
    Int128       => "i128",

    UInt8        => "i8",
    UInt16       => "i16",
    UInt32       => "i32",
    UInt64       => "i64",
    UInt128      => "i128",

    #Float16     => "half",
    Float32      => "float",
    Float64      => "double",
)
# Add the Ptr translations
foreach(x -> (d[Ptr{x}] = d[Int]), collect(keys(d)))

# LT = LLVM Type (scalar and vectors), we keep type names intentionally short
# to make the signatures smaller
const LT{T} = Union{LVec{<:Any, T}, T}

suffix(N::Integer, ::Type{Ptr{T}}) where {T} = "v$(N)p0$(T<:IntegerTypes ? "i" : "f")$(8*sizeof(T))"
suffix(N::Integer, ::Type{T}) where {T}      = "v$(N)$(T<:IntegerTypes   ? "i" : "f")$(8*sizeof(T))"
suffix(::Type{T}) where {T}                  = "$(T<:IntegerTypes        ? "i" : "f")$(8*sizeof(T))"

dotit(f) = replace(string(f), "_" => ".")
llvm_name(llvmf, N, T)                           = string("llvm", ".", dotit(llvmf), ".", suffix(N, T))
llvm_name(llvmf, ::Type{LVec{N, T}}) where {N,T} = string("llvm", ".", dotit(llvmf), ".", suffix(N, T))
llvm_name(llvmf, ::Type{T}) where {T}            = string("llvm", ".", dotit(llvmf), ".", suffix(T))

llvm_type(::Type{T}) where {T}            = d[T]
llvm_type(::Type{LVec{N, T}}) where {N,T} = "< $N x $(d[T])>"

############
# FastMath #
############

module FastMath
    const nnan     = 1 << 0
    const ninf     = 1 << 1
    const nsz      = 1 << 2
    const arcp     = 1 << 3
    const contract = 1 << 4
    const afn      = 1 << 5
    const reassoc  = 1 << 6
    const fast     = 1 << 7
end

struct FastMathFlags{T} end
Base.@pure FastMathFlags(T::Int) = FastMathFlags{T}()

function fp_str(::Type{FastMathFlags{T}}) where {T}
    flags = String[]
    (T & FastMath.nnan     != 0) && push!(flags, "nnan")
    (T & FastMath.ninf     != 0) && push!(flags, "ninf")
    (T & FastMath.nsz      != 0) && push!(flags, "nsz")
    (T & FastMath.arcp     != 0) && push!(flags, "arcp")
    (T & FastMath.contract != 0) && push!(flags, "contract")
    (T & FastMath.afn      != 0) && push!(flags, "afn")
    (T & FastMath.reassoc  != 0) && push!(flags, "reassoc")
    (T & FastMath.fast     != 0) && push!(flags, "fast")
    return join(flags, " ")
end
fp_str(::Type{Nothing}) = ""

const FPFlags{T} = Union{Nothing, FastMathFlags{T}}

####################
# Unary operators  #
####################

const UNARY_INTRINSICS_FLOAT = [
    :sqrt
    :sin
    :cos
    :exp
    :trunc
    :exp2
    :log
    :log10
    :log2
    :fabs
    :floor
    :ceil
    :rint
    :nearbyint
    :round
]

const UNARY_INTRINSICS_INT = [
    :bitreverse
    :bswap
    :ctpop
    :ctlz
    :cttz
    :fshl
    :fshr
]
for (fs, c) in zip([UNARY_INTRINSICS_FLOAT, UNARY_INTRINSICS_INT],
                   [FloatingTypes,          IntegerTypes])
    for f in fs
        @eval begin
            @generated function $(f)(x::T) where T<:LT{<:$c}
                ff = llvm_name($(QuoteNode(f)), T)
                return :(
                    $(Expr(:meta, :inline));
                    ccall($ff, llvmcall, T, (T,), x)
                )
            end
        end
    end
end

# fneg (not an intrinsic so cannot use `ccall)
@generated function fneg(x::T, ::F=nothing) where {T<:LT{<:FloatingTypes}, F<:FPFlags}
    fpflags = fp_str(F)
    s = """
    %2 = fneg $fpflags $(llvm_type(T)) %0
    ret $(llvm_type(T)) %2
    """
    return :(
        $(Expr(:meta, :inline));
        Base.llvmcall($s, T, Tuple{T}, x)
    )
end

#####################
# Binary operators  #
#####################

const BINARY_OPS_FLOAT = [
    :fadd
    :fsub
    :fmul
    :fdiv
    :frem
]

const BINARY_OPS_INT = [
    :add
    :sub
    :mul
    :sdiv
    :udiv
    :srem
    :urem
    :shl
    :ashr
    :lshr
    :and
    :or
    :xor
]

for f in BINARY_OPS_FLOAT
    @eval @generated function $f(x::T, y::T, ::F=nothing) where {T<:LT{<:FloatingTypes}, F<:FPFlags}
        fpflags = fp_str(F)
        ff = $(QuoteNode(f))
        s = """
        %3 = $ff $fpflags $(llvm_type(T)) %0, %1
        ret $(llvm_type(T)) %3
        """
        return :(
            $(Expr(:meta, :inline));
            Base.llvmcall($s, T, Tuple{T, T}, x, y)
        )
    end
end

for f in BINARY_OPS_INT
    @eval @generated function $f(x::T, y::T) where T<:LT{<:IntegerTypes}
        ff = $(QuoteNode(f))
        s = """
        %3 = $ff $(llvm_type(T)) %0, %1
        ret $(llvm_type(T)) %3
        """
        return :(
            $(Expr(:meta, :inline));
            Base.llvmcall($s, T, Tuple{T, T}, x, y)
        )
    end
end

const BINARY_INTRINSICS_FLOAT = [
    :minnum
    :maxnum
    :minimum
    :maximum
    :copysign
    :pow
    :floor
    :ceil
    :trunc
    :rint
    :nearbyint
    :round
]

const BINARY_INTRINSICS_INT = [
    :sadd_sat
    :uadd_sat
    :ssub_sat
    :usub_sat
]

for (fs, c) in zip([BINARY_INTRINSICS_FLOAT, BINARY_INTRINSICS_INT],
                   [FloatingTypes,           IntegerTypes])
    for f in fs
        @eval @generated function $(f)(x::T, y::T) where T<:LT{<:$c}
            ff = llvm_name($(QuoteNode(f)), T,)
            return :(
                $(Expr(:meta, :inline));
                ccall($ff, llvmcall, T, (T, T), x, y)
            )
        end
    end
end

# pow, powi
for (f, c) in [(:pow, FloatingTypes), (:powi, IntegerTypes)]
    @eval @generated function $(f)(x::T, y::T2) where {T <: LT{<:FloatingTypes}, T2 <: $c}
        ff = llvm_name($(QuoteNode(f)), T)
        return :(
            $(Expr(:meta, :inline));
            ccall($ff, llvmcall, T, (T, T2), x, y)
        )
    end
end

# Overflow
const OVERFLOW_INTRINSICS = [
    :sadd_with_overflow
    :uadd_with_overflow
    :ssub_with_overflow
    :usub_with_overflow
    :smul_with_overflow
    :umul_with_overflow
]

const SUPPORTS_VEC_OVERFLOW = Base.libllvm_version >= v"9"
for f in OVERFLOW_INTRINSICS
    @eval @generated function $f(x::LVec{N, T}, y::LVec{N, T}) where {N, T <: IntegerTypes}
        if !SUPPORTS_VEC_OVERFLOW
            return :(error("LLVM version 9.0 or greater required (Julia 1.5 or greater)"))
        end
        ff = llvm_name($(QuoteNode(f)), N, T)
        if $(QuoteNode(f)) == :smul_with_overflow && Sys.ARCH == :i686 && T == Int64
            str = "this intrinsic ($ff) is broken on i686"
            return :(error($str))
        end
        decl = "declare {<$N x $(d[T])>, <$N x i1>} @$ff(<$N x $(d[T])>, <$N x $(d[T])>)"

        # Julia passes Tuple{[U]Int8, Bool} as [2 x i8] so we need to special case that scenario
        ret_type = sizeof(T) == 1 ? "[2 x <$N x i8>]" : "{<$N x $(d[T])>, <$N x i8>}"

        s = """
        %res = call {<$N x $(d[T])>, <$N x i1>} @$ff(<$N x $(d[T])> %0, <$N x $(d[T])> %1)
        %plus     = extractvalue {<$N x $(d[T])>, <$N x i1>} %res, 0
        %overflow = extractvalue {<$N x $(d[T])>, <$N x i1>} %res, 1
        %overflow_ext = zext <$(N) x i1> %overflow to <$(N) x i8>
        %new_tuple   = insertvalue $ret_type undef,      <$N x $(d[T])> %plus,         0
        %new_tuple_2 = insertvalue $ret_type %new_tuple, <$N x i8>      %overflow_ext, 1
        ret $ret_type %new_tuple_2
        """
        return :(
            $(Expr(:meta, :inline));
            Base.llvmcall(($decl, $s), Tuple{LVec{N, T}, LVec{N, Bool}}, Tuple{LVec{N, T}, LVec{N, T}}, x, y)
        )
    end
end

# Comparisons
const CMP_FLAGS_FLOAT = [
    :false
    :oeq
    :ogt
    :oge
    :olt
    :ole
    :one
    :ord
    :ueq
    :ugt
    :uge
    :ult
    :ule
    :une
    :uno
    :true
]

const CMP_FLAGS_INT = [
    :eq
    :ne
    :sgt
    :sge
    :slt
    :sle
    :ugt
    :uge
    :ult
    :ule
]

for flag in CMP_FLAGS_FLOAT
    ftot = Symbol(string("fcmp_", flag))
    @eval @generated function $ftot(x::LVec{N, T}, y::LVec{N, T}, ::F=nothing) where {N, T <: FloatingTypes, F<:FPFlags}
        fpflags = fp_str(F)
        fflag = $(QuoteNode(flag))
        s = """
        %res = fcmp $(fpflags) $(fflag) <$(N) x $(d[T])> %0, %1
        %resb = zext <$(N) x i1> %res to <$(N) x i8>
        ret <$(N) x i8> %resb
        """
        return :(
            $(Expr(:meta, :inline));
            Base.llvmcall($s, LVec{N, Bool}, Tuple{LVec{N, T}, LVec{N, T}}, x, y)
        )
    end
end

for flag in CMP_FLAGS_INT
    ftot = Symbol(string("icmp_", flag))
    @eval @generated function $ftot(x::LVec{N, T}, y::LVec{N, T}) where {N, T <: IntegerTypes}
        fflag = $(QuoteNode(flag))
        s = """
        %res = icmp $(fflag) <$(N) x $(d[T])> %0, %1
        %resb = zext <$(N) x i1> %res to <$(N) x i8>
        ret <$(N) x i8> %resb
        """
        return :(
            $(Expr(:meta, :inline));
            Base.llvmcall($s, LVec{N, Bool}, Tuple{LVec{N, T}, LVec{N, T}}, x, y)
        )
    end
end

#####################
# Ternary operators #
#####################

@generated function select(cond::LVec{N, Bool}, x::LVec{N, T}, y::LVec{N, T}) where {N, T}
    s = """
    %cond = trunc <$(N) x i8> %0 to <$(N) x i1>
    %res = select <$N x i1> %cond, <$N x $(d[T])> %1, <$N x $(d[T])> %2
    ret <$N x $(d[T])> %res
    """
    return :(
        $(Expr(:meta, :inline));
        Base.llvmcall($s, LVec{N, T}, Tuple{LVec{N, Bool}, LVec{N, T}, LVec{N, T}}, cond, x, y)
    )
end

const MULADD_INTRINSICS = [
    :fmuladd,
    :fma,
]

for f in MULADD_INTRINSICS
    @eval @generated function $(f)(a::LVec{N, T}, b::LVec{N, T}, c::LVec{N, T}) where {N, T<:FloatingTypes}
        ff = llvm_name($(QuoteNode(f)), N, T)
        return :(
            $(Expr(:meta, :inline));
            ccall($ff, llvmcall, LVec{N, T}, (LVec{N, T}, LVec{N, T}, LVec{N, T}), a, b, c)
        )
    end
end

################
# Load / store #
################

# These alignment numbers feels a bit dubious
n_align(align, N, T) = align ? N * sizeof(T) : sizeof(T)
temporal_str(temporal) = temporal ? ", !nontemporal !{i32 1}" : ""

@generated function load(x::Type{LVec{N, T}}, ptr::Ptr{T},
                         ::Val{Al}=Val(false), ::Val{Te}=Val(false)) where {N, T, Al, Te}
    s = """
    %ptr = inttoptr $(d[Int]) %0 to <$N x $(d[T])>*
    %res = load <$N x $(d[T])>, <$N x $(d[T])>* %ptr, align $(n_align(Al, N, T)) $(temporal_str(Te))
    ret <$N x $(d[T])> %res
    """
    return :(
        $(Expr(:meta, :inline));
        Base.llvmcall($s, LVec{N, T}, Tuple{Ptr{T}}, ptr)
    )
end

@generated function maskedload(ptr::Ptr{T}, mask::LVec{N,Bool},
                               ::Val{Al}=Val(false), ::Val{Te}=Val(false)) where {N, T, Al, Te}
    # TODO: Allow setting the passthru
    decl = "declare <$N x $(d[T])> @llvm.masked.load.$(suffix(N, T))(<$N x $(d[T])>*, i32, <$N x i1>, <$N x $(d[T])>)"
    s = """
    %mask = trunc <$(N) x i8> %1 to <$(N) x i1>
    %ptr = inttoptr $(d[Int]) %0 to <$N x $(d[T])>*
    %res = call <$N x $(d[T])> @llvm.masked.load.$(suffix(N, T))(<$N x $(d[T])>* %ptr, i32 $(n_align(Al, N, T)), <$N x i1> %mask, <$N x $(d[T])> zeroinitializer)
    ret <$N x $(d[T])> %res
    """
    return :(
        $(Expr(:meta, :inline));
        Base.llvmcall(($decl, $s), LVec{N, T}, Tuple{Ptr{T}, LVec{N,Bool}}, ptr, mask)
    )
end

@generated function store(x::LVec{N, T}, ptr::Ptr{T},
                          ::Val{Al}=Val(false), ::Val{Te}=Val(false)) where {N, T, Al, Te}
    s = """
    %ptr = inttoptr $(d[Int]) %1 to <$N x $(d[T])>*
    store <$N x $(d[T])> %0, <$N x $(d[T])>* %ptr, align $(n_align(Al, N, T)) $(temporal_str(Te))
    ret void
    """
    return :(

        $(Expr(:meta, :inline));
        Base.llvmcall($s, Cvoid, Tuple{LVec{N, T}, Ptr{T}}, x, ptr)
    )
end

@generated function maskedstore(x::LVec{N, T}, ptr::Ptr{T}, mask::LVec{N,Bool},
                               ::Val{Al}=Val(false), ::Val{Te}=Val(false)) where {N, T, Al, Te}
    # TODO: Allow setting the passthru
    decl = "declare <$N x $(d[T])> @llvm.masked.store.$(suffix(N, T))(<$N x $(d[T])>, <$N x $(d[T])>*, i32, <$N x i1>)"
    s = """
    %mask = trunc <$(N) x i8> %2 to <$(N) x i1>
    %ptr = inttoptr $(d[Int]) %1 to <$N x $(d[T])>*
    %res = call <$N x $(d[T])> @llvm.masked.store.$(suffix(N, T))(<$N x $(d[T])> %0, <$N x $(d[T])>* %ptr, i32 $(n_align(Al, N, T)), <$N x i1> %mask)
    ret void
    """
    return :(
        $(Expr(:meta, :inline));
        Base.llvmcall(($decl, $s), Cvoid, Tuple{LVec{N, T}, Ptr{T}, LVec{N,Bool}}, x, ptr, mask)
    )
end

@generated function maskedexpandload(ptr::Ptr{T}, mask::LVec{N,Bool}) where {N, T}
    # TODO: Allow setting the passthru
    decl = "declare <$N x $(d[T])> @llvm.masked.expandload.$(suffix(N, T))($(d[T])*, <$N x i1>, <$N x $(d[T])>)"
    s = """
    %mask = trunc <$(N) x i8> %1 to <$(N) x i1>
    %ptr = inttoptr $(d[Int]) %0 to $(d[T])*
    %res = call <$N x $(d[T])> @llvm.masked.expandload.$(suffix(N, T))($(d[T])* %ptr, <$N x i1> %mask, <$N x $(d[T])> zeroinitializer)
    ret <$N x $(d[T])> %res
    """
    return :(
        $(Expr(:meta, :inline));
        Base.llvmcall(($decl, $s), LVec{N, T}, Tuple{Ptr{T}, LVec{N, Bool}}, ptr, mask)
    )
end

@generated function maskedcompressstore(x::LVec{N, T}, ptr::Ptr{T},
                                        mask::LVec{N,Bool}) where {N, T}
    decl = "declare <$N x $(d[T])> @llvm.masked.compressstore.$(suffix(N, T))(<$N x $(d[T])>, $(d[T])*, <$N x i1>)"
    s = """
    %mask = trunc <$(N) x i8> %2 to <$(N) x i1>
    %ptr = inttoptr $(d[Int]) %1 to $(d[T])*
    call <$N x $(d[T])> @llvm.masked.compressstore.$(suffix(N, T))(<$N x $(d[T])> %0, $(d[T])* %ptr, <$N x i1> %mask)
    ret void
    """
    return :(
        $(Expr(:meta, :inline));
        Base.llvmcall(($decl, $s), Cvoid, Tuple{LVec{N, T}, Ptr{T}, LVec{N, Bool}}, x, ptr, mask)
    )
end

####################
# Gather / Scatter #
####################

@generated function maskedgather(ptrs::LVec{N,Ptr{T}},
                                 mask::LVec{N,Bool}, ::Val{Al}=Val(false)) where {N, T, Al}
    # TODO: Allow setting the passthru
    decl = "declare <$N x $(d[T])> @llvm.masked.gather.$(suffix(N, T))(<$N x $(d[T])*>, i32, <$N x i1>, <$N x $(d[T])>)"
    s = """
    %mask = trunc <$(N) x i8> %1 to <$(N) x i1>
    %ptrs = inttoptr <$N x $(d[Int])> %0 to <$N x $(d[T])*>
    %res = call <$N x $(d[T])> @llvm.masked.gather.$(suffix(N, T))(<$N x $(d[T])*> %ptrs, i32 $(n_align(Al, N, T)), <$N x i1> %mask, <$N x $(d[T])> zeroinitializer)
    ret <$N x $(d[T])> %res
    """
    return :(
        $(Expr(:meta, :inline));
        Base.llvmcall(($decl, $s), LVec{N, T}, Tuple{LVec{N, Ptr{T}}, LVec{N, Bool}}, ptrs, mask)
    )
end

@generated function maskedscatter(x::LVec{N, T}, ptrs::LVec{N, Ptr{T}},
                                  mask::LVec{N,Bool}, ::Val{Al}=Val(false)) where {N, T, Al}

    decl = "declare <$N x $(d[T])> @llvm.masked.scatter.$(suffix(N, T))(<$N x $(d[T])>, <$N x $(d[T])*>, i32, <$N x i1>)"
    s = """
    %mask = trunc <$(N) x i8> %2 to <$(N) x i1>
    %ptrs = inttoptr <$N x $(d[Int])> %1 to <$N x $(d[T])*>
    call <$N x $(d[T])> @llvm.masked.scatter.$(suffix(N, T))(<$N x $(d[T])> %0, <$N x $(d[T])*> %ptrs, i32 $(n_align(Al, N, T)), <$N x i1> %mask)
    ret void
    """
    return :(
        $(Expr(:meta, :inline));
        Base.llvmcall(($decl, $s), Cvoid, Tuple{LVec{N, T}, LVec{N, Ptr{T}}, LVec{N, Bool}}, x, ptrs, mask)
    )
end

######################
# LVector Operations #
######################

@generated function extractelement(x::LVec{N, T}, i::I) where {N, T, I <: IntegerTypes}
    s = """
    %3 = extractelement <$N x $(d[T])> %0, $(d[I]) %1
    ret $(d[T]) %3
    """
    return :(
        $(Expr(:meta, :inline));
        Base.llvmcall($s, T, Tuple{LVec{N, T}, $i}, x, i)
    )
end

@generated function insertelement(x::LVec{N, T}, v::T, i::IntegerTypes) where {N, T}
    s = """
    %4 = insertelement <$N x $(d[T])> %0, $(d[T]) %1, $(d[i]) %2
    ret <$N x $(d[T])> %4
    """
    return :(
        $(Expr(:meta, :inline));
        Base.llvmcall($s, LVec{N, T}, Tuple{LVec{N, T}, T, typeof(i)}, x, v, i)
    )
end

_shuffle_vec(I) = join((string("i32 ", i == :undef ? "undef" : Int32(i::Integer)) for i in I), ", ")
@generated function shufflevector(x::LVec{N, T}, y::LVec{N, T}, ::Val{I}) where {N, T, I}
    shfl = _shuffle_vec(I)
    M = length(I)
    s = """
    %res = shufflevector <$N x $(d[T])> %0, <$N x $(d[T])> %1, <$M x i32> <$shfl>
    ret <$M x $(d[T])> %res
    """
    return :(
        $(Expr(:meta, :inline));
        Base.llvmcall($s, LVec{$M, T}, Tuple{LVec{N, T}, LVec{N, T}}, x, y)
    )
end

@generated function shufflevector(x::LVec{N, T}, ::Val{I}) where {N, T, I}
    shfl = _shuffle_vec(I)
    M = length(I)
    s = """
    %res = shufflevector <$(N) x $(d[T])> %0, <$N x $(d[T])> undef, <$M x i32> <$shfl>
    ret <$M x $(d[T])> %res
    """
    return :(
        $(Expr(:meta, :inline));
        Base.llvmcall($s, LVec{$M, T}, Tuple{LVec{N, T}}, x)
    )
end

@generated function constantvector(v::T, y::Type{LVec{N, T}}) where {N, T}
    s = """
    %2 = insertelement <$N x $(d[T])> undef, $(d[T]) %0, i32 0
    %res = shufflevector <$N x $(d[T])> %2, <$N x $(d[T])> undef, <$N x i32> zeroinitializer
    ret <$N x $(d[T])> %res
    """
    return :(
        $(Expr(:meta, :inline));
        Base.llvmcall($s, LVec{N, T}, Tuple{T}, v)
    )
end

#########################
# Conversion Operations #
#########################

const CAST_SIZE_CHANGE_FLOAT = [
    (:fptrunc, >)
    (:fpext, <)
]

const CAST_SIZE_CHANGE_INT = [
    (:trunc, >)
    (:zext, <)
    (:sext, <)
]

for (fs, c) in zip([CAST_SIZE_CHANGE_FLOAT, CAST_SIZE_CHANGE_INT],
                   [FloatingTypes,          IntegerTypes])
    for (f, criteria) in fs
        @eval @generated function $f(::Type{LVec{N, T2}}, x::LVec{N, T1}) where {N, T1 <: $c, T2 <: $c}
            sT1, sT2 = sizeof(T1) * 8, sizeof(T2) * 8
            # Not changing size is not allowed
            if !$criteria(sT1, sT2)
                str = "size of conversion type ($T2: $sT2) must be $($criteria) than the element type ($T1: $sT1)"
                return :(throw(ArgumentError($str)))
            end
            ff = $(QuoteNode(f))
            s = """
            %2 = $ff <$(N) x $(d[T1])> %0 to <$(N) x $(d[T2])>
            ret <$(N) x $(d[T2])> %2
            """
            return :(
                $(Expr(:meta, :inline));
                Base.llvmcall($s, LVec{N, T2}, Tuple{LVec{N, T1}}, x)
            )
        end
    end
end

const CONVERSION_FLOAT_TO_INT = [
    :fptoui,
    :fptosi
]

const CONVERSION_INT_TO_FLOAT = [
    :uitofp,
    :sitofp
]

for (fs, (from, to)) in zip([CONVERSION_FLOAT_TO_INT,       CONVERSION_INT_TO_FLOAT],
                           [(FloatingTypes, IntegerTypes), (IntegerTypes, FloatingTypes)])
    for f in fs
        @eval @generated function $f(::Type{LVec{N, T2}}, x::LVec{N, T1}) where {N, T1 <: $from, T2 <: $to}
            ff = $(QuoteNode(f))
            s = """
            %2 = $ff <$(N) x $(d[T1])> %0 to <$(N) x $(d[T2])>
            ret <$(N) x $(d[T2])> %2
            """
            return :(
                $(Expr(:meta, :inline));
                Base.llvmcall($s, LVec{N, T2}, Tuple{LVec{N, T1}}, x)
            )
        end
    end
end

###########
# Bitcast #
###########

@generated function bitcast(::Type{T1}, x::T2) where {T1<:LT, T2<:LT}
    sT1, sT2 = sizeof(T1) * 8, sizeof(T2) * 8
    if sT1 != sT2
        return :(throw(ArgumentError(("size of conversion type ($($T1): $($sT1)) must be equal to the vector type ($($T2): $($sT2))"))))
    end
    s = """
    %2 = bitcast $(llvm_type(T2)) %0 to $(llvm_type(T1))
    ret $(llvm_type(T1)) %2
    """
    return :(
        $(Expr(:meta, :inline));
        Base.llvmcall($s, T1, Tuple{T2}, x)
    )
end

##################################
# Horizontal reductions (LLVM 9) #
##################################

const HORZ_REDUCTION_OPS_FLOAT = [
    :fmax
    :fmin
]

const HORZ_REDUCTION_OPS_INT = [
    :and
    :or
    :mul
    :add
    :smax
    :umax
    :smin
    :umin
]

for (fs, c) in zip([HORZ_REDUCTION_OPS_FLOAT, HORZ_REDUCTION_OPS_INT],
                   [FloatingTypes,            IntegerTypes])
    for f in fs
        f_red = Symbol("reduce_", f)
        @eval @generated function $f_red(x::LVec{N, T}) where {N,T<:$c}
            ff = llvm_name(string("experimental.vector.reduce.", $(QuoteNode(f))), N, T)
            decl = "declare $(d[T]) @$ff(<$N x $(d[T])>)"
            s2 = """
            %res = call $(d[T]) @$ff(<$N x $(d[T])> %0)
            ret $(d[T]) %res
            """
            return quote
                $(Expr(:meta, :inline));
                Base.llvmcall($(decl, s2), T, Tuple{LVec{N, T},}, x)
            end
        end
    end
end

# The fadd and fmul reductions take an initializer
horz_reduction_version = Base.libllvm_version < v"9" ? "" : "v2."
for (f, neutral) in [(:fadd, "0.0"), (:fmul, "1.0")]
    f_red = Symbol("reduce_", f)
    @eval @generated function $f_red(x::LVec{N, T}) where {N,T<:FloatingTypes}
        ff = llvm_name(string("experimental.vector.reduce.$horz_reduction_version", $(QuoteNode(f))), N, T)
        decl = "declare $(d[T]) @$ff($(d[T]), <$N x $(d[T])>)"
        s2 = """
        %res = call $(d[T]) @$ff($(d[T]) $($neutral), <$N x $(d[T])> %0)
        ret $(d[T]) %res
        """
        return quote
            $(Expr(:meta, :inline));
            Base.llvmcall($(decl, s2), T, Tuple{LVec{N, T},}, x)
        end
    end
end

# See: https://llvm.org/docs/LangRef.html#id839
@generated function reduce_add(x::LVec{N, Bool}) where {N}
    native_bit_width = sizeof(Int) * 8
    if N < native_bit_width
        ret = """
        %res = zext i$(N) %maskipopcnt to i$(native_bit_width)
        ret i$(native_bit_width) %res
        """
    elseif N == native_bit_width
        ret = "ret i$(native_bit_width) %maskipopcnt"
    else
        ret = """
        %res = trunc i$(N) %maskipopcnt to i$(native_bit_width)
        ret i$(native_bit_width) %res
        """
    end
    decl = "declare i$(N) @llvm.ctpop.i$(N)(i$(N))"
    s = """
    %mask = trunc <$(N) x i8> %0 to <$(N) x i1>
    %maski = bitcast <$(N) x i1> %mask to i$(N)
    %maskipopcnt = call i$(N) @llvm.ctpop.i$(N)(i$(N) %maski)
    $(ret)
    """
    return :(
        $(Expr(:meta, :inline));
        Base.llvmcall(($decl, $s), Int, Tuple{LVec{N, Bool}}, x)
    )
end

end


text:
pyQt = pyimport("PyQt5.QtWidgets")
pyQtCore = pyimport("PyQt5.QtCore")
include("py_sdQTTable.jl") #sd_QTTable
include("py_sdQTButton.jl") #sd_QTButton
include("copyFolder.jl")

include("py_messageBox.jl")
@pydef  mutable struct sd_mainForm <: pyQt.QDialog
    function  __init__(self,parent = nothing )
        # name1 = string( name ,"Julia")
        pyQt.QDialog.__init__(self,parent)
        # self.setContextMenuPolicy(pyQtCore.ContextMenuPolicy.CustomContextMenu)
        # self.customContextMenuRequested.connect(self.contextMenuShow)
        # super(WidgetGallery, self).__init__(parent)

        self.originalPalette = pyQt.QApplication.palette()
        self.taskcp =nothing

        self.createTopLeftGroupBox()
        self.createBottomLeftTabWidget()
        topLayout = pyQt.QHBoxLayout()

        mainLayout = pyQt.QGridLayout()
        mainLayout.addLayout(topLayout, 2, 1,-1,-1)
        mainLayout.addWidget(self.topLeftGroupBox, 1, 1)
        # mainLayout.addWidget(self.topRightGroupBox, 1, 1)
        mainLayout.addWidget(self.bottomLeftTabWidget, 2, 1,-1,-1)
        mainLayout.setRowStretch(1, 1)
        mainLayout.setRowStretch(2, 10)
        mainLayout.setRowMinimumHeight(2,300)
        mainLayout.setColumnMinimumWidth(1,200)
        # mainLayout.setColumnStretch(2, 1)
        mainLayout.setColumnStretch(1, 1)
        self.setLayout(mainLayout)

        self.changeStyle("Windows")

        self.loadPairs()
    end

    function readinPairs(self)
        pairs=[]
        open("folderpairs.his","r") do file
              while !eof(file)
                 !eof(file) && (a1=readline(file))
               !eof(file) && (a2=readline(file))
            #    println(a1=>a2)
               push!(pairs ,a1 => a2)
              end
        end
        return pairs
    end
    function loadPairs(self)
        pairs =self.readinPairs()
        # println( pairs)
        self.tableWidget.setFilePairs(pairs)
    end
    function changeStyle(self, styleName)
        print( styleName)
        pyQt.QApplication.setStyle(pyQt.QStyleFactory.create(styleName))
        self.changePalette()
    end

    function changePalette(self)
        pyQt.QApplication.setPalette(pyQt.QApplication.style().standardPalette())
    end

    function showFolderLst(self)
        self.bottomLeftTabWidget.setTabPosition(1)
    end

    function showCpResult(self)

        self.bottomLeftTabWidget.setCurrentWidget( self.tab2 )
        self.textEdit.setReadOnly(true )
    end

    function appendText(self,text)
        self.textEdit.insertPlainText( text )
    end
    function clearResultText(self)
        self.textEdit.setText("")
    end
    function createTopLeftGroupBox(self)
        self.bottomLeftTabWidget = pyQt.QTabWidget()
        self.bottomLeftTabWidget.setSizePolicy(pyQt.QSizePolicy.Preferred,
                pyQt.QSizePolicy.Ignored)

        tab1 = pyQt.QWidget()
        self.tableWidget = sd_QTTable(10, 1,self)

        tab1hbox = pyQt.QHBoxLayout()
        tab1hbox.setContentsMargins(5, 5, 5, 5)
        tab1hbox.addWidget(self.tableWidget)
        tab1.setLayout(tab1hbox)

        self.tab2 = pyQt.QWidget()
        self.textEdit = pyQt.QTextEdit()

        self.textEdit.setPlainText("\n\n\n\n")

        tab2hbox = pyQt.QHBoxLayout()
        tab2hbox.setContentsMargins(5, 5, 5, 5)
        tab2hbox.addWidget(self.textEdit)
        self.tab2.setLayout(tab2hbox)

        self.bottomLeftTabWidget.addTab(tab1, "源、目标文件对照表")
        self.bottomLeftTabWidget.addTab(self.tab2, "执行结果：")
    end

    function getfolderPairs(self)
        pairs= self.tableWidget.getFilePairs()
        return pairs
        # pairs= self.getFilePairs()
    end
    function createBottomLeftTabWidget(self)
        self.topLeftGroupBox = pyQt.QGroupBox("")

        button1 = sd_QTButton("开始拷贝")
        button2 = sd_QTButton("保存对照表")
        button3 = sd_QTButton("退出")
        layout = pyQt.QHBoxLayout()
        layout.addWidget(button1)
        layout.addWidget(button2)
        layout.addWidget(button3)
        # layout.addWidget(checkBox)
        layout.addStretch(1)
        self.topLeftGroupBox.setLayout(layout)

        button3.clicked.connect( self.exitOutEventProc )
        button2.clicked.connect( self.save2logEventProc )
        button1.clicked.connect( self.cpFoldersEventProc )
    end

    function exitOutEventProc( self,idx)#退出
        self.save2logEventProc( idx)
        self.close()
    end

    function save2logEventProc( self,idx) #开始拷贝
        pairs= self.tableWidget.getFilePairs()
        length(pairs ) < 1 && return

        open("folderpairs.his","w") do file
            for idx =1:length( pairs )
                write(file,pairs[idx][1] *"\n")
                write(file,pairs[idx][2] *"\n" )
            end
        end
        # close(file)
    end
    function cpFoldersEventProc( self,idx)#保存对照表+
        if !isnothing( self.taskcp )
                println( "have task running")
                return
        end
        self.showCpResult()
        self.clearResultText()
        pairs= self.tableWidget.getFilePairs()
        # println(pairs )
        for (a1,a2 ) in pairs
          self.taskcp =  @async proCopyFolder( a1,a2)
            while !istaskdone( self.taskcp )
                a= getProgress()
                # println(a )
                self.appendText( a )
                sleep(0.01)
            end
            self.appendText( "!!!! " )
            # println("Copy Completed!!!! ")
        end
        showresult()
        self.taskcp  = nothing
    end
end


text:
###
### Experimental API
###

# This is very similar to Rust's `Try`` trait and `ControlFlow` enum.
# https://doc.rust-lang.org/std/ops/trait.Try.html
# https://doc.rust-lang.org/std/result/enum.Result.html
# https://doc.rust-lang.org/std/ops/enum.ControlFlow.html

struct Break{T}
    result::T
end

struct Continue{T}
    result::T
end

function branch end
function resultof end
function valueof end

###
### Implementation
###

branch(ok::Ok) = Continue(ok)
branch(err::Err) = Break(err)
branch(result::AbstractResult) =
    if Try.isok(result)
        Continue(result)
    else
        Break(result)
    end

resultof(br) = br.result

valueof(br::Continue{<:AbstractResult}) = Try.unwrap(br.result)
valueof(br::Break{<:AbstractResult}) = Try.unwrap_err(br.result)

branch(some::Some) = Continue(some)
branch(::Nothing) = Break(nothing)

valueof(br::Continue{<:Some}) = something(br.result)
valueof(::Break{Nothing}) = nothing

const var"@or_return" = var"@?"
macro or_return(ex)  # aka @?
    quote
        br = branch($(esc(ex)))
        if br isa Break
            return br.result
        else
            valueof(br)
        end
    end
end

macro and_return(ex)
    quote
        br = branch($(esc(ex)))
        if br isa Continue
            return br.result
        else
            valueof(br)
        end
    end
end

const var"@_return" = Try.var"@return"
macro _return(ex)
    quote
        br = branch($(esc(ex)))
        if br isa Continue
            return valueof(br)
        else
            valueof(br)
        end
    end
end

function Try.and_then(f, result)
    br = branch(result)
    if br isa Continue
        f(valueof(br))
    else
        br.result
    end
end

function Try.or_else(f, result)
    br = branch(result)
    if br isa Break
        f(valueof(br))
    else
        br.result
    end
end

function Try.unwrap_or_else(f, result)
    br = branch(result)
    if br isa Break
        f(valueof(br))
    else
        valueof(br)
    end
end

Try.and(result) = result

function Try.and(a, b)
    br = branch(a)
    if br isa Break
        br.result
    else
        b
    end
end

Try.and(a, b, c, rest...) = Try.and(Try.and(a, b), c, rest...)

Try.or(result) = result

function Try.or(a, b)
    br = branch(a)
    if br isa Continue
        br.result
    else
        b
    end
end

Try.or(a, b, c, rest...) = Try.or(Try.or(a, b), c, rest...)

macro and(ex, rest...)
    exprs = map(esc, Any[ex, rest...])
    foldr(exprs; init = pop!(exprs)) do result, ex
        br = esc(gensym(:br))
        quote
            $br = branch($result)
            if $br isa Break
                $br.result
            else
                $ex
            end
        end
    end
end

macro or(ex, rest...)
    exprs = map(esc, Any[ex, rest...])
    foldr(exprs; init = pop!(exprs)) do result, ex
        br = esc(gensym(:br))
        quote
            $br = branch($result)
            if $br isa Continue
                $br.result
            else
                $ex
            end
        end
    end
end

###
### Currying
###

# TODO: Automate currying?

function Try.and_then(f::F) where {F}
    function and_then_closure(result)
        Try.and_then(f, result)
    end
end

function Try.or_else(f::F) where {F}
    function or_else_closure(result)
        Try.or_else(f, result)
    end
end

function Try.unwrap_or_else(f::F) where {F}
    function unwrap_or_else(result)
        Try.unwrap_or_else(f, result)
    end
end


text:
using FisherMetrics
using Test

@testset "FisherMetrics.jl" begin
    # Write your tests here.
end

max_stars_repo_name: ap-browning/FisherMetrics

max_issues_repo_name: ap-browning/FisherMetrics

max_forks_repo_name: ap-browning/FisherMetrics


text:
"""
    kahanSumOneValue(a, b, c)

`a`+`b`=`sum` is done, where the lost information in that `sum` is saved in `c` (according to Kahan Summation Algorithm). Each input only represents a number.

# Args

* `a`: it is a number, where `a`>`b` should be true for the appropiate work of the algorithm.
* `b`: the value that will be added to `a`.
* `c`: the first time this function is called in an iteration, the value of `c` should be 0. After that, this value should contain the lost information of the last sum, for its proper function.

# Returns

* `sum`: `a` + `b`.
* `c`: information lost in the sum `a`+`b`.
"""
function kahanSumOneValue(a, b, c)
    sum = a
    y = b - c          # So far, so good: c is zero.
    t = sum + y        # Alas, sum is big, y small, so low-order digits of y are lost.
    c = (t - sum) - y  # (t - sum) cancels the high-order part of y; subtracting y recovers negative (low part of y)
    sum = t            # Algebraically, c should always be zero. Beware overly-aggressive optimizing compilers!
    return sum, c
end

"""
    kahanSumOneValue(a, b, c)

`a`+`b`=`sum` is done, where the lost information in that `sum` is saved in `c` (according to Kahan Summation Algorithm). Each input can represent an array or a matrix.

# Args

* `a`: it is a number, where `a`>`b` should be true for the appropiate work of the algorithm. `sum` values are overwritten here.
* `b`: the value that will be added to `a`.
* `c`: the first time this function is called in an iteration, the value of `c` should be 0. After that, this value should contain the lost information of the last sum, for its proper function. `c` values are overwritten here.
"""
function kahanSum!(a, b, c)

    if length(size(a)) == 1 # Vector
        for i in 1:length(a)
            y = b[i] - c[i]        # So far, so good: c is zero.
            t = a[i] + y           # Alas, a is big, y small, so low-order digits of y are lost.
            c[i] = (t - a[i]) - y  # (t - a) cancels the high-order part of y; subtracting y recovers negative (low part of y)
            a[i] = t               # Algebraically, c should always be zero. Beware overly-aggressive optimizing compilers!
        end
    elseif length(size(a)) == 2 # Matrix
        for i in 1:size(a)[1]
            for j in 1:size(a)[2]
                y = b[i,j] - c[i,j]        # So far, so good: c is zero.
                t = a[i,j] + y           # Alas, a is big, y small, so low-order digits of y are lost.
                c[i,j] = (t - a[i,j]) - y  # (t - a) cancels the high-order part of y; subtracting y recovers negative (low part of y)
                a[i,j] = t               # Algebraically, c should always be zero. Beware overly-aggressive optimizing compilers!
            end
        end
    else
        printf("Warning: kahanSum! works only for arrays and matrices.")
    end

end

max_stars_repo_name: salanueva/KeplerFlow

max_issues_repo_name: salanueva/KeplerFlow

max_forks_repo_name: salanueva/KeplerFlow


text:
struct InstructionList
    instructions::Vector{Tuple{Symbol,Tuple{Symbol,Any,Any}}}
    var::Symbol
    n::Base.RefValue{Int}
end

function InstructionList(; var::Symbol = :_ι, n::Base.RefValue{Int} = Ref(0))
    InstructionList(Tuple{Symbol,Tuple{Symbol,Any,Any}}[], var, n)
end

function instruction_list(
    exprs::Vector{Expression};
    subexpressions::Dict{Expression,Expression} = Dict{Expression,Expression}(),
    perform_cse::Bool = true,
)
    v = InstructionList()
    PSE = Set{Symbol}()
    if perform_cse
        exprs, CSE = cse(exprs)
        merge!(CSE, subexpressions)
    else
        CSE = subexpressions
    end
    v, map(ex -> flat_expr!(v, ex, CSE, PSE), exprs)
end

Base.length(v::InstructionList) = v.n[]
function Base.show(io::IO, v::InstructionList)
    for (id, arg) in v.instructions
        println(io, :($id = $(Expr(:call, arg...))))
    end
end

function Base.push!(v::InstructionList, x)
    id = Symbol(v.var, v.n[] += 1)
    push!(v.instructions, (id, x))
    id
end
function Base.push!(v::InstructionList, x::Pair{Symbol,<:Tuple{Symbol,Any,Any}})
    push!(v.instructions, (first(x), last(x)))
    first(x)
end

function flat_expr!(
    v,
    ex::Basic,
    # common sub expression
    CSE::Dict{Expression,Expression},
    # processed sub expression
    PSE::Set{Symbol} = Set{Symbol}(),
)
    t = class(ex)
    if t == :Symbol
        s = Symbol(to_string(ex))
        if haskey(CSE, ex) && s ∉ PSE
            push!(PSE, s)
            val = flat_expr!(v, CSE[ex], CSE, PSE)
            if val isa Symbol
                v.instructions[end] = (s, last(v.instructions[end]))
            else
                push!(v, s => val)
            end
        end
        return s
    elseif t == :Pow
        x, k = ModelKit.args(ex)
        return push!(v, (:^, flat_expr!(v, x, CSE, PSE), convert(Int, k)))
    elseif t == :Add || t == :Mul
        op = t == :Add ? :+ : :*
        vec = ModelKit.args(ex)
        op_arg = flat_expr!(v, vec[1], CSE, PSE)
        for i = 2:length(vec)
            op_arg = push!(v, (op, op_arg, flat_expr!(v, vec[i], CSE, PSE)))
        end
        return op_arg
    elseif t == :Div
        x, y = ModelKit.args(ex)
        return push!(
            v,
            (:/, flat_expr!(v, x, CSE, PSE), flat_expr!(v, y, CSE, PSE)),
        )
    else
        return to_number(ex)
    end
end

function Base.diff(
    list::InstructionList,
    vars::Vector{Symbol},
    f::Vector{Symbol},
)
    seed = Dict{Tuple{Symbol,Int},Any}()
    for (i, v) in enumerate(vars)
        seed[(v, i)] = 1
    end
    dlist = diff!(list, length(vars), seed)

    J = Matrix{Union{Nothing,Symbol,Number}}(undef, length(f), length(vars))
    for (j, v) in enumerate(vars), (i, fi) in enumerate(f)
        if haskey(seed, (fi, j))
            J[i, j] = seed[(fi, j)]
        else
            J[i, j] = nothing
        end
    end
    dlist, J
end

"""
    diff!(list::InstructionList, N::Int, diff_map)

Returns an `InstructionList` list computing the derivatives of `N` variables.
The derivatives of the instructions are stored in `diff_map`.
"""
function diff!(list::InstructionList, N::Int, diff_map)
    n = length(list)
    v = InstructionList(n = Ref(n))
    for (id, el) in list.instructions
        (op, arg1, arg2) = el

        if op == :^
            let
                p1 = p2 = :NONE
                instr_added = false
                for ∂i = 1:N
                    exp::Int = arg2
                    if haskey(diff_map, (arg1, ∂i))
                        if p2 == :NONE
                            if exp == 2
                                p2 = push!(v, (:*, 2, arg1))
                            else
                                p1 = push!(v, (:^, arg1, exp - 1))
                                p2 = push!(v, (:*, exp, p1))
                            end
                        end
                        if !instr_added
                            if exp == 2
                                push!(v, id => (:^, arg1, 2))
                            else
                                push!(v, id => (:*, p1, arg1))
                            end
                            instr_added = true
                        end
                        ∂el = diff_map[(arg1, ∂i)]
                        if ∂el != 1
                            diff_map[(id, ∂i)] = push!(v, (:*, p2, ∂el))
                        else
                            diff_map[(id, ∂i)] = p2
                        end
                    elseif p2 != :NONE && !instr_added
                        push!(v, id => el)
                        instr_added = true
                    end

                end
                if !instr_added
                    push!(v, id => el)
                end
            end
        elseif op == :*
            let
                for ∂i = 1:N
                    ∂i == 1 && push!(v, id => el)

                    has_∂1 = haskey(diff_map, (arg1, ∂i))
                    has_∂2 = haskey(diff_map, (arg2, ∂i))

                    if has_∂2
                        a2::Symbol = arg2
                        ∂arg2 = diff_map[(a2, ∂i)]
                        if ∂arg2 != 1
                            e1 = push!(v, (:*, arg1, ∂arg2))
                        else
                            e1 = arg1
                        end
                    end

                    if has_∂1
                        a1::Symbol = arg1
                        ∂arg1 = diff_map[(a1, ∂i)]
                        if ∂arg1 != 1
                            e2 = push!(v, (:*, ∂arg1, arg2))
                        else
                            e2 = arg2
                        end
                    end

                    if has_∂1 && has_∂2
                        diff_map[(id, ∂i)] = push!(v, (:+, e1, e2))
                    elseif has_∂1
                        diff_map[(id, ∂i)] = e2
                    elseif has_∂2
                        diff_map[(id, ∂i)] = e1
                    end
                end
            end
        elseif op == :/
            let
                a1::Symbol
                a2::Symbol
                for ∂i = 1:N
                    ∂i == 1 && push!(v, id => el)

                    has_∂1 = haskey(diff_map, (arg1, ∂i))
                    has_∂2 = haskey(diff_map, (arg2, ∂i))

                    if has_∂1 && has_∂2
                        a1 = arg1
                        a2 = arg2
                        ∂arg1 = diff_map[(a1, ∂i)]
                        ∂arg2 = diff_map[(a2, ∂i)]
                        e1 = push!(v, (:*, ∂arg1, a2))
                        e2 = push!(v, (:*, a1, ∂arg2))
                        e3 = push!(v, (:-, e1, e2))
                        e4 = push!(v, :(:^, arg2, 2))
                        diff_map[(id, ∂i)] = push!(v, :(:/, e3, e4))
                    elseif has_∂1
                        a1 = arg1
                        ∂arg1 = diff_map[(a1, ∂i)]
                        diff_map[(id, ∂i)] = push!(v, (:/, ∂arg1, arg2))
                    elseif has_∂2
                        a2 = arg2
                        ∂arg2 = diff_map[(a2, ∂i)]
                        e1 = push!(v, (:*, arg1, ∂arg2))
                        e2 = push!(v, (:*, -1, e1))
                        e3 = push!(v, :(:^, arg2, 2))
                        diff_map[(id, ∂i)] = push!(v, :(:/, e2, e3))
                    end
                end
            end
        elseif op == :+
            let
                for ∂i = 1:N
                    ∂i == 1 && push!(v, id => el)

                    has_∂1 = haskey(diff_map, (arg1, ∂i))
                    has_∂2 = haskey(diff_map, (arg2, ∂i))

                    if has_∂1 && has_∂2
                        diff_map[(id, ∂i)] = push!(
                            v,
                            (:+, diff_map[(arg1, ∂i)], diff_map[(arg2, ∂i)]),
                        )
                    elseif has_∂1
                        diff_map[(id, ∂i)] = diff_map[(arg1, ∂i)]
                    elseif has_∂2
                        diff_map[(id, ∂i)] = diff_map[(arg2, ∂i)]
                    end
                end
            end
        end
    end
    v
end

## Higher order AD. This follows Chapter 13 of Griewank and Walther
"""
    univariate_diff!(list::InstructionList, K::Int, diff_map)

Returns an `InstructionList` list computing the first `K` derivatives.
The derivatives of the instructions are stored in `diff_map`.
"""
function univariate_diff!(list::InstructionList, K::Int, diff_map)
    n = length(list)
    v = InstructionList(n = Ref(n))
    for (id, el) in list.instructions
        (op, arg1, arg2) = el

        if op == :^
            r::Int = arg2
            if haskey(diff_map, (arg1, 1)) && r == 2
                push!(v, id => el)
                for k = 2:K
                    w_k = nothing
                    for j = 0:k
                        u_j = j == 0 ? arg1 : get(diff_map, (arg1, j), nothing)
                        u_kj = j == k ? arg1 :
                               get(diff_map, (arg1, k - j), nothing)
                        if j < k - j && u_j !== nothing && u_kj !== nothing
                            s = push!(v, (:*, u_j, u_kj))
                            if w_k === nothing
                                w_k = s
                            else
                                w_k = push!(v, (:+, w_k, s))
                            end
                        end
                        if j ≥ k - j && w_k !== nothing
                            w_k = push!(v, (:*, 2, w_k))
                        end
                        if j == k - j
                            s = push!(v, (:*, u_j, u_kj))
                            if w_k === nothing
                                w_k = s
                            else
                                w_k = push!(v, (:+, w_k, s))
                            end
                        end
                        if j ≥ k - j
                            break
                        end
                    end
                    if w_k !== nothing
                        diff_map[(id, k)] = w_k
                    end
                end
            elseif haskey(diff_map, (arg1, 1)) && r > 2
                w_0 = push!(v, id => (:^, arg1, r))
                u_0_inv = push!(v, (:inv_not_zero, arg1, nothing))
                # w_k is v_k in Griewank
                w = Symbol[]
                for k = 1:K
                    s = nothing
                    for j = 1:k
                        u_j = get(diff_map, (arg1, j), nothing)
                        if u_j !== nothing
                            ũ_j = j == 1 ? u_j : push!(v, (:*, j, u_j))
                        else
                            ũ_j = nothing
                        end

                        w_kj = k == j ? w_0 : w[k-j]
                        if ũ_j == 1
                            s_j = w_kj
                        elseif ũ_j !== nothing
                            s_j = push!(v, (:*, w_kj, ũ_j))
                        else
                            s_j = nothing
                        end
                        if s_j !== nothing
                            s = isnothing(s) ? s_j : push!(v, (:+, s, s_j))
                        end
                    end
                    if s !== nothing
                        s = push!(v, (:*, r, s))
                    end

                    t = nothing
                    for j = 1:k-1
                        u_kj = get(diff_map, (arg1, k - j), nothing)
                        w̃_j = j == 1 ? w[1] : push!(v, (:*, j, w[j]))

                        if u_kj == 1
                            t_j = w̃_j
                        elseif u_kj !== nothing
                            t_j = push!(v, (:*, u_kj, w̃_j))
                        else
                            t_j = nothing
                        end

                        if t_j !== nothing
                            t = isnothing(t) ? t_j : push!(v, (:+, t, t_j))
                        end
                    end
                    if t == nothing && s == nothing
                        continue
                    elseif t == nothing
                        w̃_k = push!(v, (:*, u_0_inv, s))
                    else # s cannot be nothing
                        w̃_k = push!(v, (:*, u_0_inv, push!(v, (:-, s, t))))
                    end

                    if k == 1
                        w_k = w̃_k
                    else
                        w_k = push!(v, (:/, w̃_k, k))
                    end

                    push!(w, w_k)
                    diff_map[(id, k)] = w_k
                end
            else
                push!(v, id => el)
            end
        elseif op == :*
            push!(v, id => el)

            for k = 1:K
                c_k = nothing
                for j = 0:k
                    if j == 0
                        a = arg1
                        b = get(diff_map, (arg2, k - j), nothing)
                    elseif j == k
                        a = get(diff_map, (arg1, k), nothing)
                        b = arg2
                    else
                        a = get(diff_map, (arg1, j), nothing)
                        b = get(diff_map, (arg2, k - j), nothing)
                    end

                    if a === nothing || b === nothing
                        continue
                    end

                    if c_k === nothing
                        c_k = push!(v, (:*, a, b))
                    else
                        c_k = push!(v, (:+, c_k, push!(v, (:*, a, b))))
                    end
                end

                if c_k === nothing
                    break
                else
                    diff_map[(id, k)] = c_k
                end
            end
        elseif op == :+
            push!(v, id => el)
            for k = 1:K
                a = get(diff_map, (arg1, k), nothing)
                b = get(diff_map, (arg2, k), nothing)
                if a !== nothing && b !== nothing
                    diff_map[(id, k)] = push!(v, (:+, a, b))
                elseif a !== nothing
                    diff_map[(id, k)] = a
                elseif b !== nothing
                    diff_map[(id, k)] = b
                end
            end
        elseif op == :/
            error("Not implemented")
        end
    end
    v
end

@inline sqr(x) = x^2
@inline function sqr(z::Complex)
    x, y = reim(z)
    Complex((x + y) * (x - y), (x + x) * y)
end

"""
    inv_not_zero(x)

Invert x unless it is 0, then return 0.
"""
@inline inv_not_zero(x) = ifelse(iszero(x), x, inv(x))

function unroll_pow(var, n)
    n == 0 && return :(one($var))
    n == 1 && return var
    n == 2 && return :(sqr($var))
    n == 3 && return :($var * sqr($var))
    n == 4 && return :(sqr(sqr($var)))
    n == 5 && return :($var * sqr(sqr($var)))

    # base to expansion shows up which power it is needed to compute
    d = digits(n, base = 2)
    exprs = map(2:length(d)) do i
        :($(Symbol(:x, 1 << (i - 1))) = sqr($(Symbol(:x, 1 << (i - 2)))))
    end
    prods = Symbol[]
    for (i, di) in enumerate(d)
        if !iszero(di)
            push!(prods, Symbol(:x, 1 << (i - 1)))
        end
    end
    if length(prods) > 1
        push!(exprs, Expr(:call, :*, prods...))
    end
    Expr(:let, :(x1 = $var), Expr(:block, exprs...))
end

function to_expr(
    list::InstructionList,
    var_map = Dict{Symbol,Union{Expr,Symbol}}(),
    assignements = Dict{Symbol,Expr}(),
)
    exprs = Expr[]
    for (id, (op, arg1, arg2)) in list.instructions
        if op == :^
            x::Symbol = arg1
            k::Int = arg2
            if k < 0
                push!(
                    exprs,
                    :($id = inv($(unroll_pow(get(var_map, x, x), -k)))),
                )
            else
                push!(exprs, :($id = $(unroll_pow(get(var_map, x, x), k))))
            end
        else
            a = get(var_map, arg1, arg1)
            if arg2 !== nothing
                b = get(var_map, arg2, arg2)
                push!(exprs, :($id = $(Expr(:call, op, a, b))))
            else
                push!(exprs, :($id = $(Expr(:call, op, a))))
            end
        end

        if haskey(assignements, id)
            push!(exprs, assignements[id])
        end
    end
    Expr(:block, exprs...)
end


text:
function dump_bandstructure( evals, kpath; filename="TEMP_bands.dat")
    Nkpt = size(kpath)[2]
    Nstates = size(evals)[1]

    Nkspin = size(evals)[2]
    @assert( Nkpt == Nkspin )

    Xcoords = Array{Float64}(undef,Nkpt)
    Xcoords[1] = 0.0
    for ik = 2:Nkpt
        dk = kpath[:,ik] - kpath[:,ik-1]
        Xcoords[ik] = Xcoords[ik-1] + norm( dk )
    end

    f = open(filename,"w")

    for ik = 1:Nkpt
        @printf(f, "%18.10f ", Xcoords[ik])
        for ist = 1:Nstates
            @printf(f, "%18.10f ", evals[ist,ik])
        end
        @printf(f, "\n")
    end
    close(f)
end


text:
using Test
using VariantsNMF

X1 = [0.21 0.40 0.4; 0.01 0.4 0.45]
Y1 = [0.12 0.48 0.4; 0.1 0.55 0.45]

α_values = [-0.5, 0, 0.33, 1, 1.5, 2, 2.8]
α_divg_values = [0.283580, 0.195769, 0.162286, 0.124183, 0.109648, 0.101371, 0.095786]
β_values = [-0.5, 0, 0.33, 1, 1.5, 2, 2.8]
β_divg_values = [7.81200, 1.65435, 0.64633, 0.12418, 0.04732, 0.02255, 0.00952]

# test α-divergence function
@testset "α_divergence" begin
    for (α, α_divg) in zip(α_values, α_divg_values)
        @test α_divergence(X1, X1, α) ≈ 0  atol = 1e-9
        @test α_divergence(Y1, Y1, α) ≈ 0  atol = 1e-9
        @test α_divergence(X1, Y1, α) ≈ α_divg atol = 1e-5
    end
end

# test β-divegence function
@testset "β_divergence" begin
    for (β, β_divg) in zip(β_values, β_divg_values)
        @test β_divergence(X1, X1, β) ≈ 0  atol = 1e-9
        @test β_divergence(Y1, Y1, β) ≈ 0  atol = 1e-9
        @test β_divergence(X1, Y1, β) ≈ β_divg atol = 1e-5
    end
end


text:
function initWB!(
    cLayer::Conv1D,
    p::Type{T} = Float64::Type{Float64};
    He::Bool = true,
    coef::AbstractFloat = 0.01,
    zro::Bool = false,
) where {T}

    f = cLayer.f
    cl = cLayer.channels
    cl_1 = cLayer.prevLayer.channels

    inputS = cLayer.inputS
    s_H = cLayer.s
    f_H = cLayer.f

    # ## note this is the input to the matmult after padding
    # if cLayer.padding == :same
    #     p_H = Integer(ceil((s_H * (n_Hi - 1) - n_Hi + f_H) / 2))
    #     n_H = n_Hi + 2p_H
    # elseif cLayer.padding == :valid
    #     n_H = n_Hi
    # end

    if He
        coef = sqrt(2 / cl_1)
    end

    if zro
        W = zeros(T, f..., cl_1, cl)
    else
        W = T.(randn(f..., cl_1, cl) .* coef)
    end
    B = zeros(T, repeat([1], length(f) + 1)..., cl)

    cLayer.W, cLayer.B = W, B
    cLayer.K = Matrix{T}(undef,0,0)
    #unroll(cLayer, (n_H, ci, m))
    cLayer.dW = zeros(T, f..., cl_1, cl)
    cLayer.dK = Matrix{T}(undef,0,0)
    #unroll(cLayer, (n_H, ci, m), :dW)
    cLayer.dB = deepcopy(B)
    return nothing
end #initWB

function initWB!(
    cLayer::Conv2D,
    p::Type{T} = Float64::Type{Float64};
    He::Bool = true,
    coef::AbstractFloat = 0.01,
    zro::Bool = false,
) where {T}

    f = cLayer.f
    cl = cLayer.channels
    cl_1 = cLayer.prevLayer.channels

    inputS = cLayer.inputS
    s_H, s_W = cLayer.s
    f_H, f_W = cLayer.f

    ## note this is the input to the matmult after padding
    # if cLayer.padding == :same
    #     p_H = Integer(ceil((s_H * (n_Hi - 1) - n_Hi + f_H) / 2))
    #     p_W = Integer(ceil((s_W * (n_Wi - 1) - n_Wi + f_W) / 2))
    #     n_H, n_W = n_Hi + 2p_H, n_Wi + 2p_W
    # elseif cLayer.padding == :valid
    #     n_H = n_Hi
    #     n_W = n_Wi
    # end

    if He
        coef = sqrt(2 / cl_1)
    end

    if zro
        W = zeros(T, f..., cl_1, cl)
    else
        W = T.(randn(f..., cl_1, cl) .* coef)
    end
    B = zeros(T, repeat([1], length(f) + 1)..., cl)

    cLayer.W, cLayer.B = W, B
    cLayer.K = Matrix{T}(undef,0,0)
    #unroll(cLayer, (n_H, n_W, ci, m))
    cLayer.dW = zeros(T, f..., cl_1, cl)
    cLayer.dK = Matrix{T}(undef,0,0)
    #unroll(cLayer, (n_H, n_W, ci, m), :dW)
    cLayer.dB = deepcopy(B)
    return nothing
end #initWB

function initWB!(
    cLayer::Conv3D,
    p::Type{T} = Float64::Type{Float64};
    He::Bool= true,
    coef::AbstractFloat = 0.01,
    zro::Bool = false,
) where {T}

    f = cLayer.f
    cl = cLayer.channels
    cl_1 = cLayer.prevLayer.channels

    inputS = cLayer.inputS
    s_H, s_W, s_D = cLayer.s
    f_H, f_W, f_D = cLayer.f

    ## note this is the input to the matmult after padding
    # if cLayer.padding == :same
    #     p_H = Integer(ceil((s_H * (n_Hi - 1) - n_Hi + f_H) / 2))
    #     p_W = Integer(ceil((s_W * (n_Wi - 1) - n_Wi + f_W) / 2))
    #     p_D = Integer(ceil((s_D * (n_Di - 1) - n_Di + f_D) / 2))
    #     n_H, n_W, n_D = n_Hi + 2p_H, n_Wi + 2p_W, n_Di + 2p_D
    # elseif cLayer.padding == :valid
    #     n_H = n_Hi
    #     n_W = n_Wi
    #     n_D = n_Di
    # end

    if He
        coef = sqrt(2 / cl_1)
    end

    if zro
        W = zeros(T, f..., cl_1, cl)
    else
        W = T.(randn(f..., cl_1, cl) .* coef)
    end
    B = zeros(T, repeat([1], length(f) + 1)..., cl)

    cLayer.W, cLayer.B = W, B
    cLayer.K = Matrix{T}(undef,0,0)
    #unroll(cLayer, (n_H, n_W, n_D, ci, m))
    cLayer.dW = zeros(T, f..., cl_1, cl)
    cLayer.dK = Matrix{T}(undef,0,0)
    #unroll(cLayer, (n_H, n_W, n_D, ci, m), :dW)
    cLayer.dB = deepcopy(B)
    return nothing
end #initWB

###Pooling Layers

function initWB!(
    cLayer::P,
    p::Type{T} = Float64::Type{Float64};
    He::Bool = true,
    coef::AbstractFloat = 0.01,
    zro::Bool = false,
) where {T,P<:PoolLayer}

    return nothing
end

export initWB!

### initVS!

function initVS!(cLayer::CL, optimizer::Symbol) where {CL <: ConvLayer}

    if optimizer == :adam || optimizer == :momentum
        cLayer.V[:dw] = deepcopy(cLayer.dW)
        cLayer.V̂dk = deepcopy(cLayer.dK)
        cLayer.V[:db] = deepcopy(cLayer.dB)
        if optimizer == :adam
            cLayer.S = deepcopy(cLayer.V)
            cLayer.Ŝdk = deepcopy(cLayer.V̂dk)
        end #if optimizer == :adam
    end #if optimizer == :adam || optimizer == :momentum

    return nothing
end #function initVS!

initVS!(cLayer::P, optimizer::Symbol) where {P<:PoolLayer} = nothing

export initVS!


text:
using Pkg
ENV["PYTHON"] = "/root/miniconda3/envs/condy/bin/python3"
#Pkg.add("anotherPackageHere")
Pkg.add("PyCall")
Pkg.build("PyCall")
Pkg.instantiate()
Pkg.precompile()


text:
export HS

# References:
# * Lefkimmiatis, Stamatios, John Paul Ward, and Michael Unser. "Hessian Schatten-norm regularization for linear inverse problems." IEEE transactions on image processing 22.5 (2013): 1873-1888.
# * Lefkimmiatis, Stamatios, and Michael Unser. "Poisson image reconstruction with Hessian Schatten-norm regularization." IEEE transactions on image processing 22.11 (2013): 4314-4327.

function Δr1r1(x)
    return @tullio res[i, j] := x[i+2, j+0] - 2 * x[i+1, j] + x[i, j] (i in 1:size(x)[1]-2, j in 1:size(x)[2]-2)
end

function Δr2r2(x)
    return @tullio res[i, j] := x[i+0, j+2] - 2 * x[i, j+1] + x[i, j] (i in 1:size(x)[1]-2, j in 1:size(x)[2]-2)
end

function Δr1r2(x)
    return @tullio res[i, j] := x[i+1, j+1] - x[i+1, j] - x[i, j+1]  + x[i, j] (i in 1:size(x)[1]-2, j in 1:size(x)[2]-2)
end

"""
    HS(; p=1)

Hessian Schatten norm. `p` determines which Schatten norm is used.

This regularizer only works with 2D arrays at the moment.
"""
function HS(;p=1)
    if isone(p)
        return HS1
    end

    f(x) = HSp(x, p=p)
    return f
end

"""
Hessian schatten norm for p=1 efficiently with Tullio.
"""
function HS1(arr)
    H11 = Δr1r1(arr)
    H22 = Δr2r2(arr)
    return schatten_norm_1(H11, H22)
end

function schatten_norm_1(a, d)
    @tullio A[i, j] := a[i, j] + d[i, j]
    @tullio res = abs(1f-8 + A[i, j])
end

"""
Hessian schatten norm for p.
But not as fast as p=1
"""
function HSp(arr; p=1)
    H11 = Δr1r1(arr)
    H22 = Δr2r2(arr)
    H12 = Δr1r2(arr)

    res = schatten_norm_tullio(H11, H12, H22, p)
    return sum(res)
end

function schatten_norm(H11, H12, H22, p)
    λ₁, λ₂ = eigvals_symmetric(H11, H12, H22)
    return (λ₁^p + λ₂^p )^(1/p)
end

function schatten_norm_tullio(H11, H12, H22, p)
    λ₁, λ₂ = eigvals_symmetric_tullio(H11, H12, H22)
    return @tullio res = abs((1f-8 + λ₁[i, j]^p + λ₂[i, j]^p))^(1/p)
end

"""
    eigvals_symmetric(a,b,c)

Calculate the eigenvalues of the matrix
[a b; b d] analytically.
"""
function eigvals_symmetric(a, b, d)
    A = a+d
    B = sqrt((a-d)^2+4*b^2)
    λ₁ = 0.5 * (A + B)
    λ₂ = 0.5 * (A - B)
    return λ₁, λ₂
end

function eigvals_symmetric_tullio(a, b, d)
    @tullio A[i, j] := a[i, j] + d[i, j]
    @tullio B[i, j] := sqrt(1f-8 + (a[i, j]-d[i, j])^2+4*b[i, j]^2)
    @tullio λ₁[i, j] := 0.5 * (A[i, j] + B[i, j])
    @tullio λ₂[i, j] := 0.5 * (A[i, j] - B[i, j])
    return λ₁, λ₂
end


text:
# Error type and decoding functionality

export ZeError

struct ZeError <: Exception
    code::ze_result_t
end

Base.convert(::Type{ze_result_t}, err::ZeError) = err.code

Base.showerror(io::IO, err::ZeError) =
    print(io, "ZeError: ", description(err), " (code $(reinterpret(Int32, err.code)), $(name(err)))")

Base.show(io::IO, ::MIME"text/plain", err::ZeError) = print(io, "ZeError($(err.code))")

name(err::ZeError) = string(err.code)

## COV_EXCL_START
function description(err::ZeError)
    if err.code == ZE_RESULT_SUCCESS
        "success"
    elseif err.code == ZE_RESULT_NOT_READY
        "synchronization primitive not signaled"
    elseif err.code == ZE_RESULT_ERROR_DEVICE_LOST
        "device hung, reset, was removed, or driver update occurred"
    elseif err.code == ZE_RESULT_ERROR_OUT_OF_HOST_MEMORY
        "insufficient host memory to satisfy call"
    elseif err.code == ZE_RESULT_ERROR_OUT_OF_DEVICE_MEMORY
        "insufficient device memory to satisfy call"
    elseif err.code == ZE_RESULT_ERROR_MODULE_BUILD_FAILURE
        "error occurred when building module, see build log for details"
    elseif err.code == ZE_RESULT_ERROR_INSUFFICIENT_PERMISSIONS
        "access denied due to permission level"
    elseif err.code == ZE_RESULT_ERROR_NOT_AVAILABLE
        "resource already in use and simultaneous access not allowed"
    elseif err.code == ZE_RESULT_ERROR_UNINITIALIZED
        "driver is not initialized"
    elseif err.code == ZE_RESULT_ERROR_UNSUPPORTED_VERSION
        "generic error code for unsupported versions"
    elseif err.code == ZE_RESULT_ERROR_UNSUPPORTED_FEATURE
        "generic error code for unsupported features"
    elseif err.code == ZE_RESULT_ERROR_INVALID_ARGUMENT
        "generic error code for invalid arguments"
    elseif err.code == ZE_RESULT_ERROR_INVALID_NULL_HANDLE
        "handle argument is not valid"
    elseif err.code == ZE_RESULT_ERROR_HANDLE_OBJECT_IN_USE
        "object pointed to by handle still in-use by device"
    elseif err.code == ZE_RESULT_ERROR_INVALID_NULL_POINTER
        "pointer argument may not be nullptr"
    elseif err.code == ZE_RESULT_ERROR_INVALID_SIZE
        "size argument is invalid (e.g., must not be zero)"
    elseif err.code == ZE_RESULT_ERROR_UNSUPPORTED_SIZE
        "size argument is not supported by the device (e.g., too large)"
    elseif err.code == ZE_RESULT_ERROR_UNSUPPORTED_ALIGNMENT
        "alignment argument is not supported by the device (e.g., too small)"
    elseif err.code == ZE_RESULT_ERROR_INVALID_SYNCHRONIZATION_OBJECT
        "synchronization object in invalid state"
    elseif err.code == ZE_RESULT_ERROR_INVALID_ENUMERATION
        "enumerator argument is not valid"
    elseif err.code == ZE_RESULT_ERROR_UNSUPPORTED_ENUMERATION
        "enumerator argument is not supported by the device"
    elseif err.code == ZE_RESULT_ERROR_UNSUPPORTED_IMAGE_FORMAT
        "image format is not supported by the device"
    elseif err.code == ZE_RESULT_ERROR_INVALID_NATIVE_BINARY
        "native binary is not supported by the device"
    elseif err.code == ZE_RESULT_ERROR_INVALID_GLOBAL_NAME
        "global variable is not found in the module"
    elseif err.code == ZE_RESULT_ERROR_INVALID_KERNEL_NAME
        "kernel name is not found in the module"
    elseif err.code == ZE_RESULT_ERROR_INVALID_FUNCTION_NAME
        "function name is not found in the module"
    elseif err.code == ZE_RESULT_ERROR_INVALID_GROUP_SIZE_DIMENSION
        "group size dimension is not valid for the kernel or device"
    elseif err.code == ZE_RESULT_ERROR_INVALID_GLOBAL_WIDTH_DIMENSION
        "global width dimension is not valid for the kernel or device"
    elseif err.code == ZE_RESULT_ERROR_INVALID_KERNEL_ARGUMENT_INDEX
        "kernel argument index is not valid for kernel"
    elseif err.code == ZE_RESULT_ERROR_INVALID_KERNEL_ARGUMENT_SIZE
        "kernel argument size does not match kernel"
    elseif err.code == ZE_RESULT_ERROR_INVALID_KERNEL_ATTRIBUTE_VALUE
        "value of kernel attribute is not valid for the kernel or device"
    elseif err.code == ZE_RESULT_ERROR_INVALID_COMMAND_LIST_TYPE
        "command list type does not match command queue type"
    elseif err.code == ZE_RESULT_ERROR_OVERLAPPING_REGIONS
        "copy operations do not support overlapping regions of memory"
    elseif err.code == ZE_RESULT_ERROR_UNKNOWN
        "unknown or internal error"
    else
        "no description for this error"
    end
end
## COV_EXCL_STOP

@enum_without_prefix _ze_result_t ZE_

## API call wrapper

# outlined functionality to avoid GC frame allocation
@noinline function throw_api_error(res)
    throw(ZeError(res))
end


text:
function same_size(arrs...)
    s1 = size(arrs[1])
    return all(x -> size(x) == s1, arrs)
end


text:
using Test
import Varthaka

@test cmp(Varthaka.test_me(), "Varthaka") == 0


text:
using MySQL
using FTPClient
using BioMedQuery.PubMed
using BioMedQuery.DBUtils
using LightXML
using DataFrames
using Distributed

"""
    load_medline(db_con, output_dir; start_file=1, end_file=972, year=2019, test=false)

Given a MySQL connection and optionally the start and end files, fetches the medline files, parses the xml, and loads into a MySQL DB (assumes tables already exist). The raw (xml.gz) and parsed (csv) files will be stored in the output_dir.

## Arguments
* `db_con` : A MySQL Connection to a db (tables must already be created - see `PubMed.create_tables!`)
* `output_dir` : root directory where the raw and parsed files should be stored
* `start_file` : which medline file should the loading start at
* `end_file` : which medline file should the loading end at (default is last file in 2018 baseline)
* `year` : which year medline is (current is 2018)
* `test` : if true, a sample file will be downloaded, parsed, and loaded instead of the baseline files
"""
function load_medline!(db_con::MySQL.Connection, output_dir::String; start_file::Int=1, end_file::Int=972, year::Int=2019, test::Bool=false)

    ftp_con = init_medline(output_dir, test)

    set_innodb_checks!(db_con,0,0,0)

    if test
        start_file = 1
        end_file = 1
    end

    @info "Getting files from Medline"
    @sync for n = start_file:end_file
        @async get_ml_file(get_file_name(n, year, test), ftp_con, output_dir, test)
    end

    @info "Parsing files into CSV"
    pmap(x -> parse_ml_file(get_file_name(x, year, test), output_dir), start_file:end_file)

    @info "Loading CSVs into MySQL"
    for n = start_file:end_file

        fname = get_file_name(n, year, test) ::String
        println("Loading file: ", fname)

        csv_prefix = "$(fname[1:end-7])_"
        csv_path = joinpath(output_dir,"medline","parsed_files")

        db_insert!(db_con, csv_path, csv_prefix)
    end

    set_innodb_checks!(db_con)

    @info "All files processed - closing FTP connection"
    close_cons(ftp_con)

    return nothing
end

"""
    init_medline(output_dir, test=false)

Sets up environment (folders), and connects to medline FTP Server and returns the connection.
"""
function init_medline(output_dir::String, test::Bool=false)
    ## SET UP ENVIRONMENT
    @info "======Setting up folders and creating FTP Connection======"

    ml_dir = joinpath(output_dir,"medline")

    if isdir(ml_dir)
        println("medline directory already exists")
    else
        mkdir(ml_dir)
    end

    if isdir(joinpath(ml_dir, "raw_files"))
        println("raw files directory already exists")
    else
        mkdir(joinpath(ml_dir, "raw_files"))
    end

    if isdir(joinpath(ml_dir, "parsed_files"))
        println("parsed files directory already exists")
    else
        mkdir(joinpath(ml_dir, "parsed_files"))
    end

    # Initialize FTP
    ftp_init()

    ftp_con = get_ftp_con(test)

    return ftp_con
end

"""
    get_file_name(fnum::Int, year::Int = 2018, test = false)
Returns the medline file name given the file number and year.
"""
function get_file_name(fnum::Int, year::Int, test::Bool=false)
    nstr = lpad(string(fnum), 4, "0") # pad iterator with leading zeros so total length is 4
    y2 = string(year)[3:4]
    if test
        y2 = "sample" * y2
    end
    return "pubmed$(y2)n$nstr.xml.gz"
end

"""
    get_ml_file(fname::String, conn::ConnContext, output_dir)

Retrieves the file with fname and puts in medline/raw_files.  Returns the HTTP response.
"""
function get_ml_file(fname::String, conn::ConnContext, output_dir::String, test::Bool = false)
    println("Getting file: ", fname)
    # get file

    if test
        # copy pubmedsample to right folder
        source = joinpath(@__DIR__, "..", "..","test",fname)
        target = joinpath(output_dir,"medline","raw_files",fname)
        cp(source, target, force=true)
    end

    path = joinpath(output_dir,"medline","raw_files",fname)
    if isfile(path)
        resp = "File already exists, using local file"
    else
        # conn = get_ftp_con()
        resp = ftp_get(conn, fname, path)
        # ftp_close_connection(conn)
    end

    return resp
end

"""
    get_ftp_con(test = false)
Get an FTP connection
"""
function get_ftp_con(test::Bool = false)
    options = RequestOptions("ftp://ftp.ncbi.nlm.nih.gov/pubmed/baseline/")

    try
        global conn = ftp_connect(options) # returns connection and response
    catch e
        if test
            @warn "FTP error during package test"
            return ConnContext(options)
        else
            rethrow(e)
        end
    end

    return conn[1]# get ConnContext object
end

"""
    parse_ml_file(fname::String, output_dir::String)

Parses the medline xml file into a dictionary of dataframes. Saves the resulting CSV files to medline/parsed_files.
"""
function parse_ml_file(fname::String, output_dir::String)
    println("Parsing file: ", fname)

    parsed_path = joinpath(output_dir,"medline","parsed_files","$(fname[1:end-7])_basic.csv")
    if isfile(parsed_path)
        resp = "File already exists, using local file"
    else
        path = joinpath(output_dir,"medline","raw_files",fname)
        doc = parse_file(path)
        raw_articles = root(doc)

        dfs = PubMed.parse_articles(raw_articles)

        dfs_to_csv(dfs, joinpath(output_dir,"medline","parsed_files"), "$(fname[1:end-7])_")

        free(doc)
    end

    return nothing
end

"""
    close_cons(ftp_con)
closes connection and cleans up
"""
function close_cons(ftp_con::ConnContext)
    # Close FTP Connection
    ftp_close_connection(ftp_con)
    ftp_cleanup()

    return nothing
end

max_stars_repo_name: bcbi/NLM


text:
abstract type SubstitutionModel end

const SM = SubstitutionModel

"""
`NucleicAcidSubstitutionModel` is an abstract type that contains all models
describing a substitution process impacting biological sequences of `DNA` or
`RNA` with continous time Markov models.
"""
abstract type NucleicAcidSubstitutionModel <: SubstitutionModel end

const NASM = NucleicAcidSubstitutionModel

const Qmatrix = SMatrix{4, 4, Float64}

const Pmatrix = SMatrix{4, 4, Float64}


text:
using Agents
using Random
using DataFrames
using Distributions
using DrWatson
using StatsPlots
using BenchmarkTools

function rate_to_proportion(r::Float64,t::Float64)
    1-exp(-r*t)
end;

mutable struct Person <: AbstractAgent
    id::Int64
    status::Symbol
end

function init_model(β::Float64,c::Float64,γ::Float64,N::Int64,I0::Int64)
    properties = @dict(β,c,γ)
    model = ABM(Person; properties=properties)
    for i in 1:N
        if i <= I0
            s = :I
        else
            s = :S
        end
        p = Person(i,s)
        p = add_agent!(p,model)
    end
    return model
end;

function agent_step!(agent, model)
    transmit!(agent, model)
    recover!(agent, model)
end;

function transmit!(agent, model)
    # If I'm not susceptible, I return
    agent.status != :S && return
    ncontacts = rand(Poisson(model.properties[:c]))
    for i in 1:ncontacts
        # Choose random individual
        alter = random_agent(model)
        if alter.status == :I && (rand() ≤ model.properties[:β])
            # An infection occurs
            agent.status = :I
            break
        end
    end
end;

function recover!(agent, model)
    agent.status != :I && return
    if rand() ≤ model.properties[:γ]
            agent.status = :R
    end
end;

susceptible(x) = count(i == :S for i in x)
infected(x) = count(i == :I for i in x)
recovered(x) = count(i == :R for i in x);

δt = 0.1
nsteps = 400
tf = nsteps*δt
t = 0:δt:tf;

β = 0.05
c = 10.0*δt
γ = rate_to_proportion(0.25,δt);

N = 1000
I0 = 10;

Random.seed!(1234);

abm_model = init_model(β,c,γ,N,I0)

to_collect = [(:status, f) for f in (susceptible, infected, recovered)]
abm_data, _ = run!(abm_model, agent_step!, nsteps; adata = to_collect);

abm_data[!,:t] = t;

plot(t,abm_data[:,2],label="S",xlab="Time",ylabel="Number")
plot!(t,abm_data[:,3],label="I")
plot!(t,abm_data[:,4],label="R")

max_stars_repo_name: medewitt/covidresearch

max_issues_repo_name: medewitt/covidresearch

max_forks_repo_name: medewitt/covidresearch


text:
using AlphaGo
using AlphaGo: stone_features
using AlphaGo: GoPosition, PlayerMove, BLACK, WHITE, play_move!
using Base.Test

include("test_utils.jl")
env = GoEnv(9)
@testset "features" begin
  EMPTY_ROW = repeat(".", env.N) * '\n'
  TEST_BOARD = load_board("""
                          .X.....OO
                          X........
                          """ * repeat(EMPTY_ROW, 7), env)
  TEST_POSITION = Position(env,
      board = TEST_BOARD,
      n = 3,
      komi = 6.5,
      caps = (1, 2),
      ko = nothing,
      recent = [PlayerMove(BLACK, (1, 2)),
              PlayerMove(WHITE, (1, 9)),
              PlayerMove(BLACK, (2, 1))],
      to_play = BLACK);

  TEST_BOARD2 = load_board("""
                          .XOXXOO..
                          XO.OXOX..
                          XXO..X...
                          """ * repeat(EMPTY_ROW, 6), env)

  TEST_POSITION2 = Position(env,
      board = TEST_BOARD2,
      n = 0,
      komi = 6.5,
      caps = (0, 0),
      ko = nothing,
      recent = Vector{PlayerMove}(),
      to_play = BLACK);

  TEST_POSITION3 = Position(env);
  for coord in ((1, 1), (1, 2), (1, 3), (1, 4), (2, 2))
    play_move!(TEST_POSITION3, coord, mutate = true)
  end
  # resulting position should look like this:
  # X.XO.....
  # .X.......
  # .........

  @testset "stone_features" begin
    f = stone_features(TEST_POSITION3)
    @test TEST_POSITION3.to_play == WHITE
    @test size(f) == (9, 9, 16)
    @test assertEqualArray(f[:, :, 1], load_board("""
                                              ...X.....
                                              .........""" * repeat(EMPTY_ROW, 7), env))

    @test assertEqualArray(f[:, :, 2], load_board("""
                                              X.X......
                                              .X.......""" * repeat(EMPTY_ROW, 7), env))

    @test assertEqualArray(f[:, :, 3], load_board("""
                                              .X.X.....
                                              .........""" * repeat(EMPTY_ROW, 7), env))

    @test assertEqualArray(f[:, :, 4], load_board("""
                                              X.X......
                                              .........""" * repeat(EMPTY_ROW, 7), env))

    @test assertEqualArray(f[:, :, 5], load_board("""
                                              .X.......
                                              .........""" * repeat(EMPTY_ROW, 7), env))

    @test assertEqualArray(f[:, :, 6], load_board("""
                                              X.X......
                                              .........""" * repeat(EMPTY_ROW, 7), env))

    for i = 11:16
      @test assertEqualArray(f[:, :, i], zeros(env.N, env.N))
    end
  end
end


text:
using JuMP
using Base, GLPK, Clp, DataStructures, MathOptInterface
import MathOptInterface: VariablePrimalStart
import DataStructures: PriorityQueue, enqueue!, dequeue!
import Base.Order.Reverse

export solve_sp, find_w

## utilities

function find_w(f::Function,df::Function,l::Number,u::Number,z::Number)
    if l >= z
        w = l
    elseif u <= z
        w = u
    else
        g = w -> df(w)*(w-l) - (f(w)-f(l))
        if g(u) >= 0
            w = u
        else
            w = bisection( g, z, u )
        end
    end
    w
end

function bisection(f, a, b, tol=1e-9, maxiters=1000)
    @assert(round(f(a), digits=9)>=0)
    @assert(f(b)<0)
    for i=1:maxiters
        mid = a + (b-a)/2
        fmid = f(mid)
        if abs(fmid) < tol
            return mid
        end
        if f(mid) > 0
            a = mid
        else
            b = mid
        end
    end
    @warn "hit maximum iterations in bisection search"
    return (b+a)/2
end

## calculate convex hull
function calculate_hull(x_val, w, l, fs)
    nvar = size(w)[1]
    t = zeros(nvar)
    for i=1:nvar
        xi = x_val[i]
        if xi >= w[i]
            t[i] = fs[i](xi)
        else
            slopeatl = (fs[i](w[i]) - fs[i](l[i]))/(w[i] - l[i])
            offsetatl = fs[i](l[i])
            t[i] = offsetatl + slopeatl*(xi - l[i])
        end
    end
    return t
end

function model_problem(
    l,
    u,
    w,
    problem::SigmoidalProgram,
    m = Model(optimizer_with_attributes(() -> GLPK.Optimizer(want_infeasibility_certificates=false), "tm_lim" => 10000, "msg_lev" => GLPK.GLP_MSG_OFF))
)
   # Clp solver also works with the following syntax
   # m = Model(optimizer_with_attributes(Clp.Optimizer, "LogLevel" => 0))
   # m = Model(optimizer_with_attributes(GLPK.Optimizer,"tm_lim" => 60000, "msg_lev" => GLPK.GLP_MSG_OFF))

   nvar = length(l)
   fs,dfs = problem.fs, problem.dfs

   # Define our variables to be inside a box
   @variable(m, x[i=1:nvar])
   for i=1:nvar
       set_lower_bound(x[i], l[i])
       set_upper_bound(x[i], u[i])
   end
   # epigraph variable
   @variable(m, t[i=1:nvar])

   # Require that t be in the hypograph of fhat, approximating as pwl function
   # At first, we add only the bit of fhat from l to w, and the tangent at u
   for i=1:nvar
       if w[i] > l[i]
           slopeatl = (fs[i](w[i]) - fs[i](l[i]))/(w[i] - l[i])
           offsetatl = fs[i](l[i])
           @constraint(m, t[i] <= offsetatl + slopeatl*(x[i] - l[i]))
       else
           @constraint(m, t[i] <= fs[i](l[i]) + dfs[i](l[i])*(x[i] - l[i]))
       end
       @constraint(m, t[i] <= fs[i](u[i]) + dfs[i](u[i])*(x[i] - u[i]))
   end
   # Add other problem constraints
   addConstraints!(m, x, problem)

   @objective(m, Max, sum(t))
   return m
end

## maximize concave hull
function maximize_fhat(l, u, w, problem::SigmoidalProgram,
                       m = model_problem(l, u, w, problem);
                       maxiters = 10, TOL = 1e-6, verbose=0)
                       # init_x=zeros(length(l))) # GLPK doesn't allow initialization, but could be useful for other solvers

    nvar = length(l)
    maxiters *= nvar
    fs,dfs = problem.fs, problem.dfs
    x = m[:x]
    t = m[:t]

    # Now solve and add hypograph constraints until the solution stabilizes
    start = time()
    try
        optimize!(m)
    catch y
    end
    elapsed_time = time() - start
    status = termination_status(m)

    if elapsed_time >= 10
        status = "time_limit"
        return zeros(1, nvar), zeros(1, nvar), status
    end

    for i=1:maxiters
        #println("maximize_fhat iteration: ", i)
        if status == MathOptInterface.OPTIMAL
            x_val = value.(x)
            t_val = value.(t)

            solved = true

            for i=1:length(x_val)
                # Check if t is in the hypograph of f, allowing some tolerance
                xi, ti = x_val[i], t_val[i]
                if xi > w[i]
                    if ti > fs[i](xi) + TOL
                        solved = false
                        @constraint(m, t[i] <= fs[i](xi) + dfs[i](xi)*(x[i] - xi))
                    end
                end
            end
            if solved
                if verbose>=2 println("solved problem to within $TOL in $i iterations") end
                break
            else
                start = time()
                try
                    optimize!(m)
                catch y
                end
                elapsed_time = time() - start
                status = termination_status(m)
                if elapsed_time >= 10
                    status = "time_limit"
                    break
                end
            end
        else
            break
        end
    end
    # refine t a bit to make sure it's really on the convex hull
    if status == MathOptInterface.OPTIMAL
        x_val = value.(x)
        t = calculate_hull(x_val, w, l, fs)
        return x_val, t, status
    else
        return zeros(1, nvar), zeros(1, nvar), status
        #return -Inf, -Inf, status
    end
end

## Nodes of the branch and bound tree
struct Node
    l::Array{Float64,1}
    u::Array{Float64,1}
    w::Array{Float64,1}
    x #::Array{Float64,1} no, it's a jump dict
    lb::Float64
    ub::Float64
    maxdiff_index::Int64
    m # JUMP model
    function Node(l,u,w,problem,init_x=Nothing,
                  m = model_problem(l, u, w, problem);
                  kwargs...)
        nvar = length(l)
        # find upper and lower bounds
        if init_x != Nothing
            x = init_x
            s = Float64[problem.fs[i](x[i]) for i=1:size(x)[1]]
            lb = sum(s)
            t = calculate_hull(x, w, l, problem.fs)
            ub = sum(t)
            maxdiff_index = argmax(t-s)
        else
            x, t, status = maximize_fhat(l, u, w, problem, m; kwargs...)
            if status==MathOptInterface.OPTIMAL
                x = max.(x, l)
                s = Float64[problem.fs[i](x[i]) for i=1:nvar]
                ub = sum(t)
                lb = sum(s)
                maxdiff_index = argmax(t-s)
            elseif status=="time_limit"
                ub = -Inf; lb = -Inf; maxdiff_index = -1
            else
                ub = -Inf; lb = -Inf; maxdiff_index = 1
            end
        end
        new(l,u,w,x,lb,ub,maxdiff_index,m)
    end
end

function Node(l,u,problem::SigmoidalProgram, init_x; kwargs...)
    nvar = length(l)
    # find w
    w = zeros(nvar)
    for i=1:nvar
        w[i] = find_w(problem.fs[i],problem.dfs[i],l[i],u[i],problem.z[i])
    end
    Node(l,u,w,problem,init_x; kwargs...)
end

## Branching rule
function split(n::Node, problem::SigmoidalProgram, verbose=0; kwargs...)
    i = n.maxdiff_index
    # split at x for x < z; otherwise split at z
    # (this achieves tighter fits on both children when z < x < w)
    splithere = min(n.x[i], problem.z[i])
    if verbose>=2 println("split on coordinate $i at $(n.x[i])") end

    # this does not correctly copy the model; left and right node contaminate each other
    # if model can only be reused for one node, it's more valuable to reuse it for the right node, since left node only has one constraint
    # right_m = copy(n.m)
    # set_optimizer(right_m, optimizer_with_attributes(GLPK.Optimizer,"tm_lim" => 60000, "msg_lev" => GLPK.GLP_MSG_OFF))
    # left_m = copy(n.m)
    # set_optimizer(left_m, optimizer_with_attributes(GLPK.Optimizer,"tm_lim" => 60000, "msg_lev" => GLPK.GLP_MSG_OFF))
    # left_m = n.m

    # left child
    #println("left child")
    left_u = copy(n.u)
    left_u[i] = splithere
    left_w = copy(n.w)
    left_w[i] = find_w(problem.fs[i],problem.dfs[i],n.l[i],left_u[i],problem.z[i])
    left = Node(n.l, left_u, left_w, problem; kwargs...)
    if left.maxdiff_index==-1
        return left, left
    end

    # right child
    #println("right child")
    right_l = copy(n.l)
    right_l[i] = splithere
    right_w = copy(n.w)
    right_w[i] = find_w(problem.fs[i],problem.dfs[i],right_l[i],n.u[i],problem.z[i])
    #right = Node(right_l, n.u, right_w, problem, Nothing, n.m; kwargs...)
    right = Node(right_l, n.u, right_w, problem; kwargs...)
    return left, right
end

## Branch and bound
function solve_sp(l, u, problem::SigmoidalProgram, init_x=Nothing;
                  TOL = 1e-2, maxiters = 100, verbose = 0, maxiters_noimprovement = Inf)
    subtol = TOL/length(l)/10
    root = Node(l, u, problem, init_x; TOL=subtol)
    bestnodes = Node[]
    ubs = Float64[]
    lbs = Float64[]
    pq = PriorityQueue{Node, Float64}(Reverse)
    if isnan(root.ub)
        println("Problem infeasible")
        return pq, bestnodes, lbs, ubs, 4
    elseif root.maxdiff_index==-1
        println("Solver time limit reached")
        return pq, bestnodes, lbs, ubs, 5
    end
    push!(bestnodes,root)
    push!(ubs,root.ub)
    push!(lbs,root.lb)
    enqueue!(pq, root, root.ub)
    #println("split maxiters: ", maxiters)
    for i=1:maxiters
        if verbose>=1
            println("iteration: ", i)
        end

        if ubs[end] - lbs[end] < TOL * lbs[end]
            println("found solution within tolerance $(ubs[end] - lbs[end]) in $i iterations")
            return pq, bestnodes, lbs, ubs, 0
        end
        if length(pq) > 0
            node = dequeue!(pq)
        else
            println("Stop iteration as node queue is empty at iteration $i.")
            return pq, bestnodes, lbs, ubs, 1
        end
        push!(ubs,min(node.ub, ubs[end]))
        left, right = split(node, problem; TOL=subtol)

        if left.maxdiff_index==-1 || right.maxdiff_index==-1
            println("Solver time limit reached")
            return pq, bestnodes, lbs, ubs, 5
        end

        if left.lb > lbs[end] && left.lb >= right.lb
            push!(lbs,left.lb)
            push!(bestnodes,left)
        elseif right.lb > lbs[end]
            push!(lbs,right.lb)
            push!(bestnodes,right)
        else
            push!(lbs,lbs[end])
        end
        if verbose>=2
            println("(lb, ub) = ($(lbs[end]), $(ubs[end]))")
        end

        # prune infeasible or obviously suboptimal nodes
        if !isnan(left.ub) && left.ub >= lbs[end]
            enqueue!(pq, left, left.ub)
            if verbose>=2 println("enqueued left") end
        else
            if verbose>=2 println("pruned left") end
        end
        if !isnan(right.ub) && right.ub >= lbs[end]
            enqueue!(pq, right, right.ub)
            if verbose>=2 println("enqueued right") end
        else
            if verbose>=2 println("pruned right") end
        end

        # break if no improvement for too long
        if length(lbs) > maxiters_noimprovement && lbs[end]==lbs[end-maxiters_noimprovement]
            println("Break after exceeding max iterations with no improvement at iteration $i.")
            return pq, bestnodes, lbs, ubs, 2
        end
    end
    println("Max iterations reached")
    return pq, bestnodes, lbs, ubs, 3
end