text:
module Replay

using Crayons

include("FakePTYs.jl")
using .FakePTYs: open_fake_pty

const CTRL_C = '\x03'
const UP_ARROW = "\e[A"
const DOWN_ARROW = "\e[B"
const RIGHT_ARROW = "\e[C"
const LEFT_ARROW = "\e[D"

export CTRL_C, UP_ARROW, DOWN_ARROW, RIGHT_ARROW, LEFT_ARROW
export replay

function clearline(; move_up::Bool = false)
    buf = IOBuffer()
    print(buf, "\x1b[2K") # clear line
    print(buf, "\x1b[999D") # rollback the cursor
    move_up && print(buf, "\x1b[1A") # move up
    print(buf |> take! |> String)
end

function clearlines(H::Integer)
    for i = 1:H
        clearline(move_up = true)
    end
end

function type_with_ghost(line::AbstractString)
    juliaprompt = "julia> "
    clearline()
    for index in collect(eachindex(line))
        print(crayon"green bold", juliaprompt)
        print(crayon"reset")
        println(join(line[begin:index]))
        clearlines(1)
        duration = if 30 < length(line)
            0.0125
        elseif 15 < length(line) < 30
            0.05
        else
            0.1
        end
        sleep(duration)
    end
end

function setup_pty(color = :yes; julia_project = "@."::AbstractString)
    if color in [:yes, true]
        color = "yes"
    else
        color = "no"
    end
    pts, ptm = open_fake_pty()
    blackhole = Sys.isunix() ? "/dev/null" : "nul"
    julia_exepath = joinpath(Sys.BINDIR::String, Base.julia_exename())
    replproc = withenv(
        "JULIA_HISTORY" => blackhole,
        "JULIA_PROJECT" => "$julia_project",
        "CI" => get(ENV, "CI", "false"),
        "TERM" => ""
    ) do
        # MyNote
        # -e 'if ENV["CI"] != ... ' is required to pass Pkg.test(), or "IOError: stream is closed or unusable" will happen.
        # idk why we need print("\x1b[?25l") inside of -e '...'
        # It seems print("\x1b[?25l") inside of the replay function does not work???
        run(
            ```$(julia_exepath)
                --cpu-target=native --startup-file=no --color=$(color)
                -e 'if ENV["CI"] != "true"; using Pkg; Pkg.instantiate(); print("\x1b[?25l"); end'
                -i ```,
            pts, pts, pts; wait = false
        )
    end
    Base.close_stdio(pts)
    return replproc, ptm
end

function replay(repl_lines::Vector{T}, buf::IO = stdout; color = :yes, use_ghostwriter = false, julia_project = "@.") where {T<:AbstractString}
    # c.f. MyNote above
    print("\x1b[?25l") # hide cursor
    replproc, ptm = setup_pty(color; julia_project)
    # Prepare a background process to copy output from process until `pts` is closed
    output_copy = Base.BufferStream()
    tee = @async try
        while !eof(ptm)
            l = readavailable(ptm)
            write(buf, l)
            Sys.iswindows() && (sleep(0.1); yield(); yield()) # workaround hang - probably a libuv issue?
            write(output_copy, l)
        end
        close(output_copy)
        close(ptm)
    catch ex
        close(output_copy)
        close(ptm)
        if !(ex isa Base.IOError && ex.code == Base.UV_EIO)
            rethrow() # ignore EIO on ptm after pts dies
        end
    end
    # wait for the definitive prompt before start writing to the TTY
    readuntil(output_copy, "julia>")
    sleep(0.1)
    readavailable(output_copy)

    # let's replay!
    for line in repl_lines
        sleep(1)
        bytesavailable(output_copy) > 0 && readavailable(output_copy)

        use_ghostwriter && type_with_ghost(line)

        if endswith(line, "\x03")
            write(ptm, line)
        else
            write(ptm, line, "\n")
        end
        readuntil(output_copy, "\n")
        # wait for the next prompt-like to appear
        # NOTE: this is rather inaccurate because the Pkg REPL mode is a special flower
        readuntil(output_copy, "\n")
        readuntil(output_copy, "> ")
    end

    sleep(1)
    use_ghostwriter && type_with_ghost("exit()")
    write(ptm, "exit()\n")
    sleep(1)
    wait(tee)
    success(replproc) || Base.pipeline_error(replproc)
    close(ptm)
    print("\x1b[?25h") # unhide
    return buf
end

replay(repl_script::String, buf::IO = stdout; color = :yes, use_ghostwriter = false, julia_project = "@.") = replay(split(repl_script::String, '\n'; keepempty = false), buf; color, use_ghostwriter, julia_project)

end # module


text:
module TestC

    import ..TestB.c as y
    @show y
end

ext: jl


text:
#src execute = true

# # [Parsing Data](@id parsing)
#
# `PowerSystems.jl` supports the creation of a `System` from a variety of common data formats:
#
# - [MATPOWER](@ref pm_data) (code copied with permission from [`PowerModels.jl`](https://github.com/lanl-ansi/PowerModels.jl))
# - [PSS/e RAW Files](@ref pm_data) (code copied with permission from [`PowerModels.jl`](https://github.com/lanl-ansi/PowerModels.jl))
# - [PSS/e DYR Files](@ref dyr_data)
# - [PowerSystems table data (CSV Files)](@ref table_data)
#
# ## [MATPOWER / PSS/e](@id pm_data)
#
# The following code will create a System from a MATPOWER or PSS/e file:
#
# ```julia
# sys = System(joinpath(data_dir, "case5.m"))
# ```
#
# ## [PSS/e dynamic data parsing](@id dyr_data)
#
# PSS/e's dynamic model library is extensive, we currently support parsing a limited amount
# of models out of the box.
#
# | Machine models | AVR Models | Prime Movers | PSS models |
# |----------------|------------|--------------|------------|
# | GENSAE         | IEEET1     | HYGOV        | IEEEST     |
# | GENSAL         | ESDC1A     | IEEEG1       |            |
# | GENROE         | ESAC1A     | GGOV1        |            |
# | GENCLS         | ESST4B     |              |            |
# | GENROU         | EXAC2      |              |            |
# |                | EXPIC1     |              |            |
# |                | ESAC6A     |              |            |
# |                | EXAC1      |              |            |
# |                | SCRX       |              |            |
# |                | ESDC2A     |              |            |
#
# ### Creating a Dynamic System using .RAW and .DYR data
#
# A `PowerSystems.jl` system can be created using a .RAW and a .DYR file. In this example we will create the following three bus system using the following RAW file:
#
# ```raw
# 0, 100, 33, 0, 0, 60  / 24-Apr-2020 19:28:39 - MATPOWER 7.0.1-dev
#
#
#      101, 'BUS 1       ',       138, 3,    1,    1, 1,           1.02,        0,  1.1,  0.9,  1.1,  0.9
#      102, 'BUS 2       ',       138, 2,    1,    1, 1,           1.0142,           0,  1.1,  0.9,  1.1,  0.9
#      103, 'BUS 3       ',       138, 2,    1,    1, 1,           1.0059,           0,  1.1,  0.9,  1.1,  0.9
# 0 / END OF BUS DATA, BEGIN LOAD DATA
#      101,  1, 1,    1,    1,       100,       20, 0, 0, 0, 0, 1, 1, 0
#      102,  1, 1,    1,    1,       70,       10, 0, 0, 0, 0, 1, 1, 0
#      103,  1, 1,    1,    1,       50,       10, 0, 0, 0, 0, 1, 1, 0
# 0 / END OF LOAD DATA, BEGIN FIXED SHUNT DATA
# 0 / END OF FIXED SHUNT DATA, BEGIN GENERATOR DATA
#      101,  1,       20,         0,       100,      -100,    1.02, 0,     100, 0, 0, 0, 0, 1, 1, 100,       318,         0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1
#      102,  1,       100,         0,       100,      -100,   1.0142, 0,     100, 0, 0.7, 0, 0, 1, 1, 100,       318,         0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1
#      103,  1,       100,         0,       100,      -100,   1.0059, 0,     100, 0, 0.2, 0, 0, 1, 1, 100,       318,         0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1
# 0 / END OF GENERATOR DATA, BEGIN BRANCH DATA
#      101,      103, 1,  0.01000,     0.12,      0.0,     250,     250,     250, 0, 0, 0, 0, 1, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1
#      101,      102, 1,  0.01000,     0.12,      0.0,     250,     250,     250, 0, 0, 0, 0, 1, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1
#      102,      103, 1,  0.01000,     0.12,      0.0,     250,     250,     250, 0, 0, 0, 0, 1, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1
# 0 / END OF BRANCH DATA, BEGIN TRANSFORMER DATA
# 0 / END OF TRANSFORMER DATA, BEGIN AREA DATA
# 0 / END OF AREA DATA, BEGIN TWO-TERMINAL DC DATA
# 0 / END OF TWO-TERMINAL DC DATA, BEGIN VOLTAGE SOURCE CONVERTER DATA
# 0 / END OF VOLTAGE SOURCE CONVERTER DATA, BEGIN IMPEDANCE CORRECTION DATA
# 0 / END OF IMPEDANCE CORRECTION DATA, BEGIN MULTI-TERMINAL DC DATA
# 0 / END OF MULTI-TERMINAL DC DATA, BEGIN MULTI-SECTION LINE DATA
# 0 / END OF MULTI-SECTION LINE DATA, BEGIN ZONE DATA
# 0 / END OF ZONE DATA, BEGIN INTER-AREA TRANSFER DATA
# 0 / END OF INTER-AREA TRANSFER DATA, BEGIN OWNER DATA
# 0 / END OF OWNER DATA, BEGIN FACTS CONTROL DEVICE DATA
# 0 / END OF FACTS CONTROL DEVICE DATA, BEGIN SWITCHED SHUNT DATA
# 0 / END OF SWITCHED SHUNT DATA, BEGIN GNE DEVICE DATA
# 0 / END OF GNE DEVICE DATA, BEGIN INDUCTION MACHINE DATA
# 0 / END OF INDUCTION MACHINE DATA
# Q
# ```
#
# This system is a three bus system with three generators, three loads and three branches. The dynamic data for the generators is provided in the DYR file:
#
# ```raw
#   101 'GENROE' 1   8.000000  0.030000  0.400000  0.050000  6.500000  0.000000  1.800000
#   1.700000  0.300000  0.550000  0.250000  0.200000  0.039200  0.267200  /
#   101 'ESST1A' 1   1  1  0.01  99  -99  1  10  1  1  200  0  4  -4  4  -4  0  0  1  0  3  /
#   102 'GENCLS' 1   0.0   0.0 /
#   103 'GENCLS' 1   3.1   2.0 /
# ```
#
# That assigns a GENROU generator and a ESST1A voltage regulator at the generator located at bus 101, while classic machine models for the generators located at bus 102 and 103.
#
# To create the system we can do it passing both files directories:
#
# ```@example raw_dyr_system
# using PowerSystems
# data_dir = "../../../data" # hide
# RAW_dir = joinpath(data_dir, "PSSE_test/ThreeBusNetwork.raw")
# DYR_dir = joinpath(data_dir, "PSSE_test/TestGENCLS.dyr")
# dyn_system = System(RAW_dir, DYR_dir, runchecks = false)
# ```
#
# ## [PowerSystems Table Data](@id table_data)
#
# This is a custom format that allows users to define power system component data
# by category and column with custom names, types, and units.
#
# ### Categories
#
# Components for each category must be defined in their own CSV file. The
# following categories are currently supported:
#
# - branch.csv
# - bus.csv (required)
#   - columns specifying `area` and `zone` will create a corresponding set of `Area` and `LoadZone` objects.
#   - columns specifying `max_active_power` or `max_reactive_power` will create `PowerLoad` objects when nonzero values are encountered and will contribute to the `peak_active_power` and `peak_reactive_power` values for the
#   corresponding `LoadZone` object.
# - dc_branch.csv
# - gen.csv
# - load.csv
# - reserves.csv
# - storage.csv
#
# These must reside in the directory passed when constructing PowerSystemTableData.
#
# ### [Adding Time Series Data](@id parsing_time_series)
#
# PowerSystems requires a metadata file that maps components to their time series
# data in order to be able to automatically construct time_series from raw data
# files. The following fields are required for each time array:
#
# - `simulation`:  User description of simulation
# - `resolution`:  Resolution of time series in seconds
# - `module`:  Module that defines the abstract type of the component
# - `category`:  Type of component. Must map to abstract types defined by the "module"
#   entry (Bus, ElectricLoad, Generator, LoadZone, Reserve)
# - `component_name`:  Name of component
# - `name`:  User-defined name for the time series data.
# - `normalization_factor`:  Controls normalization of the data. Use 1.0 for
#   pre-normalized data. Use 'Max' to divide the time series by the max value in the
#   column. Use any float for a custom scaling factor.
# - `scaling_factor_multiplier_module`:  Module that defines the accessor function for the
# scaling factor
# - `scaling_factor_multiplier`:  Accessor function of the scaling factor
# - `data_file`:  Path to the time series data file
#
# Notes:
#
# - The "module", "category", and "component_name" entries must be valid arguments to retrieve
# a component using `get_component(${module}.${category}, sys, $name)`.
# - The "scaling_factor_multiplier_module" and the "scaling_factor_multiplier" entries must
# be sufficient to return the scaling factor data using
# `${scaling_factor_multiplier_module}.${scaling_factor_multiplier}(component)`.
#
# PowerSystems supports this metadata in either CSV or JSON formats. Refer to
# [RTS_GMLC](https://github.com/GridMod/RTS-GMLC/blob/master/RTS_Data/FormattedData/SIIP/timeseries_pointers.json)
# for an example.
#
# #### Performance considerations
#
# By default PowerSystems stores time series data in HDF5 files. It does not keep
# all of the data in memory. This means that every time you access a time_series
# PowerSystems will have to read the data from storage, which will add latency. If
# you know ahead of time that all of your data will fit in memory then you can
# change this behavior by passing `time_series_in_memory = true` when you create
# the System.
#
# If the time series data is stored in HDF5 then PowerSystems will use the tmp filesystem by
# default. You can change this by passing `time_series_directory = X` when you create the
# System. This is required if the time series data is larger than the amount of tmp space
# available. You can also override the location by setting the environment variable
# SIIP_TIME_SERIES_DIRECTORY to another directory.
#
# ### Customization
#
# The tabular data parser in `PowerSystems.jl` can be customized to read a variety of
# datasets by configuring:
#
#  - [which type of generator (`<:Generator`) to create based on the fuel and prime mover specifications](@ref csv_genmap)
#  - [property names](@ref csv_columns), [units](@ref csv_units), and per units conversions](@ref csv_per_unit) in *.csv files
#
# Here is an example of how to construct a System with all customizations listed in this section:
#
# ```julia
# data_dir = "/data/my-data-dir"
# base_power = 100.0
# descriptors = "./user_descriptors.yaml"
# timeseries_metadata_file = "./timeseries_pointers.json"
# generator_mapping_file = "./generator_mapping.yaml"
# data = PowerSystemTableData(
#     data_dir,
#     base_power,
#     descriptors;
#     timeseries_metadata_file = timeseries_metadata_file,
#     generator_mapping_file = generator_mapping_file,
# )
# sys = System(data, time_series_in_memory = true)
# ```
#
# Examples configuration files can be found in the [RTS-GMLC](https://github.com/GridMod/RTS-GMLC/) repo:
#
#  - [user_descriptors.yaml](https://github.com/GridMod/RTS-GMLC/blob/master/RTS_Data/FormattedData/SIIP/user_descriptors.yaml)
#  - [generator_mapping.yaml](https://github.com/GridMod/RTS-GMLC/blob/master/RTS_Data/FormattedData/SIIP/generator_mapping.yaml)
#
# #### [CSV Data Configurations](@id csv_data)
#
# ##### [Custom construction of generators](@id csv_genmap)
#
# PowerSystems supports custom construction of subtypes of the abstract type Generator based
# on `fuel` and `type`. The parsing code detects these fields in the raw data and then
# constructs the concrete type listed in the passed generator mapping file. The default file
# is `src/parsers/generator_mapping.yaml`. You can override this behavior by specifying your
# own file when constructing PowerSystemTableData.
#
# ##### [Column names](@id csv_columns)
#
# PowerSystems provides am input mapping capability that allows you to keep your own
# column names.
#
# For example, when parsing raw data for a generator the code expects a column
# called `name`. If the raw data instead defines that column as `GEN UID` then
# you can change the `custom_name` field under the `generator` category to
# `GEN UID` in your YAML file.
#
# To enable the parsing of a custom set of csv files, you can generate a configuration
# file (such as `user_descriptors.yaml`) from the defaults, which are stored
# in `src/descriptors/power_system_inputs.json`.
#
# ```python
# python ./bin/generate_config_file.py ./user_descriptors.yaml
# ```
#
# Next, edit this file with your customizations.
#
# Note that the user-specific customizations are stored in YAML rather than JSON
# to allow for easier editing. The next few sections describe changes you can
# make to this YAML file.  Do not edit the default JSON file.
#
# A list of `PowerSystemTableData` accepted columns for each CSV is available
# [here](@ref tabledata_input_config).
#
# ##### [Per-unit conversion](@id csv_per_unit)
#
# For more info on the per-unit conventions in `PowerSystems.jl`, refer to the [per-unit
# section of the system documentation](@ref per_unit).
#
# PowerSystems defines whether it expects a column value to be per-unit system base,
# per-unit device base, or natural units in `power_system_inputs.json`
# ([see default values](@ref tabledata_input_config)). If it expects a
# per-unit convention that differs from your values then you can set the `unit_system` in
# `user_descriptors.yaml` and PowerSystems will automatically convert the values. For
# example, if you have a `max_active_power` value stored in natural units (MW), but
# `power_system_inputs.json` specifies `unit_system: device_base`, you can enter
# `unit_system: natural_units` in `user_descriptors.yaml` and PowerSystems will divide
# the value by the value of the corresponding entry in the column identified by the
# `base_reference` field in `power_system_inputs.json`. You can also override the
# `base_reference` setting by adding `base_reference: My Column` to make device base
# per-unit conversion by dividing the value by the entry in `My Column`. System base
# per-unit conversions always divide the value by the system `base_power` value
# instantiated when constructing a `System`.
#
# ##### [Unit conversion](@id csv_units)
#
# PowerSystems provides a limited set of unit conversions. For example, if
# `power_system_inputs.json` indicates that a value's unit is degrees but
# your values are in radians then you can set `unit: radian` in
# your YAML file. Other valid `unit` entries include `GW`, `GWh`, `MW`, `MWh`, `kW`,
# and `kWh`.
#

