text:
# MIT License
#
# Copyright (c) 2020 Joshua E Gentry

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

# =================================================================================================
"""
    tess_image(
            inst::TessInst,
            pix::Pix
        )::Nothing

Set the image to be OCRed.

__Arguments:__

| T | Name | Default | Description
|:--| :--- | :------ | :----------
| R | inst |         | The instance to load the image into.
| R | pix  |         | The image to load into Tesseract to perform OCR.

__Details:__

Only 1 image can be set, calling this multiple times will replace the last image.

__Example:__

```jldoctest
julia> using Tesseract

julia> download_languages()
true

julia> instance = TessInst()
Allocated Tesseract instance.

julia> pix = sample_pix()
Image (500, 600) at 32ppi

julia> tess_image(instance, pix)
```

See also: [`tess_resolution`](@ref), [`pix_read`](@ref)
"""
function tess_image(
            inst::TessInst,
            pix::Pix
        )::Nothing

    if is_valid(inst) == true && is_valid(pix) == true
        ccall(
            (:TessBaseAPISetImage2, TESSERACT),
            Cvoid,
            (Ptr{Cvoid},Ptr{Cvoid}),
            inst,
            pix.ptr
        )
    elseif is_valid(inst) == false
        @error "Instance has been freed."
    else
        @error "Pix object has been freed."
    end
    nothing
end

# =================================================================================================
"""
    tess_resolution(
        inst::TessInst, # The instance to configure.
        ppi             # The PPI of the source image.
    )::Nothing

Set the resolution of the image in ppi.

__Arguments:__

| T | Name | Default | Description
|:--| :--- | :------ | :----------
| R | inst |         | The instance of load the image into.
| R | ppi  |         | The PPI (pixels per inch) of the source image.

__Example:__

```jldoctest
julia> using Tesseract

julia> download_languages()
true

julia> instance = TessInst()
Allocated Tesseract instance.

julia> pix = sample_pix()
Image (500, 600) at 32ppi

julia> tess_image(instance, pix)

julia> tess_resolution(instance, 72)
```

See also: [`tess_image`](@ref)
"""
function tess_resolution(
            inst::TessInst,
            ppi::Integer
        )::Nothing
    # ---------------------------------------------------------------------------------------------
    # Make sure the object is still valid.
    if is_valid(inst) == false
        @error "Instance has been freed."
        return nothing
    end

    # ---------------------------------------------------------------------------------------------
    # Make the C call.
    ccall(
        (:TessBaseAPISetSourceResolution, TESSERACT),
        Cvoid,
        (Ptr{Cvoid}, Cint),
        inst,
        ppi
    )
    nothing
end


text:
"""
	Obstacle

Abstract type representing obstacle.
"""
abstract type Obstacle end

"""
	CircularObstacle <: Obstacle

Circular obstacle with center c and radius r
"""
struct CircularObstacle <: Obstacle
    c::Vector{Float64}
    r::Float64
end

"""
	ControlBarrierFunction(O::CircularObstacle, α)

Construct a CBF for a ciruclar obstacle.
"""
function ControlBarrierFunction(O::CircularObstacle, α)
    h(x) = (x[1] - O.c[1])^2 + (x[2] - O.c[2])^2 - O.r^2
    return ControlBarrierFunction(h, α)
end

ext: jl


text:
"""
This module contains functiosn for generating datasets suitable for testing density estimation
"""
module GenerateDatasets

using Distributions
using StatsBase
using StaticArrays
import DensityEstimationML: sample_from_cdf

export approximate_support, original_sample

struct RectangularInterval{T}
    lb::T
    ub::T
end

Base.minimum(iv::RectangularInterval) = iv.lb
Base.maximum(iv::RectangularInterval) = iv.ub

"""
    approximate_support(d)

Returns the effective support for a distribution,
This can be overloaded to provide a restricted support for distribitions that have infinitely large supports,
but can be approximated by bounding them.
Falls back to `support(d)` + 10% to either side
"""
function approximate_support(d)
    sp = support(d)
    a=minimum(sp)
    b=maximum(sp)
    len = b-a
    typeof(sp)(a-0.1*len, b+0.1*len)
end

"""
    original_sample(d)
Returns a sample of a distribution of the size that was used in the corresponding paper.
Falls back to a default of 5000 elements
"""
original_sample(d) = rand(d, 5000)

"""
The distribution used by
Magdon-Ismail and Atiya, in their 1998 NIPS paper "Neural Networks for Denity Estimation."
This is a 1D GMM with two peaks, without any conditioning on any properties.
for their SIC and SLC method they tested with  n=100,
for their SIC they additionally tested with n=200

This dataset is not conditional.
"""
struct MagdonIsmailAndAtiya <: UnivariateDistribution{Continuous}
end

function backing_distribition(::MagdonIsmailAndAtiya)
    comp1 = Normal(-30, 3)
    comp2 = Normal(9, 9)
    MixtureModel([comp1, comp2], [0.3, 0.7])
end

approximate_support(::MagdonIsmailAndAtiya) = RealInterval(-50.0, 50.0)
original_sample(d::MagdonIsmailAndAtiya) = rand(backing_distribition(d), 200)

"""
The distribution used by example 1 of
Likas A, 2001, "Probability density estimation using artificial neural networks"

They only estimate it between -12 and 12

This dataset is not conditional.
"""
struct Likas1 <: UnivariateDistribution{Continuous}
end

function backing_distribition(::Likas1)
    MixtureModel([Normal(-7,0.5), Uniform(-3,-1), Uniform(1,3), Normal(7,0.5)])
end
approximate_support(::Likas1) = RealInterval(-12.0,12.0)
original_sample(d::Likas1) = rand(backing_distribition(d), 5000)

for wrapper in [:MagdonIsmailAndAtiya, :Likas1]
    for (mod,funs) in [(:Distributions, [:sampler, :mode, :modes]),
                       (:Base, [:minimum, :maximum, :mean, :var]),
                       (:StatsBase, [:skewness])
                      ]
        for ff in funs
            @eval $mod.$ff(d::$wrapper)=$mod.$ff(backing_distribition(d))
        end
    end
    for (mod,funs) in [(:Distributions, [:pdf, :logpdf, :cdf, :insupport]),
                       (:Base, [:quantile]),
                       (:StatsBase, [:entropy])
                      ]
        for ff in funs
            @eval $mod.$ff(d::$wrapper, x::Real)=$mod.$ff(backing_distribition(d), x)
        end
    end

    for (mod,funs) in [(:Distributions, [:mgf, :cf])]
        for ff in funs
            @eval $mod.$ff(d::$wrapper, x)=$mod.$ff(backing_distribition(d), x)
        end
    end

    @eval StatsBase.kurtosis(d::$wrapper, x::Bool)=StatsBase.kurtosis(backing_distribition(d),x)
end

"""
Generates a dataset with the same parameters as used in example 2 of tests by
Likas A, 2001, "Probability density estimation using artificial neural networks"
and also in
Modha, D. S. & Fainman, 1994 Y. A learning law for density estimation

The PDF is given by:

\[f(x)=\begin{cases}
    \frac{2-\frac{x}{2}}{6.5523}    & 0\le x<2\\
    \frac{2-(x-3)^{2}}{6.5523}      & 0\le x<3+\sqrt{2}\\
    0                   & otherwise
\end{cases}
\]

Thus the CDF is given by
\[
F(x)=\dfrac{1}{6.5523}\begin{cases}
    2x-\frac{x^{2}}{4} & 0\le x<2\\
    \frac{-x^{3}}{3}+3x^{2}-7x+\frac{23}{3} & 2\le x<3+\sqrt{2}\\
    0 & x<0
    1 & x>3+\sqrt{2}
\end{cases}
\]

Both works use the same 5000 samples
This dataset is not conditional.
"""
struct Likas2<:UnivariateDistribution{Continuous}
end

function Distributions.pdf(::Likas2, x::Real)
    1/6.5523 * if (x ≤ 0)
    0.0
    elseif (0<x ≤ 2)
    2.0-x/2
    elseif (2<x≤3+√2)
    2.0-(x-3.0)^2
    else
    @assert(3+√2 < x)
    0.0
    end
end

function Distributions.cdf(::Likas2, x::Real)
    1/6.5523 * if (x ≤ 0)
    0.0
    elseif (0<x ≤ 2)
    2x-x^2/4
    elseif (2<x≤3+√2)
    -x^3/3 + 3x^2 -7x + 23/3
    else
    @assert(3+√2 < x)
    6.5523
    end
end

function Base.rand(::Likas2)
    likas_2_cdf(x) = cdf(Likas2(), x[1])
    likas_2_pdf(x) = pdf(Likas2(), x[1])
    sample_from_cdf(likas_2_cdf, likas_2_pdf, [2.5])[1]
end

original_sample(d::Likas2) = rand(d, 5000)

Base.minimum(::Likas2) = -1.0
Base.maximum(::Likas2) = 6.0

const ModhaAndFainman=Likas2

"""
Generates a dataset with the same parameters as used in example 3 of tests by
Likas A, 2001, "Probability density estimation using artificial neural networks"

It is a 2D uniform retangle between 0 and 0.2 on both axis

Likas used 5000 samples.

This dataset is not conditional.
"""
struct Likas3 <: ContinuousMultivariateDistribution
end

Distributions.support(::Likas3)=RectangularInterval(@SVector([0., 0.]), @SVector([0.2, 0.2]))
approximate_support(::Likas3)=RectangularInterval(@SVector([-0.1, -0.1]), @SVector([0.3, 0.3]))
Base.rand(::Likas3) = 0.2 .* rand(2)
Base.length(::Likas3) = 2
original_sample(::Likas3, n=5000) = 0.2 .* rand((2, n))

function Distributions._pdf(d::Likas3, X::AbstractVector)
    length(X)==2 || ArgumentError("Likas3 only supports 2D inputs")
    if all(0.<X.<0.2)
        1/(0.2^2)
    else
        0.0
    end
end

Distributions._logpdf(d::Likas3, X::AbstractVector) = log(Distributions._pdf(d,X))

end #module


text:
using BivMatFun
using Test

@testset "BivMatFun.jl" begin
    if isempty(ARGS) || "all" in ARGS
        all_tests = true
    else
        all_tests = false
    end

    if all_tests || "experiment1" in ARGS
        @test include("experiment1.jl");
    end

    if all_tests || "experiment2" in ARGS
        @test include("experiment2.jl");
    end

    if all_tests || "experiment3" in ARGS
        @test include("experiment3.jl");
    end
end

max_stars_repo_name: numpi/BivMatFun

max_issues_repo_name: numpi/BivMatFun

max_forks_repo_name: numpi/BivMatFun


text:
"""
    @qe [expression]

Equivalent to:

    quote
        \$(MacroTools.striplines(esc(expression)))
    end
end
"""
macro qe(ex)
    Expr(:quote, MacroTools.striplines(esc(ex)))
end

macro test_inferred(ex)
    @qe begin
        $Test.@test (($Test.@inferred $ex); true)
    end
end

macro test_noalloc(ex)
    @qe begin
        local nbytes = $memory($BenchmarkTools.@benchmark $ex samples = 1 evals = 1)
        $iszero(nbytes) ? true : $error("Allocated $(BenchmarkTools.prettymemory(nbytes))")
    end
end

macro includetests(testname)
    if __source__.file !== nothing
        _dirname = dirname(String(__source__.file::Symbol))
        root = isempty(_dirname) ? pwd() : abspath(_dirname)
    else
        root = pwd()
    end
    ex = quote
        if $isempty(ARGS)
            local testfiles = $filter(f -> $match(r"^test_.*\.jl$", $basename(f)) !== nothing, $flattendir($root, dirs=false, join=false))
        else
            local testfiles = $map(f -> $endswith(f, ".jl") ? f : "$f.jl", ARGS)
        end
        @includetests $testname testfiles
    end
    esc(ex)
end

macro includetests(testname, testfiles)
    ex = quote
        local testname = $testname
        $Test.@testset "$testname" begin
            for file in $testfiles
                $Test.@testset "$file" begin
                    $Base.include($__module__, file)
                end
            end
        end
    end
    esc(ex)
end


text:
function gradient!(G::AbstractMatrix, P::AbstractMatOrFac, Q::AbstractMatOrFac,
                   Q²::AbstractMatOrFac, Y::AbstractMatrix, Z::Real, β::Real = 0)
    @. G *= β
    attractive_force!(G, P, Q, Y) # for t-SNE
    repulsive_force!(G, P, Q², Y, Z)
    return G
end

# attractive force can be calculated quickly just based on sparse approximation to P
# Z is normalization constant of Q
# function attractive_force!(G::AbstractMatrix, P::SparseMatrixCSC, Q::AbstractMatOrFac, Y::AbstractMatrix)
#     y = [c for c in eachcol(Y)]
#     I, J, _ = findnz(P)
#     for (i, j) in zip(I, J) # IDEA: parallelize over columns?
#          @. G[:, i] += 4 * P[i, j] * Q[i, j] * (y[i] - y[j])
#     end
#     return G
# end

# generic fallback
function attractive_force!(G::AbstractMatrix, P::AbstractMatOrFac, Q::AbstractMatOrFac, Y::AbstractMatrix)
    PQ = hadamard_product(P, Q) # fast through sparsity in P, check if special case is necessary
    PQY = PQ * Y'
    sumPQ = sum(PQ, dims = 2)
    @. G += 4(sumPQ' * Y - PQY')
    return G
end

