text:
@testset "Network" begin
    @testset "Peer" begin
        @testset "IP -> Bytes" begin
            want = hex2bytes("00000000000000000000ffff0a000001")
            @test Bitcoin.ip2bytes(ip"10.0.0.1") == want
        end
        @testset "Serialize" begin
            peer = Bitcoin.Peer(1, ip"10.0.0.1", UInt16(8333))
            want = hex2bytes("010000000000000000000000000000000000ffff0a000001208d")
            @test Bitcoin.serialize(peer, true) == want
        end
    end
    @testset "Envelope" begin
        @testset "Parse" begin
            msg = hex2bytes("f9beb4d976657261636b000000000000000000005df6e0e2")
            envelope = Bitcoin.io2envelopes(msg)[1]
            @test envelope.command == "verack"
            @test envelope.payload == UInt8[]
            msg = hex2bytes("f9beb4d976657273696f6e0000000000650000005f1a69d2721101000100000000000000bc8f5e5400000000010000000000000000000000000000000000ffffc61b6409208d010000000000000000000000000000000000ffffcb0071c0208d128035cbc97953f80f2f5361746f7368693a302e392e332fcf05050001")
            envelope = Bitcoin.io2envelopes(msg)[1]
            @test envelope.command == "version"
            @test envelope.payload == msg[25:end]
        end
        @testset "Serialize" begin
            msg = hex2bytes("f9beb4d976657261636b000000000000000000005df6e0e2")
            envelope = Bitcoin.io2envelopes(msg)[1]
            @test Bitcoin.serialize(envelope) == msg
            msg = hex2bytes("f9beb4d976657273696f6e0000000000650000005f1a69d2721101000100000000000000bc8f5e5400000000010000000000000000000000000000000000ffffc61b6409208d010000000000000000000000000000000000ffffcb0071c0208d128035cbc97953f80f2f5361746f7368693a302e392e332fcf05050001")
            envelope = Bitcoin.io2envelopes(msg)[1]
            @test Bitcoin.serialize(envelope) == msg
            want = hex2bytes("f9beb4d976657261636b000000000000000000005df6e0e2")
            envelope = Bitcoin.NetworkEnvelope("verack", UInt8[])
            @test Bitcoin.serialize(envelope) == want

        end
    end
    @testset "Message" begin
        @testset "Serialize" begin
            @testset "Version" begin
                v = VersionMessage(zero(UInt64),zero(UInt64))
                @test bytes2hex(Bitcoin.serialize(v)) == "7f11010000000000000000000000000000000000000000000000000000000000000000000000ffff00000000208d000000000000000000000000000000000000ffff00000000208d0000000000000000102f626974636f696e2e6a6c3a302e312f0000000001"
            end
            @testset "GetHeaders" begin
                block_hex = "0000000000000000001237f46acddf58578a37e213d2a6edc4884a2fcad05ba3"
                gh = GetHeadersMessage(hex2bytes(block_hex))
                @test bytes2hex(Bitcoin.serialize(gh)) == "7f11010001a35bd0ca2f4a88c4eda6d213e2378a5758dfcd6af437120000000000000000000000000000000000000000000000000000000000000000000000000000000000"
            end
            @testset "GetData" begin
                hex_msg = "020300000030eb2540c41025690160a1014c577061596e32e426b712c7ca00000000000000030000001049847939585b0652fba793661c361223446b6fc41089b8be00000000000000"
                get_data = GetDataMessage()
                block1 = hex2bytes("00000000000000cac712b726e4326e596170574c01a16001692510c44025eb30")
                append!(get_data, Bitcoin.FILTERED_BLOCK_DATA_TYPE, block1)
                block2 = hex2bytes("00000000000000beb88910c46f6b442312361c6693a7fb52065b583979844910")
                append!(get_data, Bitcoin.FILTERED_BLOCK_DATA_TYPE, block2)
                @test bytes2hex(Bitcoin.serialize(get_data)) == hex_msg
            end
        end
        @testset "Parse" begin
            @testset "Headers" begin
                hex_msg = "0200000020df3b053dc46f162a9b00c7f0d5124e2676d47bbe7c5d0793a500000000000000ef445fef2ed495c275892206ca533e7411907971013ab83e3b47bd0d692d14d4dc7c835b67d8001ac157e670000000002030eb2540c41025690160a1014c577061596e32e426b712c7ca00000000000000768b89f07044e6130ead292a3f51951adbd2202df447d98789339937fd006bd44880835b67d8001ade09204600"
                s = hex2bytes(hex_msg)
                headers = Bitcoin.PARSE_PAYLOAD["headers"](s)
                @test length(headers.headers) == 2
            end
        end
        @testset "MerkleBlock" begin
            hex_msg = "0100000082bb869cf3a793432a66e826e05a6fc37469f8efb7421dc880670100000000007f16c5962e8bd963659c793ce370d95f093bc7e367117b3c30c1f8fdd0d9728776381b4d4c86041b554b852907000000043612262624047ee87660be1a707519a443b1c1ce3d248cbfc6c15870f6c5daa2019f5b01d4195ecbc9398fbf3c3b1fa9bb3183301d7a1fb3bd174fcfa40a2b6541ed70551dd7e841883ab8f0b16bf04176b7d1480e4f0af9f3d4c3595768d06820d2a7bc994987302e5b1ac80fc425fe25f8b63169ea78e68fbaaefa59379bbf011d"
            @testset "Parse" begin
                s = hex2bytes(hex_msg)
                msg = Bitcoin.PARSE_PAYLOAD["merkleblock"](s)
                @test length(msg.hashes) == 4
                @test msg.flags == [true, false, true, true, true, false, false, false]
            end
            @testset "Is Valid" begin
                s = hex2bytes(hex_msg)
                msg = Bitcoin.PARSE_PAYLOAD["merkleblock"](s)
                @test Bitcoin.is_valid(msg)
            end
        end
    end
end


text:
"""
    critical_density(λ_laser)

Critical density of a plasma [cm^-3]

# Arguments
 - `λ_laser`: laser light wavelength, usually in nm
"""
function critical_density(λ_laser)
    return ((ElectronMass * π * SpeedOfLightInVacuum^2) /
            (λ_laser^2 * e_statC^2)) |> u"cm^-3"
end

"""

Refractive index of a plasma

# Arguments
 - `n_e`: electron number density, usually in cm^-3
 - `λ_laser`: laser light wavelength, usually in nm
"""
function refractive_index(n_e, λ_laser)
    η = 1 - (n_e / critical_density(λ_laser)) |> NoUnits
    return η
end

ext: jl


text:
@doc raw"""
    ExponentiatedKernel([a=1])

The exponentiated kernel is a Mercer kernel given by:

```math
\kappa(\mathbf{x},\mathbf{y}) = \exp\left(a \mathbf{x}^\intercal \mathbf{y} \right)
\qquad a > 0
```
"""
struct ExponentiatedKernel{T<:AbstractFloat} <: MercerKernel{T}
    alpha::HyperParameter{T}
    ExponentiatedKernel{T}(α::Real) where {T<:AbstractFloat} = new{T}(
        HyperParameter(convert(T,α), interval(OpenBound(zero(T)), nothing))
    )
end
ExponentiatedKernel(α::T1 = 1.0) where {T1<:Real} = ExponentiatedKernel{floattype(T1)}(α)

@inline exponentiatedkernel(z::T, α::T) where {T<:AbstractFloat} = exp(α*z)

@inline basefunction(::ExponentiatedKernel) = ScalarProduct()
@inline kappa(κ::ExponentiatedKernel{T}, z::T) where {T} = exponentiatedkernel(z, getvalue(κ.alpha))


text:
function getConcurrencyLevel(conf, counts)
	allPredictions  = useLockModel(conf.lock_conf, counts, conf.workloadName);

	if allPredictions.TimeSpentWaiting[1] < 1e-5
		cl = 0; # as if there's not that many people in the system!
	else
		cl = allPredictions.M_total;
	end
	cl
end

max_stars_repo_name: barzan/dbseer

max_issues_repo_name: barzan/dbseer

max_forks_repo_name: barzan/dbseer


text:
using Test
using LazIO

@testset "Conversion" begin
    # test with maximum fields where applicable
    lazp = LazIO.LazPoint(
        X = Int32(5000),
        Y = Int32(6000),
        Z = Int32(7000),
        intensity = UInt16(65535),

        # Evil types
        # U8 return_number : 3;
        # U8 number_of_returns : 3;
        # U8 scan_direction_flag : 1;
        # U8 edge_of_flight_line : 1;  8
        # U8 classification : 5;
        # U8 synthetic_flag : 1;
        # U8 keypoint_flag  : 1;
        # U8 withheld_flag : 1;      8
        return_number = UInt8(255),
        # number_of_returns = UInt8(0),
        # scan_direction_flag = UInt8(0),
        # edge_of_flight_line = UInt8(0),
        classification = UInt8(255),
        # synthetic_flag = UInt8(0),
        # keypoint_flag = UInt8(0),
        # withheld_flag = UInt8(0),
        scan_angle_rank = Int8(127),
        user_data = UInt8(255),
        point_source_ID = UInt16(65535),

        # Another evil type
        # I16 extended_scan_angle;
        # U8 extended_point_type : 2;
        # U8 extended_scanner_channel : 2;
        # U8 extended_classification_flags : 4;  8
        # U8 extended_classification;
        # U8 extended_return_number : 4;
        # U8 extended_number_of_returns : 4;  8
        extended_scan_angle = Int16(32767),
        extended_point_type = UInt8(255),
        # extended_scanner_channel = UInt8(0),
        # extended_classification_flags = UInt8(0),
        extended_classification = UInt8(255),
        extended_return_number = UInt8(25),
        # extended_number_of_returns = UInt8(0),
        dummy = ntuple(i -> UInt8(255), 7),
        gps_time = Float64(8791823.25),
        rgb = ntuple(i -> UInt16(65535), 4),
        wave_packet = ntuple(i -> UInt8(255), 29),
        num_extra_bytes = Int32(24244),
        extra_bytes = pointer("")
    )

    las0 = convert(LazIO.LasPoint0, lazp)
    las1 = convert(LazIO.LasPoint1, lazp)
    las2 = convert(LazIO.LasPoint2, lazp)
    las3 = convert(LazIO.LasPoint3, lazp)

    # test with all zero fields
    lazp = LazIO.LazPoint()
    las0 = convert(LazIO.LasPoint0, lazp)
    las1 = convert(LazIO.LasPoint1, lazp)
    las2 = convert(LazIO.LasPoint2, lazp)
    las3 = convert(LazIO.LasPoint3, lazp)
end


text:
struct TelegramError{T} <: Exception
    msg::T
end

mutable struct TelegramClient
    token::String
    chat_id::String
    parse_mode::String
    ep::String
end

