text:
module gf

using LinearAlgebra

include("DivTable8.jl")
include("MulTable8.jl")

export GF, rank, rref_with_pivots!

struct GF <: Number
    x::Int64
    function GF(x::Integer)
        if x > 255 || x < 0
            error("the field only supports 0-255")
        end
        return new(x)
    end
end

# unary plus
Base.:(+)(a::GF) = a
# unary minus
Base.:-(a::GF) = a
# binary plus
Base.:(+)(a::GF, b::GF)::GF = GF(a.x ⊻ b.x)
# binary minus
Base.:-(a::GF, b::GF)::GF = GF(a.x ⊻ b.x)

# abs
Base.:abs(a::GF)::Integer = a == GF(0) ? 0 : 1

# # FIXME: this is a trivial norm, perhaps should implement a general p-norm?
LinearAlgebra.:norm(a::GF) = a == GF(0) ? 0 : 1

function LinearAlgebra.:rank(A::AbstractArray{GF,2})
    isempty(A) && return 0 # 0-dimensional case
    # s = svdvals(A)
    # tol = max(atol, rtol*s[1])
    # count(x -> x > tol, s)
    B, p = rref_with_pivots!(A)
    return length(p)
end

# # binary multiplication
Base.:*(a::GF, b::GF)::GF = GF(mulTable8[a.x << 8 | b.x + 1])
Base.:/(a::GF, b::GF)::GF = GF(divTable8[a.x << 8 | b.x + 1])

# # modified from https://github.com/blegat/RowEchelon.jl/blob/master/src/RowEchelon_with_pivots.jl
function rref_with_pivots!(A::Matrix{T}, ɛ=T <: Union{Rational,Integer} ? 0 : eps(norm(A,Inf))) where T
    nr, nc = size(A)
    pivots = Vector{Int64}()
    i = j = 1
    while i <= nr && j <= nc
        (m, mi) = findmax(abs.(A[i:nr,j]))
        mi = mi+i - 1
        if m <= ɛ
            if ɛ > 0
                A[i:nr,j] .= zero(T)
            end
            j += 1
        else
            for k=j:nc
                A[i, k], A[mi, k] = A[mi, k], A[i, k]
            end
            d = A[i,j]
            for k = j:nc
                A[i,k] /= d
            end
            for k = 1:nr
                if k != i
                    d = A[k,j]
                    for l = j:nc
                        A[k,l] -= d*A[i,l]
                    end
                end
            end
            append!(pivots,j)
            i += 1
            j += 1
        end
    end
    return A, pivots
end

end # module

max_stars_repo_name: tavimori/gf

max_issues_repo_name: tavimori/gf

max_forks_repo_name: tavimori/gf


text:
n_draws = 100
grid = 0:10:500
resampling = BinnedResampling(grid, n_draws)

@test resampling isa BinnedResampling
@test typeof(resampling) <: AbstractBinnedUncertainValueResampling

# Create the resampling scheme. Use kernel density estimates to distribution
# in each bin.
resampling = BinnedResampling(grid, n_draws, bin_repr = UncertainScalarKDE)
@test resampling isa BinnedResampling{UncertainScalarKDE}
resampling = BinnedResampling(UncertainScalarKDE, grid, n_draws)
@test resampling isa BinnedResampling{UncertainScalarKDE}

# Represent each bin as an equiprobably population
resampling = BinnedResampling(grid, n_draws, bin_repr = UncertainScalarPopulation)
@test resampling isa BinnedResampling{UncertainScalarPopulation}
resampling = BinnedResampling(UncertainScalarPopulation, grid, n_draws)
@test resampling isa BinnedResampling{UncertainScalarPopulation}

# Keep raw values for each bin (essentially the same as UncertainScalarPopulation,
# but avoids storing an additional vector of weights for the population members).
resampling = BinnedResampling(grid, n_draws, bin_repr = RawValues)
@test resampling isa BinnedResampling{RawValues}
resampling = BinnedResampling(RawValues, grid, n_draws)
@test resampling isa BinnedResampling{RawValues}

ext: jl


text:
function symeigen(A::AbstractMatrix)
    E, U = eigen(A)
    E, Matrix(U)
end

