text:
using Test
using Logging
using NetworkOptions
using NetworkOptions: BUNDLED_CA_ROOTS

const TEST_URLS = [
    "" # not a valid host name
    "com"
    ".com" # not a valid host name
    "example.com"
    "user@example.com"
    "https://example.com"
    "xyz://example.com" # protocol doesn't matter
    "\x02H2ݼtOٲ\0RƆW9\e]B>#ǲR" # not a valid host name
]

const TRANSPORTS = [nothing, "ssl", "tls", "ssh", "xyz"]

host_variants(host::AbstractString) = Dict(
    "$host"              => true,
    ".$host"             => false,
    "$host."             => false,
    "/$host"             => false,
    "$host/"             => false,
    "user@$host"         => true,
    "user@$host:path"    => true,
    "user@$host:/path"   => true,
    "user@$host/"        => false,
    "user@$host/:"       => false,
    "user@$host/path"    => false,
    "user@$host/:path"   => false,
    "user@$host."        => false,
    "user@.$host"        => false,
    "@$host"             => false,
    "$host@"             => false,
    "://$host"           => false,
    "https://$host"      => true,
    "https://$host/"     => true,
    "https://$host/path" => true,
    "https://$host."     => false,
    "https://.$host"     => false,
    "https:/$host"       => false,
    "https:///$host"     => false,
    "xyz://$host"        => true,
    "xyz://$host."       => false,
    "xyz://.$host"       => false,
)

function clear_vars!(ENV)
    delete!(ENV, "JULIA_NO_VERIFY_HOSTS")
    delete!(ENV, "JULIA_SSL_NO_VERIFY_HOSTS")
    delete!(ENV, "JULIA_SSH_NO_VERIFY_HOSTS")
end

function without_warnings(body::Function)
    log_level = Logging.min_enabled_level(current_logger())
    disable_logging(Logging.Warn)
    try body()
    finally
        disable_logging(log_level-1)
    end
end


