text:
using ProjAssistant
@quickactivate

## ------------------------------------------------------
@time begin
    import Chemostat_InSilico
    const Dyn = Chemostat_InSilico.Dynamic
    import Chemostat_InSilico.Dynamic:
        maxent, normalizeP!, load_simdat, get_status,
        UNDONE_SIM_STATUS, DEAD_SIM_STATUS,
        EXPLODED_SIM_STATUS, NITERS_SIM_STATUS,
        STST_SIM_STATUS, collect_ts

    using Base.Threads
    using ExtractMacro

end

## ----------------------------------------------------------------------------
# ME_AVZ_EXPECTED_AVUG_BOUNDED
let
    MEmode = "ME_AVZ_EXPECTED_AVUG_BOUNDED"

    # batch = (; push_frec, Xts, sgts, z_avts, ug_avts, cgD_Xts, Pzts, Pugts)
    simid = :SimD3
    params = lglob(Dyn, simid, :params, :finite_cg)
    @extract params: Ds ϵs cg

    _simparams_pool = Iterators.product(Ds, ϵs, cg)
    tot_sims = length(_simparams_pool)
    Ch = Channel() do Ch_
        for (simi, (D, ϵ, cg)) in enumerate(_simparams_pool)
            put!(Ch_, (simi, (D, ϵ, cg)))
        end
    end

    @threads for _ in 1:nthreads()
        thid = threadid()

        for (simi, (D, ϵ, cg)) in Ch
            simparams = (;D, ϵ, cg)

            # ---------------------------------------------------------------
            # Welcome && status && cache
            status = Dyn.get_status(simid, (;D, ϵ, cg))
            datfile = Dyn.maxent_file(simid, MEmode, simparams)

            @info("At", simid, simparams, thid, status, relpath(datfile))
            isfile(datfile) && continue
            (status != STST_SIM_STATUS) && continue

            # ---------------------------------------------------------------
            # Sim globals
            S = Dyn.load_simdat(simid, simparams)
            X = S.X
            cgD_X = S.cg * S.D / S.X
            V = S.V
            S = nothing; GC.gc()

            # ---------------------------------------------------------------
            # maxent
            me_out = maxent(V, D, cgD_X;
                # Top loop
                top_maxiter = 100,

                # GD globals
                gd_maxiter = 1000,
                gdth = 1e-2,
                verbose = true,

                # stst
                stst_th = 1e-2,
                stst_w = 10,

                # info
                further_info = (;
                    progress = string(simi, " of ", tot_sims, " (", div(simi * 100, tot_sims), "%)")
                )
            )

            sdat(Dyn, me_out, datfile)

        end # for (simi, (D, ϵ, cg))

    end # for thid
end


text:
function add_loads_or_constraint_transforms!(m::Model, variablesin)

end

function logdensity_model(modelin::Model, variabledescription)
    m = deepcopy(modelin)
    @unpack vars = m
    nvars = length(vars)
    θ = addvar!(m, Symbol("##θ##")) # parameter vector
    for i ∈ 1:nvars
        var = vars[i]
        descript = variabledescription[i]
        if istracked(descript)# load from Symbol("##θ##")
            # Plan is to define func `constrain` taking in type of the variable, model, and Symbol("##θ##"), returning that variable.

            var.initialized = false
        else# load from Symbol("##DATA##")

        end
    end

end

ext: jl


text:
import Blosc

# port of https://github.com/Blosc/c-blosc/blob/3a668dcc9f61ad22b5c0a0ab45fe8dad387277fd/hdf5/blosc_filter.c (copyright 2010 Francesc Alted, license: MIT/expat)

typealias H5T_class_t Cint # C enum

# The following definitions mirror those in H5Zpublic.h for HDF5 1.8.x:
const H5Z_FLAG_REVERSE = 0x0100
const H5Z_FLAG_OPTIONAL = 0x0001
const H5Z_CLASS_T_VERS = 1
typealias H5Z_filter_t Cint
immutable H5Z_class2_t
    version::Cint # = H5Z_CLASS_T_VERS
    id::H5Z_filter_t # Filter ID number
    encoder_present::Cuint # Does this filter have an encoder?
    decoder_present::Cuint # Does this filter have a decoder?
    name::Ptr{UInt8} # Comment for debugging
    can_apply::Ptr{Void} # The "can apply" callback
    set_local::Ptr{Void} # The "set local" callback
    filter::Ptr{Void} # The filter callback
end

const FILTER_BLOSC_VERSION = 2
const FILTER_BLOSC = 32001 # Filter ID registered with the HDF Group for Blosc
const blosc_name = "blosc"

const blosc_flags_ = Array(Cuint, 1)
const blosc_nelements_ = Array(Csize_t, 1)
const blosc_values = Array(Cuint, 8)
const blosc_chunkdims = Array(Hsize, 32)
function blosc_set_local(dcpl::Hid, htype::Hid, space::Hid)
    blosc_nelements_[1] = 8
    if ccall((:H5Pget_filter_by_id2,libhdf5), Herr,
             (Hid, H5Z_filter_t, Ptr{Cuint}, Ptr{Csize_t}, Ptr{Cuint},
              Csize_t, Ptr{UInt8}, Ptr{Cuint}),
             dcpl, FILTER_BLOSC, blosc_flags_, blosc_nelements_, blosc_values,
             0, C_NULL, C_NULL) < 0
        return convert(Herr, -1)
    end
    flags = blosc_flags_[1]
    nelements = max(blosc_nelements_[1], 4)
    # Set Blosc info in first two slots
    blosc_values[1] = FILTER_BLOSC_VERSION;
    blosc_values[2] = Blosc.VERSION_FORMAT;
    ndims = ccall((:H5Pget_chunk,libhdf5), Cint, (Hid, Cint, Ptr{Hsize}),
                  dcpl, 32, blosc_chunkdims);
    chunksize = prod(resize!(blosc_chunkdims, ndims))
    if ndims < 0 || ndims > 32 || chunksize > Blosc.MAX_BUFFERSIZE
        return convert(Herr, -1)
    end
    htypesize = ccall((:H5Tget_size,libhdf5), Csize_t, (Hid,), htype)
    htypesize == 0 && returnconvert(Herr, -1)
    if ccall((:H5Tget_class,libhdf5),H5T_class_t,(Hid,),htype) == H5T_ARRAY
        hsuper = ccall((:H5Tget_super,libhdf5), Hid, (Hid,), htype)
        basetypesize = ccall((:H5Tget_size,libhdf5),Csize_t,(Hid,), hsuper)
        ccall((:H5Tclose,libhdf5), Herr, (Hid,), hsuper)
    else
        basetypesize = htypesize
        end
    # Limit large typesizes (they are pretty inefficient to shuffle
    # and, in addition, Blosc does not handle typesizes larger than
    # blocksizes).
    if basetypesize > Blosc.MAX_TYPESIZE
        basetypesize = 1
    end
    blosc_values[3] = basetypesize
    blosc_values[4] = chunksize * htypesize # size of the chunk
    if ccall((:H5Pmodify_filter,libhdf5), Herr,
             (Hid, H5Z_filter_t, Cuint, Csize_t, Ptr{Cuint}),
             dcpl, FILTER_BLOSC, flags, nelements, blosc_values) < 0
        return convert(Herr, -1)
        end
    return convert(Herr, 1)
end

function blosc_filter(flags::Cuint, cd_nelmts::Csize_t,
                      cd_values::Ptr{Cuint}, nbytes::Csize_t,
                      buf_size::Ptr{Csize_t}, buf::Ptr{Ptr{Void}})
    typesize = unsafe_load(cd_values, 3) # The datatype size
    outbuf_size = unsafe_load(cd_values, 4)
    # Compression level:
    clevel = cd_nelmts >= 5 ? unsafe_load(cd_values, 5) : convert(Cuint, 5)
    # Do shuffle:
    doshuffle = cd_nelmts >= 6 ? unsafe_load(cd_values, 6) != 0 : true
    # to do: set compressor based on compcode in unsafe_load(cd_values, 7)?

    if (flags & H5Z_FLAG_REVERSE) == 0 # compressing
        # Allocate an output buffer exactly as long as the input data; if
        # the result is larger, we simply return 0. The filter is flagged
        # as optional, so HDF5 marks the chunk as uncompressed and proceeds.
        outbuf_size = unsafe_load(buf_size)
        outbuf = Libc.malloc(outbuf_size)
        outbuf == C_NULL && return convert(Csize_t, 0)
        status = Blosc.blosc_compress(clevel, doshuffle, typesize, nbytes,
                                      unsafe_load(buf), outbuf, nbytes)
        status < 0 && (Libc.free(outbuf); return convert(Csize_t, 0))
    else # decompressing
        # Extract the exact outbuf_size from the buffer header.
        #
        # NOTE: the guess value got from "cd_values" corresponds to the
        # uncompressed chunk size but it should not be used in a general
        # cases since other filters in the pipeline can modify the buffer
        # size.
        outbuf_size, cbytes, blocksize = Blosc.cbuffer_sizes(unsafe_load(buf))
        outbuf = Libc.malloc(outbuf_size)
        outbuf == C_NULL && return convert(Csize_t, 0)
        status = Blosc.blosc_decompress(unsafe_load(buf), outbuf, outbuf_size)
        status <= 0 && (Libc.free(outbuf); return convert(Csize_t, 0))
    end

    if status != 0
        Libc.free(unsafe_load(buf));
        unsafe_store!(buf, outbuf)
        unsafe_store!(buf_size, outbuf_size)
        return convert(Csize_t, status) # size of compressed/decompressed data
    end
    Libc.free(outbuf); return convert(Csize_t, 0)
end

# register the Blosc filter function with HDF5
function register_blosc()
    c_blosc_set_local = cfunction(blosc_set_local, Herr, (Hid,Hid,Hid))
    c_blosc_filter = cfunction(blosc_filter, Csize_t,
                                 (Cuint, Csize_t, Ptr{Cuint}, Csize_t,
                                  Ptr{Csize_t}, Ptr{Ptr{Void}}))
    if ccall((:H5Zregister, libhdf5), Herr, (Ptr{H5Z_class2_t},),
             &H5Z_class2_t(H5Z_CLASS_T_VERS,
                           FILTER_BLOSC,
                           1, 1,
                           pointer(blosc_name),
                           C_NULL,
                           c_blosc_set_local,
                           c_blosc_filter)) < 0
        error("can't register Blosc filter")
    end
end

const _set_blosc_values = Cuint[0,0,0,0, 5,1,0]
function h5p_set_blosc(p::HDF5Properties, level::Integer=5)
    0 <= level <= 9 || throw(ArgumentError("blosc compression $level not in [0,9]"))
    _set_blosc_values[5] = level
    status = ccall((:H5Pset_filter,libhdf5), Herr,
          (Hid, H5Z_filter_t, Cuint, Csize_t, Ptr{Cuint}),
          p.id, FILTER_BLOSC, H5Z_FLAG_OPTIONAL, 7, _set_blosc_values)
    status < 0 && error("Error setting blosc compression level")
    nothing
end


text:
import .YaoBlocks: _apply!
include("shortcuts.jl")
include("FSimGate.jl")


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# src/eago_optimizer/relax.jl
# Defines routines used construct the relaxed subproblem.
#############################################################################

"""
$(FUNCTIONNAME)

Applies the safe cut checks detailed in Khajavirad, 2018 [Khajavirad, Aida,
and Nikolaos V. Sahinidis. "A hybrid LP/NLP paradigm for global optimization
relaxations." Mathematical Programming Computation 10.3 (2018): 383-421] to
ensure that only numerically safe affine relaxations are added. Checks that
i) ``|b| <= safe b`, ii) `safe_l <= abs(ai) <= safe u`, and iii) violates
`safe_l <= abs(ai/aj) <= safe_u`.
"""
function is_safe_cut!(m::GlobalOptimizer, f::SAF)

    safe_l = m._parameters.cut_safe_l
    safe_u = m._parameters.cut_safe_u
    safe_b = m._parameters.cut_safe_b

    (abs(f.constant) > safe_b) && return false # violates |b| <= safe_b
    term_count = length(f.terms)
    @inbounds for i = 1:term_count
        ai = f.terms[i].coefficient
        if ai !== 0.0
            if !(safe_l <= abs(ai) <= safe_u)  # violates safe_l <= abs(ai) <= safe_u
                return false
            end
            @inbounds for j = i:term_count     # violates safe_l <= abs(ai/aj) <= safe_u
                aj = f.terms[j].coefficient
                if aj !== 0.0
                    if !(safe_l <= abs(ai/aj) <= safe_u)
                        return false
                    end
                end
            end
        end
    end
    return true
end

function add_affine_relaxation!(m::GlobalOptimizer{R,S,Q}, f::SAF, check_safe::Bool) where {R,S,Q<:ExtensionType}
    valid_cut_flag = !check_safe || is_safe_cut!(m, f)
    if valid_cut_flag
        s = LT(-f.constant + _constraint_tol(m))
        f.constant = 0.0
        ci = MOI.add_constraint(_relaxed_optimizer(m), f, s)::CI{SAF,LT}
        push!(m._affine_relax_ci, ci)
    end
    return valid_cut_flag
end

"""
$(FUNCTIONNAME)

Relaxs the constraint by adding an affine constraint to the model.
"""
function relax! end