# ## [`PowerSystemTableData` Accepted CSV Columns](@id tabledata_input_config)
#
# The following tables describe default CSV column definitions accepted by the
# `PowerSystemeTableData` parser defined by `src/descriptors/power_system_inputs.json`:
#

using PowerSystems # hide
function create_md() # hide
    descriptor = PowerSystems._read_config_file(joinpath( # hide
        dirname(pathof(PowerSystems)), # hide
        "descriptors", # hide
        "power_system_inputs.json", # hide
    )) # hide
    columns = [ # hide
        "name", # hide
        "description", # hide
        "unit", # hide
        "unit_system", # hide
        "base_reference", # hide
        "default_value", # hide
        "value_options", # hide
        "value_range", # hide
    ] # hide
    header = "| " * join(columns, " | ") * " |\n" * repeat("|----", length(columns)) * "|\n" # hide
    s = "" # hide
    for (cat, items) in descriptor # hide
        csv = "" # hide
        for name in PowerSystems.INPUT_CATEGORY_NAMES # hide
            if name[2] == cat # hide
                csv = name[1] # hide
                break # hide
            end # hide
        end # hide
        csv == "" && continue # hide
        s = string(s, "### $csv.csv:\n\n") # hide
        s = string(s, header) # hide
        for item in items # hide
            extra_cols = setdiff(keys(item), columns) # hide
            if !isempty(extra_cols) # hide
                throw(@error "config file fields not included in header" extra_cols) # hide
            end # hide
            row = [] # hide
            for col in columns # hide
                val = string(get(item, col, " ")) # hide
                if col == "default_value" && val == " " # hide
                    val = "*REQUIRED*" # hide
                end # hide
                push!(row, val) # hide
            end # hide
            s = string(s, "|" * join(row, "|") * "|\n") # hide
        end # hide
        s = string(s, "\n") # hide
    end # hide
    s = replace(s, r"[_$]" => s"\\\g<0>"); # hide
    return s # hide
end # hide
txt = create_md(); # hide
fname = joinpath(dirname(dirname(pathof(PowerSystems))), "docs", "src", "modeler_guide", "generated_inputs_tables.md") # hide
open(fname, "w") do io # hide
      write(io, txt) # hide
end # hide
println(" "); # hide

# APPEND_MARKDOWN("docs/src/modeler_guide/generated_inputs_tables.md")


text:
export Tokenizer

using ..Units
using TextAnalysis

struct Tokenizer <: AbstractTokenizer end;

function (tokenizer::Tokenizer)(sentence::StringDocument)
  tokens(sentence)
end

max_stars_repo_name: LarsWl/DependencyParser

max_issues_repo_name: LarsWl/DependencyParser

max_forks_repo_name: LarsWl/DependencyParser


text:
using Distributed
using Statistics
using LinearAlgebra
using DelimitedFiles
using Dates
import Base.ImmutableDict

using Random

@everywhere include("../src/VMC.jl")

function nsgd(L, P, N, N_sgd; hX=3.0)

    outfile = "../data/2d-tfim-" * string(Dates.today())
    mkpath(outfile)

    θ = 0.1 * ones(3P) .+ P*rand(3P)
    θ[3:3:end] .= zero(P)
    h = Base.ImmutableDict(:ZZ => 1.0, :X => hX, :Z => 0.0)

    lr₀ = 0.1
    lr = lr₀
    num_chain = nworkers()

    println("Initial parameters are :", θ)
    println("Threads are ", num_chain)

    M = VMC.Model(L, θ, h)
    Obs = [VMC.Thermo(P) for i in 1:num_chain]
    Grad = [VMC.Gradient(Obs[i], P) for i in 1:num_chain]
    Sto = [VMC.StochasticMatrix(Obs[i], P) for i in 1:num_chain]
    C = [VMC.Config(size(M)) for i in 1:num_chain]

    Es = zeros(20)
    norms = zeros(20)
    grad_list = zeros(3P, 20)

    @time E, grad, mag, σ², sgrad = VMC.pmc(C, M, Obs, Grad, 10^4, Sto, num_chain)

    T = 0

    for i in 1:N_sgd
        M = VMC.Model(L, θ, h)
        Obs = [VMC.Thermo(P) for i in 1:num_chain]
        @time E, grad, mag, σ², sgrad = VMC.pmc(C, M, Obs, Grad, N, Sto, num_chain)

        normalize!(grad)
        θ.+= lr * grad
        θ.= abs.(θ)

        popfirst!(Es)
        popfirst!(norms)
        grad_list .= circshift(grad_list, (0, 1))

        push!(Es, E)
        push!(norms, norm(grad))
        grad_list[:, end] .= grad

        T0 = 6.222652893066406e-7 * N * L^2
        # println("Expected time is ", T0, " seconds")
        T += T0
        # println("total time is ", T/60/60, " hours")

        N = Int(floor(10^3 * 2^floor(i / 20)))
        lr = lr₀  * 0.4^floor(i / 20)

        # println("\nIteration :", i)
        # println("θ is :", θ)
        # println("N is :", N)
        # println("Δ is :", lr)
        # println("E is :", mean(Es)/L^2)
        # println("σ² is :", σ²/L^2)
        # println("M is :", mag/L^2)
        # println("gradient is :", grad)

        if N > 50 && mod(N, 5) ≡ 0
            res = vcat(L, h[:X], sum(θ), θ, mean(Es)/L^2, mag/L^2, std(Es)/L^2)'

            FILE = open(outfile * "/L-$L-p-$P-h-$hX.txt", "a")
            writedlm(FILE, res)
            close(FILE)
        end
    end

    return E
end

# nsgd(2, 1, 1, 10)
# @time nsgd(12, 1, 10^3, 160) # 20 minutes

# 200 is 20min
# 125 is 3 hours
# 160 is 7.3 hours


text:
export Chain
export popitems

"""
    Chain(blocks...)

Chain multiple blocks / functions together, so that they are called in sequence on a given input.
`Chain` also supports indexing and slicing, e.g. `c[2]` or `c[1:3]` or c[[1,3]].

# Examples
    julia> m = Chain(Dense(256,128), lstm(128,64), Dense(64,10));
"""
mutable struct Chain <: Block
    blocks::Vector
    function Chain(sequence::Vector)
        nblocks = length(sequence)
        blocks = Vector(undef,nblocks)
        for i = 1:nblocks
            blocks[i] = sequence[i]
        end
        new(blocks)
    end
    function Chain(sequence...)
        nblocks = length(sequence)
        blocks = Vector(undef,nblocks)
        for i = 1:nblocks
            blocks[i] = sequence[i]
        end
        new(blocks)
    end
    function Chain(nblocks::Int)
        new(Vector(undef,nblocks))
    end