function repulsive_force!(G::AbstractMatrix, P::AbstractMatOrFac, Q²::AbstractMatOrFac, Y::AbstractMatrix, Z::Real)
    n, m = size(P) # is square
    # we compute in one product
    # 1. row-wise sum of Q²
    # 2. for each dimension, run fast transform once -> d * n * log(n) complexity
    Y1 = hcat(Y', ones(m)) # IDEA: lazy or re-allocate?
    Q²Y1 = Q²*Y1 # does all relevant computations in one matrix-matrix product
    Q²Y = @view Q²Y1[:, 1:end-1]
    sumQ² = @view Q²Y1[:, end]
    @. G += 4/Z * (Q²Y' - sumQ²' * Y) # can this be re-written as a single product?
    return G
end

################################## ADAM ########################################
# implementation of adam optimizer
# gradient first and second moment increment using
function gradient!(G::AbstractMatrix, G²::AbstractMatrix,
                    P::AbstractMatOrFac, Q::AbstractMatOrFac, Q²::AbstractMatOrFac,
                    Y::AbstractMatrix, Z::Real, β₁::Real = 0, β₂::Real = 0)
    @. G *= β₁
    @. G² *= β₂
    attractive_force!(G, G², P, Q, Y, β₁, β₂) # for t-SNE
    repulsive_force!(G, G², P, Q², Y, Z, β₁, β₂)
    return G, G²
end

# generic fallback
function attractive_force!(G::AbstractMatrix, G²::AbstractMatrix, P::AbstractMatOrFac,
                           Q::AbstractMatOrFac, Y::AbstractMatrix, β₁::Real, β₂::Real)
    PQ = hadamard_product(P, Q) # fast through sparsity in P,  check if special case is necessary
    PQY = PQ * Y'
    sumPQ = sum(PQ, dims = 2)
    ∇ = @. 4(sumPQ' * Y - PQY')
    @. G += (1-β₁) * ∇ # can this be re-written as a single product?
    @. G² += (1-β₂) * ∇^2
    return G
end

function repulsive_force!(G::AbstractMatrix, G²::AbstractMatrix, P::AbstractMatOrFac,
                Q²::AbstractMatOrFac, Y::AbstractMatrix, Z::Real, β₁::Real, β₂::Real)
    n, m = size(P) # is square
    # we compute in one product
    # 1. row-wise sum of Q²
    # 2. for each dimension, run fast transform once -> d * n * log(n) complexity
    Y1 = hcat(Y', ones(m)) # IDEA: lazy or re-allocate?
    Q²Y1 = Q²*Y1 # does all relevant computations in one matrix-matrix product
    Q²Y = @view Q²Y1[:, 1:2]
    sumQ² = @view Q²Y1[:, 3]
    # @. G += 4/Z * (Q²Y' - sumQ²' * Y) # can this be re-written as a single product?
    # @. G² += (4/Z * (Q²Y' - sumQ²' * Y))^2
    ∇ = @. 4/Z * (Q²Y' - sumQ²' * Y)
    @. G += (1-β₁) * ∇ # can this be re-written as a single product?
    @. G² += (1-β₂) * ∇^2
    return G
end

# stepsize to debias moments for Adam
stepsize(α::Real, β₁::Real, β₂::Real, t::Int) = α * sqrt(1-β₂^t)/(1-β₁^t)


text:
using NAssets

config = ( seed = 123
            ,ctl_model= NAssets.GraphModel(1)
            ,ntw_topo = NAssets.GraphModel(4)
            ,size=16
            ,n_steps=40
            ,drop_proportion=0.1
            ,prob_random_walks = 0.1
            ,benchmark = false
            ,animation = true
            ,k=0
            ,Β = 0.
            ,ctl_k=0
            ,ctl_Β = 0.
            ,custom_topo = nothing
            ,ctl_custom_topo = nothing
            ,mnt_policy = 1
            ,ntw_services = [(1,7),(4,1),(5,14),(12,8)]
            ,mnt_wc_duration = 20
            ,mnt_bc_duration = 10
            ,mnt_wc_cost = 150.
            ,mnt_bc_cost = 100.
            ,traffic_dist_params = [1,0.05]
            ,data_dir = "data/testing/"
            ,plots_dir = "plots/testing/"
        )
single_run_with_logging(config)


text:
# Adapted from
# https://github.com/JuliaData/AbstractTables.jl/blob/e5afc569504ecf08ec769fd52a78da4027eab35f/src/AbstractTable/show.jl

let
    local io = IOBuffer(Array(UInt8, 80), true, true)
    global ourstrwidth
    function ourstrwidth(x::Any)::Int
        truncate(io, 0)
        ourshowcompact(io, x)
        return position(io)
    end
    ourstrwidth(x::AbstractString) = strwidth(x) + 2 # -> Int
    # myconv = VERSION < v"0.4-" ? convert : Base.unsafe_convert
    ourstrwidth(s::Symbol) =
        @compat Int(
            ccall(:u8_strwidth, Csize_t, (Ptr{UInt8}, ),
            Base.unsafe_convert(Ptr{UInt8}, s))
        )
end

ourshowcompact(io::IO, x::Any)::Void = showcompact(io, x)
ourshowcompact(io::IO, x::AbstractString)::Void = showcompact(io, x)
ourshowcompact(io::IO, x::Symbol)::Void = print(io, x)

function getmaxwidths(tbl::AbstractTable, rowlabel, limit, offset)::Vector{Int}
    ncols = ncol(tbl)
    widths = [ Vector{Int}() for j in 1:ncols ]
    maxwidths = Array{Int}(ncol(tbl) + 1)
    undefstrwidth = ourstrwidth(Base.undef_ref_str)

    rows = eachrow(tbl)
    st = start(rows)
    i = 1
    while (i <= offset) & (!done(rows, st))
        row, st = next(rows, st)
        i += 1
    end
    while (i <= offset + limit) & (!done(rows, st))
        row, st = next(rows, st)
        for (j, v) in enumerate(row)
            try
                push!(widths[j], ourstrwidth(v))
            catch
                push!(widths[j], undefstrwidth)
            end
        end
        i += 1
    end
    flds = fields(tbl)
    for j in 1:ncols
        if isempty(widths[j])
            maxwidths[j] = ourstrwidth(flds[j])
        else
            maxwidths[j] = max(maximum(widths[j]), ourstrwidth(flds[j]))
        end
    end
    maxwidths[end] = max(ourstrwidth(rowlabel), ndigits(limit)+1)
    return maxwidths
end

function getprintedwidth(maxwidths::Vector{Int})::Int
    # Include length of line-initial |
    totalwidth = 1
    for i in 1:length(maxwidths)
        # Include length of field + 2 spaces + trailing |
        totalwidth += maxwidths[i] + 3
    end
    return totalwidth
end

pad(io, padding)::Void = print(io, " "^padding)

function print_bounding_line(io, maxwidths, j_left, j_right)::Void
    rowmaxwidth = maxwidths[end]
    print(io, '├')
    print(io, "─"^(rowmaxwidth + 2))
    print(io, '┼')
    for j in j_left:j_right
        print(io, "─"^(maxwidths[j] + 2))
        if j < j_right
            print(io, '┼')
        else
            print(io, '┤')
        end
    end
    print(io, '\n')
end

# NOTE: returns a Bool indicating whether or not there are more rows than what
# has printed
function print_tbl_rows(
    io, tbl, maxwidths, j_left, j_right, rowlabel, limit, offset
)::Bool
    rowmaxwidth = maxwidths[end]
    flds = fields(tbl)
    rows = eachrow(tbl, flds[j_left:j_right]...)
    st = start(rows)
    i = 1
    while (i <= offset) & (!done(rows, st))
        row, st = next(rows, st)
        i += 1
    end
    while (i <= offset + limit) & (!done(rows, st))
        row, st = next(rows, st)
        @printf(io, "│ %d", i)
        pad(io, rowmaxwidth - ndigits(i))
        print(io, " │ ")
        # print table entry
        for j in j_left:j_right
            v = row[j]
            strlen = ourstrwidth(v)
            ourshowcompact(io, v)
            pad(io, maxwidths[j] - strlen)
            if j == j_right
                if i == (limit+offset)
                    print(io, " │")
                else
                    print(io, " │\n")
                end
            else
                print(io, " │ ")
            end
        end
        i += 1
    end
    return (i > (limit+offset)) & (!done(rows, st))
end

function print_tbl_footer(
    io, tbl::AbstractTable, ::RowDimUnknown, more_rows, j_right, limit, offset
)::Void
    println(io, "\n⋮")
    ncols = ncol(tbl)
    if j_right < ncols
        flds = fields(tbl)
        @printf(io, "with more rows and %d more columns: ", ncols - j_right)
        for j in j_right+1:(ncol-1)
            field = flds[j]
            print(io, "$field, ")
        end
        print(io, "$field.")
    else
        print(io, "with more rows.")
    end
end

function print_tbl_footer(
    io, tbl::AbstractTable, ::HasRowDim, more_rows, j_right, limit, offset
)::Void
    ncols = ncol(tbl)
    if more_rows
        println(io, "\n⋮")
        if j_right < ncols
            if offset > 0
                @printf(
                    io,
                    "with %d more rows (skipped the first %d rows) and %d more columns: ",
                    nrow(tbl)-(limit+offset),
                    offset,
                    ncols - j_right
                )
            else
                @printf(
                    io,
                    "with %d more rows and %d more columns: ",
                    nrow(tbl)-limit,
                    ncols - j_right
                )
            end
            flds = fields(tbl)
            for j in j_right+1:(ncols-1)
                field = flds[j]
                print(io, "$field, ")
            end
            field = flds[ncols]
            print(io, "$field.")
        else
            if offset > 0
                @printf(
                    io,
                    "with %d more rows (skipped the first %d rows).",
                    nrow(tbl)-(limit+offset),
                    offset
                )
            else
                @printf(io, "with %d more rows.", nrow(tbl)-limit)
            end
        end
    else
        if j_right < ncols
            println(io, "\n⋮")
            @printf(io, "with %d more columns: ", ncols - j_right)
            flds = fields(tbl)
            for j in j_right+1:(ncols-1)
                field = flds[j]
                print(io, "$field, ")
            end
            field = flds[ncols]
            print(io, "$field.")
        end
    end
    return
end

function print_tbl_header(io, tbl, maxwidths, j_left, j_right, rowlabel)::Void
    rowmaxwidth = maxwidths[end]
    flds = fields(tbl)
    @printf(io, "│ %s", rowlabel)
    pad(io, rowmaxwidth - ourstrwidth(rowlabel))
    print(io, " │ ")
    for j in j_left:j_right
        field = flds[j]
        ourshowcompact(io, field)
        pad(io, maxwidths[j] - ourstrwidth(field))
        j == j_right ? print(io, " │\n") : print(io, " │ ")
    end
    print_bounding_line(io, maxwidths, j_left, j_right)
    return
end

function getchunkbounds(
    io, maxwidths::Vector{Int}, splitchunks::Bool,
    availablewidth::Int=displaysize(io)[2]
)::Vector{Int}
    ncols = length(maxwidths) - 1
    rowmaxwidth = maxwidths[end]
    if splitchunks
        chunkbounds = [0]
        # Include 2 spaces + 2 | characters for row/col label
        totalwidth = rowmaxwidth + 4
        for j in 1:ncols
            # Include 2 spaces + | character in per-column character count
            totalwidth += maxwidths[j] + 3
            if totalwidth > availablewidth
                push!(chunkbounds, j - 1)
                totalwidth = rowmaxwidth + 4 + maxwidths[j] + 3
            end
        end
        push!(chunkbounds, ncols)
    else
        chunkbounds = [0, ncols]
    end
    return chunkbounds
end

# 1 space for line-initial | + length of field + 2 spaces + trailing |
printedwidth(maxwidths)::Void = foldl((x,y)->x+y+3, 1, maxwidths)

function show_ci_tbl(
    io, tbl::AbstractTable, rowlabel = :Row, displaysummary = true,
    splitchunks = true, showall = false, limit = 10, offset = 0
)::Void
    ncols = ncol(tbl)
    if ncols > 0
        displaysummary && println(io, summary(tbl))
        maxwidths = getmaxwidths(tbl, rowlabel, limit, offset)
        chunkbounds = getchunkbounds(io, maxwidths, splitchunks | !showall)
        if !showall
            j_left, j_right = 1, chunkbounds[2]
            print_tbl_header(io, tbl, maxwidths, j_left, j_right, rowlabel)
            more_rows = print_tbl_rows(
                io, tbl, maxwidths, j_left, j_right, rowlabel, limit, offset
            )
            print_tbl_footer(
                io, tbl, tblrowdim(tbl), more_rows, j_right, limit, offset
            )
        else
            nchunks = length(chunkbounds) - 1
            more_rows = false
            for r in 1:nchunks
                j_left = chunkbounds[r] + 1
                j_right = chunkbounds[r + 1]

                print_tbl_header(io, tbl, maxwidths, j_left, j_right, rowlabel)
                more_rows |= print_tbl_rows(
                    io, tbl, maxwidths, j_left, j_right, rowlabel, limit, offset
                )
            end
            print_tbl_footer(
                io, tbl, tblrowdim(tbl), more_rows, false, limit, offset
            )
        end
    else
        @printf(io, "An empty %s", typeof(tbl))
    end
    return
end

show_ci_tbl(tbl::AbstractTable, limit, offset)::Void =
    show(STDOUT, tbl, :Row, true, limit, offset)


text:
using ArgParse

include("functions_closest.jl")

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table! s begin
        "--target"
            help = "The alignment to search for closest matches, in fasta format"
            required = true
        "--query"
            help = "The sequence(s) to find matches for, in fasta format"
            required = true
        "--reference", "-r"
            help = "The reference sequence, in fasta format"
            required = false
        "--csv-out", "-o"
            help = "The name of csv file containing closest matches to write"
            default = "out.csv"
        "--fast"
            help = "Restrict comparison to sites that are different from the reference.
                    Requires --reference to be provided; doesn't calculate a per-site distance"
            action = :store_true
    end

    return parse_args(s)
end

function main()

    parsed_args = parse_commandline()

    if parsed_args["fast"]
        if parsed_args["reference"] == nothing
            e = error("need to provide --reference if you use --fast")
            throw(e)
        end

        fast_closest(parsed_args["target"], parsed_args["query"], parsed_args["reference"], parsed_args["csv-out"])
    else
        closest(parsed_args["target"], parsed_args["query"], parsed_args["csv-out"])
    end
end

main()

max_stars_repo_name: cov-ert/julialign

max_issues_repo_name: cov-ert/julialign

max_forks_repo_name: cov-ert/julialign


text:
using ProgressMeter

  """
      bootstrapPosteriorEstimation( errMat::Matrix{Float64}, samplingBatchSize::Int64, nrRuns::Int64 )

      compute posterior p( h* = h | S ).
      #Arguments
      - `errMat::Matrix{Float64}}`: each column is the prediction error of one hypothesis.
      - `samplingBatchSize::Int64`: sample size per main iteration.
      - `nrRuns::Int64`:            number of passes over predictions.
      #Return
      - `Vector{Float64}`:          Distribution p( h* = h | S ).
  """
  @views function bootstrapPosteriorEstimation( errMat::Matrix{Float64}, samplingBatchSize::Int64, nrRuns::Int64 )
    len                = size( errMat, 1 );
    width              = size( errMat, 2 );
    res                = zeros( Int, width );
    samplingIndexCache = zeros( Int, samplingBatchSize );
    @showprogress 1 "Computing..." for i in 1:1:nrRuns
      samplingIndexCache    .= rand( 1:size( errMat, 1 ), samplingBatchSize ) ;
      @inbounds label        = argminUniProb.argminProb( mean( errMat[samplingIndexCache,:], dims=1 )[1,:] );
      @inbounds res[label]  += 1;
    end
    return  res ./ nrRuns;
  end

  """
      bootstrapPosteriorEstimation!( errMat::Matrix{Float64}, samplingBatchSize::Int64, nrRuns::Int64, p::Array{Float64} )

      compute posterior p( h* = h | S ).
      #Arguments
      - `errMat::Matrix{Float64}}`: each column is the prediction error of one hypothesis.
      - `samplingBatchSize::Int64`: sample size per main iteration.
      - `nrRuns::Int64`:            number of passes over predictions.
      - `p::Vector{Float64}`:       resulting posterior p( h* = h | S ).
      #Return
      - `nothing`:                  nothing.
  """
  function bootstrapPosteriorEstimation!( errMat::Matrix{Float64}, samplingBatchSize::Int64, nrRuns::Int64, p::Vector{Float64} )
    pB = bootstrapPosteriorEstimation( Matrix( errMat ), samplingBatchSize, nrRuns );
    for (i,val) in enumerate( pB )
      p[i] = val;
    end
  end


text:
module ManyBody

include("term.jl")

end


text:
using MicroLogging

function checkerror(err)
    if err[] != C_NULL
        funcname = unsafe_string(rs2_get_failed_function(err[]))
        funcargs = unsafe_string(rs2_get_failed_args(err[]))
        message = unsafe_string(rs2_get_error_message(err[]))
        @error "rs_error was raised when calling $funcname($funcargs):\n    $message"
    end
end

function deviceinfo(dev)
    err = Ref{Ptr{rs2_error}}(0)
    device_info = rs2_get_device_info(dev, RS2_CAMERA_INFO_NAME, err)
    checkerror(err)
    @info "This device is an $(unsafe_string(device_info))"
    device_info = rs2_get_device_info(dev, RS2_CAMERA_INFO_SERIAL_NUMBER, err)
    checkerror(err)
    @info "  Serial number: $(unsafe_string(device_info))"
    device_info = rs2_get_device_info(dev, RS2_CAMERA_INFO_FIRMWARE_VERSION, err)
    checkerror(err)
    @info "  Firmware version: $(unsafe_string(device_info))"
    device_info = rs2_get_device_info(dev, RS2_CAMERA_INFO_PHYSICAL_PORT, err)
    checkerror(err)
    @info "  Physical port: $(unsafe_string(device_info))"
    device_info = rs2_get_device_info(dev, RS2_CAMERA_INFO_ADVANCED_MODE, err)
    checkerror(err)
    @info "  Advanced mode: $(unsafe_string(device_info))"
    device_info = rs2_get_device_info(dev, RS2_CAMERA_INFO_PRODUCT_ID, err)
    checkerror(err)
    @info "  Product ID: $(unsafe_string(device_info))"
end

function get_depth_frame(frames)
    # get the number of frames embedded within the composite frame
    num_of_frames = rs2_embedded_frames_count(frames, err)
    checkerror(err)
    for i = 0:num_of_frames-1
        # extract frame
        frame = rs2_extract_frame(frames, i, err)
        checkerror(err)
        # extract stream profile
        profile = rs2_get_frame_stream_profile(frame, err)
        checkerror(err)
        stream, format = Ref{rs2_stream}(0), Ref{rs2_format}(0)
        index, unique_id, framerate = Ref{Cint}(0), Ref{Cint}(0), Ref{Cint}(0)
        rs2_get_stream_profile_data(profile, stream, format, index, unique_id, framerate, err)
        checkerror(err)
        stream[] == RS2_STREAM_DEPTH && return frame
    end
    return C_NULL
end


text:
"""
Library of electrical models.
This library contains electrical components to build up analog circuits.
"""
module Electrical

using ModelingToolkit, Symbolics, IfElse, OrdinaryDiffEq
using OffsetArrays

@parameters t
D = Differential(t)

using ..Blocks: RealInput, RealOutput

include("utils.jl")
include("Analog/ideal_components.jl")
include("Analog/sensors.jl")
include("Analog/sources.jl")
# include("Digital/components.jl")
# include("Digital/gates.jl")
# include("Digital/tables.jl")
# include("Digital/sources.jl")

# TODO:
# - digital
# - machines
# - multi-phase

export #Interface
       Pin,
       # Analog Components
       Capacitor, Ground, Inductor, Resistor, Conductor,
       Short, IdealOpAmp,
       # Analog Sensors
       CurrentSensor, PotentialSensor, VoltageSensor,
       PowerSensor, MultiSensor,
       # Analog Sources
       Voltage, Current

       # # Digital Gates
       # And, Or, Not, Xor, Nand, Nor, Xnor,
       # # Digital components
       # HalfAdder, FullAdder, MUX, DEMUX, Encoder, Decoder,
       # # Digital Sources
       # DigitalPin, Pulse, PulseDiff

end


text:
"""
    Lloyd <: AbstractKMeansAlg

Basic algorithm for k-means calculation.
"""
struct Lloyd <: AbstractKMeansAlg end

"""
    Kmeans!(alg::AbstractKMeansAlg, containers, design_matrix, k; n_threads = nthreads(), k_init="k-means++", max_iters=300, tol=1e-6, verbose=true)

Mutable version of `kmeans` function. Definition of arguments and results can be
found in `kmeans`.

Argument `containers` represent algorithm specific containers, such as labels, intermidiate
centroids and so on, which are used during calculations.
"""
function kmeans!(alg::Lloyd, containers, X, k;
                n_threads = Threads.nthreads(),
                k_init = "k-means++", max_iters = 300,
                tol = 1e-6, verbose = false, init = nothing)
    nrow, ncol = size(X)
    centroids = init == nothing ? smart_init(X, k, n_threads, init=k_init).centroids : deepcopy(init)

    converged = false
    niters = 1
    J_previous = 0.0

    # Update centroids & labels with closest members until convergence
    while niters <= max_iters
        @parallelize n_threads ncol chunk_update_centroids(alg, containers, centroids, X)
        collect_containers(alg, containers, centroids, n_threads)
        J = sum(containers.J)

        if verbose
            # Show progress and terminate if J stopped decreasing.
            println("Iteration $niters: Jclust = $J")
        end

        # Check for convergence
        if (niters > 1) & (abs(J - J_previous) < (tol * J))
            converged = true
            break
        end

        J_previous = J
        niters += 1
    end

    @parallelize n_threads ncol sum_of_squares(containers, X, containers.labels, centroids)
    totalcost = sum(containers.sum_of_squares)

    # Terminate algorithm with the assumption that K-means has converged
    if verbose & converged
        println("Successfully terminated with convergence.")
    end

    # TODO empty placeholder vectors should be calculated
    # TODO Float64 type definitions is too restrictive, should be relaxed
    # especially during GPU related development
    return KmeansResult(centroids, containers.labels, Float64[], Int[], Float64[], totalcost, niters, converged)
end

kmeans(design_matrix, k;
    n_threads = Threads.nthreads(),
    k_init = "k-means++", max_iters = 300, tol = 1e-6,
    verbose = false, init = nothing) =
        kmeans(Lloyd(), design_matrix, k; n_threads = n_threads, k_init = k_init, max_iters = max_iters, tol = tol,
            verbose = verbose, init = init)

"""
    create_containers(::Lloyd, k, nrow, ncol, n_threads)

Internal function for the creation of all necessary intermidiate structures.

- `centroids_new` - container which holds new positions of centroids
- `centroids_cnt` - container which holds number of points for each centroid
- `labels` - vector which holds labels of corresponding points
"""
function create_containers(::Lloyd, k, nrow, ncol, n_threads)
    lng = n_threads + 1
    centroids_new = Vector{Array{Float64,2}}(undef, lng)
    centroids_cnt = Vector{Vector{Int}}(undef, lng)

    for i in 1:lng
        centroids_new[i] = Array{Float64, 2}(undef, nrow, k)
        centroids_cnt[i] = Vector{Int}(undef, k)
    end

    labels = Vector{Int}(undef, ncol)

    J = Vector{Float64}(undef, n_threads)

    # total_sum_calculation
    sum_of_squares = Vector{Float64}(undef, n_threads)

    return (centroids_new = centroids_new, centroids_cnt = centroids_cnt,
            labels = labels, J = J, sum_of_squares = sum_of_squares)
end

function chunk_update_centroids(::Lloyd, containers, centroids, X, r, idx)
    # unpack containers for easier manipulations
    centroids_new = containers.centroids_new[idx]
    centroids_cnt = containers.centroids_cnt[idx]
    labels = containers.labels

    centroids_new .= 0.0
    centroids_cnt .= 0
    J = 0.0
    @inbounds for i in r
        min_dist = distance(X, centroids, i, 1)
        label = 1
        for j in 2:size(centroids, 2)
            dist = distance(X, centroids, i, j)
            label = dist < min_dist ? j : label
            min_dist = dist < min_dist ? dist : min_dist
        end
        labels[i] = label
        centroids_cnt[label] += 1
        for j in axes(X, 1)
            centroids_new[j, label] += X[j, i]
        end
        J += min_dist
    end

    containers.J[idx] = J
end

function collect_containers(alg::Lloyd, containers, centroids, n_threads)
    if n_threads == 1
        @inbounds centroids .= containers.centroids_new[1] ./ containers.centroids_cnt[1]'
    else
        @inbounds containers.centroids_new[end] .= containers.centroids_new[1]
        @inbounds containers.centroids_cnt[end] .= containers.centroids_cnt[1]
        @inbounds for i in 2:n_threads
            containers.centroids_new[end] .+= containers.centroids_new[i]
            containers.centroids_cnt[end] .+= containers.centroids_cnt[i]
        end

        @inbounds centroids .= containers.centroids_new[end] ./ containers.centroids_cnt[end]'
    end
end


text:
# This is part of Calendars.jl. See the copyright note there.
# =================== Interactive calendar ==================

# Interactively compute calendar dates in the REPL.

using Dates: yearmonthday, now

"""
Interactively query the dates for all calendars in the REPL.

```julia
julia> IDate()
```
"""
function IDate()

    msg = "You should have entered a numeric value, bye ..."
    noi = "This was not a valid input! Try again ..."
    nod = "This was not a valid date! Try again ..."

    date = yearmonthday(now())
    cd = CalendarDates("EC", date, true)
    dy = DayOfYear("EC", date)
    wd = WeekDay(cd)
    println("\nIt's ", wd, " the ", dy, "-th day of the EC-year.")
    println()

    while true
    @label start
        println("Enter a calendar specifier:")
        println("    EC => European")
        println("    CE => Common  ")
        println("    JD => Julian  ")
        println("    AM => Hebrew  ")
        println("    AH => Islamic ")
        println("    ID => IsoDate ")
        println("or choose a historical event:")
        println("    TV => Treaty of Verdun            ")
        println("    SC => Sack of Constantinople      ")
        println("    LV => Birth of Leonardo da Vinci  ")
        println("    SR => Start of the Reformation    ")
        println("    FL => Great Fire of London        ")
        println("    SB => The Storming of the Bastille")
        println("or Ctrl+C to quit.")
        println()

        ct = try
            readline()
        catch e
            println("... bye")
            return
        end

        if ct in ["EC", "CE", "JD", "AM", "AH", "ID"]
            cs = ct

            print("Enter year      (>0): ")
            y = try
                parse(Int, readline())
            catch e
                println(msg)
                return
            end
            if cs == "ID"
                print("Enter week   (1..53): ")
            else
                print("Enter month  (1..12): ")
            end
            m = try
                parse(Int, readline())
            catch e
                println(msg)
                return
            end
            if cs == "ID"
                print("Enter day     (1..7): ")
            else
                print("Enter day    (1..31): ")
            end
            d = try
                parse(Int, readline())
            catch e
                println(msg)
                return
            end
            println()
            date = (StringToSymbol(cs), y, m, d)

        elseif ct in ["TV", "SC", "LV", "SR", "FL", "SB"]
            ct == "TV" && (date = (EC, 843,8,10))  # Treaty of Verdun
            ct == "SC" && (date = (EC,1204,4,12))  # Sack of Constantinople
            ct == "LV" && (date = (EC,1452,4,15))  # Birth of Leonardo da Vinci
            ct == "SR" && (date = (EC,1517,10,31)) # Start of the Reformation
            ct == "FL" && (date = (EC,1666,9,2))   # Great Fire of London
            ct == "SB" && (date = (EC,1789,7,14))  # The Storming of the Bastille
            cs = "EC"
        else
            println(noi)
            @goto start
        end

        if ! isValidDate(date)
           println(nod)
        else
            wd = " "
            try
                cd = CalendarDates(date, true)
                wd = WeekDay(cd)
            catch e
                println(e)
                return
            end
            println("\n", wd, "the ",
                   DayOfYear(date), "-th day of the ",
                   cs, "-year.")
        end
        println()
    end
end


text:
struct MouseEvent
    x::Real
    y::Real
    button::Cuint
    state::Cuint
end

struct ScrollEvent
    x::Real
    y::Real
    direction::GEnum
    state::Cuint
    dx::Float64
    dy::Float64
end

struct CanvasEvents
    mouseup::Observable{MouseEvent}
    mousedown::Observable{MouseEvent}
    mousemove::Observable{MouseEvent}
    mousedrag::Observable{MouseEvent}
    scroll::Observable{ScrollEvent}
    ids::Vector{Culong}
    widget::GtkWidget
    function CanvasEvents(canvas::GtkCanvas)
        event = MouseEvent(-1, -1, 0, 0)
        scroll = ScrollEvent(-1, -1, 0, 0, 0, 0)
        ids = Culong[]
        this = new(event, event, event, event, scroll, ids, canvas)
        push!(ids, Gtk.on_signal_button_press(mousedown_cb, canvas, false, this))
        push!(ids, Gtk.on_signal_button_release(mouseup_cb, canvas, false, this))
        push!(ids, Gtk.on_signal_motion(mousemove_cb, canvas, 0, 0, false, this))
        push!(ids, Gtk.on_signal_scroll(mousescroll_cb, canvas, false, this))
        return this
    end
end

mutable struct Canvas <: Widget{GtkCanvas}
    widget::GtkCanvas
    mouse::CanvasEvents
    width::Real
    height::Real
    tickcb::Cuint
    function Canvas(width::Integer = -1, height::Integer = -1; props...)
        gcanvas = @widget GtkCanvas(width, height)
        for id in gcanvas.mouse.ids
            disconnect(gcanvas, id)
        end
        empty!(gcanvas.mouse.ids)
        mouse = CanvasEvents(gcanvas)
        canvas = new(gcanvas, mouse, -1, -1, 0)
        canvas["is-focus"] = true
        gcpreserve(gcanvas, canvas)
        return canvas
     end
end

sized(canvas::Canvas) = canvas.widget.is_sized

function Base.getproperty(canvas::Canvas, prop::Symbol)
    prop == :width && return width(getfield(canvas, 1))
    prop == :height && return height(getfield(canvas, 1))
    return getfield(canvas, prop)
end

events(canvas::Canvas) = getfield(canvas, :mouse)

function Canvas(update::Function; loop::Bool = true, width = -1, height = -1, props...)
    canvas = Canvas(width, height; props...)
    if loop
        onupdate(update, canvas)
    else
        ondraw(update, canvas)
    end
    return canvas
end

function mousedown_cb(ptr::Ptr, eventp::Ptr, this::CanvasEvents)
    evt = unsafe_load(eventp)
    this.mousedown[] = MouseEvent(evt.x, evt.y, evt.button, evt.state)
    return Cint(false)
end

function mouseup_cb(ptr::Ptr, eventp::Ptr, this::CanvasEvents)
    evt = unsafe_load(eventp)
    this.mouseup[] = MouseEvent(evt.x, evt.y, evt.button, evt.state)
    return Cint(false)
end

function mousemove_cb(ptr::Ptr, eventp::Ptr, this::CanvasEvents)
    evt = unsafe_load(eventp)
    button = 0
    this.mousemove[] = MouseEvent(evt.x, evt.y, 0, evt.state)
    if evt.state & Gtk.GdkModifierType.BUTTON1 != 0
        button = 1
    elseif evt.state & Gtk.GdkModifierType.BUTTON2 != 0
        button = 2
    elseif evt.state & Gtk.GdkModifierType.BUTTON3 != 0
        button = 3
    end
    if button != 0
        this.mousedrag[] = MouseEvent(evt.x, evt.y, button, evt.state)
    end
    Cint(false)
end

function mousescroll_cb(ptr::Ptr, eventp::Ptr, this::CanvasEvents)
    evt = unsafe_load(eventp)
    this.scroll[] = ScrollEvent(
        evt.x,
        evt.y,
        evt.direction == 1 ? -1 : 1,
        evt.state,
        evt.delta_x,
        evt.delta_y
    )
    Cint(false)
end


text:
#!/usr/bin/env julia

#https://github.com/joshbode/zmq-example/blob/master/julia-server.jl

using JSON
using ZMQ
using Logging

BASE_PORT = int(ARGS[1])
WORKER = int(ARGS[2])
PORT = BASE_PORT + WORKER

Logging.configure(level=Logging.DEBUG)
context = ZMQ.Context()
socket = ZMQ.Socket(context, ZMQ.REP)
ZMQ.bind(socket, "tcp://*:$PORT")
Logging.info("Ready")
while true
    data =takebuf_string(convert(IOStream, ZMQ.recv(socket)))
    Logging.info("Received: [", join(data, ", "), "]")
    result = ["ok"]
    ZMQ.send(socket, JSON.json(result))
    Logging.info("Sent: $result")
end

max_stars_repo_name: lidiapopescu/JuliaDT

max_issues_repo_name: lidiapopescu/JuliaDT

max_forks_repo_name: JuliaComputing/JuliaDT


text:
using Causal

# Deifne model
@defmodel model begin
    @nodes begin
        gen = SinewaveGenerator(amplitude=1., frequency=1/2π)
        adder = Adder(signs=(+, -))
        ds = ContinuousLinearSystem(A=fill(-1., 1, 1), state=[1.])
        writer = Writer(input=Inport(2))
    end
    @branches begin
        gen[1] => adder[1]
        adder[1] => ds[1]
        ds[1] => adder[2]
        ds[1] => writer[1]
        gen[1] => writer[2]
    end
end

# Simulate the model
tinit, tsample, tfinal = 0., 0.01, 10.
sim = simulate!(model, tinit, tsample, tfinal)

# Read and plot data
t, x = read(getnode(model, :writer).component)
t, x = read(getnode(model, :writer).component)
using Plots
plot(t, x[:, 1], label="r(t)", xlabel="t")
plot!(t, x[:, 2], label="y(t)", xlabel="t")
plot!(t, 6 / 5 * exp.(-2t) + 1 / 5 * (2 * sin.(t) - cos.(t)), label="Analytical Solution")


text:
using DFWannier
using ProfileView
T = Float32
x = WannierModel{T}("/Users/ponet/Documents/Fysica/PhD/GeTe/colin/paperxsf/test1","/Users/ponet/Documents/Fysica/PhD/GeTe/fullrel/GeTe_bands.out",[[Atom(T[0.0,0.0,-0.0239129,-0.155854]...) for i=1:4]...,[Atom(T[0.0,0.0,5.5540692,0.318205]...) for i=1:4]...]);
@time tbbandssoc = calculate_eig_cm_angmom_soc(x,90:0.2:110);
Profile.clear()
using BenchmarkTools
@benchmark tbbandssoc = calculate_eig_cm_angmom_soc(x,90:0.2:110)
@profile tbbandssoc1 = calculate_eig_cm_angmom_soc(x,50:0.2:110);
@time calculate_angmom(x.wfcs[1],x.wfcs[3])
Profile.clear()
Profile.init(1000000000, 0.00100000)
ProfileView.view()
@code_warntype read_xsf_file("/Users/ponet/Documents/Fysica/PhD/GeTe/colin/paperxsf/test1")

struct perfTest
  val::Float64
  function perfTest()
    new(rand(Float64))
  end
end

struct perfTest2
  val::Float64
  t::Bool
  function perfTest2()
    new(rand(Float64),rand(Bool))
  end
end

const test_array = [perfTest() for i=1:10000000]
const test_array2 = [perfTest2() for i=1:10000000]

function bench_test(t)
  out = 0.0
  @inbounds for i=1:10000000
    out += t[i].val
  end
  out
end
# using StaticArrays
function bench_test2(t)
  out = 0.0
  v = Array{perfTest2,1}(12)
  @inbounds for i=1:10:10000000
    unsafe_copy!(v,1,t,i,12)
    @inbounds for j = 1:12
    out += v[j].val
    end
  end
  out
end
function bench_test2(t)
  out = 0.0
  @inbounds for i=1:10000000
    out += t[i].val
  end
  out
end

using BenchmarkTools
@code_native bench_test(test_array)
@code_native bench_test2(test_array2)


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: mediatailor
using AWS.Compat
using AWS.UUIDs

"""
    delete_playback_configuration(name)
    delete_playback_configuration(name, params::Dict{String,<:Any})

Deletes the playback configuration for the specified name.

# Arguments
- `name`: The identifier for the playback configuration.

"""
delete_playback_configuration(Name; aws_config::AbstractAWSConfig=global_aws_config()) = mediatailor("DELETE", "/playbackConfiguration/$(Name)"; aws_config=aws_config)
delete_playback_configuration(Name, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = mediatailor("DELETE", "/playbackConfiguration/$(Name)", params; aws_config=aws_config)

"""
    get_playback_configuration(name)
    get_playback_configuration(name, params::Dict{String,<:Any})

Returns the playback configuration for the specified name.

# Arguments
- `name`: The identifier for the playback configuration.

"""
get_playback_configuration(Name; aws_config::AbstractAWSConfig=global_aws_config()) = mediatailor("GET", "/playbackConfiguration/$(Name)"; aws_config=aws_config)
get_playback_configuration(Name, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = mediatailor("GET", "/playbackConfiguration/$(Name)", params; aws_config=aws_config)

"""
    list_playback_configurations()
    list_playback_configurations(params::Dict{String,<:Any})

Returns a list of the playback configurations defined in AWS Elemental MediaTailor. You can
specify a maximum number of configurations to return at a time. The default maximum is 50.
Results are returned in pagefuls. If MediaTailor has more configurations than the specified
maximum, it provides parameters in the response that you can use to retrieve the next
pageful.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"MaxResults"`: Maximum number of records to return.
- `"NextToken"`: Pagination token returned by the GET list request when results exceed the
  maximum allowed. Use the token to fetch the next page of results.
"""
list_playback_configurations(; aws_config::AbstractAWSConfig=global_aws_config()) = mediatailor("GET", "/playbackConfigurations"; aws_config=aws_config)
list_playback_configurations(params::AbstractDict{String, Any}; aws_config::AbstractAWSConfig=global_aws_config()) = mediatailor("GET", "/playbackConfigurations", params; aws_config=aws_config)

"""
    list_tags_for_resource(resource_arn)
    list_tags_for_resource(resource_arn, params::Dict{String,<:Any})

Returns a list of the tags assigned to the specified playback configuration resource.

# Arguments
- `resource_arn`: The Amazon Resource Name (ARN) for the playback configuration. You can
  get this from the response to any playback configuration request.

"""
list_tags_for_resource(ResourceArn; aws_config::AbstractAWSConfig=global_aws_config()) = mediatailor("GET", "/tags/$(ResourceArn)"; aws_config=aws_config)
list_tags_for_resource(ResourceArn, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = mediatailor("GET", "/tags/$(ResourceArn)", params; aws_config=aws_config)

"""
    put_playback_configuration()
    put_playback_configuration(params::Dict{String,<:Any})

Adds a new playback configuration to AWS Elemental MediaTailor.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"AdDecisionServerUrl"`: The URL for the ad decision server (ADS). This includes the
  specification of static parameters and placeholders for dynamic parameters. AWS Elemental
  MediaTailor substitutes player-specific and session-specific parameters as needed when
  calling the ADS. Alternately, for testing you can provide a static VAST URL. The maximum
  length is 25,000 characters.
- `"AvailSuppression"`: The configuration for Avail Suppression. Ad suppression can be used
  to turn off ad personalization in a long manifest, or if a viewer joins mid-break.
- `"Bumper"`: The configuration for bumpers. Bumpers are short audio or video clips that
  play at the start or before the end of an ad break.
- `"CdnConfiguration"`: The configuration for using a content delivery network (CDN), like
  Amazon CloudFront, for content and ad segment management.
- `"ConfigurationAliases"`: Predefined aliases for dynamic variables.
- `"DashConfiguration"`: The configuration for DASH content.
- `"LivePreRollConfiguration"`: The configuration for pre-roll ad insertion.
- `"ManifestProcessingRules"`: The configuration for manifest processing rules. Manifest
  processing rules enable customization of the personalized manifests created by MediaTailor.
- `"Name"`: The identifier for the playback configuration.
- `"PersonalizationThresholdSeconds"`: The maximum duration of underfilled ad time (in
  seconds) allowed in an ad break.
- `"SlateAdUrl"`: The URL for a high-quality video asset to transcode and use to fill in
  time that's not used by ads. AWS Elemental MediaTailor shows the slate to fill in gaps in
  media content. Configuring the slate is optional for non-VPAID configurations. For VPAID,
  the slate is required because MediaTailor provides it in the slots that are designated for
  dynamic ad content. The slate must be a high-quality asset that contains both audio and
  video.
- `"TranscodeProfileName"`: The name that is used to associate this playback configuration
  with a custom transcode profile. This overrides the dynamic transcoding defaults of
  MediaTailor. Use this only if you have already set up custom profiles with the help of AWS
  Support.
- `"VideoContentSourceUrl"`: The URL prefix for the master playlist for the stream, minus
  the asset ID. The maximum length is 512 characters.
- `"tags"`: The tags to assign to the playback configuration.
"""
put_playback_configuration(; aws_config::AbstractAWSConfig=global_aws_config()) = mediatailor("PUT", "/playbackConfiguration"; aws_config=aws_config)
put_playback_configuration(params::AbstractDict{String, Any}; aws_config::AbstractAWSConfig=global_aws_config()) = mediatailor("PUT", "/playbackConfiguration", params; aws_config=aws_config)

"""
    tag_resource(resource_arn, tags)
    tag_resource(resource_arn, tags, params::Dict{String,<:Any})

Adds tags to the specified playback configuration resource. You can specify one or more
tags to add.

# Arguments
- `resource_arn`: The Amazon Resource Name (ARN) for the playback configuration. You can
  get this from the response to any playback configuration request.
- `tags`: A comma-separated list of tag key:value pairs. For example:
 {
 \"Key1\":
  \"Value1\",
 \"Key2\": \"Value2\"
 }

"""
tag_resource(ResourceArn, tags; aws_config::AbstractAWSConfig=global_aws_config()) = mediatailor("POST", "/tags/$(ResourceArn)", Dict{String, Any}("tags"=>tags); aws_config=aws_config)
tag_resource(ResourceArn, tags, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = mediatailor("POST", "/tags/$(ResourceArn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tags"=>tags), params)); aws_config=aws_config)

"""
    untag_resource(resource_arn, tag_keys)
    untag_resource(resource_arn, tag_keys, params::Dict{String,<:Any})

Removes tags from the specified playback configuration resource. You can specify one or
more tags to remove.

# Arguments
- `resource_arn`: The Amazon Resource Name (ARN) for the playback configuration. You can
  get this from the response to any playback configuration request.
- `tag_keys`: A comma-separated list of the tag keys to remove from the playback
  configuration.

"""
untag_resource(ResourceArn, tagKeys; aws_config::AbstractAWSConfig=global_aws_config()) = mediatailor("DELETE", "/tags/$(ResourceArn)", Dict{String, Any}("tagKeys"=>tagKeys); aws_config=aws_config)
untag_resource(ResourceArn, tagKeys, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = mediatailor("DELETE", "/tags/$(ResourceArn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tagKeys"=>tagKeys), params)); aws_config=aws_config)


text:
using NewPkgEval
using Test

v = NewPkgEval.build_julia("v1.2.0")

pkgnames = ["JSON", "TimerOutputs", "Crayons", "Example"]
pkgs = NewPkgEval.read_pkgs(pkgnames)
results = NewPkgEval.run(pkgs, 2, v)
for pkg in pkgnames
    @test results[pkg] == :ok
    output = read(joinpath(@__DIR__, "..", logs, "logs-"*string(v), pkg*".log"), String)
    @test occursin("Testing $pkg tests passed", output)
end
results = NewPkgEval.run(pkgs, 2, v; time_limit = 0.1)
for pkg in pkgnames
    @test results[pkg] == :fail
end

mktemp() do path, io
    NewPkgEval.run_sandboxed_julia(`-e 'print(1337)'`; ver=v"1.1.1", stdout=io)
    close(io)
    @test read(path, String) == "1337"
end


text:
# EXAMPLE: reduction
# here we have 3 ranks

using MPI, LinearAlgebra

MPI.Init()
comm = MPI.COMM_WORLD
rank = MPI.Comm_rank(comm)

if rank == 0
    data = [7.0,8.0,9.0,10.0,11.0,12.0]
else
    data = Float64[]
end

# partial data
data = MPI.Scatter(data, 2, 0, comm) #each rank has 2 elements
#partial dot
small_dot = dot(data,data)

#all dot
all_dot = MPI.Reduce(small_dot, +, 0, comm)

if rank==0
    @show all_dot
end

MPI.Finalize()

## output:
# all_dot = 559.0


text:
using Powersense
using GLPK
using Ipopt
using Test

# @testset "Unit tests" begin
#     include("unittests.jl")
# end

@testset "MathOptInterface" begin
    include("opt/MOI_wrapper.jl")
end

@testset "External Solver Attributes Implementation with Toy Example" begin
    include("opt/ext_solver.jl")
    @test typeof(optimizer_solver) == MOI.OptimizerWithAttributes
    @test isapprox(xsol, -1.0, rtol=1e-4)
    @test isapprox(ysol, -1.0, rtol=1e-4)
    @test status == MOI.LOCALLY_SOLVED
    include("../examples/opt/opt_example.jl")
    @test isapprox(xsol, -1.0, rtol=1e-4)
    @test isapprox(ysol, -1.0, rtol=1e-4)
end

@testset "ACOPF Formulations" begin
    include("../examples/opf/ACOPF_formulations_example.jl")
    include("opf/ACOPF_formulations.jl")
end


text:
using Compat.Test

import IJulia: helpmode, error_content, docdict

content = error_content(UndefVarError(:a))
@test "UndefVarError" == content["ename"]

@test haskey(docdict("import"), "text/plain")
@test haskey(docdict("sum"), "text/plain")


text:
using Test
using ModelingToolkit
@variables t

@connector function TwoPhaseFluidPort(;name, P=0.0, m_flow=0.0, h_outflow=0.0)
    vars = @variables h_outflow(t)=h_outflow [connect=Stream] m_flow(t)=m_flow [connect=Flow] P(t)=P
    ODESystem(Equation[], t, vars, []; name=name)
end

function MassFlowSource_h(;name,
        h_in=420e3,
        m_flow_in=-0.01,
    )

    pars = @parameters begin
        h_in=h_in
        m_flow_in=m_flow_in
    end

    vars = @variables begin
        P(t)
    end

    @named port = TwoPhaseFluidPort()

    subs = [port]

    eqns = Equation[]

    push!(eqns, port.P ~ P)
    push!(eqns, port.m_flow ~ -m_flow_in)
    push!(eqns, port.h_outflow ~ h_in)

    compose(ODESystem(eqns, t, vars, pars; name=name), subs)
end

# Simplified components.
function AdiabaticStraightPipe(;name,
        kwargs...,
    )

    vars = []
    pars = []

    @named port_a = TwoPhaseFluidPort()
    @named port_b = TwoPhaseFluidPort()

    subs = [port_a; port_b]

    eqns = Equation[]

    #=
    push!(eqns, port_a.P ~ port_b.P)
    push!(eqns, 0 ~ port_a.m_flow + port_b.m_flow)
    push!(eqns, port_b.h_outflow ~ instream(port_a.h_outflow))
    push!(eqns, port_a.h_outflow ~ instream(port_b.h_outflow))
    =#

    push!(eqns, connect(port_a, port_b))
    sys = ODESystem(eqns, t, vars, pars; name=name)
    sys = compose(sys, subs)
end

function SmallBoundary_Ph(;name,
        P_in=1e6,
        h_in=400e3,
    )

    vars = []

    pars = @parameters begin
        P=P_in
        h=h_in
    end

    @named port1 = TwoPhaseFluidPort()

    subs = [port1]

    eqns = Equation[]

    push!(eqns, port1.P ~ P)
    push!(eqns, port1.h_outflow ~ h)

    compose(ODESystem(eqns, t, vars, pars; name=name), subs)
end

# N1M1 model and test code.
function N1M1(;name,
        P_in=1e6,
        h_in=400e3,
        kwargs...,
    )

    @named port_a = TwoPhaseFluidPort()
    @named source = SmallBoundary_Ph(P_in=P_in, h_in=h_in)

    subs = [port_a; source]

    eqns = Equation[]

    push!(eqns, connect(source.port1, port_a))

    sys = ODESystem(eqns, t, [], [], name=name)
    sys = compose(sys, subs)
end

@named n1m1 = N1M1()
@named pipe = AdiabaticStraightPipe()
@named sink = MassFlowSource_h(m_flow_in=-0.01, h_in=400e3)

streams_a = [n1m1.port_a, pipe.port_a]
streams_b = [pipe.port_b, sink.port]

eqns = [
        connect(n1m1.port_a, pipe.port_a)
        connect(pipe.port_b, sink.port)
       ]

@named sys = ODESystem(eqns, t)
@named n1m1Test = compose(sys, n1m1, pipe, sink)
@test_nowarn structural_simplify(n1m1Test)

# N1M2 model and test code.
function N1M2(;name,
        P_in=1e6,
        h_in=400e3,
        kwargs...,
    )

    @named port_a = TwoPhaseFluidPort()
    @named port_b = TwoPhaseFluidPort()

    @named source = SmallBoundary_Ph(P_in=P_in, h_in=h_in)

    subs = [port_a; port_b; source]

    eqns = Equation[]

    push!(eqns, connect(source.port1, port_a))
    push!(eqns, connect(source.port1, port_b))

    sys = ODESystem(eqns, t, [], [], name=name)
    sys = compose(sys, subs)
end

@named n1m2 = N1M2()
@named sink1 = MassFlowSource_h(m_flow_in=-0.01, h_in=400e3)
@named sink2 = MassFlowSource_h(m_flow_in=-0.01, h_in=400e3)

eqns = [
        connect(n1m2.port_a, sink1.port)
        connect(n1m2.port_b, sink2.port)
       ]

@named sys = ODESystem(eqns, t)
@named n1m2Test = compose(sys, n1m2, sink1, sink2)
@test_nowarn structural_simplify(n1m2Test)

@named n1m2 = N1M2()
@named pipe1 = AdiabaticStraightPipe()
@named pipe2 = AdiabaticStraightPipe()
@named sink1 = MassFlowSource_h(m_flow_in=-0.01, h_in=400e3)
@named sink2 = MassFlowSource_h(m_flow_in=-0.01, h_in=400e3)

eqns = [
        connect(n1m2.port_a, pipe1.port_a)
        connect(pipe1.port_b, sink1.port)

        connect(n1m2.port_b, pipe2.port_a)
        connect(pipe2.port_b, sink2.port)
       ]

@named sys = ODESystem(eqns, t)
@named n1m2AltTest = compose(sys, n1m2, pipe1, pipe2, sink1, sink2)
@test_nowarn structural_simplify(n1m2AltTest)

# N2M2 model and test code.
function N2M2(;name,
        kwargs...,
    )

    @named port_a = TwoPhaseFluidPort()
    @named port_b = TwoPhaseFluidPort()
    @named pipe = AdiabaticStraightPipe()

    streams_a = [port_a, pipe.port_a]
    streams_b = [pipe.port_b, port_b]

    subs = [port_a; port_b; pipe]

    eqns = Equation[]

    push!(eqns, connect(port_a, pipe.port_a))
    push!(eqns, connect(pipe.port_b, port_b))

    sys = ODESystem(eqns, t, [], [], name=name)
    sys = compose(sys, subs)
end

@named n2m2 = N2M2()
@named source = MassFlowSource_h(m_flow_in=-0.01, h_in=400e3)
@named sink = SmallBoundary_Ph(P_in=1e6, h_in=400e3)

eqns = [
        connect(source.port, n2m2.port_a)
        connect(n2m2.port_b, sink.port1)
       ]

@named sys = ODESystem(eqns, t)
@named n2m2Test = compose(sys, n2m2, source, sink)
@test_nowarn structural_simplify(n2m2Test)

# array var
@connector function VecPin(;name)
    sts = @variables v[1:2](t)=[1.0,0.0] i[1:2](t)=1.0 [connect = Flow]
    ODESystem(Equation[], t, [sts...;], []; name=name)
end

@named vp1 = VecPin()
@named vp2 = VecPin()
@named vp3 = VecPin()

@named simple = ODESystem([connect(vp1, vp2, vp3)], t)
sys = expand_connections(compose(simple, [vp1, vp2, vp3]))
@test sort(equations(sys), by=string) == sort([
                          vp1.v[1] ~ vp2.v[1]
                          vp1.v[2] ~ vp2.v[2]
                          vp1.v[1] ~ vp3.v[1]
                          vp1.v[2] ~ vp3.v[2]
                          0 ~ -vp1.i[1] - vp2.i[1] - vp3.i[1]
                          0 ~ -vp1.i[2] - vp2.i[2] - vp3.i[2]
                        ], by=string)


text:
using LinearAlgebra, ForwardDiff

function NewtonProjeto(f, x)
    ∇f(x) = ForwardDiff.gradient(f, x)
    H(x) = ForwardDiff.hessian(f, x)
    m, n = size(H(x))
    η = 1e-2
    Β = 1.0e-3
    num_backtrack = 0
    iter = 0
    while norm(∇f(x)) > 1e-4
        h = H(x)

        if minimum(diag(h)) > 0
            ρ = 0.0
        else
            ρ = -minimum(diag(h)) + Β
        end

        k = 0
        while   issuccess(cholesky(h + ρ*I, check=false)) == false
                ρ = max(2ρ, Β)
                k = k+1
        end

        F = cholesky(h + ρ*I)
        J = F.L
        y = - J \ ∇f(x)
        M = J'
        d = M \ y

        # Armijo
        α = 1.0
        while f(x + α * d) ≥ f(x) + η * α * dot(d, ∇f(x))
            α = α / 2
            num_backtrack = num_backtrack + 1
            if α < 1e-8
                error("Erro no backtracking")
            end
        end
        x = x + α * d

        iter += 1
        if iter > 10000
            error("Nao converge")
        end
    end

    return x, iter, num_backtrack
end


text:
using LinearAlgebra
using Test

using SDPA

include("example1.jl")
include("example2.jl")

include("MOI_wrapper.jl")


text:
#Define constants for US76 atmosphere model
g0us76(::Type{Earth}) = 9.80665
r0us76(::Type{Earth}) = 6356.766
rsus76(::Type{Earth}) = 8314.32
naus76(::Type{Earth}) = Float64(6.022169*10^26)
m0us76(::Type{Earth}) = 28.9644
r0mus76(::Type{Earth}) = 6356755.
ρ0us76(::Type{Earth}) = 1.225 #kg/m^3
Rus76(::Type{Earth}) = 287 #J/Kg K For Earth Air
γus76(::Type{Earth}) = 1.4 #For Earth Air
#const r = rs/m0

#Starting Temperature per Interval (per airlayer)
tmsus76(::Type{Earth}) = [288.15, 216.65, 216.65, 228.65, 270.65, 270.65, 214.65, 186.87]

#dT/dG (Thermal Lapse rate) per intercal (tropo-, strato-, meso-, and thermosphere)
tgus76(::Type{Earth}) = [-6.5, 0., 1., 2.8, 0., -2.8, -2., 0.]

#Starting Pressure per interval
psus76(::Type{Earth}) = [101325., 22632., 5474.87, 868.014, 110.905, 66.938, 3.9564, .39814]
m_m0us76(::Type{Earth}) = [1., .999996, .999989, .999971, .999941, .999909, .99987, .999829,
                             .999786, .999741, .999694, .999641, .999579]

#Volume Percentage of element in Air, He, ...., N2, O2, Ar
fius76(::Type{Earth}) = [5.24e-6, 0.0, .78084, .209476, .00934]
zssus76(::Type{Earth}) = [86.,87.,88.,89.,90.,91.,92.,93.,94.,95.,96.,
                97.,98.,99.,100.,101.,102.,103.,104.,105.,106.,107.,108.,
                109.,110.,111.,112.,113.,114.,115.,116.,117.,118.,119.,120.]
dnssus76(::Type{Earth}) = [1.447e20,1.212e20,1.014e20,8.496e19,
                7.116e19,5.962e19,4.993e19,4.178e19,3.494e19,2.92e19,
                2.44e19,2.038e19,1.702e19,1.422e19,1.189e19,9.99e18,
                8.402e18,7.071e18, 5.956e18,5.021e18,4.237e18,3.578e18,
                3.023e18,2.552e18,2.144e18,1.8e18,1.524e18,1.301e18,
                1.118e18,9.681e17,8.43e17,7.382e17,
                6.498e17,5.748e17,5.107e17]
mssus76(::Type{Earth}) = [28.95,28.95,28.94,28.93,28.91,28.89,28.86,
              28.82,28.78,28.73,28.68,28.62,28.55,28.48,28.4,28.3,28.21,28.1,
              28.,27.88,27.77,27.64,27.52,27.39,27.27,27.14,27.02,26.9,26.79,
              26.68,26.58,26.48,26.38,26.29,26.2]
zsus76(::Type{Earth}) = [0.,11.019,20.063,32.162,47.35,51.413,71.802, 86.]
c_b5us76(::Type{Earth}) = 0.0
c_b6us76(::Type{Earth}) = 1e-10

z7us76(::Type{Earth}) = 86.
t7us76(::Type{Earth}) = 186.8673
z8us76(::Type{Earth}) = 91.
z9us76(::Type{Earth}) = 110.
t9us76(::Type{Earth}) = 240.
z10us76(::Type{Earth}) = 120.
t10us76(::Type{Earth}) = 360.
tinfus76(::Type{Earth}) = 1000.
tlk9us76(::Type{Earth}) = 12.
tcus76(::Type{Earth}) = 263.1905
taus76(::Type{Earth}) = -76.3232
saus76(::Type{Earth}) = -19.9426
rlamus76(::Type{Earth}) = 0.01875

import AstroDynBase: CelestialBody

function US76Atmos(Body::Type{<:CelestialBody}, PosRsph::Vector{Float64})
  zm = Altitude(Body, PosRsph)
  zKm = zm/1000;
  if (zKm < 86.)
    for i=1:7
      if (zsus76(Body)[i] <= zKm) && (zKm < zsus76(Body)[i + 1])
        zl = Int(round(zsus76(Body)[i]))
        zlm = zl*1000
        m = m0us76(Body);
        if zKm > 80.
          #Linear Interpolation on molecular weight
          iz = Int(round(((zKm - 80.)*2) + 1))
          zz1 = (iz - 1)*0.5 + 80.
          zz2 = zz1 + 0.5
          aux = (zKm - zz1)/(zz2 - zz1)
          m = m0us76(Body)*(m_m0us76(Body)[iz] + aux*(m_m0us76(Body)[iz + 1] - m_m0us76(Body)[iz]))
        end #if zKm > 80.
        ht = r0us76(Body)*zKm/(r0us76(Body) + zKm)
        hm = ht*1000

        #Pressure
        if abs(tgus76(Body)[i] - c_b5us76(Body) <= c_b6us76(Body))
          p = psus76(Body)[i]*exp(-(g0us76(Body)*m0us76(Body)*(hm - zlm))/(rsus76(Body)*tmsus76(Body)[i]))
        else
          d1 = tmsus76(Body)[i]/(tmsus76(Body)[i] + tgus76(Body)[i] * (ht - zl))
          d2 = g0us76(Body)*m0us76(Body)/(rsus76(Body)*tgus76(Body)[i]*0.001)
          p = psus76(Body)[i]*d1^d2
        end

        #molecular scale temperature
        tm = tmsus76(Body)[i] + tgus76(Body)[i]*(ht-zl)

        #Kinetic temperature
        tkin = tm*m/m0us76(Body)

        break;
      end #if (zs[i] <= zKm) && (zKm < zs[i + 1])
    end #for i=1:7
  elseif (zKm >= 86.) && (zKm < 120.)
    tkin = tempr_(zKm)

    #linear inteprolation on molecular weight
    iz = Int(round((zKm - 86) + 1))
    aux = (zKm - zssus76(Body)[iz])/(zssus76(Body)[iz+1] - zssus76(Body)[iz])
    m = mssus76(Body)[iz] + aux*(mssus76(Body)[iz + 1] - mssus76(Body)[iz])

    #linear interpolation on logarithm of totanl number Density
    aux = log(dnssus76(Body)[iz + 1]/dnssus76(Body)[iz])/(zssus76(Body)[iz+1] - zssus76(Body)[iz])
    d = dnssus76(Body)[iz] *exp(aux*(zKm - zssus76(Body)[iz]))

    #Molecular scale temperature
    tm = tkin*m0us76(Body)/m

    #Pressure
    p = d*rsus76(Body)*tkin/naus76(Body)

  end #if (zKm < 86.)
  #check whether altitude is in the range of 0 - 120 km
  if (zKm < 0) || (zKm > 120)
    p = 0.0
    p = 0.0
    tkin = 0.0
    a = 0.0
  else
    ρ = m*p/rsus76(Body)/tkin
    a = sqrt(γus76(Body)*rsus76(Body)*tkin/m0us76(Body))
  end #if
  return [ρ, p, tkin, a]
end #US76Atmos

function tempr_(Body::Type{<:CelestialBody},zKm::Float64)
  #= This function calculates the kinetic temperature for altitudes
   greater than 86 km, according to the US Standard Atmosphere 1976. =#
  if (zKm >= z7us76(Body)) && (zKm < z8us76(Body))
    return t = t7us76(Body)
  elseif (zKm >= z8us76(Body)) && (zKm < z9us76(Body))
    return t = tcus76(Body) + taus76(Body)*sqrt(1-(zKm - z8us76(Body))/saus76(Body)*((zKm - z8us76(Body))/saus76(Body)))
  elseif (zKm >= z9us76(Body)) && (zKm < z10us76(Body))
    return t = t9us76(Body) + tlk9us76(Body)*(zKm - z9us76(Body))
  elseif zKm >= z10us76(Body)
    zeta = (zKm - z10us76(Body))*(r0us76(Body) + z10us76(Body))/(r0us76(Body) + zKm)
    return t = tinfus76(Body) - (tinfus76(Body) - t10us76(Body))*exp(-rlamus76(Body)*zeta)
  end #if
end #tempr_


text:
module ch12
using Printf
# ------
# Tuples
# https://benlauwens.github.io/ThinkJulia.jl/latest/book.html#chap12
# ------

# --------------------
# Tuples Are Immutable
# https://benlauwens.github.io/ThinkJulia.jl/latest/book.html#_tuples_are_immutable
# --------------------
x = ('a', 'b')
typeof(x)

# single element tuple
x = ('a')
typeof(x) # Char
x = ('a', )
typeof(x) # Tuple{Char}

# tuple()
x = tuple(1, 'a', pi)
x[2]
x[2:3]

# compare
y = tuple(1, 'a', pi)

x == y
x > y

# ----------------
# Tuple Assignment
# ----------------
a, b = 1, 2, 3

a = 5
b = 7
a, b = b, a

# split
addr = "julius.caesar@rome"
uname, domain = split(addr, '@')

# -----------------------
# Tuples as Return Values
# -----------------------
x = divrem(7, 3)
q,r = divrem(7,3)
# -------------------------------
# Variable-length Argument Tuples
# -------------------------------

# gather (slurp)
function printall(args...)
    println(args)
end
printall(1, 2, '3')

# scatter (splat)
x = (7,3)
# divrem(x) error
divrem(x...)

# Exercise 12-1
function sumall(args...)
    sum = 0
    for value in args
        sum += value
    end

    return sum
end

sumall(1,2,3)

# -----------------
# Arrays and Tuples
# -----------------
str = "abc"
x = [1, 2, 3]
# zip
z = zip(str, x)
for p in z
    println(p)
end
#println(collect(z)) # [('a', 1), ('b', 2), ('c', 3)]

for (char, num) in z
    @printf("%g %s\n", num, char)
end

# enumerate
for (ind, char) in enumerate("abc")
    #println(typeof(ind))
    #println(typeof(char))
    @printf("%i %s\n", ind, char)
end

# -----------------------
# Dictionaries and Tuples
# -----------------------
dict = Dict('a'=>1, 'b'=>2, 'c'=>3)
for (key,value) in dict
    @printf("%s %i\n", key, value)
end

z = zip(str, x)
dict = Dict(z)
for (key,value) in dict
    @printf("%s %i\n", key, value)
end

dict = Dict()
dict['a', 1] = 3
dict['b', 1] = 4
for ((well, stage), num) in dict
    @printf("%s,%i: %i\n", well, stage, num)
end

# ----------------------
# Sequences of Sequences
# ----------------------
println(sort([13, 7, 5]))
println(reverse([13, 7, 5]))
end


text:
const INITIAL_READ_SIZE = 1024
const MAX_READ_SIZE     = 16*1024*1024

@assert MAX_READ_SIZE%INITIAL_READ_SIZE==0 && ispow2(div(MAX_READ_SIZE,INITIAL_READ_SIZE)) "MAX_READ_SIZE should be a power of two multiple of INITIAL_READ_SIZE"

struct Tag
	t::UInt16 # two-character representation
end
Tag(a, b) = Tag(UInt16(a)&0xff | (UInt16(b)<<8))
function Tag(s::String)
	@assert length(s)==2 "Length of tag string must be 2"
	Tag(s[1],s[2])
	#Tag(UInt16(s[1])&0xff | (UInt16(s[2])&0xff)<<8)
end

# allow writing tag"NH"
macro tag_str(p)
  Tag(p)
end

print(io::IO,tag::Tag) = print(io,Char(tag.t&0xff),Char(tag.t>>8))
show(io::IO,tag::Tag) = print(io,tag)

mutable struct BamRead
	file::BamFile # file containing the read

	# fields of known size
	block_size::Int32
	refID::Int32
	pos::Int32
	bin_mq_nl::UInt32
	flag_nc::UInt32
	l_seq::Int32
	next_refID::Int32
	next_pos::Int32
	tlen::Int32

	# remaining data of unknown size
	buf::Array{UInt8,1}

	# information to make indexing into buf quicker
	#cigar_start::Int
	seq_start::Int
	qual_start::Int
	#tags_start::Int

	# tags_parsed::Bool
	# tagDict::Dict{Tag,Int} # TAG->offset to start of tag (type byte) in buf
	tagsParsePos::Int # byte offset of next tag to parse
	tagInfo::Array{Tuple{Tag,Int},1}
end

#BamRead() = BamRead(0,0,0,0,0,0,0,0,0,Array{UInt8,1}(INITIAL_READ_SIZE),-1,-1)
#BamRead() = BamRead(0,0,0,0,0,0,0,0,0,Array{UInt8,1}(INITIAL_READ_SIZE),-1,-1,false,Dict{UInt16,Int}())
# BamRead(f::BamFile) = BamRead(f,0,0,0,0,0,0,0,0,0,Array{UInt8,1}(INITIAL_READ_SIZE),-1,-1,false,Dict{UInt16,Int}())
BamRead(f::BamFile) = BamRead(f,0,0,0,0,0,0,0,0,0,zeros(UInt8,INITIAL_READ_SIZE),-1,-1,0,Tuple{Tag,Int}[])

block_size(r::BamRead) = Int(r.block_size)

refID(r::BamRead) = Int(r.refID)
pos(r::BamRead) = Int(r.pos)+1
bin_mq_nl(r::BamRead) = UInt(r.bin_mq_nl)
flag_nc(r::BamRead) = UInt(r.flag_nc)
l_seq(r::BamRead) = Int(r.l_seq)
next_refID(r::BamRead) = Int(r.next_refID)
next_pos(r::BamRead) = Int(r.next_pos)+1
tlen(r::BamRead) = Int(r.tlen)

# access to things encoded in the fields
bin(r::BamRead)         = Int(r.bin_mq_nl>>16)
mapq(r::BamRead)        = Int((r.bin_mq_nl>>8)&0x000000ff)
l_read_name(r::BamRead) = Int(r.bin_mq_nl&0x000000ff)
flag(r::BamRead)        = Int(r.flag_nc>>16)
n_cigar_op(r::BamRead)  = Int(r.flag_nc&0x0000ffff)

flagset(r::BamRead,f::Int64) = (flag(r)&f) != 0

# access to things using BamFile info
ref_name(r::BamRead)      = sequence_name(r.file,refID(r))
next_ref_name(r::BamRead) = sequence_name(r.file,next_refID(r))

getuint32(r::BamRead, byteOffs::Int) = (UInt32(r.buf[byteOffs])) | (UInt32(r.buf[byteOffs+1])<<8) | (UInt32(r.buf[byteOffs+2])<<16) | (UInt32(r.buf[byteOffs+3])<<24)
getuint16(r::BamRead, byteOffs::Int) = (UInt16(r.buf[byteOffs])) | (UInt16(r.buf[byteOffs+1])<<8)
getuint8( r::BamRead, byteOffs::Int) = r.buf[byteOffs]

getint32(r::BamRead, byteOffs::Int) = reinterpret(Int32, getuint32(r,byteOffs))
getint16(r::BamRead, byteOffs::Int) = reinterpret(Int16, getuint16(r,byteOffs))
getint8( r::BamRead, byteOffs::Int) = reinterpret(Int8, getuint8(r,byteOffs))

function compute_offsets!(r::BamRead)
	#r.cigar_offs = l_read_name(r)
	#r.seq_offs   = r.cigar_offs + 4*n_cigar_op(r)
	r.seq_start   = l_read_name(r)+1 + 4*n_cigar_op(r)
	r.qual_start  = r.seq_start + ((l_seq(r)+1)>>1)
	#r.tags_start  = r.qual_start + l_seq(r)
	# r.tags_parsed = false
	# isempty(r.tagDict) || (r.tagDict = Dict{Tag,Int}()) # TODO: avoid creating a new dict for every read, this causes too many allocations.

	r.tagsParsePos = r.qual_start + l_seq(r)
	resize!(r.tagInfo,0)

	return
end

# overwrite read with data from stream
#function nextread!(io::IO,r::BamRead)
function nextread!(io,r::BamRead)
	r.block_size = read(io, Int32)
	r.refID      = read(io, Int32)
	r.pos        = read(io, Int32)
	r.bin_mq_nl  = read(io, UInt32)
	r.flag_nc    = read(io, UInt32)
	r.l_seq      = read(io, Int32)
	r.next_refID = read(io, Int32)
	r.next_pos   = read(io, Int32)
	r.tlen       = read(io, Int32)

	reqsize = Int(r.block_size)-8*4
	#reqsize > length(r.buf) && error("Read buffer too small!") # we could expand it, but if the file is corrupt, we could allocate huge things here...
	@assert reqsize>=0 "Corrupt Bam Read, block_size too small"

	# expand read buffer within reasonable limits
	# println("pointer: ", pointer(r.buf), ", size: ", length(r.buf))
	if reqsize>length(r.buf)
		sz = nextpow2(reqsize)
		@assert sz<=MAX_READ_SIZE "Read does not fit in maximum read buffer size ($MAX_READ_SIZE bytes)."
		r.buf = Array{UInt8,1}(sz)
	end

	# println("reqsize: ", reqsize)

	n = readbytes!(io,r.buf,reqsize)
	n!=reqsize && error("Unexpected end of stream when loading read.")
	compute_offsets!(r)
	r
end

# CIGAR and Seq needs to be (immutable) types rather than simple arrays so that we can interpret them for the user

# TODO: should CIGAR be a subtype of AbstractArray???
struct CIGAR{T<:AbstractArray{UInt8}}
	# TODO: Probably make the array type a type parameter.
	# This way, it also makes sure we don't need to wrap the array into a subarray when it's not necessary.
	# buf::SubArray{UInt8,1,Array{UInt8,1},Tuple{UnitRange{Int}},true}
	buf::T
end

# allow user to construct CIGAR from CIGAR string
function CIGAR(str::AbstractString)
	# len = 0 # compute length of CIGAR
	# for c in str
	# 	len += isalpha(c)
	# end

	# compute cigar length and preallocate array
	cigarLength = length(str) - mapreduce(isnumeric,+,str)
	buf = zeros(UInt32,cigarLength)

	# now interpret cigar string
	len = 0
	ops = Dict([('M',CIGAR_M), ('I',CIGAR_I), ('D',CIGAR_D), ('N',CIGAR_N), ('S',CIGAR_S), ('H',CIGAR_H), ('P',CIGAR_P), ('=',CIGAR_Eq), ('X',CIGAR_X)])

	i = 1
	for c in str
		if isnumeric(c)
			len = len*10 + (c-'0')
			continue
		end
		@assert len>0 "Malformed CIGAR string."
		buf[i] = (len<<4) | ops[c]
		len = 0
		i = i+1
	end
	CIGAR(view(reinterpret(UInt8,buf),1:cigarLength*4))
end

struct SeqElement
	e::Int
end

# TODO: should Seq be a subtype of AbstractArray???
struct Seq <: AbstractArray{SeqElement,1}
	#buf::SubArray{UInt8,1,Array{UInt8,1},Tuple{UnitRange{Int}},1}
	buf::Array{UInt8,1} # the entire buffer with the read (same as BamRead.buf)
	#start::Int          # offset into read buffer - TODO: change to half-bytes to allow substrings???
	start::Int         # offset in half-bytes into read buffer
	len::Int            # length in half-bytes!
end

# access to variable length fields
# read_name(r::BamRead) = SubString(String(r.buf), 1, l_read_name(r)-1)
read_name(r::BamRead) = String(r.buf[1:l_read_name(r)-1])
cigar(r::BamRead)     = (s=l_read_name(r)+1; CIGAR(view(r.buf, s:s+n_cigar_op(r)*4-1)))
seq(r::BamRead)       = Seq(r.buf, (r.seq_start-1)*2, l_seq(r))
qual(r::BamRead)      = view(r.buf, r.qual_start:r.qual_start+l_seq(r)-1)

# TODO: make sure tags are properly copied if we allow the read to not start at the beginning of buf
function copy(r::BamRead)
	BamRead(r.file, r.block_size, r.refID, r.pos, r.bin_mq_nl, r.flag_nc,
			r.l_seq, r.next_refID, r.next_pos, r.tlen,
			copy(view(r.buf,1:r.block_size-8*4)),
			r.seq_start, r.qual_start, r.tagsParsePos, copy(r.tagInfo))
end

copy(c::CIGAR) = CIGAR(view(copy(c.buf),1:length(c.buf))) # copy to array and make SubArray covering the whole
#copy(s::Seq)   = Seq(copy(view(s.buffer,s.start:s.start+((s.len+1)>>1))), 1, s.len)
copy(s::Seq)   = Seq(copy(view(s.buf,s.start>>1+1:(s.start+s.len+1)>>1)), mod(s.start,2), s.len)

struct CIGARElement
	e::UInt32
end
CIGARElement(c::CIGAR,i::Int) = CIGARElement((UInt32(c.buf[i])) | (UInt32(c.buf[i+1])<<8) | (UInt32(c.buf[i+2])<<16) | (UInt32(c.buf[i+3])<<24))

op(e::CIGARElement) = Int(e.e&0xf)
const op_chars = Char['M','I','D','N','S','H','P','=','X']
op_char(e::CIGARElement) = op_chars[op(e)+1]
len(e::CIGARElement) = Int(e.e>>4)

print(io::IO,e::CIGARElement) = print(io,len(e),op_char(e))

# CIGAR iterators
# Iterators
# start(::CIGAR) = 1
# next(c::CIGAR, i::Int) = (CIGARElement(c,i), i+4)
# done(c::CIGAR, i::Int) = i > length(c.buf)

function iterate(c::CIGAR, i::Int=1)
	i>length(c.buf) && return nothing
	(CIGARElement(c,i), i+4)
end
eltype(::Type{T}) where {T<:CIGAR} = CIGARElement
length(c::CIGAR) = div(length(c.buf),4)

function alignedlength(cig::CIGAR)
	l = 0
	for c in cig
		o = op(c)
		if o==CIGAR_M || o==CIGAR_X || o==CIGAR_Eq || o==CIGAR_D || o==CIGAR_N
			l += len(c)
		end
	end
	l
end

endpos(r::BamRead) = pos(r) + alignedlength(cigar(r)) - 1

function SeqElement(s::Seq, i::Int)
	b = Int(s.buf[(i+1)>>1]) # get the correct byte
	b = i&1==0 ? b&0x0f : b>>4 # select low or high nybble depending on index
	SeqElement(b)
end

const base_indices = [0,1,2,0,3,0,0,0,4,0,0,0,0,0,0,0]
const base_chars = UInt8['=','A','C','M','G','R','S','V','T','W','Y','H','K','D','B','N'];
base_mask(e::SeqElement) = e.e
base_index(e::SeqElement) = base_indices[e.e+1]
base_char(e::SeqElement) = base_chars[e.e+1] # TODO: should this be renamed, return Char or change in some other way?

# Iterators
# start(::Seq) = 1
# next(s::Seq, i::Int) = (SeqElement(s,i+s.start), i+1)
# done(s::Seq, i::Int) = i > s.len
function iterate(s::Seq, i::Int=1)
	i > s.len && return nothing
	(SeqElement(s,i+s.start), i+1)
end
eltype(::Type{Seq}) = SeqElement

size(s::Seq) = (s.len,)
length(s::Seq) = s.len
# linearindexing(::Type{Seq}) = Base.LinearFast()
Base.IndexStyle(::Type{Seq}) = Base.IndexStyle(Type{Array})

# TODO: is this the proper way to do bounds checking???
checkbounds(s::Seq, i::Int) = 1<=i<=s.len || throw(BoundsError("$(s.len)-element Seq",i))
checkbounds(::Type{Bool}, s::Seq, i::Int) = 1<=i<=s.len

# random access
#getindex(s::Seq, i::Int) = (checkbounds(s,i); SeqElement(s,i+(s.start-1)*2))
getindex(s::Seq, i::Int) = (checkbounds(s,i); SeqElement(s,i+s.start))

# construct view into seq
#view(s::Seq, rng::UnitRange{Int}) = (checkbounds(s,rng[1]); checkbounds(s,rng[end]); Seq(s.buf,s.start+rng[1]-1,length(rng)))
function view(s::Seq, rng::UnitRange{Int})
	length(rng)==0 && return Seq(s.buf,0,0) # anything with length 0 will do
	checkbounds(s,rng[1])
	checkbounds(s,rng[end])
	Seq(s.buf,s.start+rng[1]-1,length(rng))
end

function print(io::IO,cig::CIGAR)
	for c in cig
		print(io,c)
	end
end
show(io::IO,cig::CIGAR) = print(io,cig)

# TODO: implement this in a better (still typestable) way.
function string(cig::CIGAR)
	io = IOBuffer()
	print(io,cig)
	String(take!(io)) # for type stability
end
cigar_str(r::BamRead) = string(cigar(r)) # convenience access

function print(io::IO,seq::Seq)
	for (i,n) in enumerate(seq)
		print(io,Char(base_char(n)))
	end
end
show(io::IO,seq::Seq) = print(io,seq)
function string(seq::Seq)
	str = zeros(UInt8,seq.len) # allocate buffer
	for (i,n) in enumerate(seq)
		str[i] = base_char(n)
	end
	String(str)
end
seq_str(r::BamRead) = string(seq(r)) # convenience access

# tags - array implementation

read_tag_name(r::BamRead, offs::Int) = Tag(r.buf[offs], r.buf[offs+1])
read_tag_type(r::BamRead, offs::Int) = Char(r.buf[offs])

# sizes for fixed-size tags
const tagTypeSizes = Dict{Char,Int}('A'=>1,'c'=>1,'C'=>1,'s'=>2,'S'=>2,'i'=>4,'I'=>4,'f'=>4)

function parsenexttag!(r::BamRead)
	tag = read_tag_name(r,r.tagsParsePos) # two-character tag
	r.tagsParsePos += 2

	# store offset in array
	push!(r.tagInfo, (tag, r.tagsParsePos))

	# now skip rest of tag
	tagType = read_tag_type(r,r.tagsParsePos)
	r.tagsParsePos += 1

	if tagType=='Z'
		# iterate until 0 character
		while r.buf[r.tagsParsePos] != 0x0
			r.tagsParsePos += 1
		end
		r.tagsParsePos += 1
	elseif tagType=='B'
		error("Tags with array type (B) are not yet supported.")
	elseif tagType=='H' # is this tag possible in BAM???
		error("Tags with type H are not yet supported.")
	else
		# fixed-size tag
		r.tagsParsePos += tagTypeSizes[tagType]
	end
	tag
end

# find tag and parse offset until it is found
function findtag!(r::BamRead, tag::Tag)

	# first check in already parsed tags
	for (currTag, offs) in r.tagInfo
		currTag == tag && return offs
	end

	# and if it wasn't found, continue looking in unparsed part
	while r.tagsParsePos < r.block_size-8*4 # used part of buf
		currTag = parsenexttag!(r)
		currTag == tag && return r.tagInfo[end][2]
	end

	-1
end

hastag(r::BamRead, tag::Tag) = findtag!(r,tag) != -1
hastag(r::BamRead, tag::String) = hastag(r,Tag(tag))

# TODO: support default value if tag doesn't exist?
function tagint(r::BamRead, tag::Tag)
	offs = findtag!(r,tag)
	tagType = read_tag_type(r,offs)
	offs += 1

	#println("tagtype: $tagType")

	# get value from byte buffer
	if tagType=='i'
		v = Int(getint32(r,offs))
	elseif tagType=='I'
		v = Int(getuint32(r,offs))
	elseif tagType=='s'
		v = Int(getint16(r,offs))
	elseif tagType=='S'
		v = Int(getuint16(r,offs))
	elseif tagType=='c'
		v = Int(getint8(r,offs))
	elseif tagType=='C'
		v = Int(getuint8(r,offs))
	else
		error("Cannot convert tag type \"$tagType\" to Int")
	end
	v
end
tagint(r::BamRead, tag::String) = tagint(r,Tag(tag))

# TODO: support default value if tag doesn't exist?
function tagchar(r::BamRead, tag::Tag)
	offs = findtag!(r,tag)
	tagType = read_tag_type(r,offs)
	offs += 1

	@assert tagType=='A' "Cannot convert tag type \"$tagType\" to Char"
	Char(getuint8(r,offs))
end
tagchar(r::BamRead, tag::String) = tagchar(r,Tag(tag))

# TODO: support default value if tag doesn't exist?
function tagfloat32(r::BamRead, tag::Tag)
	offs = findtag!(r,tag)
	tagType = read_tag_type(r,offs)
	offs += 1

	@assert tagType=='f' "Cannot convert tag type \"$tagType\" to Float64"
	reinterpret(Float32,getuint32(r,offs))
end
tagfloat32(r::BamRead, tag::String) = tagfloat(r,Tag(tag))

tagfloat(r::BamRead, tag) = Float(tagfloat(r,tag))

# TODO: support default value if tag doesn't exist?
function tagstring(r::BamRead, tag::Tag)
	offs = findtag!(r,tag)
	tagType = read_tag_type(r,offs)
	offs += 1

	@assert tagType=='Z' "Cannot convert tag type \"$tagType\" to string"

	start = offs

	# find end of string - TODO: no need to do this twice, store end pos in dict somehow?
	while r.buf[offs] != 0x0
		offs += 1
	end
	offs -= 1 # do not include NULL terminator

	# SubString(String(r.buf),start,offs)
	String(r.buf[start:offs])
end
tagstring(r::BamRead, tag::String) = tagstring(r, Tag(tag))

function printtag(io::IO, r::BamRead, tag::Tag, offs::Int)
	tagType = read_tag_type(r,offs)
	if tagType=='i' || tagType=='s' || tagType=='c' || tagType=='I' || tagType=='S' || tagType=='C'
		print(io, '\t', tag, ":i:", tagint(r,tag))
	elseif tagType=='A'
		print(io, '\t', tag, ":A:", tagchar(r,tag))
	elseif tagType=='f'
		print(io, '\t', tag, ":f:", tagfloat(r,tag)) # TODO: check exactly how floats are normally printed to SAM
	elseif tagType=='Z'
		print(io, '\t', tag, ":Z:", tagstring(r,tag))
	else
		error("Tags with array type (", tagType, ") are not currently supported.")
	end
end

function printtags(io::IO, r::BamRead)
	# first loop over already parsed tags
	for (tag, offs) in r.tagInfo
		printtag(io, r, tag, offs)
	end

	# and continue with previously unparsed tags
	while r.tagsParsePos < r.block_size-8*4 # used part of buf
		tag = parsenexttag!(r)
		printtag(io, r, tag, r.tagInfo[end][2])
	end
end

# # tags - dict implementation

# #read_tag_name(r::BamRead, offs::Int) = UInt16(r.buf[offs]) | UInt16(r.buf[offs+1])<<8
# read_tag_name(r::BamRead, offs::Int) = Tag(r.buf[offs], r.buf[offs+1])
# read_tag_type(r::BamRead, offs::Int) = Char(r.buf[offs])

# # sizes for fixed-size tags
# const tagTypeSizes = Dict{Char,Int}('A'=>1,'c'=>1,'C'=>1,'s'=>2,'S'=>2,'i'=>4,'I'=>4,'f'=>4)

# function parsetags(r::BamRead)
# 	offs = r.qual_start + l_seq(r)

# 	while offs < r.block_size-8*4 # used part of buf
# 		tag = read_tag_name(r,offs) # two-character tag
# 		offs += 2

# 		# store offset in dict
# 		r.tagDict[tag] = offs

# 		# now skip rest of tag
# 		tagType = read_tag_type(r,offs)
# 		offs += 1

# 		if tagType=='Z'
# 			#error("Tags with string type are not yet supported.")
# 			# println(SubString(String(r.buf),offs,offs+10-1))
# 			# println(r.buf[offs+10])

# 			# iterate until 0 character
# 			while r.buf[offs] != 0x0
# 				offs += 1
# 			end
# 			offs += 1
# 		elseif tagType=='B'
# 			error("Tags with array type (B) are not yet supported.")
# 		elseif tagType=='H' # is this tag possible in BAM???
# 			error("Tags with type H are not yet supported.")
# 		else
# 			# fixed-size tag
# 			offs += tagTypeSizes[tagType]
# 		end
# 	end

# 	r.tags_parsed = true
# 	return
# end

# function hastag(r::BamRead, tag::Tag)
# 	r.tags_parsed || parsetags(r)
# 	#tagname2uint16(tag) ∈ keys(r.tagDict)
# 	haskey(r.tagDict,tag)
# end
# hastag(r::BamRead, tag::String) = hastag(r,Tag(tag))

# # TODO: support default value if tag doesn't exist?
# function tagint(r::BamRead, tag::Tag)
# 	r.tags_parsed || parsetags(r)
# 	offs = r.tagDict[tag]
# 	tagType = read_tag_type(r,offs)
# 	offs += 1

# 	#println("tagtype: $tagType")

# 	# get value from byte buffer
# 	if tagType=='i'
# 		v = Int(getint32(r,offs))
# 	elseif tagType=='I'
# 		v = Int(getuint32(r,offs))
# 	elseif tagType=='s'
# 		v = Int(getint16(r,offs))
# 	elseif tagType=='S'
# 		v = Int(getuint16(r,offs))
# 	elseif tagType=='c'
# 		v = Int(getint8(r,offs))
# 	elseif tagType=='C'
# 		v = Int(getuint8(r,offs))
# 	else
# 		error("Cannot convert tag type \"$tagType\" to Int")
# 	end
# 	v
# end
# tagint(r::BamRead, tag::String) = tagint(r,Tag(tag))

# # TODO: support default value if tag doesn't exist?
# function tagchar(r::BamRead, tag::Tag)
# 	r.tags_parsed || parsetags(r)
# 	offs = r.tagDict[tag]
# 	tagType = read_tag_type(r,offs)
# 	offs += 1

# 	@assert tagType=='A' "Cannot convert tag type \"$tagType\" to Char"

# 	Char(getuint8(r,offs))
# end
# tagchar(r::BamRead, tag::String) = tagchar(r,Tag(tag))

# # TODO: support default value if tag doesn't exist?
# function tagstring(r::BamRead, tag::Tag)
# 	r.tags_parsed || parsetags(r)
# 	offs = r.tagDict[tag]
# 	tagType = read_tag_type(r,offs)
# 	offs += 1

# 	@assert tagType=='Z' "Cannot convert tag type \"$tagType\" to string"

# 	start = offs

# 	# find end of string - TODO: no need to do this twice, store end pos in dict somehow?
# 	while r.buf[offs] != 0x0
# 		offs += 1
# 	end
# 	offs -= 1 # do not include NULL terminator

# 	SubString(String(r.buf),start,offs)
# end
# tagstring(r::BamRead, tag::String) = tagstring(r, Tag(tag))

# create SAM string representation of BamRead
function print(io::IO, r::BamRead)
	# QNAME
	# FLAG
	# RNAME
	# POS
	# MAPQ
	# CIGAR
	# RNEXT
	# PNEXT
	# TLEN
	# SEQ
	# QUAL
	print(io, read_name(r), '\t')
	print(io, flag(r), '\t')
	print(io, ref_name(r), '\t')
	print(io, pos(r), '\t')
	print(io, mapq(r), '\t')
	cig = string(cigar(r))
	print(io, length(cig)==0 ? '*' : cig, '\t')
	#print(io, next_ref_name(r), '\t')
	print(io, next_refID(r)==refID(r) ? "=" : next_ref_name(r), '\t')
	print(io, next_pos(r), '\t')
	print(io, tlen(r), '\t')
	print(io, seq(r), '\t')
	print(io, String(qual(r).+0x21)) # NO tab, will be added by tag printing if necessary

	printtags(io, r)
end
show(io::IO, r::BamRead) = print(io,r)

function sam(r::BamRead)
	io = IOBuffer()
	print(io,r)
	String(take!(io))
end


text:
module HelloCExecutable

using HelloCExecutable_jll
export hello_world

function hello_world()
    cmd = HelloCExecutable_jll.hello_world(identity)
    msg = `$cmd` |> pipeline |> read |> String # "Hello, Goma!\n"
    # Remove a newline from `msg`
    msg = msg |> chomp # Should be `"Hello, Goma!"`
    return msg
end

end


text:
include("../emi.jl")
using EMI
using EMI.Draw, EMI.Gmsh

# Something like this with circular
#   ^
#  / \
# [   ]
#  \ /
#   U

const R, dx, dy = 1, 0.2, 0.3
loop = Loop([Line(Point(0, R+dy+dx), Point(0, R+dy)),
             Line(Point(0, R+dy), Point(dy, R+dy)),
             CircleArc(Point(dy, R+dy), Point(R+dy, R+dy), Point(R+dy, dy)),
             Line(Point(R+dy, dy), Point(R+dy, 0)),
             Line(Point(R+dy, 0), Point(R+dy+dx, 0)),
             Line(Point(R+dy+dx, 0), Point(R+dy+dx, dy)),
             CircleArc(Point(R+dy+dx, dy), Point(R+dy+dx, R+dy), Point(2R+dy+dx, R+dy)),
             Line(Point(2R+dy+dx, R+dy), Point(2R+2*dy+dx, R+dy)),
             Line(Point(2R+2*dy+dx, R+dy), Point(2R+2*dy+dx, R+dy+dx)),
             Line(Point(2R+2*dy+dx, R+dy+dx), Point(2R+dy+dx, R+dy+dx)),
             CircleArc(Point(2R+dy+dx, R+dy+dx), Point(R+dy+dx, R+dy+dx), Point(R+dy+dx, 2R+dy+dx)),
             Line(Point(R+dy+dx, 2R+dy+dx), Point(R+dy+dx, 2R+2*dy+dx)),
             Line(Point(R+dy+dx, 2R+2*dy+dx), Point(R+dy, 2R+2*dy+dx)),
             Line(Point(R+dy, 2R+2*dy+dx), Point(R+dy, 2R+dy+dx)),
             CircleArc(Point(R+dy, 2R+dy+dx), Point(R+dy, R+dy+dx), Point(dy, R+dy+dx)),
             Line(Point(dy, R+dy+dx), Point(0, R+dy+dx))])

canvas = Canvas()
canvas = canvas + loop

set_bbox!(canvas, 0.2, 0.3)
gmsh_script(canvas, 0.2)

max_stars_repo_name: MiroK/emi-cylinders

max_issues_repo_name: MiroK/emi-cylinders

max_forks_repo_name: MiroK/emi-cylinders


text:
function mm_p(t::SimplexSplitting.Triangulation)

    const n_simplices = size(t.simplex_inds, 1)
    const dim = size(t.points, 2)

    # Mismatch for the markovity of the Markov matrix is at most delta if volume_tolerance
    # tolerance as below
    const delta::Float64 = 1/10^5
    const voltol::Float64 = delta/n_simplices

    # Tolerance for similary of convex expansion coefficients of simplex vertices in simplexintersection function.
    const convex_params_tol::Float64 = 1/10^12

    intvols = SharedArray{Float64}(n_simplices, n_simplices) # intersecting volumes

    @sync @parallel for i in 1:n_simplices
        imvol = t.volumes_im[i]
        #println("\tImage #", i, "/", n_simplices)
        for j in 1:n_simplices
            vol = t.volumes[j]
            if vol * imvol > 0 && (vol/imvol) > voltol
                intvol = Simplices.simplexintersection(t.points[t.simplex_inds[j, :], :].', t.impoints[t.simplex_inds[i, :], :].') / imvol
                intvols[i, j] = intvol
            end
        end
    end

    return Array(intvols)
end

function mm_ij()
end


text:
include("./imports.jl")

include("crop.jl")
include("utils.jl")


text:
function aggregate(x::AbstractArray{T, 1}, by::AbstractArray{T2, 1}, FUN::Function = mean) where {
    T<:Real, T2<:Real }

    grps = unique(by)
    # vals = ones(length(grps));
    vals = []
    @views for grp in grps
        val = FUN(x[by .== grp])
        push!(vals, val)
    end
    Dict(grps .=> vals)
end


text:
const O₂ₐ = 0.209L*L^-1 # Volume of O₂ in the air
const Dₗᵢ = 3.17 # Diffusion coefficient of substrate in liquid phase, dimensionless
const Dₒₐ = 1.67 # Diffusion coefficient of oxygen in air, dimensionless
const pₛₓ = 0.024 # Fraction of soil carbon that is considered soluble
const Eₛ = 10cm # Effective soil septh (assume only 10 first cm produce CO2)


text:
using KNITRO, JuMP, Base.Test

m = Model(solver=KnitroSolver())
@variable(m, x)
setValue(x, 5)
@NLobjective(m, Min, -log(x)+x^2)
solve(m)
@test_approx_eq_eps getobjectivevalue(m) 0.84657 1e-5
@test_approx_eq_eps getvalue(x) 0.707107 1e-5

ktrmod = internalmodel(m)
MathProgBase.freemodel!(ktrmod)


text:
using Test, Spirograph

spirograph(10, 3, 5)
@test true


text:
using Neurodevo
using Random

include("cgp.jl")
include("darwin.jl")
include("datasets.jl")

function classify(m::Model, X::Array{Float64}, Y::Array{Int64})
    fit = 0.0
    total_time = 0
    nout = length(unique(Y))
    for epoch in 1:m.cfg["n_epochs"]
        labels = zeros(Int64, length(Y))
        for i in eachindex(Y)
            outputs = zeros(nout)
            for s in 1:m.cfg["n_step"]
                outputs, t, bytes, gctime, memallocs = @timed Neurodevo.step!(m, X[:, i])
                total_time += t
            end
            if total_time >= m.cfg["time_max"]
                return [fit - 1.0]
            end
            labels[i] = argmax(outputs)
            if labels[i] == Y[i]
                Neurodevo.reward!(m, [1.0])
            end
        end
        fit = sum(labels .== Y) / length(Y)
    end
    [fit]
end

function data_evaluation(ind::NeurodevoInd)
    cfg = cgp_cfg(Config("cfg/evo.yaml"), ind.genes[1])
    c = cgp_controller(cfg, ind.genes[2:end]; cinds=[4, 5, 9, 10])
    m = Model(cfg, c)

    X, Y = ind.func(ind.seed)
    nin = size(X, 1)
    nout = length(unique(Y))
    layered_init!(m, nin, nout; nreward=1, nhidden=cfg["n_hidden"])
    classify(m, X, Y)
end

function NeurodevoInd(genes::Array{Array{Float64}}, fitness::Array{Float64})
    NeurodevoInd(genes, fitness, get_iris, 1)
end

function generation(e::Evolution)
    fits = [get_iris, get_diabetes, get_glass]
    for i in e.population
        i.seed = floor(Int64, (e.gen - 1) / 10)
        func = mod(floor(Int64, (e.gen - 1) / e.cfg["goal_gen"]), length(fits)) + 1
        i.func = fits[func]
    end
end


text:
# specify an Exasim version to run
version = "Version0.1";

# External packages
using Revise, DelimitedFiles, SymPy

# Add Exasim to Julia search path
cdir = pwd(); ii = findlast("Exasim", cdir);
include(cdir[1:ii[end]] * "/Installation/setpath.jl");

# Exasim packages
using Preprocessing, Mesh, Gencode, Postprocessing

# create pde structure and mesh structure
pde, mesh = Preprocessing.initializeexasim(version);

# Define PDE model: governing equations, initial solutions, and boundary conditions
pde.model = "ModelW";            # ModelC, ModelD, ModelW
include("pdemodel.jl");          # include the PDE model file

# Set discretization parameters, physical parameters, and solver parameters
pde.porder = 3;                  # polynomial degree
pde.torder = 3;                  # time-stepping order of accuracy
pde.nstage = 2;                  # time-stepping number of stages
pde.physicsparam = [1.0 10.0 0.0 1.0]; # relative permitivity, wave number vector, wave speed
pde.tau = [1.0];                 # DG stabilization parameter
pde.dt = 0.1*ones(100);   # time step sizes
pde.soltime = collect(1:length(pde.dt)); # steps at which solution are collected
pde.visdt = 0.1; # visualization timestep size
pde.GMRESrestart=40;
pde.NLiter=3;

# Choose computing platform and set number of processors
#pde.platform = "gpu";           # choose this option if NVIDIA GPUs are available
pde.mpiprocs = 4;                # number of MPI processors

# read a grid from a file
mesh.p, mesh.t = Mesh.readmesh("grid.bin",0);
# expressions for disjoint boundaries
mesh.boundaryexpr = [p -> (p[2,:] .< -12.0+1e-3), p -> (p[1,:] .> 12.0-1e-3), p -> (p[2,:] .> 12.0-1e-3), p -> (p[1,:] .< -12.0+1e-3), p -> (p[1,:] .< 20)];
mesh.boundarycondition = [1 1 1 1 2]; # Set boundary condition for each disjoint boundary
# expressions for curved boundaries
mesh.curvedboundary = [0 0 0 0 1];
mesh.curvedboundaryexpr = [p -> 0, p -> 0, p -> 0, p -> 0, p -> sqrt.(p[1,:].*p[1,:]+p[2,:].*p[2,:]).-1.0];

# call exasim to generate and run C++ code to solve the PDE model
sol, pde, mesh,~,~,~,~  = Postprocessing.exasim(pde,mesh);

# visualize the numerical solution of the PDE model using Paraview
pde.visscalars = ["velocity", 1, "displacement", 4];  # list of scalar fields for visualization
pde.visvectors = ["displacement gradient", [2, 3]]; # list of vector fields for visualization
Postprocessing.vis(sol,pde,mesh); # visualize the numerical solution
print("Done!");

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim