text:
export iauPas
"""
Position-angle from spherical coordinates.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  vector/matrix support function.

Given:
   al     double     longitude of point A (e.g. RA) in radians
   ap     double     latitude of point A (e.g. Dec) in radians
   bl     double     longitude of point B
   bp     double     latitude of point B

Returned (function value):
         double     position angle of B with respect to A

Notes:

   1. The result is the bearing (position angle), in radians, of point
   B with respect to point A.  It is in the range -pi to +pi.  The
   sense is such that if B is a small distance "east" of point A,
   the bearing is approximately +pi/2.

   2. Zero is returned if the two points are coincident.

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# double iauPas(double al, double ap, double bl, double bp)

function iauPas(al::Real, ap::Real, bl::Real, bp::Real)
   return ccall((:iauPas, libsofa_c), Cdouble,
       (Cdouble, Cdouble, Cdouble, Cdouble),
       convert(Float64, al),
       convert(Float64, ap),
       convert(Float64, bl),
       convert(Float64, bp))
end


text:
println("Hellow world!")
2

ext: jl

max_stars_repo_name: Timothee-ANNE/testGit

max_issues_repo_name: Timothee-ANNE/testGit

max_forks_repo_name: Timothee-ANNE/testGit


text:
"""
    parse_and_number(treestring::S)::GeneralNode where S<:AbstractString

--- INTERNAL ---
Parse a newick string and then call set_binary and number_nodes! on the
resulting tree
"""
function parse_and_number(treestring::S)::GeneralNode where S<:AbstractString
    p_tree2 = parsing_newick_string(string(treestring))
    set_binary!(p_tree2)
    MCPhyloTree.number_nodes!(p_tree2)
    p_tree2
end # parse_and_number

"""
    ASDSF(args::String...; freq::Int64=1, check_leaves::Bool=true,
          min_splits::Float64=0.1, show_progress::Bool=true)::Vector{Float64}

Calculate the average standard deviation of split frequencies for two or more
files containing newick representations of trees. Default frequency is 1 and
by default only trees with the same leafsets are supported. The default minimal
splits threshold is 0.1. The progress bar is activated by default.
"""
function ASDSF(args::String...; freq::Int64=1, check_leaves::Bool=true,
               min_splits::Float64=0.1, show_progress::Bool=true
               )::Vector{Float64}

    length(args) < 2 && throw(ArgumentError("At least two input files are needed."))
    splitsQueue = [Accumulator{Tuple{String, String}, Int64}()]
    splitsQueues = [Vector{Accumulator{Tuple{String, String}, Int64}}()]
    for arg in args
        push!(splitsQueues[1], Accumulator{Tuple{Set{String}, Set{String}}, Int64}())
    end # for
    iter = zip([eachline(arg) for arg in args]...)
    nsams = countlines(args[1]) / freq
    asdsf_size = ceil(Int,nsams)
    ASDSF_vals = [zeros(ceil(Int,asdsf_size))]

    # tree_dims is hardcoded to 1 because there are no multiple dims in files
    res = ASDSF_int(splitsQueue, splitsQueues, iter, 1, ASDSF_vals, freq, check_leaves,
              min_splits, show_progress, basic=true)[1][1]
    asdsf_size > nsams ? res[1:end-1] : res
end # ASDSF

"""
    ASDSF(args::Vector{String}...; freq::Int64=1, check_leaves::Bool=true,
          min_splits::Float64=0.1, show_progress::Bool=true)::Vector{Float64}

Calculate the average standard deviation of split frequencies for two or more
Vectors containing newick representations of trees. Default frequency is 1 and
by default only trees with the same leafsets are supported. The default minimal
splits threshold is 0.1. The progress bar is activated by default.
"""
function ASDSF(args::Vector{String}...; freq::Int64=1, check_leaves::Bool=true,
               min_splits::Float64=0.1, show_progress::Bool=true
               )::Vector{Float64}

    length(args) < 2 && throw(ArgumentError("At least two input arrays are needed."))
    splitsQueue = [Accumulator{Tuple{String, String}, Int64}()]
    splitsQueues = [Vector{Accumulator{Tuple{String, String}, Int64}}()]
    for arg in args
        push!(splitsQueues[1], Accumulator{Tuple{Set{String}, Set{String}}, Int64}())
    end # for
    iter = zip(args...)
    ASDSF_vals = [zeros(Int(length(iter) / freq))]
    # tree_dims is hardcoded to 1 because there are no multiple dims in Vectors
    ASDSF_int(splitsQueue, splitsQueues, iter, 1, ASDSF_vals, freq, check_leaves,
              min_splits, show_progress; basic=true)[1]
end # ASDSF

"""
    ASDSF(model::ModelChains; freq::Int64=1, check_leaves::Bool=true,
          min_splits::Float64=0.1, show_progress::Bool=true
          )::Vector{Vector{Float64}}

Calculate the average standard deviation of split frequencies for the trees in
different chains in a ModelChains object. Default frequency is 1 and by default
only trees with the same leafsets are supported. The default minimal splits
threshold is 0.1. The progress bar is activated by default.
"""
function ASDSF(model::ModelChains; freq::Int64=1, check_leaves::Bool=true,
               min_splits::Float64=0.1, show_progress::Bool=true
               )::Vector{Vector{Float64}}

    tree_dims::UnitRange{Int64} = 1:size(model.trees, 2)
    splitsQueue = [Accumulator{Tuple{String, String}, Int64}() for x in tree_dims]
    splitsQueues = [Vector{Accumulator{Tuple{String, String}, Int64}}() for x in tree_dims]
    nchains = size(model.trees, 3)
    for i in 1:nchains
        for j in tree_dims
            push!(splitsQueues[j], Accumulator{Tuple{Set{String}, Set{String}}, Int64}())
        end # for
    end # for
    trees = Array{Vector{AbstractString}, 2}(undef, size(model.trees, 1), nchains)
    if length(tree_dims) > 1
        trees = Array{Vector{AbstractString}, 2}(undef, size(model.trees, 1), nchains)
        for i in 1:size(model.trees, 1)
            for j in 1:nchains
                trees[i, j] = model.trees[i,:,j]
            end # for
        end # for
    end # if
    iter = zip([trees[:,c] for c in 1:nchains]...)
    ASDSF_vals::Vector{Vector{Float64}} = [zeros(Int(floor(length(iter) / freq))) for x in tree_dims]
    ASDSF_int(splitsQueue, splitsQueues, iter, tree_dims, ASDSF_vals, freq,
              check_leaves, min_splits, show_progress)
end # ASDSF

"""
    ASDSF(r_channels::Vector{RemoteChannel}, n_trees::Int64,
          tree_dims::UnitRange{Int64}, min_splits::Float64
          )::Tuple{Vector{Vector{Float64}}, ConvergenceStorage}

--- INTERNAL ---
Calculates - on-the-fly - the average standard deviation of split frequencies
for the trees generated by MCMC draws from a model. Takes a vector of remote
channels (where the generated trees are stored during the mcmc simulation) and
the total number of trees in each chain as arguments. The default minimal splits
threshold is 0.1.
"""
function ASDSF(r_channels::Vector{RemoteChannel{Channel{Array{AbstractString,1}}}},
               n_trees::Int64, tree_dims::UnitRange{Int64}, min_splits::Float64;
               cs::Union{Nothing, ConvergenceStorage}=nothing
               )::Tuple{Vector{Vector{Float64}}, ConvergenceStorage}

    iter = 1:n_trees
    ASDSF_vals::Vector{Vector{Float64}} = [zeros(Int(n_trees)) for x in tree_dims]
    nchains = length(r_channels)
    if isnothing(cs)
        splitsQueue = [Accumulator{Tuple{String, String}, Int64}() for x in tree_dims]
        splitsQueues = [Vector{Accumulator{Tuple{String, String}, Int64}}() for x in tree_dims]
        for i in 1:nchains
            for j in tree_dims
                push!(splitsQueues[j], Accumulator{Tuple{Set{String}, Set{String}}, Int64}())
            end # for
        end # for
        ASDSF_int(splitsQueue, splitsQueues, iter, tree_dims, ASDSF_vals, 1, false,
                  min_splits, false; r_channels=r_channels)
    else
        splitsQueue = cs.splitsQueue
        splitsQueues = cs.splitsQueues
        run = cs.run
        ASDSF_int(splitsQueue, splitsQueues, iter, tree_dims, ASDSF_vals, 1,
                  false, min_splits, false; r_channels=r_channels,
                  total_runs=run)
    end # if/else
end # ASDSF

"""
    ASDSF_int(splitsQueue, splitsQueues, iter, tree_dims, ASDSF_vals, freq,
              check_leaves, min_splits, show_progress; r_channels=nothing,
              run::Int64=1, basic=false
              )::Tuple{Vector{Vector{Float64}}, ConvergenceStorage}

--- INTERNAL ---
Handles the computation of the Average Standard Deviation of Split Frequencies.
"""
function ASDSF_int(splitsQueue, splitsQueues, iter, tree_dims, ASDSF_vals, freq,
                   check_leaves, min_splits, show_progress; r_channels=nothing,
                   total_runs::Int64=1, basic=false
                   )::Tuple{Vector{Vector{Float64}}, ConvergenceStorage}

    all_keys = [Set{Tuple{String,String}}() for x in tree_dims]

    if show_progress
        prog = ProgressMeter.Progress(length(ASDSF_vals[1]),"Computing ASDSF: ")
    end # if
    gen = 1

    for (i, line) in enumerate(iter)
        if mod(i, freq) == 0
            if !isnothing(r_channels)
                line = [take!(rc) for rc in r_channels]
            end
            for td in tree_dims
                trees = basic ? [parse_and_number(tree) for tree in line] :
                                [parse_and_number(tree[td]) for tree in line]
                check_leaves && check_leafsets(trees)

                # get all bipartitions
                cmds = Accumulator.(countmap.(get_bipartitions.(trees)))

                for (ind,acc) in enumerate(cmds)
                    merge!(splitsQueues[td][ind], acc)
                end
                new_splits = merge(cmds...)

                all_keys[td] = union(all_keys[td], keys(new_splits))
                merge!(splitsQueue[td], new_splits)

                tmp = 0.0
                M = 0.0
                for split in all_keys[td]
                    tmp_i = 0.0
                    keep = false
                    ova = splitsQueue[td][split]/(length(trees)*total_runs)
                    for r in 1:length(trees)
                        if splitsQueues[td][r][split]/total_runs > min_splits
                            keep = true
                        end
                        tmp_i += (splitsQueues[td][r][split]/total_runs - ova)^2

                    end
                    tmp_i /= (length(trees)-1)
                    tmp += keep ? sqrt(tmp_i) : 0.0
                    M += keep ? 1 : 0
                end
                ASDSF_vals[td][gen] = tmp / M

            end # for
            gen += 1
            total_runs += 1
            show_progress && ProgressMeter.next!(prog)
        end # if
    end # for
    show_progress && ProgressMeter.finish!(prog)
    conv_storage = ConvergenceStorage(splitsQueue, splitsQueues, total_runs)
    ASDSF_vals, conv_storage
end # ASDSF_int

max_stars_repo_name: erathorn/MCPhylo

max_issues_repo_name: erathorn/MCPhylo

max_forks_repo_name: erathorn/MCPhylo


text:
module Pipe
#Reflow piped things replacing the _ in the next section

export @pipe

const PLACEHOLDER = :_

function replace(arg::Any, target)
    arg #Normally do nothing
end

function replace(arg::Symbol, target)
    if arg==PLACEHOLDER
        target
    else
        arg
    end
end

function replace(arg::Expr, target)
    rep = copy(arg)
    rep.args = map(x->replace(x, target), rep.args)
    rep
end

function rewrite(ff::Expr, target)
    rep_args = map(x->replace(x, target), ff.args)
    if ff.args != rep_args
        #_ subsitution
        ff.args = rep_args
        return ff
    end

    #No subsitution was done (no _ found)
    #Apply to a function that is being returned by ff,
    #(ff could be a function call or something more complex)
    rewrite_apply(ff,target)
end

function rewrite_broadcasted(ff::Expr, target)
    temp_var = gensym()
    rep_args = map(x->replace(x, temp_var), ff.args)
    if ff.args != rep_args
        #_ subsitution
        ff.args = rep_args
        return :($temp_var->$ff)
    end

    #No subsitution was done (no _ found)
    #Apply to a function that is being returned by ff,
    #(ff could be a function call or something more complex)
    rewrite_apply_broadcasted(ff,target)
end

function rewrite_apply(ff, target)
    :($ff($target)) #function application
end

function rewrite_apply_broadcasted(ff, target)
    temp_var = gensym()
    :($temp_var->$ff($temp_var))
end

function rewrite(ff::Symbol, target)
    if ff==PLACEHOLDER
        target
    else
        rewrite_apply(ff,target)
    end
end

function rewrite_broadcasted(ff::Symbol, target)
    if ff==PLACEHOLDER
        target
    else
        rewrite_apply_broadcasted(ff,target)
    end
end

function funnel(ee::Any) #Could be a Symbol could be a literal
    ee #first (left most) input
end

function funnel(ee::Expr)
    if (ee.args[1]==:|>)
        target = funnel(ee.args[2]) #Recurse
        rewrite(ee.args[3],target)
    elseif (ee.args[1]==:.|>)
        target = funnel(ee.args[2]) #Recurse
        rewritten = rewrite_broadcasted(ee.args[3],target)
        ee.args[3] = rewritten
        ee
    else
        #Not in a piping situtation
        ee #make no change
    end
end

macro pipe(ee)
    esc(funnel(ee))
end

end


text:
# Autogenerated wrapper script for jq_jll for aarch64-apple-darwin
export jq

JLLWrappers.@generate_wrapper_header("jq")
JLLWrappers.@declare_executable_product(jq)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        jq,
        "bin/jq",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
function main_func()
    ands::Array{Bool} = []
    ors::Array{Bool} = []
    xors::Array{Bool} = []
    for a in [false, true]
        for b in [false, true]
            push!(ands, (a & b))
            push!(ors, (a | b))
            push!(xors, (a ⊻ b))
        end
    end
    @assert(ands == [false, false, false, true])
    @assert(ors == [false, true, true, true])
    @assert(xors == [false, true, true, false])
    println(join(["OK"], " "))
end

function main()
    main_func()
end

main()


text:
using Underscores
using Utils

const Edge = Tuple{Int,String}
const Edges = Array{Edge,1}

function parse_file(f)
    graph = Dict{String,Edges}()
    for line in eachline(f)
        outer, inner = split(line[1:end-1], " bags contain ")
        @assert !haskey(graph, outer)
        edges = Edges()
        if !occursin("no other", inner)
            for target in split(inner[1:end-1], ", ")
                m = match(r"(\d+) (\w+ \w+)", target)
                push!(edges, (parse(Int, m[1]), m[2]))
            end
        end
        graph[outer] = edges
    end
    graph
end

function load()
    open(joinpath(@__DIR__, "07_color-bags.txt"), "r") do f
        parse_file(f)
    end
end

bags = load()
const mine = "shiny gold"

function find_path(bags, color, what)
    color == what && return true

    for edge in bags[color]
         if find_path(bags, edge[2], what)
            return true
         end
    end

    return false
end

function ex1(bags)
    count(@_ filter(_ != mine, keys(bags))) do rule
        find_path(bags, rule, mine)
    end
end

function ex2(bags, color)
    edges = bags[color]
    isempty(edges) && return 1

    1 + sum(edges) do (count, inner)
        count * ex2(bags, inner)
    end
end

ex2(bags) = ex2(bags, mine) - 1

println("At least one $mine bag: ", ex1(bags))
println("Bags required inside $mine: ", ex2(bags))

using Test

@testset "Handy Haversacks" begin
    input = """
        light red bags contain 1 bright white bag, 2 muted yellow bags.
        dark orange bags contain 3 bright white bags, 4 muted yellow bags.
        bright white bags contain 1 shiny gold bag.
        muted yellow bags contain 2 shiny gold bags, 9 faded blue bags.
        shiny gold bags contain 1 dark olive bag, 2 vibrant plum bags.
        dark olive bags contain 3 faded blue bags, 4 dotted black bags.
        vibrant plum bags contain 5 faded blue bags, 6 dotted black bags.
        faded blue bags contain no other bags.
        dotted black bags contain no other bags.
        """
    example = parse_file(IOBuffer(input))

    @testset "example 1" begin
        @test length(example) == 9
        @test ex1(example) == 4
    end

    @testset "example 2" begin
        @test ex2(example) == 32

        input2 = """
            shiny gold bags contain 2 dark red bags.
            dark red bags contain 2 dark orange bags.
            dark orange bags contain 2 dark yellow bags.
            dark yellow bags contain 2 dark green bags.
            dark green bags contain 2 dark blue bags.
            dark blue bags contain 2 dark violet bags.
            dark violet bags contain no other bags.
            """
        example2 = parse_file(IOBuffer(input2))

        @test ex2(example2) == 126
    end

    @testset "results" begin
        @test ex1(bags) == 229
        @test ex2(bags) == 6683
    end
end

max_stars_repo_name: CmdQ/AoC


text:
using Documenter, DistributedArrays

makedocs(
    modules = [DistributedArrays],
    format = Documenter.HTML(),
    sitename = "DistributedArrays.jl",
    pages = [
        "Introduction" => "index.md"
        "API" => "api.md"
    ],
    doctest = true
)

deploydocs(
    repo = "github.com/JuliaParallel/DistributedArrays.jl.git",
)


text:
module Opsimath # module

using Printf

abstract type AbstractOpsimath{T} end

include("./units.jl")
include("./pyplot.jl")
include("./materials/materials.jl")
include("./spectra/spectra.jl")

end # module


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

# deep copying

# Note: deepcopy_internal(::Any, ::IdDict) is
#       only exposed for specialization by libraries

"""
    deepcopy(x)