"""
$(FUNCTIONNAME)

Default routine for relaxing quadratic constraint `func < 0.0` on node `n`.
Takes affine bounds of convex part at point `x0` and secant line bounds on
concave parts.
"""
function affine_relax_quadratic!(m::GlobalOptimizer, func::SQF, buffer::Dict{Int,Float64}, saf::SAF) where {R,S,Q<:ExtensionType}

    quadratic_constant = func.constant

    # Affine terms only contribute coefficients, so the respective
    # values do not contribute to the cut. Since all quadratic terms
    # are considered to be branch variables we exclude any potential
    # need to retrieve variable bounds from locations other than
    # the node.
    for term in func.quadratic_terms
        a = term.coefficient
        i = term.variable_1.value
        j = term.variable_2.value
        x0 = _lower_solution(FullVar(), m, i)
        xL = _lower_bound(FullVar(), m, i)
        xU = _upper_bound(FullVar(), m, i)
        if i == j
            if a > 0.0
                buffer[i] += a*x0
                quadratic_constant -= 0.5*a*x0*x0
            else
                if !isinf(xL) && !isinf(xU)
                    buffer[i] += 0.5*a*(xL + xU)
                    quadratic_constant -= 0.5*a*xL*xU
                else
                    return false
                end
            end
        else
            y0 = _lower_solution(FullVar(), m, j)
            yL = _lower_bound(FullVar(), m, j)
            yU = _upper_bound(FullVar(), m, j)
            if a > 0.0
                if (!isinf(xL) && !isinf(yL)) && ((xU - xL)*y0 + (yU - yL)*x0 <= xU*yU - xL*yL)
                    buffer[i] += a*yL
                    buffer[j] += a*xL
                    quadratic_constant -= a*xL*yL

                elseif !isinf(xU) && !isinf(yU)
                    buffer[i] += a*yU
                    buffer[j] += a*xU
                    quadratic_constant -= a*xU*yU
                else
                    return false
                end
            else
                if (!isinf(xU) && !isinf(yL)) && ((xU - xL)*y0 - (yU - yL)*x0 <= xU*yL - xL*yU)
                    buffer[i] += a*yL
                    buffer[j] += a*xU
                    quadratic_constant -= a*xU*yL
                elseif !isinf(xL) && !isinf(yU)
                    buffer[i] += a*yU
                    buffer[j] += a*xL
                    quadratic_constant -= a*xL*yU
                else
                    return false
                end
            end
        end
    end

    for t in func.affine_terms
        buffer[t.variable.value] += t.coefficient
    end
    count = 1
    for (key, value) in buffer
        saf.terms[count] = SAT(value, VI(key))
        buffer[key] = 0.0
        count += 1
    end
    saf.constant = quadratic_constant
    return
end

"""
$(TYPEDSIGNATURES)
"""
function relax!(m::GlobalOptimizer, f::BQI, k::Int, check_safe::Bool)
    affine_relax_quadratic!(m, f.func, f.buffer, f.saf)
    valid_cut_flag = add_affine_relaxation!(m, f.saf, check_safe)
    return valid_cut_flag
end

"""
$(TYPEDSIGNATURES)
"""
function relax!(m::GlobalOptimizer, f::BQE, i::Int, check_safe::Bool)
    affine_relax_quadratic!(m, f.func, f.buffer, f.saf)
    valid_cut_flag = add_affine_relaxation!(m, f.saf, check_safe)
    affine_relax_quadratic!(m, f.minus_func, f.buffer, f.saf)
    valid_cut_flag &= add_affine_relaxation!(m, f.saf, check_safe)
    return valid_cut_flag
end

"""
$(TYPEDSIGNATURES)
"""
function check_set_affine_nl!(m::GlobalOptimizer{R,S,Q}, f::BufferedNonlinearFunction{V,N,T}, finite_cut::Bool, check_safe::Bool) where {V,R,S,N,T<:RelaxTag,Q<:ExtensionType}
    valid_cut_flag = finite_cut && (!check_safe || is_safe_cut!(m, f.saf))
    if valid_cut_flag
        lt = LT(_constraint_tol(m) - f.saf.constant)
        f.saf.constant = 0.0
        push!(m._affine_relax_ci, MOI.add_constraint(_relaxed_optimizer(m), f.saf, lt))
    end
    return valid_cut_flag
end

"""
$(TYPEDSIGNATURES)
"""
function relax!(m::GlobalOptimizer{R,S,Q}, f::BufferedNonlinearFunction{V,N,T}, k::Int, check_safe::Bool) where {V,R,S,N,T<:RelaxTag,Q<:ExtensionType}
    d = m._working_problem._relaxed_evaluator
    x = d.variable_values.x
    d.pass_number = k
    forward_pass!(d, f)
    valid_cut_flag = true
    num_feasible_cut = true

    grad_sparsity = sparsity(f)
    if is_num(f)
        f.saf.constant = num(f)
        for i = 1:length(grad_sparsity)
            f.saf.terms[i] = SAT(0.0, VI(grad_sparsity[i]))
        end
        if !(lower_bound(f) <= num(f) <= upper_bound(f))
            num_feasible_cut = false
        end
    else
        v = set(f)
        if !isempty(v)
            # if has less than or equal to bound (<=)
            if isfinite(upper_bound(f))
                lower_cut_valid = !isnan(v.cv) && isfinite(v.cv)
                if lower_cut_valid
                    f.saf.constant = v.cv - upper_bound(f)
                    for (i, k) in enumerate(grad_sparsity)
                        c = v.cv_grad[i]
                        f.saf.terms[i] = SAT(c, VI(k))
                        f.saf.constant -= c*x[k]
                    end
                    valid_cut_flag = check_set_affine_nl!(m, f, lower_cut_valid, check_safe)
                end
            end
            # if has greater than or equal to bound (>=)
            if isfinite(lower_bound(f))
                upper_cut_valid = !isnan(v.cc) && isfinite(v.cc)
                if upper_cut_valid
                    f.saf.constant = -v.cc + lower_bound(f)
                    for (i, k) in enumerate(grad_sparsity)
                        c = -v.cc_grad[i]
                        f.saf.terms[i] = SAT(c, VI(k))
                        f.saf.constant -= c*x[k]
                    end
                    valid_cut_flag &= check_set_affine_nl!(m, f, upper_cut_valid, check_safe)
                end
            end
        else
            num_feasible_cut = false
        end
    end
    return valid_cut_flag, num_feasible_cut
end

relax!(m::GlobalOptimizer, f::Union{Nothing, VI, AFI}, k::Int, b::Bool) = true

max_stars_repo_name: MatthewStuber/EAGO

max_issues_repo_name: MatthewStuber/EAGO

max_forks_repo_name: MatthewStuber/EAGO


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

export
    abspath,
    basename,
    dirname,
    expanduser,
    contractuser,
    homedir,
    isabspath,
    isdirpath,
    joinpath,
    normpath,
    realpath,
    relpath,
    splitdir,
    splitdrive,
    splitext,
    splitpath

if Sys.isunix()
    const path_separator    = "/"
    const path_separator_re = r"/+"
    const path_directory_re = r"(?:^|/)\.{0,2}$"
    const path_dir_splitter = r"^(.*?)(/+)([^/]*)$"
    const path_ext_splitter = r"^((?:.*/)?(?:\.|[^/\.])[^/]*?)(\.[^/\.]*|)$"

    splitdrive(path::String) = ("",path)
elseif Sys.iswindows()
    const path_separator    = "\\"
    const path_separator_re = r"[/\\]+"
    const path_absolute_re  = r"^(?:[A-Za-z]+:)?[/\\]"
    const path_directory_re = r"(?:^|[/\\])\.{0,2}$"
    const path_dir_splitter = r"^(.*?)([/\\]+)([^/\\]*)$"
    const path_ext_splitter = r"^((?:.*[/\\])?(?:\.|[^/\\\.])[^/\\]*?)(\.[^/\\\.]*|)$"

    function splitdrive(path::String)
        m = match(r"^([^\\]+:|\\\\[^\\]+\\[^\\]+|\\\\\?\\UNC\\[^\\]+\\[^\\]+|\\\\\?\\[^\\]+:|)(.*)$"s, path)
        String(m.captures[1]), String(m.captures[2])
    end
else
    error("path primitives for this OS need to be defined")
end

"""
    splitdrive(path::AbstractString) -> (AbstractString, AbstractString)

On Windows, split a path into the drive letter part and the path part. On Unix systems, the
first component is always the empty string.
"""
splitdrive(path::AbstractString)

"""
    homedir() -> String

Return the current user's home directory.

!!! note
    `homedir` determines the home directory via `libuv`'s `uv_os_homedir`. For details
    (for example on how to specify the home directory via environment variables), see the
    [`uv_os_homedir` documentation](http://docs.libuv.org/en/v1.x/misc.html#c.uv_os_homedir).
"""
function homedir()
    buf = Base.StringVector(AVG_PATH - 1) # space for null-terminator implied by StringVector
    sz = RefValue{Csize_t}(length(buf) + 1) # total buffer size including null
    while true
        rc = ccall(:uv_os_homedir, Cint, (Ptr{UInt8}, Ptr{Csize_t}), buf, sz)
        if rc == 0
            resize!(buf, sz[])
            return String(buf)
        elseif rc == Base.UV_ENOBUFS
            resize!(buf, sz[] - 1) # space for null-terminator implied by StringVector
        else
            uv_error("homedir()", rc)
        end
    end
end

if Sys.iswindows()
    isabspath(path::AbstractString) = occursin(path_absolute_re, path)
else
    isabspath(path::AbstractString) = startswith(path, '/')
end

"""
    isabspath(path::AbstractString) -> Bool

Determine whether a path is absolute (begins at the root directory).

# Examples
```jldoctest
julia> isabspath("/home")
true

julia> isabspath("home")
false
```
"""
isabspath(path::AbstractString)

"""
    isdirpath(path::AbstractString) -> Bool

Determine whether a path refers to a directory (for example, ends with a path separator).

# Examples
```jldoctest
julia> isdirpath("/home")
false

julia> isdirpath("/home/")
true
```
"""
isdirpath(path::String) = occursin(path_directory_re, splitdrive(path)[2])

"""
    splitdir(path::AbstractString) -> (AbstractString, AbstractString)

Split a path into a tuple of the directory name and file name.

# Examples
```jldoctest
julia> splitdir("/home/myuser")
("/home", "myuser")
```
"""
function splitdir(path::String)
    a, b = splitdrive(path)
    _splitdir_nodrive(a,b)
end

# Common splitdir functionality without splitdrive, needed for splitpath.
_splitdir_nodrive(path::String) = _splitdir_nodrive("", path)
function _splitdir_nodrive(a::String, b::String)
    m = match(path_dir_splitter,b)
    m === nothing && return (a,b)
    cs = m.captures
    getcapture(cs, i) = cs[i]::AbstractString
    c1, c2, c3 = getcapture(cs, 1), getcapture(cs, 2), getcapture(cs, 3)
    a = string(a, isempty(c1) ? c2[1] : c1)
    a, String(c3)
end

"""
    dirname(path::AbstractString) -> AbstractString

Get the directory part of a path. Trailing characters ('/' or '\\') in the path are
counted as part of the path.

# Examples
```jldoctest
julia> dirname("/home/myuser")
"/home"

julia> dirname("/home/myuser/")
"/home/myuser"
```

See also: [`basename`](@ref)
"""
 dirname(path::AbstractString) = splitdir(path)[1]

"""
    basename(path::AbstractString) -> AbstractString

Get the file name part of a path.

!!! note
    This function differs slightly from the Unix `basename` program, where trailing slashes are ignored,
    i.e. `\$ basename /foo/bar/` returns `bar`, whereas `basename` in Julia returns an empty string `""`.

# Examples
```jldoctest
julia> basename("/home/myuser/example.jl")
"example.jl"

julia> basename("/home/myuser/")
""
```

See also: [`dirname`](@ref)
"""
basename(path::AbstractString) = splitdir(path)[2]

"""
    splitext(path::AbstractString) -> (AbstractString, AbstractString)

If the last component of a path contains a dot, split the path into everything before the
dot and everything including and after the dot. Otherwise, return a tuple of the argument
unmodified and the empty string. "splitext" is short for "split extension".

# Examples
```jldoctest
julia> splitext("/home/myuser/example.jl")
("/home/myuser/example", ".jl")

julia> splitext("/home/myuser/example")
("/home/myuser/example", "")
```
"""
function splitext(path::String)
    a, b = splitdrive(path)
    m = match(path_ext_splitter, b)
    m === nothing && return (path,"")
    a*m.captures[1], String(m.captures[2])
end

# NOTE: deprecated in 1.4
pathsep() = path_separator

"""
    splitpath(path::AbstractString) -> Vector{String}

Split a file path into all its path components. This is the opposite of
`joinpath`. Returns an array of substrings, one for each directory or file in
the path, including the root directory if present.

!!! compat "Julia 1.1"
    This function requires at least Julia 1.1.

# Examples
```jldoctest
julia> splitpath("/home/myuser/example.jl")
4-element Vector{String}:
 "/"
 "home"
 "myuser"
 "example.jl"
```
"""
splitpath(p::AbstractString) = splitpath(String(p))

function splitpath(p::String)
    drive, p = splitdrive(p)
    out = String[]
    isempty(p) && (pushfirst!(out,p))  # "" means the current directory.
    while !isempty(p)
        dir, base = _splitdir_nodrive(p)
        dir == p && (pushfirst!(out, dir); break)  # Reached root node.
        if !isempty(base)  # Skip trailing '/' in basename
            pushfirst!(out, base)
        end
        p = dir
    end
    if !isempty(drive)  # Tack the drive back on to the first element.
        out[1] = drive*out[1]  # Note that length(out) is always >= 1.
    end
    return out
end

joinpath(path::AbstractString)::String = path

if Sys.iswindows()

