text:
#############################################################################################
#                             COPY TO BLOCK MATRICES                                       #
#############################################################################################

"Copy vectors from a vector of vectors to a concatenated vector"
function copy_inds(dest, src, inds)
    for i in eachindex(inds)
        dest[inds[i]] = src[i]
    end
end

function copyback_inds(dest, src, inds)
    for i in eachindex(inds)
        src[i] = dest[inds[i]]
    end
end

"Copy constraints to a single concatenated vector"
function copy_constraints!(d, solver::ConstrainedSolver)
    conSet = get_constraints(solver)
    for (i,con) in enumerate(conSet.errvals)
        copy_inds(d, con.vals, solver.con_inds[i])
    end
    return nothing
end

function copy_multipliers!(λ, solver::ConstrainedSolver)
    conSet = get_constraints(solver)
    for (i,con) in enumerate(conSet.errvals)
        copy_inds(λ, con.λ, solver.con_inds[i])
    end
    return nothing
end

function copyback_multipliers!(λ, solver::ConstrainedSolver)
    conSet = get_constraints(solver)
    for (i,con) in enumerate(conSet.errvals)
        copyback_inds(λ, con.λ, solver.con_inds[i])
    end
    return nothing
end

"Copy active set to a single concatenated vector"
function copy_active_set!(a, solver::ConstrainedSolver)
    conSet = get_constraints(solver)
    for i = 1:length(conSet.errvals)
        copy_inds(solver.active_set, conSet.convals[i].active, solver.con_inds[i])
    end
end

"""Copy constraint Jacobians to given indices in a sparse array
Dispatches on bandedness of the constraint"""
function copy_jacobian!(D, con::ALConVal{<:TO.StageConstraint}, cinds, xinds, uinds) where T
    for (i,k) in enumerate(con.inds)
        zind = [xinds[k]; uinds[k]]
        D[cinds[i], zind] .= con.jac[i]
    end
end

function copy_jacobian!(D, con::ALConVal{<:TO.StateConstraint}, cinds, xinds, uinds) where T
    for (i,k) in enumerate(con.inds)
        D[cinds[i], xinds[k]] .= con.jac[i]
    end
end

function copy_jacobian!(D, con::ALConVal{<:TO.ControlConstraint}, cinds, xinds, uinds) where T
    for (i,k) in enumerate(con.inds)
        D[cinds[i], uinds[k]] .= con.jac[i]
    end
end

# function copy_jacobian!(D, con::ConstraintVals{T,Coupled}, cinds, xinds, uinds) where T
#     for (i,k) in enumerate(con.inds)
#         zind = [xinds[k]; uinds[k]; xinds[k+1]; uinds[k+1]]
#         D[cinds[i], zind] .= con.∇c[i]
#     end
# end

function copy_jacobian!(D, con::TO.AbstractConstraintValues{<:DynamicsConstraint{<:Explicit}},
		cinds, xinds, uinds)
	N = length(xinds)
    for (i,k) in enumerate(con.inds)
        zind = [xinds[k]; uinds[k]; xinds[k+1]]
		zind1 = [xinds[k]; uinds[k]]
		zind2 = [xinds[k+1]; uinds[k+1]]
        D[cinds[i], zind1] .= con.jac[i,1]
		D[cinds[i], xinds[k+1]] .= con.jac[i,2][:,xinds[1]]
		# D[cinds[i], zind2] .= con.jac[i,2]
    end
end

"Copy all constraint Jacobians to a sparse matrix"
function copy_jacobians!(D, solver::ConstrainedSolver)
    conSet = get_constraints(solver)
    xinds, uinds = primal_partition(solver)
    cinds = solver.con_inds

    for i = 1:length(conSet.errvals)
        copy_jacobian!(D, conSet.convals[i], cinds[i], xinds, uinds) # TODO: switch to errvals
    end
    return nothing
end

"Copy constraint Jacobians to linear indices of a vector"
function copy_jacobian!(d::AbstractVector{<:Real}, con::TO.AbstractConstraintValues, linds)
	for (j,k) in enumerate(con.inds)
		inds = linds[j]
		d[inds] = con.jac[j]
	end
end

"Copy all constraint Jacobians to linear indices of a vector"
function copy_jacobians!(jac::AbstractVector{<:Real}, solver::ConstrainedSolver)
    conSet = get_constraints(solver)
    xinds, uinds = primal_partition(solver)
    cinds = solver.con_inds
    linds = jacobian_linear_inds(solver)

    for (i,con) in enumerate(conSet.constraints)
		copy_jacobian!(jac, con, linds[i])
    end
    return nothing
end

# function copy_gradient!(grad_f, E::Vector{<:CostExpansion}, xinds, uinds)
#     for k = 1:length(uinds)
#         grad_f[xinds[k]] = E[k].x
#         grad_f[uinds[k]] = E[k].u
#     end
#     if length(xinds) != length(uinds)
#         grad_f[xinds[end]] = E[end].x
#     end
# end


text:
using JSON, Colors, Memoize

# Grabbing the open league rankings from PvPoke. These are common enough that
# automatically downloading these and making them constant makes sense.
const gamemaster =     JSON.parsefile(download(
    "https://raw.githubusercontent.com/pvpoke/pvpoke/master/src/data/" *
    "gamemaster.json"))

"""
    get_cp_limit(league)

Given a league ("great", "ultra", "master") return its cp limit
(1500, 2500, 10000). Masters has no cp limit but is signified in PvPoke with the
arbitrarily large 10000.

# Examples
```jldoctest
julia> get_cp_limit("great")
1500
"""
get_cp_limit(league::String) = league == "master" ? 10_000 :
    league == "ultra" ? 2_500 : 1_500

"""
    get_rankings(cup; league = "great")

Download the PvPoke rankings for a particular meta, optionally within a league
that is not great league. This function is memoized to avoid downloading the
same file multiple times.
"""
@memoize function get_rankings(cup::String; league = "great")
    return JSON.parsefile(download(
        "https://raw.githubusercontent.com/pvpoke/pvpoke/master/src/" *
        "data/rankings/$(cup)/overall/rankings-$(get_cp_limit(league)).json"))
end

"""
    get_overrides(cup; league = "great")

Download the PvPoke overrides used to generate rankings for a particular meta,
optionally within a league that is not great league. This function is memoized
to avoid downloading the same file multiple times.
"""
@memoize function get_overrides(cup::String; league = "great")
    cup == "great" && return get_overrides("overall", league = "great")
    cup == "ultra" && return get_overrides("overall", league = "ultra")
    cup == "master" && return get_overrides("overall", league = "master")
    return JSON.parsefile(download(
        "https://raw.githubusercontent.com/pvpoke/pvpoke/master/src/data/" *
        "overrides/$(cup)/overall/$(get_cp_limit(league)).json"))
end

@memoize get_silph_usage(cup::String) =
    JSON.parsefile(download("https://silph.gg/api/cup/" * cup * "/stats/.json"))

"""
    get_gamemaster_mon_id(name)

Given a mon's PvPoke id, get its location within the gamemaster JSON file.
"""
function get_gamemaster_mon_id(name::String)
    for i = 1:length(gamemaster["pokemon"])
        gamemaster["pokemon"][i]["speciesId"] == name && return i
    end
end

"""
    get_gamemaster_move_id(name)

Given a move's PvPoke id, get its location within the gamemaster JSON file.
"""
function get_gamemaster_move_id(name::String)
    for i = 1:length(gamemaster["moves"])
        gamemaster["moves"][i]["moveId"] == name && return i
    end
end

"""
    get_rankings_mon_id(name; league = "great", cup = "all")

Given a mon's PvPoke id, find its location within a particular meta's rankings
"""
function get_rankings_mon_id(name::String;
    league::String = "great", cup::String = "all")
    rankings = get_rankings(cup, league = league)
    for i = 1:length(rankings)
        rankings[i]["speciesId"] == name && return i
    end
    return 0
end

const typings = ("normal", "fighting", "flying", "poison", "ground",
    "rock", "bug", "ghost", "steel", "fire", "water", "grass",
    "electric", "psychic", "ice", "dragon", "dark", "fairy", "none")

# CP multipliers from PvPoke
const cpm = Dict(
    1.0 => 0.0939999967813492, 1.5 => 0.135137432089339,
    2.0 => 0.166397869586945, 2.5 => 0.192650913155325,
    3.0 => 0.215732470154762, 3.5 => 0.236572651424822,
    4.0 => 0.255720049142838, 4.5 => 0.273530372106572,
    5.0 => 0.290249884128571, 5.5 => 0.306057381389863,
    6.0 => 0.321087598800659, 6.5 => 0.335445031996451,
    7.0 => 0.349212676286697, 7.5 => 0.362457736609939,
    8.0 => 0.375235587358475, 8.5 => 0.387592407713878,
    9.0 => 0.399567276239395, 9.5 => 0.4111935532161,
    10.0 => 0.422500014305115, 10.5 => 0.432926420512509,
    11.0 => 0.443107545375824, 11.5 => 0.453059948165049,
    12.0 => 0.46279838681221, 12.5 => 0.472336085311278,
    13.0 => 0.481684952974319, 13.5 => 0.490855807179549,
    14.0 => 0.499858438968658, 14.5 => 0.5087017489616,
    15.0 => 0.517393946647644, 15.5 => 0.525942516110322,
    16.0 => 0.534354329109192, 16.5 => 0.542635753803599,
    17.0 => 0.550792694091797, 17.5 => 0.558830584490385,
    18.0 => 0.566754519939423, 18.5 => 0.57456912814537,
    19.0 => 0.582278907299042, 19.5 => 0.589887907888945,
    20.0 => 0.597400009632111, 20.5 => 0.604823648665171,
    21.0 => 0.61215728521347, 21.5 => 0.619404107958234,
    22.0 => 0.626567125320435, 22.5 => 0.633649178748576,
    23.0 => 0.6406529545784, 23.5 => 0.647580971386554,
    24.0 => 0.654435634613037, 24.5 => 0.661219265805859,
    25.0 => 0.667934000492095, 25.5 => 0.674581885647492,
    26.0 => 0.681164920330048, 26.5 => 0.687684901255373,
    27.0 => 0.694143652915955, 27.5 => 0.700542901033063,
    28.0 => 0.706884205341339, 28.5 => 0.713169074873823,
    29.0 => 0.719399094581604, 29.5 => 0.725575586915154,
    30.0 => 0.731700003147125, 30.5 => 0.734741038550429,
    31.0 => 0.737769484519958, 31.5 => 0.740785579737136,
    32.0 => 0.743789434432983, 32.5 => 0.746781197247765,
    33.0 => 0.749761044979095, 33.5 => 0.752729099732281,
    34.0 => 0.75568550825119, 34.5 => 0.758630370209851,
    35.0 => 0.761563837528229, 35.5 => 0.76448604959218,
    36.0 => 0.767397165298462, 36.5 => 0.770297293677362,
    37.0 => 0.773186504840851, 37.5 => 0.776064947064992,
    38.0 => 0.778932750225067, 38.5 => 0.781790050767666,
    39.0 => 0.784636974334717, 39.5 => 0.787473608513275,
    40.0 => 0.790300011634827, 40.5 => 0.792803968023538,
    41.0 => 0.795300006866455, 41.5 => 0.797803898371622,
    42.0 => 0.800300002098083, 42.5 => 0.802803871877596,
    43.0 => 0.805299997329711, 43.5 => 0.807803850847053,
    44.0 => 0.81029999256134, 44.5 => 0.812803835179168,
    45.0 => 0.815299987792968, 45.5 => 0.817803806620319,
    46.0 => 0.820299983024597, 46.5 => 0.822803778631297,
    47.0 => 0.825299978256225, 47.5 => 0.827803750922782,
    48.0 => 0.830299973487854, 48.5 => 0.832803753381377,
    49.0 => 0.835300028324127, 49.5 => 0.837803755931569,
    50.0 => 0.840300023555755
)

# Color scheme to match PvPoke
const colors = [RGBA(153/255, 159/255, 161/255, 1.0),
                RGBA(213/255,  63/255,  91/255, 1.0),
                RGBA(148/255, 171/255, 225/255, 1.0),
                RGBA(193/255,  98/255, 212/255, 1.0),
                RGBA(212/255, 141/255,  91/255, 1.0),
                RGBA(208/255, 196/255, 142/255, 1.0),
                RGBA(158/255, 195/255,  49/255, 1.0),
                RGBA( 89/255, 107/255, 181/255, 1.0),
                RGBA( 82/255, 142/255, 160/255, 1.0),
                RGBA(254/255, 163/255,  84/255, 1.0),
                RGBA( 86/255, 158/255, 222/255, 1.0),
                RGBA( 94/255, 189/255,  91/255, 1.0),
                RGBA(246/255, 215/255,  75/255, 1.0),
                RGBA(245/255, 126/255, 121/255, 1.0),
                RGBA(120/255, 212/255, 192/255, 1.0),
                RGBA( 14/255, 104/255, 184/255, 1.0),
                RGBA( 86/255,  86/255,  99/255, 1.0),
                RGBA(240/255, 152/255, 228/255, 1.0)]

# Purple-ish shield color to match game
const shieldColor = RGBA(235/255,13/255,199/255, 1.0)

"""
    silph_to_pvpoke(name)

Given the name of a pokemon from the Silph API, return the pvpoke id

# Examples
```jldoctest
julia> silph_to_pvpoke("Mr. Mime")
mr_mime
"""
function silph_to_pvpoke(name::String)
    name = lowercase(name)
    name = replace(name, "-" => "_")
    name = replace(name, "_normal" => "")
    name = replace(name, "mr. " => "mr_")
    name = replace(name, "alola" => "alolan")
    name = replace(name, "galar" => "galarian")
    name = replace(name, "sunshine" => "sunny")
    name = replace(name, "basculin_red" => "basculin")
    name = replace(name, "basculin_blue" => "basculin")
    name = replace(name, "darmanitan" => "darmanitan_standard")
    name = replace(name, "gastrodon_blue" => "gastrodon_west_sea")
    name = replace(name, "gastrodon_pink" => "gastrodon_east_sea")
    name = replace(name, "sirfetch'd" => "sirfetchd")
    return name
end

ext: jl


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

using BAT
using Test

using LinearAlgebra, Random, StableRNGs
using ArraysOfArrays, ValueShapes, Distributions, PDMats, StatsBase

struct _TestDensityStruct <: AbstractDensity
end

struct _UniformDensityStruct <: AbstractDensity
end

struct _DeltaDensityStruct <: AbstractDensity
end

struct _ShapeDensityStruct <: AbstractDensity
end

@testset "abstract_density" begin
    mvn = MvNormal(ones(3), PDMat(Matrix{Float64}(I,3,3)))
    ValueShapes.totalndof(td::_TestDensityStruct) = Int(3)
    BAT.sampler(td::_TestDensityStruct) = BAT.sampler(mvn)
    BAT.eval_logval_unchecked(density::_TestDensityStruct, v::Any) = Distributions.logpdf(mvn, v)

    td = _TestDensityStruct()
    tds = BAT.DensityWithShape(td, ArrayShape{Real}(3))

    @test @inferred(isequal(@inferred(BAT.var_bounds(td)), missing))
    @test @inferred(isequal(@inferred(varshape(td)), missing))

    x = rand(3)
    @test_throws ArgumentError BAT.eval_logval(td, [Inf, Inf, Inf], BAT.default_dlt())
    @test_throws BAT.DensityEvalException BAT.eval_logval(tds, [Inf, Inf, Inf], BAT.default_dlt())
    @test_throws BAT.DensityEvalException BAT.eval_logval(tds, [NaN, NaN, NaN], BAT.default_dlt())
    @test_throws ArgumentError BAT.eval_logval(tds, rand(length(mvn)+1), BAT.default_dlt())
    @test_throws ArgumentError BAT.eval_logval(tds, rand(length(mvn)-1), BAT.default_dlt())

    @test @inferred(BAT.eval_logval(tds, x, BAT.default_dlt())) == @inferred(logpdf(mvn, x))

    mvu = product_distribution([Uniform() for i in 1:3])
    BAT.varshape(ud::_UniformDensityStruct) = varshape(mvu)
    BAT.eval_logval_unchecked(ud::_UniformDensityStruct, v::Any) = logpdf(mvu, v)
    BAT.var_bounds(ud::_UniformDensityStruct) = BAT.HyperRectBounds(BAT.HyperRectVolume(zeros(3), ones(3)))
    ValueShapes.totalndof(ud::_UniformDensityStruct) = Int(3)

    x = [-Inf, 0, Inf]
    ud_shape_1 = NamedTupleShape(a=ArrayShape{Real}(1), b=ArrayShape{Real}(1), c=ArrayShape{Real}(1))
    ud_shape_2 = NamedTupleShape(a=ArrayShape{Real}(3))
    ud = _UniformDensityStruct()

    @test @inferred(BAT.eval_logval(ud, x, BAT.default_dlt())) == -Inf
    @test @inferred(BAT.eval_logval(ud_shape_1(ud), ud_shape_1(x), BAT.default_dlt())) == @inferred(logpdf(mvu, x))
    @test @inferred(BAT.eval_logval(ud_shape_2(ud), ud_shape_2(x), BAT.default_dlt())) == @inferred(logpdf(mvu, x))

    @test_throws ArgumentError BAT.eval_logval(ud, vcat(x,x), BAT.default_dlt())

    @test BAT.eval_logval(ud, x .- eps(1.0), BAT.default_dlt()) == -Inf

    @test_throws ArgumentError @inferred(BAT.eval_logval(ud, [0 0 0], BAT.default_dlt()))

    ntshape = NamedTupleShape(a=ScalarShape{Real}(), b=ScalarShape{Real}(), c=ScalarShape{Real}())
    shapedasnt = ShapedAsNT(x, ntshape)

    cvd = ConstValueDist(0)
    ValueShapes.totalndof(dd::_DeltaDensityStruct) = Int(1)
    BAT.eval_logval_unchecked(dd::_DeltaDensityStruct, v::Any) = Distributions.logpdf(cvd, v)

    dd = _DeltaDensityStruct()
    dds = BAT.DensityWithShape(dd, ScalarShape{Real}())
    @test_throws ArgumentError BAT.eval_logval(dd, 0, BAT.default_dlt())
    @test_throws BAT.DensityEvalException BAT.eval_logval(dds, 0, BAT.default_dlt())

    ntdist = NamedTupleDist(a=mvn, b=mvu)
    ValueShapes.varshape(sd::_ShapeDensityStruct) = varshape(ntdist)
    BAT.eval_logval_unchecked(sd::_ShapeDensityStruct, v) = logpdf(ntdist, v)

    x1_for_sd = rand(length(ntdist.a))
    x2_for_sd = rand(length(ntdist.b))
    x_for_sd = vcat(x1_for_sd, x2_for_sd)

    correct_shape_of_sd = NamedTupleShape(a=ArrayShape{Real}(length(ntdist.a)), b=ArrayShape{Real}(length(ntdist.b)))
    incorrect_shape_of_sd = NamedTupleShape(a=ArrayShape{Real}(length(ntdist.a)-1), b=ArrayShape{Real}(length(ntdist.b)+1))

    x_for_sd_good_shape = correct_shape_of_sd(x_for_sd)
    x_for_sd_bad_shape = incorrect_shape_of_sd(x_for_sd)

    sd = _ShapeDensityStruct()

    @test BAT.eval_logval(sd, x_for_sd_good_shape, BAT.default_dlt()) == logpdf(mvn, x1_for_sd) + logpdf(mvu, x2_for_sd)
    @test_throws ArgumentError BAT.eval_logval(sd, x_for_sd_bad_shape, BAT.default_dlt())

    @testset "rand" begin
        td = _TestDensityStruct()
        @test rand(StableRNG(7002), BAT.bat_sampler(td)) ≈ [2.386799038, 1.072161895, 0.791486531]
    end
end


text:
using DataFrames
using Test
using Random: randstring

@testset "sort it" begin
M = 100_000
str_base = [randstring(8) for i in 1:1_000]
df = DataFrame(int = rand(Int32, M), float=rand(M), str = rand(str_base, M))

@time df1 = sort(df, :int);
@time df2 = fsort(df, :int);
@test df1 == df2
@test df != df2

@time df1 = sort(df, :str);
@time df2 = fsort(df, :str);
@test df1 == df2
@test df != df2

@time df1 = sort(df, [:str, :float]);
@time df2 = fsort(df, [:str, :float]);
@test df1 == df2
@test df != df2
end

max_stars_repo_name: xiaodaigh/DataConvenience

max_issues_repo_name: xiaodaigh/DataConvenience

max_forks_repo_name: xiaodaigh/DataConvenience


text:
## A SortedDict is a wrapper around balancedTree with
## methods similiar to those of Julia container Dict.

@compat type SortedDict{K, D, Ord <: Ordering} <: Associative{K,D}
    bt::BalancedTree23{K,D,Ord}

## Zero-argument constructor, or possibly one argument to specify order.

    function (::Type{SortedDict{K,D,Ord}}){K, D, Ord <: Ordering}(o::Ord=Forward)
        bt1 = BalancedTree23{K,D,Ord}(o)
        new{K,D,Ord}(bt1)
    end

end

## external constructor to take an associative and infer
## argument types

function SortedDict{K, D, Ord <: Ordering}(d::Associative{K,D}, o::Ord=Forward)
    h = SortedDict{K,D,Ord}(o)
    for (k,v) in d
        h[k] = v
    end
    h
end

## More constructors based on those in dict.jl:
## Take pairs and infer argument
## types.  Note:  this works only for the Forward ordering.

function SortedDict{K,D}(p1::Pair{K,D}, ps::Pair{K,D}...)
    h = SortedDict{K,D,ForwardOrdering}()
    h[p1.first] = p1.second
    for p in ps
        h[p.first] = p.second
    end
    h
end

## Take pairs and infer argument
## types.  Ordering parameter must be explicit first argument.

function SortedDict{K,D, Ord <: Ordering}(o::Ord, ps::Pair{K,D}...)
    h = SortedDict{K,D,Ord}(o)
    for p in ps
        h[p.first] = p.second
    end
    h
end

## This one takes an iterable; ordering type is optional.

SortedDict{Ord <: Ordering}(kv, o::Ord=Forward) =
sorteddict_with_eltype(kv, eltype(kv), o)

function sorteddict_with_eltype{K,D,Ord}(kv, ::Type{Pair{K,D}}, o::Ord)
    h = SortedDict{K,D,Ord}(o)
    for (k,v) in kv
        h[k] = v
    end
    h
end

const SDSemiToken = IntSemiToken

const SDToken = Tuple{SortedDict,IntSemiToken}

## This function implements m[k]; it returns the
## data item associated with key k.

@inline function getindex(m::SortedDict, k_)
    i, exactfound = findkey(m.bt, convert(keytype(m),k_))
    !exactfound && throw(KeyError(k_))
    return m.bt.data[i].d
end

## This function implements m[k]=d; it sets the
## data item associated with key k equal to d.

@inline function setindex!{K, D, Ord <: Ordering}(m::SortedDict{K,D,Ord}, d_, k_)
    insert!(m.bt, convert(K,k_), convert(D,d_), false)
    m
end

## push! is an alternative to insert!; it returns the container.

@inline function push!{K, D, Ord <: Ordering}(m::SortedDict{K,D,Ord}, pr::Pair{K,D})
    insert!(m.bt, convert(K,pr[1]), convert(D,pr[2]), false)
    m
end

## This function looks up a key in the tree;
## if not found, then it returns a marker for the
## end of the tree.

@inline function find(m::SortedDict, k_)
    ll, exactfound = findkey(m.bt, convert(keytype(m),k_))
    IntSemiToken(exactfound? ll : 2)
end

## This function inserts an item into the tree.
## Unlike m[k]=d, it also returns a bool and a token.
## The bool is true if the inserted item is new.
## It is false if there was already an item
## with that key.
## The token points to the newly inserted item.

@inline function insert!{K,D, Ord <: Ordering}(m::SortedDict{K,D,Ord}, k_, d_)
    b, i = insert!(m.bt, convert(K,k_), convert(D,d_), false)
    b, IntSemiToken(i)
end

@inline eltype{K,D,Ord <: Ordering}(m::SortedDict{K,D,Ord}) =  Pair{K,D}
@inline eltype{K,D,Ord <: Ordering}(::Type{SortedDict{K,D,Ord}}) =  Pair{K,D}
@inline function in{K,D,Ord <: Ordering}(pr::Pair, m::SortedDict{K,D,Ord})
    i, exactfound = findkey(m.bt,convert(K,pr[1]))
    return exactfound && isequal(m.bt.data[i].d,convert(D,pr[2]))
end

@inline in(::Tuple{Any,Any}, ::SortedDict) =
    throw(ArgumentError("'(k,v) in sorteddict' not supported in Julia 0.4 or 0.5.  See documentation"))

@inline keytype{K,D,Ord <: Ordering}(m::SortedDict{K,D,Ord}) = K
@inline keytype{K,D,Ord <: Ordering}(::Type{SortedDict{K,D,Ord}}) = K
@inline valtype{K,D,Ord <: Ordering}(m::SortedDict{K,D,Ord}) = D
@inline valtype{K,D,Ord <: Ordering}(::Type{SortedDict{K,D,Ord}}) = D
@inline ordtype{K,D,Ord <: Ordering}(m::SortedDict{K,D,Ord}) = Ord
@inline ordtype{K,D,Ord <: Ordering}(::Type{SortedDict{K,D,Ord}}) = Ord

## First and last return the first and last (key,data) pairs
## in the SortedDict.  It is an error to invoke them on an
## empty SortedDict.

@inline function first(m::SortedDict)
    i = beginloc(m.bt)
    i == 2 && throw(BoundsError())
    return Pair(m.bt.data[i].k, m.bt.data[i].d)
end

@inline function last(m::SortedDict)
    i = endloc(m.bt)
    i == 1 && throw(BoundsError())
    return Pair(m.bt.data[i].k, m.bt.data[i].d)
end

@inline orderobject(m::SortedDict) = m.bt.ord

@inline function haskey(m::SortedDict, k_)
    i, exactfound = findkey(m.bt,convert(keytype(m),k_))
    exactfound
end

@inline function get{K,D,Ord <: Ordering}(m::SortedDict{K,D,Ord}, k_, default_)
    i, exactfound = findkey(m.bt, convert(K,k_))
   return exactfound? m.bt.data[i].d : convert(D,default_)
end

function get!{K,D,Ord <: Ordering}(m::SortedDict{K,D,Ord}, k_, default_)
    k = convert(K,k_)
    i, exactfound = findkey(m.bt, k)
    if exactfound
        return m.bt.data[i].d
    else
        default = convert(D,default_)
        insert!(m.bt,k, default, false)
        return default
    end
end

function getkey{K,D,Ord <: Ordering}(m::SortedDict{K,D,Ord}, k_, default_)
    i, exactfound = findkey(m.bt, convert(K,k_))
    exactfound? m.bt.data[i].k : convert(K, default_)
end

## Function delete! deletes an item at a given
## key

@inline function delete!(m::SortedDict, k_)
    i, exactfound = findkey(m.bt,convert(keytype(m),k_))
    !exactfound && throw(KeyError(k_))
    delete!(m.bt, i)
    m
end

@inline function pop!(m::SortedDict, k_)
    i, exactfound = findkey(m.bt,convert(keytype(m),k_))
    !exactfound && throw(KeyError(k_))
    d = m.bt.data[i].d
    delete!(m.bt, i)
    d
end

## Check if two SortedDicts are equal in the sense of containing
## the same (K,D) pairs.  This sense of equality does not mean
## that semitokens valid for one are also valid for the other.

function isequal(m1::SortedDict, m2::SortedDict)
    ord = orderobject(m1)
    if !isequal(ord, orderobject(m2)) || !isequal(eltype(m1), eltype(m2))
        throw(ArgumentError("Cannot use isequal for two SortedDicts unless their element types and ordering objects are equal"))
    end
    p1 = startof(m1)
    p2 = startof(m2)
    while true
        if p1 == pastendsemitoken(m1)
            return p2 == pastendsemitoken(m2)
        end
        if p2 == pastendsemitoken(m2)
            return false
        end
        k1,d1 = deref((m1,p1))
        k2,d2 = deref((m2,p2))
        if !eq(ord,k1,k2) || !isequal(d1,d2)
            return false
        end
        p1 = advance((m1,p1))
        p2 = advance((m2,p2))
    end
end

function mergetwo!{K,D,Ord <: Ordering}(m::SortedDict{K,D,Ord},
                                        m2::Associative{K,D})
    for (k,v) in m2
        m[convert(K,k)] = convert(D,v)
    end
end

function packcopy{K,D,Ord <: Ordering}(m::SortedDict{K,D,Ord})
    w = SortedDict(Dict{K,D}(),orderobject(m))
    mergetwo!(w,m)
    w
end

function packdeepcopy{K,D,Ord <: Ordering}(m::SortedDict{K,D,Ord})
    w = SortedDict(Dict{K,D}(),orderobject(m))
    for (k,v) in m
        newk = deepcopy(k)
        newv = deepcopy(v)
        w[newk] = newv
    end
    w
end

function merge!{K,D,Ord <: Ordering}(m::SortedDict{K,D,Ord},
                                     others::Associative{K,D}...)
    for o in others
        mergetwo!(m,o)
    end
end

function merge{K,D,Ord <: Ordering}(m::SortedDict{K,D,Ord},
                                    others::Associative{K,D}...)
    result = packcopy(m)
    merge!(result, others...)
    result
end

similar{K,D,Ord<:Ordering}(m::SortedDict{K,D,Ord}) =
SortedDict{K,D,Ord}(orderobject(m))


text:
using caffeCalJuli
using Test

@testset "caffeCalJuli.jl" begin
    # Write your tests here.
end

max_stars_repo_name: Kaikei-e/caffeCalJuli

max_issues_repo_name: Kaikei-e/caffeCalJuli

max_forks_repo_name: Kaikei-e/caffeCalJuli


text:
import Flatten: flattenable
@metadata initial_value nothing
import FieldMetadata: @units, units
import FieldMetadata: @limits, limits
import FieldMetadata: @prior, prior
import FieldMetadata: @description, description
@metadata bounds nothing
import FieldMetadata: @logscaled, logscaled
@metadata reference nothing

"""
    AbstractParameters{T} <: AbstractVector{T}

An abstract type for AIBECS model parameters.

Parameters in AIBECS use the following convenience packages:

- Parameters
- FieldMetadata
- FieldDefaults
- Flatten
- Unitful
- DataFrames
- Distributions

These aim to allow for some nice features, which include

- nice syntax for unpacking parameters in functions via the `@unpack` macro (fron UnPack.jl)
- additional metadata on parameters
- easy conversion to and from vectors
- use of units and automatic conversions if necessary
- pretty table-format displays
- loading and saving to and from CSV files
- prior estimates for bayesian inference and optimization

See the list of examples to get an idea of how to generate parameters for your model.

# Examples

Generate a simple parameter type via

```jldoctest
julia> struct SimpleParams{T} <: AbstractParameters{T}
           α::T
           β::T
           γ::T
       end
SimpleParams
```

To create an instance of the `SimpleParams(Float64)` type, you can do

```jldoctest
julia> p = SimpleParams(1.0, 2.0, 3.0)
SimpleParams{Float64}
│ Row │ Symbol │ Value   │
│     │ Symbol │ Float64 │
├─────┼────────┼─────────┤
│ 1   │ α      │ 1.0     │
│ 2   │ β      │ 2.0     │
│ 3   │ γ      │ 3.0     │
```

One of the core features from Parameters is unpacking in functions, e.g.,

```jldoctest
julia> function simplef(p)
           @unpack α, γ = p
           return α + γ
       end
simplef (generic function with 1 method)

julia> simplef(p) # 1.0 + 3.0
4.0
```

More complex examples are permitted by adding metadata (thanks to FieldMetadata.jl).
You can add units

```jldoctest
julia> @units struct UnitParams{T} <: AbstractParameters{T}
           α::T | u"km"
           β::T | u"hr"
           γ::T | u"m/s"
       end ;

julia> p = UnitParams(1.0, 2.0, 3.0)
UnitParams{Float64}
│ Row │ Symbol │ Value   │ Unit     │
│     │ Symbol │ Float64 │ Unitful… │
├─────┼────────┼─────────┼──────────┤
│ 1   │ α      │ 1.0     │ km       │
│ 2   │ β      │ 2.0     │ hr       │
│ 3   │ γ      │ 3.0     │ m s^-1   │
```

Note that when adding units to your parameters, they will be converted to SI
when unpacked, as in, e.g.,

```jldoctest
julia> function speed(p)
           @unpack α, β, γ = p
           return α / β + γ
       end
speed (generic function with 1 method)

julia> speed(p) # (1.0 km / 2.0 hr + 3 m/s) in m/s
3.138888888888889
```

Another example for optimizable/flattenable parameters

```jldoctest
julia> @initial_value @units @flattenable struct OptParams{T} <: AbstractParameters{T}
           α::T | 3.6 | u"km"  | true
           β::T | 1.0 | u"hr"  | false
           γ::T | 1.0 | u"m/s" | true
       end ;

julia> p = OptParams(initial_value(OptParams)...)
OptParams{Float64}
│ Row │ Symbol │ Value   │ Initial value │ Unit     │ Optimizable │
│     │ Symbol │ Float64 │ Float64       │ Unitful… │ Bool        │
├─────┼────────┼─────────┼───────────────┼──────────┼─────────────┤
│ 1   │ α      │ 3.6     │ 3.6           │ km       │ 1           │
│ 2   │ β      │ 1.0     │ 1.0           │ hr       │ 0           │
│ 3   │ γ      │ 1.0     │ 1.0           │ m s^-1   │ 1           │
```

Thanks to the FieldMetaData interface, you can chain the following preloaded metadata:

- initial_value
- units (from Unitful.jl)
- prior (from Distributions.jl)
- description (`String`)
- bounds (2-element `Tuple`)
- logscaled (`Bool`)
- flattenable (to convert to vectors of optimizable parameters only)
- reference (`String`)

Here is an example of parameter with all the possible metadata available in AIBECS:

```jldoctest
julia> @initial_value @units @prior @description @bounds @logscaled @flattenable @reference struct FullParams{T} <: AbstractParameters{T}
           α::T | 1.0 | u"km"  | Normal(0,1)    | "The distance"   | (-Inf, Inf) | false | false | "Jean et al., 2042"
           β::T | 2.0 | u"hr"  | LogNormal(0,1) | "The time"       | (   0, Inf) | true  | true  | "Claude et al. 1983"
           γ::T | 3.0 | u"mol" | Normal(1,2)    | "The # of moles" | (  -1,   1) | false | true  | "Dusse et al. 2000"
       end ;

julia> FullParams(4.0, 5.0, 6.0)
FullParams{Float64}
│ Row │ Symbol │ Value   │ Initial value │ Unit     │ Prior                            │ Description    │ Bounds      │ Logscaled │ Optimizable │ Reference          │
│     │ Symbol │ Float64 │ Float64       │ Unitful… │ Distribu…                        │ String         │ Tuple…      │ Bool      │ Bool        │ String             │
├─────┼────────┼─────────┼───────────────┼──────────┼──────────────────────────────────┼────────────────┼─────────────┼───────────┼─────────────┼────────────────────┤
│ 1   │ α      │ 4.0     │ 1.0           │ km       │ Normal{Float64}(μ=0.0, σ=1.0)    │ The distance   │ (-Inf, Inf) │ 0         │ 0           │ Jean et al., 2042  │
│ 2   │ β      │ 5.0     │ 2.0           │ hr       │ LogNormal{Float64}(μ=0.0, σ=1.0) │ The time       │ (0, Inf)    │ 1         │ 1           │ Claude et al. 1983 │
│ 3   │ γ      │ 6.0     │ 3.0           │ mol      │ Normal{Float64}(μ=1.0, σ=2.0)    │ The # of moles │ (-1, 1)     │ 0         │ 1           │ Dusse et al. 2000  │
```

Note that there is no check that the metadata you give is consistent.
These metadata will hopefully be useful for advanced usage of AIBECS, e.g., using prior information and/or bounds for optimization.
"""
abstract type AbstractParameters{T} <: AbstractVector{T} end

"""
    symbols(p)

Returns the symbols in `p`.

Can also be used directly on the type of `p`
(because the symbols of `p::T` are contained in the type `T`).
"""
symbols(::Type{T}) where {T <: AbstractParameters} = fieldnames(T)
symbols(::T) where {T <: AbstractParameters} = symbols(T)

"""
    flattenable_symbols(p)

Returns the flattenable symbols in `p`.

The flattenable symbols are those symbols that are kepth when using `p` as a vector,
e.g., when doing `vec(p)`.
(Useful when passing parameters to an optimization routine expeting a vector of optimizable parameters.)

Can also be used directly on the type of `p`
(because the flattenable symbols of `p::T` are contained in the type `T`).
"""
flattenable_symbols(::T) where {T <: AbstractParameters} = flattenable_symbols(T)
flattenable_symbols(::Type{T}) where {T <: AbstractParameters} = fieldnames(T)[collect(flattenable(T))]

"""
    values(p::T) where {T <: AbstractParameters}

Returns a vector of **all** the values of `p`.

Note that `values(p)` is different from `vec(p)`.
"""
Base.values(p::T) where {T <: AbstractParameters} = [getfield(p, s) for s in symbols(T)]

"""
    flattenable_values(p::T) where {T <: AbstractParameters}

Returns a vector of the **flattenable** values of `p`.

Note that `vec(p)` is different from `values(p)`.
"""
flattenable_values(p::T) where {T <: AbstractParameters} = [getfield(p, s) for s in flattenable_symbols(T)]

"""
    vec(p::T) where {T <: AbstractParameters}

Returns a **SI-unit-converted** vector of flattenable values of `p`.

Note that `vec(p) ≠ flattenable_values(p)` if `p` has units.
"""
Base.vec(p::T) where {T <: AbstractParameters} = [UnPack.unpack(p, Val(s)) for s in flattenable_symbols(T)]

"""
    length(p::AbstractParameter)

Returns the length of the **flattened/optimzable** vector of `p`.

May be different from the number of parameters.
Can also be used directly on the type of `p`.
"""
Base.length(::T) where {T <: AbstractParameters} = sum(flattenable(T))
Base.length(::Type{T}) where {T <: AbstractParameters} = sum(flattenable(T))

"""
    size(p::AbstractParameter)

Returns the size of the **flattened/optimzable** vector of `p`.

May be different from the number of parameters.
Can also be used directly on the type of `p`.
"""
Base.size(::T) where {T <: AbstractParameters} = (length(T),)
Base.size(::Type{T}) where {T <: AbstractParameters} = (length(T),)

function Base.show(io::IO, p::T) where T <: AbstractParameters
    print(T)
    t = table(p)
    show(io, t; summary=false)
end
function Base.show(io::IO, m::MIME"text/plain", p::T) where T <: AbstractParameters
    print(T)
    t = table(p)
    show(io, m, t; summary=false)
end

function table(p::T, nf::NamedTuple) where {T <: AbstractParameters}
    t = DataFrame(Symbol = collect(symbols(p)), Value = values(p))
    for (n,f) in zip(keys(nf), nf)
        setproperty!(t, n, collect(f(p)))
    end
    return t
end

"""
    table(p)

Returns a `DataFrame` (a table) of `p`.
Useful for printing and saving into an actual text/latex table.
"""
function table(p::AbstractParameters)
    ks, vs = Symbol[], Function[]
    all(isnothing, initial_value(p)) || (push!(ks, Symbol("Initial value")); push!(vs, initial_value))
    all(isequal(1), units(p)) || (push!(ks, :Unit); push!(vs, units))
    all(isnothing, prior(p)) || (push!(ks, Symbol("Prior")); push!(vs, prior))
    all(isempty, description(p)) || (push!(ks, Symbol("Description")); push!(vs, description))
    all(isnothing, bounds(p)) || (push!(ks, Symbol("Bounds")); push!(vs, bounds))
    any(logscaled(p)) && (push!(ks, Symbol("Logscaled")); push!(vs, logscaled))
    all(flattenable(p)) || (push!(ks, Symbol("Optimizable")); push!(vs, flattenable))
    all(isnothing, reference(p)) || (push!(ks, Symbol("Reference")); push!(vs, reference))
    nf = (; zip(ks, vs)...)
    t = table(p, nf)