end

Base.length(c::Chain)     = length(c.blocks)
Base.lastindex(c::Chain)  = length(c.blocks)
Base.firstindex(c::Chain) = 1
Base.getindex(c::Chain, k...)     =  c.blocks[k...]
Base.setindex!(c::Chain, v, k...) = (c.blocks[k...] = v)
Base.iterate(c::Chain, i=firstindex(c)) = i>length(c) ? nothing : (c[i], i+1)

function clone(this::Chain; type::Type=Array{Float32})
    nblocks = length(this.blocks)
    cloned  = Chain(nblocks)
    for i = 1:nblocks
        cloned[i] = clone(this[i], type=type)
    end
    return cloned
end

function Base.show(io::IO, c::Chain)
    print(io,  "Chain(\n")
    join(io,c.blocks, "\n")
    println(io, "\n)")
end

function popitems(blocks::Vector{T}, list) where T
    lenb = length(blocks)
    lenl = length(list)
    @assert(lenb>lenl)
    newblocks = Vector(undef,0)
    for i = 1:lenb
        if i ∉ list
            push!(newblocks,blocks[i])
        end
    end
    return newblocks
end

"""
    unbiasedof(model::Chain)

unbiased weights of Chain block
"""
function unbiasedof(model::Chain)
    weights = Vector(undef, 0)
    for m in model
        append!(weights, unbiasedof(m))
    end
    return weights
end

function paramsof(c::Chain)
    params = Vector{Variable}(undef,0)
    for i = 1:length(c)
        p = paramsof(c[i])
        if p ≠ nothing
            append!(params, p)
        end
    end
    return params
end

function xparamsof(c::Chain)
    xparams = Vector{XVariable}(undef,0)
    for i = 1:length(c)
        p = xparamsof(c[i])
        if p ≠ nothing
            append!(xparams, p)
        end
    end
    return xparams
end

function nparamsof(c::Chain)
    nparams = 0
    for i = 1:length(c)
        nparams += nparamsof(c[i])
    end
    return nparams
end

function bytesof(model::Chain, unit::String="MB")
    n = nparamsof(model) * elsizeof(model[1].w)
    return blocksize(n, uppercase(unit))
end

function resethidden(c::Chain)
    for i = 1:length(c)
        if typeof(c[i]) ∈ RNNLIST
            resethidden(c[i])
        end
    end
end

function forward(c::Chain, x::Variable)
    for i = 1:length(c)
        x = forward(c[i], x)
    end
    return x
end

function predict(c::Chain, x)
    for i = 1:length(c)
        x = predict(c[i], x)
    end
    return x
end

function to(type::Type, c::Chain)
    for m in c
        m = to(type, m)
    end
    return c
end

function to!(type::Type, c::Chain)
    for m in c
        m = to(type, m)
    end
    return nothing
end

max_stars_repo_name: sonosole/Delta

max_issues_repo_name: sonosole/Delta

max_forks_repo_name: sonosole/Delta


text:
########################################################################################
#
# Creates the DamerauLevenshtein metric
# The DamerauLevenshtein distance is the minimum number of operations (consisting of insertions,
# deletions or substitutions of a single character, or transposition of two adjacent characters)
# required to change one string into the other.
#
# Initial implementation based on http://blog.softwx.net/2015/01/optimizing-damerau-levenshtein_15.html
#
# This implementation based on https://github.com/matthieugomez/StringDistances.jl
#

function common_prefix(s1::String32, s2::String32)
    l = 1
    ls1 = length(s1)
    ls2 = length(s2)
    while l <= ls1 && l <= ls2
        s1[l] != s2[l] && break
        l += 1
    end
    return l - 1
end

function evaluate!(t1::String32, t2::String32, max_dist::Int, v0::Vector{Int}, v2::Vector{Int})
    if length(t1) > length(t2)
        t1, t2 = t2, t1
    end
    len1, len2 = length(t1), length(t2)
    len2 - len1 > max_dist && return max_dist + 1
    # prefix common to both strings can be ignored
    k = common_prefix(t1, t2)
    k == len1 && return len2 - k
    for i in 1:(len2 - k)
        v0[i] = i
    end

    offset = 1 + max_dist - (len2 - len1)
    i2_start = 1
    i2_end = max_dist

    i1 = 1
    current = i1
    prevch1 = t1[k + 1]
    for i in k+1:len1
        ch1 = t1[i]
        left = (i1 - 1)
        current = i1
        nextTransCost = 0
        prevch2 = t2[k + 1]
        i2_start += (i1 > offset) ? 1 : 0
        i2_end = min(i2_end + 1, len2)
        i2 = 1
        for j in k+1:len2
            @inbounds ch2 = t2[j]
            if i2_start <= i2 <= i2_end
                above = current
                thisTransCost = nextTransCost
                @inbounds nextTransCost = v2[i2]
                # cost of diagonal (substitution)
                @inbounds v2[i2] = current = left
                # left now equals current cost (which will be diagonal at next iteration)
                @inbounds left = v0[i2]
                if ch1 != ch2
                    # insertion
                    if left < current
                        current = left
                    end
                    # deletion
                    if above < current
                        current = above
                    end
                    current += 1
                    if (i1 != 1) & (i2 != 1) & (ch1 == prevch2) & (prevch1 == ch2)
                        thisTransCost += 1
                        if thisTransCost < current
                            current = thisTransCost
                        end
                    end
                end
                v0[i2] = current
            end
            i2 += 1
            prevch2 = ch2
        end
        @inbounds v0[i1 + len2 - len1] > max_dist && return max_dist + 1
        i1 += 1
        prevch1 = ch1
    end

    return current
end


text:
# NOTE: tests are sensitive to the line number at which statements appear
function f1(x, y)
    # A comment
    return x + y
end

@noinline function throws()
    x = nothing
    error("oops")
end
@inline inlined() = throws()
call_throws() = inlined()

f2(x, y) = x + y

@inline function multilinesig(x::Int,
                              y::String)
    z = x + 1
    return z
end

function f50()   # issue #50
    todB(x) = 10*log10(x)
    println("100x is $(todB(100)) dB.")
end

func_1st_nokwarg() = true
func_2nd_kwarg(; kw=2) = true

module Foo
module Bar
function fit end
end
end

function Foo.Bar.fit(m)
    return m
end

Foo.Bar.fit(a, b) = a + b

# Issue #81
function hasrettype(x::Real)::Float32
    return x*x + x
end


text:
export jats2pdf

using EzXML

function jats2pdfs(path::String)
    for dir in readdir(path)
        println(dir)
        filepath = joinpath(path, dir, dir)
        jats2pdf(filepath)
    end
end

function jats2pdf(path::String)
    xml = readxml("$path.xml")
    front = findfirst("article/front", xml)
    body = findfirst("article/body", xml)

    texts = readpdf("$path.pdf")
    texts = tokenize(texts)
    doc = PDDocument(texts)
    function f(xpath, node, label)
        nodes = findall(xpath, node)
        s = 1
        for n in nodes
            str = replace(nodecontent(n), " " => "")
            r = annotate!(doc, str, label, s)
            r == nothing || (s = last(r)+1)
        end
    end
    #f("front/journal-meta/journal-title-group/journal-title", "journal")
    f("article-meta/title-group/article-title", front, "title")
    f("article-meta/abstract", front, "abstract")
    f("sec/title", body, "section")

    open(GzipCompressorStream,"$path.pdf.txt.gz","w") do io
        lines = toconll(doc)
        for line in lines
            write(io, line, "\n")
        end
    end
end


text:
import Base.*
import Base.size
import Base.eltype
import LinearAlgebra.issymmetric

import LinearAlgebra.mul!

struct SqLinOp{Tv,Ti}
    issym::Bool
    value::Tv
    n::Ti
    multFn::Function #this is bad in terms of types?

    SqLinOp{Tv,Ti}(issym,value,n,multFn) where {Tv,Ti} = new(issym,value,n,multFn)
end

SqLinOp(issym,value::Tv,index::Ti,multFn) where {Tv,Ti} = SqLinOp{Tv,Ti}(issym,value,index,multFn)

eltype(A::SqLinOp{Tv,Ti}) where {Tv,Ti} = Tv

size(A::SqLinOp{Tv,Ti}, d::Ti) where {Tv,Ti} = A.n

size(A::SqLinOp{Tv,Ti}) where {Tv,Ti} = (A.n,A.n)

issymmetric(A::SqLinOp{Tv,Ti}) where {Tv,Ti} = A.issym

function *(A::SqLinOp{Tv,Ti}, b) where {Tv,Ti}
    return A.multFn(b)
end

# put in for compat with Arpack in Julia v0.7
function mul!(Y, A::SqLinOp{Tv,Ti}, B) where {Tv,Ti}
    Y1 = A*B
    for i in 1:A.n
        Y[i] = Y1[i]
    end
end

function testId(n::Ti) where Ti
    return M = SqLinOp{Float64,Int64}(true,1.0,n,x -> x)
end


text:
using Test
using MPI

MPI.Init()

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

# assuming there are at least two processes
x = rank == 0 ? collect(reshape(1.0:16.0, 4, 4)) : zeros(4, 4)

subarray = MPI.Type_Create_Subarray(2,
                                    Cint[4, 4],
                                    Cint[2, 2],
                                    Cint[0, 0],
                                    MPI.MPI_ORDER_FORTRAN,
                                    Float64)
MPI.Type_Commit!(subarray)

# test blocking send
if rank == 0
    MPI.Send(x, 1, subarray, 1, 0, comm)
elseif rank == 1
    MPI.Recv!(x, 1, subarray, 0, 0, comm)
    @test x == [1 5 0 0;
                2 6 0 0;
                0 0 0 0;
                0 0 0 0]
end

# test non blocking send
if rank == 0
    MPI.Isend(x, 1, subarray, 1, 0, comm)
elseif rank == 1
    req = MPI.Irecv!(x, 1, subarray, 0, 0, comm)
    MPI.Wait!(req)
    @test x == [1 5 0 0;
                2 6 0 0;
                0 0 0 0;
                0 0 0 0]
end

MPI.Finalize()


text:
function CalTrack(Depth,cal,bit;DepthFrom=minimum(Depth),DepthTo=maximum(Depth),Dtick=false)

    calax=gca()
    calax.plot(cal,Depth, color="red",linewidth=1, linestyle="--")
    calax.plot(bit,Depth, color="black",linewidth=0.5)
    calax.set_ylim([DepthTo,DepthFrom])
    calax.set_xlim([5,17])
    calax.set_xlabel("Diameter Inches")
    calax.set_yticks(range(DepthFrom, stop=DepthTo, length=51),minor=true)
    calax.set_yticks(range(DepthFrom, stop=DepthTo, length=11))
    calax.xaxis.tick_top()
    calax.xaxis.set_label_position("top")
    calax.tick_params("both",labelsize=8)
    setp(calax.get_yticklabels(),visible=Dtick)
    calax.fill_betweenx(Depth,cal,bit,where=(cal .< bit),color="yellow")
    calax.fill_betweenx(Depth,cal,bit,where=(cal .> bit),color=(0.2,0.4,0.8))
end


text:
## Root finding for Chebyshev expansions
#
#  Contains code that is based in part on Chebfun v5's chebfun/@chebteck/roots.m,
# which is distributed with the following license:

# Copyright (c) 2015, The Chancellor, Masters and Scholars of the University
# of Oxford, and the Chebfun Developers. All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
#     * Redistributions of source code must retain the above copyright
#       notice, this list of conditions and the following disclaimer.
#     * Redistributions in binary form must reproduce the above copyright
#       notice, this list of conditions and the following disclaimer in the
#       documentation and/or other materials provided with the distribution.
#     * Neither the name of the University of Oxford nor the names of its
#       contributors may be used to endorse or promote products derived from
#       this software without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR
# ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
# ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
# (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

function complexroots{C<:Chebyshev}(f::Fun{C})
    if ncoefficients(f)==0 || (ncoefficients(f)==1 && isapprox(f.coefficients[1],0))
        warn("Tried to take roots of a zero function.  Returning [].")
        Complex128[]
    elseif ncoefficients(f)==1
        Complex128[]
    elseif ncoefficients(f)==2
        Complex128[-f.coefficients[1]/f.coefficients[2]]
    else
        fromcanonical(f,colleague_eigvals(f.coefficients))
    end
end