function joinpath(path::AbstractString, paths::AbstractString...)::String
    result_drive, result_path = splitdrive(path)

    local p_drive, p_path
    for p in paths
        p_drive, p_path = splitdrive(p)

        if startswith(p_path, ('\\', '/'))
            # second path is absolute
            if !isempty(p_drive) || !isempty(result_drive)
                result_drive = p_drive
            end
            result_path = p_path
            continue
        elseif !isempty(p_drive) && p_drive != result_drive
            if lowercase(p_drive) != lowercase(result_drive)
                # different drives, ignore the first path entirely
                result_drive = p_drive
                result_path = p_path
                continue
            end
        end

        # second path is relative to the first
        if !isempty(result_path) && result_path[end] ∉ ('\\', '/')
            result_path *= "\\"
        end

        result_path = result_path * p_path
    end

    # add separator between UNC and non-absolute path
    if !isempty(p_path) && result_path[1] ∉ ('\\', '/') && !isempty(result_drive) && result_drive[end] != ':'
        return result_drive * "\\" * result_path
    end

    return result_drive * result_path
end

else

function joinpath(path::AbstractString, paths::AbstractString...)::String
    for p in paths
        if isabspath(p)
            path = p
        elseif isempty(path) || path[end] == '/'
            path *= p
        else
            path *= "/" * p
        end
    end
    return path
end

end # os-test

"""
    joinpath(parts::AbstractString...) -> String

Join path components into a full path. If some argument is an absolute path or
(on Windows) has a drive specification that doesn't match the drive computed for
the join of the preceding paths, then prior components are dropped.

Note on Windows since there is a current directory for each drive, `joinpath("c:", "foo")`
represents a path relative to the current directory on drive "c:" so this is equal to "c:foo",
not "c:\\foo". Furthermore, `joinpath` treats this as a non-absolute path and ignores the drive
letter casing, hence `joinpath("C:\\A","c:b") = "C:\\A\\b"`.

# Examples
```jldoctest
julia> joinpath("/home/myuser", "example.jl")
"/home/myuser/example.jl"
```
"""
joinpath

"""
    normpath(path::AbstractString) -> String

Normalize a path, removing "." and ".." entries.

# Examples
```jldoctest
julia> normpath("/home/myuser/../example.jl")
"/home/example.jl"
```
"""
function normpath(path::String)
    isabs = isabspath(path)
    isdir = isdirpath(path)
    drive, path = splitdrive(path)
    parts = split(path, path_separator_re)
    filter!(x->!isempty(x) && x!=".", parts)
    while true
        clean = true
        for j = 1:length(parts)-1
            if parts[j] != ".." && parts[j+1] == ".."
                deleteat!(parts, j:j+1)
                clean = false
                break
            end
        end
        clean && break
    end
    if isabs
        while !isempty(parts) && parts[1] == ".."
            popfirst!(parts)
        end
    elseif isempty(parts)
        push!(parts, ".")
    end
    path = join(parts, path_separator)
    if isabs
        path = path_separator*path
    end
    if isdir && !isdirpath(path)
        path *= path_separator
    end
    string(drive,path)
end

"""
    normpath(path::AbstractString, paths::AbstractString...) -> String

Convert a set of paths to a normalized path by joining them together and removing
"." and ".." entries. Equivalent to `normpath(joinpath(path, paths...))`.
"""
normpath(a::AbstractString, b::AbstractString...) = normpath(joinpath(a,b...))

"""
    abspath(path::AbstractString) -> String

Convert a path to an absolute path by adding the current directory if necessary.
Also normalizes the path as in [`normpath`](@ref).
"""
abspath(a::String) = normpath(isabspath(a) ? a : joinpath(pwd(),a))

"""
    abspath(path::AbstractString, paths::AbstractString...) -> String

Convert a set of paths to an absolute path by joining them together and adding the
current directory if necessary. Equivalent to `abspath(joinpath(path, paths...))`.
"""
abspath(a::AbstractString, b::AbstractString...) = abspath(joinpath(a,b...))

if Sys.iswindows()

function longpath(path::AbstractString)
    p = cwstring(path)
    buf = zeros(UInt16, length(p))
    while true
        n = ccall((:GetLongPathNameW, "kernel32"), stdcall,
            UInt32, (Ptr{UInt16}, Ptr{UInt16}, UInt32),
            p, buf, length(buf))
        windowserror(:longpath, n == 0)
        x = n < length(buf) # is the buffer big enough?
        resize!(buf, n) # shrink if x, grow if !x
        x && return transcode(String, buf)
    end
end

end # os-test

"""
    realpath(path::AbstractString) -> String

Canonicalize a path by expanding symbolic links and removing "." and ".." entries.
On case-insensitive case-preserving filesystems (typically Mac and Windows), the
filesystem's stored case for the path is returned.

(This function throws an exception if `path` does not exist in the filesystem.)
"""
function realpath(path::AbstractString)
    req = Libc.malloc(_sizeof_uv_fs)
    try
        ret = ccall(:uv_fs_realpath, Cint,
                    (Ptr{Cvoid}, Ptr{Cvoid}, Cstring, Ptr{Cvoid}),
                    C_NULL, req, path, C_NULL)
        if ret < 0
            ccall(:uv_fs_req_cleanup, Cvoid, (Ptr{Cvoid},), req)
            uv_error("realpath($(repr(path)))", ret)
        end
        path = unsafe_string(ccall(:jl_uv_fs_t_ptr, Cstring, (Ptr{Cvoid},), req))
        ccall(:uv_fs_req_cleanup, Cvoid, (Ptr{Cvoid},), req)
        return path
    finally
        Libc.free(req)
    end
end

if Sys.iswindows()
# on windows, ~ means "temporary file"
expanduser(path::AbstractString) = path
contractuser(path::AbstractString) = path
else
function expanduser(path::AbstractString)
    y = iterate(path)
    y === nothing && return path
    c, i = y::Tuple{eltype(path),Int}
    c != '~' && return path
    y = iterate(path, i)
    y === nothing && return homedir()
    y[1]::eltype(path) == '/' && return homedir() * path[i:end]
    throw(ArgumentError("~user tilde expansion not yet implemented"))
end
function contractuser(path::AbstractString)
    home = homedir()
    if path == home
        return "~"
    elseif startswith(path, home)
        return joinpath("~", relpath(path, home))
    else
        return path
    end
end
end

"""
    expanduser(path::AbstractString) -> AbstractString

On Unix systems, replace a tilde character at the start of a path with the current user's home directory.
"""
expanduser(path::AbstractString)

"""
    contractuser(path::AbstractString) -> AbstractString

On Unix systems, if the path starts with `homedir()`, replace it with a tilde character.
"""
contractuser(path::AbstractString)

"""
    relpath(path::AbstractString, startpath::AbstractString = ".") -> AbstractString

Return a relative filepath to `path` either from the current directory or from an optional
start directory. This is a path computation: the filesystem is not accessed to confirm the
existence or nature of `path` or `startpath`.

On Windows, case sensitivity is applied to every part of the path except drive letters. If
`path` and `startpath` refer to different drives, the absolute path of `path` is returned.
"""
function relpath(path::String, startpath::String = ".")
    isempty(path) && throw(ArgumentError("`path` must be specified"))
    isempty(startpath) && throw(ArgumentError("`startpath` must be specified"))
    curdir = "."
    pardir = ".."
    path == startpath && return curdir
    path_drive, path_without_drive = splitdrive(path)
    startpath_drive, startpath_without_drive = splitdrive(startpath)
    path_arr  = split(abspath(path_without_drive),      path_separator_re)
    start_arr = split(abspath(startpath_without_drive), path_separator_re)
    if Sys.iswindows()
        lowercase(path_drive) != lowercase(startpath_drive) && return abspath(path)
    end
    i = 0
    while i < min(length(path_arr), length(start_arr))
        i += 1
        if path_arr[i] != start_arr[i]
            i -= 1
            break
        end
    end
    pathpart = join(path_arr[i+1:something(findlast(x -> !isempty(x), path_arr), 0)], path_separator)
    prefix_num = something(findlast(x -> !isempty(x), start_arr), 0) - i - 1
    if prefix_num >= 0
        prefix = pardir * path_separator
        relpath_ = isempty(pathpart)     ?
            (prefix^prefix_num) * pardir :
            (prefix^prefix_num) * pardir * path_separator * pathpart
    else
        relpath_ = pathpart
    end
    return isempty(relpath_) ? curdir :  relpath_
end
relpath(path::AbstractString, startpath::AbstractString) =
    relpath(String(path), String(startpath))

for f in (:isdirpath, :splitdir, :splitdrive, :splitext, :normpath, :abspath)
    @eval $f(path::AbstractString) = $f(String(path))
end


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

include("weighting.jl")
include("bat_sample.jl")
include("mcmc/mcmc.jl")


text:
"""
    Binion{T <: Real} = Elliptic{T}

A binion is an elliptic Cayley-Dickson construct with real numbers.
"""
const Binion{T <: Real} = Elliptic{T}

function Binion(a::T, b::T) where T <: Real
    Binion{T}(a, b)
end

function Binion(a::T) where T <: Real
    Binion{T}(a, zero(T))
end

function Binion(a::Real, b::Real)
    Binion(promote(a, b)...)
end

"""
    unreal(z::Binion)

The unreal part of a complex number is a 1-vector.
"""
function unreal(z::Binion{T}) where T <: Real
    T[z.r]
end

function show(io::IO, z::Binion)
    print(io, "[1: ")
    print(io, z.l)
    print(io, ", i: ")
    print(io, z.r)
    print(io, "]")
end

function random(::Type{Binion{T}}) where T <: Real
    Binion{T}(rand(T), rand(T))
end

function random(::Type{Binion{BigInt}})
    Binion{BigInt}(randomBigInt(), randomBigInt())
end

function random(::Type{Binion{BigFloat}})
    Binion{BigFloat}(randomBigFloat(), randomBigFloat())
end

function random(::Type{Binion{Rational{BigInt}}})
    Binion{Rational{BigInt}}(randomBigRational(), randomBigRational())
end

function convert(::Type{Binion{T}}, a::Real) where T <: Real
    Binion{T}(a, zero(T))
end

function convert(::Type{Binion{T}}, z::Binion) where T <: Real
    Binion{T}(z.l, z.r)
end

function promote_rule(::Type{Binion{T}}, ::Type{S}) where {T <: Real, S <: Real}
    Binion{promote_type(T, S)}
end

function promote_rule(::Type{Binion{T}}, ::Type{Binion{S}}) where {T <: Real, S <: Real}
    Binion{promote_type(T, S)}
end

function selfstar(z::Binion)
    (z + star(z)) / 2
end

function anti_selfstar(z::Binion)
    (z - star(z)) / 2
end


text:
a = 20

if a > 10
a = a - 5
end

println(a)


text:
import VideoIO
using GLPlot, React, GLAbstraction

device  = VideoIO.DEFAULT_CAMERA_DEVICE
format  = VideoIO.DEFAULT_CAMERA_FORMAT
camera  = VideoIO.opencamera(device, format)
img   = VideoIO.read(camera)

# Just for fun, lets apply a laplace filter:

#async=true, for REPL use. Then you don't have to call renderloop(window)
window  = createdisplay(#=async=true =#)
img     = glplot(Texture(img, 3))

#Get Gpu memory object
glimg   = img.uniforms[:image]

#Asynchronous updating with React:
lift(Timing.fpswhen(window.inputs[:open], 30.0)) do x
  newframe = VideoIO.read(camera)
  update!(glimg,  mapslices(reverse, newframe, 3)) # needs to be mirrored :(
end

renderloop(window)


text:
# Copyright(c) 2020: Matthew Wilhelm & Matthew Stuber.
# This work is licensed under the Creative Commons Attribution-NonCommercial-
# ShareAlike 4.0 International License. To view a copy of this license, visit
# http://creativecommons.org/licenses/by-nc-sa/4.0/ or send a letter to Creative
# Commons, PO Box 1866, Mountain View, CA 94042, USA.
#############################################################################
# Dynamic Bounds - pODEs Discrete
# A package for discretize and relax methods for bounding pODEs.
# See https://github.com/PSORLab/DynamicBoundspODEsDiscrete.jl
#############################################################################
# src/DynamicBoundspODEsDiscrete.jl
# Defines main module.
#############################################################################

module DynamicBoundspODEsDiscrete

using McCormick, DocStringExtensions, DynamicBoundsBase,
      Reexport, LinearAlgebra, StaticArrays, ElasticArrays, Polynomials

using ForwardDiff: Chunk, Dual, Partials, construct_seeds, single_seed,
      JacobianConfig, vector_mode_dual_eval, value, vector_mode_jacobian!,
      jacobian!

using DiffEqSensitivity: extract_local_sensitivities, ODEForwardSensitivityProblem
using DiffEqBase: remake, AbstractODEProblem, AbstractContinuousCallback, solve
using Sundials

using DiffResults: JacobianResult, MutableDiffResult

import DynamicBoundsBase: relax!, set!, setall!, get, getall!, getall, relax!,
                          integrate!, supports
import Base: setindex!, getindex, copyto!, literal_pow, copy

import Base.MathConstants.golden

export DiscretizeRelax, AdamsMoulton, BDF, LohnerContractor, HermiteObreschkoff

export StepParams, StepResult, ExistStorage, ContractorStorage, reinitialize!,
       existence_uniqueness!, improvement_condition, single_step!, set_xX!,
       state_contractor_steps, state_contractor_γ, state_contractor_k, excess_error,
       set_Δ!, compute_X0!, set_P!, contains, calc_alpha, mul_split!, μ!, ρ!

const DBB = DynamicBoundsBase

abstract type AbstractStateContractor end

"""
AbstractStateContractorName

The subtypes of `AbstractStateContractorName` are used
to specify the manner of contractor method to be used
by `DiscretizeRelax` in the discretize and relax scheme.
"""
abstract type AbstractStateContractorName end

"""
state_contractor_k(d::AbstractStateContractorName)

Retrieves the order of the existence test to
be used with
"""
function state_contractor_k(d::AbstractStateContractorName)
    error("No method with AbstractStateContractorName $d defined.")
end

"""
state_contractor_γ(d::AbstractStateContractorName)
"""
function state_contractor_γ(d::AbstractStateContractorName)
    error("No method with AbstractStateContractorName $d defined.")
end

"""
state_contractor_steps(d::AbstractStateContractorName)
"""
function state_contractor_steps(d::AbstractStateContractorName)
    error("No method with AbstractStateContractorName $d defined.")
end

"""
state_contractor_integrator(d::AbstractStateContractorName)
"""
function state_contractor_integrator(d::AbstractStateContractorName)
    error("No method with AbstractStateContractorName $d defined.")
end

"""
μ!(out,xⱼ,x̂ⱼ,η)

Used to compute the arguments of Jacobians (`x̂ⱼ + η(xⱼ - x̂ⱼ)`) used by the parametric Mean Value
Theorem. The result is stored to `out`.
"""
function μ!(out::Vector{Interval{Float64}}, xⱼ::Vector{Interval{Float64}}, x̂ⱼ::Vector{Float64}, η::Interval{Float64})
    out .= xⱼ
    return nothing
end
function μ!(out::Vector{MC{N,T}}, xⱼ::Vector{MC{N,T}}, x̂ⱼ::Vector{Float64}, η::Interval{Float64}) where {N, T<:RelaxTag}
    @__dot__ out = x̂ⱼ + η*(xⱼ - x̂ⱼ)
    return nothing
end

"""
ρ!(out,p,p̂ⱼ,η)

Used to compute the arguments of Jacobians (`p̂ⱼ + η(p - p̂ⱼ)`) used by the parametric Mean Value Theorem.
The result is stored to `out`.
"""
function ρ!(out::Vector{Interval{Float64}}, p::Vector{Interval{Float64}}, p̂::Vector{Float64}, η::Interval{Float64})
    out .= p
    return nothing
end
function ρ!(out::Vector{MC{N,T}}, p::Vector{MC{N,T}}, p̂::Vector{Float64}, η::Interval{Float64}) where {N, T<:RelaxTag}
    @__dot__ out = p̂ + η*(p - p̂)
    return nothing
end

include("StaticTaylorSeries/StaticTaylorSeries.jl")
using .StaticTaylorSeries

include("DiscretizeRelax/utilities/fixed_buffer.jl")
include("DiscretizeRelax/utilities/mul_split.jl")
include("DiscretizeRelax/utilities/fast_set_index.jl")
include("DiscretizeRelax/utilities/qr_utilities.jl")
include("DiscretizeRelax/utilities/coeff_calcs.jl")

include("DiscretizeRelax/utilities/taylor_functor.jl")
include("DiscretizeRelax/utilities/jacobian_functor.jl")
include("DiscretizeRelax/utilities/single_step.jl")

include("DiscretizeRelax/method/higher_order_enclosure.jl")
include("DiscretizeRelax/method/lohners_qr.jl")
include("DiscretizeRelax/method/hermite_obreschkoff.jl")
include("DiscretizeRelax/method/wilhelm_2019.jl")
include("DiscretizeRelax/method/pilms.jl")

include("DiscretizeRelax/utilities/discretize_relax.jl")
include("DiscretizeRelax/utilities/relax.jl")
include("DiscretizeRelax/utilities/access_functions.jl")

end # module


text:
#!/usr/bin/env julia

"""
This program generates batch sricpt to run IGV
Loads BAM files from localhost tunnel
Reads Variant positions from VEP output files for each chromosome
"""

fh = ARGS[1] #      File with list of .vep file names
file = open(fh)
vep = readlines(file)

fP = ARGS[2] # GBM-pairs.txt
fil2 = open(fP)
bamList = readlines(fil2)

here = pwd() # presently working directory

for line in bamList[1:5]
    out = []
    target = chomp(line)
    bam = split(target, ',')
    tumor, normal = bam[1], bam[2]
    #rm("$here/PNG/$target", recursive=true)
    mkdir("$here/PNG/$target")  # Creating snapshot directory
    push!(out, "new")
    push!(out, "genome hg19")
    push!(out, "load http://localhost:8000/bam/$tumor.bam,$here/somatic.wig,$here/germline.wig.tdf,http://localhost:8000/bam/$normal.bam")
    push!(out, "maxPanelHeight 1000")
    push!(out, "snapshotDirectory $here/PNG/$target")
    for file in vep
        file = chomp(file)
        ln = open(file)
        lines = readlines(ln)
        for po in lines
            if ismatch(r"^\#", po)
                continue
            else
                pos = split(po, '\t')
                allele = split(pos[3], ':')
                if allele == "0:0" || parse(Int, allele[2] < 3
                    continue
                else
                    position = "goto "*pos[2]
                    push!(out, position)
                    push!(out, "sort position")
                    push!(out, "snapshot")
                end
            end
        end
    end
    push!(out, "exit")
    writedlm("$target.igv.txt", out)
end


text:
using Test, SHA, SparseArrays
using SBML
import Pkg

@testset "SBML test suite" begin
    include("version.jl")
    include("ecoli_flux.jl")
    include("loadmodels.jl")
end


text:
using MPI
using CLIMA
using CLIMA.Mesh.Topologies
using CLIMA.Mesh.Grids
using CLIMA.DGBalanceLawDiscretizations
using Printf
using LinearAlgebra
using Logging
using CLIMA.MPIStateArrays

if !@isdefined integration_testing
    const integration_testing = parse(
        Bool,
        lowercase(get(ENV, "JULIA_CLIMA_INTEGRATION_TESTING", "false")),
    )
end

@inline function auxiliary_state_initialization!(aux, x, y, z)
    @inbounds begin
        r = hypot(x, y, z)
        aux[1] = r
        aux[2] = x / r
        aux[3] = y / r
        aux[4] = z / r
    end
end

using Test
function run(mpicomm, topl, ArrayType, N, FT)
    grid = DiscontinuousSpectralElementGrid(
        topl,
        FloatType = FT,
        DeviceArray = ArrayType,
        polynomialorder = N,
        meshwarp = Topologies.cubedshellwarp,
    )

    spacedisc = DGBalanceLaw(
        grid = grid,
        length_state_vector = 0,
        flux! = (x...) -> (),
        numerical_flux! = (x...) -> (),
        numerical_boundary_flux! = (x...) -> (),
        auxiliary_state_length = 4,
        auxiliary_state_initialization! = (x...) ->
            auxiliary_state_initialization!(x...),
    )

    exact_aux = copy(spacedisc.auxstate)

    DGBalanceLawDiscretizations.grad_auxiliary_state!(spacedisc, 1, (2, 3, 4))

    euclidean_distance(exact_aux, spacedisc.auxstate)
end

let
    CLIMA.init()
    ArrayTypes = (CLIMA.array_type(),)

    mpicomm = MPI.COMM_WORLD

    numelem = (5, 5, 1)

    polynomialorder = 4

    base_Nhorz = 4
    base_Nvert = 2
    Rinner = 1 // 2
    Router = 1

    expected_result = [
        2.0924087890777517e-04
        1.3897932154337201e-05
        8.8256018429045312e-07
        5.5381072850485303e-08
    ]

    lvls = integration_testing ? length(expected_result) : 1

    @testset "$(@__FILE__)" for ArrayType in ArrayTypes
        for FT in (Float64,) #Float32)
            err = zeros(FT, lvls)
            for l in 1:lvls
                @info (ArrayType, FT, "sphere", l)
                Nhorz = 2^(l - 1) * base_Nhorz
                Nvert = 2^(l - 1) * base_Nvert
                Rrange = range(FT(Rinner); length = Nvert + 1, stop = Router)
                topl = StackedCubedSphereTopology(mpicomm, Nhorz, Rrange)
                err[l] = run(mpicomm, topl, ArrayType, polynomialorder, FT)

                @test expected_result[l] ≈ err[l]
            end
            if integration_testing
                @info begin
                    msg = ""
                    for l in 1:(lvls - 1)
                        rate = log2(err[l]) - log2(err[l + 1])
                        msg *= @sprintf("\n  rate for level %d = %e\n", l, rate)
                    end
                    msg
                end
            end
        end
    end
end

nothing

max_stars_repo_name: yairchn/CLIMA

max_issues_repo_name: yairchn/CLIMA

max_forks_repo_name: yairchn/CLIMA


text:
function mark_perimeter!(r)
    side1 = Sud
    side2 = West
    h, w = go_to_corner!(r, side1, side2)
    for s in (Ost, Nord, West, Sud)
        mark_move!(r, s)
    end
    move_back(r, inverse_side(side1), inverse_side(side2), h, w)
end

max_stars_repo_name: dmitrii-a-piskunov/mirea-progs

max_issues_repo_name: dmitrii-a-piskunov/mirea-progs

max_forks_repo_name: dmitrii-a-piskunov/mirea-progs


text:
export read_data, read_data!

function read_data_seisio!(S::SeisData, filestr::String, v::Int64)
  S_in = rseis(filestr)
  L = length(S_in)
  for i in 1:L
    if typeof(S_in[i]) == SeisChannel
      push!(S, S_in[i])
      break
    elseif typeof(S_in[i]) <: GphysChannel
      C = convert(SeisChannel, S_in[i])
      push!(S, C)
      break
    elseif typeof(S_in[i]) == SeisData
      append!(S, S_in[i])
      break
    elseif typeof(S_in[i]) <: GphysData
      S1 = convert(SeisData, S_in[i])
      append!(S, S1)
      break
    elseif typeof(S_in[i]) == SeisEvent
      (v > 0) && @warn(string("Obj ", i, " is type SeisEvent: only reading :data field"))
      S1 = convert(SeisData, getfield(S_in[i], :data))
      append!(S, S1)
      break
    else
      (v > 0) && @warn(string("Obj ", i, " skipped (Type incompatible with read_data)"))
    end
  end
  return nothing
end

@doc """
    read_data

Generic wrapper for reading file data.

    S = read_data(fmt, filestr [, keywords])

Read data in file format `fmt` matching file pattern `filestr` into a new
SeisData object `S`.

    S = read_data(filestr [, keywords])

Read from files matching file pattern `filestr` into a new SeisData object `S`.
Calls `guess(filestr)` to identify the file type based on the first file
matching pattern `filestr`. Much slower than manually specifying file type.

This single-string syntax will not work with SACPZ or SEED RESP files; `guess`
cannot detect them.

    read_data!(S, fmt, filestr [, keywords])

Read data in file format `fmt` matching file pattern `filestr` into an existing
SeisData object `S`.

    read_data!(S, filestr [, keywords])

As above, but calls `guess(filestr)` to set the format based on the first file
matching pattern `filestr`.

| Format                    | String          |
| :---                      | :---            |
| AH-1                      | ah1             |
| AH-2                      | ah2             |
| Bottle                    | bottle          |
| GeoCSV, time-sample pair  | geocsv          |
| GeoCSV, sample list       | geocsv.slist    |
| Lennartz SLIST            | lennartz        |
| Mini-SEED                 | mseed           |
| PASSCAL SEG Y             | passcal         |
| PC-SUDS                   | suds            |
| SAC                       | sac             |
| SEG Y (rev 0 or rev 1)    | segy            |
| SEISIO                    | seisio          |
| SLIST (ASCII sample list) | slist           |
| UW                        | uw              |
| Win32                     | win32           |

|KW      | Used By  | Type    | Default   | Meaning                         |
|:---    |:---      |:---     |:---       |:---                             |
|cf      | win32    | String  | ""        | win32 channel info filestr      |
|full    | sac      | Bool    | false     | read full header into `:misc`?  |
|        | segy     |         |           |                                 |
|        | suds     |         |           |                                 |
|        | uw       |         |           |                                 |
|        | ah       |         |           |                                 |
|nx_add  | mseed    | Int64   | 360000    | min. increase when resizing `:x`|
|        | bottle   |         |           |                                 |
|        | win32    |         |           |                                 |
|nx_new  | mseed    | Int64   | 86400000  | `length(C.x)` for new channels  |
|        | bottle   |         |           |                                 |
|        | win32    |         |           |                                 |
|jst     | win32    | Bool    | true      | are sample times JST (UTC+9)?   |
|swap    | mseed    | Bool    | true      | byte swap?                      |
|v       | all      | Int64   | 0         | verbosity                       |

### Format-Dependent Performance Tips

#### `mseed`, `win32`, `bottle`
Adjust `nx_new` and `nx_add` based on expected data vector lengths. If the
largest data vector has `Nmax` samples, and the smallest has `Nmin`, we
recommend `nx_new`≈`Nmin` and `nx_add`≈`Nmax - Nmin`.

Default values can be changed in SeisIO keywords, e.g.,
```julia
SeisIO.KW.nx_new = 60000
SeisIO.KW.nx_add = 360000
```

The system-wide defaults are `nx_new=86400000` and `nx_add=360000`. Using these
values with very small jobs will greatly decrease performance.

#### SeisIO native format

`read_data("seisio", ...)` largely exists as a convenience wrapper; it reads
only the first SeisIO object from each file that can be converted to a SeisData
structure. For more complicated read operations on SeisIO files, use `rseis`.

## Examples
1. `S = read_data("uw", "99011116541W", full=true)`
    + Read UW-format data file `99011116541W`
    + Store full header information in `:misc`
2. `read_data!(S, "sac", "MSH80*.SAC")`
    + Read SAC-format files matching string pattern `MSH80*.SAC`
    + Read into existing SeisData object `S`
3. `S = read_data("win32", "20140927*.cnt", cf="20140927*ch", nx_new=360000)`
    + Read win32-format data files with names matching pattern `2014092709*.cnt`
    + Use ASCII channel information filenames that match pattern `20140927*ch`
    + Assign new channels an initial size of `nx_new` samples

See also: SeisIO.KW, get_data, guess, rseis
""" read_data!
function read_data!(S::GphysData, fmt::String, filestr::String;
  full    ::Bool    = false,              # full SAC/SEGY hdr
  cf      ::String  = "",                 # win32 channel info file
  jst     ::Bool    = true,               # are sample times JST (UTC+9)?
  nx_add  ::Int64   = KW.nx_add,          # append nx_add to overfull channels
  nx_new  ::Int64   = KW.nx_new,          # new channel samples
  swap    ::Bool    = false,              # do byte swap?
  v       ::Int64   = KW.v                # verbosity level
  )

  one_file = safe_isfile(filestr)

  if fmt == "sac"
    fv = getfield(BUF, :sac_fv)
    iv = getfield(BUF, :sac_iv)
    cv = getfield(BUF, :sac_cv)
    checkbuf_strict!(fv, 70)
    checkbuf_strict!(iv, 40)
    checkbuf_strict!(cv, 192)
    if one_file
      read_sac_file!(S, filestr, fv, iv, cv, full)
    else
      files = ls(filestr)
      for fname in files
        read_sac_file!(S, fname, fv, iv, cv, full)
      end
    end

  elseif fmt == "seisio"
    if one_file
      read_data_seisio!(S, filestr, v)
    else
      files = ls(filestr)
      for fname in files
        read_data_seisio!(S, fname, v)
      end
    end

  elseif ((fmt == "miniseed") || (fmt == "mseed"))
    setfield!(BUF, :swap, swap)
    if one_file
      read_mseed_file!(S, filestr, v, nx_new, nx_add)
    else
      files = ls(filestr)
      for fname in files
        read_mseed_file!(S, fname, v, nx_new, nx_add)
      end
    end

# ============================================================================
# Data formats that aren't SAC, SEISIO, or SEED begin here and are alphabetical
# by first KW

  elseif fmt == "ah1"
    if one_file
      read_ah1!(S, filestr, v=v, full=full)
    else
      files = ls(filestr)
      for fname in files
        read_ah1!(S, fname, v=v, full=full)
      end
    end

  elseif fmt == "ah2"
    if one_file
      read_ah2!(S, filestr, v=v, full=full)
    else
      files = ls(filestr)
      for fname in files
        read_ah2!(S, fname, v=v, full=full)
      end
    end

  elseif fmt == "bottle"
    read_bottle!(S, filestr, v, nx_new, nx_add)

  elseif fmt in ("geocsv", "geocsv.tspair", "geocsv", "geocsv.slist")
    tspair = (fmt == "geocsv" || fmt == "geocsv.tspair")
    if one_file
      read_geocsv_file!(S, filestr, tspair)
    else
      files = ls(filestr)
      for fname in files
        read_geocsv_file!(S, fname, tspair)
      end
    end

  elseif fmt == "lennartz"
    if one_file
      read_slist!(S, filestr, lennartz=true)
    else
      files = ls(filestr)
      for fname in files
        read_slist!(S, fname, lennartz=true)
      end
    end

  elseif fmt == "segy" || fmt == "passcal"
    passcal = fmt == "passcal"
    buf     = getfield(BUF, :buf)
    shorts  = getfield(BUF, :int16_buf)
    ints    = getfield(BUF, :int32_buf)
    checkbuf!(buf, 240)

    if one_file
      append!(S, read_segy_file(filestr, buf, shorts, ints, passcal, swap, full))
    else
      files = ls(filestr)
      for fname in files
        append!(S, read_segy_file(fname, buf, shorts, ints, passcal, swap, full))
      end
    end

  elseif fmt == "slist"
    if one_file
      read_slist!(S, filestr)
    else
      files = ls(filestr)
      for fname in files
        read_slist!(S, fname)
      end
    end

  elseif fmt == "suds"
    if one_file
      append!(S, SUDS.read_suds(filestr, full=full, v=v))
    else
      files = ls(filestr)
      for fname in files
        append!(S, SUDS.read_suds(fname, full=full, v=v))
      end
    end

  elseif fmt == "uw"
    if one_file
      append!(S, UW.uwdf(filestr, v=v, full=full))
    else
      files = ls(filestr)
      for fname in files
        append!(S, UW.uwdf(fname, v=v))
      end
    end

  elseif fmt == "win32" || fmt =="win"
    readwin32!(S, filestr, cf, jst=jst, nx_new=nx_new, nx_add=nx_add, v=v)

  else
    error("Unknown file format!")
  end

  return nothing
end

@doc (@doc read_data!)
function read_data(fmt::String, filestr::String;
  full    ::Bool    = false,              # full SAC/SEGY hdr
  cf      ::String  = "",                 # win32 channel info file
  jst     ::Bool    = true,               # are sample times JST (UTC+9)?
  nx_add  ::Int64   = KW.nx_add,          # append nx_add to overfull channels
  nx_new  ::Int64   = KW.nx_new,          # new channel samples
  swap    ::Bool    = false,              # do byte swap?
  v       ::Int64   = KW.v                # verbosity level
  )

  S = SeisData()
  read_data!(S, fmt, filestr,
    full    = full,
    cf      = cf,
    jst     = jst,
    nx_add  = nx_add,
    nx_new  = nx_new,
    swap    = swap,
    v       = v
    )
  return S
end

function read_data(filestr::String;
  full    ::Bool    = false,              # full SAC/SEGY hdr
  cf      ::String  = "",                 # win32 channel info file
  jst     ::Bool    = true,               # are sample times JST (UTC+9)?
  nx_add  ::Int64   = KW.nx_add,          # append nx_add to overfull channels
  nx_new  ::Int64   = KW.nx_new,          # new channel samples
  v       ::Int64   = KW.v                # verbosity level
  )

  if safe_isfile(filestr)
    g = guess(filestr)
  else
    files = ls(filestr)
    g = guess(files[1])
  end
  S = SeisData()
  read_data!(S, g[1], filestr,
    full    = full,
    cf      = cf,
    jst     = jst,
    nx_add  = nx_add,
    nx_new  = nx_new,
    swap    = g[2],
    v       = v
    )
  return S
end

function read_data!(S::GphysData, filestr::String;
  full    ::Bool    = false,              # full SAC/SEGY hdr
  cf      ::String  = "",                 # win32 channel info file
  jst     ::Bool    = true,               # are sample times JST (UTC+9)?
  nx_add  ::Int64   = KW.nx_add,          # append nx_add to overfull channels
  nx_new  ::Int64   = KW.nx_new,          # new channel samples
  v       ::Int64   = KW.v                # verbosity level
  )

  if safe_isfile(filestr)
    g = guess(filestr)
  else
    files = ls(filestr)
    g = guess(files[1])
  end

  read_data!(S, g[1], filestr,
    full    = full,
    cf      = cf,
    jst     = jst,
    nx_add  = nx_add,
    nx_new  = nx_new,
    swap    = g[2],
    v       = v
    )
  return S
end


text:
#!/usr/bin/env julia

using ArgParse

s = ArgParseSettings()
@add_arg_table s begin
    "-M"
    help = "Print out the list of files this program will create."
    action = :store_true
    "config"
    help = "a YAML configuration file"
    default = "config.yaml"
end

parsed_args = parse_args(ARGS, s)

import YAML
config = YAML.load(open(parsed_args["config"]))

using DiskTracer.model
using DiskTracer.constants

species = config["species"]
transition = config["transition"]
lam0 = lam0s[species*transition]
model = config["model"]
pars = convert_dict(config["parameters"], config["model"])

mol = mols[species*transition]

import PyPlot.plt
using LaTeXStrings

# Make a 2D plot of the source function and Upsilon
# density structure
function plot_S(pars::AbstractParameters)

    nr = 128
    nz = 128
    rs = linspace(0.1, 600 * AU, nr)
    zs = linspace(0.1, 600 * AU, nz)

    rr = rs./AU
    zz = zs./AU

    xx = Array{Float64}(nz, nr)
    yy = Array{Float64}(nz, nr)

    Ss = Array{Float64}(nz, nr)

    for i=1:nz
        xx[i, :] = rr
    end

    for j=1:nr
        yy[:, j] = zz
    end

    for i=1:nz
        for j=1:nr
            Ss[i,j] = DiskTracer.model.S_nu(rs[j], zs[i], mol.nu_0, pars)
        end
    end

    fig = plt[:figure]()
    ax = fig[:add_subplot](111)

    ax[:set_ylabel](L"$z$ [AU]")
    ax[:set_xlabel](L"$r$ [AU]")

    img = ax[:contourf](xx, yy, log10.(Ss)) #
    # img = ax[:contourf](xx, yy, Ss, levels=levels)

    # ax[:set_xscale]("log")

    fig[:subplots_adjust](left=0.15, bottom=0.15, right=0.77)

    cax = fig[:add_axes]([0.82, 0.22, 0.03, 0.65])
    cb = fig[:colorbar](img, cax=cax)

    plt[:savefig]("source_function.png")
end

function plot_Upsilon(pars::AbstractParameters)

    nr = 128
    nz = 128
    rs = linspace(0.1, 600 * AU, nr)
    zs = linspace(0.1, 600 * AU, nz)

    rr = rs./AU
    zz = zs./AU

    xx = Array{Float64}(nz, nr)
    yy = Array{Float64}(nz, nr)

    Us = Array{Float64}(nz, nr)

    for i=1:nz
        xx[i, :] = rr
    end

    for j=1:nr
        yy[:, j] = zz
    end

    for i=1:nz
        for j=1:nr
            Us[i,j] = DiskTracer.model.Upsilon(rs[j], zs[i], mol.nu_0, pars, mol)
        end
    end

    fig = plt[:figure]()
    ax = fig[:add_subplot](111)

    ax[:set_ylabel](L"$z$ [AU]")
    ax[:set_xlabel](L"$r$ [AU]")

    ext = (rr[1], rr[end], zz[1], zz[end])
    img = ax[:imshow](log10.(Us), vmin=-20, origin="lower", extent=ext) #
    # img = ax[:contourf](xx, yy, Ss, levels=levels)

    # ax[:set_xscale]("log")

    fig[:subplots_adjust](left=0.15, bottom=0.15, right=0.77)

    cax = fig[:add_axes]([0.82, 0.22, 0.03, 0.65])
    cb = fig[:colorbar](img, cax=cax)

    plt[:savefig]("Upsilon.png")
end

# plot_S(pars)
plot_Upsilon(pars)


text:
# test/test_all_mirt.jl

using Test

srclist = (
"../data",
"algorithm",
"fbp",
"io",
"plot",
"regularize",
"system",
"utility",
)

"""
test_all_mirt()

run all MIRT tests
"""
function test_all_mirt()

	@testset "all MIRT" begin
		for root in srclist
		#	@show root
			@testset "$root" begin
				include("../src/" * root * "/z-test.jl")
			end
		end
	end

	true
end


text:
" Implements a CoupledConvection operator "
type ADCoupledConvection <: Kernel
    u::Variable
    other_var::Variable
end

import MOOSE.computeQpResidual

@inline function computeQpResidual(kernel::ADCoupledConvection, qp::Int64, i::Int64)
    u = kernel.u

    other_var = kernel.other_var

    return (other_var.grad[qp] ⋅ u.grad[qp]) * u.phi[qp][i]
end

# Specifically leaving out computeQpJacobian() to test Automatic Differentiation
# Note: Any Kernel can be used with AD, even if it defines computeQpJacobian()
# We're just leaving it out here to make _doubly_ sure that there is no way
# computeQpJacobian() is getting called!


text:
# procedure for SDP, *in this order*:
# instantiate SDPSession(nconstraints, nmoments)
# call sdp_obj! a bunch of times in order to set objective
# call sdp_con! a bunch of times (in any order) to set constraints)
# call sdp_solve to return a solution

# solution structure & methods
type SDPSolution
    primalobj :: Float64
    dualobj :: Float64
    primalmatrix :: Array{Float64,2}
    dualmatrix :: Array{Float64,2}
    status :: Symbol # :Normal, :Infeasible, :Unbounded, :Warning, :Error
end

function sdp_sol_infeasible(maximize :: Bool)
    obj = maximize ? -Inf : Inf
    SDPSolution(obj, obj, zeros(0,0), zeros(0,0), :Infeasible)
end

function sdp_sol_unbounded(maximize :: Bool)
    obj = maximize ? Inf : -Inf
    SDPSolution(obj, obj, zeros(0,0), zeros(0,0), :Unbounded)
end

type SDPSession
    maximize :: Bool
    nmoments :: Int64
    initval :: SparseMatrixCSC{Float64,Int64}
    nconstraints :: Int64
    constraints_missing_objective :: Int64
    objective :: Array{Float64,1}
    offset :: Float64
    fname :: AbstractString
    io :: IO
end
SDPSession(maximize, nconstraints, nmoments, offset) = SDPSession(maximize, nmoments, spzeros(nmoments,nmoments), nconstraints, nconstraints, Float64[], offset, sdp_init(nconstraints,nmoments)...)

function sdp_init(nconstraints, nmoments)
    fname,io = mktemp()
    println(io, nconstraints)
    println(io, 1) # number of blocks
    println(io, nmoments)
    fname,io
end

function sdp_obj!(sess :: SDPSession, val :: Float64)
    push!(sess.objective, val)
    @printf(sess.io, "%.18e ", sess.maximize ? -val : val)
    if (sess.constraints_missing_objective -= 1) == 0
        println(sess.io)
    end
end

function sdp_con!(sess :: SDPSession, con :: Int64, i :: Int64, j :: Int64, val :: Float64)
    if con > sess.nconstraints || con < 0
        throw(ArgumentError("invalid constraint number $con ($(sess.nconstraints) expected"))
    end
    if sess.nconstraints == 0
        sess.initval[i,j] = val
        sess.initval[j,i] = val
    end
    @printf(sess.io, "%d 1 %d %d %.18e\n", con, i, j, val)
end

retcode_to_status = [
    :Unbounded, #"Success, but the problem is primal infeasible.",
    :Infeasible, #"Success, but the problem is dual infeasible.",
    :Warning, #"Partial success: full accuracy was not achieved.",
    :Error, #"Failure: maximum iterations reached.",
    :Error, #"Failure: stuck at edge of primal infeasibility.",
    :Error, #"Failure: stuck at edge of dual infeasibility.",
    :Error, #"Failure: lack of progress.",
    :Error, #"Failure: X, Z, or O was singular.",
    :Error, #"Failure: detected NaN on Inf values."
]

function sdp_solve(sess :: SDPSession; call="csdp")
    close(sess.io)

    # First handle totally unconstrained SDPs -- these throw an error in CSDP.
    # The dual matrix is totally constrained to minus `initval` -- check if this is PSD.
    if sess.nconstraints == 0
        eigmin = eigs(-sess.initval; nev=1, which=:SR)[1][1] # least eigenvalue
        if eigmin > -1e-9 # initial value is PSD, feasible
            return SDPSolution(sess.offset,sess.offset,full(-sess.initval),zeros(sess.nmoments,sess.nmoments),0)
        else # infeasible
            println("The SDP is infeasible.")
            return sdp_unbounded(sess.maximize)
        end
    end

    outfname, outio = mktemp()

    primalobj = NaN
    dualobj = NaN
    retcode = 0

    # note: our primal is CSDP's dual, so everything here is flipped that way
    for l in eachline(`$call $(sess.fname) $outfname`)
        print(l)

        if startswith(l, "Failure: return code is ")
            retcode = int(strip(split(l, "is ")[2]))
            if retcode > 0 && retcode <= length(csdp_messages)
                println(csdp_messages[retcode])
            else
                println("Unrecognized return value.")
            end

        elseif startswith(l, "Success: SDP is primal infeasible")
            println("The SDP is infeasible.")
            return sdp_sol_infeasible(sess.maximize)
        elseif startswith(l, "Success: SDP is dual infeasible")
            println("The SDP is unbounded.")
            return sdp_sol_unbounded(sess.maximize)
        end
    end

    dualvector = map(float, split(chomp(readline(outio))))

    primalmatrix = zeros(sess.nmoments, sess.nmoments)
    dualmatrix = zeros(sess.nmoments, sess.nmoments)
    for l in eachline(outio)
        toks = split(chomp(l), " ")
        i = parse(Int,toks[3])
        j = parse(Int,toks[4])
        v = parse(Float64,toks[5])
        if parse(Int,toks[1]) == 1
            primalmatrix[i,j] = v
            primalmatrix[j,i] = v
        else
            dualmatrix[i,j] = v
            dualmatrix[j,i] = v
        end
    end

    # compute objective
    if retcode != 1 && retcode != 2
        primalobj = dot(dualvector, sess.objective) + sess.offset
        dualobj = primalobj # XXX this is no good -- want to know the duality gap.
    end

    SDPSolution(primalobj, dualobj, primalmatrix, dualmatrix, retcode==0 ? :Normal : retcode_to_status[retcode])
end

max_stars_repo_name: willperry/SumOfSquaresOptimization


text:
## Preamble
using Plots

include("../AcousticPropagation.jl")

## n²-Linear Profile
c₀ = 1550
z₀ = 1e3
f = 2e3
src = AcousticPropagation.Source(AcousticPropagation.Position(0, z₀), AcousticPropagation.Signal(f))
ocn = AcousticPropagation.Medium((r, z) -> c₀/sqrt(1 + 2.4z/c₀), 3.5e3)
bty = AcousticPropagation.Boundary(z₀)
θ₀ = -acos(ocn.c(0, z₀)/ocn.c(0, 150))
θ₀s = θ₀.*(0:0.5:1)

ray = AcousticPropagation.Ray.(θ₀, src, ocn, bty)

plot(yaxis = :flip)
plot!(ray.Sol, vars = (1, 2))

## n²-Linear Profile (a)
c₀ = 1550
z₀ = 1e3
f = 2e3
src = AcousticPropagation.Source(AcousticPropagation.Position(0, z₀), AcousticPropagation.Signal(f))
ocn = AcousticPropagation.Medium((r, z) -> c₀/sqrt(1 + 2.4z/c₀), 3.5e3)
bty = AcousticPropagation.Boundary(z₀)
θ₀ = -acos(ocn.c(0, z₀)/ocn.c(0, 150))

δθ₀ = π/40
# ray = AcousticPropagation.Ray(θ₀, src, ocn, bty)
beam = AcousticPropagation.Beam(θ₀, src, ocn, bty)
TL(s, n) = -20log10(abs(beam.b(s, n)))

s = range(0, beam.S, length = 400)
n = 300.0.*range(-1, 1, length = 400)

heatmap(s, n, TL,
	seriescolor = cgrad(:jet, rev = true))


text:
@testset "disjunctive.jl" begin
    @testset "getEST(task::Task)" begin
        trailer = SeaPearl.Trailer()
        task1 = SeaPearl.IntVar(1, 3, "t1", trailer)
        task2 = SeaPearl.IntVar(2, 3, "t2", trailer)
        p1 = 2
        p2 = 3
        tasks = Vector{SeaPearl.IntVar}([task1, task2])
        processing_time = [p1,p2]
        constraint = SeaPearl.Disjunctive(tasks, processing_time, trailer)

        @test constraint.active.value
        @test SeaPearl.getEST(constraint.tasks[1]) == 1
        @test SeaPearl.getEST(constraint.tasks[2]) == 2
    end

    @testset "getLCT(task::Task)" begin
        trailer = SeaPearl.Trailer()
        task1 = SeaPearl.IntVar(1, 3, "t1", trailer)
        task2 = SeaPearl.IntVar(2, 3, "t2", trailer)
        p1 = 2
        p2 = 3
        tasks = Vector{SeaPearl.IntVar}([task1, task2])
        processing_time = [p1,p2]
        constraint = SeaPearl.Disjunctive(tasks, processing_time, trailer)

        @test constraint.active.value
        @test SeaPearl.getLCT(constraint.tasks[1]) == 5
        @test SeaPearl.getLCT(constraint.tasks[2]) == 6
    end

    @testset "getECT(task::Task)" begin
        trailer = SeaPearl.Trailer()
        task1 = SeaPearl.IntVar(1, 3, "t1", trailer)
        task2 = SeaPearl.IntVar(2, 3, "t2", trailer)
        p1 = 2
        p2 = 3
        tasks = Vector{SeaPearl.IntVar}([task1, task2])
        processing_time = [p1,p2]
        constraint = SeaPearl.Disjunctive(tasks, processing_time, trailer)

       @test constraint.active.value
        @test SeaPearl.getECT(constraint.tasks[1]) == 3
        @test SeaPearl.getECT(constraint.tasks[2]) == 5
    end

    @testset "getLST(task::Task)" begin
        trailer = SeaPearl.Trailer()
        task1 = SeaPearl.IntVar(1, 3, "t1", trailer)
        task2 = SeaPearl.IntVar(2, 3, "t2", trailer)
        p1 = 2
        p2 = 3
        tasks = Vector{SeaPearl.IntVar}([task1, task2])
        processing_time = [p1,p2]
        constraint = SeaPearl.Disjunctive(tasks, processing_time, trailer)

        @test constraint.active.value
        @test SeaPearl.getLST(constraint.tasks[1]) == 3
        @test SeaPearl.getLST(constraint.tasks[2]) == 3
    end

    @testset "function Disjunctive(earliestStartingTime::Array{IntVar},
                    processingTime::Array{Int}, trailer)" begin

        trailer = SeaPearl.Trailer()
        task1 = SeaPearl.IntVar(1, 1, "task1", trailer)
        task2 = SeaPearl.IntVar(8, 8, "task2", trailer)
        task3 = SeaPearl.IntVar(2, 3, "task3", trailer)
        task4 = SeaPearl.IntVar(1, 9, "task4", trailer)

        p1 = 2
        p2 = 1
        p3 = 3
        p4 = 3

        tasks = Vector{SeaPearl.IntVar}([task1, task2, task3, task4])
        processing_time = [p1, p2, p3, p4]
        constraint = SeaPearl.Disjunctive(tasks, processing_time, trailer)

        @test size(constraint.tasks)[1] == 4
        @test constraint.active.value
        @test constraint in task1.onDomainChange
        @test constraint in task2.onDomainChange
        @test constraint in task3.onDomainChange
        @test constraint in task4.onDomainChange

     end

    @testset "propagate example from article" begin
        trailer = SeaPearl.Trailer()
        task1 = SeaPearl.IntVar(0, 1, "task1", trailer)
        task2 = SeaPearl.IntVar(1, 5, "task2", trailer)
        p1 = 4
        p2 = 1
        tasks = Vector{SeaPearl.IntVar}([task1, task2])
        processing_time = [p1,p2]
        constraint = SeaPearl.Disjunctive(tasks, processing_time, trailer)

        toPropagate = Set{SeaPearl.Constraint}([constraint])
        prunedDomains = SeaPearl.CPModification()

        res = SeaPearl.propagate!(constraint, toPropagate, prunedDomains)

        @test res
        @test length(task1.domain) == 2
        @test length(task2.domain) == 2
        @test 0 in task1.domain
        @test 1 in task1.domain
        @test 4 in task2.domain
        @test 5 in task2.domain
    end

    @testset "propagate! with 3 tasks" begin
        trailer = SeaPearl.Trailer()
        task1 = SeaPearl.IntVar(1, 3, "task1", trailer)
        task2 = SeaPearl.IntVar(2, 3, "task2", trailer)
        task3 = SeaPearl.IntVar(1, 5, "task3", trailer)

        p1 = 2
        p2 = 3
        p3 = 3
        tasks = Vector{SeaPearl.IntVar}([task1, task2, task3])
        processing_time = [p1, p2, p3]
        constraint = SeaPearl.Disjunctive(tasks, processing_time, trailer)

        toPropagate = Set{SeaPearl.Constraint}([constraint])
        prunedDomains = SeaPearl.CPModification()

        res = SeaPearl.propagate!(constraint, toPropagate, prunedDomains)

        @test res
        @test length(task1.domain) == 3
        @test length(task2.domain) == 2
        @test length(task3.domain) == 1
        @test 1 in task1.domain
        @test 2 in task1.domain
        @test 3 in task1.domain
        @test 2 in task2.domain
        @test 3 in task2.domain
        @test 5 in task3.domain
    end

    @testset "propagate! with 4 tasks" begin
        trailer = SeaPearl.Trailer()
        task1 = SeaPearl.IntVar(1, 1, "task1", trailer)
        task2 = SeaPearl.IntVar(8, 8, "task2", trailer)
        task3 = SeaPearl.IntVar(2, 3, "task3", trailer)
        task4 = SeaPearl.IntVar(1, 9, "task4", trailer)

        p1 = 2
        p2 = 1
        p3 = 3
        p4 = 3

        tasks = Vector{SeaPearl.IntVar}([task1, task2, task3, task4])
        processing_time = [p1, p2, p3, p4]
        constraint = SeaPearl.Disjunctive(tasks, processing_time, trailer)

        toPropagate = Set{SeaPearl.Constraint}([constraint])
        prunedDomains = SeaPearl.CPModification()

        res = SeaPearl.propagate!(constraint, toPropagate, prunedDomains)

        @test res
        @test length(task1.domain) == 1
        @test length(task2.domain) == 1
        @test length(task3.domain) == 1
        @test length(task4.domain) == 1

        @test 1 in task1.domain
        @test 8 in task2.domain
        @test 3 in task3.domain
        @test 9 in task4.domain
        @test !constraint.active.value
    end

    @testset "propagate! failed" begin
        trailer = SeaPearl.Trailer()
        task1 = SeaPearl.IntVar(1, 3, "task1", trailer)
        task2 = SeaPearl.IntVar(2, 3, "task2", trailer)
        task3 = SeaPearl.IntVar(1, 4, "task3", trailer)
        p1 = 2
        p2 = 3
        p3 = 3
        tasks = Vector{SeaPearl.IntVar}([task1, task2, task3])
        processing_time = [p1, p2, p3]
        constraint = SeaPearl.Disjunctive(tasks, processing_time, trailer)

        toPropagate = Set{SeaPearl.Constraint}([constraint])
        prunedDomains = SeaPearl.CPModification()

        res = SeaPearl.propagate!(constraint, toPropagate, prunedDomains)

        @test !res
    end

    @testset "propagate! failed" begin
        trailer = SeaPearl.Trailer()

        task1 = SeaPearl.IntVar(1, 1, "task1", trailer)
        task2 = SeaPearl.IntVar(3, 3, "task2", trailer)
        task3 = SeaPearl.IntVar(5, 5, "task3", trailer)

        p1 = 2
        p2 = 3
        p3 = 3
        tasks = Vector{SeaPearl.IntVar}([task1, task2, task3])
        processing_time = [p1, p2, p3]

        constraint = SeaPearl.Disjunctive(tasks, processing_time, trailer)

        toPropagate = Set{SeaPearl.Constraint}([constraint])
        prunedDomains = SeaPearl.CPModification()

        res = SeaPearl.propagate!(constraint, toPropagate, prunedDomains)
        @test !res
    end

    @testset "propagate! full solving with fixed 3 tasks and Infeasible status" begin
        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)

        task1 = SeaPearl.IntVar(1, 1, "task1", trailer)
        task2 = SeaPearl.IntVar(3, 3, "task2", trailer)
        task3 = SeaPearl.IntVar(5, 5, "task3", trailer)

        p1 = 2
        p2 = 3
        p3 = 3
        SeaPearl.addVariable!(model, task1)
        SeaPearl.addVariable!(model, task2)
        SeaPearl.addVariable!(model, task3)

        tasks = Vector{SeaPearl.IntVar}([task1, task2, task3])
        processing_time = [p1, p2, p3]
        SeaPearl.addConstraint!(model,SeaPearl.Disjunctive(tasks, processing_time, trailer))

        variableSelection = SeaPearl.MinDomainVariableSelection{false}()
        status = SeaPearl.solve!(model; variableHeuristic=variableSelection)

        @test status == :Infeasible
        @test length(model.statistics.solutions) == 1 # one infeasible solution
    end

    @testset "propagate! full solving with 3 tasks and Infeasible status" begin
        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)

        task1 = SeaPearl.IntVar(1, 1, "task1", trailer)
        task2 = SeaPearl.IntVar(2, 3, "task2", trailer)
        task3 = SeaPearl.IntVar(1, 5, "task3", trailer)

        p1 = 2
        p2 = 3
        p3 = 3
        tasks = Vector{SeaPearl.IntVar}([task1, task2, task3])
        processing_time = [p1, p2, p3]
        SeaPearl.addConstraint!(model,SeaPearl.Disjunctive(tasks, processing_time, trailer))
        SeaPearl.addVariable!(model, task1)
        SeaPearl.addVariable!(model, task2)
        SeaPearl.addVariable!(model, task3)

        variableSelection = SeaPearl.MinDomainVariableSelection{false}()
        status = @time SeaPearl.solve!(model; variableHeuristic=variableSelection)

        @test status == :Infeasible
        @test model.statistics.numberOfSolutions == 0

    end

    @testset "propagate! chain" begin
        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)

        task1 = SeaPearl.IntVar(1, 3, "task1", trailer)
        task2 = SeaPearl.IntVar(2, 3, "task2", trailer)
        task3 = SeaPearl.IntVar(1, 7, "task3", trailer)

        p1 = 2
        p2 = 3
        p3 = 3
        tasks = Vector{SeaPearl.IntVar}([task1, task2, task3])
        processing_time = [p1, p2, p3]

        constraint = SeaPearl.Disjunctive(tasks, processing_time, trailer)

        toPropagate = Set{SeaPearl.Constraint}([constraint])
        prunedDomains = SeaPearl.CPModification()
        res = SeaPearl.propagate!(constraint, toPropagate, prunedDomains)

        @test res
        @test length(task1.domain) == 3
        @test length(task2.domain) == 2
        @test length(task3.domain) == 3
        @test 1 in task1.domain
        @test 2 in task1.domain
        @test 3 in task1.domain
        @test 2 in task2.domain
        @test 3 in task2.domain
        @test 5 in task3.domain
        @test 6 in task3.domain
        @test 7 in task3.domain

        SeaPearl.remove!(task1.domain, 3)
        res = SeaPearl.propagate!(constraint, toPropagate, prunedDomains)

        @test res
        @test length(task1.domain) == 2
        @test length(task2.domain) == 1
        @test length(task3.domain) == 2
        @test 1 in task1.domain
        @test 2 in task1.domain
        @test 3 in task2.domain
        @test 6 in task3.domain
        @test 7 in task3.domain

        SeaPearl.remove!(task1.domain, 2)
        res = SeaPearl.propagate!(constraint, toPropagate, prunedDomains)
        @test res
        @test length(task1.domain) == 1
        @test length(task2.domain) == 1
        @test length(task3.domain) == 2
        @test 1 in task1.domain
        @test 3 in task2.domain
        @test 6 in task3.domain
        @test 7 in task3.domain
    end

    @testset "full solving with 3 tasks and optimal status" begin
        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)

        task1 = SeaPearl.IntVar(1, 1, "task1", trailer)
        task2 = SeaPearl.IntVar(2, 3, "task2", trailer)
        task3 = SeaPearl.IntVar(1, 6, "task3", trailer)
        SeaPearl.addVariable!(model, task1)
        SeaPearl.addVariable!(model, task2)
        SeaPearl.addVariable!(model, task3)

        p1 = 2
        p2 = 3
        p3 = 3
        tasks = Vector{SeaPearl.IntVar}([task1, task2, task3])
        processing_time = [p1, p2, p3]
        SeaPearl.addConstraint!(model,SeaPearl.Disjunctive(tasks, processing_time, trailer))

        variableSelection = SeaPearl.MinDomainVariableSelection{false}()
        status = @time SeaPearl.solve!(model; variableHeuristic=variableSelection)

        @test status == :Optimal
        @test model.statistics.numberOfSolutions == 1
        @test 1 in task1.domain
        @test 3 in task2.domain
        @test 6 in task3.domain
        @test 1 == length(task1.domain)
        @test 1 == length(task2.domain)
        @test 1 == length(task3.domain)
    end

    @testset "full solving with 3 tasks, optimal status and multiple solution" begin
        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)

        task1 = SeaPearl.IntVar(1, 1, "task1", trailer)
        task2 = SeaPearl.IntVar(2, 3, "task2", trailer)
        task3 = SeaPearl.IntVar(1, 7, "task3", trailer)
        SeaPearl.addVariable!(model, task1)
        SeaPearl.addVariable!(model, task2)
        SeaPearl.addVariable!(model, task3)

        p1 = 2
        p2 = 3
        p3 = 3
        tasks = Vector{SeaPearl.IntVar}([task1, task2, task3])
        processing_time = [p1, p2, p3]
        SeaPearl.addConstraint!(model,SeaPearl.Disjunctive(tasks, processing_time, trailer))

        variableSelection = SeaPearl.MinDomainVariableSelection{false}()
        status = @time SeaPearl.solve!(model; variableHeuristic=variableSelection)

        @test status == :Optimal
        @test model.statistics.numberOfSolutions == 2
        @test 1 in task1.domain
        @test 3 in task2.domain
        @test 6 in task3.domain
        @test 7 in task3.domain
        @test 1 == length(task1.domain)
        @test 1 == length(task2.domain)
        @test 2 == length(task3.domain)
    end

    @testset "job shop test wiht two tasks" begin
        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)
        nbMachine = 2
        nbTask = 2
        timeLimit = 10
        tasks = Matrix(undef, nbTask, nbMachine)
        objectif = SeaPearl.IntVar(0, timeLimit, "obj", trailer)
        SeaPearl.addVariable!(model, objectif)
        processingTime = [1 4; 4 1]
        endTask = Matrix(undef, nbTask, nbMachine)

        for i in 1:nbTask
            for j in 1:nbMachine
                tasks[i,j] = SeaPearl.IntVar(0, timeLimit, "task_"*string(i)*"_"*string(j), trailer)
                endTask[i,j] = SeaPearl.IntVarViewOffset(tasks[i,j], processingTime[i,j], "end_"*string(i)*"_"*string(j))
                SeaPearl.addVariable!(model, tasks[i,j])
                SeaPearl.addConstraint!(model, SeaPearl.LessOrEqual(endTask[i,j], objectif, trailer))
            end
        end

        for i in 1:nbTask
            SeaPearl.addConstraint!(model, SeaPearl.Disjunctive([tasks[i,j] for j in 1:nbMachine],
                                                            [processingTime[i,j] for j in 1:nbMachine],
                                                            trailer))
        end
        for j in 1:nbMachine
            SeaPearl.addConstraint!(model, SeaPearl.Disjunctive([tasks[i,j] for i in 1:nbTask],
                                                            [processingTime[i,j] for i in 1:nbTask],
                                                            trailer))
        end

        SeaPearl.addObjective!(model, objectif)
        variableSelection = SeaPearl.MinDomainVariableSelection{false}()
        status = @time SeaPearl.solve!(model; variableHeuristic=variableSelection,)

        @test status == :Optimal
        @test (model.statistics.solutions[end-1]["obj"]) == 5
    end