end

"""
    latex(p)

Returns a LaTeX-formatted table of the parameters.
"""
latex(p::AbstractParameters) = show(stdout, MIME("text/latex"), table(p))

"""
    unpack(p <: AbstractParameters, s)

Unpacks the parameter `s` from `p`.

Note this is specialized and will convert the parameter value to SI units.
"""
@inline UnPack.unpack(p::T, ::Val{f}) where {T<:AbstractParameters,f} = ustrip(upreferred(getproperty(p, f) * units(T, f)))

"""
    +(p::T, v::Vector) where {T <: AbstractParameters}

Adds the flattened vector `v` to `p`.

**Warning:** This method for `+` is implemented only for differentiation
using dual and hyperdual numbers.
If you want to change the values of `p`, you should do so explicitly
rather than use this `+` method.
"""
Base.:+(p::T, v::Vector) where {T <: AbstractParameters} = reconstruct(T, vec(p) + v)

"""
    reconstruct(T, v)

Reconstructs the parameter of type `T` from flattenable vector `v`.
"""
function reconstruct(::Type{T}, v::Tv) where {T <: AbstractParameters, Tv <: Vector}
    all(isnothing, initial_value(T)) && error("Can't reconstruct without initial values")
    vunits = units(T)[[flattenable(T)...]]
    v = @. v * upreferred(vunits) |> vunits |> ustrip
    reconstructed_v = convert(Tv, collect(initial_value(T)))
    reconstructed_v[collect(flattenable(T))] .= v
    return T.name.wrapper(reconstructed_v...)
end

"""
    getindex(p::T, i) where {T <: AbstractParameters}

Returns the i-th element of vec(p).

This is not efficient and only used for testing the derivatives with ForwardDiff.
"""
Base.getindex(p::T, i) where {T <: AbstractParameters} = getindex(vec(p), i)

function (::Type{T})(args::Quantity...) where {T <: AbstractParameters}
    all(isequal(1), units(T)) && error("$T needs `units` for this construction to work")
    return T([ustrip(x |> units(T, f)) for (x,f) in zip(args, fieldnames(T))]...)
end

function (::Type{T})(;kwargs...) where {T <: AbstractParameters}
    all(isnothing, initial_value(T)) && length(kwargs) ≠ length(symbols(T)) && error("$T needs `initial_value` if one of the parameters is not supplied")
    value(f::Symbol, v::Quantity) = ustrip(v |> units(T, f))
    value(f::Symbol, v) = v
    return T([f ∈ keys(kwargs) ? value(f, kwargs[f]) : initial_value(T, f) for f in fieldnames(T)]...)
end

#===============================
Writing to savable formats
===============================#

function Base.Dict(p::T, s=symbols(p)) where {T<:AbstractParameters}
    v = [getfield(p,s) for s in s]
    u = [units(p,s) for s in s]
    return Dict([(s,v*u) for (s,v,u) in zip(s,v,u)])
end
function Base.NamedTuple(p::T, s=symbols(p)) where {T<:AbstractParameters}
    v = [getfield(p,s) for s in s]
    u = [units(p,s) for s in s]
    return (; zip(s, v .* u)...)
end

#=====================
mismatch of parameters
=====================#
"""
    mismatch(p::AbstractParameters)

Returns the sum of the negative log-likelihood of each flattenable parameter.
"""
function mismatch(p::AbstractParameters)
    return sum(-logpdf(prior(p,k), UnPack.unpack(p,Val(k))) for k in flattenable_symbols(p))
end
function ∇mismatch(p::AbstractParameters)
    return transpose([-gradlogpdf(prior(p,k), UnPack.unpack(p,Val(k))) for k in flattenable_symbols(p)])
end

# The functions below is just for ForwardDiff to work with vectors instead of p
# which requires `mismatch` to know about the priors, which are containted in `T`
function generate_objective(ωs, μx, σ²x, v, ωp, ::Type{T}) where {T<:AbstractParameters}
    nt, nb = length(ωs), length(v)
    tracers(x) = state_to_tracers(x, nb, nt)
    f(x, p) = ωp * mismatch(T, p) +
        sum([ωⱼ * mismatch(xⱼ, μⱼ, σⱼ², v) for (ωⱼ, xⱼ, μⱼ, σⱼ²) in zip(ωs, tracers(x), μx, σ²x)])
    return f
end
function mismatch(::Type{T}, v) where {T<:AbstractParameters}
    return sum(-logpdf(prior(T,k), v[i]) for (i,k) in enumerate(flattenable_symbols(T)))
end

#==================
Change of variables
==================#
"""
    subfun

Returns the substitution function for the change of variables of parameters.

If the prior of parameter `pᵢ` is `LogNormal`, then the substitution function is `exp`.
If the prior is `Uniform`, then the change of variables is the logit function.
Otherwise, it's `identity`.
"""
function subfun(::Type{T}) where {T<:AbstractParameters}
    return λ -> reconstruct(T, [subfun(T, s)(λᵢ) for (λᵢ,s) in zip(λ, flattenable_symbols(T))])
end
function ∇subfun(::Type{T}) where {T<:AbstractParameters}
    return λ -> reconstruct(T, [∇subfun(T, s)(λᵢ) for (λᵢ,s) in zip(λ, flattenable_symbols(T))])
end
function ∇²subfun(::Type{T}) where {T<:AbstractParameters}
    return λ -> reconstruct(T, [∇²subfun(T, s)(λᵢ) for (λᵢ,s) in zip(λ, flattenable_symbols(T))])
end
function invsubfun(::Type{T}) where {T<:AbstractParameters}
    return p -> [invsubfun(T, s)(pᵢ) for (pᵢ,s) in zip(vec(p), flattenable_symbols(T))]
end
# substitution function (change of variables) is determined from prior distribution
subfun(::Type{T}, s::Symbol) where {T<:AbstractParameters} = subfun(prior(T,s))
∇subfun(::Type{T}, s::Symbol) where {T<:AbstractParameters} = ∇subfun(prior(T,s))
∇²subfun(::Type{T}, s::Symbol) where {T<:AbstractParameters} = ∇²subfun(prior(T,s))
invsubfun(::Type{T}, s::Symbol) where {T<:AbstractParameters} = invsubfun(prior(T,s))
# Fallback rule for change of variables is identity
subfun(::Distribution) = identity
∇subfun(::Distribution) = x -> one(x)
∇²subfun(::Distribution) = x -> zero(x)
invsubfun(::Distribution) = identity
# p = exp(λ) for LogNormal
subfun(::LogNormal) = exp
∇subfun(::LogNormal) = exp
∇²subfun(::LogNormal) = exp
invsubfun(::LogNormal) = log
# p = logistic(λ) for Uniform
subfun(d::Uniform) = λ -> d.a + (d.b - d.a) / (exp(-λ) + 1)
∇subfun(d::Uniform) = λ -> (d.b - d.a) * exp(-λ) / (exp(-λ) + 1)^2
∇²subfun(d::Uniform) = λ -> (d.a - d.b) * exp(-λ) / (exp(-λ) + 1)^2 + 2(d.b - d.a) * exp(-2λ) / (exp(-λ) + 1)^3
invsubfun(d::Uniform) = p -> -log((d.b - d.a) / (p - d.a) - 1)

export subfun, ∇subfun, ∇²subfun, invsubfun

export AbstractParameters, latex


text:
# The Layer abstract type and empty function definitions, to be expanded by subtypes.
# To be included in MERAKit.jl

"""
Abstract supertype of all layer types, e.g. `BinaryLayer` and `TernaryLayer`.
"""
abstract type Layer end

layertype(::Type{L}) where {L<:Layer} = L

"""
    baselayertype(::Type{<:Layer})

Return the unparametrised layer type, e.g. for `BinaryLayer{ComplexSpace, Float64, false}`,
return just `BinaryLayer`.
"""
baselayertype(::Type{L}) where {L<:Layer} = Base.typename(L).wrapper

"""
    randomlayer(::Type{T <: Layer}, T, Vin, Vout, Vint = Vout; random_disentangler = false)

Return a MERA layer with Haar random tensors.

`T` is the `Layer` type, and `Vin` and `Vout` are the input and output spaces. `Vint` is an
internal vector space for the layer, connecting the disentanglers to the isometries. If
`random_disentangler = true`, the disentangler is also a random unitary, if `false`
(default), it is the identity or the product of two single-site isometries, depending on if
the disentangler is supposed to be unitary or isometric.

Each subtype of `Layer` should have its own method for this function.
"""
function randomlayer end

"""
    outputspace(layer::Layer)

Return the vector space of the downwards-pointing (towards the physical level) indices of
`layer`.

See also: [`inputspace`](@ref), [`internalspace`](@ref)
"""
function outputspace end

"""
    inputspace(layer::Layer)

Return the vector space of the upwards-pointing (towards the scale invariance) indices of
`layer`.

See also: [`outputspace`](@ref), [`internalspace`](@ref)
"""
function inputspace end

"""
    internalspace(layer::Layer)

Return the internal vector space of `layer`. This need not be defined for all layer types.

See also: [`outputspace`](@ref),  [`inputspace`](@ref)
"""
function internalspace end

"""
    expand_inputspace(layer::Layer, V_new)

Return a new layer where the tensors have been padded with zeros as necessary to change the
input space to be `V_new`.

Each subtype of `Layer` should have its own method for this function.
"""
function expand_inputspace end

"""
    expand_outputspace(layer::Layer, V_new)

Return a new layer where the tensors have been padded with zeros as necessary to change the
output space to be `V_new`.

Each subtype of `Layer` should have its own method for this function.
"""
function expand_outputspace end

"""
    expand_internalspace(layer::Layer, V_new)

Return a new layer where the tensors have been padded with zeros as necessary to change the
interal vector space to be `V_new`.

Each subtype of `Layer` that has an internal vector space should have its own method for
this function.
"""
function expand_internalspace end

"""
    space_invar_intralayer(layer::Layer)

Return `true` if the indices within `layer` are compatible with each other, false otherwise.
"""
function space_invar_intralayer end

"""
    space_invar_interlayer(layer::T, next_layer::T) where {T <: Layer}

Return true if the indices between the two layers are compatible with each other,
false otherwise. `layer` is below `next_layer`.
"""
function space_invar_interlayer end

"""
    ascending_fixedpoint(layer::Layer)

Return the operator that is the fixed point of the average ascending superoperator of this
layer, normalised to have norm 1.
"""
function ascending_fixedpoint end

"""
    scalingoperator_initialguess(l::Layer, irrep)

Return an initial guess to be used in the iterative eigensolver that solves for scaling
operators of the layer `l`.
"""
function scalingoperator_initialguess end

"""
    ascend(op, layer::Layer)
Ascend a local operator `op` from the bottom of `layer` to the top.
"""
function ascend end

"""
    descend(rho, layer::Layer)

Descend a local density matrix `rho` from the top of `layer` to the bottom.
"""
function descend end

"""
    environment(op, layer::Layer, rho; vary_disentanglers = true)

Compute the environments with respect to `op` of all the tensors in the layer, and return
them as a `Layer`. `rho` is the local density matrix at the top indices of this layer.

If `vary_disentanglers = false`, only compute the environments for the isometries, and set
the environments for the disentanglers to zero.
"""
function environment end

"""
    minimize_expectation_ev(layer::Layer, env::Layer; vary_disentanglers = true)

Return a new layer that minimizes the expectation value with respect to the environment
`env`, using an Evenbly-Vidal update.

If `vary_disentanglers = false`, only update the isometries.
"""
function minimize_expectation_ev end


text:
"""
    Grids

A simple 1-dimensional uniform grid
for finite difference method.
"""
module Grids

export Grid, over_elems, over_elems_real, over_elems_ghost
export first_elem_above_surface, get_z

"""
    Grid{T}

A simple 1-dimensional uniform grid of
type `T` for finite difference method.
"""
struct Grid{T}
  """
  Number of physical elements
  """
  n_elem :: Int

  """
  Element size
  """
  dz :: T

  """
  Inverse of element size squared
  """
  dzi :: T

  """
  Inverse of element size squared
  """
  dzi2 :: T

  """
  z-coordinate
  """
  z :: Vector{T}

  """
  Number of ghost points per side
  """
  n_ghost :: Int

  """
  z-coordinate at the cell center at the top of the domain
  """
  zc_top :: T

  """
  z-coordinate at the top domain boundary
  """
  zn_top :: T

  """
  z-coordinate at the cell center at the bottom of the domain
  """
  zc_surf :: T

  """
  z-coordinate at the bottom domain boundary
  """
  zn_surf :: T
end

"""
    Grid(z_min::T, z_max::T, n_elem_real::Int, n_ghost::Int = 1)

A simple grid implementation that accepts the domain interval
`z_min` and `z_max`, the number of elements `n_elem_real` and
the number of ghost points per side `n_ghost`.
"""
function Grid(z_min::T, z_max::T, n_elem_real::Int, n_ghost::Int = 1) where T
  n_elem = n_elem_real+2*n_ghost
  dz = (z_max-z_min)/n_elem_real
  dzi = inv(dz)
  dzi2 = dzi^2
  z = [(i-n_ghost)*dz-dz/2 for i in 1:n_elem]
  zc_surf = z[1+n_ghost]
  zn_surf = (z[n_ghost]+z[1+n_ghost])/2
  zc_top = z[end-n_ghost]
  zn_top = (z[end-n_ghost]+z[end-n_ghost+1])/2
  return Grid(n_elem, dz, dzi, dzi2, z, n_ghost, zc_top, zn_top, zc_surf, zn_surf)