Create a deep copy of `x`: everything is copied recursively, resulting in a fully
independent object. For example, deep-copying an array produces a new array whose elements
are deep copies of the original elements. Calling `deepcopy` on an object should generally
have the same effect as serializing and then deserializing it.

As a special case, functions can only be actually deep-copied if they are anonymous,
otherwise they are just copied. The difference is only relevant in the case of closures,
i.e. functions which may contain hidden internal references.

While it isn't normally necessary, user-defined types can override the default `deepcopy`
behavior by defining a specialized version of the function
`deepcopy_internal(x::T, dict::IdDict)` (which shouldn't otherwise be used),
where `T` is the type to be specialized for, and `dict` keeps track of objects copied
so far within the recursion. Within the definition, `deepcopy_internal` should be used
in place of `deepcopy`, and the `dict` variable should be
updated as appropriate before returning.
"""
deepcopy(x) = deepcopy_internal(x, IdDict())::typeof(x)

deepcopy_internal(x::Union{Symbol,Core.MethodInstance,Method,GlobalRef,DataType,Union,Task},
                  stackdict::IdDict) = x
deepcopy_internal(x::Tuple, stackdict::IdDict) =
    ntuple(i->deepcopy_internal(x[i], stackdict), length(x))
deepcopy_internal(x::Module, stackdict::IdDict) = error("deepcopy of Modules not supported")

function deepcopy_internal(x::SimpleVector, stackdict::IdDict)
    if haskey(stackdict, x)
        return stackdict[x]
    end
    y = Core.svec(Any[deepcopy_internal(x[i], stackdict) for i = 1:length(x)]...)
    stackdict[x] = y
    return y
end

function deepcopy_internal(x::String, stackdict::IdDict)
    if haskey(stackdict, x)
        return stackdict[x]
    end
    y = GC.@preserve x unsafe_string(pointer(x), sizeof(x))
    stackdict[x] = y
    return y
end

function deepcopy_internal(@nospecialize(x), stackdict::IdDict)
    T = typeof(x)::DataType
    nf = nfields(x)
    (isbitstype(T) || nf == 0) && return x
    if haskey(stackdict, x)
        return stackdict[x]
    end
    y = ccall(:jl_new_struct_uninit, Any, (Any,), T)
    if T.mutable
        stackdict[x] = y
    end
    for i in 1:nf
        if isdefined(x,i)
            ccall(:jl_set_nth_field, Cvoid, (Any, Csize_t, Any), y, i-1,
                  deepcopy_internal(getfield(x,i), stackdict))
        end
    end
    return y::T
end

function deepcopy_internal(x::Array, stackdict::IdDict)
    if haskey(stackdict, x)
        return stackdict[x]
    end
    _deepcopy_array_t(x, eltype(x), stackdict)
end

function _deepcopy_array_t(@nospecialize(x), T, stackdict::IdDict)
    if isbitstype(T)
        return (stackdict[x]=copy(x))
    end
    dest = similar(x)
    stackdict[x] = dest
    for i = 1:(length(x)::Int)
        if ccall(:jl_array_isassigned, Cint, (Any, Csize_t), x, i-1) != 0
            xi = ccall(:jl_arrayref, Any, (Any, Csize_t), x, i-1)
            if !isbits(xi)
                xi = deepcopy_internal(xi, stackdict)
            end
            ccall(:jl_arrayset, Cvoid, (Any, Any, Csize_t), dest, xi, i-1)
        end
    end
    return dest
end

function deepcopy_internal(x::Union{Dict,IdDict}, stackdict::IdDict)
    if haskey(stackdict, x)
        return stackdict[x]::typeof(x)
    end

    if isbitstype(eltype(x))
        return (stackdict[x] = copy(x))
    end

    dest = empty(x)
    stackdict[x] = dest
    for (k, v) in x
        dest[deepcopy_internal(k, stackdict)] = deepcopy_internal(v, stackdict)
    end
    dest
end

max_stars_repo_name: greimel/julia

max_issues_repo_name: greimel/julia

max_forks_repo_name: greimel/julia


text:
using QuasiDefinite
using Base.Test

# generate random nxn spd matrix
function gen_spd(n,seed=0)
    srand(seed)
    A = rand(n,n)
    A = A*A' + eye(n)
    A
end

function test_qdtf2!(A,uplo::Char)
    m, n = size(A)
    LD = copy(A)
    LD, info = QuasiDefinite.qdtf2!(uplo,LD)
    @test info == 0
    D = diagm(diag(LD))
    if uplo == 'L'
        L = tril(LD,-1) + eye(m)
        U = triu(LD)
        @test_approx_eq_eps norm(L*D*L' - A,1) 0.0 1e-12
    else
        L = tril(LD)
        U = triu(LD,1) + eye(m)
        @test_approx_eq_eps norm(U'*D*U - A,1) 0.0 1e-12
    end
    @test_approx_eq_eps norm(L*U - A,1) 0.0 1e-12
end

# set size and generate matrix
n = 35
A = gen_spd(n)

# test unblocked LDL^T
test_qdtf2!(A,'L')
test_qdtf2!(A,'U')


text:
# ------------------------------------------------------------------------------------------
# # Convolutional Neural Networks
# (c) Deniz Yuret, 2019
# * Objectives: See the effect of sparse and shared weights implemented by convolutional
# networks.
# * Prerequisites: [MLP models](40.mlp.ipynb), [MNIST](20.mnist.ipynb)
# * New functions:
# [conv4](http://denizyuret.github.io/Knet.jl/latest/reference/#Knet.conv4),
# [pool](http://denizyuret.github.io/Knet.jl/latest/reference/#Knet.pool),
# [mat](http://denizyuret.github.io/Knet.jl/latest/reference/#Knet.mat)
#
# ![image](https://github.com/denizyuret/Knet.jl/blob/master/docs/src/images/le_net.png?raw=
# true)
# ([image source](http://www.dataiku.com/blog/2015/08/18/Deep_Learning.html))
#
# To improve the performance further, we can use a convolutional neural networks (CNN). See
# the [course notes](http://cs231n.github.io/convolutional-networks/) by Andrej Karpathy for
# a good introduction to CNNs. We will implement the
# [LeNet](http://yann.lecun.com/exdb/lenet) model which consists of two convolutional layers
# followed by two fully connected layers. We will describe and use the
# [conv4](http://denizyuret.github.io/Knet.jl/latest/reference/#Knet.conv4) and
# [pool](http://denizyuret.github.io/Knet.jl/latest/reference/#Knet.pool) functions provided
# by Knet for the implementation of convolutional nets.
#
# Even though MLPs and CNNs are both universal function approximators and both achieve 0
# error on the training set, we will see that a CNN converges a lot faster and generalizes a
# lot better with less overfitting achieving a 99.5% test accuracy on MNIST. The sparse
# connectivity and shared weights of a CNN give it an inductive bias appropriate for image
# features allowing it to learn better with less data.
# ------------------------------------------------------------------------------------------

# Setup display width, load packages, import symbols
ENV["COLUMNS"]=72
using Pkg; for p in ("Knet","Plots"); haskey(Pkg.installed(),p) || Pkg.add(p); end
using Base.Iterators: flatten
using Statistics: mean
using Knet: Knet, conv4, pool, mat, KnetArray, nll, zeroone, progress, sgd, param, param0, dropout, relu, Data

# ------------------------------------------------------------------------------------------
# ## Introduction to convolution
# ------------------------------------------------------------------------------------------

# Convolution operator in Knet
@doc conv4

# Convolution in 1-D
w = reshape([1.0,2.0,3.0], (3,1,1,1)); @show w
x = reshape([1.0:7.0...], (7,1,1,1)); @show x
@show y = conv4(w, x);  # size Y = X - W + 1 = 5 by default

# Padding
w = reshape([1.0,2.0,3.0], (3,1,1,1)); @show w
x = reshape([1.0:7.0...], (7,1,1,1)); @show x
@show y2 = conv4(w, x, padding=(1,0));  # size Y = X + 2P - W + 1 = 7 with padding=1
# To preserve input size (Y=X) for a given W, what padding P should we use?

# Stride
w = reshape([1.0,2.0,3.0], (3,1,1,1)); @show w
x = reshape([1.0:7.0...], (7,1,1,1)); @show x
@show y3 = conv4(w, x; padding=(1,0), stride=3);  # size Y = 1 + floor((X+2P-W)/S)

# Mode
w = reshape([1.0,2.0,3.0], (3,1,1,1)); @show w
x = reshape([1.0:7.0...], (7,1,1,1)); @show x
@show y4 = conv4(w, x, mode=0);  # Default mode (convolution) inverts w
@show y5 = conv4(w, x, mode=1);  # mode=1 (cross-correlation) does not invert w

# Convolution in more dimensions
x = reshape([1.0:9.0...], (3,3,1,1))

w = reshape([1.0:4.0...], (2,2,1,1))

y = conv4(w, x)

# Convolution with multiple channels, filters, and instances
# size X = [X1,X2,...,Xd,Cx,N] where d is the number of dimensions, Cx is channels, N is instances
x = reshape([1.0:18.0...], (3,3,2,1))

# size W = [W1,W2,...,Wd,Cx,Cy] where d is the number of dimensions, Cx is input channels, Cy is output channels
w = reshape([1.0:24.0...], (2,2,2,3));

# size Y = [Y1,Y2,...,Yd,Cy,N]  where Yi = 1 + floor((Xi+2Pi-Wi)/Si), Cy is channels, N is instances
y = conv4(w,x)

# ------------------------------------------------------------------------------------------
# See http://cs231n.github.io/assets/conv-demo/index.html for an animated example.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Introduction to Pooling
# ------------------------------------------------------------------------------------------

# Pooling operator in Knet
@doc pool

# 1-D pooling example
x = reshape([1.0:6.0...], (6,1,1,1)); @show x
@show pool(x);

# Window size
x = reshape([1.0:6.0...], (6,1,1,1)); @show x
@show pool(x; window=3);  # size Y = floor(X/W)

# Padding
x = reshape([1.0:6.0...], (6,1,1,1)); @show x
@show pool(x; padding=(1,0));  # size Y = floor((X+2P)/W)

# Stride
x = reshape([1.0:10.0...], (10,1,1,1)); @show x
@show pool(x; stride=4);  # size Y = 1 + floor((X+2P-W)/S)

# Mode (using KnetArray here; not all modes are implemented on the CPU)
x = KnetArray(reshape([1.0:6.0...], (6,1,1,1))); @show x
@show pool(x; padding=(1,0), mode=0)  # max pooling
@show pool(x; padding=(1,0), mode=1)  # avg pooling
@show pool(x; padding=(1,0), mode=2); # avg pooling excluding padded values (is not implemented on CPU)

# More dimensions
x = reshape([1.0:16.0...], (4,4,1,1))

pool(x)

# Multiple channels and instances
x = reshape([1.0:32.0...], (4,4,2,1))

# each channel and each instance is pooled separately
pool(x)  # size Y = (Y1,...,Yd,Cx,N) where Yi are spatial dims, Cx and N are identical to input X

# ------------------------------------------------------------------------------------------
# ## Experiment setup
# ------------------------------------------------------------------------------------------

# Load data (see mnist.ipynb)
include(Knet.dir("data","mnist.jl"))  # Load data
dtrn,dtst = mnistdata();              # dtrn and dtst = [ (x1,y1), (x2,y2), ... ] where xi,yi are minibatches of 100

(x,y) = first(dtst)
println.(summary.((x,y)));

# For running experiments
function trainresults(file,model; o...)
    if (print("Train from scratch? "); readline()[1]=='y')
        takeevery(n,itr) = (x for (i,x) in enumerate(itr) if i % n == 1)
        r = ((model(dtrn), model(dtst), zeroone(model,dtrn), zeroone(model,dtst))
             for x in takeevery(length(dtrn), progress(sgd(model,repeat(dtrn,100)))))
        r = reshape(collect(Float32,flatten(r)),(4,:))
        Knet.save(file,"results",r)
        Knet.gc() # To save gpu memory
    else
        isfile(file) || download("http://people.csail.mit.edu/deniz/models/tutorial/$file",file)
        r = Knet.load(file,"results")
    end
    println(minimum(r,dims=2))
    return r
end

# ------------------------------------------------------------------------------------------
# ## A convolutional neural network model for MNIST
# ------------------------------------------------------------------------------------------

# Define a convolutional layer:
struct Conv; w; b; f; p; end
(c::Conv)(x) = c.f.(pool(conv4(c.w, dropout(x,c.p)) .+ c.b))
Conv(w1::Int,w2::Int,cx::Int,cy::Int,f=relu;pdrop=0) = Conv(param(w1,w2,cx,cy), param0(1,1,cy,1), f, pdrop)

# Redefine dense layer (See mlp.ipynb):
struct Dense; w; b; f; p; end
(d::Dense)(x) = d.f.(d.w * mat(dropout(x,d.p)) .+ d.b) # mat reshapes 4-D tensor to 2-D matrix so we can use matmul
Dense(i::Int,o::Int,f=relu;pdrop=0) = Dense(param(o,i), param0(o), f, pdrop)

# Let's define a chain of layers
struct Chain
    layers
    Chain(layers...) = new(layers)
end
(c::Chain)(x) = (for l in c.layers; x = l(x); end; x)
(c::Chain)(x,y) = nll(c(x),y)
(c::Chain)(d::Data) = mean(c(x,y) for (x,y) in d)

lenet =   Chain(Conv(5,5,1,20),
                Conv(5,5,20,50),
                Dense(800,500,pdrop=0.3),
                Dense(500,10,identity,pdrop=0.3))
summary.(l.w for l in lenet.layers)

lenet(x,y)

# ------------------------------------------------------------------------------------------
# ## CNN vs MLP
# ------------------------------------------------------------------------------------------

# 1.08e-02  100.00%┣████████████████▉┫ 60000/60000 [03:50/03:50, 260.67i/s]
# [0.000135032; 0.0196918; 0.0; 0.0053]
cnn = trainresults("cnn113.jld2", lenet);

mlp = Knet.load("mlp113f.jld2","results");

using Plots; default(fmt=:png,ls=:auto)

# Comparison to MLP shows faster convergence, better generalization
plot([mlp[1,:], mlp[2,:], cnn[1,:], cnn[2,:]],ylim=(0.0,0.1),
     labels=[:trnMLP :tstMLP :trnCNN :tstCNN],xlabel="Epochs",ylabel="Loss")

plot([mlp[3,:], mlp[4,:], cnn[3,:], cnn[4,:]],ylim=(0.0,0.03),
    labels=[:trnMLP :tstMLP :trnCNN :tstCNN],xlabel="Epochs",ylabel="Error")

# ------------------------------------------------------------------------------------------
# ## Convolution vs Matrix Multiplication
# ------------------------------------------------------------------------------------------

# Convolution and matrix multiplication can be implemented in terms of each other.
# Convolutional networks have no additional representational power, only statistical efficiency.
# Our original 1-D example
@show w = reshape([1.0,2.0,3.0], (3,1,1,1))
@show x = reshape([1.0:7.0...], (7,1,1,1))
@show y = conv4(w, x);  # size Y = X - W + 1 = 5 by default

# Convolution as matrix multiplication (1)
# Turn w into a (Y,X) sparse matrix
w2 = Float64[3 2 1 0 0 0 0; 0 3 2 1 0 0 0; 0 0 3 2 1 0 0; 0 0 0 3 2 1 0; 0 0 0 0 3 2 1]

@show y2 = w2 * mat(x);

# Convolution as matrix multiplication (2)
# Turn x into a (W,Y) dense matrix (aka the im2col operation)
# This is used to speed up convolution with known efficient matmul algorithms
x3 = Float64[1 2 3 4 5; 2 3 4 5 6; 3 4 5 6 7]

@show w3 = [3.0 2.0 1.0]
@show y3 = w3 * x3;

# Matrix multiplication as convolution
# This could be used to make a fully connected network accept variable sized inputs.
w = reshape([1.0:6.0...], (2,3))

x = reshape([1.0:3.0...], (3,1))

y = w * x

# Consider w with size (Y,X)
# Treat each of the Y rows of w as a convolution filter
w2 = copy(reshape(Array(w)', (3,1,1,2)))

# Reshape x for convolution
x2 = reshape(x, (3,1,1,1))

# Use conv4 for matrix multiplication
y2 = conv4(w2, x2; mode=1)

# ------------------------------------------------------------------------------------------
# * So there is no difference between the class of functions representable with an MLP vs
# CNN.
# * Sparse connections and weight sharing give CNNs more generalization power with images.
# * Number of parameters in MLP256: (256x784)+256+(10x256)+10 = 203530
# * Number of parameters in LeNet: (5*5*1*20)+20+(5*5*20*50)+50+(500*800)+500+(10*500)+10 =
# 431080
# ------------------------------------------------------------------------------------------

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
# This program generates a 3D synthetic seismic section by convolution with
# a the Ricker wavelet defined in Seismic.jl's ShotProfileWEM. To do this
# it reads a velocity model in, calculates reflection coeficients, and uses
# a 1D convolution per trace in the frequency domain.
#
# function call vars:
# peakF is the Ricker peak frequency, and samp is the sampling interval
#
# by GRAM | 8 Apr 2016

using Seismic

function modeler(vel="", out="", peakF=100, samp=0.001, verbose=0)

	# import vel model
	vel,vel_h = SeisRead(vel)
	nz = size(vel)[1]
	nx = size(vel)[2]
	ny = size(vel)[3]

	# initialize
	dz = 10
	dx = 20
	dy = 20

	# build a wavelet for modeling
	temp = ricker(peakF, samp)
	rick = zeros(Float32, size(temp)[1])
	for i = 1 : size(temp)[1]
		rick[i] = convert(Float32, temp[i])
	end

	if verbose > 0
		print("\nThis Ricker wavelet is ",size(rick)[1]," samples long.\n\n")
	end

	#write out wavelet for use with finite difference modeling
	ex = Seismic.Extent(convert(Int32,size(rick)[1]), convert(Int32,1), convert(Int32,1),
		1, 1, 0, 0, 0, 0, 0, convert(Float32,dz), convert(Float32,1),
		convert(Float32,1), 1, 1, "Depth", "mx", "0", "", "", "", "", "",
		"", "", "")
	h_w = Array(Header,1);
	h_w[1] = Seismic.InitSeisHeader();
	h_w[1].tracenum = 1;
	h_w[1].n1 = size(rick)[1];
	h_w[1].d1 = samp;
	SeisWrite("../dat/wav", rick[:], h_w, ex);

	# calculate (normal incidence) reflection coeficients
	r = zeros(Float32,nz,nx,ny)
	d1 = d2 = 1 # acoustic assumption

	for k = 1 : ny
		for j = 1 : nx
			for i = 1 : nz - 1
				r[i,j,k] = ((d2 * vel[i+1,j,k] - d1 * vel[i,j,k]) /
					   (d2 * vel[i+1,j,k] + d1 * vel[i,j,k])) ^ 2
			end
		end
	end

#=
	# write out reflection coefficients
	h = Array(Header,nx*ny)
	ex = Seismic.Extent(convert(Int32,nz), convert(Int32,nx), convert(Int32,ny),
		1, 1, 0, 0, 0, 0, 0, convert(Float32,dz), convert(Float32,dx),
		convert(Float32,dy), 1, 1, "Depth", "mx", "my", "", "", "", "", "",
		"", "", "")

	for i = 1 : ny
		for j = 1 : nx
			h[j + nx * (i - 1)] = Seismic.InitSeisHeader()
			h[j + nx * (i - 1)].tracenum = convert(Int32, j + nx * (i - 1))
			h[j + nx * (i - 1)].o1 = convert(Float32, 0)
			h[j + nx * (i - 1)].n1 = convert(Int32, nz)
			h[j + nx * (i - 1)].d1 = convert(Float32, dz)
			h[j + nx * (i - 1)].mx = convert(Float32, j)
			h[j + nx * (i - 1)].my = convert(Float32, i)
		end
	end

	SeisWrite("../dat/refl",r,vel_h,ex)
=#

	# convolve RC with Ricker wavelet to generate 0 offset section
	nz_new = size(rick)[1] + nz - 1
	m = zeros(Float32, nz_new, nx, ny)
	for j = 1 : ny
		for i = 1 : nx
			m[:,i,j] = conv(r[:,i,j], rick)
		end
	end

	# window out middle samples of convolution for physical reality
	extra = size(m)[1] - nz
	if isodd(extra) == true
		for j = 1 : ny
			for i = 1 : nx
				vel[:,i,j] = m[Int(ceil(extra / 2)) : nz - 1 + Int(ceil(extra / 2)), i, j]
			end
		end
	else
		for j = 1 : ny
			for i = 1 : nx
				vel[:,i,j] = m[Int(extra / 2) : nz - 1 + Int(extra / 2), i, j]
			end
		end
	end

	# write out the wavefield
	h = Array(Header,nx*ny)
	ex = Seismic.Extent(convert(Int32,nz), convert(Int32,nx), convert(Int32,ny),
		1, 1, 0, 0, 0, 0, 0, convert(Float32,dz), convert(Float32,dx),
		convert(Float32,dy), 1, 1, "Depth", "mx", "my", "", "", "", "", "",
		"", "", "")

	SeisWrite(out,vel,vel_h,ex)

end # end function


text:
using Distributed
using Random

# @everywhere using JuliaRL
@everywhere include("test_tilecoder_mountaincar.jl")

function run_mountain_car_parallel(n_runs=10, n_episodes=10000)

    experiment(run) = mountain_car_test(0.5/8, 0.1, 1.0, 32, 8; n_episodes=n_episodes, seed=1002349+run)

    futures = Array{Future, 1}()

    for run in 1:n_runs
        append!(futures, [@spawn experiment(run)])
    end

    results = zeros(n_runs, n_episodes)

    @showprogress 1.0 "Run Fetched: " for run in 1:n_runs
        results[run, :] .= fetch(futures[run])
    end

    return results

end

max_stars_repo_name: mkschleg/JuliaRL

max_issues_repo_name: mkschleg/JuliaRL

max_forks_repo_name: mkschleg/JuliaRL


text:
immutable FuncExpr
    name::Symbol
    slots
    args
    body::Vector
    returntype::Type
end

function FuncExpr(f, types)
    if isa(f, Core.IntrinsicFunction)
        error("$f is an intrinsic function which can't be transpiled")
    end
    lam_typed, ret_type = get_lambda(code_typed, f, types)
    slots = slot_mapping(lam_typed)
    ast = lam_typed.code
    name = Symbol(f)
    FuncExpr(name, slots, view(slots, 2:(length(slots))), ast, ret_type)
end

function get_ast(li::CodeInfo)
    ast = li.code
    if isa(ast, Vector{UInt8})
        return Base.uncompressed_ast(li)
    end
    ast
end


text:
#= src/lapack.jl
=#

module _LAPACK
import LinearAlgebra.BLAS: @blasfunc
using LinearAlgebra: BlasFloat, BlasInt, chkstride1
import Base: Nothing
const liblapack = Base.liblapack_name

for (geqrf, gelqf, orgqr, orglq, elty) in
      ((:sgeqrf_, :sgelqf_, :sorgqr_, :sorglq_, :Float32   ),
       (:dgeqrf_, :dgelqf_, :dorgqr_, :dorglq_, :Float64   ),
       (:cgeqrf_, :cgelqf_, :cungqr_, :cunglq_, :ComplexF32 ),
       (:zgeqrf_, :zgelqf_, :zungqr_, :zunglq_, :ComplexF64))
  @eval begin
    function geqrf!(
        A::AbstractMatrix{$elty}, tau::Vector{$elty}, work::Vector{$elty})
      chkstride1(A)
      m, n  = size(A)
      k     = min(m, n)
      tau   = length(tau) < k ? Array{$elty}(undef, k) : tau
      lwork = BlasInt(-1)
      info  = Ref{BlasInt}()
      for i = 1:2
        ccall((@blasfunc($geqrf), liblapack), Nothing,
              (Ref{BlasInt}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt},
               Ptr{$elty}, Ptr{$elty}, Ref{BlasInt}, Ptr{BlasInt}),
              m, n, A, max(1,stride(A,2)),
              tau, work, lwork, info)
        if i == 1
          lwork = BlasInt(real(work[1]))
          work  = length(work) < lwork ? Array{$elty}(undef, lwork) : work
          lwork = length(work)
        end
      end
      A, tau, work
    end

    function gelqf!(
        A::AbstractMatrix{$elty}, tau::Vector{$elty}, work::Vector{$elty})
      chkstride1(A)
      m, n  = size(A)
      k     = min(m, n)
      tau   = length(tau) < k ? Array{$elty}(undef, k) : tau
      lwork = BlasInt(-1)
      info  = Ref{BlasInt}()
      for i = 1:2
        ccall((@blasfunc($gelqf), liblapack), Nothing,
              (Ref{BlasInt}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt},
               Ptr{$elty}, Ptr{$elty}, Ref{BlasInt}, Ptr{BlasInt}),
              m, n, A, max(1,stride(A,2)),
              tau, work, lwork, info)
        if i == 1
            lwork = BlasInt(real(work[1]))
            work  = length(work) < lwork ? Array{$elty}(undef, lwork) : work
            lwork = length(work)
        end
      end
      A, tau, work
    end

    function orglq!(
        A::AbstractMatrix{$elty}, tau::Vector{$elty}, k::Integer,
        work::Vector{$elty})
      chkstride1(A)
      n = size(A, 2)
      m = min(n, size(A, 1))
      0 <= k <= min(m, length(tau)) || throw(DimensionMismatch)
      lwork = BlasInt(-1)
      info  = Ref{BlasInt}()
      for i = 1:2
          ccall((@blasfunc($orglq), liblapack), Nothing,
                (Ref{BlasInt}, Ref{BlasInt}, Ref{BlasInt}, Ptr{$elty},
                 Ref{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ref{BlasInt},
                 Ptr{BlasInt}),
                m, n, k, A,
                max(1,stride(A,2)), tau, work, lwork,
                info)
          if i == 1
              lwork = BlasInt(real(work[1]))
              work  = length(work) < lwork ? Array{$elty}(undef, lwork) : work
              lwork = length(work)
          end
      end
      A, tau, work
    end

    function orgqr!(
        A::AbstractMatrix{$elty}, tau::Vector{$elty}, k::Integer,
        work::Vector{$elty})
      chkstride1(A)
      m = size(A, 1)
      n = min(m, size(A,2))
      0 <= k <= min(n, length(tau)) || throw(DimensionMismatch)
      lwork = BlasInt(-1)
      info  = Ref{BlasInt}()
      for i = 1:2
        ccall((@blasfunc($orgqr), liblapack), Nothing,
              (Ref{BlasInt}, Ref{BlasInt}, Ref{BlasInt}, Ptr{$elty},
               Ref{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ref{BlasInt},
               Ptr{BlasInt}),
              m, n, k, A,
              max(1,stride(A,2)), tau, work, lwork,
              info)
        if i == 1
          lwork = BlasInt(real(work[1]))
          work  = length(work) < lwork ? Array{$elty}(undef, lwork) : work
          lwork = length(work)
        end
      end
      A, tau, work
    end
  end
end

for (laqps, elty, relty) in ((:slaqps_, :Float32,    :Float32),
                             (:dlaqps_, :Float64,    :Float64),
                             (:claqps_, :ComplexF32,  :Float32),
                             (:zlaqps_, :ComplexF64, :Float64))
  @eval begin
    function laqps!(
        offset::BlasInt, nb::BlasInt, kb::Ref{BlasInt},
        A::AbstractMatrix{$elty},
        jpvt::AbstractVector{BlasInt}, tau::AbstractVector{$elty},
        vn1::AbstractVector{$relty}, vn2::AbstractVector{$relty},
        auxv::AbstractVector{$elty}, F::AbstractVector{$elty})
      m, n = size(A)
      ccall(
        (@blasfunc($laqps), liblapack), Nothing,
        (Ref{BlasInt}, Ref{BlasInt}, Ref{BlasInt}, Ref{BlasInt}, Ptr{BlasInt},
         Ptr{$elty}, Ref{BlasInt}, Ptr{BlasInt}, Ptr{$elty},
         Ptr{$relty}, Ptr{$relty}, Ptr{$elty}, Ptr{$elty}, Ref{BlasInt}),
        m, n, offset, nb, kb,
        A, max(1,stride(A,2)), jpvt, tau,
        vn1, vn2, auxv, F, max(1,n))
    end
  end
end

end  # module


text:
# managment of MPI tags

"""
  Manages MPI tags.  Specifically, keeps track of which tags are in use and
  provides for getting unused tags.  Currently it does not consider which tags
  are in use with which receiving processes.

  **Implementation Notes**

  The management scheme attempts to re-use tags as soon as they are freed.
  This requires an O(N) search in some cases to figure out the next tag.
  This could be optimized in various ways, but it doesn't seem worth it
  considering only a handful of tags are used at any one time.
"""
mutable struct MPITagManager
  used_tags::BitArray{1}  # true if tag in inuse, false otherwise
                          # if large numbers of tags are used, this could be
                          # a SparseVector{Bool} instead
  next_tag::Cint
  start_tag::Cint
end

"""
  Outer constrcutor for `MPITagManager`](@ref)

  **Inputs**

   * starting tag: an integer specifying the first tag the tag manager will
                   use.  Using tags lower than this value is not allowed.