end


text:
# Code for video 3_2:

Pkg.add("DataArrays")
using DataArrays

108 + NA # NA
10.5 * NA # NA

darr = @data([7.1, 3.14, NA, 5.01, 42.8])
# 5-element DataArrays.DataArray{Float64,1}:
#   7.1
#   3.14
#    NA
#   5.01
#  42.8

sum(darr) # NA
mean(darr) # NA
sum(dropna(darr)) # 58.05
mean(dropna(darr)) # 14.5125

repl = 2.5 # 2.5
sum(convert(Array, darr, repl)) # 60.55
mean(convert(Array, darr, repl)) # 12.11

# Pooling or factors (categorical data)
pdarr = @pdata(["B", "O", "AB", "AB", "A", "O"])
# 6-element DataArrays.PooledDataArray{ASCIIString,UInt32,1}:
#  "B"
#  "O"
#  "AB"
#  "AB"
#  "A"
#  "O"
pdarr = pool(["B", "O", "AB", "AB", "A", "O"])
# same result

levels(pdarr)
# 4-element Array{ASCIIString,1}:
#  "A"
#  "AB"
#  "B"
#  "O"

Pkg.add("DataFrames")
using DataFrames

df = DataFrame()
df[:sepal_length] = @data([5.1, 4.9, 5.9])
df[:sepal_width] = @data([3.5, 3.0, 3.0])
df[:petal_length] = @data([1.4, 1.4, 5.1])
df[:petal_width] = @data([0.2, 0.2, 1.8])
df[:species] = @data(["I. setosa", "I. setosa", "I. virginica"])
show(df)
3x5 DataFrames.DataFrame
| Row | sepal_length | sepal_width | petal_length | petal_width |
|-----|--------------|-------------|--------------|-------------|
| 1   | 5.1          | 3.5         | 1.4          | 0.2         |
| 2   | 4.9          | 3.0         | 1.4          | 0.2         |
| 3   | 5.9          | 3.0         | 5.1          | 1.8         |