#function roots(f::Fun)
#    irts=map(real,filter!(x->abs(x)<=1.+10eps(),filter!#(isreal,complexroots(f.coefficients))))
#
#    map!(x->x>1.?1.:x,irts)
#    map!(x->x<-1.?-1.:x,irts)
#
#    if length(irts)==0
#        Float64[]
#    else
#        fromcanonical(f,irts)
#    end
#end

function roots(f::Fun)
    f2=Fun(f,domain(f)) # default is to convert to Chebyshev/Fourier
    if space(f2)==space(f)
        error("roots not implemented for "*string(typeof(f)))
    else
        roots(f2)
    end
end

function roots{C<:Chebyshev}( f::Fun{C} )
# FIND THE ROOTS OF AN IFUN.

    d = domain(f)
    c = f.coefficients
    vscale = maxabs(values(f))
    if vscale == 0
        warn("Tried to take roots of a zero function.  Returning [].")
        return eltype(domain(f))[]
    end

    hscale = maximum( [abs(first(d)), abs(last(d))] )
    htol = eps(2000.)*max(hscale, 1)  # TODO: choose tolerance better

    if eltype(f) == BigFloat
        r = rootsunit_coeffs(convert(Vector{Float64},c./vscale), Float64(htol))
        # Map roots from [-1,1] to domain of f:
        rts = fromcanonical(d,r)
        fp = differentiate(f)
        while norm(f(rts)) > 1000eps(eltype(f))
            rts .-=f(rts)./fp(rts)
        end
    elseif eltype(f) == Complex{BigFloat}
        r = rootsunit_coeffs(convert(Vector{Complex{Float64}},c./vscale), Float64(htol))
        # Map roots from [-1,1] to domain of f:
        rts = fromcanonical(d,r)
        fp = differentiate(f)
        while norm(f(rts)) > 1000eps(eltype(f))
            rts .-=f(rts)./fp(rts)
        end
    else
        cvscale=c./vscale
        r = rootsunit_coeffs(cvscale, Float64(htol))

        # Check endpoints, as these are prone to inaccuracy
        # which can be deadly.
        if (isempty(r) || !isapprox(last(r),1.)) && abs(sum(cvscale)) < htol
            push!(r,1.)
        end
        if (isempty(r) || !isapprox(first(r),-1.)) && abs(alternatingsum(cvscale)) < htol
            insert!(r,1,-1.)
        end

        # Map roots from [-1,1] to domain of f:
        rts = fromcanonical(d,r)
    end
    return rts
end

function colleague_matrix{T<:Number}( c::Vector{T} )
#TODO: This is command isn't typed correctly
# COMPUTE THE ROOTS OF A LOW DEGREE POLYNOMIAL BY USING THE COLLEAGUE MATRIX:
    n = length(c) - 1
    A=zeros(T,n,n)

    for k=1:n-1
        A[k+1,k]=0.5
        A[k,k+1]=0.5
    end
    for k=1:n
        A[1,end-k+1]-=0.5*c[k]/c[end]
    end
    A[n,n-1]=one(T)

    # TODO: can we speed things up because A is lower-Hessenberg
    # Standard colleague matrix (See [Good, 1961]):
    A
end

function colleague_balance!(M)
    n=size(M,1)
    if n<2
        return M
    end

    d=1+abs(M[1,2])
    M[1,2]/=d;M[2,1]*=d

    for k=3:n
        M[k-1,k]*=d;M[k,k-1]/=d
        d+=abs(M[1,k])
        M[1,k]/=d;M[k-1,k]/=d;M[k,k-1]*=d
    end
    M
end

# colleague_eigvals( c::Vector )=hesseneigvals(colleague_balance!(colleague_matrix(c)))

colleague_eigvals( c::Vector )=eigvals(colleague_matrix(c))

function PruneOptions( r, htol::Float64 )
# ONLY KEEP ROOTS IN THE INTERVAL

    # Remove dangling imaginary parts:
    r = real( r[ abs(imag(r)) .< htol ] )
    # Keep roots inside [-1 1]:
    r = sort( r[ abs(r) .<= 1+htol ] )
    # Put roots near ends onto the domain:
    r = min( max( r, -1 ), 1)

    # Return the pruned roots:
    return r
end

rootsunit_coeffs{T<:Number}(c::Vector{T}, htol::Float64)=rootsunit_coeffs(c,htol,ClenshawPlan(T,Chebyshev(),length(c),length(c)))
function rootsunit_coeffs{S,T<:Number}(c::Vector{T}, htol::Float64,clplan::ClenshawPlan{S,T})
# Computes the roots of the polynomial given by the coefficients c on the unit interval.

    # If recursive subdivision is used, then subdivide [-1,1] into [-1,splitPoint] and [splitPoint,1].
    const splitPoint = -0.004849834917525

    # Simplify the coefficients by chopping off the tail:
    nrmc=norm(c, 1)
    @assert nrmc > 0  # There's an error if we make it this far
    c = chop(c,eps()*nrmc)
    n = length(c)

    if n == 0

        # EMPTY FUNCTION
        r = Float64[]

    elseif n == 1

        # CONSTANT FUNCTION
        r = ( c[1] == 0.0 ) ? Float64[0.0] : Float64[]

    elseif n == 2

        # LINEAR POLYNOMIAL
        r1 = -c[1]/c[2];
        r = ( (abs(imag(r1))>htol) | (abs(real(r1))>(1+htol)) ) ? Float64[] : Float64[max(min(real(r1),1),-1)]

    elseif n <= 70

        # COLLEAGUE MATRIX
        # The recursion subdividing below will keep going until we have a piecewise polynomial
        # of degree at most 50 and we likely end up here for each piece.

        # Adjust the coefficients for the colleague matrix
        # Prune roots depending on preferences:
        r = PruneOptions( colleague_eigvals(c), htol )::Vector{Float64}
    else

        #  RECURSIVE SUBDIVISION:
        # If n > 50, then split the interval [-1,1] into [-1,splitPoint], [splitPoint,1].
        # Find the roots of the polynomial on each piece and then concatenate. Recurse if necessary.

        # Evaluate the polynomial at Chebyshev grids on both intervals:
        #(clenshaw! overwrites points, which only makes sence if c is real)

        v1 = isa(c,Vector{Float64})?clenshaw!( c, points([-1,splitPoint],n),clplan):clenshaw( c, points([-1,splitPoint],n),clplan)
        v2 = isa(c,Vector{Float64})?clenshaw!( c, points([splitPoint,1] ,n),clplan):clenshaw( c, points([splitPoint,1] ,n),clplan)

        # Recurse (and map roots back to original interval):
        p = plan_chebyshevtransform( v1 )
        r = [ (splitPoint - 1)/2 + (splitPoint + 1)/2*rootsunit_coeffs( chebyshevtransform(v1,p), 2*htol,clplan) ;
                 (splitPoint + 1)/2 + (1 - splitPoint)/2*rootsunit_coeffs( chebyshevtransform(v2,p), 2*htol,clplan) ]

    end

    # Return the roots:
    r
end

function extremal_args{S<:PiecewiseSpace}(f::Fun{S})
    return cat(1,[extremal_args(fp) for fp in vec(f)]...)
end

function extremal_args(f::Fun)
    d=domain(f)
    if isa(d,PeriodicInterval)
        roots(differentiate(f))
    elseif isa(d,PeriodicDomain)  # avoid complex domains
        S=typeof(space(f))
        fromcanonical(f,extremal_args(setcanonicaldomain(f)))
    else
        dab=∂(domain(f))
        if ncoefficients(f) <=2 #TODO this is only relevant for Polynomial bases
            dab
        else
            [dab;roots(differentiate(f))]
        end
    end
end

for op in (:(Base.maximum),:(Base.minimum),:(Base.extrema),:(Base.maxabs),:(Base.minabs))
    @eval function $op{S<:RealSpace,T<:Real}(f::Fun{S,T})
        pts = extremal_args(f)

        $op(f(pts))
    end
end

for op in (:(Base.maxabs),:(Base.minabs))
    @eval function $op(f::Fun)
        # complex spaces/types can have different extrema
        pts = extremal_args(abs(f))

        $op(f(pts))
    end
end

for op in (:(Base.maximum),:(Base.minimum),:(Base.maxabs),:(Base.minabs))
    @eval begin
        $op{SV,DD<:UnionDomain,d,T<:Real}(f::Fun{PiecewiseSpace{SV,RealBasis,DD,d},T}) =
            $op(map($op,vec(f)))
    end
end

Base.extrema{SV,DD<:UnionDomain,d,T<:Real}(f::Fun{PiecewiseSpace{SV,RealBasis,DD,d},T}) =
    mapreduce(extrema,(x,y)->extrema([x...;y...]),vec(f))

for op in (:(Base.indmax),:(Base.indmin))
    @eval begin
        function $op{S<:RealSpace,T<:Real}(f::Fun{S,T})
            # the following avoids warning when differentiate(f)==0
            pts = extremal_args(f)
            # the extra real avoids issues with complex round-off
            pts[$op(real(f(pts)))]
        end

        function $op{S,T}(f::Fun{S,T})
            # the following avoids warning when differentiate(f)==0
            pts = extremal_args(f)
            fp=f(pts)
            @assert norm(imag(fp))<100eps()
            pts[$op(real(fp))]
        end
    end
end

for op in (:(Base.findmax),:(Base.findmin))
    @eval begin
        function $op{S<:RealSpace,T<:Real}(f::Fun{S,T})
            # the following avoids warning when differentiate(f)==0
            pts = extremal_args(f)
            ext,ind = $op(f(pts))
	        ext,pts[ind]
        end
    end
end

## Root finding for Laurent expansion

function companion_matrix{T}(c::Vector{T})
    n=length(c)-1

    if n==0
        zeros(T,0,0)
    else
        A=zeros(T,n,n)
        for k=1:n
            A[k,end]=-c[k]/c[end]
        end
        for k=2:n
            A[k,k-1]=one(T)
        end
        A
    end
end

if isdir(Pkg.dir("AMVW"))
    using AMVW
    function complexroots(coefficients::Vector)
        c=chop(coefficients,10eps())

        # Only use special routine for large roots
        if length(c)≥70
            Main.AMVW.rootsAMVW(c)
        else
            hesseneigvals(companion_matrix(c))
        end
    end
else
    complexroots(coefficients::Vector)=hesseneigvals(companion_matrix(chop(coefficients,10eps())))
end

complexroots(neg::Vector,pos::Vector)=complexroots([flipdim(chop(neg,10eps()),1);pos])
complexroots{DD}(f::Fun{Laurent{DD}})=mappoint(Circle(),domain(f),complexroots(f.coefficients[2:2:end],f.coefficients[1:2:end]))
complexroots{DD}(f::Fun{Taylor{DD}})=mappoint(Circle(),domain(f),complexroots(f.coefficients))

function roots{DD}(f::Fun{Laurent{DD}})
    irts=filter!(z->in(z,Circle()),complexroots(Fun(f.coefficients,Laurent(Circle()))))
    if length(irts)==0
        Complex{Float64}[]
    else
        rts=fromcanonical(f,tocanonical(Circle(),irts))
        if isa(domain(f),PeriodicInterval)
            sort!(real(rts))  # Make type safe?
        else
            rts
        end
    end
end

roots{D}(f::Fun{Fourier{D}})=roots(Fun(f,Laurent))

function roots{P<:PiecewiseSpace}(f::Fun{P})
    rts=mapreduce(roots,vcat,vec(f))
    k=1
    while k < length(rts)
        if isapprox(rts[k],rts[k+1])
            rts=rts[[1:k,k+2:end]]
        else
            k+=1
        end
    end

    rts
end

## Root finding for JacobiWeight expansion

# Add endpoints for JacobiWeight
# TODO: what about cancellation?
function roots{S<:JacobiWeight,T}(f::Fun{S,T})
    sp=space(f)
    d=domain(sp)
    rts=roots(Fun(f.coefficients,sp.space))
    if sp.α > 0
        rts=[first(d);rts]
    end
    if sp.β > 0
        rts=[rts;last(d)]
    end
    rts
end


text:
"""
    Pathfinding.position_delta(pathfinder::AStar{D}, from::NTuple{Int,D}, to::NTuple{Int,D})
Returns the absolute difference in coordinates between `from` and `to` taking into account
periodicity of `pathfinder`.
"""
position_delta(pathfinder::GridPathfinder{D,true}, from::Dims{D}, to::Dims{D}) where {D} =
    min.(abs.(to .- from), size(pathfinder.walkmap) .- abs.(to .- from))

position_delta(pathfinder::GridPathfinder{D,false}, from::Dims{D}, to::Dims{D}) where {D} =
    abs.(to .- from)