end

function Base.show(io::IO, grid::Grid)
  println(io, "----------------- Grid")
  println(io, "z       = ",grid.z)
  println(io, "n_elem  = ",grid.n_elem)
  println(io, "n_ghost = ",grid.n_ghost)
  println(io, "zc_top  = ",grid.zc_top)
  println(io, "zn_top  = ",grid.zn_top)
  println(io, "zc_surf = ",grid.zc_surf)
  println(io, "zn_surf = ",grid.zn_surf)
  println(io, "-----------------")
end

"""
    get_z(grid::Grid, k::Int)

Get the z-coordinate given element index
"""
@inline get_z(grid::Grid, k::Int) = grid.z[k]

"""
    over_elems(grid::Grid)

Get the range of indexes to traverse real and ghost grid elements
"""
@inline over_elems(grid::Grid) = 1:grid.n_elem

"""
    over_elems_real(grid::Grid)

Get the range of indexes to traverse only real grid elements
"""
@inline over_elems_real(grid::Grid) = 1+grid.n_ghost : grid.n_elem-grid.n_ghost

"""
    over_elems_ghost(grid::Grid)

Get the range of indexes to traverse only ghost grid elements
"""
@inline over_elems_ghost(grid::Grid) = setdiff(over_elems(grid), over_elems_real(grid))

"""
    first_elem_above_surface(grid::Grid)

Get the first element index above the surface
"""
@inline first_elem_above_surface(grid::Grid) = 1+grid.n_ghost

end

max_stars_repo_name: christophernhill-dev/CLIMA

max_issues_repo_name: christophernhill-dev/CLIMA

max_forks_repo_name: christophernhill-dev/CLIMA


text:
# most machines will be higher resolution than this, but we're playing it safe
const RESOLUTION = 1000 # 1 μs = 1000 ns

##############
# Parameters #
##############

type Parameters
    seconds::Float64
    samples::Int
    evals::Int
    overhead::Int
    gctrial::Bool
    gcsample::Bool
    time_tolerance::Float64
    memory_tolerance::Float64
end

const DEFAULT_PARAMETERS = Parameters(5.0, 10000, 1, 0, true, false, 0.05, 0.01)

function Parameters(; seconds = DEFAULT_PARAMETERS.seconds,
                    samples = DEFAULT_PARAMETERS.samples,
                    evals = DEFAULT_PARAMETERS.evals,
                    overhead = DEFAULT_PARAMETERS.overhead,
                    gctrial = DEFAULT_PARAMETERS.gctrial,
                    gcsample = DEFAULT_PARAMETERS.gcsample,
                    time_tolerance = DEFAULT_PARAMETERS.time_tolerance,
                    memory_tolerance = DEFAULT_PARAMETERS.memory_tolerance)
    return Parameters(seconds, samples, evals, overhead, gctrial,
                      gcsample, time_tolerance, memory_tolerance)
end

function Parameters(default::Parameters; seconds = nothing, samples = nothing,
                    evals = nothing, overhead = nothing, gctrial = nothing,
                    gcsample = nothing, time_tolerance = nothing,
                    memory_tolerance = nothing)
    params = Parameters()
    params.seconds = seconds != nothing ? seconds : default.seconds
    params.samples = samples != nothing ? samples : default.samples
    params.evals = evals != nothing ? evals : default.evals
    params.overhead = overhead != nothing ? overhead : default.overhead
    params.gctrial = gctrial != nothing ? gctrial : default.gctrial
    params.gcsample = gcsample != nothing ? gcsample : default.gcsample
    params.time_tolerance = time_tolerance != nothing ? time_tolerance : default.time_tolerance
    params.memory_tolerance = memory_tolerance != nothing ? memory_tolerance : default.memory_tolerance
    return params::BenchmarkTools.Parameters
end

@compat function Base.:(==)(a::Parameters, b::Parameters)
    return a.seconds == b.seconds &&
           a.samples == b.samples &&
           a.evals == b.evals &&
           a.overhead == b.overhead &&
           a.gctrial == b.gctrial &&
           a.gcsample == b.gcsample &&
           a.time_tolerance == b.time_tolerance &&
           a.memory_tolerance == b.memory_tolerance
end

Base.copy(p::Parameters) = Parameters(p.seconds, p.samples, p.evals, p.overhead, p.gctrial,
                                      p.gcsample, p.time_tolerance, p.memory_tolerance)

function loadparams!(a::Parameters, b::Parameters, fields...)
    fields = isempty(fields) ? fieldnames(Parameters) : fields
    for f in fields
        setfield!(a, f, getfield(b, f))
    end
    return a
end

################################
# RESOLUTION/OVERHEAD settings #
################################

@noinline nullfunc() = nothing

@noinline function overhead_sample(evals)
    start_time = time_ns()
    for _ in 1:evals
        nullfunc()
    end
    sample_time = time_ns() - start_time
    return Int(cld(sample_time, evals))
end

function estimate_overhead()
    x = typemax(Int)
    for _ in 1:10000
        y = overhead_sample(RESOLUTION)
        if y < x
            x = y
        end
    end
    return x
end


text:
using GtkInteract
using Plots
backend(:immerse)

## Manipulate can display text output in a label
@manipulate for n=1:10, x=1:10
   (n,x)
end

## Similary, a plot can be made. Here we use Immerse through Plots:
@manipulate for n=1:10
    plot(sin, 0, n*pi)
end

## Unicode works, as does passing along functions through Options.
@manipulate for ϕ = 0:π/16:4π, f = Dict(:sin => sin, :cos => cos)
   plot(x -> f(x + ϕ), 0, pi)
end

## Can use plot! for adding layers
@manipulate for n=1:10
    plot(sin, 0, n*pi)
    plot!(cos, 0, n*pi, color=:red)
end

## and we can combine (mimicing an example from Interact)
@manipulate for ϕ = 0:π/16:4π, f = Dict(:sin => sin, :cos => cos), both = false
    if both
        plot(θ -> sin(θ + ϕ), 0, 8)
	plot!(θ -> cos(θ + ϕ), 0, 8, color=:red)
    else
       plot(θ -> f(θ + ϕ), 0, 8)
    end
end

## immerse works as directly, as well. In this case, the immerse toolbar is available:

using Immerse
@manipulate for n=1:10
    Immerse.plot(x=rand(n*10), y=rand(n*10))
end

# Text output
## The basic text output uses a label.
## This allows for some PANGO markup
## (https://developer.gnome.org/pango/stable/PangoMarkupFormat.html)
@manipulate for n = 1:20
    x = n > 10 ? "$n" : string(n)
    "The value is $x"
end

## Sadly, no latex output widget. So this example fails to render nicely
using SymPy
x = Sym("x")
@manipulate for n=1:10
    a = diff(sin(x^2), x, n)
    a
    jprint(a)
end

# Output widgets

## There are controls for gathering input, and output widgets for display. These consist of
## cairographic, textarea, label, and progress
## These values are `push!`ed onto within the function call. The function call should return `nothing`,
## else the output will be displayed as well in either a label or plot container.

## The label can be replaced by a multiline text buffer. The syntax is
## a bit awkward.  output widgets have values `push!`ed onto them and a final
## value of `nothing` is used.
@manipulate for n=1:10, x=1:10, out=textarea()
    push!(out, (n,x))
    nothing
end

## The progres bar is another output widget
@manipulate for n=1:100, pb=progress()
    push!(pb, n)
    nothing
end

## We can have more than one output widget
@manipulate for n=1:10, cg1=cairographic(width=300, height=200), cg2=cairographic(width=300, height=200)
    push!(cg1, Plots.plot(cos, 0, n*pi))
    push!(cg2, Plots.plot(sin, 0, n*pi))
    nothing
end

## Control widgets can be used instead of their being derived from the argument. This
## gives some flexibility and allows the labels to be different than the variable
## name associated with the control.
@manipulate for rb=radiobuttons([1,2,3],label="rb"), cb=checkbox(true, label="checkbox")
    (rb, cb)
end

## An example of @vchuravy from https://github.com/JuliaLang/Interact.jl/issues/36
using Distributions
@manipulate for α in 1:100, β = 1:100
    Plots.plot(x -> pdf(Beta(α, β), x), 0, 1)
end

## An example by @stevengj (https://github.com/JuliaLang/Interact.jl/issues/36)

# n steps of Newton iteration for sqrt(a), starting at x
function newton(a, x, n)
    for i = 1:n
        x = 0.5 * (x + a/x)
    end
    return x
end

# output x as HTML, with digits matching x0 printed in bold
function matchdigits(x::Number, x0::Number)
    s = string(x)
    s0 = string(x0)
    buf = IOBuffer()
    matches = true
    i = 0
    print(buf, "")           # pango  matches HTML
    while (i += 1) <= length(s)
        i % 30 == 0 && print(buf, "\n") # not

        if matches && i <= length(s0) && isdigit(s[i])
            if s[i] == s0[i]
                print(buf, s[i])
                continue
            end
            print(buf, "")
            matches = false
        end
        print(buf, s[i])
    end
    matches && print(buf, "")
    takebuf_string(buf)
end

set_bigfloat_precision(1024)
sqrt2 = sqrt(big(2))

@manipulate for l=label("Number of steps"), n = slider(0:9, value=0, label="n")
   matchdigits(newton(big(2), 2, n), sqrt2)
end


text:
using Plots, LaTeXStrings; pyplot()

a, c, d = 2, 1, 5
next(x,y) = [a*x*(1-x) - x*y, -c*y + d*x*y]
equibPoint = [(1+c)/d ,(d*(a-1)-a*(1+c))/d]

initX = [0.8,0.05]
tEnd = 100

traj = [[] for _ in 1:tEnd]
traj[1] = initX

for t in 2:tEnd
    traj[t] = next(traj[t-1]...)
end

scatter([traj[1][1]], [traj[1][2]],
	c=:black, ms=10,
	label="Initial state")
plot!(first.(traj),last.(traj),
	c=:blue, ls=:dash, m=(:dot, 5, Plots.stroke(0)),
	label="Model trajectory")
scatter!([equibPoint[1]], [equibPoint[2]],
	c=:red, shape=:cross, ms=10, label="Equlibrium point",
	xlabel=L"X_1", ylabel=L"X_2")

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
neighbors_bags(g) = ScatteredBags([setdiff(inneighbors(g, i), i) for i in vertices(g)])

graph2bag(verticesnode) = BagNode(verticesnode, [1:nobs(verticesnode)])

function joinve(xv, xe, g)
    @assert size(xe,2) == ne(g)

    bags = [Vector{Int}() for i in 1:nv(g)]
    foreach(e -> push!(bags[e.dst], e.src), edges(g))
    o = vcat(matchvertex2edges(xv, edges(g)), xe)
    BagNode(ArrayNode(o), ScatteredBags(bags))
end

function vertices2bag(mg, vertex2array)
    vertex_properties = reduce(cat, [ArrayNode(vertex2array(mg, i)) for i in vertices(mg)])
    TreeNode((vertex_properties, BagNode(vertex_properties,  neighbors_bags(mg))))
end

function vertices2bag(mg, vertex2array, edge2array)
    vertex_properties = reduce(hcat, [vertex2array(mg, v) for v in vertices(mg)])
    edge_properties = reduce(hcat, [edge2array(mg, e) for e in edges(mg)])
    TreeNode(ArrayNode(vertex_properties), joinve(vertex_properties, edge_properties, mg))
end

get_vertex_labels(mg, get_label) = [get_label(mg, v) for v in vertices(mg)]


text:
#
# This file is part of the Simulate.jl Julia package, MIT license
#
# Paul Bayer, 2020
#
# This is a Julia package for discrete event simulation
#

"""
    RTC::Vector{ClockChannel}

Real time clocks are registered to RTC and can be accessed and operated
through it.
"""
const RTC = Vector{ClockChannel}()


text:
# imports
using AccretionFormulae
using Plots
using Measures

# initialising plot backend
multiplier = 3
font_size = multiplier*10
line_width = multiplier*2
margin = multiplier*10mm
gr(
    legendfontsize=font_size,
    margin=margin,
    right_margin=margin*0.5,
    framestyle=:box,
    grid=false,
    lw=line_width,
    xlabel="R (M)",
    xguidefontsize=font_size,
    ylabel="Temperature (10⁶ K)",
    yguidefontsize=font_size,
    tickfontsize=font_size,
    thickness_scaling=1
    )

"""
Automatically determines the labels for various lines to go into the legend
"""
function label(value, name, units)
    scale = floor(log(10, value)) - 1
    if scale > 4
        scale = Int32(scale)
        value = "\$10^{$scale}\$"
    else
        value = "\$$value\$"
    end
    label = "$name = $value $(units) "
end

# constants
M_☼ = 1.99e30
G = 6.67e-11
c = 3e8
M = 10*M_☼
a_star = 0.0

# masses and spins to be changed
masses = [10, 100, 200, 10e6]
spins = [0.0, 0.5, 0.8, 0.95, 0.998]

# legend labels, colours and styles
mass_labels = permutedims(label.(masses, "M", "\$\\mathrm{M}_{\\odot}\$"))
spin_labels = permutedims(label.(spins, "Spin", ""))
colors = [:black :blue :red :green :purple]
styles = [:solid :dash :dot :dashdot :dashdotdot]