| Row | species        |
|-----|----------------|
| 1   | "I. setosa"    |
| 2   | "I. setosa"    |
| 3   | "I. virginica" |

names = DataFrame(ID = [1, 2, 3], Name = ["Alex Smith", "Melissa Fields", "Brad Pitt"])
3x2 DataFrames.DataFrame
| Row | ID | Name             |
|-----|----|------------------|
| 1   | 1  | "Alex Smith"     |
| 2   | 2  | "Melissa Fields" |
| 3   | 3  | "Brad Pitt"      |

names[:Name]
3-element DataArrays.DataArray{ASCIIString,1}:
 "Alex Smith"
 "Melissa Fields"
 "Brad Pitt"

blood_types = DataFrame(ID = [2, 1, 3], blood_type = ["O", "AB", "B"])
3x2 DataFrames.DataFrame
| Row | ID | blood_type |
|-----|----|------------|
| 1   | 2  | "O"        |
| 2   | 1  | "AB"       |
| 3   | 3  | "B"        |

full = join(names, blood_types, on = :ID)
3x3 DataFrames.DataFrame
| Row | ID | Name             | blood_type |
|-----|----|------------------|------------|
| 1   | 1  | "Alex Smith"     | "AB"       |
| 2   | 2  | "Melissa Fields" | "O"        |
| 3   | 3  | "Brad Pitt"      | "B"        |