"""
function MPITagManager(starting_tag::Integer=1)

  used_tags = BitArray(0)
  next_tag = starting_tag

  return MPITagManager(used_tags, next_tag, starting_tag)
end

"""
  Returns the next free tag

  **Inputs**

   * mgr: an MPITagManager

  **Outputs**

   * a Cint containing the tag
"""
function getNextTag(mgr::MPITagManager)

  next_tag = mgr.next_tag
  next_tag_offset = next_tag - mgr.start_tag + 1

  resize_array = next_tag_offset > length(mgr.used_tags)
  if resize_array
    resize!(mgr.used_tags, next_tag_offset)
  end

  mgr.used_tags[next_tag_offset] = true

  # figure out what the next tag will be
  if resize_array
    new_next_tag = next_tag + 1
  else
    # find next tag using O(N) search
    # hopefully this is the rare case because tags mostly get used and rarely
    # get freed
    new_next_tag = findfirst(mgr.used_tags, false) + mgr.start_tag - 1
    if new_next_tag == mgr.start_tag - 1
      new_next_tag = length(mgr.used_tags) + mgr.start_tag
    end
  end

  mgr.next_tag = new_next_tag
  checkNextTag(new_next_tag)

  return next_tag
end

"""
  Marks a tag as being in use.  This is useful if (for some reason) the user
  wants to use a specific tag number rather than calling [`getNextTag`](@ref)

  **Inputs**

   * mgr: [`MPITagManager`](@ref)
   * tag: an integer specifying the tag to be marked.  This tag must be
          greater than the starting tag passed to the `mgr` constructor