"""
    Pathfinding.delta_cost(pathfinder::GridPathfinder{D}, metric::M, from, to) where {M<:CostMetric}
Calculate an approximation for the cost of travelling from `from` to `to` (both of
type `NTuple{N,Int}`. Expects a return value of `Float64`.
"""
function delta_cost(
    pathfinder::GridPathfinder{D,periodic,true},
    metric::DirectDistance{D},
    from::Dims{D},
    to::Dims{D},
) where {D,periodic}
    delta = collect(position_delta(pathfinder, from, to))

    sort!(delta)
    carry = 0
    hdist = 0
    for i in D:-1:1
        hdist += metric.direction_costs[i] * (delta[D+1-i] - carry)
        carry = delta[D+1-i]
    end
    return hdist
end

function delta_cost(
    pathfinder::GridPathfinder{D,periodic,false},
    metric::DirectDistance{D},
    from::Dims{D},
    to::Dims{D},
) where {D,periodic}
    delta = position_delta(pathfinder, from, to)

    return sum(delta) * metric.direction_costs[1]
end

delta_cost(
    pathfinder::GridPathfinder{D},
    metric::MaxDistance{D},
    from::Dims{D},
    to::Dims{D},
) where {D} = max(position_delta(pathfinder, from, to)...)

delta_cost(
    pathfinder::GridPathfinder{D},
    metric::PenaltyMap{D},
    from::Dims{D},
    to::Dims{D},
) where {D} =
    delta_cost(pathfinder, metric.base_metric, from, to) +
    abs(metric.pmap[from...] - metric.pmap[to...])

delta_cost(pathfinder::GridPathfinder{D}, from::Dims{D}, to::Dims{D}) where {D} =
    delta_cost(pathfinder, pathfinder.cost_metric, from, to)


text:
__precompile__()

module Base58
export base58encode, base58decode, base58checkencode, base58checkdecode

import SHA: sha256

if VERSION < v"0.7.0-DEV.3213"
    codeunits(x) = convert(Array{UInt8}, x)
end

struct NotInAlphabetException <: Exception end

const BASE58CHARS = (
    codeunits("123456789ABCDEFGHJKLMNPQRSTUVWXYZabcdefghijkmnopqrstuvwxyz")...,
)
const SPACE = UInt8(' ')
const ZEROBASE58 = UInt8('1')

const REVBASE58CHARS = if VERSION < v"0.7.0-DEV.4455"
    ntuple(i -> findfirst(BASE58CHARS, UInt8(i)) |>
                x -> x == 0 ? typemax(UInt8) : UInt8(x - 1),
           typemax(UInt8))
else
    ntuple(i -> findfirst(isequal(UInt8(i)), BASE58CHARS) |>
                x -> x == nothing ? typemax(UInt8) : UInt8(x - 1),
           typemax(UInt8))
end

const BASE = 58

encode(x::UInt8) = @inbounds return BASE58CHARS[x + 1]

function base58encode(x::T) where T <: Union{DenseArray{UInt8, 1},
                                             NTuple{N, UInt8} where N}

    if length(x) == 0
        return UInt8[]
    end

    n_zeros = 0
    for c in x
        c != 0 && break
        n_zeros += 1
    end

    if n_zeros == length(x)
        return fill(ZEROBASE58, length(x))
    end

    length_result = ceil(Int64, length(x) * log(256) / log(BASE) + 1)

    res = zeros(UInt8, (length_result, ))
    l = 0

    i_x = 1
    @inbounds while i_x <= length(x)

        carry = UInt64(x[i_x])

        i = 0
        i_res = length_result
        while (carry != 0 || i < l) && i_res > 1
            carry += UInt64(0x100) * UInt64(res[i_res])
            carry, res[i_res] = divrem(carry, UInt64(BASE))

            i_res -= 1
            i += 1
        end

        @assert carry == 0
        l = i
        i_x += 1
    end

    i_res = 1
    @inbounds while i_res < length_result && res[i_res] == 0
        res[i_res] = 0
        i_res += 1
    end

    res = res[i_res - n_zeros:end]
    @inbounds for i in eachindex(res)
        res[i] = encode(res[i])
    end

    return res
end

function base58decode(x::T) where T <: Union{ DenseArray{UInt8, 1},
                                              NTuple{N, UInt8} where N }
    # Leading spaces will be ignored
    i = 1
    for xx in x
        if xx == SPACE
            i += 1
        else
            break
        end
    end

    n_zeros = 0
    @inbounds while i <= length(x) && x[i] == ZEROBASE58
        n_zeros += 1
        i += 1
    end

    out_size = ceil(Int64, (length(x) - n_zeros) * log(58) / log(256) + 1)
    res = zeros(UInt8, out_size)

    l = 0
    @inbounds while i <= length(x) && x[i] != SPACE

        carry = UInt64(REVBASE58CHARS[x[i]])
        if carry == UInt64(typemax(UInt8))
            throw(ArgumentError(
                "Letter " * Char(x[i]) * " not in Base58 alphabet"
            ))
        end

        j = 0
        k = out_size
        while (carry != 0 || j < l) && k > 0

            carry +=  UInt64(BASE) * res[k]
            carry, res[k] = divrem(carry, UInt64(256))

            k -= 1
            j += 1
        end

        @assert carry == 0

        l = i
        i += 1
    end

    @inbounds while i < length(x) && x[i] == SPACE
        i += 1
    end

    # The res = ...; return res is only for julia 0.6 compat
    @inbounds res = if out_size > 0

        i = 1
        while i <= out_size && res[i] == 0
            i += 1
        end

        if n_zeros > 0
            prepend!(res[i:end], zeros(UInt8, n_zeros))
        else
            res[i:end]
        end
    elseif n_zeros > 0
        zeros(UInt8, n_zeros)
    else
        res
    end
    return res
end

function base58checkencode(payload::T) where
    T <: Union{DenseArray{UInt8, 1}, NTuple{N, UInt8} where N}

    checksum = sha256(sha256(payload))[1:4]

    base58encode([payload..., checksum...])
end

function base58checkdecode(x::T, check::Bool = true) where
    T <: Union{DenseArray{UInt8, 1}, NTuple{N, UInt8} where N}

    dec = base58decode(x)
    payload = dec[1:end-4]
    checksum = dec[end-3:end]

    if check
        if sha256(sha256(payload))[1:4] != checksum
            throw(ArgumentError("Invalid address"))
        end
    end

    return payload
end

end # module Base58


text:
@testset "oettli-präger method" begin

    A = [2..4 -2..1; -1..2 2..4]
    b = [-2..2, -2..2]

    #p = solve(A, b, NonLinearOettliPrager())

    polyhedra = solve(A, b, LinearOettliPrager())

    for pnt in [[-4, -3], [3, -4], [4, 3], [-3, 4]]
        #@test any(pnt ∈ x for x in p.boundary)
        @test sum(pnt ∈ pol for pol in polyhedra) == 1
    end

    for pnt in [[-5, 5], [5, 5], [5, -5], [-5, -5]]
        @test all(pnt ∉ pol for pol in polyhedra)
    end
end


text:
#######################################
# IO
# Right now this basically just extends FileIO to support some things as text files

const TextFile = Union{File{format"TXT"}, File{format"SHA256"}}

function loadfile(::Type{T}, file::File) where T
    T(load(file))::T # Fallback to FileIO
end

function loadfile(::Type{T}, file::TextFile) where T  # specialize on T only to prevent ambiguity
    replace(read(file.filename, String), "\r"=>"") # ignore CRLF/LF difference
end

function loadfile(::Type{<:Number}, file::File{format"TXT"})
    parse(Float64, loadfile(String, file))
end

function savefile(file::File, content)
    save(file, content) # Fallback to FileIO
end

function savefile(file::TextFile, content)
    write(file.filename, string(content))
end

function query_extended(filename)
    file, ext = splitext(filename)
    # TODO: make this less hacky
    if uppercase(ext) == ".SHA256"
        res = File{format"SHA256"}(filename)
    else
        res = query(filename)
        if res isa File{DataFormat{:UNKNOWN}}
            res = File{format"TXT"}(filename)
        end
    end
    res
end

"""
    _convert(T::Type{<:DataFormat}, x; kw...) -> out

Convert `x` to a validate content for file data format `T`.
"""
_convert(::Type{<:DataFormat}, x; kw...) = x

# plain TXT
_convert(::Type{DataFormat{:TXT}}, x; kw...) = replace(string(x), "\r"=>"") # ignore CRLF/LF difference
_convert(::Type{DataFormat{:TXT}}, x::Number; kw...) = x
function _convert(::Type{DataFormat{:TXT}}, x::AbstractArray{<:AbstractString}; kw...)
    return join(x, '\n')
end
function _convert(
    ::Type{DataFormat{:TXT}}, img::AbstractArray{<:Colorant};
    size = (20,40), kw...)

    # encode image into string
    strs = @withcolor ImageInTerminal.encodeimg(
        ImageInTerminal.SmallBlocks(),
        ImageInTerminal.TermColor256(),
        img,
        size...)[1]
    return join(strs,'\n')
end

# SHA256
_convert(::Type{DataFormat{:SHA256}}, x; kw...) = bytes2hex(sha256(string(x)))
function _convert(::Type{DataFormat{:SHA256}}, img::AbstractArray{<:Colorant}; kw...)
    # encode image into SHA256
    size_str = bytes2hex(sha256(reinterpret(UInt8,[map(Int64,size(img))...])))
    img_str = bytes2hex(sha256(reinterpret(UInt8,vec(rawview(channelview(img))))))

    return size_str * img_str
end


text:
##########################################
# Some uncertain datasets
##########################################
UI = UncertainIndexDataset(uidxs)

##########################################
# Apply functions to datasets `n` times
##########################################
n = 3
@test resample(median, UI, n) isa Vector{T} where T <: Real
@test resample(cor, UI, UI, n) isa Vector{T} where T <: Real


text:
# ------------ #
# Original API #
# ------------ #

# splidef.m -- DONE
function splidef(breaks, evennum=0, k::Int=3)
    p = SplineParams(breaks, evennum, k)
    return (
        length(p.breaks),
        minimum(p.breaks),
        maximum(p.breaks),
        Any[p.breaks, 0, k],
        p
    )
end

# splinode.m  -- DONE
splinode(breaks::AbstractVector, evennum::Int, k::Int=3) =
    nodes(SplineParams(breaks, evennum, k))

function splidop(breaks, evennum=0, k=3, order=1)
    D, p = derivative_op(SplineParams(breaks, evennum, k), order)
    n = length(breaks) + k - 1

    D, n-order, breaks[1], breaks[end], Any[breaks, evennum, k-order]
end

# splibas.m -- DONE
function splibase(breaks::AbstractVector, evennum, k=3, x=splinode(breaks, evennum, k),
                  order=0)
    B = evalbase(SplineParams(breaks, evennum, k), x, order)
    B, x
end


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

function _recpol(x::T, y::T, degrees::Bool) where {T<:AbstractFloat}
    if degrees
        return hypot(x, y), rad2deg(atan(y, x))
    else
        return hypot(x, y), atan(y, x)
    end
end

"""
    recpol(x, y[, degrees=true]) -> radius, angle

### Purpose ###

Convert 2D rectangular coordinates to polar coordinates.

### Explanation ###

This is the partial inverse function of `polrec`.

### Arguments ###

* `x`: the abscissa coordinate of the point.  It may be a scalar or an array.
* `y`: the ordinate coordinate of the point.  It may be a scalar or an array
  of the same lenth as `x`.
* `degrees` (optional boolean keyword): if `true`, the output `angle` is given
 in degrees, otherwise in radians.  It defaults to `false`.

Mandatory arguments may also be passed as the 2-tuple `(x, y)`, so that it is
possible to execute `polrec(recpol(x, y))`.

### Output ###

A 2-tuple `(radius, angle)` with the polar coordinates of the input.  The
coordinate `angle` coordinate lies in the range \$[-\\pi, \\pi]\$ if
`degrees=false`, or \$[-180, 180]\$ when `degrees=true`.

If `x` and `y` are arrays, `radius` and `angle` are arrays of the same length as
`radius` and `angle`.

### Example ###

Calculate polar coordinates \$(r, \\varphi)\$ of point with rectangular
coordinates \$(x, y) = (2.24, -1.87)\$.

```jldoctest
julia> using AstroLib

julia> r, phi = recpol(2.24, -1.87)
(2.917961617293826, -0.6956158538564537)
```

Angle \$\\varphi\$ is given in radians.

"""
recpol(x::Real, y::Real; degrees::Bool=false) =
    _recpol(promote(float(x), float(y))..., degrees)