"""
References:
    * Seeger, M., Hetzel, A., Dai, Z., Meissner, E., & Lawrence, N. D. (2018). Auto-Differentiating Linear Algebra.
"""
function symeigen_back(E::AbstractVector{T}, U, dE, dU; η=1e-40) where T
    all(x->x isa Nothing, (dU, dE)) && return nothing
    η = T(η)
    if dU === nothing
        D = Diagonal(dE)
    else
        F = E .- E'
        F .= F./(F.^2 .+ η)
        dUU = dU' * U .* F
        D = (dUU + dUU')/2
        if dE !== nothing
            D = D + Diagonal(dE)
        end
    end
    U * D * U'
end


text:
function read_ts_from_json(ts_path)
    if isfile(ts_path)
        time_series  = Dict()
        open(ts_path, "r") do f
            time_series  = JSON.parse(f)
        end
    else
        @error "no time series data is available at $(ts_path)"
    end
    return time_series
end

function set_pq_from_timeseries!(mn, ts_data, variable)
    for step in 1:steps
        network = mn["nw"]["$(step)"]
        for idx in keys(ts_data)
            network["load"][idx][variable] = ts_data[idx]["$(step-1)"] / network["baseMVA"]
        end
    end
    return mn
end


text:
module TodoGrabber

include("Todo.jl")
include("grabber.jl");
include("storingdb.jl")

function todo_db(directory::String)
    filenames = get_filenames("./tests")
    todos = find_todos(filenames)
    storetodos(todos)
end # todo_db

end # module


text:
@testset "780.reaching-points.jl" begin
    @test reaching_points(1, 4, 999999997, 4) == true
    @test reaching_points(9, 5, 12, 8) == false
    @test reaching_points(1, 1, 3, 5) == true
    @test reaching_points(3, 3, 12, 9) == true
    @test reaching_points(9, 10, 9, 19) == true
    @test reaching_points(3, 3, 3, 3) == true
end


text:
using Documenter, Swarm

makedocs(
    modules = [Swarm],
    format = :html,
    checkdocs = :exports,
    sitename = "Swarm.jl",
    pages = Any["index.md"]
)

deploydocs(
    repo = "github.com/phelipe/Swarm.jl.git",
)


text:
module TrafficAssignment

# package code goes here
using LightGraphs, Optim

include("load_network.jl")
include("frank_wolfe.jl")

export
        load_ta_network,
        ta_frank_wolfe,
        TA_Data

end # module


text:
# unitary?
v0 = [.2 .4 1]
w0 = PointVortex(v0)
@test rawData(w0) == v0

v1 = [.2 .4 1;0.7 1.5 -1;-.3 1.2 1]
w1 = PointVortex(v1)
@test rawData(w1) == v1

w3 = randPointVortex(1000)
v3 = rawData(w3)
@test v3 == rawData(w3)

# single vortex creation and detection
@test foundNear(1)
@test foundNear(30)


text:
function mnist_load(file_str::String)::Tuple
    dict = Dict()

    open(file_str, "r") do file
        global dict
        dict = JSON.parse(file)
    end

    # <- ::Vector{Tuple{Vector{String},String}}
    train = dict["training"]
    valid = dict["validation"]
    test = dict["test"]

    # -> ::Vector{Tuple{Vector{Float64},Float64}}
    for t in train
        t[1] = parse(Float64,t[1])
        for s in t[0]
            s = parse(Float64,s)
        end
    end
    for t in valid
        t[1] = parse(Float64,t[1])
        for s in t[0]
            s = parse(Float64,s)
        end
    end
    for t in test
        t[1] = parse(Float64,t[1])
        for s in t[0]
            s = parse(Float64,s)
        end
    end
    return (train,valid,test)
end

max_stars_repo_name: ian-double-u/NNfromScratch

max_issues_repo_name: ian-double-u/NNfromScratch

max_forks_repo_name: ian-double-u/NNfromScratch


text:
function StructDefinition{false}(spec::SpecStruct)
    p = Dict(
        :category => :struct,
        :decl => :(
            $(struct_name(spec)) <:
            $(spec.is_returnedonly ? :ReturnedOnly : :(VulkanStruct{$(needs_deps(spec))}))
        ),
    )
    if spec.is_returnedonly
        p[:fields] = map(x -> :($(wrap_identifier(x))::$(idiomatic_julia_type(x))), children(spec))
    else
        p[:fields] = [:(vks::$(spec.name))]
        needs_deps(spec) && push!(p[:fields], :(deps::Vector{Any}))

        foreach(parent_handles(spec)) do member
            handle_type = remove_vk_prefix(member.type)
            name = wrap_identifier(member)
            field_type = is_optional(member) ? :(OptionalPtr{$handle_type}) : handle_type
            push!(p[:fields], :($name::$field_type))
        end
    end

    StructDefinition{false}(spec, p)
end

function Constructor(def::StructDefinition{false})
    spec = def.spec
    cconverted_members = getindex(children(spec), findall(is_semantic_ptr, children(spec).type))
    cconverted_ids = map(wrap_identifier, cconverted_members)
    p = Dict(
        :category => :function,
        :name => name(def),
        :relax_signature => true,
    )
    if needs_deps(spec)
        p[:body] = quote
            $(
                (
                    :($id = cconvert($(m.type), $id)) for
                    (m, id) ∈ zip(cconverted_members, cconverted_ids)
                )...
            )
            deps = [$((cconverted_ids)...)]
            vks = $(spec.name)($(map(vk_call, children(spec))...))
            $(name(def))(vks, deps, $(wrap_identifier.(parent_handles(spec))...))
        end
    else
        p[:body] = :($(name(def))($(spec.name)($(map(vk_call, children(spec))...)), $(wrap_identifier.(parent_handles(spec))...)))
    end
    potential_args = filter(x -> x.type ≠ :VkStructureType, children(spec))
    add_func_args!(p, spec, potential_args)
    Constructor(def, p)
end

StructDefinition{true}(spec::Spec) = StructDefinition{true}(StructDefinition{false}(spec))

function StructDefinition{true}(def::StructDefinition{false})
    spec = def.spec
    @assert !spec.is_returnedonly
    p = Dict(
        :category => :struct,
        :decl => :(
            $(struct_name(spec, true)) <: HighLevelStruct
        ),
        :fields => Expr[],
    )
    for member in filter(!drop_field, children(spec))
        T = hl_type(member)
        if hl_is_optional(member)
            T = :(OptionalPtr{$T})
        end
        push!(p[:fields], :($(wrap_identifier(member))::$T))
    end
    StructDefinition{true}(spec, p)
end

drop_field(x::Spec) = drop_arg(x) || x.name == :sType

function Constructor(T::StructDefinition{false}, x::StructDefinition{true})
    spec = T.spec
    p = Dict(
        :category => :function,
        :name => name(T),
        :args => [:(x::$(name(x)))],
        :short => true,
    )
    args = Expr[]
    kwargs = Expr[]
    for member in filter(!drop_field, children(spec))
        id = wrap_identifier(member)
        id_deref = :(x.$id)
        Tsym = hl_type(member)
        _Tsym = @match innermost_type(Tsym) begin
            GuardBy(is_hl) => @match Tsym begin
                :(Vector{$T}) => :(Vector{$(Symbol(:_, T))})
                :(NTuple{$N,$T}) => :(NTuple{$N,$(Symbol(:_, T))})
                ::Symbol => Symbol(:_, Tsym)
                _ => error("Unhandled conversion for type $Tsym")
            end
            _ => nothing
        end
        ex = isnothing(_Tsym) ? id_deref : :(convert_nonnull($_Tsym, $id_deref))
        if is_optional(member)
            push!(kwargs, ex == id_deref ? ex : Expr(:kw, id, ex))
        else
            push!(args, ex)
        end
    end
    p[:body] = reconstruct_call(Dict(:name => p[:name], :args => args, :kwargs => kwargs))
    Constructor(T, p)
end

function Convert(T::StructDefinition{false}, x::StructDefinition{true})
    p = Dict(
        :category => :function,
        :name => :convert,
        :args => [:(T::Type{$(name(T))}), :(x::$(name(x)))],
        :body => :(T(x)),
        :short => true,
    )
    Convert(T, x, p)
end

function Constructor(def::StructDefinition{true})
    spec = def.spec
    p = Dict(
        :category => :function,
        :name => name(def),
        :args => [],
        :kwargs => [],
        :short => true,
        :relax_signature => true,
    )
    args = []
    for member in filter(!drop_field, children(spec))
        id = wrap_identifier(member)
        if is_optional(member)
            push!(p[:kwargs], Expr(:kw, id, hl_default(member)))
        else
            push!(p[:args], :($id::$(signature_type(hl_type(member)))))
        end
        push!(args, id)
    end
    p[:body] = reconstruct_call(Dict(:name => p[:name], :args => args))
    Constructor(def, p)
end

function hl_default(member::SpecStructMember)
    !embeds_sentinel(member) && return :C_NULL
    @match hl_type(member) begin
        :String => ""
        :(Vector{$t}) => :($t[])
        _ => default(member)
    end
end

hl_is_optional(member::SpecStructMember) = is_optional(member) && !embeds_sentinel(member) && member.name ≠ :pNext

function embeds_sentinel(member::SpecStructMember)
    @match hl_type(member) begin
        :(Vector{$_}) => return !is_optional(member) # optional means C_NULL can be provided, and may have different semantics than for an empty vector
        :String || :UInt32 || :UInt64 || :Int32 || :Int => return true
        _ => nothing
    end

    @match member.type begin
        GuardBy(in(spec_enums.name)) || GuardBy(in(spec_bitmasks.name)) || GuardBy(is_flag) => true
        _ => false
    end
end


text:
include("AbstractTypes/AbstractType.jl")

mutable struct TechnologyHistory <: Buffer
    cashFlow::Array{Real,1} # 7 columns: Revenue, subsidy,VOM, fuelCost, emissionCost ,FOM, CAPEX, NPV
    dispatchHistory::Array{Real,1} # 3 columns: profitableHours,marginalHours,unprofitableHours
    hourlyHistory::Array{Real,2} # 4 columns:  Revenue, VOM, fuelCost, emissionCost
    dailyHistory::Array{Real,2} # 4 columns:  Revenue, VOM, fuelCost, emissionCost
    yearlyHistory::Array{Real,2} # 6 columns:  Revenue, VOM, fuelCost, emissionCost, FOM, CAPEX

    day2hour::Int64
    year2day::Int64

    Reinitialize!::Function # Reinitialize the history after prefiltering

    function TechnologyHistory(systemSettings::SystemSettings)
        self = new()
        self.cashFlow = zeros(8)
        self.dispatchHistory = zeros(3)
        self.day2hour = systemSettings.simulationSettings.day2hour
        self.year2day = systemSettings.simulationSettings.year2day
        self.hourlyHistory = zeros(self.day2hour,5)
        self.dailyHistory = zeros(self.year2day,5)
        self.yearlyHistory = Array{Real,2}(undef,0,7)

        self.Reinitialize! = function()
            self.cashFlow = zeros(8)
            self.dispatchHistory = zeros(3)
            self.hourlyHistory = zeros(self.day2hour,5)
            self.dailyHistory = zeros(self.year2day,5)
            self.yearlyHistory = Array{Real,2}(undef,0,7)
        end
        return self
    end
end

mutable struct DesignProperty <: Property
    availability::Real # range 0.0 - 1.0
    constructionTimeStamp::Real # The year that the plant is built
    constructionTime::Real # lead time
    deconstructionTime::Real # required time to dismantle
    efficiency::Real # Fuel -> electricity, range 0.0 - 1.0
    emissionLevel::Real #ton/MWh
    lifeTime::Real # life expectance of the plant
    installedUnitNumber::Real # The installed unit number of this plant
    unitCapacity::Real # capacity of each unit

    installedCapacity::Real # This is not direct input, calculated use its internal function
    maximumOutPut::Real

    CalculateInstalledCapacity::Function

    function DesignProperty(availability::Real, constructionTimeStamp::Real, constructionTime::Real, deconstructionTime::Real, efficiency::Real, emissionLevel::Real,lifeTime::Real, unitCapacity::Real)
        self = new()
        self.availability = availability
        self.constructionTimeStamp = constructionTimeStamp
        self.constructionTime = constructionTime
        self.deconstructionTime = deconstructionTime
        self.efficiency = efficiency
        self.emissionLevel = emissionLevel
        self.lifeTime = lifeTime
        self.unitCapacity = unitCapacity

        self.CalculateInstalledCapacity =function(unitNumber::Real)
            self.installedUnitNumber = unitNumber
            self.installedCapacity = self.installedUnitNumber * self.unitCapacity
            self.maximumOutPut = self.installedUnitNumber * self.unitCapacity
        end

        return self
    end
end

mutable struct EconomicProperty <: Property
    #Cost of each unit
    unitConstructionCost::Real #€/unit
    unitDeconstructionCost::Real #€/unit
    unitVOMCost::Real #€/MWh
    unitFOMCost::Real #€/unit/year

    domesticLearningRate::Real
    nonDomesticLearningRate::Real

    #Cost of the plant (which is composed of several plants)
    constructionCost::Real #€
    deconstructionCost::Real #€
    VOMCost::Real  #€/MWh
    fuelCost::Real #€/MWh
    emissionCost::Real #€/MWh
    FOMCost::Real #€/year
    yearlyInstallment::Float64  #€/year
    virtualHistory::TechnologyHistory # Store information during virtual auction
    realHistory::TechnologyHistory # Store information during real auction

    CalculatePlantParameter::Function

    function EconomicProperty(unitConstructionCost::Real,unitDeconstructionCost::Real, unitVOMCost::Real,unitFOMCost::Real,systemSettings::SystemSettings)
        self = new()
        self.unitConstructionCost = unitConstructionCost
        self.unitDeconstructionCost = unitDeconstructionCost
        self.unitVOMCost = unitVOMCost
        self.unitFOMCost = unitFOMCost
        self.virtualHistory = TechnologyHistory(systemSettings)
        self.realHistory = TechnologyHistory(systemSettings)

        #Calculate plant parameters if there are more than 1 unit
        self.CalculatePlantParameter = function(unitNumber::Real)
            self.constructionCost = self.unitConstructionCost * unitNumber
            self.deconstructionCost = self.unitDeconstructionCost * unitNumber
            self.VOMCost = self.unitVOMCost
            self.FOMCost = self.unitFOMCost * unitNumber
        end

        return self
    end
end

mutable struct ConventionalPlants <: Technology
    name::Symbol
    type::Symbol
    designProperties::DesignProperty
    economicProperties::EconomicProperty
    offerPrice::Float64 # €/MWh
    offerQuantity::Float64 # MWh
    state::Array{Symbol,1} # state1: online/offline; state2: subsidized/unsubsidized

    ChangeState::Function
    UpdateCost::Function
    MultipleUnits::Function
    CalculateYearlyInstallment::Function
    UpdateHistory!::Function
    AggregateHistory!::Function

    function ConventionalPlants(type::Symbol, designProperties::DesignProperty,economicProperties::EconomicProperty,systemSettings::SystemSettings)
        self = new()
        self.type = type
        self.designProperties = designProperties
        self.economicProperties = economicProperties

        self.UpdateCost = function(newEmissionPrice::Real,newFuelPrice::Real)
            self.economicProperties.emissionCost = newEmissionPrice * self.designProperties.emissionLevel
            self.economicProperties.fuelCost = newFuelPrice / self.designProperties.efficiency
            self.offerPrice = self.economicProperties.unitVOMCost + self.economicProperties.fuelCost
        end

        # NOTE: source: https://en.wikipedia.org/wiki/Mortgage_calculator
        self.CalculateYearlyInstallment = function(systemSettings::SystemSettings)
            interestRate = systemSettings.interestRate
            capital = self.economicProperties.constructionCost *(1-systemSettings.simulationSettings.downpaymentRatio)
            if interestRate > 0.0
                payBackPeriod = self.designProperties.lifeTime
                self.economicProperties.yearlyInstallment = (interestRate*(1 + interestRate) ^ payBackPeriod) / ((1 + interestRate) ^ payBackPeriod - 1) * capital
            end
            # When interest rate is 0, the formula is not valid
            if interestRate == 0.0
                self.economicProperties.yearlyInstallment = capital/self.designProperties.lifeTime
            end
        end

        # NOTE: some necessary initilization of this plant
        self.MultipleUnits = function(unitNumber::Real,systemSettings::SystemSettings)
            self.designProperties.CalculateInstalledCapacity(unitNumber)
            self.economicProperties.CalculatePlantParameter(unitNumber)
            self.offerPrice = self.economicProperties.VOMCost + self.economicProperties.fuelCost
            self.offerQuantity = self.designProperties.maximumOutPut
            self.CalculateYearlyInstallment(systemSettings)
        end

        # NOTE: Record the evaluation history or real life auction history
        self.UpdateHistory! = function(auctionIdentify::Symbol,history::TechnologyHistory,powerFlow::PowerFlow, hour::Int64,day::Int64,subsidyAmount::Real,systemSettings::SystemSettings)
            if auctionIdentify == :real
                repeatTime = 1
            end

            if auctionIdentify == :virtual
                repeatTime = systemSettings.milestoneYear
            end
            marketPrice = powerFlow.marketPrice
            quantity = powerFlow.quantity

            representativeDayWeightFactors = systemSettings.simulationSettings.representativeData[:Weights]
            year2day = systemSettings.simulationSettings.year2day
            day2hour = systemSettings.simulationSettings.day2hour

            #NOTE: update Revenue, VOM, FOM, CAPEX
            thisRevenue = marketPrice * quantity
            thisSubsidy = subsidyAmount * quantity
            thisVOM = self.economicProperties.unitVOMCost * quantity
            thisFuelCost = self.economicProperties.fuelCost * quantity
            thisEmissionCost = self.economicProperties.emissionCost * quantity

            history.hourlyHistory[hour,1] = thisRevenue
            history.hourlyHistory[hour,2] = thisSubsidy
            history.hourlyHistory[hour,3] = thisVOM
            history.hourlyHistory[hour,4] = thisFuelCost
            history.hourlyHistory[hour,5] = thisEmissionCost

            if hour == day2hour
                history.dailyHistory[day,:] = sum(history.hourlyHistory,dims=1) .* representativeDayWeightFactors[day,2]
                history.hourlyHistory = zeros(day2hour,5)
            end

            if day == year2day && hour == day2hour #End of a year
                thisYearlyRevenue = sum(history.dailyHistory[:,1])
                thisYearlySubsidy = sum(history.dailyHistory[:,2])
                thisYearlyVOM = sum(history.dailyHistory[:,3])
                thisYearlyFuelCost = sum(history.dailyHistory[:,4])
                thisYearlyEmissionCost = sum(history.dailyHistory[:,5])
                thisYearlyFOM = self.economicProperties.FOMCost
                thisYearlyHistory = repeat([thisYearlyRevenue thisYearlySubsidy thisYearlyVOM thisYearlyFuelCost thisYearlyEmissionCost thisYearlyFOM self.economicProperties.yearlyInstallment], repeatTime)
                #println("one year",size(thisYearlyHistory))
                #println(day," ",hour)
                history.yearlyHistory = vcat(history.yearlyHistory, thisYearlyHistory)
                history.dailyHistory = zeros(year2day,5)
                #println("Yearly: ", size(history.yearlyHistory))
            end

            #NOTE:update dispatch history
            if marketPrice > self.offerPrice
                history.dispatchHistory[1] += representativeDayWeightFactors[day,2] * repeatTime
            end
            if marketPrice == self.offerPrice
                history.dispatchHistory[2] += representativeDayWeightFactors[day,2] * repeatTime
            end
            if marketPrice < self.offerPrice
                history.dispatchHistory[3] += representativeDayWeightFactors[day,2] * repeatTime
            end
        end

        #NOTE: calculate total 1-> Revenue, 2-> VOM, 3-> FOM, 4-> CAPEX, 5-> NPV
        self.AggregateHistory! = function(history::TechnologyHistory,systemSettings::SystemSettings)
            interestRate = systemSettings.interestRate
            cashFlowHorizon = size(history.yearlyHistory)[1]
            if cashFlowHorizon>=1
                yearVector = collect(0:1:(cashFlowHorizon-1))
                lifeDiscountVector = exp.(-interestRate .* yearVector)
                history.cashFlow[1:7] = lifeDiscountVector' * history.yearlyHistory
                history.cashFlow[8] = sum(history.cashFlow[1:2]) - sum(history.cashFlow[3:7]) - self.economicProperties.constructionCost * systemSettings.simulationSettings.downpaymentRatio
            end
        end
        return self
    end
end

mutable struct RenewablePlants <:Technology
    name::Symbol
    type::Symbol #:OnShoreWind,:OffShoreWind,:Solar
    designProperties::DesignProperty
    economicProperties::EconomicProperty
    offerPrice::Float64
    offerQuantity::Float64
    capacityFactors::Matrix{Float64}

    ChangeState::Function
    MultipleUnits::Function
    TimeSeries::Function
    CalculateYearlyInstallment::Function
    UpdateHistory!::Function
    AggregateHistory!::Function

    function RenewablePlants(type::Symbol,designProperties::DesignProperty,economicProperties::EconomicProperty,capacityFactorImproveRate::Float64,systemSettings::SystemSettings)
        self = new()
        self.type = type
        self.designProperties = designProperties
        self.economicProperties = economicProperties
        self.capacityFactorImproveRate = capacityFactorImproveRate
        self.capacityFactors = systemSettings.simulationSettings.representativeData[self.type]

        # NOTE: source: https://en.wikipedia.org/wiki/Mortgage_calculator
        self.CalculateYearlyInstallment = function(systemSettings::SystemSettings)
            interestRate = systemSettings.interestRate
            capital = self.economicProperties.constructionCost *(1-systemSettings.simulationSettings.downpaymentRatio)
            if interestRate > 0.0
                payBackPeriod = self.designProperties.lifeTime
                self.economicProperties.yearlyInstallment = (interestRate*(1 + interestRate) ^ payBackPeriod) / ((1 + interestRate) ^ payBackPeriod - 1) * capital
                #println(self.economicProperties.yearlyInstallment)
            end
            # When interest rate is 0, the formula is not valid
            if interestRate == 0.0
                self.economicProperties.yearlyInstallment = capital/self.designProperties.lifeTime
                #println(self.economicProperties.yearlyInstallment)
            end
        end

        # NOTE: some necessary initilization of this plant
        self.MultipleUnits = function(unitNumber::Real,systemSettings::SystemSettings)
            self.designProperties.CalculateInstalledCapacity(unitNumber)
            self.economicProperties.CalculatePlantParameter(unitNumber)
            self.offerPrice = self.economicProperties.VOMCost
            #self.offerQuantity = self.designProperties.installedCapacity
            self.CalculateYearlyInstallment(systemSettings)
        end

        self.TimeSeries = function(day::Int64,hour::Int64)
            self.offerQuantity = self.capacityFactors[day,hour] * self.designProperties.maximumOutPut
        end

        # NOTE: Record the evaluation history or real life auction history
        self.UpdateHistory! = function(auctionIdentify::Symbol, history::TechnologyHistory,powerFlow::PowerFlow, hour::Int64,day::Int64,subsidyAmount::Real,systemSettings::SystemSettings)
            if auctionIdentify == :real
                repeatTime = 1
            end

            if auctionIdentify == :virtual
                repeatTime = systemSettings.milestoneYear
            end
            marketPrice = powerFlow.marketPrice
            quantity = powerFlow.quantity

            representativeDayWeightFactors = systemSettings.simulationSettings.representativeData[:Weights]
            year2day = systemSettings.simulationSettings.year2day
            day2hour = systemSettings.simulationSettings.day2hour

            #NOTE: update Revenue, VOM, FOM, CAPEX
            thisRevenue = marketPrice * quantity
            thisSubsidy = subsidyAmount * quantity
            thisVOM = self.economicProperties.VOMCost * quantity
            thisFuelCost = 0
            thisEmissionCost = 0

            history.hourlyHistory[hour,1] = thisRevenue
            history.hourlyHistory[hour,2] = thisSubsidy
            history.hourlyHistory[hour,3] = thisVOM
            history.hourlyHistory[hour,4] = thisFuelCost
            history.hourlyHistory[hour,5] = thisEmissionCost

            if hour == day2hour
                history.dailyHistory[day,:] = sum(history.hourlyHistory,dims=1) * representativeDayWeightFactors[day,2]
                history.hourlyHistory = zeros(day2hour,5)
            end

            if day == year2day && hour == day2hour #End of a year
                thisYearlyRevenue = sum(history.dailyHistory[:,1])
                thisYearlySubsidy = sum(history.dailyHistory[:,2])
                thisYearlyVOM = sum(history.dailyHistory[:,3])
                thisYearlyFuelCost = sum(history.dailyHistory[:,4])
                thisYearlyEmissionCost = sum(history.dailyHistory[:,5])
                thisYearlyFOM = self.economicProperties.FOMCost
                thisYearlyHistory = repeat([thisYearlyRevenue thisYearlySubsidy thisYearlyVOM thisYearlyFuelCost thisYearlyEmissionCost thisYearlyFOM self.economicProperties.yearlyInstallment], repeatTime)
                #println("one year",size(thisYearlyHistory))
                #println(day," ",hour)
                history.yearlyHistory = vcat(history.yearlyHistory, thisYearlyHistory)
                history.dailyHistory = zeros(year2day,5)
                #println("Yearly: ", size(history.yearlyHistory))
            end

            #NOTE:update dispatch history
            if marketPrice > self.offerPrice
                history.dispatchHistory[1] += representativeDayWeightFactors[day,2] * repeatTime
            end
            if marketPrice == self.offerPrice
                history.dispatchHistory[2] += representativeDayWeightFactors[day,2] * repeatTime
            end
            if marketPrice < self.offerPrice
                history.dispatchHistory[3] += representativeDayWeightFactors[day,2] * repeatTime
            end
        end

        #NOTE: calculate total 1-> Revenue, 2-> VOM, 3-> FOM, 4-> CAPEX, 5-> NPV
        self.AggregateHistory! = function(history::TechnologyHistory,systemSettings::SystemSettings)
            interestRate = systemSettings.interestRate
            cashFlowHorizon = size(history.yearlyHistory)[1]
            if cashFlowHorizon>=1
                yearVector = collect(0:1:(cashFlowHorizon-1))
                lifeDiscountVector = exp.(-interestRate .* yearVector)
                history.cashFlow[1:7] = lifeDiscountVector' * history.yearlyHistory
                history.cashFlow[8] = sum(history.cashFlow[1:2]) - sum(history.cashFlow[3:7]) - self.economicProperties.constructionCost * systemSettings.simulationSettings.downpaymentRatio
            end
        end
        return self
    end
end


text:
# get random matrix for least squares problem
A     = ( sprandn(100,13,.1), sprandn(100,200,.2), randn(100,24))
i1     = (1:round(Int64,size(A[1],1)/2))
i2     = (round(Int64,size(A[1],1)/2)+1:size(A[1],1))
pFor1  = LSparam(A[1][i1,:],[])
pFor2  = LSparam(A[1][i2,:],[])

# single LS problem
pFor   = LSparam(A[2],[])
# parallel with dynamic scheduling
pForp  = [pFor1; pFor2];
# remote references
pForRef    = Array{RemoteChannel}(undef,2)
i1     = (1:round(Int64,size(A[3],1)/2))
i2     = (round(Int64,size(A[3],1)/2)+1:size(A[3],1))
A1     = A[3][i1,:]
A2     = A[3][i2,:]
workerList = workers()
nw         = length(workerList)
pForRef[1] = initRemoteChannel(LSparam,workerList[1%nw+1],A1,[])
pForRef[2] = initRemoteChannel(LSparam,workerList[2%nw+1],A2,[])

pFors = (pForp,pFor,pForRef)

for k=1:length(pFors)
	# test getSensMatVec for pFor as $(typeof(pFors[k]))...
	(mt,nt) = getSensMatSize(pFors[k])
	nd      = getNumberOfData(pFors[k])
	(mk,nk) = size(A[k])

	@test mt==nd
	@test mk==mt
	@test nt==nk

	At = getSensMat(randn(nk),pFors[k])
	if At[1] isa Future
	    At = vcat([fetch(Ai) for Ai in At]...)
	end
	@test norm(At-A[k],Inf)/norm(A[k],Inf) < 1e-14
end


text:
# This file was generated, do not modify it. # hide
@from_network mach begin
    mutable struct CompositeModel
        regressor=ridge_model
    end
end


text:
"""
    read_fasta_records(filename)

Read a FASTA file and return records.

# Returns:
- `records::Vector{FASTA.Record}`

"""
function read_fasta_records(filename)
    stream = open(FASTA.Reader, filename)
    records = FASTA.Record[]
    for entry in stream
        push!(records, entry)
    end
    return records
end

"""
    read_fasta(filename)

Read a FASTA file and convert to a given sequence type.

# Returns:
- `seqs::Vector{T}`

"""
function read_fasta(filename; seqtype=DNASequence)
    records = read_fasta_records(filename)
    return seqtype[sequence(r) for r in records]
end

"""
    write_fasta(filename, seqs; names)

Write sequences to a FASTA file.

# Arguments:
- `filename`: file into which to write
- `seqs`: sequences to write
- `names::Vector{String}`: optional list of corresponding names

"""
function write_fasta(filename, seqs; names=String[])
    if length(names) == 0
        names = ["seq_$i" for i in 1:length(seqs)]
    end
    stream = open(FASTA.Writer, filename)
    for i = 1:length(seqs)
        write(stream, FASTA.Record(names[i], seqs[i]))
    end
    close(stream)
end

"""
    read_fastq_records(filename)

Read a FASTQ file and return records.

# Returns:
- `records::Vector{FASTQ.Record}`

"""
function read_fastq_records(filename)
    stream = open(FASTQ.Reader, filename)
    records = FASTQ.Record[]
    for record in stream
        if any([q < 0 for q in FASTQ.quality(record)])
            error("$(FASTQ.identifier(record)) in $filename contains negative phred values")
        end
        push!(records, record)
    end
    return records
end

"""
    read_fastq(filename)

Read a FASTQ file and convert to a given sequence type.

# Returns:
- `seqs::Vector{T}`:
- `phreds::Vector{Vector{Phred}}`: Phred values
- `names::Vector{String}`: sequence names

"""
function read_fastq(filename; seqtype=DNASequence)
    records = read_fastq_records(filename)
    seqs = seqtype[]
    phreds = Vector{Phred}[]
    names = String[]
    for record in records
        push!(seqs, sequence(record))
        push!(phreds, FASTQ.quality(record))
        push!(names, FASTQ.identifier(record))
    end
    return seqs, phreds, names
end

"""
    write_fastq(filename, seqs, phreds; names)

Write sequences to a FASTA file.

# Arguments:
- `filename`: file into which to write
- `seqs`: sequences to write
- `phreds`: corresponding Phred scores
- `names::Vector{String}`: optional list of corresponding names

"""
function write_fastq(filename, seqs, phreds::Vector{Vector{Phred}};
                     names=String[])
    stream = open(FASTQ.Writer, filename)
    i = 0
    if length(names) != length(seqs)
        names = [string("seq_", i) for i in 1:length(seqs)]
    end
    for (s, q, n) in zip(seqs, phreds, names)
        i += 1
        write(stream, FASTQ.Record(n, s, q))
    end
    close(stream)
end


text:
# for convenience (use DISSEQT typealiases instead?)
const Seq = Tuple{String,String}
const Reference = Array{Seq,1}

function askforconfirmation(str::AbstractString)
    while true
        println(str, " (y/n)")
        response = strip(readline())
        lowercase(response) in ["y","yes"] && return true
        lowercase(response) in ["n","no"] && return false
    end
end

# Makes sure folder exists and is empty.
# [Default] Asks for confirmation before removing any files.
# return false if failed
function makecleanfolder(folder::AbstractString, askForConfirmation=true)
    if isdir(folder) && !isempty(readdir(folder))
        (askForConfirmation && askforconfirmation("\"$folder\" folder is non-empty, delete all files?")) || return false
        rm(folder, recursive=true)
    end
    isdir(folder) || mkdir(folder)
    return true
end

# expands (pattern, filename) to (pattern, filename, filepath, localfilepath)
function getreferenceinfo(syn, ref::Tuple{T,String}, referenceFolder::String) where {T}
    path = childpath(syn, referenceFolder, ref[2])
    ref[1], ref[2], path, expanduser(localpath(syn,path,downloadLocation=synapseCacheDir,ifcollision="overwrite.local"))
end
function getreferenceinfo(syn, refs::Vector, referenceFolder::String)
    [getreferenceinfo(syn,r,referenceFolder) for r in refs]
end


text:
# Exported functions

# TODO:
# Add tests for the sensibleness of mi and cmi
# Offer alternative target
# Look into standardising the documentation
# Implement (or find implementation of) Bayesian blocks in Julia,
# 	to not be dependent on the Python module

export get_shrinkage_entropy, get_shrinkage_mi, get_shrinkage_cmi, get_shrinkage_total_correlation

"""
Calculates a James-Stein shrinkage estimate for the entropy of a set
of observed values, via the following steps:
(1) convert the observed values to frequencies of discrete bins
(2) convert the frequencies to probablities using shrinkage estimation
(3) run the probabilities through the entropy formula.
Can calculate joint entropy for up to 3 sets of observed values.

Parameters:

valuesX - dxn Array{Float64,2} - The observed values, where d is
the number of dimensions and n is the number of values.

[valuesY] - dxn Array{Float64,2} - The observed values, where d is
the number of dimensions and n is the number of values.

[valuesZ] - dxn Array{Float64,2} - The observed values, where d is
the number of dimensions and n is the number of values.

Keyword arguments:

[base=2] - the base of the logarithm, which determines the units
[lambda=nothing] - the shrinkage intensity, between 0 and 1 inclusive, which
		will be calculated if this value is not given
[mode="uniformwidth"] - the discretization mode: "uniformwidth", "uniformcount" or
		"bayesianblocks"
"""
function get_shrinkage_entropy(valuesX::Array{Float64,2}; base::Number = 2, lambda = nothing, mode::ASCIIString = "uniformwidth")
	probabilities = getprobabilitiesshrinkage(getfrequencies(valuesX, mode), lambda)
	return applyentropyformula(probabilities, base)
end
function get_shrinkage_entropy(valuesX::Array{Float64,2}, valuesY::Array{Float64,2}; base::Number = 2, lambda = nothing, mode::ASCIIString = "uniformwidth")
	jointprobabilities = getprobabilitiesshrinkage(getjointfrequencies(valuesX, valuesY, mode), lambda)
	return applyentropyformula(jointprobabilities, base)
end
function get_shrinkage_entropy(valuesX::Array{Float64,2}, valuesY::Array{Float64,2}, valuesZ::Array{Float64,2}; base::Number = 2, lambda = nothing, mode::ASCIIString = "uniformwidth")
	jointprobabilities = getprobabilitiesshrinkage(getjointfrequencies(valuesX, valuesY, valuesZ, mode), lambda)
	return applyentropyformula(jointprobabilities, base)
end

"""
Calculates a James-Stein shrinkage estimate for the mutual information
between two sets of observed values.

Parameters:

valuesX - dxn Array{Float64,2} - A set of observed values, where d is
the number of dimensions and n is the number of values.

valuesY - dxn Array{Float64,2} - Another set of observed values, where
d is the number of dimensions and n is the number of values.

Keyword arguments:

[base=2] - the base of the logarithm, which determines the units
[lambda=nothing] - the shrinkage intensity, between 0 and 1 inclusive, which
		will be calculated if this value is not given
[mode="uniformwidth"] - the discretization mode: "uniformwidth", "uniformcount" or
		"bayesianblocks"
"""
# NB "mi" stands for mutual information
function get_shrinkage_mi(valuesX::Array{Float64,2}, valuesY::Array{Float64,2}; base::Number = 2, lambda = nothing, mode::ASCIIString = "uniformwidth")
	entropyX = get_shrinkage_entropy(valuesX, base = base, lambda = lambda, mode = mode)
	entropyY = get_shrinkage_entropy(valuesY, base = base, lambda = lambda, mode = mode)
	entropyXY = get_shrinkage_entropy(valuesX, valuesY, base = base, lambda = lambda, mode = mode)
	return applymutualinformationformula(entropyX, entropyY, entropyXY)
end

"""
Calculates a James-Stein shrinkage estimate for the conditional mutual
information between two sets of observed values, given a third set of
observed values.

Parameters:

valuesX - dxn Array{Float64,2} - A set of observed values, where d is
the number of dimensions and n is the number of values.

valuesY - dxn Array{Float64,2} - Another set of observed values, where
d is the number of dimensions and n is the number of values.

valuesZ - dxn Array{Float64,2} - The set of observed values, that will
be conditioned on, where d is the number of dimensions and n is the
number of values.

Keyword arguments:

[base=2] - the base of the logarithm, which determines the units
[lambda=nothing] - the shrinkage intensity, between 0 and 1 inclusive, which
		will be calculated if this value is not given
[mode="uniformwidth"] - the discretization mode: "uniformwidth", "uniformcount" or
		"bayesianblocks"
"""
# NB "cmi" stands for conditional mutual information
function get_shrinkage_cmi(valuesX::Array{Float64,2}, valuesY::Array{Float64,2}, valuesZ::Array{Float64,2}; base::Number = 2, lambda = nothing, mode::ASCIIString = "uniformwidth")
	entropyZ = get_shrinkage_entropy(valuesZ, base = base, lambda = lambda, mode = mode)
	entropyXZ = get_shrinkage_entropy(valuesX, valuesZ, base = base, lambda = lambda, mode = mode)
	entropyYZ = get_shrinkage_entropy(valuesY, valuesZ, base = base, lambda = lambda, mode = mode)
	entropyXYZ = get_shrinkage_entropy(valuesX, valuesY, valuesZ, base = base, lambda = lambda, mode = mode)
	return applyconditionalmutualinformationformula(entropyZ, entropyXZ, entropyYZ, entropyXYZ)
end

function get_shrinkage_total_correlation(valuesX::Array{Float64,2}, valuesY::Array{Float64,2}, valuesZ::Array{Float64,2}; base::Number = 2, lambda = nothing, mode::ASCIIString = "uniformwidth")
	entropyX = get_shrinkage_entropy(valuesX, base = base, lambda = lambda, mode = mode)
	entropyY = get_shrinkage_entropy(valuesY, base = base, lambda = lambda, mode = mode)
	entropyZ = get_shrinkage_entropy(valuesZ, base = base, lambda = lambda, mode = mode)
	entropyXYZ = get_shrinkage_entropy(valuesX, valuesY, valuesZ, base = base, lambda = lambda, mode = mode)
	return applytotalcorrelationformula(entropyX, entropyY, entropyZ, entropyXYZ)
end

max_stars_repo_name: Tchanders/ShrinkageEstimators

max_issues_repo_name: Tchanders/ShrinkageEstimators

max_forks_repo_name: Tchanders/ShrinkageEstimators


text:
@testset "intersect_lange!" begin
    x = [1, 3, 4, 5, 7, 9]
    y = [2, 3, 5, 6]
    @test MendelImpute.intersect_size_sorted(x, y) == 2
    MendelImpute.intersect_lange!(x, y)
    @test all(x .== [3, 5])
    @test all(y .== [2, 3, 5, 6])

    x = [1, 3, 4, 7]
    y = [2, 3, 5, 6, 7, 10]
    @test MendelImpute.intersect_size_sorted(x, y) == 2
    MendelImpute.intersect_lange!(x, y)
    @test all(x .== [3, 7])
    @test all(y .== [2, 3, 5, 6, 7, 10])

    # allow repeats, although we don't have any in MendelImpute
    x = [3, 4, 7, 7, 7, 10]
    y = [2, 3, 5, 7, 7, 10]
    @test MendelImpute.intersect_size_sorted(x, y) == 4
    MendelImpute.intersect_lange!(x, y)
    @test all(x .== [3, 7, 7, 10])
    @test all(y .== [2, 3, 5, 7, 7, 10])
end


text:
module StateVectorBackend

    include("./StateVectorFunctions.jl")
    include("../backend_Models.jl")
    export execute_backend
    using .StateVectorFunctions
    using .BackendModels

    function execute_backend(n_qregs::Int, gates, model, counts::Int64)
        if model.zero_state == true
            state = zeros(ComplexF64, 2^n_qregs)
            state[1] = 1.
        else
            if length(model.init_state) != 2^n_qregs
                error("Initial state length must be equal to 2^(n_qubits).")
            else
                state = model.init_state
            end
        end

        backend_gates = call_backend_gates(gates)
        states = []
        cregs = []
        for i in 1:counts
            _backend_gates = copy(backend_gates)
            _state = copy(state)
            creg_indices = []
            results = []
            for gate in _backend_gates
                apply!(gate, _state, n_qregs)
                if typeof(gate) == M
                    push!(creg_indices, gate._cregidx)
                    push!(results, gate._result)
                end
            end

            if length(creg_indices) > 0
                creg = zeros(Int64, max(creg_indices...))
                for (i, idx) in enumerate(creg_indices)
                    creg[idx] = results[i]
                end
            else
                creg = zeros(Int64, n_qregs)
            end
            push!(states, _state)
            push!(cregs, creg)
        end
        return states, cregs
    end

    function call_backend_gates(gates)
        backend_gates = []
        for gate in gates
            name = gate._name
            fields = fieldnames(typeof(gate))
            field_names = []
            for i in fields
                push!(field_names, string(i))
            end
            given_property = [getproperty(gate, property) for property in fields]
            code = gate._name * "("
            for prop in given_property[1:end-1]
                code *= string(prop)
                code *= ", "
            end
            code *= ")"
            new_gate = eval(Meta.parse(code))
            push!(backend_gates, new_gate)
        end
        return backend_gates
    end

end

max_stars_repo_name: terasakisatoshi/blueqat-julia


text:
using CSVFiles, CSV, DelimitedFiles
using DataFrames
using CodecZlib
using Random

flatten(x::Array{<:Array,1})= Iterators.flatten(x)|> collect|> flatten
flatten(x::Array{<:Number,1})= x

"""
random_sample(sample_size::Integer = SAMPLE_SIZE)
Generates a random sample of indices from 1 to TOTAL_DATA and bins them with a bin width of sample_size.
"""
function random_sample(sample_size::Integer = SAMPLE_SIZE)
    sample = sort!(rand(1:TOTAL_DATA, sample_size))
    binned = bin(sample, TOTAL_DATA, sample_size)
    return convert(Array{Array{Int32,1}}, binned)
end

"""
bin(data, n, Δ)
Returns data binned into n bins of width Δ
"""
function bin(data, n, Δ)
    bins = [[] for i = 1:Δ:n]
    for i ∈ eachindex(bins)
        for (j, X) ∈ enumerate(data)
            if (i - 1)*Δ < X ≤ i*Δ
                push!(bins[i], X)
            else
                data = data[j:end]
                break
            end
        end
    end
    return bins
end

"""
viewdf(df::DataFrame, row, col, T::Type)
Returns the dataframe value at the specified row and column and converts this value to type T
"""
function viewdf(df::DataFrame, row, col, T::Type)
    return convert(T, df[row, col])
end

"""
export_gz(out, dir, file)
Exports out argument as a .gzip file to "dir/file"
"""
function export_gz(out, dir, file)
    open("$dir/$file", "w") do io
        @info("Opening ", "$dir/$file")
        stream = GzipCompressorStream(io)
        @info("Writing to path...")
        CSV.write(stream, out)
        close(stream)
        @info("Data exported.")
    end
end

"""
export_csv(csv, dir::String, file_name::String)
Exports csv argument as a csv file to "dir/file"
"""
function export_csv(csv, dir::String, file::String)
    open("$dir/$file", "w") do f
        writedlm(f, csv, ',')
    end
end

"""
permute!(data)
Generates a permutation Π, permutes data by Π, and returns Π
"""
function permute!(data)
    Π = randperm(size(data,1))
    data[:] = data[Π]
    return Π
end

"""
split(data, train_percent)
Returns data split into two lists, the first containing the first (train_percent) elements of data, the second containing the last (1 - train_percent) elements of data.
"""
function split(data, train_percent)
    train_data = []
    test_data = []
    for (i,d) in enumerate(data)
        if i/size(data,1) ≤ train_percent
            push!(train_data, d)
        else
            push!(test_data, d)
        end
    end
    return train_data, test_data
end

"""
_mkpaths(path::String, subdirs::AbstractArray{<:String})

Makes paths at path/subdirs[1], ... ,path/subdirs[end]

Returns the first subdirectory path (/path/subdirs[1])
"""
function _mkpaths(path::String, subdirs::AbstractArray{<:String})
    [mkpath("$path/$sd") for sd in subdirs]
    return "$path/$n/$(subdirs[1])"
end

"""
split_export(data, onts, path; train_percent=0.8, csv=nothing)

Splits data by an amount (train_percent) and exports both lists to a .gzip file.

If csv is given, exports csv as .gzip file.
"""
function split_export(data, onts, path; train_percent=0.8, csv=nothing)
    total_data =[(d[1],o) for (d,o) ∈ zip(data, onts)]
    Π = permute!(total_data)
    train_set, test_set = split(total_data, train_percent)
    export_path = _mkpaths(path, ["data", "models/$ONTOLOGY/temp", "models/$ONTOLOGY/logs"])

    train_data = array_to_df(train_set)
    test_data = array_to_df(test_set)
    export_gz(train_data, export_path, "train.gz")
    export_gz(test_data, export_path, "test.gz")

    if !(csv ≡ nothing)
        csv = flatten(csv)[Π]
        train_idx, test_idx = split(csv, train_percent)
        export_csv(train_idx, "$export_path", "train_idx.csv")
        export_csv(test_idx, "$export_path", "test_idx.csv")
    end
end

"""
loaddf(path; train_data=false, test_data=false)
Loads training or testing dataframe
"""
function loaddf(path; train_data=false, test_data=false)
    if train_data
        return DataFrame(load(File(format"CSV", path*"/train.gz")))
    elseif test_data
        return DataFrame(load(File(format"CSV", path*"/test.gz")))
    end
end

"""
Loads dataframe specified by index
"""
loaddf(path, index) = DataFrame(load(File(format"CSV", "path/index.gz")));

"""
function import_csv(import_path::String; file = nothing)
Imports csv file at "import_path/file"
"""
function import_csv(import_path::String; file = nothing)
    csv = AbstractArray{Int}
    !isnothing(file) && (import_path = "$import_path/$file")
    open(import_path) do f
        csv = readdlm(import_path, ',', Int32)
    end
    return csv
end

"""
import_data(import_path::String, dbId::Integer; train_data=false, test_data=false)
Returns data from dataframe at import_path with database id dbID
"""
function import_data(import_path::String, dbId::Integer; train_data=false, test_data=false)
    df = loaddf(import_path, train_data=train_data, test_data=test_data)
    data = df_to_array(df, dbId=dbId, reimport=true)
    @info("IMPORTED ", size(df,1), " POLYTOPES.")
    return data
end

"""
import_split(import_path::String)
Returns training and testing data after importing from import_path
"""
function import_split(import_path::String)
    train_df = loaddf(import_path, train_data=true)
    test_df = loaddf(import_path, test_data=true)
    train_set = df_to_array(train_df, reimport=true)
    test_set = df_to_array(test_df, reimport=true)
    return train_set, test_set
end

"""
array_to_df(data)

Converts array data to dataframe for exporting
"""
function array_to_df(data)
    out = DataFrame()
    for d in data
        ## Need to transpose matrix again for columnar output
        len = convert(Int32, size(d[1]',1))
        delim = zeros(len)
        Δ = hcat(delim, d[1]')
        header = zeros(DIM + 1)
        header[1] += len
        header[2] += d[2][1]
        header[3] += d[2][2]
        header[4] += d[2][3]
        Δ = vcat(header', Δ)
        Δ = convert(DataFrame, Δ)
        append!(out, Δ)
    end
    return out
end

"""
df_to_array(df::DataFrame; dbId=nothing, required=nothing, reimport=false)

Converts dataframe to array data for training, testing, or evaluation
"""
function df_to_array(df::DataFrame; dbId=nothing, required=nothing, reimport=false)
    dbId ≡ nothing && (dbId = 1)
    dfPolyId = (dbId - 1) * SAMPLE_SIZE + 1
    data = []
    onts = []
    currentRow = 1
    skip = viewdf(df, currentRow, 1, Int32)
    dfSize = size(df,1)
    polyMat = []

    if reimport == true
        required = [x for x in 1:SAMPLE_SIZE]
        dfPolyId = 1
    elseif required ≡ nothing
        required = [x for x in ((dbId - 1) * SAMPLE_SIZE + 1 ):(dbId * SAMPLE_SIZE)]
    end

    for reqPolyId in required
        while(reqPolyId != dfPolyId && currentRow < dfSize)
            #get skip
            skip = viewdf(df, currentRow, 1, Int32)
            #increase currentRow
            currentRow += skip + 1
            #iterate polytope count
            dfPolyId += 1
        end
        if currentRow + 1 > dfSize; break; end
        #When we arrive at the correct polytope, grab the matrix and the ontology and add them to the dataset
        skip = viewdf(df, currentRow, 1, Int32)
        !NUMBERS_ONLY && (polyMat = viewdf(df, currentRow+1:currentRow+skip, 2:DIM+1, Array))
        ontology = viewdf(df, currentRow, ontology_index, Int32)
        !reimport && (labels = convert.(Int32,viewdf(df, currentRow:currentRow, 2:4, Array)))
        ############################
        #Debugging and verification
        ############################
        # if reqPolyId == 100000
        #     println(currentRow)
        #     display(df[currentRow:currentRow+skip+1,:])
        # end
        ###########################

        if reimport == false
            if NUMBERS_ONLY == false
                v = viewdf(df, currentRow + 1, 1, Int32)
                p = viewdf(df, currentRow + 2, 1, Int32)
                v_dual = viewdf(df, currentRow + 3, 1, Int32)
                p_dual = viewdf(df, currentRow + 4, 1, Int32)
                polyMat = vcat([v p v_dual p_dual], polyMat)
            else
                polyMat = viewdf(df, currentRow + 1:currentRow + 4, 1:1, Array)'
            end
        else
            NUMBERS_ONLY == true && (polyMat = viewdf(df, currentRow + 1:currentRow + 1, 2:5, Array))
        end
        push!(data, (polyMat',ontology))
        !reimport && push!(onts, labels)
    end
    println("Database: ", dbId)
    println("REQUIRED: $(length(required))")
    println("GOT: $(length(data)) ", length(data) == length(required))
    println("---------------")
    !reimport && (return data, onts)
    return data
end

"""
new_sample(;required = nothing, db_path = DB_PATH,
                     export_path = "None",
                     sample_size = SAMPLE_SIZE,
                     split = false,
                     train_percent = 0.8,
                     db_start = 1,
                     db_end = TOTAL_NUM_DB,
                     get_all=false)
Generates a new data sample from raw database of size sample_size.

Splits data into training and testing sets and exports to a .gzip file at export_path.
"""
function new_sample(;required = nothing, db_path = DB_PATH,
                     export_path = "None",
                     sample_size = SAMPLE_SIZE,
                     split = false,
                     train_percent = 0.8,
                     db_start = 1,
                     db_end = TOTAL_NUM_DB,
                     get_all=false)

    isnothing(required) && (required = random_sample(sample_size))
    get_all && (required = [x for x in 1:sample_size])

    data = []
    onts = []
    for i in db_start:db_end
        df = loaddf(db_path, i);
        dbData, labels = df_to_array(df, dbId=i, required=required[i])
        append!(data, dbData)
        append!(onts,labels)
    end
    @info("OBTAINED: $(size(data,1)) IN TOTAL")

    if export_path != "None"
        if split == true
            split_export(data, onts, export_path; train_percent=train_percent, csv=required)
        else
            n = 1; while(isdir("$export_path/$n.gz")) n+=1; end
            export_gz(data, export_path, "$n.gz")
        end
    else
        @warn("DATA NOT EXPORTED: NO export_path GIVEN. RETURNING INSTEAD.")
        return data, onts
    end
end

"""
augment!(data)
Pads data with zeros until each datapoint is of size MAX_SIZE
"""
function augment!(data)
    X = zeros(DIM,MAX_SIZE,size(data,1))
    Y = []
    augmentedX = []

    for d in data
        diff = MAX_SIZE - size(d[1],2)
        diff > 0 ? push!(augmentedX, [d[1] zeros(DIM, diff)]) : nothing
        push!(Y, d[2])
    end

    for (i,d) in enumerate(augmentedX)
        X[:,:,i] = X[:,:,i] + d
    end

    for i in 1:length(data)
        data[i] = (X[:,:,i], Y[i])
    end
end

"""
unpackage(data)
Returns data and labels as seperate lists
"""
function unpackage(data)
    X = zeros(DIM*MAX_SIZE, size(data,1))
    Y = []
    for (i,d) in enumerate(data)
        X[:,i] = X[:,i] + d[1][:,]
        append!(Y, d[2])
    end

    return X, Y
end

"""
Conversions between euler number and index.
"""
euler_to_index(x) = convert.(Int32, round.(map(y -> 0.5*(y + 960) + 1, x)))
index_to_euler(x) = convert.(Int32, round.(map(y -> 2*(y - 1) - 960, x)))


text:
using CSV
using DataFrames

lcfraw = CSV.File( "/mnt/data/lcf/1718/tab/tab/dvhh_ukanon_2017-18.tab" ) |> DataFrame
lcf = DataFrame()

lcf[!,:age_u_18]=lcfraw[!,:A020]+lcfraw[!,:A021]+lcfraw[!,:A022]+lcfraw[!,:A030]+lcfraw[!,:A031]+lcfraw[!,:A032]
lcf[!,:age_18_plus]=lcfraw[!,:A049]-lcf[!,:age_u_18]
lcf[!,:tenure_type] =  lcfraw[!,:A122]
lcf[!,:region] =  lcfraw[!,:Gorx]
lcf[!,:economic_pos] =  lcfraw[!,:A093]
lcf[!,:age_of_oldest] =  lcfraw[!,:a065p]

lcf[!,:total_consumpt] =  lcfraw[!,:P600t]
lcf[!,:food_and_drink] =  lcfraw[!,:P601t]
lcf[!,:alcohol_tobacco] =  lcfraw[!,:P602t]
lcf[!,:clothing] =  lcfraw[!,:P603t]
lcf[!,:housing] =  lcfraw[!,:P604t]
lcf[!,:household_goods] =  lcfraw[!,:P605t]
lcf[!,:health] =  lcfraw[!,:P606t]
lcf[!,:transport] =  lcfraw[!,:P607t]
lcf[!,:communication] =  lcfraw[!,:P608t]
lcf[!,:recreation] =  lcfraw[!,:P609t]
lcf[!,:education] =  lcfraw[!,:P610t]
lcf[!,:restaurants_etc] =  lcfraw[!,:P611t]
lcf[!,:miscellaneous] =  lcfraw[!,:P612t]
lcf[!,:non_consumption] =  lcfraw[!,:P620tp]
lcf[!,:total_expend] =  lcfraw[!,:P630tp]
lcf[!,:equiv_scale] =  lcfraw[!,:OECDSc]
lcf[!,:weekly_net_inc] =  lcfraw[!,:P389p]

deleterows!( lcf, (lcf[!,:total_consumpt] .< 0.0) )

deleterows!( lcf, (lcf[!,:total_expend] .<0))

deleterows!( lcf, (lcf[!,:housing] .<0))

deleterows!( lcf, (lcf[!,:weekly_net_inc] .< 0.0) )

CSV.write( "/home/graham_s/lcf017_8.tab", lcf, delim='\t' )

lcf[!,:l_total_cons] = log.(lcf[!,:total_consumpt ])
lcf[!,:l_food] = log.(lcf[!,:food_and_drink ])
lcf[!,:l_alc_tob] = log.(lcf[!,:alcohol_tobacco ])
lcf[!,:l_cloth] = log.(lcf[!,:clothing ])
lcf[!,:l_housing] = log.(lcf[!,:housing ])
lcf[!,:l_h_goods] = log.(lcf[!,:household_goods ])
lcf[!,:l_health] = log.(lcf[!,:health ])
lcf[!,:l_transport] = log.(lcf[!,:transport ])
lcf[!,:l_comms] = log.(lcf[!,:communication ])
lcf[!,:l_rec] = log.(lcf[!,:recreation ])
lcf[!,:l_educ] = log.(lcf[!,:education ])
lcf[!,:l_rest_etc] = log.(lcf[!,:restaurants_etc ])
lcf[!,:l_misc] = log.(lcf[!,:miscellaneous ])
lcf[!,:l_non_cons] = log.(lcf[!,:non_consumption ])
lcf[!,:l_total_exp] = log.(lcf[!,:total_expend ])
lcf[!,:l_net_inc] = log.(lcf[!,:weekly_net_inc ])

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
using Base.Test
using ReverseDiffSource

@printf("*********************\njulia version = %s, dir %s\n", Base.VERSION_STRING, pwd())

typealias Vec Vector{Float64}
typealias Mat Matrix{Float64}
typealias SymMat Symmetric{Float64,Matrix{Float64}}

AAt(A::Mat) = Symmetric(A*A')
AtA(A::Mat) = Symmetric(A'*A)
sumsq(x::Vec) = norm(x)^2

# Make matrix from diagonal and strict lower triangle,
# e.g. D = [d11 d22 d33 d44]
#      LT = [L21 L31 L32 L41 L42 L43]
# Outputting
#  [d11   0   0   0]
#  [L21 d12   0   0] # row r: Ls starting at sum_i=1^r
#  [L31 L32 d33   0]
#  [L41 L42 L43 d44]
function ltri_unpack(D::Vec, LT::Vec)
  d=length(D)
  make_row(r::Int, L) = hcat(reshape([ L[i] for i=1:r-1 ],1,r-1), D[r], zeros(1,d-r))
  row_start(r::Int) = div((r-1)*(r-2),2)
  inds(r) = row_start(r)+(1:r-1)
  vcat([ make_row(r, LT[inds(r)]) for r=1:d ]...)
end

ltri_unpack(D, L) = ltri_unpack([Float64(d) for d in D], [Float64(l) for l in L])

@test [11 0 0 0; 21 22 0 0 ; 31 32 33 0 ; 41 42 43 44] == ltri_unpack([11 22 33 44], [21 31 32 41 42 43])

LL = ltri_unpack([1 2 3 4], [21 31 32 41 42 43])
@printf("An example lower triangle made from diag and LT=\n%s\n", LL)

function ltri_pack{M<:AbstractMatrix}(L::M)
  d=size(L,1)

  make_row(r::Int, L) = hcat(reshape([ L[i] for i=1:r-1 ],1,r-1), D[r], zeros(1,d-r))
  row_start(r::Int) = (r-1)*(r-2)/2
  diag(L), hcat([L[r,1:r-1] for r=1:d ]...)
end

ltri_pack{T}(L::LowerTriangular{T, Matrix{T}}) = ltri_pack(full(L)) ## Until  packed storage is implemented for those

@test ltri_unpack(ltri_pack(LL)...) == LL

@printf("pak=%s\n", ltri_pack(LL))

# Conventional GMM.
# This doesn't even use inverse covariance, because as soon as you
# start down that route, you may as well go for lpGMM below.
type GMM
  n::Int           # number of Gaussians
  d::Int           # dimension of Gaussian
  alphas::Vec      # weights: n, require sum(alphas)==1
  mus::Array{Vec}  # means: n, each dx1
  sigmas::Array{SymMat}  # covariances: n, each dxd symmetric positive definite
end

function log_likelihood(g::GMM, x::Vec)
  total = 0
  for k=1:g.n
    mean = g.mus[k]
    weight =  g.alphas[k]
    Σ = g.sigmas[k]
    mahalanobis = dot(mean - x, inv(Σ) * (mean - x))
    @printf("m=%s, ", det(inv(full(Σ))))
    total += weight / sqrt(det(2pi*full(Σ))) * exp(-0.5*mahalanobis)
  end
  @printf("w=%f\n", total)
  log(total)
end

n=3
d=2
alphas=rand(n); alphas /= sum(alphas);
mus=[randn(d) for k=1:n]
sigmas=[AAt(randn(d,d)) for k=1:n]
test_gmm = GMM(n,d,alphas,mus,sigmas)
@printf("An example gmm = %s\n", test_gmm)

x = randn(d) # Test point

ll0 = log_likelihood(test_gmm, x)
@printf("ll0=%f\n", ll0)

##########################################################################
# Log-parametrized GMM.
# Weights are strictly positive, covariances are parameterized by their inverse
# square roots (lower triangular).
type lpGMM
  n::Int           # number of Gaussians
  d::Int           # dimension of Gaussian
  alphas::Vec      # log weights: n
  mus::Array{Vec}  # means: n, each dx1
  LDs::Array{Vec}  # square-root-inverse-covariances, log(diagonal): n, each d x 1
  LTs::Array{Vec}  # square-root-inverse-covariances, lower triangle: n, each d*(d-1)/2 x 1
end

# Convert simple GMM to lpGMM
function lpGMM(g::GMM)
  LTs = Array{Vec}(g.n)
  LDs = Array{Vec}(g.n)
  for k=1:g.n
    L = inv(chol(g.sigmas[k].data, Val{:L}))
    D, T = ltri_pack(L)
    LDs[k], LTs[k] = vec(log(D)), vec(T)
  end
  lpGMM(g.n,g.d,log(g.alphas),g.mus,LDs,LTs)
end

# Convert log-parameterized-GMM to simple GMM UnivariateGMM
function GMM(l::lpGMM)
  alphas::Vec = exp(l.alphas)/sum(exp(l.alphas))
  mus::Array{Vec} = l.mus
  Ls = [ltri_unpack(exp(l.LDs[i]), l.LTs[i]) for i=1:l.n]
  sigmas::Array{SymMat} = map(A->inv(Symmetric(A'*A)), Ls)
  GMM(l.n,l.d,alphas,mus,sigmas)
end

g = lpGMM(test_gmm)

#@printf("gmm=%s\n**\n", GMM(g))

########################################

const halflog2π = log(2π)/2

function log_likelihood(g::lpGMM, x::Vec)
  total = 0
  weight_normalizer = sum(exp(g.alphas))
  for k=1:g.n
    InvLowerTriangle = ltri_unpack(exp(g.LDs[k]), g.LTs[k])
    mean = g.mus[k]
    weight =  exp(g.alphas[k])/weight_normalizer # Weights parametrized as logs
    mahalanobis = sumsq(InvLowerTriangle * (mean - x))
    @printf("m=%s, ", det(InvLowerTriangle'*InvLowerTriangle))
    total += weight * det(InvLowerTriangle) * exp(-0.5*mahalanobis)
  end
  @printf("tot=%f, ", total)
  log(total) - halflog2π*g.d
end

ll1 = log_likelihood(g, x)
@printf("ll0=%f, ll1=%f, rat=%f\n", ll0, ll1, ll0/ll1)
@test_approx_eq_eps ll0 ll1 1e-12

########################################
## A better implementation, with some factorizing and logsumexp:

# logsumexp()
function logsumexp_both(x::Array{Float64,1})
  A = maximum(x);
  ema = exp(x-A);
  sema = sum(ema);
  l = log(sema) + A;
  Jacobian = ema/sema;
  return (l, Jacobian);
end
logsumexp(x::Array{Float64,1}) = logsumexp_both(x)[1]
r = rand(5);
@test_approx_eq_eps logsumexp(r) log(sum(exp(r))) 1.0e-8

@deriv_rule logsumexp(x::AbstractArray)  x  logsumexp_both(x)[2].*ds

function logsumexp_1st(x::Array{Float64,1})
  A = maximum(x);
  ema = exp(x-A);
  sema = sum(ema);
  log(sema) + A;
end
rdiff(logsumexp_1st, (zeros(3),)).code

mahal(L::Mat, mu::Vec, x::Vec) = sumsq(L*(x-mu))

function log_likelihood_2(g::lpGMM, x::Vec)

  d_mahal = [0.5*mahal(ltri_unpack(exp(g.LDs[i]), g.LTs[i]), g.mus[i], x) for i in 1:g.n]

  determinants = [sum(g.LDs[i]) for i in 1:g.n]

  logsumexp(g.alphas + determinants - d_mahal) - logsumexp(g.alphas) - halflog2π*g.d
end

ll2 = log_likelihood_2(g, x)
@printf("ll0=%f, ll2=%f, rat=%f\n", ll0, ll2, ll0/ll2)
@test_approx_eq_eps ll0 ll2 1e-12

# Matrix version, in case needed
function log_likelihood_mat(alphas::Vec, mus::Mat, LDs::Mat, LTs::Mat, x::Vec)
  n = length(alphas)
  d = size(mus,1)
  L(i::Int) = ltri_unpack(exp(LDs[:,i]), LTs[:,i])
  d_mahal = [0.5*mahal(L(i), mus[:,i], x) for i in 1:n]

  determinants = [sum(LDs[:,i]) for i in 1:n]

  logsumexp(alphas + determinants - d_mahal) - logsumexp(alphas) - halflog2π*d
end

ll3 = log_likelihood_mat(g.alphas, hcat(g.mus...), hcat(g.LDs...), hcat(g.LTs...), x)
@printf("ll0=%f, ll3=%f, rat=%f\n", ll0, ll3, ll0/ll3)

using Stats
using Distributions

function draw(g::GMM)
  rng = Distributions.Categorical(exp(alphas)/sum(exp(alphas)))
  k = rand(rng)
  rng = Distributions.MultivariateNormal(mus[k], make_L(exp(g.LDs[k]), g.LTs[k]))
  x::Vec = rand(rng)
end

#   using Gadfly
#   x = [draw(g) for k in 1:100]
#   Gadfly.plot(sin, 0,1)

###################
# Define the derivative of sumsq
rdiff(sumsq, (zeros(3),))

@deriv_rule sumsq(x)   x   2x

sumsq2(x::Number,y::Number) = x^2 + y^2
sumsq2_rdiff = rdiff(sumsq2, (0.,0.))

#d1_sumsq2(x,y) = sumsq2_rdiff(x,y)[2]
#d2_sumsq2(x,y) = sumsq2_rdiff(x,y)[3]
#@deriv_rule sumsq2(x,y)   x   d1_sumsq2(x,y)
#@deriv_rule sumsq2(x,y)   y   d2_sumsq2(x,y)

rosenbrock(x::Float32,y::Float32) = sumsq2(1 - x, 10*(x - y^2))
rosenbrock_diff = rdiff(rosenbrock, (0.,0.))
rosenbrock_diff(1.,2.)

rosenbrockN_ex = :((1 - x[1])^2 + 100*sumsq(x[2:end] - x[1:end-1].^2))
rdiff(rosenbrockN_ex, x=zeros(4))

mahal_ex= :(mahal(ltri_unpack(x[1:3], x[4:6]), x[7:9], x[10:12]))

@eval mahal_new(x) = $mahal_ex
mahal_new(rand(12))

f=rdiff(mahal, (zeros(3,3),zeros(3),zeros(3),), )
f.
f(rand(3,3), rand(3), rand(3))

# This one not so good: unmanaged type int
lla(alphas::Vec) = log_likelihood_mat(alphas, hcat(g.mus...), hcat(g.LDs...), hcat(g.LTs...), x)
dx = rdiff(lla, alphas=g.alphas)
dx([-1,-248,-240,-1,-2])

Pkg.dir("ReverseDiffSource")
#           Query(

# abstract Visitor

# type myvisitor <: Visitor
#   index::Int64

#   myvisitor() = new(0)
# end

# v = myvisitor()

# function visit(v::myvisitor, obj::Number)
#   ++v.index
#   println(v)
# end

# function visit(v::Visitor, obj::Any)
#   fields = names(obj)
#   if length(fields) == 0
#     println("Cannot handle a [", typeof(obj), "]")
#   else
#     for n in names(obj)
#       visit(v, obj.(n))
#     end
#   end
# end

# visit(v, gmm_test)

# methods(Array)

max_stars_repo_name: dsyme/ADBench

max_issues_repo_name: dsyme/ADBench

max_forks_repo_name: dsyme/ADBench


text:
const opts = Base.JLOptions()
const inline_flag = opts.can_inline == 1 ? `` : `--inline=no`
const cov_flag = (opts.code_coverage == 1) ? `--code-coverage=user` :
                 (opts.code_coverage == 2) ? `--code-coverage=all` :
                 ``

function run_test(script)
    srvrscript = joinpath(dirname(@__FILE__), script)
    srvrcmd = `$(joinpath(JULIA_HOME, "julia")) $cov_flag $inline_flag $script`
    println("Running tests from ", script, "\n", "="^60)
    ret = run(srvrcmd)
    println("Finished ", script, "\n", "="^60)
    nothing
end

ENV["NODE_META_IMPL"] = "DagScheduler.ShmemMeta.ShmemExecutorMeta"
ENV["CLUSTER_META_IMPL"] = "DagScheduler.SimpleMeta.SimpleExecutorMeta"
println("Running with ShmemMeta and SimpleMeta")
println("===================================")
run(`ipcs -a`)
println("===================================")
run_test("runtests_master_only.jl")
run_test("runtests_1node.jl")
run_test("runtests_2node.jl")
println("===================================")
run(`ipcs -a`)

ENV["NODE_META_IMPL"] = "DagScheduler.ShmemMeta.ShmemExecutorMeta"
ENV["CLUSTER_META_IMPL"] = "DagScheduler.FdbMeta.FdbExecutorMeta"
println("Running with ShmemMeta and FdbMeta")
println("===================================")
run(`ipcs -a`)
println("===================================")
run_test("runtests_1node.jl")
run_test("runtests_2node.jl")
println("===================================")
run(`ipcs -a`)

ENV["NODE_META_IMPL"] = "DagScheduler.FdbMeta.FdbExecutorMeta"
ENV["CLUSTER_META_IMPL"] = "DagScheduler.FdbMeta.FdbExecutorMeta"
println("Running with FdbMeta")
println("===================================")
run_test("runtests_fdb_queue.jl")
run_test("runtests_fdb_dict.jl")
run_test("runtests_master_only.jl")


text:
### Implements a list of RLEVectors like Bioconductor's RLEList
### Not doing compressed versions where the latter concatenates the underlying
###   data into one vector for cheap vectorization

type RLEVectorList{T1,T2 <: Integer} <: AbstractArray{T1, 1}
    els::Vector{ RLEVector{T1,T2} }
    function RLEVectorList(rles...)
        rvl = ""
        return(rvl)
    end
end


text:
#
# Computes the total potential energy of the system using r2
#

function utotal3!(ft, p, data)
    first_atom = zeros(Int64, data.nc, data.nc)
    next_atom = zeros(Int64, data.N)
    makell!(first_atom, next_atom, p, data)
    ut = 0.
    for i = 1:length(ft)
        ft[i] .= 0
    end
    for ip = 1:data.N
        icell = trunc(Int64, p[ip][1]/data.cutoff) + 1
        jcell = trunc(Int64, p[ip][2]/data.cutoff) + 1
        for ic = icell-1:icell+1
            iw = wrap_cell(data.nc, ic)
            for jc = jcell-1:jcell+1
                jw = wrap_cell(data.nc, jc)
                jp = first_atom[iw,jw]
                while jp > 0
                    if ip < jp
                        r2 = (p[ip][1]-p[jp][1])^2 + (p[ip][2]-p[jp][2])^2 # r2pbc(p[ip], p[jp], data.bdim)
                        if r2 ≤ data.cutoff2
                            up, fp = upair3(r2, p[ip], p[jp], data)
                            ut += up
                            ft[ip] .+= fp
                            ft[jp] .-= fp
                        end
                    end
                    jp = next_atom[jp]
                end
            end
        end
    end
    return ut
end

function utotal3(p, data)
    ft = [ Array{Float64}(undef, 2) for i=1:length(p) ]
    return utotal3!(ft, p, data)
end


text:
"""
    sphere(D)

The well-known D-dimensional Sphere function.
"""
function sphere(D=10)
    # Objective function
    f(x) = sum(x.*x)

    bounds = Array([-100.0ones(D) 100.0ones(D)]')

    x = zeros(D)
    return f, bounds, [generateChild(x, f(x)) ]

end

"""
    discus(D)

The well-known D-dimensional Discus function.
"""
function discus(D = 10)
    # Objective function
    f(x) = 1e6x[1].^2 + sum(x[2:end] .^2)

    bounds = Array([-10.0ones(D) 10.0ones(D)]')

    x = zeros(D)
    return f, bounds, [generateChild(x, f(x))]
end

"""
    rastrigin(D)

The well-known D-dimensional Rastrigin function.
"""
function rastrigin(D = 10)

    # Objective function
    f(x) = 10D+ sum(x.*x - 10cos.(2π*x))

    bounds = Array([-5.0ones(D) 5.0ones(D)]')

    x = zeros(D)
    return f, bounds, [generateChild(x, f(x)) ]

end


text:
using StatsBase
using RDatasets
w = WeightVec([1.,2.,3.],6.)
eltype(wv)
length(wv)
isempty(wv)
values(wv)

max_stars_repo_path: Chapter 3/2.jl

max_stars_repo_name: PacktPublishing/Julia-for-data-science

max_issues_repo_path: Chapter 3/2.jl

max_issues_repo_name: PacktPublishing/Julia-for-data-science

max_forks_repo_path: Chapter 3/2.jl

max_forks_repo_name: PacktPublishing/Julia-for-data-science


text:
# --------------------------------------------
# Process regional projections and observations
# - Read TG estimates
# - Estimate trend + acceleration using Hector
# - Extrapolate data from tide gauges
# - Read NCA5 projections
# --------------------------------------------

module RegionalProjections
using NetCDF
using DelimitedFiles
using Plots
using Statistics
using NCDatasets
using MAT
using Interpolations
using CSV
dir_code = homedir()*"/Projects/2021_NOAA_TR/Code/"
include(dir_code*"Masks.jl")
include(dir_code*"Hector.jl")
include(dir_code*"ProcessObservations.jl")
function RunRegionalProjections(settings)
    println("\nRegional observations, trajectories and scenarios...")
    local_obs  = ProcessObservations.ReadLocalObs(settings)
    region_obs = ProcessObservations.ComputeRegionObs(local_obs,settings)
    ComputeRegionalTrajectory!(region_obs,settings)
    NCA5_regional = ReadNCA5Regional(settings)
    correct_baselines!(region_obs,NCA5_regional,settings)
    SaveRegionData(region_obs,NCA5_regional,settings)
    println("Regional observations, trajectories and scenarios done\n")
    return nothing
end

function ComputeRegionalTrajectory!(region_obs,settings)
    # -----------------------------------------------
    # Compute the trajectory for regional sea level
    # and extrapolate
    # -----------------------------------------------

    println("  Computing trajectories...")
    years_estimate = intersect(settings["years_trajectory"],settings["years_tg"])
    y_acc = findall(in(years_estimate),settings["years_tg"])

    amat = ones(length(y_acc),3)
    amat[:,2] = settings["years_tg"][y_acc] .- mean(settings["years_tg"][y_acc])
    amat[:,3] = (settings["years_tg"][y_acc] .- mean(settings["years_tg"][y_acc])).^2

    amat_extend = ones(length(settings["years_trajectory"]),3)
    amat_extend[:,2] .= settings["years_trajectory"] .- mean(settings["years_tg"][y_acc])
    amat_extend[:,3] = (settings["years_trajectory"] .- mean(settings["years_tg"][y_acc])).^2
    for region in settings["regions"]
        trend_file = Hector.EstTrend(settings["years_tg"][y_acc],region_obs[region]["rsl"][y_acc];accel=true,model="Powerlaw",SA=false,SSA=false,monthly=false,tref=mean(settings["years_tg"][y_acc]))
        μ_sol = [trend_file["bias"],trend_file["trend"],trend_file["accel"]]
        σ_sol= [trend_file["bias_sigma"],trend_file["trend_sigma"],trend_file["accel_sigma"]]
        sol_arr = reshape(μ_sol,(1,3)) .+ randn(5000,3) .* reshape(σ_sol,(1,3))
        trajectory_arr = zeros(Float32,length(settings["years_trajectory"]),5000)
        [trajectory_arr[:,i] = amat_extend * @views sol_arr[i,:] for i ∈ 1:5000]
        trajectory_stats = zeros(Float32,length(settings["years_trajectory"]),3)
        [trajectory_stats[t,:] = quantile((@views trajectory_arr[t,:]),[0.17,0.50,0.83]) for t ∈ 1:length(settings["years_trajectory"])]
        region_obs[region]["rsl_trajectory"] = trajectory_stats
        region_obs[region]["trend"] = [μ_sol[2]-σ_sol[2],μ_sol[2],μ_sol[2]+σ_sol[2]]
        region_obs[region]["accel"] = 2 .*[μ_sol[3]-σ_sol[3],μ_sol[3],μ_sol[3]+σ_sol[3]]
    end
    return nothing
end

function ReadNCA5Regional(settings)
    println("  Reading NCA5 scenarios...")

    # Read mask and grid information
    mask = Masks.ReadMask(settings)
    area = ComputeGridArea(mask["ϕ"],mask["θ"])

    fn = settings["dir_NCA5"]*"NCA5_Low_grid.nc"
    lon_NCA5 = ncread(fn,"lon")
    lat_NCA5 = ncread(fn,"lat")
    years_NCA5 = convert.(Float32,ncread(fn,"years",start=[1],count=[-1]))
    pct_NCA5 = convert.(Float32,ncread(fn,"percentiles"))

    # Read and interpolate the data
    # Structure NCA5_projections[scenario][process][region]
    # Store in dictionary
    NCA5_regional = Dict()
    NCA5_regional["years"] = years_NCA5
    NCA5_regional["percentiles"] = pct_NCA5
    for scn in settings["NCA5_scenarios"]
        println("   Scenario "*scn*"...")
        NCA5_regional[scn] = Dict()
        fn = settings["dir_NCA5"]*"NCA5_"*scn*"_grid.nc"
        for prc in settings["processes"]
            NCA5_regional[scn][prc] = Dict()
            NCA5_prc_raw = convert.(Float32,ncread(fn,prc,start=[1,1,1,1],count=[-1,-1,-1,-1]));
            NCA5_prc_int = LinearInterpolation((lon_NCA5, lat_NCA5,years_NCA5,pct_NCA5),NCA5_prc_raw,extrapolation_bc = Line())(mask["ϕ"],mask["θ"],years_NCA5,pct_NCA5);
            for region in settings["regions"]
                NCA5_regional[scn][prc][region] = (sum((@. NCA5_prc_int*area*mask[region]),dims=(1,2)) / sum(area .* mask[region]))[1,1,:,:]
            end
        end
    end
    return NCA5_regional
end

function ComputeGridArea(lon,lat)
    gridsize = abs(lat[2]-lat[1])
    area = @. deg2rad(gridsize) * (sind(lat+gridsize/2)-sind(lat-gridsize/2)) * 6371000^2
    return repeat(area',size(lon,1))
end

function correct_baselines!(region_obs,NCA5_regional,settings)
    # Correct baseline: baseline for everything is trajectory value in 2000
    traj_baseline = findall(in(settings["years_baseline"]),settings["years_trajectory"])
    traj_nca5_start = findfirst(settings["years_trajectory"].==2005)

    # TG and altimetry average
    traj_tg_bl = findall(in(intersect(settings["years_trajectory"],settings["years_tg"])),settings["years_trajectory"])
    obs_tg_bl  = findall(in(intersect(settings["years_trajectory"],settings["years_tg"])),settings["years_tg"])

    obs_alt_bl = findall(isfinite,region_obs["USA"]["alt"])
    traj_alt_bl = findall(in(intersect(settings["years_trajectory"],settings["years_tg"][obs_alt_bl])),settings["years_trajectory"])
    for region in settings["regions"]
        region_obs[region]["rsl_trajectory"] .-= mean(region_obs[region]["rsl_trajectory"][traj_baseline,2])
        region_obs[region]["rsl"] = region_obs[region]["rsl"] .- mean(filter(isfinite,region_obs[region]["rsl"][obs_tg_bl])) .+ mean(region_obs[region]["rsl_trajectory"][traj_tg_bl,2])
        region_obs[region]["alt"] = region_obs[region]["alt"] .- mean(region_obs[region]["alt"][obs_alt_bl]) .+ mean(region_obs[region]["rsl_trajectory"][traj_alt_bl,2])
        for scenario in settings["NCA5_scenarios"]
            NCA5_regional[scenario]["total"][region] .+= region_obs[region]["rsl_trajectory"][traj_nca5_start,2]
            NCA5_regional[scenario]["verticallandmotion"][region] .+= region_obs[region]["rsl_trajectory"][traj_nca5_start,2]
        end
    end
    return nothing
end

function SaveRegionData(region_obs,NCA5_regional,settings)
    println("  Saving...")
    fh = Dataset(settings["fn_proj_reg"],"c")
    fh.attrib["title"] = "Regional projections"
    fh.attrib["description"] = "Regional projections, trajectory, and observations for the interagency report: Global and Regional Sea Level Rise Scenarios for the United States: Updated Mean Projections and Extreme Water Level Probabilities Along U.S. Coastlines"
    fh.attrib["processes"] = "AIS: Antarctic Ice Sheet, GIS: Greenland, glaciers: Glaciers and Ice Caps, landwaterstorage: Liquid water storage changes on land, oceandynamics: Ocean dynamics and global thermosteric expansion, verticallandmotion: Vertical land motion, total: All processes combined."
    fh.attrib["regions"] = "USA: Contiguous United States, NE: Northeast Atlantic Coast, SE: Southeast Atlantic Coast, EGOM: Eastern Gulf of Mexico, WGOM: Western Gulf of Mexico, SW: Southwest Pacific Coast, NW: Northwest Pacific Coast, PAC: Hawaiian Islands, CAR: Caribbean Islands, NAL: North Alaska, SAL: South Alaska."
    defDim(fh,"years", length(settings["years"]))
    defDim(fh,"percentiles",3)
    defDim(fh,"region", length(settings["regions"]))

    defVar(fh,"region",settings["regions"],("region",),deflatelevel=5, attrib = Dict("description" => "Region abbreviation for which the projection has been made.", "units" => "-"))
    defVar(fh,"years",settings["years"],("years",),deflatelevel=5,attrib = Dict("units" => "years"))
    defVar(fh,"percentiles",convert.(Int32,[17,50,83]),("percentiles",),deflatelevel=5, attrib = Dict("description" => "The percentile of the projection. The 50th percentile shows the median projection, and the 17th and 83rd show the upper- and lower bound on the 1 sigma level.", "units" => "-"))

    # Write observations
    rsl_obs_array = zeros(Float32,length(settings["regions"]),length(settings["years"])) .* NaN32
    alt_obs_array = zeros(Float32,length(settings["regions"]),length(settings["years"])) .* NaN32

    for (region_idx,region) ∈ enumerate(settings["regions"])
        rsl_obs_array[region_idx,:] = LinearInterpolation((convert.(Float32,settings["years_tg"])),region_obs[region]["rsl"],extrapolation_bc=NaN32)(settings["years"])
        alt_obs_array[region_idx,:] = LinearInterpolation((convert.(Float32,settings["years_tg"])),region_obs[region]["alt"],extrapolation_bc=NaN32)(settings["years"])
    end
    defVar(fh,"rsl_obs",rsl_obs_array,("region","years"),deflatelevel=5,attrib = Dict("description" => "Observed regional sea-level changes, relative to year 2000. (tide-gauge observations)", "units" => "mm"))
    defVar(fh,"alt_obs",alt_obs_array,("region","years"),deflatelevel=5,attrib = Dict("description" => "Observed regional sea-level changes, relative to year 2000. (altimetry observations)", "units" => "mm"))

    # Write trajectory
    trajectory = zeros(Float32,length(settings["regions"]),length(settings["years"]),3) .* NaN32
    [trajectory[region_idx,:,:] = (LinearInterpolation((settings["years_trajectory"],[1.0f0:3.0f0...]),region_obs[region]["rsl_trajectory"],extrapolation_bc=NaN32)[settings["years"],[1.0f0:3.0f0...]]) for (region_idx,region) ∈ enumerate(settings["regions"])]
    defVar(fh,"rsl_trajectory",trajectory,("region","years","percentiles"),deflatelevel=5, attrib = Dict("description" => "Estimated regional sea-level trajectory, relative to year 2000.", "units" => "mm"))

    # Write trends and accelerations
    trends = zeros(Float32,length(settings["regions"]),3) .* NaN32
    accels = zeros(Float32,length(settings["regions"]),3) .* NaN32
    [trends[region_idx,:] = region_obs[region]["trend"] for (region_idx,region) ∈ enumerate(settings["regions"])]
    [accels[region_idx,:] = region_obs[region]["accel"] for (region_idx,region) ∈ enumerate(settings["regions"])]
    defVar(fh,"rsl_trend",trends,("region","percentiles"),deflatelevel=5, attrib = Dict("description" => "Linear trend of estimated trajectory", "units" => "mm yr-1"))
    defVar(fh,"rsl_accel",accels,("region","percentiles"),deflatelevel=5, attrib = Dict("description" => "Acceleration (half quadratic) of estimated trajectory", "units" => "mm yr-2"))

    # Write projections
    scn_array = zeros(Float32,length(settings["regions"]),length(settings["years"]),3) .* NaN32
    for prc in settings["processes"]
        for scn ∈ settings["NCA5_scenarios"]
            for (region_idx,region) ∈ enumerate(settings["regions"])
                scn_array[region_idx,:,:] = (LinearInterpolation((NCA5_regional["years"],[1.0f0:3.0f0...]),NCA5_regional[scn][prc][region],extrapolation_bc=NaN32)[settings["years"],[1.0f0:3.0f0...]])
            end
            defVar(fh,"rsl_"*prc*"_"*scn,scn_array,("region","years","percentiles",),deflatelevel=5, attrib = Dict("description" => "Projected regional sea level for scenario "*scn*" and process "*prc*", relative to year 2000.", "units" => "mm"))
        end
    end
    close(fh)
end

end


text:
#######
# Demo setup

cool(x) = x + 1
cool(x, y) = x + y + 1

# a rule we define so we can test rules
dummy_identity(x) = x

@scalar_rule(dummy_identity(x), One())

#######

_second(t) = Base.tuple_type_head(Base.tuple_type_tail(t))

@testset "frule and rrule" begin
    @test frule(cool, 1) === nothing
    @test frule(cool, 1; iscool=true) === nothing
    @test rrule(cool, 1) === nothing
    @test rrule(cool, 1; iscool=true) === nothing

    ChainRulesCore.@scalar_rule(Main.cool(x), one(x))
    @test hasmethod(rrule, Tuple{typeof(cool),Number})
    ChainRulesCore.@scalar_rule(Main.cool(x::String), "wow such dfdx")
    @test hasmethod(rrule, Tuple{typeof(cool),String})
    # Ensure those are the *only* methods that have been defined
    cool_methods = Set(m.sig for m in methods(rrule) if _second(m.sig) == typeof(cool))
    only_methods = Set([Tuple{typeof(rrule),typeof(cool),Number},
                        Tuple{typeof(rrule),typeof(cool),String}])
    @test cool_methods == only_methods

    frx, fr = frule(cool, 1)
    @test frx == 2
    @test fr(1) == 1
    rrx, rr = rrule(cool, 1)
    @test rrx == 2
    @test rr(1) == 1
end


text:
using BenchmarkTools
using StreamModels
using DataStreams
using DataFrames

using StreamModels: @formula

n = 100_000
source = DataFrame(a = collect(1:n),
                   b = rand(n),
                   c = repeat(["a", "b"], inner=n>>1))

f = @formula( ~ a)

@benchmark StreamModels.parse!($f)

b = @benchmarkable modelmatrix(source, f)

# tune!(b)
# BenchmarkTools.save("params.jld", "b", params(b))

loadparams!(b, BenchmarkTools.load("params.jld", "b"))
source2 = copy(source)
pool!(source2, :c)
ff = DataFrames.Formula(f.lhs, f.rhs)

b2 = @benchmark ModelMatrix(ModelFrame($ff, $source2))

run(b)


text:
printstyled("  time test with synthetic timeseries data\n", color=:light_green)

ts = round(Int64, sμ*d2u(DateTime("2020-03-01T00:00:00")))
fs = 100.0
gain = 32.0
id = "XX.STA..BHZ"
loc = GeoLoc( lat = 45.560121,
              lon = -122.617068,
              el = 53.04 )
Δ = 10000
files = ["1.sac", "2.sac", "3.sac", "4.sac", "5.sac", "6.sac", "7.sac"]
tos   = [    0.0,    10.0,   20.01,   20.02,    30.0,    40.0,    60.0]
nx    = [   1000,     999,       1,     100,       1,    1000,    1000]

C = SeisChannel(id = id,
                gain = gain,
                loc = loc,
                fs = fs )
S0 = SeisData()

printstyled("    creating files\n", color=:light_green)
for i = 1:length(files)
  fname = files[i]
  C.x = randn(Float32, nx[i])
  C.t = [1 ts+round(Int64, tos[i]*1000000); nx[i] 0]
  push!(S0, C)
  writesac(C, fname=fname)
end

# Read back in
printstyled("    reading to test all channel-extension cases\n", color=:light_green)
S = SeisData()
read_data!(S, "sac", files[1])
xi = cumsum(nx)
δt = round.(Int64, tos.*sμ) .+ (nx.-1).*Δ

# Expect:
# no gaps, precise read
@test S.t[1] == [1 ts; nx[1] 0]
@test endtime(S.t[1], S.fs[1]) == Δ*(nx[1]-1)+ts
@test S0.x[1] == S.x[1]
@test S0.gain[1] == S.gain[1]

# append file 2 to create case 1
read_data!(S, "sac", files[2])
@test S.t[1] == [1 ts; xi[2] 0]
@test S0.x[2] == S.x[1][nx[1]+1:xi[2]]

# append file 3 to create case 4
read_data!(S, "sac", files[3])
@test S.t[1] == [1 ts; xi[3] 20000]
@test S0.x[3] == S.x[1][xi[2]+1:xi[3]]

# append file 4 to create case 2
read_data!(S, "sac", files[4])
@test S.t[1] == [1 ts; xi[3] 20000; xi[4] 0]
@test S0.x[4] == S.x[1][xi[3]+1:xi[4]]

# append file 5 to create case 3
δt0 = round(Int64, sμ*(tos[5]-tos[4]-nx[4]/fs)) # I canceled a +Δ and a -Δ
read_data!(S, "sac", files[5])
@test S.t[1] == [1 ts; xi[3] 20000; xi[5] δt0]
@test S0.x[5] == S.x[1][xi[4]+1:xi[5]]

# append file 6 to create case 6
t_old = copy(S.t[1])
δt0 = round(Int64, sμ*(tos[6]-tos[5]-nx[5]/fs)) # I canceled a +Δ and a -Δ
read_data!(S, "sac", files[6])
@test S.t[1] == vcat(t_old, [xi[5]+1 δt0; xi[6] 0])
@test S0.x[6] == S.x[1][xi[5]+1:xi[6]]

# append file 7 to create case 5
t_old = copy(S.t[1])
read_data!(S, "sac", files[7])
δt0 = round(Int64, sμ*(tos[7]-tos[6]-nx[7]/fs))
@test S.t[1] == vcat(t_old[1:end-1,1:2], [t_old[end,1]+1 δt0], [xi[7] 0])

# finally, did we get it all right?
printstyled("    final checks\n", color=:light_green)
T = t_collapse(vcat([t_expand(i, fs) for i in S0.t]...), fs)
X = vcat([i for i in S0.x]...)
@test T == S.t[1]
@test X == S.x[1]
@test gain == S.gain[1]
@test fs == S.fs[1]
@test id == S.id[1]
@test ≈(loc.lat, S.loc[1].lat, rtol=2eps(Float32))
@test ≈(loc.lon, S.loc[1].lon, rtol=2eps(Float32))
@test ≈(loc.el, S.loc[1].el, rtol=2eps(Float32))

printstyled("    cleanup\n", color=:light_green)
for i in files
  safe_rm(i)
end


text:
using Documenter, ImageSegmentationEvaluation

makedocs(
  format = :html,
    sitename = "Image Segmentation Evaluation",
    pages = [
        "Index"=> "index.md",
        "Subsection" => [
            "Supervised methods" => "supervised.md",
            "Unsupervised methods" => "unsupervised.md"
        ]
    ]
)


text:
struct IdSet{T} <: AbstractSet{T}
  dict::ObjectIdDict
  IdSet{T}() where T = new(ObjectIdDict())
end

Base.eltype{T}(::IdSet{T}) = T

IdSet() = IdSet{Any}()

Base.push!{T}(s::IdSet{T}, x::T) = (s.dict[x] = nothing; s)
Base.delete!{T}(s::IdSet{T}, x::T) = (delete!(s.dict, x); s)
Base.in(x, s::IdSet) = haskey(s.dict, x)

(::Type{IdSet{T}}){T}(xs) = push!(IdSet{T}(), xs...)

IdSet(xs) = IdSet{eltype(xs)}(xs)

Base.collect(s::IdSet) = Base.collect(keys(s.dict))
Base.similar(s::IdSet, T::Type) = IdSet{T}()

@forward IdSet.dict Base.length

Base.start(s::IdSet) = start(keys(s.dict))
Base.next(s::IdSet, st) = next(keys(s.dict), st)
Base.done(s::IdSet, st) = done(keys(s.dict), st)


text:
using Test
using LIKWID
using LinearAlgebra

# julia must have been started with `likwid-perfctr -C ... -g ... -m`
@assert Marker.isactive()

# init
@test isnothing(Marker.init())

LIKWID.PerfMon.init() # required for getregion below

N = 100_000_000
a = 3.141f0
z = zeros(Float32, N)
x = rand(Float32, N)
y = rand(Float32, N)

saxpy!(z,a,x,y) = z .= a .* x .+ y

# regular workflow
saxpy!(z,a,x,y)
@test Marker.startregion("saxpy!")
saxpy!(z,a,x,y)
@test Marker.stopregion("saxpy!")

# nextgroup
@test isnothing(Marker.nextgroup())

# registerregion
A = rand(100,100)
B = rand(100,100)
@test Marker.registerregion("mul")
@test Marker.startregion("mul")
for _ in 1:10
    A * B
end
@test Marker.stopregion("mul")

# getregion
x = Marker.getregion("mul")
@test typeof(x) == Tuple{Int32, Vector{Float64}, Float64, Int32}

# resetregion
@test Marker.resetregion("mul")
y = Marker.getregion("mul")
@test x != y
@test y[3] == 0
@test y[4] == 0

# close
@test isnothing(Marker.close())


text:
using Unitful
using PlutoUI
using PhysicalConstants.CODATA2018
using Plots
using Printf

k_B = BoltzmannConstant
c_0 = SpeedOfLightInVacuum
h = PlanckConstant

function planck_spectrum(lambdas, T)
	res = []
	for λ in lambdas
		push!(
			res,
			ustrip(
				uconvert(
					u"W/m^3", 2*pi*h*c_0^2 / λ.^5 / (exp(h*c_0/λ/k_B/T) - 1))
				)
			)
	end
	return res
end

function add_visible_light(lambda_unit, height, alpha)
	bar!(
		[ustrip(uconvert(lambda_unit, 400u"nm"))],
		[height],
		fillcolor=[:violet],
		bar_width=ustrip(uconvert(lambda_unit, 40u"nm")),
		fillalpha=alpha,
		label="",
	)
	bar!(
		[ustrip(uconvert(lambda_unit, 455u"nm"))],
		[height],
		fillcolor=[:blue],
		bar_width=ustrip(uconvert(lambda_unit, 70u"nm")),
		fillalpha=alpha,
		label="",
	)
	bar!(
		[ustrip(uconvert(lambda_unit, 532.5u"nm"))],
		[height],
		fillcolor=[:green],
		bar_width=ustrip(uconvert(lambda_unit, 85u"nm")),
		fillalpha=alpha,
		label="",
	)
	bar!(
		[ustrip(uconvert(lambda_unit, 580u"nm"))],
		[height],
		fillcolor=[:yellow],
		bar_width=ustrip(uconvert(lambda_unit, 10u"nm")),
		fillalpha=alpha,
		label="",
	)
	bar!(
		[ustrip(uconvert(lambda_unit, 617.5u"nm"))],
		[height],
		fillcolor=[:orange],
		bar_width=ustrip(uconvert(lambda_unit, 65u"nm")),
		fillalpha=alpha,
		label="",
	)
	bar!(
		[ustrip(uconvert(lambda_unit, 700u"nm"))],
		[height],
		fillcolor=[:red],
		bar_width=ustrip(uconvert(lambda_unit, 100u"nm")),
		fillalpha=alpha,
		label="",
	)
end

function plot_stuff(T)
    ls3 = 10 .^ (range(-2, stop=4, length=1000))
    y = planck_spectrum(ls3.*1u"nm", T*1u"K")
    add_visible_light(u"nm", maximum(y), 0.7)
    plot!(
        ls3,
        #ls3,
        y,
        #xscale=:log10,
        #yscale=:log10,
        xlim=[100, 1500],
        xlabel="Wellenlänge [nm]",
        ylabel="Intensität [W/m³]",
        label="$T K",
        color=:black,
        grid=true,
        lw=3
    )
end

for T in [3000, 6000, 12000]
    pgfplotsx()
    pl = plot()
    plot_stuff(T)
    savefig(pl, "planck_$T.tex")

    gr()
    pl = plot()
    plot_stuff(T)
    savefig(pl, "planck_$T.png")

    pl = plot()
    plot_stuff(T)
    savefig(pl, "planck_$T.svg")
end


text:
export learn_chow_liu_tree_circuit, learn_vtree_from_clt, compile_sdd_from_clt
using Graphs: outneighbors
using MetaGraphs: get_prop

"""
Learning from data a structured-decomposable circuit with several structure learning algorithms
"""
function learn_chow_liu_tree_circuit(data;
        pseudocount = 1.0, algo_kwargs=(α=1.0, clt_root="graph_center"), vtree_kwargs=(vtree_mode="balanced",))
    clt = learn_chow_liu_tree(data; algo_kwargs...)
    vtree = learn_vtree_from_clt(clt; vtree_kwargs...)
    lc = compile_sdd_from_clt(clt, vtree)
    pc = ProbCircuit(lc)
    estimate_parameters!(pc, data; pseudocount=pseudocount)
    pc, vtree
end

#############
# Learn PlainVtree from CLT
#############

"
Learn a vtree from clt,
with strategy (close to) `linear` or `balanced`
"
function learn_vtree_from_clt(clt::CLT; vtree_mode::String)::PlainVtree
    roots = [i for (i, x) in enumerate(parent_vector(clt)) if x == 0]
    rootnode = construct_children(Var.(roots), clt, vtree_mode)

    return rootnode
end

function construct_node(v::Var, clt::CLT, strategy::String)::PlainVtree
    children = Var.(outneighbors(clt, v))
    if isempty(children) # leaf node
        return PlainVtreeLeafNode(v)
    else
        right = construct_children(children, clt, strategy)
        return add_parent(v, right)
    end
end

function construct_children(children::Vector{Var}, clt::CLT, strategy::String)::PlainVtree
    sorted_vars = sort(collect(children))
    children_nodes = Vector{PlainVtree}()
    foreach(x -> push!(children_nodes, construct_node(x, clt, strategy)), sorted_vars)

    if strategy == "linear"
        construct_children_linear(children_nodes, clt)
    elseif strategy == "balanced"
        construct_children_balanced(children_nodes, clt)
    else
        throw("Unknown type of strategy")
    end
end

function construct_children_linear(children_nodes::Vector{PlainVtree}, clt::CLT)::PlainVtree
    children_nodes = Iterators.Stateful(reverse(children_nodes))

    right = popfirst!(children_nodes)
    for left in children_nodes
        right = PlainVtreeInnerNode(left, right)
    end
    return right
end

function construct_children_balanced(children_nodes::Vector{PlainVtree}, clt::CLT)::PlainVtree
    if length(children_nodes) == 1
        return children_nodes[1]
    elseif length(children_nodes) == 2
        return PlainVtreeInnerNode(children_nodes[1], children_nodes[2])
    else
        len = trunc(Int64, length(children_nodes) / 2)
        left = construct_children_balanced(children_nodes[1 : len], clt)
        right = construct_children_balanced(children_nodes[len + 1 : end], clt)
        return PlainVtreeInnerNode(left, right)
    end
end

function add_parent(parent::Var, children::PlainVtree)
    return PlainVtreeInnerNode(PlainVtreeLeafNode(parent), children)
end

#####################
# Compile PSDD from CLT and vtree
#####################

"Compile a psdd circuit from clt and vtree"
function compile_sdd_from_clt(clt::CLT, vtree::PlainVtree)::PlainStructLogicCircuit

    parent_clt = Var.(parent_vector(clt))
    v2p = Dict{PlainVtree, Vector{PlainStructLogicCircuit}}()

    function add_mapping!(v::PlainVtree, circuits)
        if !haskey(v2p, v); v2p[v] = Vector{PlainStructLogicCircuit}(); end
        foreach(c -> if !(c in v2p[v]) push!(v2p[v], c);end, circuits)
    end

    # compile vtree leaf node to terminal/true node
    function compile_from_vtree_node(v::PlainVtreeLeafNode)
        var = v.var
        children = Var.(outneighbors(clt, var))
        cpt = get_prop(clt, var, :cpt)
        parent = parent_clt[var]
        if isequal(children, [])
            circuit = compile_true_nodes(var, v; num=length(cpt) ÷ 2)
        else
            circuit = compile_canonical_literals(var, v)
        end
        add_mapping!(v, circuit)
        nothing
    end

    # compile to decision node
    function compile_from_vtree_node(v::PlainVtreeInnerNode)
        left_var = left_most_descendent(v.left).var
        right_var = left_most_descendent(v.right).var
        left_circuit = v2p[v.left]
        right_circuit = v2p[v.right]

        if parent_clt[left_var] == parent_clt[right_var] # two nodes are independent, compile to seperate decision nodes
            circuit = [compile_decision_node([l], [r], v) for (l, r) in zip(left_circuit, right_circuit)]
        elseif left_var == parent_clt[right_var] # conditioned on left
            cpt = get_prop(clt, left_var, :cpt)
            circuit = compile_decision_nodes(left_circuit, right_circuit, v; num=length(cpt) ÷ 2)
        else
            throw("PlainVtree are not learned from the same CLT")
        end
        add_mapping!(v, circuit)
        nothing
    end

    foreach(compile_from_vtree_node, vtree)

    v2p[vtree][end]
end

#####################
# Construct circuit node
#####################
"Construct decision nodes given `primes` and `subs`"
function compile_decision_node(primes::Vector{<:PlainStructLogicCircuit}, subs::Vector{<:PlainStructLogicCircuit}, vtree::PlainVtreeInnerNode)
    elements = [conjoin(prime, sub; use_vtree=vtree) for (prime, sub) in zip(primes, subs)]
    return disjoin(elements; use_vtree=vtree)
end

"Construct literal nodes given variable `var`"
function compile_canonical_literals(var::Var, vtree::PlainVtreeLeafNode)
    return [PlainStructLiteralNode( var2lit(var), vtree), PlainStructLiteralNode(-var2lit(var), vtree)]
end

"Construct true nodes given variable `var`"
function compile_true_nodes(var::Var, vtree::PlainVtreeLeafNode; num)
    pos, neg = compile_canonical_literals(var, vtree)
    return [disjoin([pos, neg]; use_vtree = vtree) for _ in 1 : num]
end

"Construct decision nodes conditiond on different distribution"
function compile_decision_nodes(primes::Vector{<:PlainStructLogicCircuit}, subs::Vector{<:PlainStructLogicCircuit}, vtree::PlainVtreeInnerNode; num)
    return [compile_decision_node(primes, subs, vtree) for _ in 1 : num]
end