rename!(blood_types, :blood_type, :BloodType)
3x2 DataFrames.DataFrame
| Row | ID | BloodType |
|-----|----|-----------|
| 1   | 2  | "O"       |
| 2   | 1  | "AB"      |
| 3   | 3  | "B"       |

fname = "iris.csv"
data = readtable(fname, separator = ',')
150x5 DataFrames.DataFrame
| Row | sepal_length | sepal_width | petal_length | petal_width |
|-----|--------------|-------------|--------------|-------------|
| 1   | 5.1          | 3.5         | 1.4          | 0.2         |
| 2   | 4.9          | 3.0         | 1.4          | 0.2         |
| 3   | 4.7          | 3.2         | 1.3          | 0.2         |
| 4   | 4.6          | 3.1         | 1.5          | 0.2         |
| 5   | 5.0          | 3.6         | 1.4          | 0.2         |
...
| Row | species        |
|-----|----------------|
| 1   | "I. setosa"    |
| 2   | "I. setosa"    |
| 3   | "I. setosa"    |
| 4   | "I. setosa"    |
| 5   | "I. setosa"    |
typeof(data) # DataFrames.DataFrame
size(data) # (150, 5)

show(data[2])
# [3.5,3.0,3.2,3.1,3.6,3.9,3.4,3.4,2.9,3.1,3.7,3.4,3.0,3.0,4.0,4.4,3.9,3.5,3.8,3.8,
# 3.4,3.7,3.6,3.3,3.4,3.0,3.4,3.5,3.4,3.2,3.1,3.4,4.1,4.2,3.1,3.2,3.5,3.1,3.0,3.4,3
# .5,2.3,3.2,3.5,3.8,3.0,3.8,3.2,3.7,3.3,3.2,3.2,3.1,2.3,2.8,2.8,3.3,2.4,2.9,2.7,2.
# 0,3.0,2.2,2.9,2.9,3.1,3.0,2.7,2.2,2.5,3.2,2.8,2.5,2.8,2.9,3.0,2.8,3.0,2.9,2.6,2.4
# ,2.4,2.7,2.7,3.0,3.4,3.1,2.3,3.0,2.5,2.6,3.0,2.6,2.3,2.7,3.0,2.9,2.9,2.5,2.8,3.3,
# 2.7,3.0,2.9,3.0,3.0,2.5,2.9,2.5,3.6,3.2,2.7,3.0,2.5,2.8,3.2,3.0,3.8,2.6,2.2,3.2,2
# .8,2.8,2.7,3.3,3.2,2.8,3.0,2.8,3.0,2.8,3.8,2.8,2.8,2.6,3.0,3.4,3.1,3.0,3.1,3.1,3.
# 1,2.7,3.2,3.3,3.0,2.5,3.0,3.4,3.0]
show(data[:sepal_width])
# [3.5,3.0,3.2,3.1,3.6,3.9,3.4,3.4,2.9,3.1,3.7,3.4,3.0,3.0,4.0,4.4,3.9,3.5,3.8,3.8,
# 3.4,3.7,3.6,3.3,3.4,3.0,3.4,3.5,3.4,3.2,3.1,3.4,4.1,4.2,3.1,3.2,3.5,3.1,3.0,3.4,3
# .5,2.3,3.2,3.5,3.8,3.0,3.8,3.2,3.7,3.3,3.2,3.2,3.1,2.3,2.8,2.8,3.3,2.4,2.9,2.7,2.
# 0,3.0,2.2,2.9,2.9,3.1,3.0,2.7,2.2,2.5,3.2,2.8,2.5,2.8,2.9,3.0,2.8,3.0,2.9,2.6,2.4
# ,2.4,2.7,2.7,3.0,3.4,3.1,2.3,3.0,2.5,2.6,3.0,2.6,2.3,2.7,3.0,2.9,2.9,2.5,2.8,3.3,
# 2.7,3.0,2.9,3.0,3.0,2.5,2.9,2.5,3.6,3.2,2.7,3.0,2.5,2.8,3.2,3.0,3.8,2.6,2.2,3.2,2
# .8,2.8,2.7,3.3,3.2,2.8,3.0,2.8,3.0,2.8,3.8,2.8,2.8,2.6,3.0,3.4,3.1,3.0,3.1,3.1,3.
# 1,2.7,3.2,3.3,3.0,2.5,3.0,3.4,3.0]