"""
    TelegramClient(token; chat_id, parse_mode, ep, use_globally = true)

Creates telegram client, which can be used to run telegram commands.

# Arguments
- `token::String`: telegram token, which can be obtained in telegram from @BotFather.
- `chat_id`: if set, used as default `chat_id` argument for all chat related commands. To get specific `chat_id`, send message to the bot in telegram application and use [`getUpdates`](@ref) command.
- `parse_mode`: if set, used as default for text messaging commands.
- `use_globally::Bool`: default `true`. If set to `true` then it current client can be used as default in all telegram commands.
- `ep`: endpoint for telegram api. By default `https://api.telegram.org/bot` is used, but you may change it, if you are using proxy or for testing purposes.
"""
function TelegramClient(token = get(ENV, "TELEGRAM_BOT_TOKEN", ""); chat_id = get(ENV, "TELEGRAM_BOT_CHAT_ID" ,""), parse_mode = "", ep = "https://api.telegram.org/bot", use_globally = true)
    client = TelegramClient(token, chat_id, parse_mode, ep)
    if use_globally
        useglobally!(client)
    end

    return client
end

mutable struct TelegramOpts
    client::TelegramClient
    upload_kw::Vector{Symbol}
    timeout::Int
end

const DEFAULT_OPTS = TelegramOpts(TelegramClient(use_globally = false), [:photo, :audio, :thumb, :document, :video, :animation, :voice, :video_note, :sticker, :png_sticker, :tgs_sticker, :certificate, :files], 30)

"""
    useglobally!(tg::TelegramClient)

Set `tg` as default client in all `Telegram.API` functions.

# Example
```julia
tg = TelegramClient(ENV["TG_TOKEN"])
useglobally!(tg)

getMe() # new client is used in this command by default.
```
"""
function useglobally!(client::TelegramClient)
    DEFAULT_OPTS.client = client

    return client
end

function token(client::TelegramClient)
    isempty(client.token) || return client.token
    token = get(ENV, "TELEGRAM_BOT_TOKEN", "")
    isempty(token) && @warn "Telegram token is empty"
    return token
end

function chatid(client::TelegramClient)
    isempty(client.chat_id) || return client.chat_id
    chat_id = get(ENV, "TELEGRAM_BOT_CHAT_ID", "")
    return chat_id
end

ioify(elem::IO) = elem
ioify(elem) = IOBuffer(elem)
function ioify(elem::IOBuffer)
    seekstart(elem)
    return elem
end

pushfiles!(body, keyword, file) = push!(body, keyword => file)
function pushfiles!(body, keyword, file::Pair)
    push!(body, keyword => HTTP.Multipart(file.first, ioify(file.second)))
end

process_params(x::AbstractString) = x
process_params(x) = JSON3.write(x)

function query(client::TelegramClient, method; params = Dict())
    req_uri = client.ep * token(client) * "/" * method
    intersection = intersect(keys(params), DEFAULT_OPTS.upload_kw)
    if isempty(intersection)
        headers = ["Content-Type" => "application/json", "Connection" => "Keep-Alive"]
        json_params = JSON3.write(params)
        res = HTTP.post(req_uri, headers, json_params; readtimeout = DEFAULT_OPTS.timeout + 1)
    else
        body = Pair[]
        for (k, v) in params
            k in intersection && continue
            push!(body, k => process_params(v))
        end
        for k in intersection
            pushfiles!(body, k, params[k])
        end
        res = HTTP.post(req_uri, ["Connection" => "Keep-Alive"], HTTP.Form(body); readtimeout = DEFAULT_OPTS.timeout + 1)
    end

    response = JSON3.read(res.body)

    if response.ok
        return response.result
    else
        throw(TelegramError(response))
    end
end

"""
    apiquery(method, client; kwargs...)

Sends `method` request to the telegram. It is recommended to use only if some of Telegram API function is not wrapped already.

# Arguments
- `method::String`: method name from the Telegram API, for example "getMe"
- `client::TelegramClient`: telegram client, can be omitted, in this case default client is used.
"""
function apiquery(method, client::TelegramClient = DEFAULT_OPTS.client; kwargs...)
    params = Dict{Symbol, Any}(kwargs)
    params[:chat_id] = get(params, :chat_id, chatid(client))
    params[:parse_mode] = get(params, :parse_mode, client.parse_mode)
    query(client, method, params = params)
end


text:
using Test
using ShuffleProofs: Node, Leaf, Tree, encode, decode

### Testing parser

@test decode("01 00000002 2D52") == Leaf(UInt8[0x2d, 0x52])
@test decode("00 00000002 01 00000001 AF 01 00000002 03E1") == Node([Leaf(UInt8[0xaf]), Leaf(UInt8[0x03, 0xe1])])
@test decode("00 00000002 00 00000002 01 00000001 AF 01 00000002 03E1 01 00000002 2D52") == Node([Node([Leaf(UInt8[0xaf]), Leaf(UInt8[0x03, 0xe1])]), Leaf(UInt8[0x2d, 0x52])])

#@test convert(Tree, "01 00000002 2D52") == Leaf(UInt8[0x2d, 0x52])
#@test convert(Tree, "00 00000002 01 00000001 AF 01 00000002 03E1") == Node([Leaf(UInt8[0xaf]), Leaf(UInt8[0x03, 0xe1])])
#@test convert(Tree, "00 00000002 00 00000002 01 00000001 AF 01 00000002 03E1 01 00000002 2D52") == Node([Node([Leaf(UInt8[0xaf]), Leaf(UInt8[0x03, 0xe1])]), Leaf(UInt8[0x2d, 0x52])])

y = hex2bytes(replace("00 00000002 00 00000002 01 00000001 AF 01 00000002 03E1 01 00000002 2D52", " "=>""))
tree = decode(y)
@test convert(Tuple{Tuple{Int32, Int32}, String}, tree) == ((175, 993), "-R")

@test Tree(((UInt8(175), UInt16(993)), "-R")) == tree

### Real life example for prime group specification:

group_spec = "00000000020100000020636f6d2e766572696669636174756d2e61726974686d2e4d6f645047726f757000000000040100000041009a91c3b704e382e0c772fa7cf0e5d6363edc53d156e841555702c5b6f906574204bf49a551b695bed292e0218337c0861ee649d2fe4039174514fe2c23c10f6701000000404d48e1db8271c17063b97d3e7872eb1b1f6e29e8ab7420aaab8162db7c832ba1025fa4d2a8db4adf69497010c19be0430f7324e97f201c8ba28a7f1611e087b3010000004100300763b0150525252e4989f51e33c4e6462091152ef2291e45699374a3aa8acea714ff30260338bddbb48fc7446b273aaada90e3ee8326f388b582ea8a073502010000000400000001"

@test encode(decode(group_spec)) == hex2bytes(group_spec);

tree = decode(group_spec)

(group_name, (p, q, g, e)) = convert(Tuple{String, Tuple{BigInt, BigInt, BigInt, UInt32}}, tree)

@test p == 2*q + 1
@test powermod(g, q + 1, p) == g


text:
function nlsfitworker(model::Function, y::Matrix{T}, x::Vector{T}, p0::Vector{T}; kwargs...) where T
  # performs Levenberg-Marquardt fitting
  nparams = length(p0)
  n = size(y,2)
  params = zeros(nparams, n)
  resids = similar(y)
  for j in 1:n
    fit = curve_fit(model, x, y[:,j], p0)
    params[:,j] = fit.param
    resids[:,j] = fit.resid
  end
  (params, resids)
end

function nlsfit(f::Function, y::Matrix{T}, idxs::Vector{Int}, x::Vector{T},
  p0::Vector{T}; kwargs...) where T
  timerStart = time_ns()
  nt, n = size(y)
  nw = nworkers()
  nidxs = length(idxs)
  params = zeros(length(p0), n)
  resids = zeros(nt, n)
  dof = nt - length(p0)
  if nw > 1
    reflist = Any[]
    idxperworker = ceil(Int, nidxs / nw)
    workeridxs = Any[idxs[(w-1)*idxperworker+1:min(w*idxperworker,nidxs)] for w in 1:nw]
    ni = map(length, workeridxs)
    @dprint "work distribution: $nt x $ni points"
    workerids = workers()
    for w in 1:nw
      r = @spawnat workerids[w] nlsfitworker(f, y[:,workeridxs[w]], x, p0; kwargs...)
      push!(reflist, r)
    end
    for w in 1:nw
      params[:,workeridxs[w]], resids[:,workeridxs[w]] = fetch(reflist[w])
    end
  else   # run on one CPU core
    @dprint "running $nt x $nidxs points on one CPU core"
    params[:,idxs], resids[:,idxs] = nlsfitworker(f, y[:,idxs], x, p0; kwargs...)
  end
  t = (time_ns()-timerStart)/1e9
  vps = nidxs / t
  @dprint @sprintf "processed %d voxels in %.1f s (%.1f vox/s)\n" nidxs t vps
  (params, resids, dof)
end

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

function cumtrapz(x::Array{Float64, 1}, y::Array{Float64}, dim::Integer=1)
    # Cumulative trapezoidal integration
    # Obtained from:
    # https://github.com/jdtuck/ElasticFDA.jl/blob/master/src/misc_funcs.jl
    perm = [dim:max(length(size(y)),dim); 1:dim-1];
    y = permutedims(y, perm);
    if ndims(y) == 1
        n = 1;
        m = length(y);
    else
        m, n = size(y);
    end

    if n == 1
        dt = diff(x)/2.0;
        z = [0; cumsum(dt.*(y[1:(m-1)] + y[2:m]))];
    else
        dt = repmat(diff(x)/2.0,1,n);
        z = [zeros(1,n); cumsum(dt.*(y[1:(m-1), :] + y[2:m, :]),1)];
        z = ipermutedims(z, perm);
    end

    return z
end


text:
function toDataFrame(s::Matrix, p::Matrix, z::Matrix, x::Matrix = zeros(size(s)) )
    df = DataFrame();
    J = size(s,2);
    for j = 0:J-1
        df[!, "shares$j"] =  s[:,j+1];
        df[!, "prices$j"] =  p[:,j+1];
        df[!, "demand_instruments$j"] =  z[:,j+1];
        df[!, "x$j"] =  x[:,j+1];
    end
    return df;
end


text:
module TestTapir

using Base.Experimental: Tapir
using TapirSchedulers
using TapirSchedulersBenchmarks.BenchFib: WSFib, DFFib, CPFib, RPFib
using Test

@noinline produce(x) = Base.inferencebarrier(x)::typeof(x)

function simple_spawn(tgf)
    Tapir.@output a b
    Tapir.@sync tgf() begin
        Tapir.@spawn a = produce(111)
        b = produce(222)
    end
    return a + b
end

function test_simple_spawn()
    @test simple_spawn(WorkStealingTaskGroup) == 333
    @test simple_spawn(DepthFirstTaskGroup) == 333
end

function simple_spawn_macro()
    Tapir.@output a b
    @sync_ws begin
        Tapir.@spawn a = produce(111)
        b = produce(222)
    end
    return a + b
end

function test_simple_spawn_macro()
    @test simple_spawn_macro() == 333
end

test_fib_work_stealing() = check_fib(WSFib.fib)
test_fib_depth_first() = check_fib(DFFib.fib)
test_fib_constant_priority() = check_fib(CPFib.fib)
test_fib_random_priority() = check_fib(RPFib.fib)