"""
function markTagUsed(mgr::MPITagManager, tag::Integer)

  if tag < mgr.start_tag
    error("cannot use tags less than start tag: $(mgr.start_tag)")
  end

  tag_offset = tag - mgr.start_tag + 1
  if tag_offset > length(mgr.used_tags)
    resize!(mgr.used_tags, tag_offset)
  end

  if mgr.used_tags[tag_offset]
    error("attempting to mark already-used tag $tag as used")
  end

  mgr.used_tags[tag_offset] = true

  return nothing
end

"""
  Marks a tag to be freed.  It is an error to free a tag that is not used.

  **Inputs**

   * tag: integer specifying the tag to be freed
"""
function freeTag(mgr::MPITagManager, tag::Integer)

  tag_offset = tag - mgr.start_tag + 1
  if tag < mgr.start_tag
    error("cannot use tags less than start tag: $(mgr.start_tag)")
  end

  if !mgr.used_tags[tag_offset]
    error("cannot free tag $tag, which is not in use")
  end

  mgr.used_tags[tag_offset] = false

  # figure out what tag to use next.  If possible, reuse this tag, but only
  # if the next tag is greater
  # Basically, try to use the first available tag
  if mgr.next_tag > tag
    mgr.next_tag = tag
  end

  checkNextTag(mgr.next_tag)

  return nothing
end

"""
  Internal function, throws an error if the maximum tag is reached

  **Inputs**

    * next_tag: next tag to be used