recpol(xy::Tuple{Real, Real}; degrees::Bool=false) =
    recpol(xy..., degrees=degrees)

function recpol(x::AbstractArray{X}, y::AbstractArray{Y};
                degrees::Bool=false) where {X<:Real, Y<:Real}
    @assert length(x) == length(y)
    typex = float(X)
    r = similar(x, typex)
    a = similar(x, typex)
    for i in eachindex(x)
        r[i], a[i] = recpol(x[i], y[i], degrees=degrees)
    end
    return r, a
end


text:
"""
    eigenSys!(h, qcols, neq)

Compute the roots and the left eigenvectors of the companion
matrix, sort the roots from large-to-small, and sort the
eigenvectors conformably.  Map the eigenvectors into the real
domain. Count the roots bigger than uprbnd.
"""
function eigenSys!(aa::Array{Float64,2}, upperbound::Float64, rowsLeft::Int64)

    roots = Array{Complex{Float64},1}(eigvals(copy(transpose(aa))))
    ww = eigvecs(copy(transpose(aa)))

    # sort eigenvalues in descending order of magnitude
    magnitudes = abs.(roots)
    highestToLowestMag = sortperm(magnitudes, rev = true) # reverse order
    roots = roots[highestToLowestMag]

    # sort eigenvecs in order found above
    ww = ww[:, highestToLowestMag]

    #  Given a complex conjugate pair of vectors W = [w1,w2], there is a
    #  nonsingular matrix D such that W*D = real(W) + imag(W).  That is to
    #  say, W and real(W)+imag(W) span the same subspace, which is all
    #  that AMA cares about.

    ww = ( real(ww) + imag(ww) )

    # count how many roots are above the upperbound threshold
    lgroots = mapreduce((mag->mag > upperbound ? 1 : 0), +, magnitudes)

    return (ww, roots, lgroots)

end # eigenSys!


text:
using PerfChecker
using Test

using PatternFolds

@testset "PatternFolds.jl" begin
    # Title of the alloc check (for logging purpose)
    title = "Basic intervals and vectors folds operation"

    # Dependencies needed to execute pre_alloc and alloc
    dependencies = [PatternFolds]

    # Target of the alloc check
    targets = [PatternFolds]

    function alloc() # 0.2.x
        # Intervals
        itv = Interval{Open,Closed}(0.0, 1.0)
        i = IntervalsFold(itv, 2.0, 1000)

        @info "Checking IntervalsFold" i pattern(i) gap(i) folds(i) size(i) length(i)

        unfold(i)
        collect(i)
        reverse(collect(i))

        # rand(i, 1000)

        # Vectors
        vf = make_vector_fold([0, 1], 2, 1000)
        @info "Checking VectorFold" vf pattern(vf) gap(vf) folds(vf) length(vf)

        unfold(vf)
        collect(vf)
        reverse(collect(vf))

        rand(vf, 1000)

        return nothing
    end

    # function alloc() #0.1.1-0.1.5
    #     i = IntervalsFold(Interval((0.0, true), (1.0, false)), 2.0, 1000)

    #     @info "Checking IntervalsFold" i pattern(i) gap(i) folds(i) size(i) length(i)

    #     unfold(i)
    #     collect(i)
    #     reverse(collect(i))

    #     # rand(i, 1000)

    #     # Vectors
    #     vf = VectorFold([0, 1], 2, 1000)
    #     @info "Checking VectorFold" vf pattern(vf) gap(vf) folds(vf) length(vf)

    #     unfold(vf)
    #     collect(vf)
    #     reverse(collect(vf))

    #     [rand(vf) for _ in 1:1000]
    # end

    # Actual call to PerfChecker
    alloc_check(title, dependencies, targets, alloc, alloc; path=@__DIR__)
end


text:
importall JuMP, Gurobi, LightGraphs, Polynomials, Roots, Combinatorics
using Plots
pyplot()

function extract_data(model::JuMP.Model)
    c = JuMP.prepAffObjective(model)

    xlb = copy(model.colLower)
    xub = copy(model.colUpper)

    return c, xlb, xub
end

function lagrange(m::JuMP.Model)
  c, xlb, xub = extract_data(m)
  c = float(c)
  c_ = copy(c)
  number_of_edges = length(c)
  x = variable(Int)
  n = Int(fzero(x^2 - x - 2*number_of_edges, 1.0))

  ind_fixas_0 = find(xub + xlb .== 0)
  ind_fixas_1 = find(xub + xlb .== 2)

  c[ind_fixas_0] = Inf
  c[ind_fixas_1] = -Inf

  distmx = cost_matrix(c,n)
  distmx_ = cost_matrix(c_,n)

  g = Graph()
  add_vertices!(g,n)
  e = create_edges(n,number_of_edges)

  for i=1:number_of_edges
    if i ∈ ind_fixas_1
      add_edge!(g,e[i,1],e[i,2])
    end
  end

  for i=1:n
    δ = degree(g,i)
    if δ > 2
      m.objVal = :Infeasible
      m.colVal = NaN
      return nothing
    end
  end
  #=
 if ind_fixas_1 !=[]
   #S = vcat([collect(combinations(2:n,i)) for i=1:n]...)
   Subsets = subsets(distmx_)
   s = size(Subsets,2)
   for i=1:s
     if Subsets[1,i] == 0
       S = pertence(Subsets[:,i])
       if length(S)>3
         sub_graph = induced_subgraph(g,S)[1]
         if ne(sub_graph) > nv(sub_graph) - 1
           m.objVal = :Infeasible
           m.colVal = NaN
           return nothing
         end
       end
     end
   end
 end
=#
  for i=1:number_of_edges
    if i ∉ ind_fixas_0
      add_edge!(g,e[i,1],e[i,2])
    end
  end

  V = copy(neighbors(g,1))
  for v in V
    rem_edge!(g,1,v)
  end

  ω = 0.8(initial_tour(n,distmx)[1])
  u = zeros(n)
  iter = 1
  max_iter = 500
  x = []
  z = []
  converge = zeros(max_iter + 1)
  y = Array{Float64,1}(n)
  c_e = Array{Float64,1}(number_of_edges)
  c_e_ = Array{Float64,1}(number_of_edges)
  while iter <= max_iter
    for i=1:number_of_edges
        c_e[i] = c[i] - u[e[i,1]] - u[e[i,2]]
        c_e_[i] = c[i] - u[e[i,1]] - u[e[i,2]]
    end

    distmx = cost_matrix(c_e,n)
    distmx_ = cost_matrix(c_e_,n)

    tree = kruskal_mst(g,distmx)
    cost = []
    for i=2:n
      cost_ = (distmx_[1,i],[1,i])
      cost = copy([cost;[cost_]])
    end
    sort!(cost, by = x -> x[1])

    g_tree = Graph()

    add_vertices!(g_tree, nv(g))

    for i in 1:length(tree)
      add_edge!(g_tree,tree[i])
    end

    k = 0
    for j=1:n - 1
      if j ∈ ind_fixas_1
         add_edge!(g_tree, Edge(cost[j][2][1],cost[j][2][2]))
         k = k +1
         if k == 2
           break
         end
      end
    end
    if k !=2
      for j=1:n - 1
        if j ∉ ind_fixas_0
           add_edge!(g_tree, Edge(cost[j][2][1],cost[j][2][2]))
           k = k +1
           if k == 2
             break
           end
        end
      end
    end

    for i=1:n
      y[i] = 2 - degree(g_tree,i)
    end
    z = (sum(full(adjacency_matrix(g_tree)).*distmx_))/2 +  2*sum(u)
    x = vcat([diag(full(adjacency_matrix(g_tree)),i) for i=1:n-1]...)
    if y == zeros(n)
      break
    end
    #=
    μ = 1/iter
    for i=1:n
      u[i] = minimum([u[i] + μ*y[i];0])
    end
    =#
    u =  u + ((ω - z)/(sum(y[i]^2 for i=1:n)))*y
    iter = iter + 1
    converge[iter] = z

  end
  m.objVal = copy(z)
  return nothing
end

function create_edges(n,number_of_edges) # Cria as arestas a partir da Matriz de Adjacencia

  e = zeros(number_of_edges,2)
  k=1;
  for i=1:n
    for j=i+1:n
      e[k,1]=i
      e[k,2]=j
      k = k+1
    end
  end
  return round(Int64,e)
end

function initial_tour(n,distmx)
  tour = Graph()
  add_vertices!(tour,n)
  t  = ne(tour)
  i = 1
  j = []
  visit = []
  while t < n - 1
    push!(visit,i)
    cost_ = copy(distmx[i,:])
    cost = copy(cost_)
    deleteat!(cost_,sort!(unique(visit)))
    val = findmin(cost_)[1]
    j = find(x->x==val,cost)[1]
    add_edge!(tour,i,j)
    i = copy(j)
    t = ne(tour)
  end
  add_edge!(tour,j,1)
   ω = sum(distmx.*full(adjacency_matrix(tour)))/2
   x = vcat([diag(full(adjacency_matrix(tour)),i) for i=1:n-1]...)
   return ω, x
end

function cost_matrix(c,n)
distmx = zeros(n,n)
  t = 0
  for i=1:n
      for j=i+1:n
        t = t+1
        distmx[i,j]=c[t]
      end
  end
  distmx = distmx + distmx'
  return distmx
end

function subsets(Matrix) # Retorna todos os subconjuntos de vértices de uma Matriz de Adjacência
  number_of_nodes = size(Matrix,1)
  Subset_S = digits(1,2,number_of_nodes)
  for i = 2:2^(number_of_nodes)-2
    aux = digits(i,2,number_of_nodes)
    Subset_S = hcat(Subset_S,aux)
  end
  return Subset_S
end
function pertence(x)
  n = length(x)
  nodes = []
  for i=1:n
    if x[i] == 1
      nodes = copy([nodes;[i]])
    end
  end
  return nodes
end


text:
if !isdir("ccv")
    run(`git clone https://github.com/liuliu/ccv`)
end
cd("ccv")
run(`git pull`)
cd("lib")
mv("makefile", "makefile.original", remove_destination=true) # we have our own
cp("../../makefile", "./makefile")
run(`make clean`)
run(`./configure`)
run(`make`)


text:
module metric

function accuracy(y, pred_proba)::Float64
    pred_proba = Array(pred_proba)
    y = Array(y)
    y_pred = mapslices(argmax, pred_proba, dims = 1)
    y_true = mapslices(argmax, y, dims = 1)

    return sum(y_true .== y_pred) / length(y_true)
end

end

max_stars_repo_name: prenc/MiraculousNN

max_issues_repo_name: prenc/MiraculousNN

max_forks_repo_name: prenc/MiraculousNN


text:
"""
`module LinearMapsAA`

Provides `AbstractArray` (actually `AbstractMatrix`)
or "Ann Arbor" version of `LinearMap` objects
"""
module LinearMapsAA

export LinearMapAA

# Indexer = AbstractVector{Int}

include("lm-aa.jl")

end # module


text:
import JSON
import Base: writemime
import MCTS: node_tag, tooltip_tag

function MCTS.create_json{P<:RobustMCTSPlanner}(visualizer::TreeVisualizer{P})
    root_id = -1
    next_id = 1
    node_dict = Dict{Int, Dict{UTF8String, Any}}()
    s_dict = Dict{Any, Int}()
    sa_dict = Dict{Any, Int}()
    san_dict = Dict{Any, Int}()
    for (s, sn) in visualizer.policy.tree
        # create state node
        node_dict[next_id] = sd = Dict("id"=>next_id,
                                       "type"=>:state,
                                       "children_ids"=>Array(Int,0),
                                       "tag"=>node_tag(s),
                                       "tt_tag"=>tooltip_tag(s),
                                       "N"=>sn.N
                                       )
        if s == visualizer.init_state
            root_id = next_id
        end
        s_dict[s] = next_id
        next_id += 1

        # create action nodes
        for (a, san) in sn.A
            node_dict[next_id] = sad = Dict("id"=>next_id,
                                            "type"=>:action,
                                            "children_ids"=>Array(Int,0),
                                            "tag"=>node_tag(a),
                                            "tt_tag"=>tooltip_tag(a),
                                            "N"=>san.N,
                                            "Q"=>san.Q
                                            )
            push!(sd["children_ids"], next_id)
            sa_dict[(s,a)] = next_id
            next_id += 1

            # create nature nodes
            for (mdp, nnode) in san.nature_nodes
                node_dict[next_id] = Dict("id"=>next_id,
                                          "type"=>:action,
                                          "children_ids"=>Array(Int,0),
                                          "tag"=>node_tag(mdp),
                                          "tt_tag"=>tooltip_tag(mdp),
                                          "N"=>nnode.N,
                                          "Q"=>nnode.Q
                                          )
                push!(sad["children_ids"], next_id)
                san_dict[(s,a,mdp)] = next_id
                next_id += 1
            end
        end
    end

    if root_id < 0
        error("""
                MCTS tree visualization: Policy does not have a node for the specified state.
            """)
    end

    # go back and refill transitions
    for (s, sn) in visualizer.policy.tree
        for (a, san) in sn.A
            for (mdp, nnode) in san.nature_nodes
                for sp in keys(nnode.transitions)
                    nd = node_dict[san_dict[(s,a,mdp)]]
                    if haskey(s_dict, sp)
                        push!(nd["children_ids"], s_dict[sp])
                    else
                        node_dict[next_id] = Dict("id"=>next_id,
                                                  "type"=>:state,
                                                  "children_ids"=>Array(Int,0),
                                                  "tag"=>node_tag(sp),
                                                  "tt_tag"=>tooltip_tag(sp),
                                                  "N"=>0
                                                  )
                        push!(nd["children_ids"], next_id)
                        next_id += 1
                    end
                end
            end
        end
    end
    json = JSON.json(node_dict)
    return (json, root_id)