# changing the mass
temperatures_mass = []
for mass in masses
    mass *= M_☼
    R_isco = AccretionFormulae.r_isco(a_star, mass)
    r_vals = LinRange(R_isco, 50*mass, 1000)
    temperature_vals = AccretionFormulae.temperature.(r_vals, a_star, mass)
    @show(mass/M_☼)
    @show(maximum(temperature_vals))
    temperature_vals ./= 1e6
    push!(temperatures_mass, temperature_vals)
end

# changing the radius units to gravitational radii
R_isco_mass = AccretionFormulae.r_isco(a_star, M)
r_vals_mass = LinRange(R_isco_mass, 50*M, 1000)
r_vals_mass /= M

# changing the spin
temperatures_spin = []
r_vals_spins = []
for spin in spins
    R_isco = AccretionFormulae.r_isco(spin, M)
    r_vals = LinRange(R_isco, 30*M, 1000)
    temperature_vals = AccretionFormulae.temperature.(r_vals, spin, M)
    @show(spin)
    @show(maximum(temperature_vals))
    temperature_vals ./= 1e6
    # every line has a ISCO, so need to create r_vals in gravitational radii individually
    r_vals = r_vals ./ M
    push!(r_vals_spins, r_vals)
    push!(temperatures_spin, temperature_vals)
end

# plots
changing_mass_plot = plot(
    r_vals_mass,
    temperatures_mass,
    label=mass_labels,
    linecolor = :black,
    linestyle = styles
    )

changing_spin_plot = plot(
    r_vals_spins,
    temperatures_spin,
    label=spin_labels,
    linecolor = :black,
    linestyle = styles
    )

combined_plot_final = plot(changing_mass_plot, changing_spin_plot, layout=(1,2), size=multiplier.*(1200,400))

# saving image
png(combined_plot_final, "emission_profiles.png")
display(combined_plot_final)


text:
@testset "py2jl mod gen" begin
    py2jl"""
def sumBy(f, seq):
    s = 0
    for each in seq:
        s = s + f(each)
    return s

result = sumBy(lambda x: 1 + x, [100, 200])
    """
    @testset "transpiled codes run in the pure Julia runtime[1]" begin
        @test result == 302
        @test sumBy(identity, [1, 2, 3]) === 6
    end

    py2jl"""
def sum_by(f, seq):
    s = 0
    for e in seq:
        s = s + f(e)
    return s
result = sum_by(lambda x: x * 10, [1, 2, 3])
    """

     @testset "transpiled codes run in the pure Julia runtime[2]" begin
        @test result == 60
        @test sum_by(identity, [1, 2, 3]) === 6
    end

    py2jl"""
def __BOND_repl():
    SENTINEL = 1
    while SENTINEL>1:
        SENTINEL = 2
    return SENTINEL
"""
    @testset "while & compare" begin
        @test __BOND_repl() == 1
    end

    py2jl"""
    def __BOND_repl():
        SENTINEL = 10
        while True:
            break
        return SENTINEL
    """
    @testset "True" begin
            @test __BOND_repl() == 10
    end
end


text:
export compose_basic, compose_epsilon

function compose_basic(a::Wfst, b::Wfst)
    # Assert that the WFSTs are using the same semiring.
    @assert a.semiring == b.semiring

    # The input and output alphabets are those of a and b respectively
    @assert a.output_alphabet == b.input_alphabet
    input_alphabet = a.input_alphabet
    output_alphabet = b.output_alphabet
    debug("a.output_alphabet == b.input_alphabet")

    # Initial states are the cross product of the states of the two FSTs
    initial_states =
            Set([(x,y) for x in a.initial_states, y in b.initial_states])
    debug(string("initial_states: ", initial_states))

    # Initialize the weights for the initial states
    initial_weights = Dict()
    for state in initial_states
        if haskey(a.initial_weights, state[1]) &&
                haskey(b.initial_weights, state[2])
            initial_weights[state] = a.semiring.multiplication(
                    a.initial_weights[state[1]], b.initial_weights[state[2]])
        end
    end
    debug(string("initial_weights: ", initial_weights))
    # States start with the initial states
    states = Set([(x,y) for x in a.initial_states, y in b.initial_states])
    # Our queue also starts with the initial states
    queue = [(x,y) for x in a.initial_states, y in b.initial_states]
    queue = reshape(queue, length(queue))
    # Let the final states start empty
    final_states = Set()
    final_weights = Dict()
    arcs::Set{Arc} = Set{Arc}()

    # From the starting states, propagate through the Fst along possible
    # composed paths adding states and weights as we go.
    while queue != []
        debug(string("queue: ", queue))
        state = shift!(queue)
        if state in [(x,y) for x in a.final_states, y in b.final_states]
            debug(string("state ", state, " in final state cross product."))
            final_states = union(final_states, Set([state]))
            final_weights[state] = a.semiring.multiplication(
                    a.final_weights[state[1]], b.final_weights[state[2]])
        end
        # For each arc combination that travels from the states of 'a' and 'b'
        for arc_combo in [(x,y) for x in a.arcs, y in b.arcs]
            if arc_combo[1].from == state[1] && arc_combo[2].from == state[2]
                # If the output of the 'a' arc == the input of the 'b' arc
                if arc_combo[1].output == arc_combo[2].input
                    next_state = (arc_combo[1].to, arc_combo[2].to)
                    if !(next_state in states)
                        states = union(states, Set([next_state]))
                        push!(queue, next_state)
                    end
                    product = a.semiring.multiplication(
                            arc_combo[1].weight, arc_combo[2].weight)
                    new_arc = Arc(state,next_state,
                            arc_combo[1].input, arc_combo[2].output,
                            product)
                    debug(string("new_arc: ", new_arc))
                    arcs = union(arcs, Set{Arc}([new_arc]))
                end
            end
        end
    end

    return Wfst(states, input_alphabet, output_alphabet, initial_states,
           final_states, arcs, initial_weights, final_weights, a.semiring)
end

function compose_epsilon(a::Wfst, b::Wfst)
    @assert a.output_alphabet == b.input_alphabet

    for arc in a.arcs
        if arc.output == "<eps>"
            arc.output = "<eps2>"
        end
    end
    for arc in b.arcs
        if arc.input == "<eps>"
            arc.input = "<eps1>"
        end
    end
    for state in a.states
        add_arc!(a, state, state, "<eps>", "<eps1>", 1.0)
    end
    for state in b.states
        add_arc!(b, state, state, "<eps2>", "<eps>", 1.0)
    end
    f = create_filter(a.output_alphabet)
    c = compose_basic(compose_basic(a,f), b)
    return c
end

# Creates a filter that is used to filter out redundant paths when composing
# WFSTs with epsilon transitions.
function create_filter(alphabet::Set{String})
    f = Wfst()
    for symbol in alphabet
        add_arc!(f, 0, 0, symbol, symbol, 1.0)
    end
    add_arc!(f, 0, 0, "<eps2>", "<eps1>", 1.0)
    add_arc!(f, 0, 1, "<eps1>", "<eps1>", 1.0)
    add_arc!(f, 1, 1, "<eps1>", "<eps1>", 1.0)
    for symbol in alphabet
        add_arc!(f, 1, 0, symbol, symbol, 1.0)
    end
    add_arc!(f, 0, 2, "<eps2>", "<eps2>", 1.0)
    add_arc!(f, 2, 2, "<eps2>", "<eps2>", 1.0)
    for symbol in alphabet
        add_arc!(f, 2, 0, symbol, symbol, 1.0)
    end
    add_initial_state!(f, 0, 1.0)
    add_final_state!(f, 0, 1.0)
    add_final_state!(f, 1, 1.0)
    add_final_state!(f, 2, 1.0)
    return f
end


text:
begin

	# Instead of maps we now reduce the maps to only the important bits
	# This reduces the memory/storage used dramatically.
	# Also calculation of binary-results change dramatically.
	# Output would be like:
	# 11110101000000000000000100000000000000000001
	# By counting the digits you simply find out by what number it is divideable.
	# In this case 24 is divideable through 1,2,3,4,6,8,12 and 24.
	# ToDo: Update main.cpp (use big number library there)
	#       I do not prefer to go back to lists or anything like that, maybe need to wait for updates.

	function normalize_num(v,KL)
		R = 0
		i = 0
		n = 0
		x = 0
		y = 0

		while true
			if v & BigInt(2)^i != 0
				R |= BigInt(2)^x
			end
			x += 1
			i += 1
			n += 1
			if (KL & BigInt(2)^n) != 0
				if x == 1
					R |= 1
					break
				end
				#println(x," ",i)
				x = 0
				n += 1
			end
		end
		R
	end
	function calc_fact_num(pos,fnum)
		R = BigInt(0)
		M = C[pos]
		N = K[fnum]
		i = 0
		for x in 1:SIZE
			for y in 1:SIZE
				if N[x,y]
					if M[x,y]
						R |= BigInt(2)^i
					end
					i += 1
				end
			end
		end
		R
	end
	function calc_kl(KL_X)
		R = BigInt(0)
		i = 0
		done = false
		for y in 1:SIZE
			for x in 1:SIZE
				if !KL_X[x,y]
					if x == 1
						done = true
					end
					R |= BigInt(2)^i
					if done
						break
					end
					i += 1
					break
				end
				i += 1
			end
			if done
				break
			end
		end
		R
	end
	function not(v)
		v $ 0xFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF
	end

	startT = Dates.unix2datetime(time())
	SIZE=2^5
	A = reshape(collect(0:SIZE*SIZE-1) % SIZE + 1,SIZE,SIZE)
	function F(x,y)
			A[y,x]
	end
	B = [F(x,y) for x=collect(1:SIZE),y=collect(1:SIZE)]

	C = Array(Array,6)
	C[1] = ((A.*B)&(2^0)).==2^0
	C[2] = ((A.*B)&(2^1)).==2^1
	C[3] = ((A.*B)&(2^2)).==2^2
	C[4] = ((A.*B)&(2^3)).==2^3
	C[5] = ((A.*B)&(2^4)).==2^4
	C[6] = ((A.*B)&(2^5)).==2^5

	K = Array(Array,6)
	K[1] = (A.*B).<=2^1
	K[2] = (A.*B).<=2^2
	K[3] = (A.*B).<=2^3
	K[4] = (A.*B).<=2^4
	K[5] = (A.*B).<=2^5
	K[6] = (A.*B).<=2^6

	#We could actually factorize with just those numbers instead of the big maps
	J1_1 = calc_fact_num(1,1)
	J1_2 = calc_fact_num(1,2)
	J1_3 = calc_fact_num(1,3)
	J1_4 = calc_fact_num(1,4)
	J1_5 = calc_fact_num(1,5)
	J1_6 = calc_fact_num(1,6)

	J2_1 = calc_fact_num(2,1)
	J2_2 = calc_fact_num(2,2)
	J2_3 = calc_fact_num(2,3)
	J2_4 = calc_fact_num(2,4)
	J2_5 = calc_fact_num(2,5)
	J2_6 = calc_fact_num(2,6)

	J3_1 = calc_fact_num(3,1)
	J3_2 = calc_fact_num(3,2)
	J3_3 = calc_fact_num(3,3)
	J3_4 = calc_fact_num(3,4)
	J3_5 = calc_fact_num(3,5)
	J3_6 = calc_fact_num(3,6)

	J4_1 = calc_fact_num(4,1)
	J4_2 = calc_fact_num(4,2)
	J4_3 = calc_fact_num(4,3)
	J4_4 = calc_fact_num(4,4)
	J4_5 = calc_fact_num(4,5)
	J4_6 = calc_fact_num(4,6)

	J5_1 = calc_fact_num(5,1)
	J5_2 = calc_fact_num(5,2)
	J5_3 = calc_fact_num(5,3)
	J5_4 = calc_fact_num(5,4)
	J5_5 = calc_fact_num(5,5)
	J5_6 = calc_fact_num(5,6)

	KL_3 = calc_kl(K[3])
	KL_4 = calc_kl(K[4])
	KL_5 = calc_kl(K[5])

	endT = Dates.unix2datetime(time())

	# Numbers would need to be generated only ONCE(then stored somewhere) and loaded when needed, wich should not take that long
	println("Time used for generating basic numbers: ",endT-startT)

	# binary 111 -> 7
	v = J1_3 & J2_3 & J3_3
	v = normalize_num(v,KL_3)
	println(reverse(bin(v)))

	println()
	# binary 1001 -> 9
	v = J1_4 & not(J2_4) & not(J3_4) & J4_4
	v = normalize_num(v,KL_4)
	println(reverse(bin(v)))

	println()

	println("J1_5:",J1_5)
	println("J2_5:",J2_5)
	println("J3_5:",J3_5)
	println("J4_5:",J4_5)
	println("J5_5:",J5_5)
	println("KL_5:",KL_5)

	println()

	startT = Dates.unix2datetime(time())
	# binary 11000 -> 24
	v = not(J1_5) & not(J2_5) & not(J3_5) & J4_5 & J5_5
	v = normalize_num(v,KL_5)
	println(reverse(bin(v)))
	endT = Dates.unix2datetime(time())
	println("Time actually needed to calculate from precalculated numbers: ",endT-startT)

end

max_stars_repo_name: TheCrazyT/factorisation-maps

max_issues_repo_name: TheCrazyT/factorisation-maps

max_forks_repo_name: TheCrazyT/factorization-maps