"""
function checkNextTag(next_tag::Integer)

  # this triggers one tag early, but that will have to do
  if next_tag == typemax(Cint)
    error("reached maximum MPI tag.  Consider freeing some tags (also, how did you manage to use so many tags in the first place?)")
  end
end


text:
using Parameters, NLsolve, SpecialFunctions, Distributions, Roots, ArbNumerics, StatsBase, LsqFit, PoissonRandom, SparseArrays, Distributed, CSV, SharedArrays, JLD2, DataFrames, ProgressMeter

import GZip: open
import Parsers: parse
import OrderedCollections: OrderedDict
import FastaIO: readfasta

include("/home/jmurga/.julia/dev/Analytical/src/parameters.jl")
include("/home/jmurga/.julia/dev/Analytical/src/fixations.jl")
include("/home/jmurga/.julia/dev/Analytical/src/polymorphism.jl")
include("/home/jmurga/.julia/dev/Analytical/src/summaryStatistics.jl")
include("/home/jmurga/.julia/dev/Analytical/src/rates.jl")
include("/home/jmurga/.julia/dev/Analytical/src/inferTools.jl")
include("/home/jmurga/.julia/dev/Analytical/src/readFasta.jl")
include("/home/jmurga/.julia/dev/Analytical/src/methods.jl")

using Distributed
addprocs()
@everywhere using Analytical, DataFrames, CSV
# Set up model
adap = Analytical.parameters(N=500,n=500,gam_neg=-457, gL=10,gH=500,Lf=10^5,B=0.999,alTot=0.4,alLow=0.2)
Analytical.	binomOp!(adap);

sfs = convert(Array,DataFrame!(CSV.File("/home/jmurga/mkt/202004/rawData/simulations/noDemog/noDemog_0.4_0.1_0.999/sfs.tsv")))
rSfs = convert.(Int64,Analytical.reduceSfs(Analytical.cumulativeSfs(sfs[:,2:end]),100))'
sCumu = Analytical.cumulativeSfs(sfs[:,2:end])

sfsPos   = sCumu[:,1] + sCumu[:,2]
sfsNopos = sCumu[:,4] + sCumu[:,2]

divergence = convert(Array,DataFrame!(CSV.File("/home/jmurga/mkt/202004/rawData/simulations/noDemog/noDemog_0.4_0.1_0.999/div.tsv")))
d = [convert(Float64,sum(divergence[1:2]))]

alpha = @. 1 - (divergence[2]/divergence[1] * rSfs[:,1]/rSfs[:,2])

inputAbc = DataFrame(alpha')

@time df = Analytical.summaryStats(adap,d,sfsPos,100,3);
@time df = Analytical.summaryStats(adap,583,d,sfsPos);

CSV.write("/home/jmurga/mkt/202004/rawData/summStat/noDemog/noDemog_0.4_0.2_0.999/noDemog_0.4_0.2_0.999.tsv", df, delim='\t',header=false);

# fac       = rand(-2:0.05:2,iterations)
# afac      = @. 0.184*(2^fac)
# bfac      = @. 0.000402*(2^fac)
#
# alTot       = rand(collect(0.01:0.01:0.6),iterations)
# lfac        = rand(collect(0.1:0.1:0.9),iterations)
# alLow       = @. round(alTot * lfac,digits=5)
# nParam      = [param for i in 1:iterations]
# ndivergence = [divergence for i in 1:iterations]
# nSfs        = [sfs for i in 1:iterations]
# nBins        = [bins for i in 1:iterations]
#
# for i in zip(afac,bfac,alTot,alLow)
#
#     param.al = i[1]; param.be = i[2];
#     param.alLow = i[3]; param.alTot = i[3]
#
#     # Solve probabilites without B effect to get to achieve α value
#     param.B = 0.999
#     Analytical.set_theta_f!(param)
#     @time Analytical.setPpos!(param)
# end

## Open empirical data
# path= "/home/jmurga/mkt/202004/rawData/";suffix="txt";
# files = path .* filter(x -> occursin(suffix,x), readdir(path))[1]

# pol,sfs,div = Analytical.parseSfs(param=adap,data=files,output="/home/jmurga/data",sfsColumns=[3,5],divColumns=[6,7],bins=100)
# sfs = sfs[:,1]
# div = [div[1]]

# CSV.write("/home/jmurga/mkt/202004/rawData/summStat/noDemog/noDemog_0.4_0.2_0.999/sfsnoDemog.tsv", inputAbc, delim='\t',header=false);
#
# function summStats(param::Analytical.parameters,iterations::Int64,div::Array,sfs::Array)
#
#     fac       = rand(-2:0.05:2,iterations)
#     afac      = @. 0.184*(2^fac)
#     bfac      = @. 0.000402*(2^fac)
#
#     alTot     = rand(collect(0.01:0.01:0.6),iterations)
#     lfac      = rand(collect(0.1:0.1:0.9),iterations)
#     alLow     = @. round(alTot * lfac,digits=5)
#     nParam = [param for i in 1:iterations]
#     nDiv = [div for i in 1:iterations]
#     nSfs = [sfs for i in 1:iterations]
#
#     wp = CachingPool(workers())
#     b = pmap(bgsIter,wp,nParam,afac,bfac,alTot,alLow,nDiv,nSfs);
# 	return(b)
# 	#=return(reduce(vcat,b))=#
# end
#
# @everywhere function bgsIter(param::Analytical.parameters,afac::Float64,bfac::Float64,alTot::Float64,alLow::Float64,div::Array,sfs::Array)
#
#     r = Array{Float64}(undef, 17, 100 + 3	)
#     # r = zeros(1,103)
#     param.al = afac; param.be = bfac;
#     param.alLow = alLow; param.alTot = alTot;
#     iter = 1
# 	for j in param.bRange
#         param.B = j
#
#         Analytical.set_theta_f!(param)
#         theta_f = param.theta_f
#         param.B = 0.999
#         Analytical.set_theta_f!(param)
#         Analytical.setPpos!(param)
#         param.theta_f = theta_f
#         param.B = j
#         x,y,z = Analytical.alphaByFrequencies(param,div,sfs,100,0.9)
#
#         r[iter,:] = z
#         iter = iter + 1;
#         # println(z)
#     end
#     # return(reduce(vcat,r))
#     return(r)
# end

# function test(p,d,s)
#     r=[]
#     for i in 1:1000
#         push!(r,bgsIter(p,d,s))
#     end
#     return(r)
# end

# reduce(vcat,reduce(vcat,a))

# Analytical.summaryStatistics(output, z)


text:
@testset "icd10" begin
    DiagnosisClass  = DiagnosisClassification.DiagnosisClass
    construct_icd10 = DiagnosisClassification.construct_icd10
    construct_ccs10 = DiagnosisClassification.construct_ccs10
    @test construct_icd10("W55.29XA") == DiagnosisClass(;
        system = "ICD10",
        value = "W55.29XA",
        )
    @test construct_icd10("W5529XA")  == DiagnosisClass(;
        system = "ICD10",
        value = "W55.29XA",
        )
    @test construct_ccs10("helloworld123")  == DiagnosisClass(;
        system = "CCS10",
        value = "helloworld123",
    )
end


text:
#=

  Zoo, buses and kids problem in Julia ConstraintSolver.jl

  https://dmcommunity.org/challenge/challenge-july-2018/
  """
  300 kids need to travel to the London zoo. The school may rent 40 seats and
  30 seats buses for 500 and 400 £. How many buses of each to minimize cost?
  """

  Also, see Alex Fleisher
  "Optimization : simply do more with less, zoo, buses and kids (Part2, python, java, C++…)"
  ""
  https://alexfleischer-84755.medium.com/optimization-simply-do-more-with-less-zoo-buses-and-kids-part2-python-java-c-cc04558e49b5

  This model was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
const CS = ConstraintSolver
include("constraints_utils.jl")

function zoo(problem,print_solutions=true,all_solutions=true)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   # "all_solutions"=> all_solutions,
                                                            "all_optimal_solutions"=>all_solutions,
                                                            "logging"=>[],

                                                            "traverse_strategy"=>:BFS,
                                                            # "traverse_strategy"=>:DFS,
                                                            # "traverse_strategy"=>:DBFS,

                                                            # "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            "branch_split"=>:InHalf,

                                                            # https://wikunia.github.io/ConstraintSolver.jl/stable/options/#branch_strategy-(:Auto)
                                                            "branch_strategy" => :IMPS, # default
                                                            # "branch_strategy" => :ABS, # Activity Based Search
                                                            # "activity.decay" => 0.999, # default 0.999
                                                            # "activity.max_probes" => 10, # default, 10
                                                            # "activity.max_confidence_deviation" => 20, # default 20

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>6,

                                                            # "backtrack" => false, # default true
                                                            # "backtrack_sorting" => false, # default true

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            # "lp_optimizer" => glpk_optimizer,
                                                            # "lp_optimizer" => ipopt_optimizer,
                                        ))
    num_kids  = problem[:num_kids]
    num_seats = problem[:num_seats]
    price     = problem[:price]
    n         = problem[:n]

    @variable(model, 1 <= x[1:n] <= 100, Int)
    @variable(model, 0 <= z <= sum(price.*num_seats), Int)

    scalar_product(model, x,num_seats,:(>=), num_kids)
    @constraint(model, z == sum(x.*price))

    @objective(model,Min,z)

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                println("solution #$sol")
                x_val = convert.(Integer,JuMP.value.(x; result=sol))
                z_val = convert.(Integer,JuMP.value.(z; result=sol))
                println("x:$x_val")
                println("z:$z_val")
                println()

            end
        end
    else
        println("status:$status")
    end

    return status
end

zoo_problem = Dict(
   :num_kids  => 300,
   :num_seats => [30,40],
   :price     => [400,500],
   :n         => 2,

)

@time zoo(zoo_problem)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
using LinearAlgebra;
using SparseArrays;
using IterativeSolvers;
using Krylov;
using LinearOperators;
using LinearMaps;
using QDLDL;
using LDLFactorizations;
using MKLSparse;

function LaLdlInit(vX, mP, vQ, mA, ρ, ρ¹, σ, numElements, numConstraints)

    hDL = ldlt([mP + sparse(σ * I, numElements, numElements) mA'; mA sparse(-ρ¹ * I, numConstraints, numConstraints)]);
    vV  = zeros(numConstraints + numElements); #<! Buffer

    vXX = @view vV[1:numElements];
    vZZ = @view vV[(numElements + 1):end];

    return vXX, vZZ, [hDL, vV];

end

function LaLdl!(tuSolver, vXX, vZZ, vX, mP, vQ, mA, vZ, vY, ρ, ρ¹, σ, numElements, numConstraints, changedΡ)

    if (changedΡ)
        tuSolver[1] = ldlt([mP + sparse(σ * I, numElements, numElements) mA'; mA sparse(-ρ¹ * I, numConstraints, numConstraints)]);
    end

    hDL = tuSolver[1];
    vV  = tuSolver[2];

    @. vXX = σ * vX - vQ;
    @. vZZ = vZ - ρ¹ * vY;
    vV .= hDL \ vV; #<! Like copy!() / copyto!()
    @. vZZ = vZ + ρ¹ * (vZZ - vY);

    return;

end

function QDLdlInit(vX, mP, vQ, mA, ρ, ρ¹, σ, numElements, numConstraints)

    hDL = QDLDL.qdldl([mP + sparse(σ * I, numElements, numElements) mA'; mA sparse(-ρ¹ * I, numConstraints, numConstraints)]);
    vV  = zeros(numConstraints + numElements); #<! Buffer

    vXX = @view vV[1:numElements];
    vZZ = @view vV[(numElements + 1):end];

    return vXX, vZZ, [hDL, vV];

end

function QDLdl!(tuSolver, vXX, vZZ, vX, mP, vQ, mA, vZ, vY, ρ, ρ¹, σ, numElements, numConstraints, changedΡ)

    if (changedΡ)
        tuSolver[1] = QDLDL.qdldl([mP + sparse(σ * I, numElements, numElements) mA'; mA sparse(-ρ¹ * I, numConstraints, numConstraints)]);
    end

    hDL = tuSolver[1];
    vV  = tuSolver[2];

    @. vXX = σ * vX - vQ;
    @. vZZ = vZ - ρ¹ * vY;
    QDLDL.solve!(hDL, vV);
    @. vZZ = vZ + ρ¹ * (vZZ - vY);

    return;

end

function FacLdlInit(vX, mP, vQ, mA, ρ, ρ¹, σ, numElements, numConstraints)

    # Using LDLFactorizations
    hDL = ldl([mP + sparse(σ * I, numElements, numElements) mA'; mA sparse(-ρ¹ * I, numConstraints, numConstraints)]);
    vV  = zeros(numConstraints + numElements); #<! Buffer

    vXX = @view vV[1:numElements];
    vZZ = @view vV[(numElements + 1):end];

    return vXX, vZZ, [hDL, vV];

end

function FacLdl!(tuSolver, vXX, vZZ, vX, mP, vQ, mA, vZ, vY, ρ, ρ¹, σ, numElements, numConstraints, changedΡ)

    if (changedΡ)
        tuSolver[1] = ldl([mP + sparse(σ * I, numElements, numElements) mA'; mA sparse(-ρ¹ * I, numConstraints, numConstraints)]);
    end

    hDL = tuSolver[1];
    vV  = tuSolver[2];

    @. vXX = σ * vX - vQ;
    @. vZZ = vZ - ρ¹ * vY;
    ldiv!(hDL, vV); #<! Like copy!() / copyto!()
    @. vZZ = vZ + ρ¹ * (vZZ - vY);

    return;

end

function ItrSolCgInit(vX, mP, vQ, mA, ρ, ρ¹, σ, numElements, numConstraints)

    mAA = mA' * mA;
    mPI = mP + sparse(σ * I, numElements, numElements);
    mL  = mPI + ρ * mAA;
    vT  = zeros(numElements); #<! Buffer

    vXX = zeros(numElements);
    vZZ = zeros(numConstraints);

    return vXX, vZZ, [mL, mPI, mAA, vT];

end

function ItrSolCg!(tuSolver, vXX, vZZ, vX, mP, vQ, mA, vZ, vY, ρ, ρ¹, σ, numElements, numConstraints, changedΡ; ϵPcg = 1e-6, numItrPcg = 1000)

    if (changedΡ)
        tuSolver[1] = tuSolver[2] + ρ * tuSolver[3];
    end

    mL = tuSolver[1];
    vT = tuSolver[4];

    @. vZZ = ρ * vZ - vY; #<! Using `vZZ` as a buffer
    mul!(vT, mA', vZZ);
    @. vT = σ * vX - vQ + vT;
    IterativeSolvers.cg!(vXX, mL, vT, abstol = ϵPcg, maxiter = numItrPcg);
    # cg!(vXX, mL, σ * vX - vQ + mA' * (ρ * vZ - vY), abstol = ϵPcg, maxiter = numItrPcg);
    mul!(vZZ, mA, vXX);

end

function LinOpCgInit(vX, mP, vQ, mA, ρ, ρ¹, σ, numElements, numConstraints)

    vXX = zeros(numElements);
    vZZ = zeros(numConstraints);

    vT  = zeros(numElements); #<! Buffer

    mL = LinearOperator(Float64, numElements, numElements, true, true, (vU, vW, α, β) -> begin
        mul!(vZZ, mA, vW);
        mul!(vU, mA', vZZ);
        mul!(vU, mP, vW, one(Float64), ρ);
        vU .= vU .+ (σ .* vW);
    end);

    return vXX, vZZ, [mL, vT];

end

function LinOpCg!(tuSolver, vXX, vZZ, vX, mP, vQ, mA, vZ, vY, ρ, ρ¹, σ, numElements, numConstraints, changedΡ; ϵPcg = 1e-6, numItrPcg = 1000)

    if (changedΡ)
        tuSolver[1] = LinearOperator(Float64, numElements, numElements, true, true, (vU, vW, α, β) -> begin
            mul!(vZZ, mA, vW);
            mul!(vU, mA', vZZ);
            mul!(vU, mP, vW, one(Float64), ρ);
            vU .= vU .+ (σ .* vW);
        end);
    end

    mL = tuSolver[1];
    vT = tuSolver[2];

    @. vZZ = ρ * vZ - vY; #<! Using `vZZ` as a buffer
    mul!(vT, mA', vZZ);
    @. vT = σ * vX - vQ + vT;
    IterativeSolvers.cg!(vXX, mL, vT, abstol = ϵPcg, maxiter = numItrPcg);
    # cg!(vXX, mL, σ * vX - vQ + mA' * (ρ * vZ - vY), abstol = ϵPcg, maxiter = numItrPcg);
    mul!(vZZ, mA, vXX);

end

function LinMapsCgInit(vX, mP, vQ, mA, ρ, ρ¹, σ, numElements, numConstraints)

    vXX = zeros(numElements);
    vZZ = zeros(numConstraints);

    vT  = zeros(numElements); #<! Buffer

    mL = LinearMap{Float64}(((vU, vW) -> begin
        mul!(vZZ, mA, vW);
        mul!(vU, mA', vZZ);
        mul!(vU, mP, vW, one(Float64), ρ);
        vU .= vU .+ (σ .* vW);
    end), numElements, numElements; issymmetric = true, isposdef = true, ismutating = true);

    return vXX, vZZ, [mL, vT];

end

function LinMapsCg!(tuSolver, vXX, vZZ, vX, mP, vQ, mA, vZ, vY, ρ, ρ¹, σ, numElements, numConstraints, changedΡ; ϵPcg = 1e-6, numItrPcg = 1000)

    if (changedΡ)
        tuSolver[1] = LinearMap{Float64}(((vU, vW) -> begin
            mul!(vZZ, mA, vW);
            mul!(vU, mA', vZZ);
            mul!(vU, mP, vW, one(Float64), ρ);
            vU .= vU .+ (σ .* vW);
        end), numElements, numElements; issymmetric = true, isposdef = true, ismutating = true);
    end

    mL = tuSolver[1];
    vT = tuSolver[2];

    @. vZZ = ρ * vZ - vY; #<! Using `vZZ` as a buffer
    mul!(vT, mA', vZZ);
    @. vT = σ * vX - vQ + vT;
    IterativeSolvers.cg!(vXX, mL, vT, abstol = ϵPcg, maxiter = numItrPcg);
    # cg!(vXX, mL, σ * vX - vQ + mA' * (ρ * vZ - vY), abstol = ϵPcg, maxiter = numItrPcg);
    mul!(vZZ, mA, vXX);

end

max_stars_repo_name: RoyiAvital/QuadraticProgramSolver

max_issues_repo_name: RoyiAvital/QuadraticProgramSolver

max_forks_repo_name: RoyiAvital/QuadraticProgramSolver


text:
## Functions to pertube a given gene tree
## - change the root position
## - perform a random NNI move
## Necessary file for gene-tree-perturbations.jl
## Claudia (June 2022)

## function to change root of gene tree `tre`
## to a random position
## `s` is the random seed
## out=true (default): we change the root based on a taxon outgroup
## out=false: we change the root to any edge (so whole clades can be outgroup)
function changeRandomRoot!(tre, s::Int; out=true::Bool)
    Random.seed!(s)
    if out
        tips = tipLabels(tre)
        n = length(tips)
        ind = sample(1:n,1)
        rootatnode!(tre,tips[ind][1])
    else
        ind = sample(1:length(tre.edge),1)
        rootonedge!(tre,ind[1], index=true)
    end
end

## function to perform a random NN move
## on the gene tree `tre`
## `s` is the random seed
function randomNNI!(tre, s::Int)
    Random.seed!(s)
    intedge = tre.edge[Bool.([e.istIdentifiable for e in tre.edge])]
    n = length(intedge)
    ind = sample(1:n,1)
    nni!(tre, intedge[ind][1])
end

max_stars_repo_name: gtiley/bad-diamonds

max_issues_repo_name: gtiley/diamond-identifiability

max_forks_repo_name: gtiley/diamond-identifiability


text:
module FloatyBgTileModule

using ..Ahorn, Maple

@mapdef Entity "HonlyHelper/FloatyBgTile" FloatyBgTile(x::Integer, y::Integer, width::Integer=2, height::Integer=2,  tiletype::String="3", disableSpawnOffset::Bool=false)

const placements = Ahorn.PlacementDict(
    "Floaty BG Tiles (HonlyHelper)" => Ahorn.EntityPlacement(
        FloatyBgTile,
        "rectangle",
        Dict{String, Any}()
    )
)

function bgTiletypeEditingOptions()
    validTiles = Ahorn.validTiles("bgTiles")
    tileNames = Ahorn.tileNames("bgTiles")

    return Dict{String, String}(
        tileNames[mat] => string(mat)
        for mat in validTiles if mat != Maple.tile_bg_names["Air"]
    )
end

function drawBgTileEntity(ctx::Ahorn.Cairo.CairoContext, room::Maple.Room, entity::Maple.Entity; material::Union{Char, Nothing}=nothing, alpha::Number=Ahorn.getGlobalAlpha(), blendIn::Bool=false)
    x = Int(get(entity.data, "x", 0))
    y = Int(get(entity.data, "y", 0))

    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    blendIn = get(entity.data, "blendin", blendIn)

    if material === nothing
        key = Ahorn.materialTileTypeKey(entity)
        tile = get(entity.data, key, "3")
        material = isa(tile, Number) ? string(tile) : tile
    end

    # Don't draw air versions, even though they shouldn't exist
    if material[1] in Ahorn.validTileEntityTiles("bgTiles")
        fakeTiles = Ahorn.createFakeTiles(room, x, y, width, height, material[1], blendIn=blendIn)
        objTiles = Ahorn.getObjectTiles(room, x, y, width, height)

        Ahorn.drawFakeTiles(ctx, room, fakeTiles, objTiles, false, x, y, alpha=alpha, clipEdges=true)
    end
end

Ahorn.editingOptions(entity::FloatyBgTile) = Dict{String, Any}(
    "tiletype" => bgTiletypeEditingOptions()
)

Ahorn.minimumSize(entity::FloatyBgTile) = 8, 8
Ahorn.resizable(entity::FloatyBgTile) = true, true

Ahorn.selection(entity::FloatyBgTile) = Ahorn.getEntityRectangle(entity)

Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::FloatyBgTile, room::Maple.Room) = drawBgTileEntity(ctx, room, entity, alpha=0.4)

end

max_stars_repo_name: HollyMagala/HonlyHelper

max_issues_repo_name: HollyMagala/HonlyHelper

max_forks_repo_name: HollyMagala/HonlyHelper


text:
# SimLynx/src/history.jl
# Licensed under the MIT License. See LICENSE.md file in the project root for
# full license information.

abstract type History{T<:Real} end

# Accumulated History

"""
    AccumulatedHistory{T<:Real} <: History{T}