end


text:
@testset "583.delete-operation-for-two-strings.jl" begin
    @test min_distance583("sea", "eat") == 2
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMBase.jl/blob/master/LICENSE

"""
    FEMBase.Test

FEMBase.Test contains utilities to test JuliaFEM extensions, including:

- all test macros from Base.Test (@test, @testset, ...)
- test problems `Poisson` and `Dirichlet`.
- test linear system solver `LUSolver`.
- test analysis `Static`.
- read ABAQUS matrices

"""
module Test

# For convenience, one can also get all @test macros and so on just by typing
# `using FEMBase.Test`, they are the very same than in `Base.Test`
using Base.Test
export @test, @test_throws, @test_broken, @test_skip,
       @test_warn, @test_nowarn
# export @test_logs, @test_deprecated
export @testset, @inferred

using FEMBase
import FEMBase: assemble_elements!, get_unknown_field_name, solve!, run!

mutable struct Poisson <: FieldProblem end

function get_unknown_field_name(::Problem{Poisson})
    return "u"
end

function assemble_elements!(problem::Problem{Poisson},
                            assembly::Assembly,
                            elements::Vector{Element{E}},
                            time::Float64) where E

    bi = BasisInfo(E)
    ndofs = length(bi)
    Ke = zeros(ndofs, ndofs)
    fe = zeros(ndofs)

    for element in elements
        fill!(Ke, 0.0)
        fill!(fe, 0.0)
        for ip in get_integration_points(element)
            J, detJ, N, dN = element_info!(bi, element, ip, time)
            k = 1.0
            if haskey(element, "coefficient")
                k = element("coefficient", ip, time)
            end
            Ke += ip.weight * k*dN'*dN * detJ
            if haskey(element, "source")
                f = element("source", ip, time)
                fe += ip.weight * N'*f * detJ
            end
        end
        gdofs = get_gdofs(problem, element)
        add!(assembly.K, gdofs, gdofs, Ke)
        add!(assembly.f, gdofs, fe)
    end
end