text:
type Decoder
    decoderFunctions

    Decoder() = new(Dict{ASCIIString,Function}(
                        "_"  => (x -> nothing),
                        ":"  => symbol,
                        "\$" => (x -> TSymbol(x)),
                        "?"  => (x -> x == "t"),
                        "b"  => base64decode,
                        "c"  => (x -> x[1]),
                        "i"  => (x -> Base.parse(Int64, x)),
                        "d"  => (x -> Base.parse(Float64, x)),
                        "f"  => (x -> decimal(lowercase(x))),
                        "r"  => (x -> TURI(x)),
                        "n"  => (x -> Base.parse(BigInt, x)),
                        "u"  => (x -> Base.Random.UUID(x)), # only string case so far
                        "t"  => parsedatetime,
                        "m"  => (x -> Dates.unix2datetime(Base.parse(Float64, x) / 1000.)), # maybe not sufficient
                        "z"  => (x -> if (x == "NaN")
                                          NaN
                                      elseif (x == "INF")
                                          Inf
                                      elseif (x == "-INF")
                                          -Inf
                                      else
                                          throw(string("Don't know how to encode: ", x))
                                      end),

                        # tag decoders
                        "'"  => (x -> x),
                        "set" => (x -> TSet(x)),
                        "link" => (x -> Link(x...)),
                        "list" => tolist,
                        "ratio" => (x -> x[1]//x[2]),
                        "cmap" => (x -> [a[1] => a[2]
                                         for a in zip(x[1:2:end], x[2:2:end])])
                    ))
end

function getindex(d::Decoder, k::AbstractString)
    d.decoderFunctions[k]
end

function add_decoder(e::Decoder, tag::AbstractString, f::Function)
    e.decoderFunctions[tag] = f
end

function decode(e::Decoder, node::Any, cache::Cache=RollingCache(), as_map_key::Bool=false)
    decode_value(e, node, cache, as_map_key)
end

function decode_value(e::Decoder, node::Any, cache::Cache, as_map_key::Bool=false)
    node
end

function decode_value(e::Decoder, node::Bool)
    node ? true : false # where we may have to add TTrue, TFalse for set issue
end

function decode_value(e::Decoder, node::Array{Any,1}, cache::Cache, as_map_key::Bool=false)
    if !isempty(node)
        if node[1] == MAP_AS_ARRAY
            returned_dict = Dict()
            for i in 2:2:length(node)
                key = decode_value(e, node[i], cache, true)
                value = decode_value(e, node[i+1], cache, as_map_key)
                returned_dict[key] = value
            end
            return returned_dict
        else
            decoded = decode_value(e, node[1], cache, as_map_key)
            if isa(decoded, Tag)
                return decode_value(e, decoded, node[2], cache, as_map_key)
            end
        end
    end

    [decode_value(e, x, cache, as_map_key) for x in node]
end

function decode_value(e::Decoder, hash::Dict, cache::Cache, as_map_key::Bool=false)

    if length(hash) != 1
        h = Dict{Any,Any}()
        for kv in hash
            key = decode_value(e, kv[1], cache, true)
            val = decode_value(e, kv[2], cache, false)
            h[key] = val
        end
        return h
    else
        for (k,v) in hash
            key = decode_value(e, k, cache, true)
            if isa(key, Tag)
                return decode_value(e, key, v, cache, as_map_key)
            end
            return Dict{Any,Any}(key => decode_value(e, v, cache, false))
        end
    end
end

function decode_value(e::Decoder, s::AbstractString, cache::Cache, as_map_key::Bool=false)
    if iscachekey(s)
        return decode_value(e, read(cache, s), cache, as_map_key)
    end

    if iscacheable(s, as_map_key)
        write!(cache, s)
    end

    if !startswith(s, ESC)
        s
    elseif startswith(s, TAG)
        Tag(s[3:end])
    elseif startswith(s, ESC_ESC) ||  startswith(s, ESC_SUB) || startswith(s, ESC_RES)
        s[2:end]
    elseif startswith(s, ESC)
        #2:2 is necessary to get str instead of char
        tag = s[2:2]
        if haskey(e.decoderFunctions, tag)
            e[tag](s[3:end])
        else
	    TaggedValue(tag, s[3:end])
        end
    end
end

function decode_value(e::Decoder, tag::Tag, value, cache::Cache, as_map_key::Bool=false)
    if haskey(e.decoderFunctions, tag.rep)
        e[tag.rep](decode_value(e, value, cache, false))
    else
        TaggedValue(tag.rep, decode_value(e, value, cache, false))
    end
end

ext: jl


text:
module ManifestUtilities

import LightGraphs
import MetaGraphs
import Pkg
import TOML

export prune_manifest

include("types.jl")

include("pkg.jl")

include("prune_manifest.jl")

end # module


text:
```
ДАНО: Робот - в произвольной клетке ограниченного прямоугольного поля без внутренних перегородок

РЕЗУЛЬТАТ: Робот - в исходном положении, в клетке с роботом стоит маркер, и все остальные клетки поля промаркированы в шахматном порядке

```

using HorizonSideRobots
originRobot = Robot(10,7,animate = false)
include("../functions.jl")

function mod_moveToBorder!(robot::Robot, side::HorizonSide, counter::Integer)
    while true
        if !isborder(robot, side)
            move!(robot, side)
            if (counter % 2 == 0)
                putmarker!(robot)
            end
            counter += 1
            println(counter)
        else
            break
        end
    end
    return counter
end

function mod_moveLikeSnake!(robot, moveSides::NTuple{2,HorizonSide}, borderSide::HorizonSide)::Nothing
    border_borderSide = false
    counter = 0
    while !border_borderSide
        for side in moveSides
            counter = mod_moveToBorder!(robot, side, counter)
            if !(isborder(robot, borderSide))
                move!(robot, borderSide)
                if (counter % 2 == 0)
                    putmarker!(robot)
                end
                counter += 1
                println(counter)
            else
                border_borderSide = true
                break
            end
        end
    end
end

#=function get_one_marked_line!(robot::Robot, side, counter::Integer)::Nothing
    while (!(isborder(robot,side)))
        if (counter % 2 == 1)
            putmarker!(robot)
            move!(robot,side)
        else
            move!(robot,side)
        end
        counter += 1
    end
end=#

# master! главная функция программы.
function master7!(robot::Robot)::Nothing
    back_path = moving_in_angle!(robot)

    mod_moveLikeSnake!(robot, (Nord, Sud), Ost)

    moving_back_to_start!(robot, back_path)
end

master7!(originRobot)
show!(originRobot)

max_stars_repo_name: Igetsuga/AlgorithmicsHomework

max_issues_repo_name: Igetsuga/AlgorithmicsHomework

max_forks_repo_name: Igetsuga/AlgorithmicsHomework


text:
using WAVI, Test

@testset "Version Updates" begin
    @info "testing the current version of WAVI against previous version...."

    function version_update_test()
        #Grid and boundary conditions
        nx = 80
        ny = 10
        nσ = 4
        x0 = 0.0
        y0 = -40000.0
        dx = 8000.0
        dy = 8000.0
        h_mask=trues(nx,ny)
        u_iszero = falses(nx+1,ny); u_iszero[1,:].=true
        v_iszero=falses(nx,ny+1); v_iszero[:,1].=true; v_iszero[:,end].=true

        #alternative bc notations
        u_iszero = ["north"]
        v_iszero = ["west", "east"]
        grid = Grid(nx = nx,
                    ny = ny,
                    nσ = nσ,
                    x0 = x0,
                    y0 = y0,
                    dx = dx,
                    dy = dy,
                    h_mask = h_mask,
                    u_iszero = u_iszero,
                    v_iszero = v_iszero)

        #Bed
        bed = WAVI.mismip_plus_bed #function definition

        #solver parameters
        maxiter_picard = 1
        solver_params = SolverParams(maxiter_picard = maxiter_picard)

        #Physical parameters
        default_thickness = 100.0 #set the initial condition this way
        accumulation_rate = 0.3
        params = Params(default_thickness = default_thickness,
                        accumulation_rate = accumulation_rate)

        #make the model
        model = Model(grid = grid,
                        bed_elevation = bed,
                        params = params,
                        solver_params = solver_params)

        #timestepping parameters
        niter0 = 0
        dt = 0.1
        end_time = 100.
        timestepping_params = TimesteppingParams(niter0 = niter0,
                                                dt = dt,
                                                end_time = end_time)

        #nb no output parameters

        simulation = Simulation(model = model,
                            timestepping_params = timestepping_params)

        #perform the simulation
        run_simulation!(simulation)
        return simulation
    end

    simulation = version_update_test();

    if  VERSION == v"1.6.2"
        filename = joinpath(dirname(@__FILE__), "v1_6_2_MISMIP_100yr_output_8kmres_maxiter1_timesteppt1.jld2")
        example_output = load(filename)
    elseif VERSION == v"1.6.1"
        filename = joinpath(dirname(@__FILE__), "v1_6_1_MISMIP_100yr_output_8kmres_maxiter1_timesteppt1.jld2")
        example_output = load(filename)
    elseif  VERSION == v"1.5.1"
        filename = joinpath(dirname(@__FILE__), "v1_5_1_MISMIP_100yr_output_8kmres_maxiter1_timesteppt1.jld2")
        example_output = load(filename)
    else
    example_output = Dict("h" => NaN, "u" => NaN, "v" => NaN, "viscosity" => NaN, "grounded_fraction" => NaN, "bed_speed" => NaN)
    end
    @test simulation.model.fields.gh.h == example_output["h"]
    @test simulation.model.fields.gu.u == example_output["u"]
    @test simulation.model.fields.gv.v == example_output["v"]
    @test simulation.model.fields.gh.ηav == example_output["viscosity"]
    @test simulation.model.fields.gh.grounded_fraction == example_output["grounded_fraction"]
    @test simulation.model.fields.gh.bed_speed == example_output["bed_speed"]
end


text:
#==============================================================================#
#
#    Structs responsible for determining the restrictions imposed on the
#    state space of a diffusion process
#
#==============================================================================#

"""
    UnboundedStateSpace <: DiffusionStateSpace

No restrictions imposed on the state-space of the process (i.e. ℝᵈ)
"""
struct UnboundedStateSpace <: DiffusionStateSpace end

"""
    LowerBoundedStateSpace{T,S,N} <: DiffusionStateSpace

Lower bounds imposed on the state-space of a diffusion process. `T` is used
to list the indices that have lower-bound restrictions, `S` indicates the values
of the lower-bounds, `N` is the total number of coordinates with lower-bound
restrictions
"""
struct LowerBoundedStateSpace{T,S,N} <: DiffusionStateSpace

    function LowerBoundedStateSpace(
            coords::NTuple{N,Integer},
            bounds::NTuple{N,Number},
        ) where N
        new{coords, bounds, N}()
    end

    function LowerBoundedStateSpace(
            coords,
            bounds,
        )
        N = length(bounds)
        @assert length(coords) == N
        @assert all( map(c->(typeof(c)<:Integer), coords) )
        @assert all( map(b->(typeof(b)<:Number), bounds) )
        new{Tuple(coords), Tuple(bounds), N}()
    end

end

bound_info(::LowerBoundedStateSpace{T,S,N}) where {T,S,N} = T,S,N

"""
    UpperBoundedStateSpace{T,S,N} <: DiffusionStateSpace

Upper bounds imposed on the state-space of a diffusion process. `T` is used
to list the indices that have upper-bound restrictions, `S` indicates the values
of the upper-bounds, `N` is the total number of coordinates with upper-bound
restrictions
"""
struct UpperBoundedStateSpace{T,S,N} <: DiffusionStateSpace
    function UpperBoundedStateSpace(coords, bounds)
        T,S,N = bound_info(LowerBoundedStateSpace(coords, bounds))
        new{T,S,N}()
    end
end

"""
    BoundedStateSpace{L,U} <: DiffusionStateSpace

Upper and lower bounds imposed on the state-space of a diffusion process.
`L` corresponds to lower bounds, `U` corresponds to upper bounds.
"""
struct BoundedStateSpace{L,U} <: DiffusionStateSpace
    function BoundedStateSpace(
            (coords_lower, bounds_lower),
            (coords_upper, bounds_upper)
        )
        L = LowerBoundedStateSpace(coords_lower, bounds_lower)
        U = UpperBoundedStateSpace(coords_upper, bounds_upper)
        new{L,U}()
    end

    function BoundedStateSpace(
        L::LowerBoundedStateSpace,
        U::UpperBoundedStateSpace
        )
        new{L,U}()
    end
end

"""
    bound_satisfied(::UnboundedStateSpace, x)

No restrictions, bounds satisfied by default
"""
@inline _bound_satisfied(::UnboundedStateSpace, x) = true

"""
    bound_satisfied(::LowerBoundedStateSpace{T,S,N}, x) where {T,S,N}

Checks if all coordinates adhere to lower bound restrictions
"""
@generated function _bound_satisfied(
        ::LowerBoundedStateSpace{T,S,N},
        x
    ) where {T,S,N}
    ex = :(true)
    for i = 1:N
        ex = :(x[T[$i]] > S[$i] ? $ex : false)
    end
    return ex
end

"""
    bound_satisfied(::UpperBoundedStateSpace{T,S,N}, x) where {T,S,N}

Checks if all coordinates adhere to upper bound restrictions
"""
@generated function _bound_satisfied(
        ::UpperBoundedStateSpace{T,S,N},
        x
    ) where {T,S,N}
    ex = :(true)
    for i = 1:N
        ex = :(x[T[$i]] < S[$i] ? $ex : false)
    end
    return ex
end

"""
    bound_satisfied(::BoundedStateSpace{L,U}, x) where {L,U}

Checks if all coordinates adhere to lower and upper bound restrictions
"""
function _bound_satisfied(::BoundedStateSpace{L,U}, x) where {L,U}
    _bound_satisfied(L, x) && _bound_satisfied(U, x)
end

function bound_info(::DiffusionProcess{T,DP,DW,SS}) where {T,DP,DW,SS}
    bound_info(SS)
end

bound_satisfied(P::DiffusionProcess, x) = _bound_satisfied(P, x)

function _bound_satisfied(::DiffusionProcess{T,DP,DW,SS}, x) where {T,DP,DW,SS}
    _bound_satisfied(SS, x)
end


text:
module JSExpr

export @js, @js_str
export JSString

using JSON

include("./jsstring.jl")
include("./ast.jl")

"""
    crawl(expr)

Crawl a given Julia expression and convert it into a `JSNode`.

# Examples
```julia-repl
julia> JSExpr.crawl(:(foo = "bar"))
:(JSAST(:(=), JSTerminal(:foo), JSTerminal("bar")))
```
"""
function crawl(ex::Expr)::Expr
    crawl(Val(ex.head), ex.args...)
end

"""
    crawl(Val(head), args...)::Expr

Enables multiple dispatch on expressions using `Val` types.
The expectation is that each dispatched crawl function returns an expression
that yields a `JSNode` by calling the crawl-function recursively on nested
expressions.

# Examples
```julia
function crawl(::Val{:+}, lhs, rhs)
    :(JSAST(:+, \$(crawl(lhs)), \$(crawl(rhs))))
end
```
"""
function crawl(::Val{T}, args...)::Expr where {T}
    error("Expression type ($(QuoteNode(T))) not supported.")
end

"""
    deparse(jsnode)

Convert a `JSNode` to `JSString`.
"""
function deparse(ex::JSAST)::JSString
    deparse(Val(ex.head), ex.args...)::JSString
end

"""
    deparse(Val(head), args...)

Transform a JSAST with the specified `head` and `args` into a `JSString`.
The expectation is that each dispatched deparse function returns a bare
`JSString` literal, formed by appropriate ordering and concatenation of the
output of recursive calls to the `deparse` function.
"""
function deparse(::Val{H}, args...) where {H}
    # This should only happen if there is an asymmetry between crawl and
    # deparse; for example, if crawl is implemented but deparse was forgotten,
    # or crawl generates a JSAST with more arguments than the deparse method
    # supports.
    error("JSAST expression type ($(QuoteNode(H))) cannot be deparsed.")
end

"""
    @js(ex)

A macro to convert a Julia expression into a `JSString`.
"""
macro js(ex)
    return Expr(:call, :deparse, crawl(ex))
end

"""
    @crawl(ex)

A macro to generate a `JSNode` for the given expression.
This is useful for making assertions about the generated AST structure, which
is in turn useful for testing, but most users should just use `@js`.
"""
macro crawl(ex)
    return crawl(ex)
end

crawl(::LineNumberNode) = :(nothing)

# For QuoteNode's, if they wrap symbols (and they usually do), we implicitly
# convert them into strings. This is what we want in a lot of places, but there
# are a few places where we need to be careful about this (namely, while
# translating dot expressions).
function crawl(node::QuoteNode)
    if isa(node.value, Symbol)
        return crawl(string(node.value))
    end
    return crawl(node.value)
end

# All other terminals
crawl(ex::T) where {T} = :(JSTerminal($(esc(ex))))

include("./literals.jl")
include("./control.jl")
include("./call.jl")
include("./infix.jl")
include("./macrocall.jl")
include("./arrays.jl")
include("./objects.jl")
include("./jskeywords.jl")
include("./interpolation.jl")
include("./function.jl")
include("./ref.jl")
include("./juliaisms.jl")
include("./compat.jl")

end


text:
import Avalon.walkstruct

conv3x3(inplanes, outplanes; stride=1, dilation=1) =
    Conv2d(inplanes, outplanes, 3; stride=stride, padding=dilation, dilation=dilation, bias=false)

conv1x1(inplanes, outplanes; stride=1) =
    Conv2d(inplanes, outplanes, 1; stride=stride, bias=false)

################################################################################
#                                BasicBlock                                    #
################################################################################

mutable struct BasicBlock
    conv1::Conv2d       # 3x3
    bn1::BatchNorm2d
    conv2::Conv2d       # 3x3
    bn2::BatchNorm2d
    downsample::Union{Sequential, Nothing}
end

function BasicBlock(inplanes::Int, planes::Int; downsample=nothing, stride::Int=1,
                    base_width::Int=64, dilation=1)
    return BasicBlock(
        conv3x3(inplanes, planes; stride=stride),
        BatchNorm2d(planes),
        conv3x3(planes, planes),
        BatchNorm2d(planes),
        downsample
    )
end

expansion(m::Type{BasicBlock}) = 1

function (m::BasicBlock)(x::AbstractArray{T, 4}) where T
    identity = x
    out = m.conv1(x)
    out = m.bn1(out)
    out = relu.(out)
    out = m.conv2(out)
    out = m.bn2(out)
    if m.downsample != nothing
        identity = m.downsample(x)
    end
    out = out .+ identity
    out = relu.(out)
    return out
end

################################################################################
#                                Bottleneck                                    #
################################################################################

mutable struct Bottleneck
    conv1::Conv2d
    bn1::BatchNorm2d
    conv2::Conv2d
    bn2::BatchNorm2d
    conv3::Conv2d
    bn3::BatchNorm2d
    downsample::Union{Sequential, Nothing}
    stride
end

function Bottleneck(inplanes, planes; stride=1, downsample=nothing, base_width=64, dilation=1)
    width = Int(planes * (base_width / 64.0))
    expansion = 4
    # Both self.conv2 and self.downsample layers downsample the input when stride != 1
    return Bottleneck(
        conv1x1(inplanes, width),
        BatchNorm2d(width),
        conv3x3(width, width; stride=stride, dilation=dilation),
        BatchNorm2d(width),
        conv1x1(width, planes * expansion),
        BatchNorm2d(planes * expansion),
        downsample,
        stride
    )
end

expansion(block::Type{Bottleneck}) = 4

function (m::Bottleneck)(x::AbstractArray{T, 4}) where T
    identity = x
    out = m.conv1(x)
    out = m.bn1(out)
    out = relu.(out)
    out = m.conv2(out)
    out = m.bn2(out)
    out = relu.(out)
    out = m.conv3(out)
    out = m.bn3(out)
    if m.downsample != nothing
        identity = m.downsample(x)
    end
    out = out .+ identity
    out = relu.(out)
    return out
end

################################################################################
#                                  ResNet                                      #
################################################################################

function make_layer(block, inplanes::Int, planes::Int, blocks::Int, base_width::Int; stride=1, dilation=1)
    norm_layer = BatchNorm2d
    # in the original code there's also `dilate` parameter which leads to pretty convoluted
    # stateful logic, thus we ignore it in the initial implementation
    # ResNet should still behave as usual given default parameters
    downsample = nothing
    if stride != 1 || inplanes != planes * expansion(block)
        downsample = Sequential(
            conv1x1(inplanes, planes * expansion(block); stride=stride),
            norm_layer(planes * expansion(block))
        )
    end
    layers = []
    push!(layers, block(inplanes, planes; stride=stride, downsample=downsample,
                        base_width=base_width, dilation=dilation))
    # expand input to the next set of blocks
    inplanes = planes * expansion(block)
    for _=2:blocks
        push!(layers, block(inplanes, planes; base_width=base_width, dilation=dilation))
    end
    return Sequential(layers...), inplanes
end

mutable struct ResNet
    conv1::Conv2d
    bn1::BatchNorm2d
    layer1::Sequential
    layer2::Sequential
    layer3::Sequential
    layer4::Sequential
    fc::Linear
end

function ResNet(block::Union{Type{BasicBlock}, Type{Bottleneck}}, block_nums::Vector{Int};
                num_classes=1000, zero_init_residual=false, width_per_group=64)
    norm_layer = BatchNorm2d
    inplanes = 64
    dilation = 1
    base_width = width_per_group
    conv1 = Conv2d(3, inplanes, 7; stride=2, padding=3, bias=false)
    bn1 = norm_layer(inplanes)
    layer1, inplanes = make_layer(block, inplanes, 64, block_nums[1], base_width)
    layer2, inplanes = make_layer(block, inplanes, 128, block_nums[2], base_width; stride=2)
    layer3, inplanes = make_layer(block, inplanes, 256, block_nums[3], base_width; stride=2)
    layer4, inplanes = make_layer(block, inplanes, 512, block_nums[3], base_width; stride=2)
    fc = Linear(512 * expansion(block), num_classes)
    resnet = ResNet(conv1, bn1, layer1, layer2, layer3, layer4, fc)
    walkstruct(resnet) do s
        if s isa Conv2d
            init_kaiming_normal!(s.W; mode=:fan_out, nonlinearity=:relu)
        elseif s isa BatchNorm2d
            init_constant!(s.gamma, 1.0)
            init_constant!(s.beta, 0.0)
        end
    end
    if zero_init_residual
        walkstruct(resnet) do s
            if s isa Bottleneck
                init_constant!(s.bn3.gamma, 0.0)
            elseif s isa BasicBlock
                init_constant!(s.bn2.gamma, 0.0)
            end
        end
    end
    return resnet
end

function (m::ResNet)(x::AbstractArray{T, 4}) where T
    x = m.conv1(x)
    x = m.bn1(x)
    x = relu.(x)
    x = maxpool2d(x, 3, stride=2, padding=1)

    x = m.layer1(x)
    x = m.layer2(x)
    x = m.layer3(x)
    x = m.layer4(x)

    # using maxpool2d instead of AdaptiveAvgPool since later isn't implemented yet
    # AdaptiveXXXPool computes kernel size automatically based on output size
    # (which is (1,1) in PyTorch implementation)
    # instead we compute kernel size ourselves as size(x)[1:2] .+ (1, 1) .- (1, 1) ./ (1, 1)
    # which is just (4, 4)
    x = maxpool2d(x, (4, 4))
    x = dropdims(x, dims=(1,2))
    x = m.fc(x)

    return x
end


text:
import J4Org: evaluate
import J4Org: with_hash_evaluate

@testset "extract_links" begin
    r=evaluate("4+5",Main)

end

# "#!" support
@testset "extract_links_full" begin
    s="# some text\n#\n#!x=5+6\n#!ones(5)\n#!y=5;\n# some text \n"
    r=with_hash_evaluate(s,"BoxingModule")
    @test r=="# some text\n#\n# #+BEGIN_SRC julia\n# x=5+6\n# ones(5)\n# y=5;\n# #+END_SRC\n# #+BEGIN_SRC julia\n# 11\n# 5-element Array{Float64,1}:\n#  1.0\n#  1.0\n#  1.0\n#  1.0\n#  1.0\n# #+END_SRC\n# some text \n"
end

# "# !" support
@testset "extract_links_full" begin
    s="# some text\n#\n# !x=5+6\n# !ones(5)\n# !y=5;\n# some text \n"
    r=with_hash_evaluate(s,"BoxingModule")
    @test r=="# some text\n#\n# #+BEGIN_SRC julia\n# x=5+6\n# ones(5)\n# y=5;\n# #+END_SRC\n# #+BEGIN_SRC julia\n# 11\n# 5-element Array{Float64,1}:\n#  1.0\n#  1.0\n#  1.0\n#  1.0\n#  1.0\n# #+END_SRC\n# some text \n"
end


text:
mutable struct FactorTreeGibbsOperatorXMLElement <: OperatorXMLElement
    el::XMLOrNothing
    lfm::LatentFactorModelXMLElement
    tdl::TraitLikelihoodXMLElement
    weight::Real

    function FactorTreeGibbsOperatorXMLElement(lfm::LatentFactorModelXMLElement,
                    tdl::TraitLikelihoodXMLElement)
        return new(nothing, lfm, tdl, 1.0)
    end
end

function make_xml(op::FactorTreeGibbsOperatorXMLElement)
    el = new_element(bn.FACTOR_TREE_GIBBS_OP)
    set_attribute(el, bn.WEIGHT, op.weight)
    make_xml(op.lfm)
    make_xml(op.tdl)
    add_ref_el(el, op.lfm.el)
    add_ref_el(el, op.tdl.el)
    op.el = el
    return el
end


text:
using Documenter, BlockArrays

# Build documentation.
# ====================

makedocs(
    # options
    modules = [BlockArrays],
    clean   = true,
    doctest = false,
    )

# Deploy built documentation from Travis.
# =======================================

deploydocs(
    # options
    deps = Deps.pip("mkdocs", "python-markdown-math"),
    repo = "github.com/KristofferC/BlockArrays.jl.git"
)


text:
function pWaves(ngrid,H,Np,Br,m,T)

# Evaluates the frequency and spatial structure of
# polar magnetic Rossby waves in a spherical geometry.
# The wave equation is formulated in terms of the
# meridional component of the magnetic field.
# The meridional coordinate is y = cos(colat), although
# a second grid array is stored as latitude for plotting.
# The radial magnetic field is assumed to be constant
# over the surface of the CMB. The same governing equation
# describes equatorial MAC waves in a spherical geometry. The
# two types of waves are distinguished by the sign of the frequency.
# The equatorial waves travel eastward (positive freq.) and the
# high latitude waves travel westward (negative freq.)
#
# input
#
# ngrid - grid used to represent the waves
# H    - layer thickness (km)
# Np   - dimensionless buoyancy frequency (N/Omega)
# Br   - radial magnetic field  (mT)
# m    - angular order m
# T    - initial estimate of period  (years, -'ve' for westward)
#
# output
#
# lat  -  array of latitude for plotting eigenfunction by(y)
# y    -  array of coordinates "y"
# byp   -  eigenfunction by(y)  (modified according to Eq. 13)
# omega - frequency of wave
# pd    - period of wave (years)
# q     - quality factor

# initial estimate of frequency
sectoyear = 365.25 * 24 * 60 * 60;
omega0 = 2 * pi / (T*sectoyear) + 0.0im;
B = Br * 1.0e-3;    # convert to Tesla

# set up grid
y,lat = mgrid(ngrid);

# solve eigenvalue problem for byp
omega,byp = modes_y(y,H,Np,B,m,omega0);

#  compute b_x
bxp = modes_x(y,byp,H,Np,B,m,omega);

# transform back to original variables
ynew,bx,by = transform(y,bxp,byp);

# evaluate period and quality factor of wave
pd = 2 * pi / (sectoyear * real(omega));
q = abs(0.5*real(omega)/imag(omega));

return pd,q,bx,by,ynew

end

max_stars_repo_name: bbuffett/MagRossby

max_issues_repo_name: bbuffett/MagRossby

max_forks_repo_name: bbuffett/MagRossby


text:
using BinaryProvider
using CxxWrap

prefix() =  BinaryProvider.Prefix(dirname(dirname(CxxWrap.jlcxx_path)))

products = Product[]
for basename in ["jlcxx_containers", "except", "extended", "functions", "hello", "basic_types", "inheritance", "parametric", "pointer_modification", "types"]
  fullname = "lib"*basename
  push!(products, LibraryProduct(prefix(), fullname, Symbol(fullname)))
end

if any(!satisfied(p; verbose=true) for p in products)
  error("libcxxwrap-julia was not built with the example libraries required for testing")
end

write_deps_file(joinpath(@__DIR__, "deps.jl"), products)


text:
Base.@pure determine_chunksize(u,alg::DiffEqBase.DEAlgorithm) = determine_chunksize(u,get_chunksize(alg))
Base.@pure function determine_chunksize(u,CS)
  if CS != 0
    return CS
  else
    return ForwardDiff.pickchunksize(length(u))
  end
end

struct NLSOLVEJL_SETUP{CS,AD} end
  Base.@pure NLSOLVEJL_SETUP(;chunk_size=0,autodiff=true) = NLSOLVEJL_SETUP{chunk_size,autodiff}()
  (::NLSOLVEJL_SETUP)(f,u0; kwargs...) = (res=NLsolve.nlsolve(f,u0; kwargs...); res.zero)
  function (p::NLSOLVEJL_SETUP{CS,AD})(::Type{Val{:init}},f,u0_prototype) where {CS,AD}
    AD ? autodiff = :forward : autodiff = :central
    OnceDifferentiable(f, u0_prototype, u0_prototype, autodiff,
                       ForwardDiff.Chunk(determine_chunksize(u0_prototype,CS)))
  end

# wrapper for non-autonomous functions
mutable struct NonAutonomousFunction{F,autonomous}
  f::F
  t
  p
end
(p::NonAutonomousFunction{F,true})(res,u) where F = p.f(res, u)
(p::NonAutonomousFunction{F,false})(res,u) where F = p.f(res, u, p.p, p.t)

mutable struct ManifoldProjection{autonomous,F,NL,O}
  g::F
  nl_rhs
  nlsolve::NL
  nlopts::O

  function ManifoldProjection{autonomous}(g, nlsolve, nlopts) where {autonomous}
    # replace residual function if it is time-dependent
    # since NLsolve only accepts functions with two arguments
    _g = NonAutonomousFunction{typeof(g),autonomous}(g, 0,0)
    new{autonomous,typeof(_g),typeof(nlsolve),typeof(nlopts)}(_g, _g, nlsolve, nlopts)
  end
end

# Now make `affect!` for this:
function (p::ManifoldProjection{autonomous,NL})(integrator) where {autonomous,NL}
  # update current time if residual function is time-dependent
  if !autonomous
    p.g.t = integrator.t
  end
  p.g.p = integrator.p

  integrator.u .= p.nlsolve(p.nl_rhs, integrator.u; p.nlopts...)
end

function Manifold_initialize(cb,u::Number,t,integrator)
  cb.affect!.nl_rhs = cb.affect!.nlsolve(Val{:init}, cb.affect!.g, [u])
  u_modified!(integrator,false)
end

function Manifold_initialize(cb,u,t,integrator)
  cb.affect!.nl_rhs = cb.affect!.nlsolve(Val{:init}, cb.affect!.g, u)
  u_modified!(integrator,false)
end

function ManifoldProjection(g; nlsolve=NLSOLVEJL_SETUP(), save=true,
                            autonomous=DiffEqBase.numargs(g)==2, nlopts=Dict{Symbol,Any}())
  affect! = ManifoldProjection{autonomous}(g, nlsolve, nlopts)
  condtion = (u,t,integrator) -> true
  save_positions = (false,save)
  DiscreteCallback(condtion, affect!;
                   initialize = Manifold_initialize,
                   save_positions=save_positions)
end

export ManifoldProjection


text:
function least_squares(A::AbstractArray{T}, b::AbstractArray{U}) where {T,U}
    if (T<:ForwardDiff.Dual && U<:ForwardDiff.Dual)
        m = size(A,1)
        n = size(A,2)

        Av = map(ForwardDiff.value, A)
        bv = map(ForwardDiff.value, b)

        P = pinv(Av)
        xv = P * bv

        Ap = map(ForwardDiff.partials, A)
        bp = map(ForwardDiff.partials, b)
        Apt = Ap[reshape(1:m*n,(m,n))'] # taking the transpose
        bpt = reshape(bp,m,1)

        Pp = P*P'*Apt*(I - Av*P) + (I - P*Av)*Apt*P'*P - P*Ap*P
        xp = Pp*bv + P*bpt

        x = map(T,xv,xp)
    elseif (T<:ForwardDiff.Dual)
        m = size(A,1)
        n = size(A,2)

        Av = map(ForwardDiff.value, A)

        P = pinv(Av)
        xv = P * b

        Ap = map(ForwardDiff.partials, A)
        Apt = Ap[reshape(1:m*n,(m,n))'] # taking the transpose

        Pp = P*P'*Apt*(I - Av*P) + (I - P*Av)*Apt*P'*P - P*Ap*P
        xp = Pp*bv

        x = map(T,xv,xp)
    elseif (U<:ForwardDiff.Dual)
        m = size(A,1)
        n = size(A,2)

        bv = map(ForwardDiff.value, b)

        P = pinv(A)
        xv = P * bv

        bp = map(ForwardDiff.partials, b)
        bpt = reshape(bp,m,1)

        xp = P*bpt

        x = map(T,xv,xp)
    else
        x = pinv(A) * b
    end

    # TODO other cases

    x
end


text:
# ------------------------------------------------------------------------------------------
# ## Model complexity
#
# In the last notebook, we saw that we could customize a model by adding a parameter. Doing
# so, we were able to fit that model to a data point. This fit was perfect, insofar as
# numerics would allow.
#
# In the next notebook, we'll see that as we add more data to our data set, fitting a model
# to our data usually becomes more challenging and the result will be less perfect.
#
# For one thing, we will find that we can add complexity to our model to capture added
# complexity in the data. We can do this by adding more parameters to our model. We'll see
# that for a data set with two data points, we can again get a "perfect" fit to our model by
# adding a second parameter to our model.
#
# However, we can't simply add a parameter to our model every time we add a data point to
# our data set, since this will lead to a phenomenon called **overfitting**.
#
# In the image below, we depict a data set that is close to linear, and models that exhibit
# underfitting, fitting well, and overfitting, from left to right:
#
# <img src="data/model_fitting.png" alt="Drawing" style="width: 800px;"/>
#
#
# In the first image, the model accounts for the slope along which the data falls, but not
# the offset.
#
# In the second image, the model accounts for both the slope and offset of the  data. Adding
# this second parameter (the offset) to the model creates a much better fit.
#
# However, we can imagine that a model can have too many parameters, where we begin to fit
# not only the high level features of the data, but also the noise. This overfitting is
# depicted in the third image.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Our aim will be to fit the data well, but avoiding *over*fitting the data!
# ------------------------------------------------------------------------------------------

max_stars_repo_path: .nbexports/introductory-tutorials/broader-topics-and-ecosystem/intro-to-ml/07. ML - Model complexity.jl

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
"""
```
features = haar_features(img, top_left, bottom_right, feat)
features = haar_features(img, top_left, bottom_right, feat, coordinates)
```

Returns an array containing the Haar-like features for the given Integral Image in the region specified by the points top_left and bottom_right.

The caller of the haar_features function specifies a rectangular region in the image by its top_left and bottom_right points and a particular
feature type, e.g. two rectangles side by side (:x2). This function then computes the values of all haar rectangular features for the particular
rectangle configuration for all possible locations, and all possible rectangle widths & heights if coordinates of certain specific features are not
provided. If they are provided then it only computes the values of those features whose coordinates are given. Calculating value of haar rectangular
feature corresponds to finding the difference of sums of all points in the different rectangles comprising the feature.

Parameters:

 -  img               = The Integral Image for which the Haar-like features are to be found
 -  top_left          = The top and left most point of the region where the features are to be found
 -  bottom_right      = The bottom and right most point of the region where the features are to be found
 -  feat              = A symbol specifying the type of the Haar-like feature to be found

```
    Currently, feat can take 5 values:

    :x2 = Two rectangles along horizontal axis
        +---------------------------------------+
        |                                       |
        | +-----+-----+                         |
        | |     |     |                         |
        | | -1  | +1  |  +---------->           |
        | |     |     |                         |
        | +-----+-----+                         |
        |                                       |
        |       +                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       v                               |
        |                                       |
        |                                       |
        |                                       |
        |                                       |
        |                                       |
        |                                       |
        +---------------------------------------+

    :y2 = Two rectangles along vertical axis
        +---------------------------------------+
        |                                       |
        | +------------+                        |
        | |     -1     |                        |
        | +------------+  +---------->          |
        | |     +1     |                        |
        | +------------+                        |
        |                                       |
        |       +                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       v                               |
        |                                       |
        |                                       |
        |                                       |
        |                                       |
        |                                       |
        |                                       |
        +---------------------------------------+

    :x3 = Three rectangles along horizontal axis
        +---------------------------------------+
        |                                       |
        | +-----+-----+-----+                   |
        | |     |     |     |                   |
        | | -1  | +1  | -1  |  +---------->     |
        | |     |     |     |                   |
        | +-----+-----+-----+                   |
        |                                       |
        |       +                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       v                               |
        |                                       |
        |                                       |
        |                                       |
        |                                       |
        |                                       |
        |                                       |
        +---------------------------------------+

    :y3 = Three rectangles along vertical axis
        +---------------------------------------+
        |                                       |
        | +------------+                        |
        | |     -1     |                        |
        | +------------+                        |
        | |     +1     |  +---------->          |
        | +------------+                        |
        | |     -1     |                        |
        | +------------+                        |
        |                                       |
        |       +                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       v                               |
        |                                       |
        |                                       |
        |                                       |
        |                                       |
        +---------------------------------------+

    :xy4 = Four rectangles along horizontal and vertical axes
        +---------------------------------------+
        |                                       |
        | +-----+-----+                         |
        | |     |     |                         |
        | | -1  | +1  |                         |
        | |     |     |                         |
        | +-----+-----+  +---------->           |
        | |     |     |                         |
        | | +1  | -1  |                         |
        | |     |     |                         |
        | +-----+-----+                         |
        |                                       |
        |       +                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       v                               |
        |                                       |
        |                                       |
        +---------------------------------------+

    The +1 and -1 signs show which rectangles are subtracted and which are added to evaluate the final haar feature.

```

 -  coordinates       = The user can provide the coordinates of the rectangles if only certain Haar-like features are to found in the given region. The required format is a 3 Dimensional array as (f,r,4) where f = number_of_features, r = number_of_rectangles and 4 is for the coordinates of the top_left and bottom_right point of the rectangle (top_left_y, top_left_x, bottom_right_y, bottom_right_x). The default value is nothing, where all the features are found

## References

M. Oren, C. Papageorgiou, P. Sinha, E. Osuna and T. Poggio, "Pedestrian detection using wavelet templates," Proceedings of IEEE Computer Society Conference on Computer Vision and Pattern Recognition, San Juan, 1997, pp. 193-199.

P. Viola and M. Jones, "Rapid object detection using a boosted cascade of simple features," Proceedings of the 2001 IEEE Computer Society Conference on Computer Vision and Pattern Recognition. CVPR 2001, 2001, pp. I-511-I-518 vol.1.

"""
function haar_features(img::AbstractArray{T, 2}, top_left::Array{I, 1}, bottom_right::Array{I, 1}, feat::Symbol, coordinates = nothing) where {T <: Union{Real, Color}, I <: Int}
    check_coordinates(top_left, bottom_right)
    rectangular_feature_coordinates = check_feature_type(feat)

    if coordinates == nothing
        rectangular_feature_coordinates = haar_coordinates(bottom_right[1] - top_left[1] + 1, bottom_right[2] - top_left[2] + 1, feat)
    else
        for i = 1:size(coordinates)[1]
            push!(rectangular_feature_coordinates, SMatrix{size(coordinates)[2], 4}(coordinates[i, :, :]))
        end
    end
    if length(rectangular_feature_coordinates) > 0
        num_feats = length(rectangular_feature_coordinates)
        num_rects = size(rectangular_feature_coordinates[1])[1]
    else
        println("No features of given type found in region!")
        num_feats = 0
        num_rects = 0
    end

    rectangular_features = zeros(T, num_feats, num_rects)
    for i = 1:num_feats
        a = rectangular_feature_coordinates[i]
        for j = 1:num_rects
            ytop = top_left[1] + a[j, 1] - 1
            ybot = top_left[1] + a[j, 3] - 1
            xtop = top_left[2] + a[j, 2] - 1
            xbot = top_left[2] + a[j, 4] - 1
            rectangular_features[i,j] = boxdiff(img, ytop, xtop, ybot, xbot)
        end
    end
    rectangular_feature_values = (sum(rectangular_features[:, 2:2:end], dims = 2) - sum(rectangular_features[:, 1:2:end], dims = 2))[:]
end

function haar_features(img::AbstractArray{T, 2}, top_left::Tuple{I, 2}, bottom_right::Tuple{I, 2}, feat::Symbol, coordinates = nothing) where {T <: Union{Real, Color}, I <: Int}
	haar_features(img, [top_left[1], top_left[2]], [bottom_right[1], bottom_right[2]], feat, coordinates)
end

"""
```
coordinates = haar_coordinates(height, width, feat)
```

Returns an array containing the coordinates of the Haar-like features of the specified type.

The caller of the haar_coordinates function specifies a rectangular region in the image by its height and width and a particular
feature type, e.g. two rectangles side by side (:x2). This function then computes the coordinates for the particular rectangle configuration
for all possible locations, and all possible rectangle widths & heights.

Parameters:

 -  height        = Height of the neighbourhood/window where the coordinates are computed
 -  width         = Width of the neighbourhood/window where the coordinates are computed
 -  feat          = A symbol specifying the type of the Haar-like feature to be found

```
    Currently, feat can take 5 values:

    :x2 = Two rectangles along horizontal axis
        +---------------------------------------+
        |                                       |
        | +-----+-----+                         |
        | |     |     |                         |
        | | -1  | +1  |  +---------->           |
        | |     |     |                         |
        | +-----+-----+                         |
        |                                       |
        |       +                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       v                               |
        |                                       |
        |                                       |
        |                                       |
        |                                       |
        |                                       |
        |                                       |
        +---------------------------------------+

    :y2 = Two rectangles along vertical axis
        +---------------------------------------+
        |                                       |
        | +------------+                        |
        | |     -1     |                        |
        | +------------+  +---------->          |
        | |     +1     |                        |
        | +------------+                        |
        |                                       |
        |       +                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       v                               |
        |                                       |
        |                                       |
        |                                       |
        |                                       |
        |                                       |
        |                                       |
        +---------------------------------------+

    :x3 = Three rectangles along horizontal axis
        +---------------------------------------+
        |                                       |
        | +-----+-----+-----+                   |
        | |     |     |     |                   |
        | | -1  | +1  | -1  |  +---------->     |
        | |     |     |     |                   |
        | +-----+-----+-----+                   |
        |                                       |
        |       +                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       v                               |
        |                                       |
        |                                       |
        |                                       |
        |                                       |
        |                                       |
        |                                       |
        +---------------------------------------+

    :y3 = Three rectangles along vertical axis
        +---------------------------------------+
        |                                       |
        | +------------+                        |
        | |     -1     |                        |
        | +------------+                        |
        | |     +1     |  +---------->          |
        | +------------+                        |
        | |     -1     |                        |
        | +------------+                        |
        |                                       |
        |       +                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       v                               |
        |                                       |
        |                                       |
        |                                       |
        |                                       |
        +---------------------------------------+

    :xy4 = Four rectangles along horizontal and vertical axes
        +---------------------------------------+
        |                                       |
        | +-----+-----+                         |
        | |     |     |                         |
        | | -1  | +1  |                         |
        | |     |     |                         |
        | +-----+-----+  +---------->           |
        | |     |     |                         |
        | | +1  | -1  |                         |
        | |     |     |                         |
        | +-----+-----+                         |
        |                                       |
        |       +                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       |                               |
        |       v                               |
        |                                       |
        |                                       |
        +---------------------------------------+

    The +1 and -1 signs show which rectangles are subtracted and which are added to evaluate the final haar feature.

```

"""
function haar_coordinates(h::Int, w::Int, feat::Symbol)
    rectangular_feature_coordinates = check_feature_type(feat)

    for I = 1:h
        for J = 1:w
            for i = 2:h
                for j = 2:w
                    if feat == :x2 && I + i - 2 <= h && J + 2*j - 3 <= w
                        push!(rectangular_feature_coordinates, (SMatrix{4, 2, Int}(I, J, I + i - 2, J + j - 2, I, J + j - 1, I + i - 2, J + 2*j - 3))')
                    elseif feat == :y2 && I + 2*i - 3 <= h && J + j - 2 <= w
                        push!(rectangular_feature_coordinates, (SMatrix{4, 2, Int}(I, J, I + i - 2, J + j - 2, I + i - 1, J, I + 2*i - 3, J + j - 2))')
                    elseif feat == :x3 && I + i - 2 <= h && J + 3*j - 4 <= w
                        push!(rectangular_feature_coordinates, (SMatrix{4, 3, Int}(I, J, I + i - 2, J + j - 2, I, J + j - 1, I + i - 2, J + 2*j - 3, I, J + 2*j - 2, I + i - 2, J + 3*j - 4))')
                    elseif feat == :y3 && I + 3*i - 4 <= h && J + j  - 2 <= w
                        push!(rectangular_feature_coordinates, (SMatrix{4, 3, Int}(I, J, I + i - 2, J + j - 2, I + i - 1, J, I + 2*i - 3, J + j - 2, I + 2*i - 2, J, I + 3*i - 4, J + j - 2))')
                    elseif feat == :xy4 && I + 2*i - 3 <= h && J + 2*j - 3 <= w
                        push!(rectangular_feature_coordinates, (SMatrix{4, 4, Int}(I, J, I + i - 2, J + j - 2, I, J + j - 1, I + i - 2, J + 2*j - 3, I + i - 1, J + j - 1, I + 2*i - 3, J + 2*j - 3, I + i - 1, J, I + 2*i - 3, J + j - 2))')
                    end
                end
            end
        end
    end

    return rectangular_feature_coordinates
end

function check_coordinates(top_left::Array{I, 1}, bottom_right::Array{I, 1}) where I <: Int
    if length(top_left) != 2 || length(bottom_right) != 2
        throw(ArgumentError("The top_left and bottom_right point must have only 2 coordinates."))
    end

    if bottom_right[1] < top_left[1]
        throw(ArgumentError("The bottom_right point must be lower than the top_left point."))
    end

    if bottom_right[2] < top_left[2]
        throw(ArgumentError("The bottom_right point must be towards the right of the top_left point."))
    end
end

function check_feature_type(feat::Symbol)
    if feat == :x2 || feat == :y2
        rectangular_feature_coordinates = Vector{SMatrix{2, 4, Int}}()
    elseif feat == :x3 || feat == :y3
        rectangular_feature_coordinates = Vector{SMatrix{3, 4, Int}}()
    elseif feat == :xy4
        rectangular_feature_coordinates = Vector{SMatrix{4, 4, Int}}()
    else
        throw(ArgumentError("The type of the feature must be either :x2, :y2, :x3, :y3 or :xy4."))
    end

    return rectangular_feature_coordinates
end


text:
"""
This constitutes the MAIN LOOP for all the simulation, it creates storage variables
as well as some arrays to compute the full evolution of the system. Most of the operations here
are carried out on the GPU device, except for the accumulation operations for the compressibility
factor and other things.
"""
function move!(
    positions::AbstractArray,
    forces::AbstractArray,
    syst::System,
    dynamics::Dynamic,
    potential::Potential,
    cycles::Integer;
    filenames::AbstractArray = nothing,
    thermal::Integer = Int(cycles / 2),
    msd::Bool = false,
    logging::Bool = false,
)
# * Accumulation variables
    total_energy = 0.0f0
    total_pressure = 0.0f0
    big_z = 0.0f0
    total_z = 0.0f0
    total_virial = 0.0f0
    samples = 0
# MSD
    msdval = 0.0f0
    ref_pos = Matrix{Float32}(undef, 3, syst.N)
    time_acc = 0.0f0
    hst_positions = Matrix{Float32}(undef, 3, syst.N)

# GPU variables for the kernels
    nthreads = 512
    gpu_blocks = ceil(Int, syst.N / nthreads)

# Allocate memory for random numbers
    rnd_matrix = Matrix{Float32}(undef, 3, syst.N)

# When equilibrating, always use periodic boundary conditions
    pbc = true

# ! Main loop
    @showprogress for i = 1:cycles
    # * Create array of random numbers
        Random.randn!(syst.rng, rnd_matrix)
        rnd_matrix .*= sqrt(2.0f0 * dynamics.τ)
        cu_rand_matrix = cu(rnd_matrix)

    # * Move the particles following the Ermak-McCammon algorithm
        CUDA.@sync begin
            @cuda threads = nthreads blocks = gpu_blocks gpu_ermak!(
            positions,
            forces,
            syst.N,
            syst.L,
            dynamics.τ,
            cu_rand_matrix,
            pbc
        )
        end
    # ! Always set forces to zero
        CUDA.fill!(forces, 0.0f0)

    # ! Compute the energy and forces, O(N^2) complexity
        CUDA.@sync begin
            @cuda threads = nthreads blocks = gpu_blocks gpu_energy!(
            positions,
            forces,
            syst.N,
            syst.L,
            syst.rc,
            potential.a,
            potential.b,
            potential.λ,
            potential.temp,
            syst.ener,
            syst.press,
        )
        end

    # * Start accumulating averages
        if i >= thermal
        # ! Turn off periodic boundary conditions
            pbc = false

        # ! Update the samples counter
            samples += 1

        # * Update the total energy of the system
            total_energy = sum(syst.ener)

        # * Extract the virial from the kernel
            total_virial = sum(syst.press)
            total_virial /= (3.0f0 * syst.N)
            big_z = 1.0f0 + total_virial
            total_z += big_z # This is just an accumulator
        # * Update the total pressure of the system
            total_pressure = big_z * syst.ρ

        # ! Save energy, pressure and compressibility to file
            open(filenames[1], "a") do io
                filener = total_energy / (syst.N * samples)
                println(io, "$(filener),$(total_pressure),$(big_z)")
            end

            if msd
            # * Reference positions for the MSD, at first time step
                if samples == 1
                    ref_pos = Array(positions)
                end

            # * Accumulate MSD
                if samples > 1
                    time_acc += dynamics.τ
                    hst_positions = Array(positions)
                    msdval = msd!(hst_positions, ref_pos)
                    msdval /= syst.N

                    open(filenames[2], "a") do io
                        println(io, "$(time_acc),$(msdval)")
                    end

                    savemsd(hst_positions, filenames[3])
                end
            end

            if logging
                # * Save snapshot of the system
                snapshot(Array(positions), filenames[4], i)
            end
        end
    end

    # ! Adjust results as averages
    total_energy /= syst.N
    total_z /= samples

    # * Show results
    println("Energy: $(total_energy)")
    println("Pressure: $(total_pressure)")
    println("Compressibility: $(big_z)")
    println("Compressibility (average): $(total_z)")
end


text:
@testset "Manopt Gradient Descent" begin
  # Test the gradient descent with
  # the distance function squared
  # on S1, such that we can easily also verify exp and log
  M = Circle()
  r = [-π/2,π/4,0.0,π/4];
  f = S1Point.(r);
  F(x) = 1/10*sum(distance.(Ref(M),f,Ref(x)).^2)
  ∇F(x) = 1/5*sum(-log.(Ref(M),Ref(x),f))
  x,rec = steepestDescent(M,F,∇F,f[1];
    stoppingCriterion = stopWhenAny(stopAfterIteration(200), stopWhenChangeLess(10^-16)),
    stepsize = ArmijoLinesearch(1.,exp,0.99,0.1),
    debug = [:Iteration," ",:Cost, :Stop, 100,"\n"],
    record = [:Iteration, :Cost, 1]
  )
  # after one step for local enough data -> equal to real valued data
  @test abs(getValue(x)-sum(r)/length(r)) ≈ 0 atol=5*10.0^(-14)
  # Test Fallbacks -> we can't do steps with the wrong combination
  p = SubGradientProblem(M,F,∇F)
  o = GradientDescentOptions(f[1],stopAfterIteration(20),ConstantStepsize(1.))
  @test_throws ErrorException initializeSolver!(p,o)
  @test_throws ErrorException doSolverStep!(p,o,1)
  @test_throws ErrorException getSolverResult(p,o)
end


text:
struct State_fixed{T<:AbstractVector}
    st::T
    ope::Op_fixed
    State_fixed(st,ope) = length(st) != binomial(dim(ope)+nume(ope)-1,nume(ope)) ? error("lenght of vector does not match dimension") : new{typeof(st)}(st/sqrt(st'*st),ope)
end

st(s::State_fixed) = s.st

ope(s::State_fixed) = s.ope

nume(s::State_fixed) = nume(s.ope)

dim(s::State_fixed) = dim(s.ope)

basis(s::State_fixed) = basis(s.ope)

typ(s::State_fixed) = eltype(s.st)

function rhosp(sta::State_fixed, precis=15)
    n = dim(sta)
    num = nume(sta)
    rhospv = spzeros(typ(sta),n,n)
    estate = st(sta)
    o = ope(sta)
    for i in 1:n
        for j in 1:n
            rhospv[i,j] = round(estate'*bdb(o, i, j)*estate, digits = precis)
        end
    end
    return rhospv
end


text:
__precompile__()

"""
This package contains routines for computing eigenvalues of second
order ordinary differential equations and in particular the
quasinormal modes (QNMs) of black holes in General Relativity using
the "Asymptotic Iteration Method" [1] using the implementation
based on the "improved" version of the AIM, described in [2].

References:

[1] (https://arxiv.org/abs/math-ph/0309066v1)
[2] (https://arxiv.org/abs/1111.5024)
"""
module QuasinormalModes

# ------------------------------------------------------------------
# 1. Imports
# ------------------------------------------------------------------

using SymEngine
using Polynomials
using PolynomialRoots
using TaylorSeries
using NLsolve
using Roots

# ------------------------------------------------------------------
# 2. Public API (Exports)
# ------------------------------------------------------------------

# --- Type hierarchy ---
export AIMProblem
export AnalyticAIMProblem
export NumericAIMProblem
export QuadraticEigenvalueProblem

# --- Extra (concrete) types ---
export AIMCache

# --- Steping methods (singleton types) ---
export Serial
export Threaded

# --- Mandatory methods ---
export λ0, S0, get_niter, get_x0
export get_ODEvar, get_ODEeigen

# --- AIM methods ---
export computeDelta!
export computeEigenvalues
export eigenvaluesInGrid

# ------------------------------------------------------------------
# 3. Type hierarchy
# ------------------------------------------------------------------

"""
Parent super-type of all problems that can be solved using the AIM.
"""
abstract type AIMProblem{N<:Unsigned,T<:Number} end

"""
Parent super-type of all problems that can be solved using the AIM semi-analytically.
"""
abstract type AnalyticAIMProblem{N<:Unsigned,T<:Number} <: AIMProblem{N,T} end

"""
Parent super-type of all problems that can be solved using the AIM numerically.
"""
abstract type NumericAIMProblem{N<:Unsigned,T<:Number} <: AIMProblem{N,T} end

"""
Parent super-type of all problems whose eigenvalue is a quadratic polynomial.
"""
abstract type QuadraticEigenvalueProblem{N<:Unsigned,T<:Number} <: AnalyticAIMProblem{N,T} end

# ------------------------------------------------------------------
# 4. Traits
# ------------------------------------------------------------------

"""
Super-type of traits describing the analyticity of eigenvalue problems.
"""
abstract type AnalyticityTrait end

"""
All problems with eigenvalues that *can* be described by analytic functions have this trait.
"""
struct IsAnalytic <: AnalyticityTrait end

"""
All problems with eigenvalues that *can't* be described by analytic functions have this trait.
"""
struct IsNumeric <: AnalyticityTrait end

"""
The default trait of AIMProblem(s).
"""
AnalyticityTrait(::Type{<:AIMProblem}) = IsNumeric()

"""
The trait of AnalyticAIMProblem(s).
"""
AnalyticityTrait(::Type{<:AnalyticAIMProblem}) = IsAnalytic()

"""
The trait of NumericAIMProblem(s).
"""
AnalyticityTrait(::Type{<:NumericAIMProblem}) = IsNumeric()

"""
All problem types must implement a λ0 function.
This behavior is enforced by the default implementations.
"""
λ0(x::T) where {T} = λ0(AnalyticityTrait(T), x)
λ0(::IsAnalytic, x) = error("Please implement a λ0 function for ", typeof(x))
λ0(::IsNumeric, x) = error("Please implement a λ0 function for ", typeof(x))

"""
All problem types must implement a S0 function.
This behavior is enforced by the default implementations.
"""
S0(x::T) where {T} = S0(AnalyticityTrait(T), x)
S0(::IsAnalytic, x) = error("Please implement a S0 function for ", typeof(x))
S0(::IsNumeric, x) = error("Please implement a S0 function for ", typeof(x))

"""
All problem types must implement get_niter to return the number of iterations to perform.
"""
get_niter(x::T) where {T} = get_niter(AnalyticityTrait(T), x)
get_niter(::IsAnalytic, x) = error("Please implement a get_niter function for ", typeof(x))
get_niter(::IsNumeric, x) = error("Please implement a get_niter function for ", typeof(x))

"""
All problem types must implement get_x0 to return AIM's point of evaluation.
"""
get_x0(x::T) where {T} = get_x0(AnalyticityTrait(T), x)
get_x0(::IsAnalytic, x) = error("Please implement a get_x0 function for ", typeof(x))
get_x0(::IsNumeric, x) = error("Please implement a get_x0 function for ", typeof(x))

"""
Analytic problems must implement an acessor to the variable of the ODE.
"""
get_ODEvar(x::T) where {T} = get_ODEvar(AnalyticityTrait(T), x)
get_ODEvar(::IsAnalytic, x) = error("Please implement a get_ODEvar function for ", typeof(x))
get_ODEvar(::IsNumeric, x) = error("Numeric problem", typeof(x), "cannot implement get_ODEvar")

"""
Analytic problems must implement an acessor to the eigenvalue of the ODE.
"""
get_ODEeigen(x::T) where {T} = get_ODEeigen(AnalyticityTrait(T), x)
get_ODEeigen(::IsAnalytic, x) = error("Please implement a get_ODEeigen function for ", typeof(x))
get_ODEeigen(::IsNumeric, x) = error("Numeric problem", typeof(x), "cannot implement get_ODEeigen")

# ------------------------------------------------------------------
# 5. AIM stepping methods
# ------------------------------------------------------------------

"""
Super-type of all stepping methods used internally by the AIM.
"""
abstract type AIMSteppingMethod end

"""
Perform all AIM steps sequentially
"""
struct Serial <: AIMSteppingMethod end

"""
Perform all AIM steps in parallel using threads
"""
struct Threaded <: AIMSteppingMethod end

# ------------------------------------------------------------------
# 6. Includes
# ------------------------------------------------------------------

# --- Include core functionality ---
include("./core/aim_cache.jl")
include("./core/analytic_tools.jl")
include("./core/aim_step.jl")
include("./core/compute_delta.jl")

# --- Include eigenvalue computing functionality ---
include("./eigenvalue_computing/compute_eigenvalues.jl")
include("./eigenvalue_computing/eigenvalues_in_grid.jl")

end # QuasinormalModes


text:
module Optimization

using Test

using Darwin
# import Darwin: crossover!, mutate!, selection

using LearningStrategies: learn!, MaxIter, strategy, Verbose
using Distributions

include("test_functions.jl")

const Entity = Vector{Float64}

function run_with(fitness, bounds, selection, crossover, mutation, generations)
    initial_population = Population([rand(bounds, 2) for _ in 1:128])
    model = PopulationModel(initial_population, fitness)
    strat = strategy(Verbose(GAStrategy(selection, crossover, mutation)),
                     MaxIter(generations))
    learn!(model, strat)
end

function test_with(fn, selection, crossover, mutation, generations)
    fitness = FitnessFunction{Entity}(function (x::Entity)
                                          -pack(fn, Entity)(x)
                                      end)
    Bounds = Uniform(bounds(fn)...)

    model = run_with(fitness, Bounds, selection, crossover, mutation, generations)

    @test fitness(model.fittest.genome) ≈ minimum(fn) atol = 0.01
    @test model.fittest.genome ≈ argmin(fn) atol = 0.15
end

function test_skip_with(fn, selection, crossover, mutation, generations)
    fitness = FitnessFunction{Entity}(function (x::Entity)
                                          -pack(fn, Entity)(x)
                                      end)
    Bounds = Uniform(bounds(fn)...)

    model = run_with(fitness, Bounds, selection, crossover, mutation, generations)

    @test_skip fitness(model.fittest.genome) ≈ minimum(fn) atol = 0.01
    @test_skip model.fittest.genome ≈ argmin(fn) atol = 0.15
end

@testset "Rosenbrock" begin
    rosenbrock = Rosenbrock{Float64}()

    test_with(rosenbrock,
              SoftmaxSelection{Entity, 2, 2}(ExponentialRate(1.0, 0.5, 0.999)),
              ArithmeticCrossover{Entity, 2, 2}(ConstantRate(0.7)),
              PointwiseMutation{Entity}(ConstantRate(0.2), Uniform(bounds(rosenbrock)...)),
              10_000)

    test_with(rosenbrock,
              TournamentSelection{Entity, 5, 1, 2}(),
              UniformCrossover{Entity, 2, 1}(LinearRate(0.8, 0.2, 1/20)),
              AdditiveMutation{Entity, Uniform}(ConstantRate(0.3), 2.0, -2.0, 2.0),
              10_000)

    test_skip_with(rosenbrock,
                   L1Selection{Entity, 2, 2}(),
                   ArithmeticCrossover{Entity, 2, 2}(ConstantRate(0.7)),
                   AdditiveMutation{Entity, Normal}(ConstantRate(0.3), 2.0, -2.0, 2.0),
                   10_000)
end

end # module


text:
module TextSegmentation

include("Utls.jl")
include("TextTiling.jl")
include("C99.jl")
include("TopicTiling.jl")

export Utls, TextTiling, C99, TopicTiling

end


text:
struct CutIterData #used for plotting performance of method
    time::Float64
    Cut::Int128
    obj::Float64
    bestbound::Float64
    boundgap::Float64
    status::Symbol
end

mutable struct Cut
    p::Float64
    ∇s::Array{Float64}
    status::Symbol
  end

  struct Dual
    α::Array{Float64}
    λ::Float64
    βl::Array{Float64}
    βu::Array{Float64}
    ρ::Array{Float64}
    w::Array{Float64}
    ofv::Float64
    status::Symbol
  end

  struct SparsePortfolioData
    μ::Array{Float64}
    Y::Array{Float64}
    d::Array{Float64}
    X::Array{Float64, 2}
    A::Array{Float64, 2}
    l::Array{Float64}
    u::Array{Float64}
    k::Int64
    n::Int64 # No stocks
    m::Int64 # No constraints
    f::Int64 # Rank of cov matrix
    min_investment::Array{Float64, 1}
    γ::Array{Float64}
  end

struct CardinalityConstrainedPortfolio
  indices::Array{Int}
  w::Array{Float64}
  λ::Float64
  α::Array{Float64}
  βl::Array{Float64}
  βu::Array{Float64}
  ρ::Array{Float64}
  boundGap::Float64
  boundGapSOCP::Float64
  numIters::Int128
  expectedReturn::Float64
  portfolioVariance::Float64
  solveTime::Float64
  nodeCount::Float64
end

struct PortfolioData #Old struct, not currently in use.
  μ::Array{Float64}
  X::Array{Float64, 2}
  normalization::Array{Float64, 2}
  X_normalized::Array{Float64, 2}
  min_investment::Array{Float64, 1}
end


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

export DataSet

mutable struct DataSet
    keys::Vector{FileKey}
end

const ds_ignore_line_expr = r"^(\s*#.*)?$"

function Base.read(::Type{DataSet}, filename::AbstractString)
    keys = open(filename) do input
        [FileKey(strip(l)) for l in eachline(input) if !occursin(ds_ignore_line_expr, l)]
    end
    DataSet(keys)
end

import Base.==
==(a::DataSet, b::DataSet) = a.keys == b.keys


text:
using Test

ϕ=60
a1 = [1,0].*2
a2 = [cos(π*ϕ/180),sin(π*ϕ/180)].*2
l = Generic2DLattice(a1,a2)
conf = CorrugatedMorse(5,0.4,0.1,l)

x = collect(0:0.1:10)
y = collect(0:0.1:10)
v = Array{Float64,2}(undef, length(x),length(y))

for (i1,xi) in enumerate(x), (i2,yi) in enumerate(y)
    z = 0
    v[i1,i2] = potential(conf,xi,yi,z)
end

using Plots
gr()
#plot(x,y,v,seriestype=:heatmap,aspect_ratio=:equal)
#gui()

start = 0
V0 = []
for s in 5:20:100
push!(V0, real(calc_VG([0,0],l,conf,0.2,s)))
end
plot(V0,seriestype=:line)
gui()

#start = 0
#z = 0.3
#V0 = fill(0.,7,7)
#for gx in -3:3, gy in -3:3
#    V0[gx+4,gy+4] = real(calc_VG([gx,gy],l,conf,z,10))
#end
#@show V0
#plot(V0,seriestype=:heatmap)
#gui()


text:
using HMMs
using Brobdingnag
using Ipopt
using MathProgBase
using Random

struct HMMOptimizer{Dist,Calc,Out} <: MathProgBase.AbstractNLPEvaluator
    hmm::HMMs.HMM{Dist,Calc,Out}
    current::Vector{Float64}
end

function HMMOptimizer( hmm::HMMs.HMM{Dist,Calc,Out} ) where {Dist,Calc,Out}
    m = length(hmm.initialprobabilities)
    return HMMOptimizer( hmm, fill( NaN, m*(m+2) ) )
end

function MathProgBase.initialize( ::HMMOptimizer{Dist,Calc,Out}, requested_features::Vector{Symbol} ) where {Dist,Calc,Out}
    unimplemented = setdiff( requested_features, [:Grad,:Jac] )
    if !isempty( unimplemented )
        error( "The following features aren't implemented: " * join( string.(unimplemented), ", " ) )
    end
end

MathProgBase.features_available( ::HMMOptimizer{Dist,Calc,Out} ) where {Dist,Calc,Out} = [:Grad]

function MathProgBase.eval_f( opt::HMMOptimizer{Dist,Brob,Float64}, x ) where {Dist}
    HMMs.setparameters!( opt.hmm, x )
    return HMMs.likelihood( opt.hmm )[end].log
end

function MathProgBase.eval_g( opt::HMMOptimizer{Dist,Calc,Out}, g, x ) where {Dist,Calc,Out}
    HMMs.setparameters!( opt.hmm, x )
    g[:] = sum(opt.hmm.transitionprobabilities, dims=2)
end

function MathProgBase.eval_grad_f( opt::HMMOptimizer{Dist,Brob,Float64}, grad_f, x ) where {Dist}
    HMMs.setparameters!( opt.hmm, x )
    l = HMMs.likelihood( opt.hmm )[end]
    dl = HMMs.dlikelihood( opt.hmm )[:,end]
    grad_f[:] = convert( Vector{Float64}, dl ./ l )
end

function MathProgBase.jac_structure( opt::HMMOptimizer{Dist,Calc,Out} ) where {Dist,Calc,Out}
    m = length(opt.hmm.initialprobabilities)
    rows = zeros(Int, m^2)
    cols = zeros(Int, m^2)
    inner = 1
    for i = 1:m
        for j = 1:m
            rows[inner] = i
            cols[inner] = (i-1)*m + j
            inner += 1
        end
    end
    return (rows,cols)
end

function MathProgBase.eval_jac_g( opt::HMMOptimizer{Dist,Calc,Out}, J, x ) where {Dist,Calc,Out}
    m = length(opt.hmm.initialprobabilities)
    J[:] = ones(m^2)
end

m = 3
graph = HMMs.fullyconnected( m )
hmm = HMMs.randomhmm( graph, calc=Brob, seed=1 )
y = rand( hmm, 10_000 );

mlhmm = HMMs.randomhmm( graph, calc=Brob, seed=2 )
HMMs.setobservations( mlhmm, y );

solver = IpoptSolver(
    max_iter=1_000,
    derivative_test="first-order",
    option_file_name="/home/atteson/ipopt.opt",
    print_level=5,
)
model = MathProgBase.NonlinearModel(solver)
opt = HMMOptimizer(mlhmm)

MathProgBase.jac_structure( opt )
MathProgBase.loadproblem!(model, m*(m+2), m, [zeros(m^2); fill(-Inf,m); zeros(m)], fill(Inf,m*(m+2)),
                          ones(m), ones(m), :Max, opt)
parameters = fill( NaN, m*(m+2) )
HMMs.getparameters!( opt.hmm, parameters )

grad_f = fill(NaN,m*(m+2))
MathProgBase.eval_grad_f( opt, grad_f, parameters )

errors = Float64[]
delta = 1e-2
for i = 1:length(parameters)
    parameters[i] += delta
    fp = MathProgBase.eval_f( opt, parameters )
    parameters[i] -= delta

    parameters[i] -= delta
    fm = MathProgBase.eval_f( opt, parameters )
    parameters[i] += delta

    d = (fp - fm)/(2*delta)
    push!( errors, abs(d / grad_f[i] - 1) )
    if( errors[end] > 1e-4 )
        @warn( "Mismatch at derivative $i" )
    end
end

(rows,cols) = MathProgBase.jac_structure( opt )

MathProgBase.setwarmstart!( model, parameters )
MathProgBase.optimize!(model)

HMMs.likelihood( opt.hmm )

emhmm = HMMs.randomhmm( graph, calc=Brob, seed=2 )
HMMs.setobservations( emhmm, y );
HMMs.em( emhmm, debug=2 )
HMMs.likelihood( opt.hmm )[end]
HMMs.likelihood( emhmm )[end]

HMMs.em( opt.hmm, debug=2 )
HMMs.likelihood( opt.hmm )

HMMs.getparameters!( opt.hmm, parameters )
HMMs.getparameters!( emhmm, parameters )
MathProgBase.eval_f( opt, parameters )


text:
using Plots

function cplx_heatmap(M)
    plot(
        heatmap(abs.(M)),
        heatmap(real.(M)),
        heatmap(imag.(M)),
        layout = (1,3),
        size = (1600,400)
    )
end

function interp_nd(G, P, X, Y)
iv = zeros(size(P)[1], length(X))
for i ∈ 1:length(X)
    index_x = argmin(abs.(G.-X[i]))
    index_y = argmin(abs.(G.-Y[i]))

    iv[:, i] = P[:, index_y, index_x]
end
return iv
end

function interp(G, P, X, Y)
iv = zeros(length(X))
for i ∈ 1:length(X)
    index_x = argmin(abs.(G.-X[i]))
    index_y = argmin(abs.(G.-Y[i]))

    iv[i] = P[index_y, index_x]
end
return iv
end


text:
module LotkaVolterra2dPlots

    using GeometricBase
    using GeometricEquations
    using GeometricProblems.Diagnostics

    using LaTeXStrings
    using Measures: mm
    using RecipesBase

    """
    Plots the solution of a 2D Lotka-Volterra model together with the energy error.

    Arguments:
    * `sol <: Solution`
    * `equ <: Equation`

    Keyword aguments:
    * `nplot=1`: plot every `nplot`th time step
    * `xlims=:auto`: xlims for solution plot
    * `ylims=:auto`: ylims for solution plot
    * `latex=true`: use LaTeX guides
    """
    @userplot Plot_Lotka_Volterra_2d
    @recipe function f(p::Plot_Lotka_Volterra_2d; nplot=1, nt=:auto, xlims=:auto, ylims=:auto, latex=true)
        if length(p.args) != 2 || !(typeof(p.args[1]) <: Solution) || !(typeof(p.args[2]) <: Equation)
            error("Lotka-Volterra plot should be given two arguments: a solution and an equation. Got: $(typeof(p.args))")
        end
        sol = p.args[1]
        equ = p.args[2]
        params = equ.parameters

        if nt == :auto
            nt = ntime(sol)
        end

        if nt > ntime(sol)
            nt = ntime(sol)
        end

        H, ΔH = compute_invariant_error(sol.t, sol.q, (t,q) -> equ.invariants[:h](t,q,params))

        size   := (1000,400)
        layout := (1,2)
        legend := :none

        right_margin := 10mm
        bottom_margin := 10mm

        guidefontsize := 18
        tickfontsize  := 12

        @series begin
            subplot := 1
            left_margin := 10mm

            if ntime(sol) ≤ 200
                markersize := 5
            else
                markersize  := 1
                markercolor := 1
                linecolor   := 1
                markerstrokewidth := 1
                markerstrokecolor := 1
            end

            # seriestype := :scatter
            if latex
                xguide := L"x_1"
                yguide := L"x_2"
            else
                xguide := "x₁"
                yguide := "x₂"
            end
            xlims  := xlims
            ylims  := ylims
            aspect_ratio := 1
            sol.q[1,0:nplot:nt], sol.q[2,0:nplot:nt]
        end

        @series begin
            subplot := 2
            if latex
                xguide := L"t"
                yguide := L"[H(t) - H(0)] / H(0)"
            else
                xguide := "t"
                yguide := "[H(t) - H(0)] / H(0)"
            end
            xlims  := (sol.t[0], Inf)
            ylims  := :auto
            yformatter := :scientific
            right_margin := 10mm
            sol.t[0:nplot:nt], ΔH[0:nplot:nt]
        end
    end

    """
    Plots the solution of a 2D Lotka-Volterra model.

    Arguments:
    * `sol <: Solution`
    * `equ <: Equation`

    Keyword aguments:
    * `nplot=1`: plot every `nplot`th time step
    * `xlims=:auto`: xlims for solution plot
    * `ylims=:auto`: ylims for solution plot
    * `latex=true`: use LaTeX guides
    """
    @userplot Plot_Lotka_Volterra_2d_Solution
    @recipe function f(p::Plot_Lotka_Volterra_2d_Solution; nplot=1, nt=:auto, xlims=:auto, ylims=:auto, latex=true)
        if length(p.args) != 2 || !(typeof(p.args[1]) <: Solution) || !(typeof(p.args[2]) <: Equation)
            error("Lotka-Volterra solution plot should be given two arguments: a solution and an equation. Got: $(typeof(p.args))")
        end
        sol = p.args[1]
        equ = p.args[2]
        params = equ.parameters

        if nt == :auto
            nt = ntime(sol)
        end

        if nt > ntime(sol)
            nt = ntime(sol)
        end

        if ntime(sol) ≤ 200
            markersize := 5
        else
            markersize  := 1
            markercolor := 1
            linecolor   := 1
            markerstrokewidth := 1
            markerstrokecolor := 1
        end

        legend := :none
        size := (400,400)

        # solution
        @series begin
            # seriestype := :scatter
            if latex
                xguide := L"x_1"
                yguide := L"x_2"
            else
                xguide := "x₁"
                yguide := "x₂"
            end
            xlims  := xlims
            ylims  := ylims
            aspect_ratio := 1
            guidefontsize := 18
            tickfontsize  := 12
            sol.q[1,0:nplot:nt], sol.q[2,0:nplot:nt]
        end
    end

    """
    Plots time traces of the solution of a 2D Lotka-Volterra model and its energy error.

    Arguments:
    * `sol <: Solution`
    * `equ <: Equation`

    Keyword aguments:
    * `nplot=1`: plot every `nplot`th time step
    * `latex=true`: use LaTeX guides
    """
    @userplot Plot_Lotka_Volterra_2d_Traces
    @recipe function f(p::Plot_Lotka_Volterra_2d_Traces; nplot=1, nt=:auto, latex=true)
        if length(p.args) != 2 || !(typeof(p.args[1]) <: Solution) || !(typeof(p.args[2]) <: Equation)
            error("Lotka-Volterra traces plot should be given two arguments: a solution and an equation. Got: $(typeof(p.args))")
        end
        sol = p.args[1]
        equ = p.args[2]
        params = equ.parameters

        if nt == :auto
            nt = ntime(sol)
        end

        if nt > ntime(sol)
            nt = ntime(sol)
        end

        H, ΔH = compute_invariant_error(sol.t, sol.q, (t,q) -> equ.invariants[:h](t,q,params))

        size   := (800,600)
        legend := :none
        guidefontsize := 18
        tickfontsize  := 12

        # traces

        layout := (3,1)

        if latex
            ylabels = (L"x_1", L"x_2")
        else
            ylabels = ("x₁", "x₂")
        end

        for i in 1:2
            @series begin
                subplot := i
                yguide := ylabels[i]
                xlims  := (sol.t[0], Inf)
                xaxis := false
                right_margin := 10mm
                sol.t[0:nplot:nt], sol.q[i,0:nplot:nt]
            end
        end

        @series begin
            subplot := 3
            if latex
                xguide := L"t"
                yguide := L"[H(t) - H(0)] / H(0)"
            else
                xguide := "t"
                yguide := "[H(t) - H(0)] / H(0)"
            end
            xlims  := (sol.t[0], Inf)
            yformatter := :scientific
            guidefontsize := 18
            tickfontsize  := 12
            right_margin := 10mm
            sol.t[0:nplot:nt], ΔH[0:nplot:nt]
        end
    end

end


text:
@testset "Unit tests" begin
    @testset "requests.jl" begin
        @testset "_add_trailing_slash" begin
            @test FHIRClientXML._add_trailing_slash(HTTP.URI("http://julialang.org")) == HTTP.URI("http://julialang.org/")
            @test FHIRClientXML._add_trailing_slash(HTTP.URI("http://julialang.org/")) == HTTP.URI("http://julialang.org/")
        end
    end
end


text:
using PtFEM

ProjDir = dirname(@__FILE__)

data = Dict(
  # Beam(ndim, nst, nxe, nip, direction, finite_element(nod, nodof), axisymmet
  # ?Beam for more details
  :struc_el => Beam(2, 1, 9, 1, :x, Line(2, 2), false),
  :properties => [4.0e4; 2.0e4],
  :etype => [1,1,2,2,2,2,2,2,2],
  :x_coords => [0.0, 2.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 8.0, 10.0],
  :support => [
    (1, [0 1]),
    (9, [0 1])
    ],
  :fixed_freedoms => [
    (1, 2, -0.001),
    (3, 1, -0.005)],
  :loaded_nodes => [
      (2, [-20.0 0.0])
      (3, [-1.0 -0.02083])
      (4, [-1.0 0.0])
      (5, [-1.0 0.0])
      (6, [-1.0 0.0])
      (7, [-1.0 0.0])
      (8, [-1.0 0.0])
      (9, [-3.8 -0.77917])
      (10, [-1.2 0.5333])
    ],
  :penalty => 1e19,
  :eq_nodal_forces_and_moments => [
    (3, [-1.0 -0.02083 -1.0 0.02083]),
    (4, [-1.0 -0.02083 -1.0 0.02083]),
    (5, [-1.0 -0.02083 -1.0 0.02083]),
    (6, [-1.0 -0.02083 -1.0 0.02083]),
    (7, [-1.0 -0.02083 -1.0 0.02083]),
    (8, [-1.0 -0.02083 -1.0 0.02083]),
    (9, [-2.8 -0.8 -1.2  0.5333])
  ]
)

data[:loaded_nodes]

data |> display
println()

@time fem, dis_df, fm_df = p43(data)
println()

display(dis_df)
println()
display(fm_df)

if VERSION.minor < 6
  using Plots
  gr(size=(400,600))

  p = Vector{Plots.Plot{Plots.GRBackend}}(undef, 3)
  titles = ["p43.1 y deflection", "p43.1 y shear force", "p43.1 z moment"]
  fors = vcat(
    convert(Array, fm_df[:, :xl_Force]),
    convert(Array, fm_df[:, :xr_Force])[end]
  )
  moms = vcat(
    convert(Array, fm_df[:, :xl_Moment]),
    convert(Array, fm_df[:, :xr_Moment])[end]
  )

  p[1] = plot(fem.displacements[2,:], ylim=(-0.005, 0.005),
    xlabel="node", ylabel="deflection [m]", color=:red,
    line=(:dash,1), marker=(:circle,4,0.8,stroke(1,:black)),
    title=titles[1], leg=false)
  p[2] = plot(fors, ylim=(-20.0, 25),
    xlabel="node", ylabel="y shear force [N]", color=:blue,
    line=(:dash,1), marker=(:circle,4,0.8,stroke(1,:black)),
    title=titles[2], leg=false)
  p[3] = plot(moms, ylim=(-40, 40),
    xlabel="node", ylabel="z moment [Nm]", color=:green,
    line=(:dash,1), marker=(:circle,4,0.8,stroke(1,:black)),
    title=titles[3], leg=false)

  plot(p..., layout=(3, 1))
  savefig(ProjDir*"/Ex43.1b.png")

end


text:
using LibPQ
using Intervals: AbstractInterval, Bounded, Closed, Open, Unbounded

# String representation of interval bounds
_left_bound(::Type{Open}) = '('
_left_bound(::Type{Closed}) = '['
_right_bound(::Type{Open}) = ')'
_right_bound(::Type{Closed}) = ']'

# Define string literals for intervals and vectors of intervals
_libpq_literal(x) = x    # default: identity

_libpq_literal(r::AbstractRange) = string('[', r.start, ',', r.stop, ']')
function _libpq_literal(i::AbstractInterval{T,L,R}) where {T,L<:Bounded,R<:Bounded}
    return string(_left_bound(L), first(i), ',', last(i), _right_bound(R))
end
function _libpq_literal(i::AbstractInterval{T,L,R}) where {T,L<:Bounded,R<:Unbounded}
    return string(_left_bound(L), first(i), ",)")
end
function _libpq_literal(i::AbstractInterval{T,L,R}) where {T,L<:Unbounded,R<:Bounded}
    return string("(,", last(i), _right_bound(R))
end
_libpq_literal(i::AbstractInterval{T,L,R}) where {T,L<:Unbounded,R<:Unbounded} = "(,)"
_libpq_literal(xs::Union{AbstractVector,Tuple}) = [_libpq_literal(x) for x in xs]

_build_params(q::SQLQuery, params) = Tuple(_libpq_literal(params[a]) for a in q.args)

# The following methods overload the SQLQuery call operator for LibPQ.jl

"""
    $(SIGNATURES)

Execute query `q` with `params` using connection `conn`

"""
function (q::SQLQuery{<:Union{Select,InsertUpdateDelete}})(conn::LibPQ.Connection, params=NamedTuple())
    return execute(conn, q.query, _build_params(q, params))
end

function (q::SQLQuery{<:Script})(conn)
    execute(conn, "BEGIN;")
    execute(conn, q.query)
    execute(conn, "COMMIT;")
    return nothing
end

function (q::SQLQuery{<:Union{InsertReturning,SelectOne}})(conn::LibPQ.Connection, params=NamedTuple())
    result = execute(conn, q.query, _build_params(q, params))
    if num_rows(result) == 0
        close(result)
        return missing
    elseif num_rows(result) > 1
        @warn "$(q.name) returned more than one row"
    end
    nt = NamedTuple{tuple(Symbol.(result.column_names)...)}([v for v in first(result)])
    close(result)
    return nt
end

function (q::SQLQuery{<:InsertUpdateDeleteMany})(conn::LibPQ.Connection, params=NamedTuple())
    execute(conn, "BEGIN;")
    LibPQ.load!(NamedTuple(k => _libpq_literal(v) for (k,v) in pairs(params) if k in q.args),
                conn,
                q.query)
    execute(conn, "COMMIT;")
    return nothing
end

function (q::SQLQuery{<:Union{SelectValue,InsertReturningValue}})(conn::LibPQ.Connection, params=NamedTuple())
    result = execute(conn, q.query, _build_params(q, params))
    if num_rows(result) == 0
        close(result)
        return missing
    elseif num_rows(result) > 1
        @warn "$(q.name) returned more than one row"
    end
    row = first(result)
    if length(row) > 1
        @warn "$(q.name) returned more than one value"
    end
    value = first(row)
    close(result)
    return value
end


text:
module Client
import Rock
using Logging
@Logging.configure(level=DEBUG)
type t
    host::ASCIIString
    port::Int64
end

function make(host::ASCIIString, port::Int64)
    t(host, port)
end
function _connect(c::t)
    connect(c.host, c.port)
end
function command(c::t, cmd::Rock.RPC.Command)
    @debug "Issuing Rock Command ($cmd) to $c"
    conn = _connect(c)
    try
        serialize(conn, cmd)
        flush(conn)
        deserialize(conn)
    catch err
        rethrow(err)
    finally
        close(conn)
    end
end

end

max_stars_repo_name: JeremyRubin/Rock


text:
SRC_DIR = joinpath(Base.source_dir(), "..", "..", "src")
OUT_DIR = joinpath(Base.source_dir(), "..", "..", "output")
mkpath(OUT_DIR)

#include(joinpath(SRC_DIR, "BridgeSDEInference.jl"))
#using Main.BridgeSDEInference
include(joinpath(SRC_DIR, "BridgeSDEInference_for_tests.jl"))

using Distributions # to define priors
using Random        # to seed the random number generator
using DataFrames
using CSV

DIR = "auxiliary"
include(joinpath(SRC_DIR, DIR, "read_and_write_data.jl"))
include(joinpath(SRC_DIR, DIR, "transforms.jl"))
include(joinpath(SRC_DIR, DIR, "utility_functions.jl"))

# decide if first passage time observations or partially observed diffusion
fptObsFlag = true

# pick dataset
filename = "path_fpt_simpleConjug.csv"

# fetch the data
(df, x0, obs, obs_time, fpt,
      fptOrPartObs) = readData(Val(fptObsFlag), joinpath(OUT_DIR, filename))
x0
param = :complexConjug
θ_init = [10.0, -8.0, 25.0, 0.0, 3.0]
P˟ = FitzhughDiffusion(param, θ_init...)
P̃ = std_aux_laws(FitzhughDiffusionAux, param, θ_init, obs, obs_time, 1)
L = @SMatrix [1. 0.]
Σ = @SMatrix [10^(-10)]

setup = MCMCSetup(P˟, P̃, fptOrPartObs)
set_observations!(setup, [L for _ in P̃], [Σ for _ in P̃], obs, obs_time, fpt)
set_imputation_grid!(setup, 1/5000)
set_transition_kernels!(setup,
                        [RandomWalk([3.0, 5.0, 1.0, 0.01, 0.5], 5)],
                        0.99, true, 3,
                        (MetropolisHastingsUpdt(),))
set_priors!(setup, Priors((ImproperPrior(),)), KnownStartingPt(x0))
set_mcmc_params!(setup, 3*10^4, 3*10^2, 10^2, 10^1, 100)
set_blocking!(setup, ChequeredBlocking(),
              (collect(1:length(obs)-2)[1:1:end], 10^(-10), SimpleChangePt(100)))
set_solver!(setup, Vern7(), NoChangePt())
initialise!(eltype(x0), setup)

Random.seed!(4)
out, elapsed = @timeit mcmc(setup)
display(out.accpt_tracker)

include(joinpath(SRC_DIR, DIR, "plotting_fns.jl"))
plot_chains(out; truth=[10.0, -8.0, 25.0, 0.0, 3.0])
plot_paths(out; transf=[(x,θ)->x, (x,θ)->conjugToRegular(x, θ[1], 0)],
           obs=(times=obs_time[2:end], vals=obs[2:end], indices=1))


text:
module MovingMesh

using FiniteDifferenceDerivatives
using ArrayViews
using DASSL
using Iterators

export Equation, solve

abstract Equation

function epsilon_eval(eps::Function,r,u,ur)
    return eps(r,u,ur)
end

function epsilon_eval(eps::Real,r,u,ur)
    return eps
end

function extracty(y,dy,npts,npde)
    t   =  y[1]
    dt  = dy[1]
    # r   = view( y,2:npts+1)
    # dr  = view(dy,2:npts+1)
    # u   = reshape_view( view( y,npts+2:length(y)), (npts, npde))
    # du  = reshape_view( view(dy,npts+2:length(y)), (npts, npde))
    r   =  y[2:npts+1]
    dr  = dy[2:npts+1]
    u   = reshape(  y[npts+2:length(y)], npts, npde)
    du  = reshape( dy[npts+2:length(y)], npts, npde)
    return t, dt, r, dr, u, du
end

function computederivatives(r,dr,u,du,derivative)
    dudr = zero(u)          # du/dr
    dudt = zero(u)
    tmp = Array(eltype(u),size(u,1))
    for j = 1:size(u,2)
        derivative(tmp,view(u,:,j),r)
        dudr[:,j] = tmp
        for i = 1:size(u,1)
            dudt[i,j] = du[i,j]-dr[i]*dudr[i,j]
        end
    end
    return dudr, dudt
end

# simple in-place gaussian filter to smoothen the monitor function
function smoothen!(M; passes = 4, args...)
    if passes == 0
        return
    else
        # forward pass
        for i = 1:length(M)-1
            M[i]=(M[i+1]+M[i])/2
        end
        M[end] = M[end-1]
        # backward pass
        for i = length(M):-1:2
            M[i]=(M[i]+M[i-1])/2
        end
        smoothen!(M, passes=passes-1)
    end
end

function movingmeshres(r,dr,M,gamma,gval,epsilon)
    npts = length(dr)
    dxi  = 1/(npts-1)
    resr = zero(dr)          # residua of moving mesh
    for i = 2:npts-1
        resr[i]   = (dr[i]-gamma*(dr[i-1]+dr[i+1]-2*dr[i])/dxi^2)-gval/epsilon*( (M[i+1]+M[i])*(r[i+1]-r[i])/dxi-(M[i]+M[i-1])*(r[i]-r[i-1])/dxi )/dxi
    end
    resr[ 1 ] =  r[ 1 ]   # mesh boundary conditions
    resr[end] = dr[end]
    return resr
end

function physicalres!(resu,rhsval,r,u,dudt,gval)
    npts, npde = size(u)
    for j = 1:npde
        resu[1,j] = u[1,j]          # enforce boundary conditions u(0,t)=0
        for i = 1:npts
            resu[i,j] = gval*rhsval[i,j]-dudt[i,j]
        end
    end
end

function meshinit(eqn,
                  uinit,
                  nx;
                  T = Float64,
                  epsilon = 1e-5,
                  gamma = 2,
                  derivative = (du,u,r)->fdd!(du,u,r,1,3),
                  args...)

    info("Mesh initialization: Start")
    r0   = linspace(eqn.rspan...,nx)
    dxi  = 1/(convert(T,nx)-1)

    dudr = zero(uinit(r0))      # du/dr

    function F(t,r,dr)
        u    = uinit(r)::Matrix{T}
        for j = 1:size(u,2)
            derivative(view(dudr,:,j),view(u,:,j),r)
        end
        M = eqn.monitor(r,u,dudr)::Array{T}
        smoothen!(M; args...)
        ϵ   = epsilon_eval(epsilon,r,u,dudr)
        res  = zero(dr)
        res[  1] = dr[  1]
        res[end] = dr[end]
        for i = 2:nx-1
            res[i]=ϵ*(dr[i]*dxi^2-gamma*(dr[i-1]+dr[i+1]-2*dr[i]))-( (M[i+1]+M[i])*(r[i+1]-r[i])-(M[i]+M[i-1])*(r[i]-r[i-1]) )
        end
        return res
    end

    u0    = uinit(r0)
    for j = 1:size(u0,2)
        derivative(view(dudr,:,j),view(u0,:,j),r0)
    end

    ϵ   = epsilon_eval(epsilon,r0,u0,dudr)
    tol = min(1,ϵ*1e7)
    err = Inf
    for (t,r,dr) in dasslIterator(F, r0, zero(T); reltol = dxi*tol, abstol = 1e-2*dxi*tol)
        err = norm(dr)*dxi
        if err < eps(T)^(4/5)
            info("Mesh initialization: Done! Total error = $err")
            return r
        end
    end

    warn("Mesh initialization: Failed! Total error = $err")
    return r0
end

function solve(eqn   :: Equation,
               uinit :: Function;
               npts :: Int = 100,
               T     = Float64,
               epsilon = 1e-5,
               gamma = 2,
               derivative = (du,u,r)->fdd!(du,u,r,1,3),
               args...)

    # generate initial data
    r0 = meshinit(eqn, uinit, npts; epsilon=epsilon, derivative=derivative, gamma=gamma, T=T, args...)
    u0 = uinit(r0)
    resu = zero(u0)
    y0 = vcat(eqn.t0,r0,vec(u0))

    # dae to solve
    function dae(tau,y,dy)
        t, dt, r, dr, u, du = extracty(y,dy,npts,eqn.npde)::(T,T,Vector{T},Vector{T},Matrix{T},Matrix{T})
        dudr, dudt = computederivatives(r,dr,u,du,derivative)

        M    = eqn.monitor(r,u,dudr)::Array{T}
        smoothen!(M; args...)
        gval = eqn.sundman(r,u,dudr)::T
        eps  = epsilon_eval(epsilon,r,u,dudr)::T
        rhs  = eqn.rhs(r,u)

        res  = Array(T,npts*eqn.npde+npts+1)
        res[1]          = dt-gval # Sundman transform equations
        res[2:npts+1]   = movingmeshres(r,dr,M,gamma,gval,eps) # moving mesh equations
        physicalres!(resu,rhs,r,u,dudt,gval) # physical equations
        res[npts+2:end] = vec(resu)

        return res
    end

    function extract_results(t,y,yt)
        t = y[1];
        x = y[2:npts+1]
        u = reshape(y[npts+2:end],npts,eqn.npde)
        dt = yt[1];
        dx = yt[2:npts+1]/dt
        du = reshape(yt[npts+2:end],npts,eqn.npde)/dt
        return t, x, u, dt, dx, du
    end

    sol = dasslIterator(dae, y0, zero(T); args...)

    return imap(z->extract_results(z...),sol)

end

end # module


text:
import Base: getindex, push!, empty!, length, sort!

export
    Population,
    getindex,
    push!,
    empty!,
    length,
    randomize!,
    sort!,
    release!

immutable Population{C, T<:AbstractComparer}
    ctx::C
    comparer::T
    decode::Function
    entities::Vector{Entity}
end

Population(ctx, comparer, decode) =
    Population(ctx, comparer, decode, Vector{Entity}())

length(pop::Population) = length(pop.entities)

getindex(pop::Population, idx) = pop.entities[idx]

push!(pop::Population, entity::Entity) = (push!(pop.entities, entity); entity)

push!(pop::Population, dna::DNA) = push!(pop, Entity(dna, pop.decode(dna)))

function randomize!(pop::Population, populationSize::Int, dnaSize::Int)
    empty!(pop)
    for i in 1:populationSize
        dna = DNA(pop.ctx, dnaSize)
        randomizeUniform!(dna)
        push!(pop, dna)
    end
    sort!(pop)
end

function sort!(pop::Population)
    sort!(pop.entities, by = e -> e.body, lt = (a, b) -> fn(pop.comparer)(a, b))
    first(pop.entities)
end

function release!(pop::Population)
    for entity in pop.entities
        release!(entity)
    end
    empty!(pop)
end

function empty!(pop::Population)
    empty!(pop.entities)
end

max_stars_repo_name: unbornchikken/julia-ml-proto

max_issues_repo_name: unbornchikken/julia-ml-proto

max_forks_repo_name: unbornchikken/julia-ml-proto