Accumulates history of values over time.

# Fields
- `data::Array{T,1}` - the value data
- `durations::Array{Float64,1}` - the duration data
"""
mutable struct AccumulatedHistory{T<:Real} <: History{T}
    data::Array{T,1}
    durations::Array{Float64,1}
    AccumulatedHistory{T}(x::T) where {T<:Real} =
        new([x], [0.0])
    AccumulatedHistory{T}() where {T<:Real} =
        new([], [])
end

"Update an accumulated history with a new value and duration."
function update!(hist::AccumulatedHistory{T}, x::T, t::Float64) where {T<:Real}
    if t < 0.0
        throw(ArgumentError("weight t cannot be negative, given $t"))
    end
    if isempty(hist.data)
        hist.data = [x]
        hist.durations = [t]
    elseif x == hist.data[end]
        hist.durations[end] += t
    else
        push!(hist.data, x)
        push!(hist.durations, t)
    end
    return nothing
end

"Plot an accumulated history with an optional title and save to file."
function plot_history(hist::AccumulatedHistory;
                      file::Union{String, Nothing} = nothing,
                      title::String = "History")
    println("Displaying $title plot")
    flush(stdout)
    t = 0.0
    x = []
    y = []
    for (value, duration) in zip(hist.data, hist.durations)
        push!(x, t); push!(y, value)
        t += duration
        push!(x, t); push!(y, value)
    end
    p = plot(x,  y, title = title,
             xlabel = "Time", ylabel = "Value",
             legend = false)
    display(p)
    if !isnothing(file)
        println("Writing $title plot to file \"$file\"")
        savefig(p, file)
    end
    return nothing
end

# Tallied History

"""
    TalliedHistory{T<:Real} <: History{T}

Tallies history of values.

# Fields
- `data::Array{T,1}` - the value data
"""
mutable struct TalliedHistory{T<:Real} <: History{T}
    data::Array{T,1}
    TalliedHistory{T}() where {T<:Real} =
        new(Array{T,1}(undef,0))
end

"Update a tallied history with a new value."
function update!(hist::TalliedHistory{T}, x::T) where {T<:Real}
    push!(hist.data, x)
    return nothing
end

"""
Plot a tallied history as a discrete histogram with an optional title and
save to file.
"""
function plot_history(hist::TalliedHistory{T};
                      file::Union{String, Nothing}=nothing,
                      title::String="Histogram") where {T<:Integer}
    low, high = extrema(hist.data) # minimum and maximum data values
    n = high-low+1 # range of data values
    histo = zeros(Int64,n) # empty histogram
    # Construct the histogram counts
    for item in hist.data
        histo[item-low+1] += 1
    end
    # Build the bar chart
    p = bar(range(low, stop=high),histo, title=title,
            xlabel="Value", ylabel="Count",
            legend=false)
    println("Displaying $title plot")
    flush(stdout)
    display(p)
    if !isnothing(file)
        println("Writing $title plot to file \"$file\"")
        savefig(p, file)
    end
    return nothing
end

"""
Plot a tallied history as a continuous histogram with an optional title and
save to file.
"""
function plot_history(hist::TalliedHistory{T};
                      file::Union{String, Nothing}=nothing,
                      title::String="Histogram") where {T<:Real}
    # Build the histogram
    p = histogram(hist.data, bins=100, title=title,
                  xlabel="Value", ylabel="Count",
                  legend=false)
    println("Displaying $title plot")
    flush(stdout)
    display(p)
    if !isnothing(file)
        println("Writing $title plot to file \"$file\"")
        savefig(p, file)
    end
    return nothing
end


text:
using JetPack, Jets, Test

@testset "reghost" begin
    nt = 512
    nx = 256
    G = JopReghost(symspace(JopReghost, Float64, nt, nx), 10.0, 20.0, .004, 10.0) # target depth, receiver depth, time sampling, space sampling
    m = rand(domain(G))
    d = rand(range(G))
    lhs, rhs = dot_product_test(G, m, d)
    @test lhs ≈ rhs
end


text:
const γ = GAMMA_AIR
const p0 = P0
const ρ0 = RHO0
const a0 = A0

"""
    qc2cas(qc)

Calculate calibrated airspeed from ASI (Air Speed indicator), differential
pressure between impact pressure and static pressure.
qc = p_t - p_s

# References

- [1] Ward, D. T. (1993). Introduction to flight test engineering. Elsevier
    Science Ltd. (page 13, formula 2.13)
"""
function qc2cas(qc)
    cas = sqrt((2*γ*p0) / ((γ - 1) * ρ0) * ((qc / p0 + 1) ^ ((γ - 1) / γ) - 1))
    return cas
end

"""
    qc2tas(qc, ρ, p)

Calculate true airspeed from ASI (Air Speed indicator), differential
pressure between impact pressure and static pressure (qc = p_t - p_s), rho
and p.

# References

- [1] Ward, D. T. (1993). Introduction to flight test engineering. Elsevier
    Science Ltd. (page 12, based on formula 2.11)
"""
function qc2tas(qc, ρ, p)
    tas2 = (2*γ*p) / ((γ - 1) * ρ0) * ((qc/p + 1) ^ ((γ - 1) / γ) - 1) * ρ0/ρ
    tas = sqrt(tas2)
    return tas
end

"""
    qc2eas(qc, p)

Calculate equivalent airspeed from ASI (Air Speed indicator), differential
pressure between impact pressure and static pressure (qc = p_t - p_s) and p.

# References

- [1] Ward, D. T. (1993). Introduction to flight test engineering. Elsevier
    Science Ltd.
"""
function qc2eas(qc, p)
    eas = sqrt((2*γ*p) / ((γ - 1) * ρ0) * ((qc/p + 1) ^ ((γ - 1) / γ) - 1))
    return eas
end

"""
    tas2eas(tas, ρ)

Calculate equivalent airspeed from true airspeed and density at current
altitude (ρ).

# References

- [1] Ward, D. T. (1993). Introduction to flight test engineering. Elsevier
    Science Ltd. (page 13, formula 2.15)
"""
function tas2eas(tas, ρ)
    eas = tas * sqrt(ρ / ρ0)
    return eas
end

"""
    eas2tas(qc, ρ)

Calculate true airspeed from equivalent airspeed and density at current
altitude (ρ).

# References

- [1] Ward, D. T. (1993). Introduction to flight test engineering. Elsevier
    Science Ltd. (page 13, formula 2.15)
"""
function eas2tas(eas, ρ)
    tas = eas / sqrt(ρ / ρ0)
    return tas
end

"""
    cas2eas(cas, ρ, p)

Calculate equivalent airspeed from calibrated airspeed, density (ρ) and
pressure (p) at the current altitude.
"""
function cas2eas(cas, ρ, p)
    tas = cas2tas(cas, ρ, p)
    eas = tas2eas(tas, ρ)
    return eas
end

"""
    eas2cas(eas, ρ, p)

Calculate calibrated airspeed from equivalent airspeed, density (ρ) and
pressure (p) at the current altitude.
"""
function eas2cas(eas, ρ, p)
    tas = eas2tas(eas, ρ)
    cas = tas2cas(tas, ρ, p)
    return cas
end

"""
    cas2tas(cas, ρ, p)

Calculate true airspeed from calibrated airspeed, density (ρ) and pressure (p)
at the current altitude.
"""
function cas2tas(cas, ρ, p)

    a = sqrt(γ * p / ρ)

    temp = (cas*cas * (γ - 1.0)/(2.0 * a0*a0) + 1.0) ^ (γ / (γ - 1.0))
    temp = (temp - 1.0) * (p0 / p)
    temp = (temp + 1.0) ^ ((γ - 1.0) / γ) - 1.0

    tas = sqrt(2.0 * a*a / (γ - 1.0) * temp)

    return tas
end

"""
    tas2cas(cas, ρ, p)

Calculate true airspeed from calibrated airspeed, density (ρ) and pressure (p)
at the current altitude.
"""
function tas2cas(tas, ρ, p)

    a = sqrt(γ * p / ρ)

    temp = (tas*tas * (γ - 1.0)/(2.0 * a*a) + 1.0) ^ (γ / (γ - 1.0))
    temp = (temp - 1.0) * (p / p0)
    temp = (temp + 1.0) ^ ((γ - 1.0) / γ) - 1.0

    cas = sqrt(2.0 * a0*a0 / (γ - 1) * temp)

    return cas
end

"""
    tas_alpha_beta_from_uvw(u, v, w)

Calculate true air speed (TAS), angle of attack (α) and angle of side-slip (β)
from aerodynamic velocity expressed in body axis.

# References

- [1] Etkin, B. (2005). Dynamics of atmospheric flight. Dover Publications
    (page 114, formulas 4.3,2 and 4.3,3)
"""
function tas_alpha_beta_from_uvw(u, v, w)

    tas = sqrt(u*u + v*v + w*w)
    alpha = atan(w / u)
    beta = asin(v / tas)
    return [tas, alpha, beta]
end

"""
    incompressible_qinf(tas, ρ)

Calculate incompressible dynamic pressure from true airspeed (tas) and density
(ρ) at current altitude.

# References

- [1] Ward, D. T. (1993). Introduction to flight test engineering. Elsevier
    Science Ltd. (page 13, formula 2.14)
"""
function incompressible_qinf(tas, ρ)
    return 0.5 * ρ * tas*tas
end

"""
    compressible_qinf(tas, p, a)

Calculate compressible dynamic pressure from Mach number and static
pressure (p)

Two different models are used depending on the Mach number:

- Subsonic case: Bernouilli's equation compressible form.
- Supersonic case: to be implemented.

# References