text:
# Portions translated from SLICOT-Reference distribution
# Copyright (c) 2002-2020 NICONET e.V.
function run_mb02dd(datfile, io=stdout)
    NIN = 5
    NOUT = 6
    KMAX = 20
    MMAX = 20
    NMAX = 20
    LDG = KMAX*( MMAX + NMAX )
    LDL = KMAX*( MMAX + NMAX )
    LDR = KMAX*( MMAX + NMAX )
    LDT = KMAX*( MMAX + NMAX )
    LDWORK = ( MMAX + NMAX - 1 )*KMAX
    LCS = 3*LDWORK
    tdim = KMAX*(MMAX+NMAX)
    T = Array{Float64,2}(undef, LDT,KMAX*(MMAX+NMAX))
    G = Array{Float64,2}(undef, LDG,KMAX*(MMAX+NMAX))
    R = Array{Float64,2}(undef, LDR,KMAX*(MMAX+NMAX))
    L = Array{Float64,2}(undef, LDL,KMAX*(MMAX+NMAX))
    CS = Array{Float64,1}(undef, LCS)
    DWORK = Array{Float64,1}(undef, LDWORK)
f = open(datfile,"r")
    readline(f)
    vs = split(readline(f))
    N = parse(BlasInt, vs[1])
    K = parse(BlasInt, vs[2])
    M = parse(BlasInt, vs[3])
    JOB =  vs[4][1]
    TYPET =  vs[5][1]
    S = ( N + M )*K