function assemble_elements!(problem::Problem{Poisson},
                            assembly::Assembly,
                            elements::Vector{Element{E}},
                            time::Float64) where E<:Union{Seg2,Seg3}

    bi = BasisInfo(E)
    ndofs = length(bi)
    Ce = zeros(ndofs, ndofs)
    fe = zeros(ndofs)
    ge = zeros(ndofs)

    for element in elements
        fill!(Ce, 0.0)
        fill!(fe, 0.0)
        fill!(ge, 0.0)
        for ip in get_integration_points(element)
            J, detJ, N, dN = element_info!(bi, element, ip, time)
            if haskey(element, "flux")
                f = element("flux", ip, time)
                fe += ip.weight * N'*f * detJ
            end
            if haskey(element, "fixed u")
                f = element("fixed u", ip, time)
                Ce += ip.weight * N'*N * detJ
                ge += ip.weight * N'*f * detJ
            end
        end
        gdofs = get_gdofs(problem, element)
        add!(assembly.C1, gdofs, gdofs, Ce')
        add!(assembly.C2, gdofs, gdofs, Ce)
        add!(assembly.f, gdofs, fe)
        add!(assembly.g, gdofs, ge)
    end
end

mutable struct Dirichlet <: BoundaryProblem end

function assemble_elements!(problem::Problem{Dirichlet},
                            assembly::Assembly,
                            elements::Vector{Element{E}},
                            time::Float64) where E

    name = get_parent_field_name(problem)
    dim = get_unknown_field_dimension(problem)

    data = Dict{Int64,Float64}()
    for element in elements
        for i=1:dim
            haskey(element, "$name $i") || continue
            gdofs = get_gdofs(problem, element)
            ldofs = gdofs[i:dim:end]
            xis = get_reference_element_coordinates(E)
            for (ldof, xi) in zip(ldofs, xis)
                data[ldof] = interpolate(element, "$name $i", xi, time)
            end
        end
    end

    for (dof, val) in data
        add!(assembly.C1, dof, dof, 1.0)
        add!(assembly.C2, dof, dof, 1.0)
        add!(assembly.g, dof, val)
    end

end

mutable struct LUSolver <: AbstractLinearSystemSolver
end

"""
    solve!(solver::LUSolver, ls::LinearSystem)

Solve linear system using LU factorization. If final system has zero rows,
add 1 to diagonal to make matrix non-singular.
"""
function solve!(ls::LinearSystem, ::LUSolver)
    A = [ls.K ls.C1; ls.C2 ls.D]
    b = [ls.f; ls.g]

    # add 1.0 to diagonal for any zero rows in system
    p = ones(2*ls.dim)
    p[unique(rowvals(A))] = 0.0
    A += spdiagm(p)

    # solve A*x = b using LU factorization and update solution vectors
    x = lufact(A) \ full(b)
    ls.u = x[1:ls.dim]
    ls.la = x[ls.dim+1:end]

    return nothing
end

mutable struct Static <: AbstractAnalysis
    time :: Float64
end

function Static()
    return Static(0.0)
end

"""
    get_linear_system(problems, time)

Assemble `problems` and construct LinearSystem for test problem.
"""
function get_linear_system(problems, time)
    # assemble matrices for all problems
    N = 0 # size of resulting matrix
    for problem in problems
        assemble!(problem, time)
        N = max(N, size(problem.assembly.K, 2))
    end
    # create new LinearSystem and add assemblies to that
    ls = LinearSystem(N)
    for problem in problems
        ls.M += sparse(problem.assembly.M, N, N)
        ls.K += sparse(problem.assembly.K, N, N)
        ls.Kg += sparse(problem.assembly.K, N, N)
        ls.f += sparse(problem.assembly.f, N, 1)
        ls.fg += sparse(problem.assembly.f, N, 1)
        ls.C1 += sparse(problem.assembly.C1, N, N)
        ls.C2 += sparse(problem.assembly.C2, N, N)
        ls.D += sparse(problem.assembly.D, N, N)
        ls.g += sparse(problem.assembly.g, N, 1)
    end
    return ls
end

"""
    run!(analysis::Analysis{Static})

Run static analysis for test problem. This is for testing purposes only and
should not be used in real simulations.
"""
function run!(analysis::Analysis{Static})
    info("This is from FEMBase.Test.Analysis{Static}, which is for testing ",
         "purposes only.")
    ls = get_linear_system(get_problems(analysis), analysis.properties.time)
    solve!(ls, LUSolver())
    normu = norm(ls.u)
    normla = norm(ls.la)
    info("Solution norms are |u| = $normu, |la| = $normla")
    return ls, normu, normla
end

export Poisson, Dirichlet, Static, LUSolver, get_linear_system, abaqus_read_mtx

function read_mtx(data::IO; dim=0)
    if dim == 0
        for ln in eachline(data)
            i,idof,j,jdof,value = map(parse, split(ln, ','))
            dim = max(dim, idof, jdof)
        end
    end
    seekstart(data)
    I = Int64[]
    J = Int64[]
    V = Float64[]
    for ln in eachline(data)
        i,idof,j,jdof,value = map(parse, split(ln, ','))
        if i < 1 || j < 1
            continue
        end
        push!(I, (i-1)*dim+idof)
        push!(J, (j-1)*dim+jdof)
        push!(V, value)
    end
    A = full(sparse(I, J, V))
    A += transpose(tril(A,-1))
    return A
end

"""
    read_mtx_from_file(filename::String; dim=0)

Read matrix data export from ABAQUS by using command:

*STEP
*MATRIX GENERATE, STIFFNESS, MASS, LOAD
*MATRIX OUTPUT, STIFFNESS, MASS, LOAD
*END STEP

If number of dofs / node (dim) is not given, it will be automatically
detemined from file.

"""
function read_mtx_from_file(filename::String; dim=0)
    open(filename) do fid
        return read_mtx(fid; dim=dim)
    end
end

function read_mtx_from_string(data::String; dim=0)
    return read_mtx(IOBuffer(data); dim=dim)
end

"""
    test_resource(resource::String) -> String

`@test_resource(resource)` expands to a string containing full path to the some
`resource` what is needed by a test.

# Example

Say, test file name is `test_run_model.jl` inside package called `Models`.
One needs a test mesh file `mesh.inp`. Then, function

```julia
mesh_file = @test_resource("mesh.inp")
```

expands to

    `~/.julia/v0.6/Models/test/test_run_model/mesh.inp`

Macro must be used inside the actual test file.
"""
macro test_resource(resource::String)
    #__source__.file === nothing && return nothing
    #filename = String(__source__.file)
    filename = Base.source_path()
    abs_dirname = abspath(dirname(filename))
    test_dir_name = first(splitext(basename(filename)))
    resource_ = joinpath(abs_dirname, test_dir_name, resource)
    if !isfile(resource_)
        warn("Cannot find resource $resource_")
    end
    return resource_
end

export read_mtx_from_file, read_mtx_from_string, @test_resource
end


text:
module JungleHelperDropLanternTrigger

using ..Ahorn, Maple

@mapdef Trigger "JungleHelper/DropLanternTrigger" DropLanternTrigger(x::Integer, y::Integer, width::Integer=8, height::Integer=8, oneUse::Bool=false, destroyLantern::Bool=true)

const placements = Ahorn.PlacementDict(
    "Drop Lantern (Jungle Helper)" => Ahorn.EntityPlacement(
       DropLanternTrigger,
        "rectangle"
    )
)

end

max_stars_repo_name: Jackal-Celeste/JungleHelper

max_issues_repo_name: Jackal-Celeste/JungleHelper

max_forks_repo_name: Jackal-Celeste/JungleHelper


text:
module A

using Typeclass

@class Copyable T Symbol begin
    clone(::T)::T
    obtain(::T, ::Symbol)
end

type Me
    self::Int
end

@instance Copyable Me begin
    clone(x::Me) = Me(x.self + 1)
    obtain(x::Me, s::Symbol) = s
end

export Me, clone, obtain, Copyable
end

using A
import A: clone, obtain
using Typeclass

x = Me(1)
println(clone(x))

type Me2
    self::Int
end
@instance Copyable Me2 begin
    function clone(x::Me2)
        Me(x.self + 100)
    end
    function obtain(x::Me2, s::Symbol)
        s
    end
end

x2 = Me2(1)
println(clone(x2))

println(typeof(Copyable))
println(obtain(x, :prop))
println(obtain(x2, :prop2))


text:
using Documenter, DeepViz

makedocs(;
    modules=[DeepViz],
    format=Documenter.HTML(),
    pages=[
        "About" => "index.md",
        "Getting Started" => "gs.md",
        "Insights" => "insights.md",
        "API Reference" => "api.md",
        "Utility Function" => "utility.md",
        "Usage" => "usage.md",
        "FAQs" => "faq.md"
    ],
    repo="https://github.com/AdarshKumar712/DeepViz.jl/blob/{commit}{path}#L{line}",
    sitename="DeepViz.jl",
    authors="Adarshkumar712 <Adarshkumar712.ak@gmail.com>",
    assets=String[],
)

deploydocs(;
    repo="github.com/AdarshKumar712/DeepViz.jl",
)


text:
using Pkg

Pkg.activate("../..")

using ColonyRuntime
using ColonyRuntime.Crypto

function main(args)
  colony_prvkey = args[1]
  colonyid = Crypto.id(colony_prvkey)
  server = ColonyRuntime.ColoniesServer("localhost", 8080)

  # register a runtime
  println("- registering a new runtime to colony " * colonyid)
  runtime_prvkey = Crypto.prvkey()
  println("  runtime_prvkey: ", runtime_prvkey)
  println("  runtimeid: ", Crypto.id(runtime_prvkey))

  runtime = ColonyRuntime.Runtime(Crypto.id(runtime_prvkey), "fibonacci_generator", "fibonacci_generator", colonyid, "AMD Ryzen 9 5950X (32) @ 3.400GHz", 32, 80326, "NVIDIA GeForce RTX 2080 Ti Rev. A", 1)
  runtime = ColonyRuntime.addruntime(server, runtime, colony_prvkey)

  # and approve it so that can use the api
  println("- approving runtime " * runtime.runtimeid)
  ColonyRuntime.approveruntime(server, runtime.runtimeid, colony_prvkey)

  # submit a process spec
  println("- submitting process spec, fibonacci_num=" * args[2] * ", target_runtime=fibonacci_solver")

  conditions = ColonyRuntime.Conditions(colonyid, [], "fibonacci_solver", 1, 1024, 0) # 1 core, 1024 Mib memory, no GPU
  env = Dict()
  env["fibonacci_num"] = args[2]

  processpec = ColonyRuntime.ProcessSpec(-1, -1, conditions, env)
  process = ColonyRuntime.submitprocess(server, processpec, runtime_prvkey)
  println("  processid: ", process.processid)
end

main(ARGS)


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

# The numeric type used to determine how the different
# versions of a package should be weighed
struct VersionWeight
    major::Int64
    minor::Int64
    patch::Int64
end
VersionWeight(major::Integer, minor::Integer) = VersionWeight(major, minor, 0)
VersionWeight(major::Integer) = VersionWeight(major, 0)
VersionWeight() = VersionWeight(0)
VersionWeight(vn::VersionNumber) = VersionWeight(vn.major, vn.minor, vn.patch)

Base.zero(::Type{VersionWeight}) = VersionWeight()

Base.typemin(::Type{VersionWeight}) = (x=typemin(Int64); VersionWeight(x, x, x))

Base.:(-)(a::VersionWeight, b::VersionWeight) =
    VersionWeight(a.major-b.major, a.minor-b.minor, a.patch-b.patch)

Base.:(+)(a::VersionWeight, b::VersionWeight) =
    VersionWeight(a.major+b.major, a.minor+b.minor, a.patch+b.patch)

Base.:(-)(a::VersionWeight) =
    VersionWeight(-a.major, -a.minor, -a.patch)

function Base.cmp(a::VersionWeight, b::VersionWeight)
    c = cmp(a.major, b.major); c != 0 && return c
    c = cmp(a.minor, b.minor); c != 0 && return c
    return cmp(a.patch, b.patch)
end
Base.isless(a::VersionWeight, b::VersionWeight) = cmp(a,b) < 0
Base.:(==)(a::VersionWeight, b::VersionWeight) = cmp(a,b) == 0

Base.abs(a::VersionWeight) =
    VersionWeight(abs(a.major), abs(a.minor), abs(a.patch))

Base.copy(a::VersionWeight) =
    VersionWeight(a.major, a.minor, a.patch)

# This isn't nice, but it's for debugging only anyway
function Base.show(io::IO, a::VersionWeight)
    print(io, "(", a.major)
    a == VersionWeight(a.major) && @goto done
    print(io, ".", a.minor)
    a == VersionWeight(a.major, a.minor) && @goto done
    print(io, ".", a.patch)
    @label done
    print(io, ")")
end


text:
"""
    SChunk

This struct is the standard container for Blosc 2 compressed data.
"""
mutable struct SChunk
    ptr ::Ptr{Lib.blosc2_schunk}
    function SChunk(schunk::Ptr{Lib.blosc2_schunk})
        result = new(schunk)
        finalizer(result) do sc
            Lib.blosc2_schunk_free(sc.ptr)
        end
    end
end

"""
    clevel(s::SChunk)

Compression level of schunk
"""
clevel(s::SChunk) = unsafe_load(s.ptr.clevel)
"""
    typesize(s::SChunk)

Type Size  of schunk
"""
typesize(s::SChunk) = unsafe_load(s.ptr.typesize)
"""
    chunksize(s::SChunk)

Size of chunkes of super-chunk
"""
chunksize(s::SChunk) = unsafe_load(s.ptr.chunksize)
"""
    blocksize(s::SChunk)

Block Size of schunk
"""
blocksize(s::SChunk) = unsafe_load(s.ptr.blocksize)
"""
    nchunks(s::SChunk)

Number of chunks in super-chunk.
"""
nchunks(s::SChunk) = unsafe_load(s.ptr.nchunks)
"""
    nbytes(s::SChunk)

Uncompressed size of all data in super-chunk
"""
nbytes(s::SChunk) = unsafe_load(s.ptr.nbytes)
"""
    cbytes(s::SChunk)

Compressed size of all data in super-chunk
"""
cbytes(s::SChunk) = unsafe_load(s.schunk.cbytes)

"""
    length(s::SChunk)

Number of chunks in super-chunk.
"""
Base.length(s::SChunk) = nchunks(s)

function checkchunksbounds(s::SChunk, bounds)
    checkbounds(Bool, 1:length(s), bounds) || throw(BoundsError(s, bounds))
    nothing
end

"""
    SChunk([storage::Storage])

Create a new super-chunk.

In case that storage.urlpath is not `nothing`, the data is stored on-disk.
if the data file(s) exists the file will be overwritten
"""
function SChunk(storage::Storage = Storage())
    return with_preserved_blosc2_storage(storage) do bstorage
        SChunk(Lib.blosc2_schunk_new(Ref(bstorage)))
    end
end

"""
    Base.copy(schunk::SChunk, [storage::Storage])

Create a copy of super-chunk with new storage settings
"""
function Base.copy(schunk::SChunk, storage::Storage = Storage())
    return with_preserved_blosc2_storage(storage) do bstorage
        SChunk(Lib.blosc2_schunk_copy(schunk.ptr, Ref(bstorage)))
    end
end


text:
# read_samples

"""

Read sample output files created by StanSample.jl and return in the requested `output_format`.
The default output_format is :table. Optionally the list of parameter symbols can be returned.

$(SIGNATURES)

# Extended help

### Required arguments
```julia
* `model`                         : SampleModel
* `output_format=:table`          : Requested format for samples
```

### Optional arguments
```julia
* `include_internals=false`       : Include internal Stan paramenters
* `return_parameters=false`       : Return a tuple of (output_format, parameter_symbols)
* `chains=1:m.num_chains*m.num_cpp_chains` : Chains to be included in output (forked processes)
* `start=1`                       : First sample to be included
* `kwargs...`                     : Capture all other keyword arguments
```

Currently supported output_formats are:

1. :table (DEFAULT: StanTable Tables object, all chains appended)
2. :array (3d array format - [samples, parameters, chains])
3. :namedtuple (NamedTuple object, all chains appended)
4. :namedtuples (Vector{NamedTuple} object, individual chains)
5. :tables (Vector{Tables} object, individual chains)
6. :dataframe (DataFrames.DataFrame object, all chains appended)
7. :dataframes (Vector{DataFrames.DataFrame} object, individual chains)
8. :keyedarray (KeyedArray object from AxisDict.jl)
9. :particles (Dict{MonteCarloMeasurements.Particles})
10. :dimarray (Appended chains DimensionalData.DimArray object)
11. :dimarrays (DimensionalData.DimArray object)
12. :mcmcchains (MCMCChains.Chains object)

Basically chains can be returned as an Array, a KeyedArray, a DimArray, a NamedTuple,
a StanTable, a DataFrame, a Particles or an MCMCChains.Chains object.

Options 8 to 12 are enabled by the presence of AxisKeys.jl, MonteCarloMeasurements.jl,
DimensionalData.jl or MCMCChains.jl.

For NamedTuple, StanTable, DimArray and DataFrame all chains are appended or can be returned
as a Vector{...} for each chain.

With the optional keyword argument `chains` a subset of chains can be included,
e.g. `chains = [2, 4]`.

The optional keyword argument `start` specifies which initial (warm-up) samples
should be removed.

Notes:
1. Use of the Stan `thinning` option will interfere with the value of start.
2. Start is the first sample included, e.g. with 1000 warm-up samples, start should be set to 1001.

The NamedTuple output-format will extract and combine parameter vectors, e.g.
if Stan's cmdstan returns `a.1, a.2, a.3` the NamedTuple will just contain `a`.

For KeyedArray and StanTable objects you can use the overloaded `matrix()` method to
extract a block of parametes:
```
stantable = read_samples(m10.4s, :table)
atable = matrix(stantable, "a")
```

For an appended DataFrame you can use e.g. `DataFrame(df, :log_lik)` to block a
set of variables, in this example the `log_lik.1, log_lik.2, etc.`.

Currently :table is the default chain output_format (a StanTable object).

In general it is safer to specify the desired output_format as this area
is still under heavy development in the Julia eco system. The default
has changed frequently!

"""
function read_samples(model::SampleModel, output_format=:table;
  include_internals=false,
  return_parameters=false,
  chains=1:model.num_chains,
  start=1,
  kwargs...)

  #println(chains)

  (res, names) = read_csv_files(model::SampleModel, output_format;
    include_internals, start, chains, kwargs...
  )

  if return_parameters
    return( (res, names) )
  else
    return(res)
  end

end

max_stars_repo_name: StanJulia/StanSample

max_issues_repo_name: StanJulia/StanSample

max_forks_repo_name: StanJulia/StanSample


text:
# This file is part of Kpax3. License is MIT.

function logmargliksplit!(cl::Vector{Int},
                          k::Int,
                          hi::Int,
                          priorC::AminoAcidPriorCol,
                          support::MCMCSupport)
  support.k = 0
  for l in 1:k
    if cl[l] != hi
      # this cluster is not affected by the split
      support.k += 1
      support.cl[support.k] = cl[l]
    end
  end
  support.k += 2

  support.logmlikcandidate = 0.0
  g = 0
  tmp1 = zeros(Float64, 3)
  tmp2 = zeros(Float64, 2)
  for b in 1:size(support.lp, 3)
    tmp1[1] = priorC.logγ[1]
    tmp1[2] = priorC.logγ[2]
    tmp1[3] = priorC.logγ[3]

    for l in 1:(support.k - 2)
      g = support.cl[l]

      tmp1[1] += support.lp[1, g, b]
      tmp1[2] += support.lp[2, g, b]

      tmp2[1] = priorC.logω[support.k][1] + support.lp[3, g, b]
      tmp2[2] = priorC.logω[support.k][2] + support.lp[4, g, b]

      tmp1[3] += (tmp2[1] > tmp2[2]) ?
                 (tmp2[1] + log1p(exp(tmp2[2] - tmp2[1]))) :
                 (tmp2[2] + log1p(exp(tmp2[1] - tmp2[2])))
    end

    # new cluster i
    tmp1[1] += support.lpi[1, b]
    tmp1[2] += support.lpi[2, b]

    tmp2[1] = priorC.logω[support.k][1] + support.lpi[3, b]
    tmp2[2] = priorC.logω[support.k][2] + support.lpi[4, b]

    tmp1[3] += (tmp2[1] > tmp2[2]) ?
               (tmp2[1] + log1p(exp(tmp2[2] - tmp2[1]))) :
               (tmp2[2] + log1p(exp(tmp2[1] - tmp2[2])))

    # new cluster j
    tmp1[1] += support.lpj[1, b]
    tmp1[2] += support.lpj[2, b]

    tmp2[1] = priorC.logω[support.k][1] + support.lpj[3, b]
    tmp2[2] = priorC.logω[support.k][2] + support.lpj[4, b]

    tmp1[3] += (tmp2[1] > tmp2[2]) ?
               (tmp2[1] + log1p(exp(tmp2[2] - tmp2[1]))) :
               (tmp2[2] + log1p(exp(tmp2[1] - tmp2[2])))

    support.logmlikcandidate += if (tmp1[1] >= tmp1[2]) &&
                                   (tmp1[1] >= tmp1[3])
                                  tmp1[1] + log1p(exp(tmp1[2] - tmp1[1]) +
                                                  exp(tmp1[3] - tmp1[1]))
                                elseif (tmp1[2] >= tmp1[1]) &&
                                       (tmp1[2] >= tmp1[3])
                                  tmp1[2] + log1p(exp(tmp1[1] - tmp1[2]) +
                                                  exp(tmp1[3] - tmp1[2]))
                                else
                                  tmp1[3] + log1p(exp(tmp1[1] - tmp1[3]) +
                                                  exp(tmp1[2] - tmp1[3]))
                                end
  end

  nothing
end