- [1] Ward, D. T. (1993). Introduction to flight test engineering. Elsevier
    Science Ltd. (page 12)
- [2] Fundamentals of Aerdynamics, 5th edition, J.D.Anderson Jr (page 550)
"""
function compressible_qinf(M, p)

    if M < 1
        pt = p * (1 + (γ - 1) / 2 * M*M) ^ (γ / (γ - 1))
    else    # (M >= 1)
        pt = p * ((1 - γ + (2 * γ * M^2))/(γ + 1)) * (((γ + 1)^2 * M^2)/((4 * γ * M^2) - 2*(γ - 1)))^(γ / (γ - 1))
    end
    return pt
end


text:
##### Microphysics models

"""
    compute_cloud_phys!

Define microphysics fields
 - `aux[:CF, k, i]`
 - `aux[:θ_cloudy, k, i]`
 - `aux[:t_cloudy, k, i]`
 - `aux[:q_tot_cloudy, k, i]`
 - `aux[:q_vap_cloudy, k, i]`
 - `aux[:θ_dry, k, i]`
 - `aux[:q_tot_dry, k, i]`

 for all `k` and all `i`
"""
function compute_cloud_phys! end

function compute_cloud_phys!(grid::Grid{FT}, q, aux, params) where {FT}
    gm, en, ud, sd, al = allcombinations(q)
    @unpack param_set = params
    @inbounds for k in over_elems_real(grid)
        q_tot = q[:q_tot, k, en]
        ts = ActiveThermoState(param_set, q, aux, k, en)
        T = air_temperature(ts)
        q_liq = PhasePartition(ts).liq
        q_vap = q_tot - q_liq
        θ = dry_pottemp(ts)
        if q_liq > 0
            aux[:CF, k] = FT(1)
            aux[:θ_cloudy, k] = θ
            aux[:t_cloudy, k] = T
            aux[:q_tot_cloudy, k] = q_tot
            aux[:q_vap_cloudy, k] = q_vap
        else
            aux[:CF, k] = FT(0)
            aux[:θ_dry, k] = θ
            aux[:q_tot_dry, k] = q_tot
        end
    end
end


text:
# Author: Lukasz Bratos

module Blocksys

import SparseArrays
using LinearAlgebra

export readMatrix, readVector, solveGauss, solveGaussWithChoice, lu, solveLu, gaussWithChoice, solveLuWithChoice, printXToFile, printBToFile, computeRightSideVector

"""
Funkcja wczytuje macierz z pliku

Dane wejściowe:
    file - ścieżka pliku wejściowego
Dane wyjściowe:
    krotka (A, n, l), gdzie:
        A - macierz rzadka w formacie SparseMatrixCSC
        n - rozmiar macierzy
        l - rozmiar bloku
"""
function readMatrix(file::String)
    open(file) do f
        firstLine = split(readline(f))
        n = parse(Int64, firstLine[1])
        l = parse(Int64, firstLine[2])

        rows = Int64[]
        columns = Int64[]
        values = Float64[]

        for line in eachline(f)
            line = split(line)
            push!(rows, parse(Int64, line[1]))
            push!(columns, parse(Int64, line[2]))
            push!(values, parse(Float64, line[3]))
        end

        A = SparseArrays.sparse(rows, columns, values)
        return (A, n, l)
    end
end

"""
Funkcja wczytująca wektor z pliku

Dane wejściowe:
    file - ścieżka pliku wejściowego
Dane wyjściowe:
    vector - wczytany wektor
"""
function readVector(file::String)::Vector{Float64}
    open(file) do f
        n = parse(Int64, readline(f))
        vector = []

        for line in eachline(f)
            push!(vector, parse(Float64, line))
        end

        return vector
    end
end

"""
Funkcja zapisująca wektor do pliku

Dane wejściowe:
    file - ścieżka do pliku
    x - wektor, który chcemy zapisać
"""
function printXToFile(file::String, x::Vector{Float64})
    open(file, "w") do f
        foreach(a->write(f, string(a), "\n"), x)
    end
end

"""
Funkcja zapisująca wektor do pliku wraz z błędem względnym

Dane wejściowe:
    file - ścieżka do pliku
    x - wektor, który chcemy zapisać
"""
function printBToFile(file::String, x::Vector{Float64})
    open(file, "w") do f
        y = ones(Float64, length(x))
        write(f, string(norm(y - x) / norm(x)), "\n")
        foreach(a->write(f, string(a), "\n"), x)
    end
end

"""
Funkcja wyliczająca wektor prawych stron na podstawie macierzy A

Dane wejściowe:
    A - macierz rzadka w formacie SparseMatrixCSC
    n - rozmiar macierzy
    l - rozmiar bloku
Dane wyjściowe:
    b - wektor prawych stron
"""
function computeRightSideVector(A::SparseArrays.SparseMatrixCSC{Float64,Int64}, n::Int64, l::Int64)
    x = ones(Float64, n)
    b = zeros(Float64, n)

    for i in 1:n
        for j in max(1, i - (2 + l)):min(n, i + l)
            b[i] += A[i, j]
        end
    end

    return b
end

"""
Funkcja aplikująca metodę eliminacji Gaussa na macierzy A

Dane wejściowe:
    A - macierz rzadka w formacie SparseMatrixCSC
    n - rozmiar macierzy
    l - rozmiar bloku
    b - wektor prawych stron
Dane wyjściowe:
    krotka (A, b), gdzie:
        A - zmodyfikowana macierz
        b - zmodyfikowany wektor prawych stron
"""
function gauss(A::SparseArrays.SparseMatrixCSC{Float64,Int64}, n::Int64, l::Int64, b::Vector{Float64})
    for k in 1:(n - 1)
        for i in (k + 1):min(n, k + l + 1)
            z = A[i, k] / A[k, k]
            A[i, k] = 0.0

            for j in (k + 1):min(n, k + l)
                A[i, j] -= z * A[k, j]
            end

            b[i] -= z * b[k]
        end
    end

    return (A, b)
end

"""
Funkcja rozwiązująca układ równań Ax = b metodą eliminacji Gaussa

Dane wejściowe:
    A - macierz rzadka współczynników w formacie SparseMatrixCSC
    n - rozmiar macierzy
    l - rozmiar bloku
    b - wektor prawych stron
Dane wyjściowe:
    x - wektor wyników
"""
function solveGauss(_A::SparseArrays.SparseMatrixCSC{Float64,Int64}, n::Int64, l::Int64, _b::Vector{Float64})::Vector{Float64}
    A, b = gauss(_A, n, l, _b)

    x = zeros(Float64, n)

    for i in n:-1:1
        sum = 0.0

        for j in (i + 1):min(n, i + l)
            sum += A[i, j] * x[j]
        end

        x[i] = (b[i] - sum) / A[i, i]
    end

    return x
end

"""
Funkcja aplikująca metodę eliminacji Gaussa (z częściowym wyborem elementu głównego) na macierzy A

Dane wejściowe:
    A - macierz rzadka w formacie SparseMatrixCSC
    n - rozmiar macierzy
    l - rozmiar bloku
    b - wektor prawych stron
Dane wyjściowe:
    krotka (A, b, perm), gdzie:
        A - zmodyfikowana macierz
        b - zmodyfikowany wektor prawych stron
        perm - wektor permutacji
"""
function gaussWithChoice(A::SparseArrays.SparseMatrixCSC{Float64,Int64}, n::Int64, l::Int64, b::Vector{Float64})
    perm = [1:n;]

    for k in 1:(n - 1)
        max = 0.0
        maxIdx = 0

        for i in k:min(n, k + l + 1)
            if abs(A[perm[i], k]) > max
                max = abs(A[perm[i], k])
                maxIdx = i
            end
        end

        perm[maxIdx], perm[k] = perm[k], perm[maxIdx]

        for i in (k + 1):min(n, k + l + 1)
            z = A[perm[i], k] / A[perm[k], k]
            A[perm[i], k] = 0.0

            for j in (k + 1):min(n, k + 2 * l)
                A[perm[i], j] -= z * A[perm[k], j]
            end

            b[perm[i]] -= z * b[perm[k]]
        end
    end

    return (A, b, perm)
end

"""
Funkcja rozwiązująca układ równań Ax = b metodą eliminacji Gaussa (z częściowym wyborem elementu głównego)

Dane wejściowe:
    A - macierz rzadka współczynników w formacie SparseMatrixCSC
    n - rozmiar macierzy
    l - rozmiar bloku
    b - wektor prawych stron
Dane wyjściowe:
    x - wektor wyników
"""
function solveGaussWithChoice(_A::SparseArrays.SparseMatrixCSC{Float64,Int64}, n::Int64, l::Int64, _b::Vector{Float64})
    A, b, perm = gaussWithChoice(_A, n, l, _b)

    x = zeros(Float64, n)

    for i in n:-1:1
        sum = 0.0

        for j in (i + 1):min(n, i + 2 * l)
            sum += A[perm[i], j] * x[j]
        end

        x[i] = (b[perm[i]] - sum) / A[perm[i], i]
    end

    return x
end

"""
Funkcja dokonująca rozkładu LU macierzy A metodą eliminacji Gaussa

Dane wejściowe:
    A - macierz rzadka w formacie SparseMatrixCSC
    n - rozmiar macierzy
    l - rozmiar bloku
Dane wyjściowe:
    A - macierz U
    L - macierz L
"""
function lu(A::SparseArrays.SparseMatrixCSC{Float64,Int64}, n::Int64, l::Int64)
    # A is our U
    L = SparseArrays.spzeros(n, n)

    for k in 1:(n - 1)
        L[k ,k] = 1.0

        for i in (k + 1):min(n, k + l + 1)
            z = A[i, k] / A[k, k]
            L[i, k] = z
            A[i, k] = 0.0
            for j in (k + 1):min(n, k + 2 * l)
                A[i, j] -= z * A[k, j]
            end
        end
    end
    L[n,n] = 1.0

    return (A, L)
end

"""
Funkcja rozwiązująca układ równań Ax = b korzystając z rozkładu LU

Dane wejściowe:
    U - macierz U w formacie SparseMatrixCSC
    L - macierz L w formacie SparseMatrixCSC
    n - rozmiar macierzy
    l - rozmiar bloku
    b - wektor prawych stron
Dane wyjściowe:
    x - wektor wyników
"""
function solveLu(U::SparseArrays.SparseMatrixCSC{Float64,Int64}, L::SparseArrays.SparseMatrixCSC{Float64,Int64}, n::Int64, l::Int64, b::Vector{Float64})::Vector{Float64}
    x = zeros(Float64, n)

    for k in 1:(n - 1)
        for i in (k + 1):min(n, k + l + 1)
            b[i] -= L[i, k] * b[k]
        end
    end

    for i in n:-1:1
        sum = 0.0

        for j in (i + 1):min(n, i + l)
            sum += U[i, j] * x[j]
        end

        x[i] = (b[i] - sum) / U[i, i]
    end

    return x
end

"""
Funkcja dokonująca rozkładu LU macierzy A metodą eliminacji Gaussa (z częściowym wyborem elementu głównego)

Dane wejściowe:
    A - macierz rzadka w formacie SparseMatrixCSC
    n - rozmiar macierzy
    l - rozmiar bloku
Dane wyjściowe:
    krotka (A, perm), gdzie:
        A - macierz U
        L - macierz L
        perm - wektor permutacji
"""
function luWithChoice(A::SparseArrays.SparseMatrixCSC{Float64,Int64}, n::Int64, l::Int64)
    # A is our U
    perm = [1:n;]

    L = SparseArrays.spzeros(n, n)

    for k in 1:(n - 1)

        max = 0.0
        maxIdx = 0

        for i in k:min(n, k + l + 1)
            if abs(A[perm[i], k]) > max
                max = abs(A[perm[i], k])
                maxIdx = i
            end
        end

        perm[maxIdx], perm[k] = perm[k], perm[maxIdx]

        for i in (k + 1):min(n, k + l + 1)
            z = A[perm[i], k] / A[perm[k], k]

            L[perm[i], k] = z
            A[perm[i], k] = 0.0

            for j in (k + 1):min(n, k + 2 * l)
                A[perm[i], j] -= z * A[perm[k], j]
            end
        end
    end

    return (A, L, perm)
end

"""
Funkcja rozwiązująca układ równań Ax = b korzystając z rozkładu LU (z częściowym wyborem elementu głównego)

Dane wejściowe:
    A - macierz rozkładu LU w formacie SparseMatrixCSC
    n - rozmiar macierzy
    l - rozmiar bloku
    b - wektor prawych stron
    perm - wektor permutacji
Dane wyjściowe:
    x - wektor wyników