function check_fib(fib)
    @test fib(1) == 1
    @test fib(2) == 1
    @test fib(3) == 2
    @test fib(4) == 3
    @test fib(5) == 5
    @test fib(6) == 8
    @test fib(10) == 55
end

end  # module


text:
module LibBSON

using Compat

include( "../deps/deps.jl")

LIBBSON_VERSION = VersionNumber(
    ccall((:bson_get_major_version, libbson), Cint, ()),
    ccall((:bson_get_minor_version, libbson), Cint, ()),
    ccall((:bson_get_micro_version, libbson), Cint, ()),
    )
export LIBBSON_VERSION

import Base.hash,
    Base.convert,
    Base.show,
    Base.error,
    Base.start,
    Base.next,
    Base.done,
    Base.string,
    Base.length

include("BSONOID.jl")
include("BSONError.jl")
include("BSONObject.jl")
include("BSONArray.jl")
include("BSONType.jl")
include("BSONIter.jl")

end


text:
using Documenter, CoupledODETools

makedocs(;
    modules=[CoupledODETools],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/jamesjscully/CoupledODETools.jl/blob/{commit}{path}#L{line}",
    sitename="CoupledODETools.jl",
    authors="Jack Scully",
    assets=String[],
)

deploydocs(;
    repo="github.com/jamesjscully/CoupledODETools.jl",
)


text:
module MaxHelpingHandCustomPlanets

using ..Ahorn, Maple

@mapdef Effect "MaxHelpingHand/CustomPlanets" CustomPlanets(only::String="*", exclude::String="", count::Number=32, color::String="", scrollx::Number=1.0, scrolly::Number=1.0,
    speedx::Number=0.0, speedy::Number=0.0, directory::String="MaxHelpingHand/customplanets/bigstars", animationDelay::Number=0.1)

placements = CustomPlanets

function Ahorn.canFgBg(effect::CustomPlanets)
    return true, true
end

function Ahorn.editingOptions(effect::CustomPlanets)
    return Dict{String, Any}(
        "directory" => String[
            "MaxHelpingHand/customplanets/bigstars",
            "MaxHelpingHand/customplanets/smallstars",
            "MaxHelpingHand/customplanets/bigplanets",
            "MaxHelpingHand/customplanets/smallplanets",
            "MaxHelpingHand/customplanets/rainbowstars"
        ]
    )
end

end

max_stars_repo_name: Kitty-Cats/MaxHelpingHand

max_issues_repo_name: Kitty-Cats/MaxHelpingHand

max_forks_repo_name: Kitty-Cats/MaxHelpingHand


text:
module HomotheciesTests

using Gridap.TensorValues
using Gridap.Arrays
using Gridap.Fields

origin = Point(1,1)
scaling = Point(2,2)

h = Homothecy(origin,scaling)

x1 = Point(0,0)
x2 = Point(1,1)
x3 = Point(2,1)
x = [x1,x2,x3]
r = Point{2,Int}[(-1, -1), (1, 1), (3, 1)]
∇r = TensorValue{2,2,Int,4}[(2, 0, 0, 2), (2, 0, 0, 2), (2, 0, 0, 2)]
test_field(h,x,r,grad=∇r)

end # module


text:
function mh(
    likelihood::Function,
    prior::Function,
    propdistpdf::Function,
    propdistsample::Function,
    startval::Vector,
    mcit::Int64,
    burn::Int64
    )

    sample = zeros(mcit,length(startval))
    sample[1,:] = startval

    for i in 1:mcit-1
        curr = sample[i,:]

        #create new sample
        proposalval = curr + propdistsample()

        #calculate acceptance probability
        A = (likelihood(proposalval)/likelihood(curr)*prior(proposalval)/
            prior(curr)*propdistpdf(proposalval-curr)/
            propdistpdf(curr-proposalval))[1]

        #accept or reject sample
        if A>=rand(Uniform(0,1))
            sample[i+1,:] = proposalval
        else
            sample[i+1,:] = curr
        end
    end
    return sample[burn+1:mcit,:]
end

function gibbssample(
    condist::Vector,
    startval::Vector,
    mcit::Int64,
    burn::Int64
    )

    sample = zeros(mcit, length(startval))
    sample[1,:] = startval

    for i in 1:mcit-1
        sample[i+1,:] = sample[i,:]
        for j in 1:size(sample, 2)
            #update each variable by itself
            sample[i+1,j] = condist[j](sample[i+1,:])
        end
    end
    return sample[burn+1:mcit,:]
end

function tresample(
    sset::Array,
    hatw::Array
    )
    snum = size(sset,1)
    rsset = zeros(snum, size(sset,2))
    hatwn = zeros(snum, 1)
    number = [1:snum...]

    for i in 1:snum
        #resample according to weight
        resample = StatsBase.sample(number, ProbabilityWeights(vec(hatw)))
        rsset[i,:] = sset[resample, :]
        hatwn[i] = hatw[resample]
    end
    return rsset, hatwn
end

function calcbeta(
    likelihood::Function,
    betaj::Float64,
    sset::Array,
    snum::Int64
    )

    betaj1 = betaj + (1.1 - betaj)*0.5
    likelihoodcov = zeros(snum, 1)
    step = (1.1 - betaj)*0.25
    cov = 0

    for i in 1:snum
        likelihoodcov[i] = likelihood(sset[i,:])
    end
    #beta approximates so that coefficient of variance --> 1
    while cov> 1.0001 || cov < 0.9999
        cov = std(likelihoodcov.^(betaj1-betaj))/mean(likelihoodcov.^(betaj1-betaj))
        if cov < 1
            betaj1 += step
            step = 0.5*step
        else
            betaj1 -= step
            step = 0.5*step
        end
        if betaj1>1
            break
        end
    end
    betaj1 = min(1, betaj1)
    return betaj1, likelihoodcov.^(betaj1-betaj)
end