if ( N<=0 || N>NMAX )
else
if ( K<=0 || K>KMAX )
else
if ( M<=0 || M>MMAX )
else
if ( LSAME( TYPET, 'R' ) )

    vs = String[]
    _isz,_jsz = (K,S)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       T[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end
else

    vs = String[]
    _isz,_jsz = (S,K)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       T[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end
end # if
# interp call 1

    INFO = SLICOT.mb02cd!(JOB, TYPET, K, N, T, G, R, L, CS, LCS, LDWORK)
    @test INFO == 0
    INFO == 0 || return

if ( INFO!=0 )
else
# interp output 1
    println(io, "R:")
    _nc = N*K
    _nr = N*K
    show(io, "text/plain", R[1:_nr,1:_nc])
    println(io)

if ( LSAME( JOB, 'R' ) || LSAME( JOB, 'A' ) )
if ( LSAME( TYPET, 'R' ) )
# interp output 2
    println(io, "G:")
    _nc = N*K
    _nr = 2*K
    show(io, "text/plain", G[1:_nr,1:_nc])
    println(io)

else
# interp output 3
    println(io, "G:")
    _nc = 2*K
    _nr = N*K
    show(io, "text/plain", G[1:_nr,1:_nc])
    println(io)

end # if
end # if
if ( LSAME( JOB, 'A' ) )
# interp output 4
    println(io, "L:")
    _nc = N*K
    _nr = N*K
    show(io, "text/plain", L[1:_nr,1:_nc])
    println(io)

end # if
if ( LSAME( TYPET, 'R' ) )
# interp call 2

    # FOREIGN.dlacpy!( 'All', N*K, K, R(1,(N-1)*K+1), LDR, R(K+1,N*K+1), LDR )
    for ir in 1:N*K
        for ic in 1:K
            R[K+ir,N*K+ic] = R[ir,(N-1)*K+ic]
        end
    end

# interp call 3

    INFO = SLICOT.mb02dd!(JOB, TYPET, K, M, N, view(T,1:K,N*K+1:tdim), T, G,
                          view(R,:,N*K+1:(N*M)*K), view(L,N*K+1:LDL,:), CS, LDWORK)
    @test INFO == 0
    INFO == 0 || return

else
# interp call 4

    #FOREIGN.dlacpy!( 'All', K, N*K, R((N-1)*K+1,1), LDR, R(N*K+1,K+1), LDR )
    for ir in 1:K
        for ic in 1:N*K
            R[N*K+ir,K+ic] = R[(N-1)*K+ir,ic]
        end
    end

# interp call 5

    INFO = SLICOT.mb02dd!(JOB, TYPET, K, M, N, view(T,N*K+1:(N+M)*K,:), T, G, view(R,N*K+1:(N+M)*K,:), view(L,:, N*K+1:(N+M)*K), CS, LDWORK)
    @test INFO == 0
    INFO == 0 || return

end # if
if ( INFO!=0 )
else
# interp output 5
    println(io, "R:")
    _nc = S
    _nr = S
    show(io, "text/plain", R[1:_nr,1:_nc])
    println(io)

if ( LSAME( JOB, 'R' ) || LSAME( JOB, 'A' ) )
if ( LSAME( TYPET, 'R' ) )
# interp output 6
    println(io, "G:")
    _nc = S
    _nr = 2*K
    show(io, "text/plain", G[1:_nr,1:_nc])
    println(io)

else
# interp output 7
    println(io, "G:")
    _nc = 2*K
    _nr = S
    show(io, "text/plain", G[1:_nr,1:_nc])
    println(io)

end # if
end # if
if ( LSAME( JOB, 'A' ) )
# interp output 8
    println(io, "L:")
    _nc = S
    _nr = S
    show(io, "text/plain", L[1:_nr,1:_nc])
    println(io)

end # if
end # if
end # if
end # if
end # if
end # if
    close(f)
end # run_mb02dd()


text:
@testset "IntegralForms" begin

  @testset "Discriminant" begin
    A = [i for i in -12:13]
    C = [1,1,0,0,1,1,0,0,1,1,0,0,0,1,0,0,1,1,0,0,1,1,0,0,1,1]
    @test isdiscriminant.(A) == C
  end

  @testset "Conductor" begin
    A = [5, 8, 12, 13, 17, 20, 21, 24, 28, 29, 32, 33, 37, 40, 41, 44, 45, 48,
         52, 53, 56, 57, 60, 61, 65, 68, 69, 72, 73, 76, 77, 80, 84, 85, 88,
         89, 92, 93, 96, 97, 101, 104, 105, 108, 109, 112, 113, 116, 117, 120,
         124, 125, 128, 129, 132, 133, 136, 137, 140, 141, 145, 148, 149, 152,
         153, 156, 157, 160, 161, 164, 165, 168, 172, 173, 176, 177, 180, 181,
         184, 185, 188, 189, 192, 193, 197, 200, 201, 204, 205, 208, 209, 212,
         213, 216, 217, 220, 221, 224, 228, 229, 232, 233, 236, 237, 240, 241,
         244, 245, 248, 249, 252, 253, 257, 260, 261, 264, 265, 268, 269, 272,
         273, 276, 277, 280, 281, 284, 285, 288, 292, 293, 296, 297, 300, 301,
         304, 305, 308, 309, 312, 313, 316, 317, 320, 321, 325, 328, 329, 332,
         333, 336, 337, 340, 341, 344, 345, 348, 349, 352, 353, 356, 357, 360,
         364, 365, 368, 369, 372, 373, 376, 377, 380, 381, 384, 385, 388, 389,
         392, 393, 396, 397, 401, 404, 405, 408, 409, 412, 413, 416, 417, 420,
         421, 424, 425, 428, 429, 432, 433, 436, 437, 440, 444, 445, 448, 449,
         452, 453, 456, 457, 460, 461, 464, 465, 468, 469, 472, 473, 476, 477,
         480, 481, 485, 488, 489, 492, 493, 496, 497, 500, 501, 504, 505, 508,
         509, 512, 513, 516, 517, 520, 521, 524, 525, 528, 532, 533, 536, 537,
         540, 541, 544, 545, 548, 549, 552, 553, 556, 557, 560, 561, 564, 565,
         568, 569, 572, 573, 577, 580, 581, 584, 585, 588, 589, 592, 593, 596,
         597, 600, 601, 604, 605, 608, 609, 612, 613, 616, 617, 620, 621, 624,
         628, 629, 632, 633, 636, 637, 640, 641, 644, 645, 648, 649, 652, 653,
         656, 657, 660, 661, 664, 665, 668, 669, 672, 673, 677, 680, 681, 684,
         685, 688, 689, 692, 693, 696, 697, 700, 701, 704, 705, 708, 709, 712,
         713, 716, 717, 720, 721, 724, 725, 728, 732, 733, 736, 737, 740, 741,
         744, 745, 748, 749, 752, 753, 756, 757, 760, 761, 764, 765, 768, 769,
         772, 773, 776, 777, 780, 781, 785, 788, 789, 792, 793, 796, 797, 800,
         801, 804, 805, 808, 809, 812, 813, 816, 817, 820, 821, 824, 825, 828,
         829, 832, 833, 836, 837, 840, 844, 845, 848, 849, 852, 853, 856, 857,
         860, 861, 864, 865, 868, 869, 872, 873, 876, 877, 880, 881, 884, 885,
         888, 889, 892, 893, 896, 897, 901, 904, 905, 908, 909, 912, 913, 916,
         917, 920, 921, 924, 925, 928, 929, 932, 933, 936, 937, 940, 941, 944,
         945, 948, 949, 952, 953, 956, 957, 960, 964, 965, 968, 969, 972, 973,
         976, 977, 980, 981, 984, 985, 988, 989, 992, 993, 996, 997, 1000]
    Acond = [1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 3, 2, 2, 1, 1, 1,
             1, 1, 1, 2, 1, 3, 1, 1, 1, 4, 2, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 3,
             1, 2, 1, 2, 3, 1, 1, 5, 4, 1, 2, 1, 1, 1, 1, 1, 1, 2, 1, 1, 3, 1,
             1, 2, 1, 2, 1, 1, 1, 1, 2, 1, 6, 1, 1, 1, 1, 3, 4, 1, 1, 5, 1, 1,
             1, 4, 1, 2, 1, 3, 1, 1, 1, 2, 2, 1, 1, 1, 1, 1, 2, 1, 2, 7, 1, 1,
             3, 1, 1, 2, 3, 1, 1, 1, 1, 4, 1, 2, 1, 1, 1, 1, 1, 6, 2, 1, 1, 3,
             5, 1, 2, 1, 2, 1, 1, 1, 1, 1, 8, 1, 5, 1, 1, 1, 3, 4, 1, 2, 1, 1,
             1, 1, 1, 2, 1, 2, 1, 3, 1, 1, 2, 3, 2, 1, 1, 1, 1, 1, 4, 1, 2, 1,
             7, 1, 3, 1, 1, 2, 9, 1, 1, 1, 1, 2, 1, 2, 1, 1, 5, 1, 1, 6, 1, 2,
             1, 1, 1, 1, 4, 1, 2, 1, 1, 1, 1, 1, 4, 1, 6, 1, 1, 1, 1, 3, 2, 1,
             1, 1, 1, 1, 1, 2, 1, 10, 1, 3, 1, 1, 1, 8, 3, 2, 1, 1, 1, 1, 5, 4,
             2, 1, 1, 1, 3, 1, 2, 1, 2, 3, 1, 1, 1, 1, 2, 1, 2, 1, 1, 1, 1, 1,
             1, 2, 1, 1, 3, 7, 1, 4, 1, 2, 1, 5, 1, 1, 11, 2, 1, 6, 1, 1, 1, 1,
             3, 2, 2, 1, 1, 1, 1, 7, 4, 1, 2, 1, 9, 1, 1, 1, 4, 3, 2, 1, 1, 1,
             1, 1, 2, 1, 1, 1, 1, 3, 1, 2, 1, 2, 3, 1, 1, 5, 1, 4, 1, 2, 1, 1,
             1, 1, 1, 12, 1, 2, 5, 1, 1, 1, 2, 1, 2, 1, 1, 1, 1, 1, 2, 1, 6, 1,
             1, 1, 1, 3, 8, 1, 2, 1, 1, 1, 1, 1, 1, 2, 1, 3, 1, 1, 1, 10, 3, 2,
             1, 1, 1, 1, 1, 2, 1, 2, 1, 1, 5, 3, 1, 8, 7, 2, 3, 1, 1, 13, 4, 1,
             2, 1, 1, 1, 1, 1, 6, 1, 2, 1, 1, 3, 1, 1, 2, 1, 2, 1, 1, 1, 1, 1,
             4, 1, 1, 1, 1, 1, 3, 4, 1, 2, 1, 1, 1, 1, 5, 2, 1, 2, 1, 3, 1, 1,
             1, 2, 3, 2, 1, 1, 1, 1, 1, 4, 2, 1, 11, 1, 9, 1, 4, 1, 14, 3, 1,
             1, 1, 1, 2, 1, 2, 1, 5]
    @test conductor.(A) == Acond
  end

  @testset "FundamentalDiscriminant" begin
    @test isfundamental_discriminant(12) == true
    @test isfundamental_discriminant(-12) == false
    D =	[1, 5, 8, 12, 13, 17, 21, 24, 28, 29, 33, 37, 40, 41, 44, 53, 56, 57,
         60, 61, 65, 69, 73, 76, 77, 85, 88, 89, 92, 93, 97, 101, 104, 105,
         109, 113, 120, 124, 129, 133, 136, 137]
    @test all(isfundamental_discriminant, D)
    D =	map(x -> -x, [3, 4, 7, 8, 11, 15, 19, 20, 23, 24, 31, 35, 39, 40, 43,
                      47, 51, 52, 55, 56, 59, 67, 68, 71, 79, 83, 84, 87, 88,
                      91, 95, 103])
    @test all(isfundamental_discriminant, D)

    @test fundamental_discriminant(20) == 5
    @test fundamental_discriminant(-20) == -20
  end

  @testset "DefiniteForms" begin
    @test isnegative_definite(QuadBin(-2, 3, -2)) == true
    @test ispositive_definite(QuadBin(2, -3, 2)) == true
    @test isindefinite(QuadBin(1, 3, 1)) == true
    @test ispositive_definite(QuadBin(1, 3, 1)) == false
    @test isnegative_definite(QuadBin(1, 3, 1)) == false
  end

  @testset "Reduction" begin
  @test isreduced(QuadBin(1, 2, 3)) == false
  @test isreduced(QuadBin(2, 1, 3)) == true
  @test isreduced(QuadBin(1, -1, 1)) == false
  @test isreduced(QuadBin(1, 1, 1)) == true
  @test isreduced(QuadBin(-1, 2, 2)) == true
  @test isreduced(QuadBin(1, 9, 4)) == false
  @test isreduced(QuadBin(1, 5, -1)) == true
  @test reduce(QuadBin(195751, 37615, 1807)) == QuadBin(1, 1, 1)
  @test reduce(QuadBin(33, 11, 5)) == QuadBin(5, -1, 27)
  @test reduce(QuadBin(15, 0, 15)) == QuadBin(15, 0, 15)
  end

  @testset "Composition" begin
  @test reduce(compose(QuadBin(2 ,2, 1), QuadBin(5, 4, 1))) == QuadBin(1, 0, 1)
  @test reduce(compose(QuadBin(1, 1, 6), QuadBin(1, 1, 6))) == QuadBin(1, 1, 6)
  @test reduce(compose(QuadBin(2, -1, 3), QuadBin(2, -1, 3))) == QuadBin(2, 1, 3)
  end

  @testset "Cycle" begin

    A1 = [QuadBin(1, 7, -6), QuadBin(6, 5, -2), QuadBin(2, 7, -3), QuadBin(3, 5, -4),
    QuadBin(4, 3, -4), QuadBin(4, 5, -3), QuadBin(3, 7, -2), QuadBin(2, 5, -6), QuadBin(6, 7, -1)]

    @test cycle(QuadBin(1,7,-6)) == A1

    A2 = [QuadBin(1, 8, -3), QuadBin(3, 4, -5), QuadBin(5, 6, -2), QuadBin(2, 6, -5),
    QuadBin(5, 4, -3), QuadBin(3, 8, -1)]

    @test cycle(QuadBin(1,8,-3)) == A2
  end

  @testset "CanSolve" begin
    @test can_solve(QuadBin(3,2,2), 28) == (true, (2, 2))
    @test can_solve(QuadBin(2,1,3), 3) == (true, (0, 1))
    @test can_solve(QuadBin(2,1,3), 5) == (false, (0, 0))
    @test can_solve(QuadBin(2,1,3), 6) == (true, (1, 1))
    @test can_solve(QuadBin(2,-1,3), 3) == (true, (0, 1))
    @test can_solve(QuadBin(2,-1,3), 5) == (false, (0, 0))
    @test can_solve(QuadBin(2,-1,3), 6) == (true, (-1, 1))
    @test can_solve(QuadBin(1,1,6), 3) == (false, (0, 0))
    @test can_solve(QuadBin(1,1,6), 5) == (false, (0, 0))
    @test can_solve(QuadBin(1,1,6), 6) == (true, (0, 1))
  end

  @testset "PrimeForm" begin
    f = prime_form(ZZ(17), ZZ(1237))
    @test discriminant(f) == 17 && f[1] == 1237
    f = prime_form(ZZ(12), ZZ(743))
    @test discriminant(f) == 12 && f[1] == 743
    f = prime_power_form(ZZ(117), ZZ(3), 2)
    @test discriminant(f) == 117 && f[1] == 9
  end
end


text:
module UserTypes
include("FixedTypes.jl")
import .FixedTypes: Param, Out, Signal, Compression, MyData, @define_struct
@define_struct Params VI VRADJ
@define_struct Aux aux1
@define_struct Comment MyComment PicturePath
@define_struct Outs VQmean VQstd
@define_struct Signals vq vd vro vwo
export Param, Out, Signal, Compression,
	Params, Aux, Comment, Outs, Signals, MyData,
	@define_struct
end

#data    = MyData(param, aux, comment, out, wfm)


text:
function fuck(male, female)
    println("$(male) fucks $(female) or $(female) fucks $(male)?")
end

fuck("alice", "bob")

kiss(male, female) = "$(male) pushes tongue down $(female)'s throat or $(female) pushes tongue down $(male)'s throat?"

println(kiss("alice", "bob"))

function hypotenuse(a, b)
    a = abs(a)
    b = abs(b)
    if a > b
        r = b/a
        return a*sqrt(1+r*r)
    end
    if b == 0
        return zero(a)
    end
    r = a/b
    return b*sqrt(1+r*r)
end

println(hypotenuse(3, 4))

using Dates
function Date(y::Int64, m::Int64=1, d::Int64=1)
    err = Dates.validargs(Dates.Date, y, m, d)
    err === nothing || throw(err)
    return Dates.Date(Dates.UTD(Dates.totaldays(y, m, d)))
end
println(Date(2003, 12, 26))
println(methods(Date))

max_stars_repo_name: JulialangOrgCN/howtojulia

max_issues_repo_name: JulialangOrgCN/howtojulia

max_forks_repo_name: JulialangOrgCN/howtojulia


text:
using SatelliteToolbox
using Plots
gr()

#=
設定パラメータ
=#
DataNum =100 #シミュレータ反復回数
dt = 60*5 ##シミュレータの計算間隔 [s]
t = DatetoJD(2019, 11, 14, 07, 30, 00)   #シミュレート開始時刻＠ユリウス日

#開始時のみ処理
tles = read_tle("ISS_TLE.txt") #TLE読み込み
orbp = init_orbit_propagator(Val{:sgp4}, tles[1]) #SGP4での軌道モデルの読み込み
eop_IAU2000A = get_iers_eop(:IAU2000A)  #IAU-2000Aでの地球モデルの読み込み
epoch_JD = DatetoJD(2000+tles[1].epoch_year, 1, 1, 00, 00, 00) + tles[1].epoch_day  #TLEの元期＠ユリウス日
o, r, v = propagate!(orbp, t - epoch_JD)

#ループ毎処理関数定義
function SatPlaceCal(n)
    #=  衛星位置の赤道面座標系での計算　=#
    o, r, v = step!(orbp, dt)
    #=  衛星位置の地心直交座標系への変換  =#
    xg = rECItoECEF(GCRF(), ITRF(), t + n*dt/(24*60*60), eop_IAU2000A) * r
    vg = rECItoECEF(GCRF(), ITRF(), t + n*dt/(24*60*60), eop_IAU2000A) * v
    return xg, vg
end

function plus(a,b)
    c = a+b
    return c
end

x_ecef_log = zeros(100, 3)    #衛星位置＠地心直交座標系
x_geod_log = zeros(100, 3)    #衛星位置＠測地座標系
v_ecef_log = zeros(100, 3)    #衛星速度＠地心直交座標系
x_eci_log = zeros(100, 3)
#=
#シミュレーションループ
for i=1:DataNum
    x_sat, v_sat = SatPlaceCal(i)
    x_ecef_log[i,:] = x_sat
    v_ecef_log[i,:] = v_sat
    println(x_ecef_log[i,:])
    x_geo = ECEFtoGeodetic([-1*x_sat[1]; x_sat[2]; x_sat[3]])
    x_geod_log[i,1]=x_geo[1]
    x_geod_log[i,2]=x_geo[2]
    x_geod_log[i,3]=x_geo[3]
    println(x_geod_log[i,:])
end
=#

#=
#=************************************************
地心直交座標系での衛星位置を3Dプロットしたい場合はコレ
************************************************=#
println(x_ecef_log)

plot(x_ecef_log[:,1], x_ecef_log[:, 2], x_ecef_log[:, 3])
xlabel!("x")
ylabel!("y")
#zlabel!("z")

savefig("test2.png")
=#

#=
anim = @animate for i=1:DataNum
    time=epoch_JD + (i-1)*dt/(60*60*24)
    x_sat, v_sat = SatPlaceCal(i)
    x_ecef_log[i,:] = x_sat
    v_ecef_log[i,:] = v_sat
    x_sat = rECEFtoECI(ITRF(),GCRF(), time, eop_IAU2000A)*x_sat
    x_eci_log[i,:] = x_sat
    println(x_eci_log[i,:])
    plot(x_eci_log[1:i, 1], x_eci_log[1:i, 2])
end

gif(anim, "test7_anim.gif", fps = 10)
=#

#=***********************************************
緯度経度に変換してGifでグラフを作成したい場合はコレ
***********************************************=#
anim = @animate for i=1:DataNum
    x_sat, v_sat = SatPlaceCal(i)
    x_ecef_log[i,:] = x_sat
    v_ecef_log[i,:] = v_sat
    #plot(x_ecef_log[1:i,1], x_ecef_log[1:i, 2], x_ecef_log[1:i, 3])

    x_geo = ECEFtoGeodetic([x_ecef_log[i,1]; x_ecef_log[i,2]; x_ecef_log[i,3]])
    x_geod_log[i,1]=x_geo[1]
    if x_geo[2]>0
        x_geod_log[i,2]=x_geo[2]
    else
        x_geod_log[i,2]=x_geo[2]+2π
    end
    x_geod_log[i,3]=x_geo[3]
    plot(rad2deg.(x_geod_log[1:i,2]),rad2deg.(x_geod_log[1:i,1]))

end

gif(anim, "test5_anim.gif", fps = 50)

max_stars_repo_name: yui-project/ACSim

max_issues_repo_name: yui-project/ACSim

max_forks_repo_name: yui-project/ACSim


text:
# Analytic
struct AnalyticAlg <: DistAlgorithm end
const Analytic = AnalyticAlg()

analytic(distop, x::RandVar) = distop(distributionsjl(x))

for dop in distops_names
  expr = :($dop(x::RandVar, ::AnalyticAlg) = analytic($dop, x))
  eval(expr)
end

# Sample mean
struct SampleMeanAlg <: DistAlgorithm end
const SampleMean = SampleMeanAlg()
mean(x, ::SampleMeanAlg; kwargs...) = samplemean(x; kwargs...)

"Expectation of `x` from `n` samples"
function samplemean(x::RandVar; n = 10_000, alg = Omega.Inference.RejectionSample, kwargs...)
    # @show kwargs
  mean(rand(x, n; alg = alg, kwargs...))
end

# Sample mean std
meanstd(x) = (m = mean(x); (mean = m, std = stdm(x, m)))

# Sample prob
sampleprob(x::RandVar, n; alg = FailUnsat, kwargs...) = samplemean(x, n; alg = alg, kwargs...)

# samplemeanᵣ(x, n; alg = RejectionSample) = ciid(ω -> samplemean(x(ω), n; alg = alg))
# samplemeanstd(x, n; alg, kwargs...) = meanstd(rand(x, n; alg = alg, kwargs...))
# samplemeanstdᵣ(x, n; alg = RejectionSample) = ciid(ω -> samplemeanstd(x(ω), n; alg = alg))
# const sampleprob = samplemean
# const sampleprobᵣ = samplemeanᵣ
# meanᵣ(x; alg, kwargs...) = samplemeanᵣ(x, 1000; alg = alg, kwargs...)
# zt: Default to RejectionSample is problematic

# # Sample probability

# "Probability x is true"
# prob(x::RandVar, n, israndvar::Type{Val{false}}) =
#   sum((rand(x, alg = RejectionSample) for i = 1:n)) / n
# lprob(x::RandVar, n = 1000) = lift(prob)(x, n, Val{false})
# prob(x::RandVar, n) = prob(x, n, Val{elemtype(x) <: RandVar})

# const bindistattrs = [:entropy, :mgf, :cf, :pdf, :logpdf, :loglikelihood, :cdf,
#                       :logcdf, :ccdf, :logccdf, :quantile, :cquantile,
#                       :invlogcdf, :invlogccdf]

# for func in bindistattrs
#   expr =
#   quote
#     $func(x::PrimRandVar, et::Type{<:Real}) = Djl.$func(distribution(x))
#     $func(x::RandVar, et::Type{<:RandVar}, t) = ciid(mean, x)
#     $func(x::RandVar, t) = $func(x, elemtype(x))
#   end
#   eval(expr)
# end
# Djl.entropy(x::RandVar) = entropy(distribution(x))
# Djl.entropy(x::RandVar, base) = entropy(distribution(x), base)

# Djl.mgf(x::RandVar, t) = mgf(distribution(x), t)
# Djl.cf(x::RandVar, t) = cf(distribution(x), t)

# Djl.pdf(x::RandVar, t) = pdf(distribution(x), t)
# Djl.logpdf(x::RandVar, t) = logpdf(distribution(x), t)
# Djl.loglikelihood(x::RandVar, t) = loglikelihood(distribution(x), t)

# Djl.cdf(x::RandVar, t) = cdf(distribution(x), t)
# Djl.logcdf(x::RandVar, t) = logcdf(distribution(x), t)

# Djl.ccdf(x::RandVar, t) = ccdf(distribution(x), t)
# Djl.logccdf(x::RandVar, t) = logccdf(distribution(x), t)

# Djl.quantile(x::RandVar, t) = quantile(distribution(x), t)
# Djl.cquantile(x::RandVar, t) = cquantile(distribution(x), t)

# Djl.invlogcdf(x::RandVar, t) = invlogcdf(distribution(x), t)
# Djl.invlogccdf(x::RandVar, t) = invlogccdf(distribution(x), t)

# Converions between Distributions and Omega