"""
function solveLuWithChoice(U::SparseArrays.SparseMatrixCSC{Float64,Int64}, L::SparseArrays.SparseMatrixCSC{Float64,Int64}, n::Int64, l::Int64, b::Vector{Float64}, perm::Vector{Int64})
    x = zeros(Float64, n)

    for k in 1:(n - 1)
        for i in (k + 1):min(n, k + l + 1)
            b[perm[i]] -= L[perm[i], k] * b[perm[k]]
        end
    end

    for i in n:-1:1
        sum = 0.0

        for j in (i + 1):min(n, i + 2 * l)
            sum += U[perm[i], j] * x[j]
        end

        x[i] = (b[perm[i]] - sum) / U[perm[i], i]
    end

    return x
end

end


text:
using MAT
using Distributed
using Random
include("deeper-cut.config.jl")

struct DataItem
    id::Int
    path::String
    size::Array{Float64}
    joints::Array
    is_single::Bool
    annorect::Dict{String, Any}

    function DataItem(
        id::Int,
    path::String,
    size::Array{Float64},
    joints::Array,
    is_single::Bool,
    annorect::Dict{String, Any}
        )
        return new(id, path, size, joints, is_single, annorect)
    end

    function DataItem(original_data_item::DataItem, new_h_size, new_w_size; shift_idx=true)
        copied_data_item = deepcopy(original_data_item)
        if (shift_idx)
            copied_data_item.joints[:, 1] .+= 1
        end

         scale_h = convert(Float32, new_h_size / copied_data_item.size[2])
    scale_w = convert(Float32,  new_w_size / copied_data_item.size[3])

        copied_data_item.size[2] = new_h_size
        copied_data_item.size[3] = new_w_size

        copied_data_item.joints[:, 2] = round.(Array{Float32}(copied_data_item.joints[:, 2]) .*= scale_w)
         copied_data_item.joints[:, 3] = round.(Array{Float32}(copied_data_item.joints[:, 3]) .*= scale_h)

        # Changing annorect is enough, no need to scale annopoints
        copied_data_item.annorect["x1"] *= scale_w
        copied_data_item.annorect["x2"] *= scale_w
        copied_data_item.annorect["y1"] *= scale_h
        copied_data_item.annorect["y2"] *= scale_h

       return copied_data_item
    end
end

function mirror_data_item(data_item)
    id = data_item.id
    path = data_item.path
    size = data_item.size
    image_width = size[3]
    # mirror joints coords
    mirrored_joints = deepcopy(data_item.joints)
    mirrored_joints[:,2] = data_item.size[3] .- data_item.joints[:,2] .|> Int32
    is_single = data_item.is_single
    # mirror joint ids
    for i in enumerate(mirrored_joints[:, 1])
       mirrored_joints[i[1], 1] = symmetric_joint_dict[i[2]]
    end
    # mirror annorect
    annorect = data_item.annorect |> deepcopy
    annorect["x1"]  =   data_item.size[3] .- data_item.annorect["x1"]
    annorect["x2"]  =   data_item.size[3] .- data_item.annorect["x2"]
    return DataItem(id, path, size, mirrored_joints, is_single, annorect)
end

function read_cropped_mpii_annotations(;should_shuffle=false)
    file = matopen(path_to_processed_mat)
    dataset = read(file, "dataset")
    close(file)
    reformat_dataset!(dataset, should_shuffle)
    return dataset
end

function reformat_dataset!(dataset, should_shuffle)
    zipped = zip(dataset["image"], dataset["size"], dataset["joints"])
    collected = collect(zipped)
    if should_shuffle
        shuffle!(collected)
    end
    vecced = vec(collected)
    dataset["size"] = map(e -> e[2],vecced)
    dataset["image"] = map(e -> e[1],vecced)
    dataset["joints"] = map(e -> e[3],vecced)

    image_metas = find_annorect_and_is_single_info(dataset["image"]);
    dataset["is_single"] = image_metas[:, 1]
    dataset["annorect"] = image_metas[:, 2]
    return dataset
end

function find_annorect_and_is_single_info(dataset_images)
    file = matopen(path_to_multi_person_mat);
    multi_person = read(file);
    close(file);

    file = matopen(path_to_single_person_mat);
    single_person = read(file);
    close(file);

    single_image_paths = map(x -> x["name"], single_person["annolist"]["image"])
    single_image_paths = vec(single_image_paths);

    multi_image_paths = map(x -> x["name"], multi_person["annolist"]["image"])
    multi_image_paths = vec(multi_image_paths);

    image_metas = map(image_path -> find_image_idx_in_single_multi_meta(image_path, single_image_paths, multi_image_paths, single_person, multi_person) ,dataset_images)

    vcatted_image_metas = vcat(image_metas...);
    return vcatted_image_metas
end

function find_image_idx_in_single_multi_meta(image_name, single_paths, multi_paths, single_person, multi_person)
    is_single = true;
    found_idx = findlast(x -> x == image_name ,single_paths)
    if found_idx == nothing
        found_idx = findlast(x -> x == image_name, multi_paths)
        is_single = false
    elseif found_idx == nothing
        throw(DomainError(image_name, "image path not found in netiher single nor multi"))
    end

    if is_single
        annorect = single_person["annolist"]["annorect"][1, found_idx]
    else
        annorect = multi_person["annolist"]["annorect"][1, found_idx]
    end

    if annorect == nothing
     throw(DomainError(image_name, "image path has no annorect"))
    end

    return reshape([is_single, annorect], (1 ,2))
end

function raw_data_to_data_item(indexed_raw_data)
    i = indexed_raw_data[1]
    img_path = indexed_raw_data[2][1]
    img_size = indexed_raw_data[2][2]
    joints_data = indexed_raw_data[2][3][:][1]
    joints_ids = joints_data[:, 1]
    id_check = all(id -> id < global_num_joints, joints_ids)
    if (!id_check)
        throw(DomainError(joints_data, "invalid joint id"))
    end
    is_single = indexed_raw_data[2][4]
    annorect = indexed_raw_data[2][5]
    return DataItem(i, img_path, img_size, joints_data, is_single, annorect)
end

function get_from_dataset(dataset, initial_index=1, last_index=max_image_number_to_read; should_use_pmap = true)
    sizes = dataset["size"][initial_index: last_index]
    images = dataset["image"][initial_index: last_index]
    joints = dataset["joints"][initial_index: last_index]
    is_singles = dataset["is_single"][initial_index: last_index]
    annorects = dataset["annorect"][initial_index: last_index]
    if should_use_pmap
            data_items = pmap(raw_data_to_data_item, enumerate(zip(images, sizes, joints, is_singles, annorects)))
    else
            data_items = map(raw_data_to_data_item, enumerate(zip(images, sizes, joints, is_singles, annorects)))
    end
    return data_items
end


text:
import Base: -, ==, abs, isless

mutable struct Interval
    steps::Int
end

-(interval::Interval) = Interval(-interval.steps)
==(a::Interval, b::Interval) = a.steps == b.steps
abs(interval::Interval) = Interval(abs(interval.steps))
isless(a::Interval, b::Interval) = a.steps < b.steps

const ROOT = Interval(0)
const OCTAVE = Interval(12)

const MINOR_THIRD = Interval(3)
const MAJOR_THIRD = Interval(4)
const PERFECT_FOURTH = Interval(5)
const PERFECT_FIFTH = Interval(7)
const MINOR_SIXTH = Interval(8)
const MAJOR_SIXTH = Interval(9)
const MINOR_SEVENTH = Interval(10)
const MAJOR_SEVENTH = Interval(11)


text:
# ToDo: Flush out the remainder of the outer join, including checking the reserve for hits.
# Can we safely copy into consign all the time and then just pop from reserve on first hit?

# Alternating march in until there is a hit, then send through done to next for march out
function Base.start{I, T}(iterator::Join{Outer, I, T})
	reserve = Dict{T, _JoinValue}()
	consign = Dict{T, _JoinValue}()
	sizehint!(reserve, length(iterator.small))
	for i = 1:iterator.smalllength

		# March in from the small source
		smallvalue = get!(reserve, iterator.small(i), _JoinValue(Vector{Int64}(0), Vector{Int64}(0)))
		push!(smallvalue.small, i)

		# March in from the large source
		largevalue = get!(reserve, iterator.large(i), _JoinValue(Vector{Int64}(0), Vector{Int64}(0)))
		push!(largevalue.large, i)

		# Return the hit on moving both to consign
		if iterator.small[i] == iterator.large[i] && length(smallvalue.large) > 0 && length(largevalue.small) > 0
			consign[iterator.small[i]] = pop!(reserve, iterator.small[i])
			return _JoinState(reserve, consign, i, i, length(largevalue.small), length(smallvalue.large))
		elseif length(smallvalue.large) > 0 && length(largevalue.small) > 0
			consign[iterator.small[i]] = pop!(reserve, iterator.small[i])
			consign[iterator.large[i]] = pop!(reserve, iterator.large[i])
			return _JoinState(reserve, consign, i, i, length(largevalue.small), length(smallvalue.large))
		elseif length(largevalue.small) > 0
			consign[iterator.large[i]] = pop!(reserve, iterator.large[i])
			return _JoinState(reserve, consign, i, i, length(largevalue.small), length(smallvalue.large))
		elseif length(smallvalue.large) > 0
			consign[iterator.small[i]] = pop!(reserve, iterator.small[i])
		end
	end

	# There were no hits so it is safe to tell next to just loop through both sources.
	# ToDo: patch this state so that is actually does this
	_JoinState(reserve, consign, length(iterator.large) + 1, length(iterator.small) + 1, 0, 0)

	#=
	# Both sources are the same size, so we start sending out all the data
	if length(iterator.small) == length(iterator.large)
		return _JoinState(reserve, consign, length(iterator.large) + 1, length(iterator.small) + 1, 0, 0)
	end

	# Still data left in the large source, so we start se
	for i = (length(iterator.small) + 1):length(iterator.large)
		largevalue = get(reserve, iterator.large(i), _JoinValue(Vector{Int64}(0), Vector{Int64}(0)))
		if length(largevalue.small) > 0
			return _JoinState(reserve, consign, i, length(iterator.small) + 1, length(largevalue.small), 0)
		end
	end

	# No matches found
	_JoinState(reserve, consign, length(iterator.large) + 1, length(iterator.small) + 1, 0, 0)
	=#
end

# March out the finds, small then big, until there are none, and then alternatingly march
# in until there is another hit
function Base.next{I, T}(iterator::Join{Outer, I}, state::_JoinState{T})

	# We are in next so there must be something to return, start by marching out the small
	# finds, on the left. All matches are always in consign.
	if state.smallmarchout > 0 && iterator.parity
		item = (iterator.small[state.smallmarchin], state.smallmarchin, state.consign[iterator.small[state.smallmarchin]].large[state.smallmarchout])
		largemarchout = state.largemarchout
		smallmarchout = state.smallmarchout - 1

	# March out the small finds, on the right
	elseif state.smallmarchout > 0
		item = (iterator.small[state.smallmarchin], state.consign[iterator.small[state.smallmarchin]].large[state.smallmarchout], state.smallmarchin)
		largemarchout = state.largemarchout
		smallmarchout = state.smallmarchout - 1

	# March out the large finds, on the right
	elseif state.largemarchout > 0 && iterator.parity
		ritem = (iterator.large[state.largemarchin], state.consign[iterator.large[state.largemarchin]].small[state.largemarchout], state.largemarchin)
		largemarchout = state.largemarchout - 1
		smallmarchout = state.smallmarchout

	# March out the large finds, on the left
	elseif state.largemarchout > 0
		item = (iterator.large[state.largemarchin], state.largemarchin, state.consign[iterator.large[state.largemarchin]].small[state.largemarchout])
		largemarchout = state.largemarchout
		smallmarchout = state.smallmarchout - 1

	# Opps something went wrong, maybe next was called out of loop, or sequence
	else
		error("next() was called on an invalid state")
	end

	# More left to march out after current, safe to loop through done to next
	if largemarchout > 0 || smallmarchout > 0
		return (item, _JoinState(state.reserve, state.consign, state.largemarchin, state.smallmarchin, largemarchout, smallmarchout))
	end

	# Nothing left to march out after the current item, so proceed with alternating march in
	if state.smallmarchin < length(iterator.small)
		for i = (state.smallmarchin + 1):length(iterator.small)

			# March in from the small source
			smallvalue = get(state.consign, iterator.small(i), _JoinValue(Vector{Int64}(0), Vector{Int64}(0)))
			if length(smallvalue.large) < 1
				smallvalue = get!(state.reserve, iterator.small(i), smallvalue)
			end
			push!(smallvalue.small, i)

			# March in from the large source
			largevalue = get(state.consign, iterator.large(i), _JoinValue(Vector{Int64}(0), Vector{Int64}(0)))
			if length(largevalue.small) < 1 && length(largevalue.large) < 1
				largevalue = get!(state.reserve, iterator.large(i), largevalue)
			end
			push!(largevalue.large, i)

			# Return the hit
			if length(smallvalue.large) > 0 || length(largevalue.small) > 0
				return (item, _JoinState(state.reserve, state.consign, i, i, length(largevalue.small), length(smallvalue.large)))
			end
		end
	end

	# Both sources are the same size, so we are done
	if length(iterator.small) == length(iterator.large)
		return (item, _JoinState(state.reserve, state.consign, length(iterator.small) + 1, length(iterator.large) + 1, 0, 0))
	end

	# Still data left in the large source, we only need roll call, without marching in
	for i = (length(iterator.small) + 1):length(iterator.large)
		largevalue = get(state.reserve, iterator.large(i), _JoinValue(Vector{Int64}(0), Vector{Int64}(0)))
		if length(largevalue.small) > 0
			return (item, _JoinState(state.reserve, state.consign, i, length(iterator.small) + 1, length(largevalue.small), 0))
		end
	end

	# Whew! made it all the way through and found no more matches than the current one.
	(item, _JoinState(state.reserve, state.consign, length(iterator.large) + 1, length(iterator.small) + 1, 0, 0))
end


text:
module TodosChannel

using Channels, App
@dependencies

function subscribe()
  Channels.subscribe(wsclient(@params), "/todos")
  "OK"
end

end

max_stars_repo_name: essenciary/genie-todo-mvc

max_issues_repo_name: essenciary/genie-todo-mvc

max_forks_repo_name: essenciary/genie-todo-mvc


text:
module Microbiome

export
    # Functions
    ## abundance
    abundancetable,
    filterabund,
    relativeabundance,
    relativeabundance!,
    rownormalize,
    rownormalize!,
    colnormalize,
    colnormalize!,
    nfeatures,
    getfeature,
    featurenames,
    featuretotals,
    nsamples,
    getsample,
    samplenames,
    sampletotals,
    ## Diversity
    ginisimpson,
    shannon,
    present,
    prevalence

using Reexport
@reexport using SpatialEcology

using Statistics
using StatsBase
using DataFrames

import Base: getindex, setindex, length

include("ecotranslations.jl")
include("abundances.jl")
include("distances.jl")

end  # module Microbiome