function calcsigma(
    sset::Array,
    hatw::Array,
    gamma::Real
    )
    sdim = size(sset, 2)
    snum = size(sset, 1)

    thet = zeros(sdim, 1)
    sig = zeros(sdim, sdim)

    #formula for calculating Sigma in Ching and Chen 2007
    for i in 1:snum
        thet = thet + sset[i,:].*hatw[i]
    end

    for i in 1:snum
        sig = sig + ((sset[i,:]-thet)*(sset[i,:]-thet)').*hatw[i]
    end
    sig = sig.*gamma^2

    return sig
end

function tmcmc(
    likelihood::Function,
    prior::Function,
    priorsample::Function,
    snum::Int64,
    gamma::Real
    )

    j = 0
    betaj = 0.0
    sset = priorsample(snum)

    while betaj < 1
        betaj1, w = calcbeta(likelihood, betaj, sset, snum)
        hatw = w./sum(w)#normalize

        sset, w = tresample(sset, hatw)
        hatw = w./sum(w)

        sig = calcsigma(sset, hatw, gamma)

        for i in 1:snum
            #MH step for each sample
            mvpropdist = MvNormal(zeros(size(sset, 2)), sig)
            propdistsample() = rand(mvpropdist)
            propdistpdf(x) = pdf(mvpropdist, x)
            likelihoodbeta(x) = likelihood(x)^(betaj1)
            sset[i,:] = mh(likelihoodbeta, prior, propdistpdf, propdistsample, sset[i,:], 2,1)
        end
        betaj = betaj1
    end

    return sset
end

function resample(
    sset::Array,
    hatw::Array
    )
    snum = size(sset,1)
    rsset = zeros(snum, size(sset,2))
    number = [1:snum...]

    for i in 1:snum
        #weighted resampling
        resample = StatsBase.sample(number, ProbabilityWeights(vec(hatw)))
        rsset[i,:] = sset[resample, :]
    end

    return rsset
end

function smc(
    likelihood::Function,
    prior::Function,
    priorsample::Function,
    propdistpdf::Function,
    propdistsample::Function,
    snum::Int64,
    term::Function
    )

    j = 0
    sset = priorsample(snum)
    w = zeros(snum,1)

    #calculate weights
    for i in 1:snum
        w[i] = likelihood(sset[i,:])*prior(sset[i,:])
    end

    hatw = w./sum(w)#normalize

    while term(sset, j)
        j += 1
        Neff = 1/sum(hatw.^2)#efficient sample size
        if Neff < snum/2 #test for degeneracy
            sset = resample(sset, hatw)
            hatw = ones(snum,1)./snum
        end

        rsset = zeros(snum, size(sset,2))
        w = zeros(snum, 1)

        for i in 1:snum
            #one step in a random direction
            rsset[i,:] = sset[i,:] + propdistsample()
            w[i] = hatw[i].*likelihood(rsset[i,:])/likelihood(sset[i,:])*
                prior(rsset[i,:])/prior(sset[i,:])*
                propdistpdf(rsset[i,:]-sset[i,:])/propdistpdf(sset[i,:]-rsset[i,:])
        end

        sset = rsset
        hatw = w./sum(w)
    end
    return sset
end

function calcequaltails(
    dsample::Array,
    quantile::Float64,
    )
    snum = size(dsample, 1)
    dsample = sort(dsample, dims = 1)

    #calc cutoff points
    low = floor(Int32,snum*(1-quantile)/2)
    high = floor(Int32,snum - snum*(1-quantile)/2)

    rv = zeros(2, size(dsample,2))

    for i in 1:size(dsample,2)
        rv[:,i] = [dsample[low,i]; dsample[high,i]]
    end
    return rv
end

function grconvergence(
    likelihood::Function,
    prior::Function,
    priorsample::Function,
    propdistpdf::Function,
    propdistsample::Function,
    mcit::Int64,
    nchains::Int64
    )

    startval = priorsample(nchains)
    dims = size(startval, 2)
    rhat = 0
    mcitb = floor(Int, mcit/2)
    samples = zeros(floor(Int, mcit/2), size(startval,2), nchains)

    #create n Markov chains
    for i in 1:nchains
        samples[:,:,i] = mh(likelihood, prior, propdistpdf, propdistsample,startval[i,:], mcit, floor(Int, mcit/2))
    end

    if dims == 1
        chainvariance = zeros(nchains)
        chainmean = zeros(nchains)
        overallmean = 0
        for i in 1:nchains
            chainmean[i] = mean(samples[:,1,i], dims=1)[1]
            overallmean += chainmean[i]/nchains
            for j in 1:floor(Int,mcit/2) chainvariance[i] += (samples[j,1,i] - chainmean[i])^2 end
        end
        #vector of variances for each chain
        chainvariance = chainvariance./(mcitb-1)

        #variance of all chains
        mvariance = 0
        for i in 1:nchains mvariance += chainvariance[i]/nchains end

        #variance bias estimator
        Bn = 0
        for i in 1:nchains Bn += (chainmean[i]-overallmean)^2 end
        Bn = Bn/(nchains-1)

        #over corrected variance
        sig2 = ((mcitb-1)/mcitb)*mvariance + Bn

        rhat = sqrt(sig2/mvariance)
    else
        chainvariance = zeros(dims, dims, nchains)
        chainmean = zeros(dims, nchains)
        overallmean = zeros(dims)

        for i in 1:nchains
            chainmean[:, i] = mean(samples[:,:,i], dims=1)
            overallmean += chainmean[:, i]./nchains
            for j in 1:floor(Int, mcit/2)
                v = samples[j,:,i]-chainmean[:,i]
                chainvariance[:,:,i] += v*v'
            end
        end
        #array of covariance matrices for each chain
        chainvariance = chainvariance./(mcitb-1)

        #covariance matrix of all chains
        mvariance = zeros(dims, dims)
        for i in 1:nchains mvariance += chainvariance[:,:,i]./nchains end

        #bias estimator
        Bn = zeros(dims, dims)
        for i in 1:nchains
            v = chainmean[:,i] - overallmean
            Bn += v*v'
        end

        Bn = Bn*mcitb/(nchains-1)
        mvariance = inv(mvariance)
        maxeig = maximum(eigvals(mvariance*Bn))

        rhat = sqrt((mcitb-1)/mcitb + maxeig/mcitb)
    end

    return rhat
end


text:
using UUIDs, TOML

abstract type AbstractCell end
mutable struct InputCell <: AbstractCell
    ID::UUID
    Text::String
    env::Dict
end

mutable struct OutputCell <: AbstractCell
    ID::UUID
    output::String
end

mutable struct Cell <: AbstractCell
    input::Cell
    output::Cell
    hidden::bool
    lang::Symbol
    env::UUID
end


text:
simplifyProject(P::Project) = map(P) do branch
	map(branch) do solution
		solution.data
	end
end

#TODO actually convert(T, Project)
function complicateProject(V) #::Vector{Vector{Vector{Float64}}}
	P = Project()
	branches = map(V) do bData
		branch = Branch(P)
		solutions::Vector{Solution} = map(bData) do sData
			Solution(sData, branch)
		end
		branch.solutions = solutions
		return branch
	end
	P.branches = branches
	return P
end

"""
    create(homotopy, jacobian, projection)
"""
create(v...) = begin
	ses = Session()
	ses.P = Project()
	ses.core = Galerkin(ses, v...)
	ses.cont = PC(ses)
	ses.viz = GalerkinViz(ses)
	show(ses.cont); show(ses.core); show(ses.viz)
	return ses
end

"""
    save(filename, session[, overwrite])
"""
save(filename, S::Session; overwrite=false) = begin
	if !isfile("save/$(filename)") || overwrite; open("save/$(filename)", "w") do f serialize(f, simplifyProject(S.P)) end
	else error("File already exists. Use  overwrite=true  .")end
	return Void
end

#TODO restore non-serializable stuff (figures, observer)
"""
    load(filename, homotopy, jacobian, projection)
"""
load(filename, v...) = begin
	V = open(deserialize, "save/$(filename)")
	ses = Session()
	ses.P = complicateProject(V)
	ses.core = Galerkin(ses, v...)
	ses.cont = PC(ses)
	ses.viz = GalerkinViz(ses)
	show(ses.cont); show(ses.core); show(ses.viz)
	return ses
end


text:
"""
    epot_monthly(datevec)

Use montly average values for potential evapotranspiration.
See report "The nordic HBV model" by Nils Roar Saelthun.

"""

function epot_monthly(datevec)

    # Monthly values of potential evapotranspiration (mm/day)

    epot_month = [0.7, 0.7, 0.7, 1.0, 1.3, 1.4, 1.3, 1.1, 1.0, 0.9, 0.7, 0.7]

    # Assign montly values to days

    epot = zeros(Float64, length(datevec))

    for i in eachindex(datevec)

        imonth = Dates.month(datevec[i])

        epot[i] = epot_month[imonth]

    end

    return(epot)

end

"""
    epot_zero(datevec)

Set potential evapotranspiration to zero.

"""

function epot_zero(datevec)

    epot = zeros(Float64, length(datevec))

    return(epot)

end

max_stars_repo_name: jmgnve/Vann

max_forks_repo_name: jmgnve/Vann


text:
using Test
include("plots_1.jl")
include("diffeq_1.jl")
include("diffeq_plots_1.jl")
include("diffeq_plots_2.jl")

#=
using SnoopCompileCore
list = []
@testset "Plots" begin
    push!(list, @snoopr include("plots.jl"))
end
@testset "DifferentialEquations" begin
    push!(list, @snoopr include("differential_equations.jl"))
end
using SnoopCompile

invalidations = list[1]
length(uinvalidated(invalidations))
trees = invalidation_trees(invalidations)
length(trees)
ftrees = filtermod(MechGluecode, trees)
ftrees = filtermod(MechGlueDiffEq, trees)

[(length(trees[i].backedges), i) for i in 1:37] |> sort!

method_invalidations = trees[34]
length(method_invalidations.backedges)
root = method_invalidations.backedges[1]
show(root; minchildren=0)

=#


text:
# This test file test that deserialization of method responses from Slack
# execute wihout errors.
#
# The tests are not actually present in this git repository. These tests
# come from the "Tester" tab in the Slack API pages for each method. The issue is that they
# come from actual requests done by my user, and so they contain a bunch of data which I
# might not want public, I believe including the test token provided by Slack. Therefore
# each test is located on my machine only, and is ignored by the .gitignore file.
#
# Note that these tests don't check that the deserialized fields have the correct values.
# Primarily these tests confirm those fields that are mandatory, and those that are optional.
# When creating the first of these test methods I came across many fields which are optional,
# because they're not present in the JSON response, but were marked as mandatory in our code.
# In the distant future I will add that functionality to these tests though. Each test file
# will have a companion Julia file which contains the expected deserialized object. It's
# time consuming to write tests though.

facts("Deserializing Slack Tester examples") do
    testcase_dir = "test/slacktester"
    rel_testcase_dir = "slacktester"
    istestcase = x -> isfile(joinpath(testcase_dir, string(x))) && endswith(x, ".json")
    testcase_dir_entries = []

    try
        testcase_dir_entries = readdir(testcase_dir)
    catch ex
        isa(ex, SystemError) || rethrow(ex)
        println("No test/slacktester directory found, so not testing Slack examples.")
        return
    end

    testcases = filter(istestcase, testcase_dir_entries)

    for testcase_file in testcases
        context("$(testcase_file)") do
            testcase_with_path = joinpath(testcase_dir, string(testcase_file))
            jl_file, _json_ext = splitext(testcase_file)
            jl_file_with_path = joinpath(rel_testcase_dir, jl_file * ".jl")
            datatype = eval(evalfile(jl_file_with_path))
            tc_json = JSON.parsefile(testcase_with_path)
            DandelionSlack.deserialize(datatype, tc_json)
        end
    end
end


text:
# Set of functions
#
# After matrices generation, Betti curves will be generated
# For better optimization, operations on indices should be done first

# TODO Check if the threshold of values is applied here and if it has some
# 		consequence on results
using Eirene
using Random

include("MatrixToolbox.jl")

struct PlottingData
	mat_size::Int64
	dim::Int
	src_pts_number::Int
	trgt_pts_number::Int
	src_points::Vector{Int}
	trgt_points::Array{Int}
	targets::Vector{Int}

	# Constructor for input data
	function PlottingData(mat_size::Int, dim::Int, src_pts_number::Int,
							trgt_pts_number::Int)
		steps_set = [0]
		src_points = [0]
		trgt_points = [0]

		new(mat_size::Int, dim::Int, src_pts_number::Int,
								trgt_pts_number::Int, src_points,
								trgt_points, steps_set)
	end

	function PlottingData(mat_size::Int, dim::Int, src_pts_number::Int,
							trgt_pts_number::Int, src_points::Vector{Int},
							trgt_points::Array{Int}, trgt_sptep::Int)
		if trgt_sptep == 0
			steps_set = [trgt_pts_number]
		else
	        steps_set = collect(1:trgt_sptep:trgt_pts_number)
	        # Make sure that all points are used
	        isempty(findall(x->x==trgt_pts_number, steps_set)) && push!(steps_set, trgt_pts_number)
		end

		new(mat_size::Int, dim::Int, src_pts_number::Int,
								trgt_pts_number::Int, src_points::Vector{Int},
								trgt_points::Array{Int}, steps_set)
	end
end

function get_replacing_points(mat_size, src_pts_number, trgt_pts_number)
	total_pts_number = src_pts_number + src_pts_number*trgt_pts_number
	total_pts_number > mat_size && error("Too many points to substitute!")

	elements_collection = randcycle(mat_size)
    src_points = elements_collection[1:src_pts_number]
    strat_val = src_pts_number+1
    stop_val = total_pts_number
    trgt_points = elements_collection[strat_val:stop_val]
    trgt_points = reshape(trgt_points, trgt_pts_number, src_pts_number)

    return src_points, trgt_points
end

# ===

function replace_matrix_rows(matrix, srcs, trgts)
    replacement_row = get_row(matrix, pt_src)

    new_matrix = set_row(matrix, trgt_points, replacement_row)
end

function get_row(matrix, pt_src)
    return  matrix[pt_src,:]
end

function set_row!(matrix::Array, pt_trgt::Int, replacement_row::Array)
    @debug "set_row! func"
    replacement_row[pt_trgt] = 0

    matrix[pt_trgt, :] .= replacement_row
    matrix[:, pt_trgt] .= replacement_row

    return matrix
end

function set_row(matrix::Array, pt_trgt::Int, replacement_row::Array)
    @debug "set_row func"
    new_matrix = copy(matrix)
    return set_row!(new_matrix, pt_trgt, replacement_row)
end

function set_row(matrix::Array, pt_trgt::Array, replacement_row::Array)
    @debug "set_row func"
    new_matrix = copy(matrix)
    for point in pt_trgt
        new_matrix = set_row!(new_matrix, point, replacement_row)
    end
    return new_matrix
end

function matrix_organization(matr, src_points, trgt_points)
	working_matrix = copy(matr)
	mat_size = size(matr, 1)
	src_number = size(src_points,1)

	swapping_sources = Any[]
	step = 0
	matrix_indices = CartesianIndices((1:mat_size, 1:mat_size))
	matrix_indices = findall(x->x[1]<x[2], matrix_indices)
	sorted_values = matr[matrix_indices]

	ordered_indices = sort!([1:mat_size;],
                        by=i->(sorted_values[i],matrix_indices[i]))

	sort(matr[:,1])
	# matr[src_points[src_pt],1]
	for src_pt = 1:src_number
		# find all source
		target_set = findall(x-> x==matr[src_points[src_pt],1], matr[:,1])

		# swap all equivalents
		for tragt = target_set
			swap_rows!(working_matrix, tragt, mat_size-step)
			step +=1
		end
	end
	return working_matrix
end

# matrix = ordered_matrices_collection[15]
function swap_rows!(matrix, src_row_num, trgt_row_num)
	src_backup = copy(matrix[src_row_num,:])
	trgt_backup = copy(matrix[trgt_row_num,:])

	matrix[src_row_num, trgt_row_num] = matrix[trgt_row_num, trgt_row_num]
	matrix[trgt_row_num, src_row_num] = matrix[src_row_num,src_row_num]

	matrix[src_row_num,src_row_num] = trgt_backup[src_row_num]
	matrix[trgt_row_num, trgt_row_num] = src_backup[trgt_row_num]

	src_backup = copy(matrix[src_row_num,:])
	matrix[src_row_num,:] .= matrix[trgt_row_num, :]
	matrix[trgt_row_num, :] .= src_backup

	matrix[:, src_row_num] = matrix[src_row_num,:]
	matrix[:, trgt_row_num] = matrix[trgt_row_num,:]
end

function swap_rows(matrix, src_row_num, trgt_row_num)
	new_matrix = copy(matrix)
	swap_rows!(new_matrix, src_row_num, trgt_row_num)
	return new_matrix
end

function ordering_matrix_analysis(test_data::PlottingData;generation_function=get_geom_matrix)
	mat_size = test_data.mat_size
	dim = test_data.dim
	src_pts_number = test_data.src_pts_number
	trgt_pts_number = test_data.trgt_pts_number
	trgt_steps = 0

	src_points, trgt_points = get_replacing_points(mat_size, src_pts_number, trgt_pts_number)
	distance_matrix = generation_function(mat_size, dim)

	distance_matrices_collection = get_dist_mat_collection(distance_matrix, src_points, trgt_points, trgt_steps)

	ordered_matrices_collection = get_ordered_set(distance_matrices_collection)
	bettis_collection = get_bettis_collection(ordered_matrices_collection)

	plot_data = PlottingData(mat_size, dim, src_pts_number, trgt_pts_number, src_points, trgt_points, trgt_steps)

	plotting_data = print_hmap_with_bettis(ordered_matrices_collection,
												bettis_collection, plot_data)

	return distance_matrices_collection, ordered_matrices_collection, bettis_collection, plot_data
end

# =================================
# Matrix modification functions

function make_matrix_steps!(input_matrix, step_number; step_size=2 )
	# input_matrix = ord_mat
	# step_number = 13
	rows = step_number:(step_number+step_size-1)
    cols = 1:step_number-1
	min_value = findmin(input_matrix[rows,cols])[1]
	input_matrix[rows,cols] .= min_value
	input_matrix[cols,rows] .= min_value
end

"""
function add_step_to_matrix(input_matrix, last_components)

Takes a symmetric matrix 'input_matrix' and appends 2 columns and 2 rows such that
resulting geometric object structure is bigger by 1 dimension. 'last_components'
determines which etries in the matrix are used for closing high dimensional simplices.
"""
function add_step_to_matrix(input_matrix, last_components)
	matrix_size = size(input_matrix,1)
	new_matrix = zeros(Int,matrix_size +2, matrix_size+2)
	new_matrix[1:matrix_size,1:matrix_size] .= input_matrix
	min_closing_component = findmin(input_matrix[last_components])[1]
	new_row1 = range(min_closing_component, length=matrix_size)
	new_row2 = range(findmax(new_row1)[1]+1, length=matrix_size)

	last = 2
	new_matrix[matrix_size+1,1:end-last] = new_row1
	new_matrix[matrix_size+2,1:end-last] = new_row2
	new_matrix[1:end-last,matrix_size+1] = new_row1
	new_matrix[1:end-last,matrix_size+2] = new_row2

	# Adjust last components
	max_new_matrix = findmax(new_matrix)[1]
	new_matrix[last_components].=input_matrix[last_components].+(max_new_matrix-min_closing_component+1)
	new_matrix[end-1,end  ] = findmax(new_matrix)[1]+1
	new_matrix[end,  end-1] = findmax(new_matrix)[1]

	new_max_val = findmax(new_matrix)[2]
	new_component = copy(last_components)
	push!(new_component, new_max_val)
	push!(new_component, CartesianIndex(new_max_val[2], new_max_val[1]))

	return new_matrix, new_component
end


text:
module TestExamplesTests
using LiterateTest.AssertAsTest: @assert
include("../examples/tests.jl")
end  # module


text:
import Base.empty!,
       Base.size,
       Base.convert

type Features
    x::Array{Float64,1}
    i::Int
    function Features(l::Int)
        new(zeros(Float64, l), 1)
    end
end

macro append!(v,w)
    quote
        $(esc(v)).x[$(esc(v)).i] = float($(esc(w)))
        $(esc(v)).i += 1
    end
end

function empty!(v::Features)
    v.i = 1 ; v.x[:] = 0
end

function size(v::Features)
    size(v.x)
end

function convert(::Type{Array{Float64,1}}, v::Features)
    v.x
end


text:
"""
    OlsFn(y,x,m=1)

LS of y on x; for one dependent variable

# Usage
(b,res,yhat,V,R2a) = OlsFn(y,x,m)

# Input
- `y::Array`:     Tx1, the dependent variable
- `x::Array`:     Txk matrix of regressors (including deterministic ones)
- `m::Int`:       scalar, bandwidth in Newey-West

# Output
- `b::Array`:     kx1, regression coefficients
- `u::Array`:     Tx1, residuals y - yhat
- `yhat::Array`:  Tx1, fitted values x*b
- `V::Array`:     kxk matrix, covariance matrix of sqrt(T)b
- `R2a::Number`:  scalar, R2 value

"""
function OlsFn(y,x,m=0)
    T    = size(y,1)
    b    = x\y
    yhat = x*b
    u    = y - yhat
    g    = x.*u
    S0   = NWFn(g,m)            #Newey-West covariance matrix
    Sxx  = -x'x/T
    V    = inv(Sxx)'S0*inv(Sxx)
    R2a  = 1 - var(u)/var(y)
    return b,u,yhat,V,R2a
end

max_stars_repo_name: snowdj/EmpiricalFinancePhD

max_issues_repo_name: AlTaken/EmpiricalFinancePhD

max_forks_repo_name: AlTaken/EmpiricalFinancePhD


text:
using TriMatrices
using TriMatrices: wraptri
using TriMatrices.Indexing: trinum
using LinearAlgebra

const _N = 4
const _data = 1:trinum(_N)
const m_upper = TriMatrix(TriUpper(), _N, _data)
const m_lower = TriMatrix(TriLower(), _N, _data)
const m_sym = TriMatrix(TriSymmetric(), _N, _data)

const _data_nd = 1:trinum(_N - 1)
const _diag = -1
const m_upper_nd = TriMatrix(TriUpper{false}(), _N, _data_nd; diag=_diag)
const m_lower_nd = TriMatrix(TriLower{false}(), _N, _data_nd; diag=_diag)
const m_sym_nd = TriMatrix(TriSymmetric{false}(), _N, _data_nd; diag=_diag)

@testset "Transpose" begin
	t_upper = transpose(m_upper)
	@test t_upper isa TriMatrix{Int, TriLower{true}}
	@test parent(t_upper) === parent(m_upper)

	t_lower = transpose(m_lower)
	@test t_lower isa TriMatrix{Int, TriUpper{true}}
	@test parent(t_lower) === parent(m_lower)

	@test transpose(m_sym) === m_sym

	t_upper_nd = transpose(m_upper_nd)
	@test t_upper_nd isa TriMatrix{Int, TriLower{false}}
	@test parent(t_upper_nd) === parent(m_upper_nd)
	@test t_upper_nd.diag == m_upper_nd.diag

	t_lower_nd = transpose(m_lower_nd)
	@test t_lower_nd isa TriMatrix{Int, TriUpper{false}}
	@test parent(t_lower_nd) === parent(m_lower_nd)
	@test t_lower_nd.diag == m_lower_nd.diag

	@test transpose(m_sym_nd) === m_sym_nd
end

@testset "wraptri" begin
	cases = [
		[m_upper, m_upper_nd] => UpperTriangular,
		[m_lower, m_lower_nd] => LowerTriangular,
		[m_sym, m_sym_nd] => Symmetric,
	]

	for (ms, T) in cases
		for m in ms
			wrapped = wraptri(m)
			@test wrapped isa T
			@test parent(wrapped) === m
		end
	end
end


text:
function scalar_field()
    # isoparametric approximation
    mesh = generate_grid(QuadraticQuadrilateral, (20, 20))
    f(x) = x[1]^2
    nodal_vals = [f(p.x) for p in mesh.nodes]

    ip_f = Lagrange{2,RefCube,2}() # function interpolation
    ip_g = Lagrange{2,RefCube,2}() # geometry interpolation

    # compute values in quadrature points
    qr = QuadratureRule{2, RefCube}(3) # exactly approximate quadratic field
    cv = CellScalarValues(qr, ip_f, ip_g)
    qp_vals = [Vector{Float64}(undef, getnquadpoints(cv)) for _ in 1:getncells(mesh)]
    for cellid in eachindex(mesh.cells)
        xe = getcoordinates(mesh, cellid)
        reinit!(cv, xe)
        for qp in 1:getnquadpoints(cv)
            qp_vals[cellid][qp] = f(spatial_coordinate(cv, qp, xe))
        end
    end

    # do a L2Projection for getting values in dofs
    projector = L2Projector(ip_f, mesh)
    projector_vals = project(projector, qp_vals, qr; project_to_nodes=false)

    # points where we want to retrieve field values
    points = [Vec((x, 0.52)) for x in range(0.0; stop=1.0, length=100)]

    # set up PointEvalHandler and retrieve values
    ph = PointEvalHandler(mesh, points)
    vals = get_point_values(ph, projector, projector_vals)
    @test f.(points) ≈ vals

    # alternatively retrieve vals from nodal values TODO: make this work?
    # vals = get_point_values(ph, nodal_vals)
    # @test f.(points) ≈ vals
end

function vector_field()
    ## vector field
    # isoparametric approximation
    mesh = generate_grid(QuadraticQuadrilateral, (20, 20))
    f(x) = Vec((x[1]^2, x[1]))
    nodal_vals = [f(p.x) for p in mesh.nodes]

    ip_f = Lagrange{2,RefCube,2}() # function interpolation
    ip_g = Lagrange{2,RefCube,2}() # geometry interpolation

    # compute values in quadrature points
    qr = QuadratureRule{2, RefCube}(3) # exactly approximate quadratic field
    cv = CellScalarValues(qr, ip_f, ip_g)
    qp_vals = [Vector{Vec{2,Float64}}(undef, getnquadpoints(cv)) for i=1:getncells(mesh)]
    for cellid in eachindex(mesh.cells)
        xe = getcoordinates(mesh, cellid)
        reinit!(cv, xe)
        for qp in 1:getnquadpoints(cv)
            qp_vals[cellid][qp] = f(spatial_coordinate(cv, qp, xe))
        end
    end

    # do a L2Projection for getting values in dofs
    projector = L2Projector(ip_f, mesh)
    projector_vals = project(projector, qp_vals, qr; project_to_nodes=false)
    # TODO: project_to_nodes should probably return dof values and not Vecs for vector fields
    # projector_vals = convert(Vector{Float64}, reinterpret(Float64, projector_vals))

    # points where we want to retrieve field values
    points = [Vec((x, 0.52)) for x in range(0.0; stop=1.0, length=100)]

    # set up PointEvalHandler and retrieve values
    ph = PointEvalHandler(mesh, points)
    vals = get_point_values(ph, projector, projector_vals)
    @test f.(points) ≈ vals

    # alternatively retrieve vals from nodal values# TODO
    # vals = get_point_values(ph, nodal_vals)
    # @test f.(points) ≈ vals
end

function superparametric()
    # superparametric approximation
    mesh = generate_grid(Quadrilateral, (20, 20))
    f(x) = x*x[1]
    ip_f = Lagrange{2,RefCube,2}() # function interpolation
    ip_g = Lagrange{2,RefCube,1}() # geometry interpolation

    # compute values in quadrature points
    qr = QuadratureRule{2, RefCube}(3) # exactly approximate quadratic field
    cv = CellScalarValues(qr, ip_f, ip_g)
    qp_vals = [Vector{Vec{2,Float64}}(undef, getnquadpoints(cv)) for i=1:getncells(mesh)]
    for cellid in eachindex(mesh.cells)
        xe = getcoordinates(mesh, cellid)
        reinit!(cv, xe)
        for qp in 1:getnquadpoints(cv)
            qp_vals[cellid][qp] = f(spatial_coordinate(cv, qp, xe))
        end
    end

    # do a L2Projection for getting values in dofs
    projector = L2Projector(ip_f, mesh)
    projector_vals = project(projector, qp_vals, qr; project_to_nodes=false)

    # points where we want to retrieve field values
    points = [Vec((x, 0.52)) for x in range(0.0; stop=1.0, length=100)]

    # set up PointEvalHandler and retrieve values
    ph = PointEvalHandler(mesh, points)
    vals = get_point_values(ph, projector, projector_vals)

    # can recover a quadratic field by a quadratic approximation
    @test f.(points) ≈ vals
end

function dofhandler()
    mesh = generate_grid(Quadrilateral, (2,2))
    dof_vals = [1., 2., 5., 4., 3., 6., 8., 7., 9.]
    points = [node.x for node in mesh.nodes] # same as nodes

    dh = DofHandler(mesh)
    push!(dh, :s, 1) # a scalar field
    close!(dh)

    ph = PointEvalHandler(mesh, points)
    vals = get_point_values(ph, dh, dof_vals, :s)
    @test vals ≈ 1.0:9.0

    # TODO
    # vals = get_point_values(ph, collect(1.0:9.0))
    # @test vals ≈ 1.0:9.0
end

function mixed_grid()
    ## Mixed grid where not all cells have the same fields

    # 5_______6
    # |\      |
    # |   \   |
    # 3______\4
    # |       |
    # |       |
    # 1_______2

    nodes = [Node((0.0, 0.0)),
            Node((1.0, 0.0)),
            Node((0.0, 1.0)),
            Node((1.0, 1.0)),
            Node((0.0, 2.0)),
            Node((1.0, 2.0))]

    cells = Ferrite.AbstractCell[Quadrilateral((1,2,4,3)),
            Triangle((3,4,6)),
            Triangle((3,6,5))]

    mesh = Grid(cells, nodes)
    addcellset!(mesh, "quads", Set{Int}((1,)))
    addcellset!(mesh, "tris", Set{Int}((2, 3)))

    ip_quad = Lagrange{2,RefCube,1}()
    ip_tri = Lagrange{2,RefTetrahedron,1}()

    f(x) = x[1]

    # compute values in quadrature points for quad
    qr = QuadratureRule{2, RefCube}(2)
    cv = CellScalarValues(qr, ip_quad)
    qp_vals_quads = [Vector{Float64}(undef, getnquadpoints(cv)) for cell in getcellset(mesh, "quads")]
    for (local_cellid, global_cellid) in enumerate(getcellset(mesh, "quads"))
        xe = getcoordinates(mesh, global_cellid)
        reinit!(cv, xe)
        for qp in 1:getnquadpoints(cv)
            qp_vals_quads[local_cellid][qp] = f(spatial_coordinate(cv, qp, xe))
        end
    end

    # construct projector
    projector = L2Projector(ip_quad, mesh; set=getcellset(mesh, "quads"))

    points = [Vec((x, 2x)) for x in range(0.0; stop=1.0, length=10)]

    # first alternative: L2Projection to dofs
    projector_values = project(projector, qp_vals_quads, qr; project_to_nodes = false)
    ph = PointEvalHandler(mesh, points)
    vals = get_point_values(ph, projector, projector_values)
    @test vals[1:5] ≈ f.(points[1:5])
    @test all(isnan, vals[6:end])
    # TODO
    # # second alternative: L2Projection to nodes
    # nodal_vals = project(projector, qp_vals_quads, qr; project_to_nodes = true)
    # vals = get_point_values(ph, nodal_vals)
    # @test vals[1:5] ≈ f.(points[1:5])
    # @test all(isnan.(vals[6:end]))

    # second alternative: assume a vector field :v
    dh = MixedDofHandler(mesh)
    field = Field(:v, ip_quad, 2)
    fh_quad = FieldHandler([field], getcellset(mesh, "quads"))
    push!(dh, fh_quad)
    field = Field(:v, ip_tri, 2)
    fh_tri = FieldHandler([field], getcellset(mesh, "tris"))
    push!(dh, fh_tri)
    close!(dh)

    dof_vals = [1., 1., 2., 2., 4., 4., 3., 3., 6., 6., 5., 5.]
    points = [node.x for node in mesh.nodes]
    ph = PointEvalHandler(mesh, points)
    vals = get_point_values(ph, dh, dof_vals, :v)
    @test vals == [Vec((i, i)) for i=1.0:6.0]
end

function oneD()
    # isoparametric approximation
    mesh = generate_grid(Line, (2,))
    f(x) = x[1]
    nodal_vals = [f(p.x) for p in mesh.nodes]

    ip_f = Lagrange{1,RefCube,1}() # function interpolation
    ip_g = Lagrange{1,RefCube,1}() # geometry interpolation

    # compute values in quadrature points
    qr = QuadratureRule{1, RefCube}(2)
    cv = CellScalarValues(qr, ip_f, ip_g)
    qp_vals = [Vector{Float64}(undef, getnquadpoints(cv)) for i=1:getncells(mesh)]
    for cellid in eachindex(mesh.cells)
        xe = getcoordinates(mesh, cellid)
        reinit!(cv, xe)
        for qp in 1:getnquadpoints(cv)
            qp_vals[cellid][qp] = f(spatial_coordinate(cv, qp, xe))
        end
    end

    # do a L2Projection for getting values in dofs
    projector = L2Projector(ip_f, mesh)
    projector_values = project(projector, qp_vals, qr; project_to_nodes=false)

    # points where we want to retrieve field values
    points = [Vec((x,)) for x in range(-1.0; stop=1.0, length=5)]

    # set up PointEvalHandler and retrieve values
    ph = PointEvalHandler(mesh, points)
    vals = get_point_values(ph, projector, projector_values)
    @test f.(points) ≈ vals

    # alternatively retrieve vals from nodal values
    # TODO
    # vals = get_point_values(ph, nodal_vals)
    # @test f.(points) ≈ vals
end

@testset "PointEvalHandler" begin
    scalar_field()
    vector_field()
    dofhandler()
    superparametric()
    mixed_grid()
    oneD()
end


text:
using HSA
using HSA.Emulation
using FactCheck

include("kernels.jl")

facts("The Emulation can be used to run kernels on the CPU") do
    context("vector add") do
        const alen = 10

        a = Array(Int, alen); rand!(a)
        b = Array(Int, alen); rand!(b)
        c = Array(Int, alen); rand!(c)
        c_expected = a + b

        run_cpu((alen, 1, 1), vadd, a,b,c)

        @fact c --> c_expected
    end

    context("matrix multiplication") do
        const arows = 3
        const acols = 4

        a = Array(Float64, arows, acols); rand!(a)
        b = Array(Float64, acols, arows); rand!(b)
        c = Array(Float64, arows, arows); rand!(c)
        c_expected = a * b

        run_cpu((arows, 1, 1), mmul, a,b,c,acols)

        @fact c --> c_expected
    end
end


text:
using FFTW
using DelimitedFiles

# using the convolutional theorem
function convolve_fft(signal1::Array{T}, signal2::Array{T}) where {T <: Number}
    return ifft(fft(signal1).*fft(signal2))
end

function main()

    # Random distribution in x
    x = rand(100)

    # Gaussian signals
    y = [exp(-((i-50)/100)^2/.01) for i = 1:100]

    # Normalization is not strictly necessary, but good practice
    normalize!(x)
    normalize!(y)

    # cyclic convolution via the convolutional theorem
    fft_output = convolve_fft(x, y)

    # outputting convolutions to different files for plotting in external code
    writedlm("fft.dat", fft_output)

end


text:
if !isdefined(:ExoplanetsSysSim) using ExoplanetsSysSim end
using Distributions
using Roots #if we want to use the 'find_zero()' function
using Optim

##### This code is a translation of 'MRpredict.R' at 'https://github.com/Bo-Ning/Predicting-exoplanet-mass-and-radius-relationship/blob/master/MR-predict/MRpredict.R'

struct MR_param_Ning2018
    Mass_min::Float64 #log10 scale
    Mass_max::Float64 #log10 scale
    Radius_min::Float64 #log10 scale
    Radius_max::Float64 #log10 scale
    deg::Int64
    weights_mle::Vector{Float64}
    indices_keep::Vector{Int64}
    #indices_div::Vector{Int64}
    #indices_rem::Vector{Int64}
end

function pdfnorm_beta(x::Float64, x_obs::Float64, x_sd::Float64, x_max::Float64, x_min::Float64, shape1::Int64, shape2::Int64 ; log::Bool=true)
    if log == true
        norm_dist = Normal(10^x, x_sd)
    else
        norm_dist = Normal(x, x_sd)
    end
    #beta_dist = Beta(shape1, shape2)
    return pdf.(norm_dist, x_obs) .* pdf.(Beta(shape1, shape2), (x-x_min)/(x_max-x_min)) / (x_max-x_min)
end

function fn_for_integrate(x_obs::Float64, x_sd::Float64, deg::Int, degree::Int, x_max::Float64, x_min::Float64 ; log::Bool=false, abs_tol::Float64=1e-10)
    shape1, shape2 = degree, deg-degree+1
    return quadgk(x -> pdfnorm_beta(x, x_obs, x_sd, x_max, x_min, shape1, shape2; log=log), x_min, x_max, abstol=abs_tol)[1]
end

function conditional_density(y::Float64, y_max::Float64, y_min::Float64, x_max::Float64, x_min::Float64, deg::Int64, w_hat::Vector{Float64}, indices_keep::Vector{Int64}; y_sd::Union{Void,Float64}=nothing, qtl::Vector{Float64}=[0.16, 0.84])
    deg_vec = 1:deg

    #To compute conditional mean, variance, quantile, distribution:
    y_beta_indv = zeros(deg)
    if y_sd == nothing
        #y_stdardize = (y-y_min)/(y_max-y_min)
        #y_beta_indv = map(degree -> pdf.(Beta(degree, deg-degree+1), (y-y_min)/(y_max-y_min)), deg_vec)::Array{Float64,1} / (y_max-y_min)
        map!(degree -> pdf.(Beta(degree, deg-degree+1), (y-y_min)/(y_max-y_min)), y_beta_indv, deg_vec)
        y_beta_indv /= (y_max-y_min)
    else
        #y_beta_indv = map(degree -> fn_for_integrate(y, y_sd, deg, degree, y_max, y_min), deg_vec)::Array{Float64,1}
        map!(degree -> fn_for_integrate(y, y_sd, deg, degree, y_max, y_min), y_beta_indv, deg_vec)
    end
    y_beta_pdf = kron(ones(deg), y_beta_indv)[indices_keep]
    ###y_beta_pdf = y_beta_indv[indices_rem]
    denominator = sum(w_hat .* y_beta_pdf::Vector{Float64})

    #Mean:
    mean_beta_indv = deg_vec/(deg+1)*(x_max-x_min)+x_min
    mean_beta = kron(mean_beta_indv, y_beta_indv)[indices_keep]
    ###mean_beta = mean_beta_indv[indices_div] .* y_beta_indv[indices_rem]
    mean_numerator = sum(w_hat .* mean_beta::Vector{Float64})
    mean = mean_numerator / denominator
    #mean = sum(w_hat .* kron(deg_vec/(deg+1)*(x_max-x_min)+x_min, y_beta_indv)) / denominator

    #Variance:
    var_beta_indv = deg_vec .* (deg-deg_vec+1) / ((deg+1)^2*(deg+2))*(x_max-x_min)^2
    var_beta = kron(var_beta_indv, y_beta_indv)[indices_keep]
    ###var_beta = var_beta_indv[indices_div] .* y_beta_indv[indices_rem]
    var_numerator = sum(w_hat .* var_beta::Vector{Float64})
    var = var_numerator /denominator
    #var = sum(w_hat .* kron(deg_vec .* (deg-deg_vec+1) / ((deg+1)^2*(deg+2))*(x_max-x_min)^2, y_beta_indv)) / denominator

    #Quantile:
    function pbeta_conditional_density(x::Float64)::Float64
        function mix_density(j::Float64)::Float64
            #x_indv_cdf::Vector{Float64} = map(degree -> cdf.(Beta(degree, deg-degree+1), (j-x_min)/(x_max-x_min)), deg_vec)
            ###quantile_beta = x_indv_cdf[indices_div] .* y_beta_indv[indices_rem]
            #quantile_numerator = sum(w_hat .* kron(x_indv_cdf, y_beta_indv)::Vector{Float64})
            #return quantile_numerator/denominator
            return sum(w_hat .* kron(map(degree -> cdf.(Beta(degree, deg-degree+1), (j-x_min)/(x_max-x_min)), deg_vec), y_beta_indv)[indices_keep]::Vector{Float64}) / denominator
        end

        return mix_density(x)
    end

    function conditional_quantile(q::Float64)::Float64
        function g(x::Float64)::Float64
            return pbeta_conditional_density(x) - q
        end
        function g2(x::Float64)::Float64
            return (pbeta_conditional_density(x) - q)^2
        end

        #return find_zero(g, (x_min, x_max), FalsePosition())
        #return find_zero(g, (x_min, x_max))::Float64 ###'find_zero' is from 'using Roots' (is this a good function to use?); Eric suggested brent's method in the Optim package

        result = optimize(g2, x_min, x_max, Brent(), rel_tol=1e-3, abs_tol=1e-3) #Brent(); GoldenSection()
        if Optim.converged(result)
            return Optim.minimizer(result)::Float64
        else
            return find_zero(g, (x_min, x_max))::Float64
        end
    end

    quantile = map(conditional_quantile, qtl)

    return (mean, var, quantile, denominator, y_beta_indv)
end

function predict_mass_given_radius(Radius::Float64, param::MR_param_Ning2018 ; R_sigma::Union{Void,Float64}=nothing, posterior_sample::Bool=false, qtl::Vector{Float64}=[0.16, 0.84])
    #Convert data to log scale:
    l_radius = log10.(Radius)

    #The 'posterior_sample == false' condition can deal with two cases:
    #Case I: if input data do not have measurement errors
    #Case II: if input data have measurement error
    if posterior_sample == false
        predicted_value = conditional_density(l_radius, param.Radius_max, param.Radius_min, param.Mass_max, param.Mass_min, param.deg, param.weights_mle, param.indices_keep; y_sd=R_sigma, qtl=qtl)
        predicted_mean = predicted_value[1]
        predicted_quantiles = predicted_value[3]
        #predicted_lower_quantile = predicted_value[3][1]
        #predicted_upper_quantile = predicted_value[3][2]
    elseif posterior_sample == true
        #Case III: if the input are posterior samples
        radius_sample = log10(Radius)
        k = length(radius_sample)
        denominator_sample = zeros(k)
        mean_sample = zeros(k)
        y_beta_indv_sample = zeros(k, param.deg)

        #Calculate mean:
        for i in 1:k
            results = cond_density_estimation(y=radius_sample[i], y_max=param.Radius_max, y_min=param.Radius_min, x_max=param.Mass_max, x_min=param.Mass_min, deg=param.deg, w_hat=param.weights_mle, qtl=quantile, only_output_mean=true) ###what function is this??? Is it part of some CRAN package...
            mean_sample[i] = (results[1]) ###why are there parentheses here? It makes no difference...
            denominator_sample[i] = results[2]
            y_beta_indv_sample[i,:] = results[3:57]
        end
        predicted_mean = mean(mean_sample)

        #Calculate 16% and 84% quantiles:
        #Mixture of the CDF of k conditional densities
        function pbeta_conditional_density(x)
            function mix_density(j)
                deg_vec = 1:param.deg
                x_indv_cdf = map(degree -> cdf.(Beta(degree, param.deg-degree+1), (j-x_min)/(x_max-x_min)), deg_vec)
                quantile_numerator = zeros(k)
                p_beta_sample = zeros(k)
                for ii in 1:k
                    quantile_numerator[ii] = sum(param.weights_mle .* kron(x_indv_cdf, y_beta_indv_sample[ii,:]))
                    p_beta_sample[ii] = quantile_numerator[ii]/denominator_sample[ii]
                end
                return p_beta = mean(p_beta_sample)
            end

            return map(mix_density, x)
        end

        function mixture_conditional_quantile(q, x_min, x_max)
            function g(x)
                return pbeta_conditional_density(x) - q
            end
            return find_zero(g, (x_min, x_max)) ###see note above
            #return Optim.minimum(optimize(g, x_min, x_max))
        end

        predicted_quantiles = map(q -> mixture_conditional_quantile(q, param.Mass_min, param.Mass_max), qtl)
        #predicted_lower_quantile = predicted_quantiles[1]
        #predicted_upper_quantile = predicted_quantiles[2]
    end

    #Return the output:
return (predicted_mean, predicted_quantiles)::Tuple{Float64, Vector{Float64}}
end

function draw_planet_mass_from_radius_Ning2018(Radius::Float64, param::MR_param_Ning2018)
    #This function takes in a Radius (in solar radii) and draws a mass (returning in solar masses) probabilistically
    Radius_in_earths = Radius/ExoplanetsSysSim.earth_radius
    @assert param.Radius_min < log10.(Radius_in_earths) < param.Radius_max
    q::Float64 = rand()
    l_mass = predict_mass_given_radius(Radius_in_earths, param; qtl=[q])[2][1]
    return (10^l_mass)*ExoplanetsSysSim.earth_mass
end

function generate_planet_mass_from_radius_Ning2018(Radius::Float64, sim_param::SimParam)
    global MR_param
    return draw_planet_mass_from_radius_Ning2018(Radius, MR_param)
end

##### To initialize the M-R model parameters:

weights_mle = CSV.read("weights.mle.csv")[:x]
degrees = 55

N_keep = 25
indices_keep = sortperm(weights_mle, rev=true)[1:N_keep]
#indices_div = 1+div.(indices_keep .- 1, degrees)
#indices_rem = 1+rem.(indices_keep .- 1, degrees)
weights = weights_mle[indices_keep]
MR_param = MR_param_Ning2018(-1., 3.809597, -0.302, 1.357509, degrees, weights, indices_keep) #(-1., 3.809597, -0.3, 1.357509, 55, weights_mle)

##### Examples:

#=
#Observation without measurement errors:

Radius = 5. #original scale, not log scale
predict_result = predict_mass_given_radius(Radius, MR_param)
println(predict_result)

#Observation with a measurement error:

Radius = 5. #original scale, not log scale
R_sigma = 0.1
predict_result = predict_mass_given_radius(Radius, MR_param; R_sigma=0.1)
println(predict_result)

#If want to change 16% and 84% quantiles to 5% and 95% quantiles:

Radius = 5. #original scale, not log scale
predict_result = predict_mass_given_radius(Radius, MR_param; qtl=[0.05, 0.95])
println(predict_result)
=#

#Input are posterior samples: ###currently broken because the function 'cond_density_estimation' is undefined
#=
Radius_dist = Normal(5, 0.5) #original scale, not log scale
Radius_samples = rand(Radius_dist, 100)
predict_result = predict_mass_given_radius(Radius_samples, MR_param; R_sigma=0.1, posterior_sample=true)
println(predict_result)
=#

##### For timing the function:
#=
sim_param = setup_sim_param_model()
#Radii = ones(10000)*ExoplanetsSysSim.earth_radius
Radii = (10.^(linspace(MR_param.Radius_min+0.01, MR_param.Radius_max-0.01, 10000)))*ExoplanetsSysSim.earth_radius
tic()
Masses = map(r -> generate_planet_mass_from_radius_Ning2018(r, sim_param), Radii)
t_elapsed = toc()
=#

##### For plotting and timing the model:
#=
#Pkg.add("PyPlot")
using PyPlot

N_keep_array = [3025, 25, 15]
colors = ["blue", "green", "red"]
Radii = 10.^(linspace(MR_param.Radius_min+0.01, MR_param.Radius_max-0.01, 1000))

fig, ax = subplots()
for (i,N_keep) in enumerate(N_keep_array)
    indices_keep = sortperm(weights_mle, rev=true)[1:N_keep]
    weights = weights_mle[indices_keep]
    MR_param = MR_param_Ning2018(-1., 3.809597, -0.302, 1.357509, degrees, weights, indices_keep)

    tic()
    Masses_predict = map(r -> predict_mass_given_radius(r, MR_param; qtl=[0.16, 0.84]), Radii)
    t_elapsed = toc()

    Masses_mean = 10.^[x[1] for x in Masses_predict]
    Masses_lower = 10.^[x[2][1] for x in Masses_predict]
    Masses_upper = 10.^[x[2][2] for x in Masses_predict]

    loglog(Radii, Masses_mean, color=colors[i], linewidth=2.0, linestyle="-", label="M = "*string(N_keep))
    loglog(Radii, Masses_lower, color=colors[i], linewidth=1.0, linestyle="--")
    loglog(Radii, Masses_upper, color=colors[i], linewidth=1.0, linestyle="--")
end
xlabel(L"$R_p (R_\oplus)$")
ylabel(L"$M_p (M_\oplus)$")
legend(loc="upper left")
=#


text:
#%%
using Glimpse
const Gl = Glimpse
import Glimpse: SystemKind, System, Component, ComponentData, update, Spatial, TimingData,
                valid_entities, singleton, component, RGBAf0
using LinearAlgebra

# define rotation component and system
struct Rotation <: ComponentData
    omega::Float32
    center::Point3f0
    axis::Vec3f0
end
RotationComponent(id) = Component(id, Rotation)

struct Rotator <: SystemKind end
rotator_system(dio) = System{Rotator}(dio, (Spatial, Rotation), (TimingData,))

function update(sys::System{Rotator})
    rotation = component(sys, Rotation)
    spatial  = component(sys, Spatial)
    dt       = Float32(singleton(sys, TimingData).dtime)
    for i in valid_entities(rotation, spatial)
        e_rotation = rotation[i]
        n          = e_rotation.axis
        r          = -e_rotation.center + spatial[i].position
        theta      = e_rotation.omega * dt
        nnd        = n * dot(n, r)
        spatial[i] = Spatial(Point3f0(e_rotation.center + nnd + (r - nnd) * cos(theta) + cross(r, n) * sin(theta)), spatial[i].velocity)
    end
end

begin
    dio = Diorama(; background = RGBAf0(0.0, 0.0, 0.0, 1.0))

    # add new component and system to the diorama, alongside with the already present rendering systems/components
    Gl.add_component!(dio, Rotation)
    Gl.add_system!(dio, rotator_system(dio))
    # add some sphere entities and the new extra component which will automatically engage with the new system
    geom = Gl.PolygonGeometry(Sphere(Point3f0(0.0), 1.0f0))
    progtag = Gl.ProgramTag{Gl.PeelingProgram}()
    for i in 1:2000
        Gl.new_entity!(dio;
                       separate = [Gl.UniformColor(rand(RGBAf0)),
                                   Gl.Spatial(50 * (1 .- rand(Point3f0)), rand(Vec3f0)),
                                   Gl.Material(), Gl.Shape(),
                                   Rotation(0.2f0, 50 * (1 .- rand(Point3f0)),
                                            normalize(1 .- rand(Vec3f0))), Gl.Dynamic(),
                                   progtag], shared = [geom])
    end
    dio.loop = @async Gl.renderloop(dio)
end
#%%


text:
abstract type A end
struct B <: A
    i::Int
end

expr1 = :(struct Foo1 i::B end)

class = B
expr2 = :(struct Foo2 i::$class end)

dump(expr1)

dump(expr2)


text:
const barchars =  ["", "▏","▎","▍","▌","▋","▊","▉","█"]

function barspec(progress, maxsubs, maxchars)
    maxpoints = maxchars * (maxsubs)

    @assert 0.0 <= progress <= 1.0

    points = round(Int, progress * maxpoints)
    fulls = points ÷ maxsubs
    subs = mod(points, maxsubs)
    @assert fulls <= maxchars

    (fulls, subs)
end

function baronly(progress, maxchars = 50, barchars = barchars)
    if isfinite(progress)
        progress = clamp(progress, 0.0, 1.0)

        maxsubs = length(barchars) - 1
        fulls, subs = barspec(progress, maxsubs, maxchars)
        rems = maxchars - (fulls + 1) #+1 is for the subs
        bar = barchars[end]^fulls * barchars[subs+1]
    else
        bar=string(" ", progress, " ")
    end
    @assert length(bar) <= maxchars

    rpad(bar, maxchars)
end

function format_progress(progress)
    if isfinite(progress)
        "$(round(100progress, digits=2))%"
    else
        string(progress)
    end
end

function bar(progress; maxchars = 50, barchars = barchars)
    bartext = baronly(progress, maxchars, barchars)
    "[$bartext] " * format_progress(progress)
end

function progress_message(progress_raw)
    progress_val(x) = x
    function progress_val(x::AbstractString)
        inner(::Nothing) = x
        inner(mm::RegexMatch) = parse(Float32, mm[1])/100

        inner(match(r"^\s*(\d*\.?\d*)\s*\%\s*$", x))
        # a number followed by a %, with whitespace whereever
    end

    message(x) = x
    message(x::Real) = bar(x)

    message(progress_val(progress_raw))
end


text:
module ContractionOrder

export ContractionTree

struct ContractionTree
    left
    right
end

function log2sumexp2(s)
    ms = maximum(s)
    return log2(sum(x->exp2(x - ms), s)) + ms
end

include("incidencelist.jl")
include("greedy.jl")

end

using .ContractionOrder: IncidenceList, ContractionTree, contract_pair!, MinSpaceOut, MinSpaceDiff, tree_greedy
export parse_eincode, ContractionOrder, optimize_greedy

function parse_eincode!(::Type{ET}, ::IncidenceList, tree, vertices_order, level=0) where ET
    ti = findfirst(==(tree), vertices_order)
    ti, NestedEinsum{ET}(ti)
end

function parse_eincode!(::Type{EIT}, incidence_list::IncidenceList, tree::ContractionTree, vertices_order, level=0) where EIT
    ti, codei = parse_eincode!(EIT, incidence_list, tree.left, vertices_order, level+1)
    tj, codej = parse_eincode!(EIT, incidence_list, tree.right, vertices_order, level+1)
    dummy = Dict([e=>0 for e in keys(incidence_list.e2v)])
    _, _, code = contract_pair!(incidence_list, vertices_order[ti], vertices_order[tj], dummy)
    ti, NestedEinsum((codei, codej), _eincode(EIT, [code.first...], level==0 ? incidence_list.openedges : code.second))
end
_eincode(::Type{DynamicEinCode{LT}}, ixs::Vector, iy::Vector) where LT = DynamicEinCode(ixs, iy)
_eincode(::Type{StaticEinCode}, ixs::Vector, iy::Vector) = StaticEinCode{(Tuple.(ixs)...,), (iy...,)}()

function parse_eincode(::Type{EIT}, incidence_list::IncidenceList{VT,ET}, tree::ContractionTree; vertices = collect(keys(incidence_list.v2e))) where {EIT,VT,ET}
    parse_eincode!(EIT, copy(incidence_list), tree, vertices)[2]
end

function parse_tree(ein, vertices)
    if isleaf(ein)
        vertices[ein.tensorindex]
    else
        if length(ein.args) != 2
            error("This eincode is not a binary tree.")
        end
        left, right = parse_tree.(ein.args, Ref(vertices))
        ContractionTree(left, right)
    end
end

"""
    optimize_greedy(eincode, size_dict; method=MinSpaceOut(), nrepeat=10)

Greedy optimizing the contraction order and return a `NestedEinsum` object. Methods are
* `MinSpaceOut`, always choose the next contraction that produces the minimum output tensor.
* `MinSpaceDiff`, always choose the next contraction that minimizes the total space.
"""
function optimize_greedy(code::DynamicEinCode{L}, size_dict::Dict; method=MinSpaceOut(), nrepeat=10) where {L}
    optimize_greedy(DynamicEinCode{L}, getixsv(code), getiyv(code), size_dict; method=method, nrepeat=nrepeat)
end
function optimize_greedy(code::StaticEinCode, size_dict::Dict; method=MinSpaceOut(), nrepeat=10)
    optimize_greedy(StaticEinCode, getixsv(code), getiyv(code), size_dict; method=method, nrepeat=nrepeat)
end
function optimize_greedy(::Type{ET}, ixs::AbstractVector{<:AbstractVector}, iy::AbstractVector, size_dict::Dict{L,TI}; method=MinSpaceOut(), nrepeat=10) where {ET, L, TI}
    if length(ixs) <= 2
        return NestedEinsum((1:length(ixs)...,), _eincode(ET, ixs, iy))
    end
    log2_edge_sizes = Dict{L,Float64}()
    for (k, v) in size_dict
        log2_edge_sizes[k] = log2(v)
    end
    incidence_list = ContractionOrder.IncidenceList(Dict([i=>ixs[i] for i=1:length(ixs)]); openedges=iy)
    tree, _, _ = tree_greedy(incidence_list, log2_edge_sizes; method=method, nrepeat=nrepeat)
    parse_eincode!(ET, incidence_list, tree, 1:length(ixs))[2]
end
function optimize_greedy(code::NestedEinsum, size_dict; method=MinSpaceOut(), nrepeat=10)
    isleaf(code) && return code
    args = optimize_greedy.(code.args, Ref(size_dict); method=method, nrepeat=nrepeat)
    if length(code.args) > 2
        # generate coarse grained hypergraph.
        nested = optimize_greedy(code.eins, size_dict; method=method, nrepeat=nrepeat)
        replace_args(nested, args)
    else
        NestedEinsum(args, code.eins)
    end
end

function replace_args(nested::NestedEinsum{ET}, trueargs) where ET
    isleaf(nested) && return trueargs[nested.tensorindex]
    NestedEinsum(replace_args.(nested.args, Ref(trueargs)), nested.eins)
end

export timespace_complexity
"""
    timespace_complexity(eincode, size_dict)

Return the time and space complexity of the einsum contraction.
The time complexity is defined as `log2(number of element multiplication)`.
The space complexity is defined as `log2(size of the maximum intermediate tensor)`.
"""
function timespace_complexity(ei::NestedEinsum, size_dict)
    log2_sizes = Dict([k=>log2(v) for (k,v) in size_dict])
    _timespace_complexity(ei, log2_sizes)
end

function timespace_complexity(ei::EinCode, size_dict)
    log2_sizes = Dict([k=>log2(v) for (k,v) in size_dict])
    _timespace_complexity(getixsv(ei), getiyv(ei), log2_sizes)
end

function _timespace_complexity(ei::NestedEinsum, log2_sizes::Dict{L}) where L
    isleaf(ei) && return (-Inf, -Inf)
    tcs = Float64[]
    scs = Float64[]
    for arg in ei.args
        tc, sc = _timespace_complexity(arg, log2_sizes)
        push!(tcs, tc)
        push!(scs, sc)
    end
    tc2, sc2 = _timespace_complexity(getixsv(ei.eins), getiyv(ei.eins), log2_sizes)
    tc = ContractionOrder.log2sumexp2([tcs..., tc2])
    sc = max(reduce(max, scs), sc2)
    return tc, sc
end

function _timespace_complexity(ixs::AbstractVector, iy::AbstractVector{T}, log2_sizes::Dict{L}) where {T, L}
    # remove redundant legs
    counts = Dict{L,Int}()
    for ix in ixs
        for l in ix
            if haskey(counts, l)
                counts[l] += 1
            else
                counts[l] = 1
            end
        end
    end
    for l in iy
        if haskey(counts, l)
            counts[l] += 1
        else
            counts[l] = 1
        end
    end
    loop_inds = L[]
    for ix in ixs
        for l in ix
            c = count(==(l), ix)
            if counts[l] > c && l ∉ loop_inds
                push!(loop_inds, l)
            end
        end
    end
    tc = isempty(loop_inds) ? -Inf : sum(l->log2_sizes[l], loop_inds)
    sc = isempty(iy) ? 0.0 : sum(l->log2_sizes[l], iy)
    return tc, sc
end