data[2:3, :]
# 2x5 DataFrames.DataFrame
# | Row | sepal_length | sepal_width | petal_length | petal_width | species     |
# |-----|--------------|-------------|--------------|-------------|-------------|
# | 1   | 4.9          | 3.0         | 1.4          | 0.2         | "I. setosa" |
# | 2   | 4.7          | 3.2         | 1.3          | 0.2         | "I. setosa" |

julia> data[145:150, [:petal_length, :species]]
6x2 DataFrames.DataFrame
# | Row | petal_length | species        |
# |-----|--------------|----------------|
# | 1   | 5.7          | "I. virginica" |
# | 2   | 5.2          | "I. virginica" |
# | 3   | 5.0          | "I. virginica" |
# | 4   | 5.2          | "I. virginica" |
# | 5   | 5.4          | "I. virginica" |
# | 6   | 5.1          | "I. virginica" |

julia> names(data)
5-element Array{Symbol,1}:
 :sepal_length
 :sepal_width
 :petal_length
 :petal_width
 :species

julia> head(data)
6x5 DataFrames.DataFrame
| Row | sepal_length | sepal_width | petal_length | petal_width | species     |
|-----|--------------|-------------|--------------|-------------|-------------|
| 1   | 5.1          | 3.5         | 1.4          | 0.2         | "I. setosa" |
| 2   | 4.9          | 3.0         | 1.4          | 0.2         | "I. setosa" |
| 3   | 4.7          | 3.2         | 1.3          | 0.2         | "I. setosa" |
| 4   | 4.6          | 3.1         | 1.5          | 0.2         | "I. setosa" |
| 5   | 5.0          | 3.6         | 1.4          | 0.2         | "I. setosa" |
| 6   | 5.4          | 3.9         | 1.7          | 0.4         | "I. setosa" |

julia> tail(data)
6x5 DataFrames.DataFrame
| Row | sepal_length | sepal_width | petal_length | petal_width |
|-----|--------------|-------------|--------------|-------------|
| 1   | 6.7          | 3.3         | 5.7          | 2.5         |
| 2   | 6.7          | 3.0         | 5.2          | 2.3         |
| 3   | 6.3          | 2.5         | 5.0          | 1.9         |
| 4   | 6.5          | 3.0         | 5.2          | 2.0         |
| 5   | 6.2          | 3.4         | 5.4          | 2.3         |
| 6   | 5.9          | 3.0         | 5.1          | 1.8         |

| Row | species        |
|-----|----------------|
| 1   | "I. virginica" |
| 2   | "I. virginica" |
| 3   | "I. virginica" |
| 4   | "I. virginica" |
| 5   | "I. virginica" |
| 6   | "I. virginica" |

eltypes(data)
5-element Array{Type{T},1}:
 Float64
 Float64
 Float64
 Float64
 UTF8String

fname = "iris.dat"
writetable(fname, data, separator = '/', header = false)

whos()
    data   9660 bytes  150x5 DataFrames.DataFrame : 150x5 …
pool!(data, [:species])
whos()
    data   5744 bytes  150x5 DataFrames.DataFrame : 150x5 …

max_stars_repo_path: Section 3/3_2.jl

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

max_issues_repo_path: Section 3/3_2.jl

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

max_forks_repo_path: Section 3/3_2.jl

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


text:
function do_sendrecv()

    comm = MPI.COMM_WORLD

    MPI.Barrier(comm)

    rank = MPI.Comm_rank(comm)
    size = MPI.Comm_size(comm)

    dst = mod(rank+1, size)
    src = mod(rank-1, size)

    N = 4

    send_mesg = Array(Float64, N)
    recv_mesg = Array(Float64, N)

    fill!(send_mesg, Float64(rank))

    rreq = MPI.Irecv!(recv_mesg, src,  src+32, comm)

    println("$rank: Sending   $rank -> $dst = $send_mesg")
    sreq = MPI.Isend(send_mesg, dst, rank+32, comm)

    stats = MPI.Waitall!([rreq, sreq])

    println("$rank: Receiving $src -> $rank = $recv_mesg")

    MPI.Barrier(comm)
end


text:
module FASTQ

import Automa
import Automa.RegExp: @re_str
import Automa.Stream: @mark, @markpos, @relpos, @abspos
import BioSymbols
import BioSequences
import BioGenerics: BioGenerics, isfilled
import BioGenerics.Automa: State
import TranscodingStreams: TranscodingStreams, TranscodingStream
import ..FASTX: identifier, hasidentifier,
    description, hasdescription, header,
    sequence, hassequence

include("quality.jl")
include("record.jl")
include("readrecord.jl")
include("reader.jl")
include("writer.jl")
include("fastqread.jl")

end

max_forks_repo_name: BioJulia/FASTX


text:
using BinaryBuilder

name = "HelloWorldC"
version = v"1.1.0"

# No sources, we're just building the testsuite
sources = [
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
mkdir -p ${prefix}/bin
cc -o ${prefix}/bin/hello_world${exeext} -g -O2 /usr/share/testsuite/c/hello_world/hello_world.c

# Also build with cmake
mkdir build && cd build
cmake .. -DCMAKE_INSTALL_PREFIX=${prefix} -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DCMAKE_BUILD_TYPE=Release
make -j${nproc}

install_license /usr/share/licenses/MIT
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

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

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

# Build the tarballs.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: mortenpi/Yggdrasil

max_issues_repo_name: mortenpi/Yggdrasil

max_forks_repo_name: mortenpi/Yggdrasil


text: include("imagetile.jl")


text:
function sorteigen(evals::Vector{T}, evecs::Matrix{T}) where {T<:Real}
    p = sortperm(abs.(evals), lt=isless, rev=true)
    return evals[p], evecs[:, p]
end


text:
#!/bin/bash
# -*- mode: julia -*-
#=
JULIA="${JULIA:-julia --color=yes --startup-file=no}"
export JULIA_PROJECT="$(dirname ${BASH_SOURCE[0]})"

set -ex
${JULIA} -e 'using Pkg; Pkg.instantiate();
             Pkg.develop(PackageSpec(path=pwd()))'

export JULIA_LOAD_PATH="@"
exec ${JULIA} "${BASH_SOURCE[0]}" "$@"
=#

include("utils.jl")
using Literate
transducers_literate()
transducers_makedocs()
deploydocs(;
    repo="github.com/tkf/Transducers.jl",
)


text:
#EasyPlotGrace AbstractBuilder interface implementation
#-------------------------------------------------------------------------------
#=NOTE:
 - TODO: Not using defaults.renderdpi or DEFAULT_RENDERDPI
=#

#==Types
===============================================================================#
mutable struct Builder{T} <: EasyPlot.AbstractBuilder
	dpi::Int
	postproc::Optional{Function}
	args::Tuple
	kwargs::Base.Iterators.Pairs
end
const GUIBuilder = Builder{:gui}
const HeadlessBuilder = Builder{:headless}
EasyPlot.AbstractBuilder(::DS{:Grace}) = GUIBuilder #Register builder
EasyPlot.AbstractBuilder(::DS{:Grace_headless}) = HeadlessBuilder #Register builder

#==Helper functions
===============================================================================#
isheadless(::Builder) = false
isheadless(::HeadlessBuilder) = true

#==Constructor interface
===============================================================================#
#Use same builder irrespective of target application:
function EasyPlot.getbuilder(::Target{T}, btype::Type{Builder{BT}}, args...;
		postproc=nothing, dpi=GracePlot.DEFAULT_DPI, kwargs...) where {T, BT}
	return btype(dpi, postproc, args, kwargs)
end

#But can't have a gui in headless mode:
function EasyPlot.getbuilder(::Target{:gui}, btype::Type{HeadlessBuilder}, args...;
		postproc=nothing, dpi=GracePlot.DEFAULT_DPI, kwargs...)
	error("Cannot create a GUI-enabled builder in headless mode")
end

#==I/O interface
===============================================================================#
function EasyPlot._show(io::IO, mime::MIME, opt::EasyPlot.ShowOptions, plot::GracePlot.Plot)
	if !EasyPlot.isauto(opt.dim)
		@warn("Cannot currently specify plot size. Using plot defaults.")
		#NOTE: Grace supports rescaling through DPI only.
	end
	GracePlot.show(io, mime, plot) #Uses dimension values in plot
end

EasyPlot._showable(mime::MIME, b::Builder) =
	hasmethod(show, (IO, typeof(mime), GracePlot.Plot))

#==Display interface
===============================================================================#
#Can only redraw (always visible if !headless):
EasyPlot.displaygui(plot::GracePlot.Plot) = (redraw(plot); return plot)

#==Plot building functions
===============================================================================#
#Render a particular EasyPlot.Plot:
function _build(g::GracePlot.GraphRef, eplot::EasyPlot.Plot,
	theme::EasyPlot.Theme, colormgr::ColorMgr)
	set(g, subtitle = eplot.title)
	fold = isa(eplot.xaxis, EasyPlot.FoldedAxis) ? eplot.xaxis : nothing

	builder = WfrmBuilder(g, theme, colormgr, fold)

	for (i, wfrm) in enumerate(eplot.wfrmlist)
		EasyPlot.addwfrm(builder, wfrm, i)
	end

	autofit(g)

	#x-axis properties:
	xmin = eplot.xext.min; xmax = eplot.xext.max
	xscale = Symbol(eplot.xaxis)

	ymin, ymax = (NaN, NaN)
	yscale = :lin
	ylabel = ""
	if length(eplot.ystriplist) > 0
		strip = eplot.ystriplist[1]
		ymin = strip.ext.min; ymax = strip.ext.max
		yscale = strip.scale
		ylabel = strip.axislabel
	end

	#Translate NaN -> nothing (supported by GracePlot).
	isnan(xmin) && (xmin=nothing); isnan(xmax) && (xmax=nothing)
	isnan(ymin) && (ymin=nothing); isnan(ymax) && (ymax=nothing)

	#Apply axes & scale settings:
	set(g, GracePlot.paxes(
		xscale = xscale, yscale = yscale,
		xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax,
	))
	set(g, legend(display=eplot.legend))

	#Apply x/y labels
	if eplot.xlabel != nothing
		set(g, xlabel=eplot.xlabel)
	end
	if ylabel != nothing
		set(g, ylabel=ylabel)
	end

	return g
end

function _build(gplot::GracePlot.Plot, ecoll::EasyPlot.PlotCollection)
	ecoll = EasyPlot.condxfrm_multistrip(ecoll, "EasyPlotGrace") #Emulate multi-strip plots
	ncols = ecoll.ncolumns
	nrows = div(length(ecoll.plotlist)-1, ncols)+1
	colormgr = ColorMgr(gplot)

	#Arrange basically allocates all subplots (GracePlot.Plot):
	arrange(gplot, (nrows, ncols))
	graphidx = 0

	title = ecoll.title
	for p in ecoll.plotlist
		g = graph(gplot, graphidx)
		_build(g, p, ecoll.theme, colormgr)
		graphidx += 1
	end

	if length(ecoll.plotlist) > 1
		#Grace has a centered title & subtitle for each subplot.
		#Thus, if we want a centered title with a multi-plot output,
		#we have to do it manually:
		w = get(gplot, :wview); h = get(gplot, :hview)
		vgap = 0.15 #Pick a reasonable value (cannot querry state)
		title = GracePlot.text(ecoll.title, loctype=:view, size=1.5, loc=(w/2,h-vgap/2.5), just=:centercenter)
		GracePlot.addannotation(gplot, title)
	else
		set(graph(gplot, 0), title = ecoll.title)
	end

	redraw(gplot)
	return gplot
end

function EasyPlot.build(b::Builder, ecoll::EasyPlot.PlotCollection)
	plot = GracePlot.new(b.args...; guimode=!isheadless(b), dpi=b.dpi, b.kwargs...)
	return _build(plot, ecoll)
end

#Last line


text:
using DensestCommonSubgraph
using Base.Test
using FactCheck

function test_binary_to_gray(n::Integer=5)
  prev_code = ""
  for i=0:(2^n-1)
    binary_code = bin(i, n)
    gray_code = binary_to_gray(binary_code)
    orig_code = gray_to_binary(gray_code)
    @fact orig_code => binary_code "test_binary_to_gray(): incorrect conversion of $i"
    if i >= 1
      xx = filter((j) -> gray_code[j] != prev_code[j], [1:n])
      @fact length(xx) => 1 "test_binary_to_gray(): incorrect gray coding $(i-1), $i"
    end

    # println("i: ", i, "\tbinary: ", binary_code, "\tgray_code: ", gray_code, "\toriginal", orig_code)
    prev_code = gray_code
  end
end