text:
export VarGamma

import LinearAlgebra: norm
import SpecialFunctions: gamma, digamma, polygamma
include(joinpath(@__DIR__, "../prob_operations.jl"))
include(joinpath(@__DIR__, "../graph_operations.jl"))
include(joinpath(@__DIR__, "../util.jl"))

mutable struct VarGamma
    """Variable with a Gamma recognition distribution"""

    # Factor graph properties
    id::Symbol
    time::Union{Integer,Nothing}
    block::Bool

    # Recognition distribution parameters
    marginal::Gamma{Float64}
    free_energy::Float64
    grad_free_energy::Tuple{Float64,Float64}

    function VarGamma(id::Symbol;
                      marginal::Gamma{Float64}=Gamma(1.,1e-12),
                      grad_free_energy::Tuple{Float64,Float64}=(1e12,1e12),
                      free_energy::Float64=1e12,
                      time::Union{Integer,Nothing}=0,
                      block::Bool=false)

     # Construct instance
     return new(id, time, block, marginal, free_energy, grad_free_energy)
 end
end

function params(variable::VarGamma)
    "Parameters of current belief"
    return params(variable.marginal)
end

function mean(variable::VarGamma)
    "Mean of current belief"
    return mean(variable.marginal)
end

function var(variable::VarGamma)
    "Variance of current belief"
    return var(variable.marginal)
end

function moments(variable::VarGamma)
    "First two moments of current belief"
    return mean(variable), var(variable)
end

function update!(graph::MetaGraph, variable::VarGamma)
    "Update recognition distribution as the product of messages"

    # List of incoming messages
    incoming = Any[]

    # Check each connected edge
    for edge in edges(graph, graph[variable.id, :id])

		# Mesages
		message_f2v = get_prop(graph, edge, :message_factor2var)

        # Check for incoming
        push!(incoming, message_f2v)
    end

    # Update marginal
    variable.marginal = prod(incoming)
end

function marginal(variable::VarGamma)
    "Marginal over variable"
    return variable.marginal
end

function entropy(variable::VarGamma)
    "Entropy of Gamma distribution"

    # Parameters
    α = shape(variable.marginal)
    θ = scale(variable.marginal)

    # Entropy of a univariate Gamma
    return α + log(θ) + log(gamma(α)) + (1-α)*digamma(α)
end

function grad_entropy(variable::VarGamma)
    "Gradient of entropy of Gamma evaluated for supplied parameters"

	# Parameters
    α = shape(variable.marginal)
    θ = scale(variable.marginal)

    # Partial derivative with respect to shape
    partial_shape = 1 + (1 - α)*polygamma(1, α)

    # Partial derivative with respect to scale
    partial_scale = 1 / θ

    # Return tuple of partial derivatives
    return (partial_shape, partial_scale)
end

function free_energy(graph::MetaGraph, variable::VarGamma)
    "Compute local free energy of variable and connecting factor nodes"

    # Initialize energy
    U = 0.

	# Iterate over edges connected to variable node
	for edge in edges(graph, graph[variable.id, :id])

		# Extract message from factor to variable
		message_f2v = get_prop(graph, edge, :message_factor2var)

		# Compute the expectation E_qx [log ν(x)]
		U += expectation(variable.marginal, message_f2v)
	end

	# Entropy
	H = entropy(variable)

	return U - H
end

function grad_free_energy(graph::MetaGraph, variable::VarGamma)
    "Compute gradient of free energy w.r.t. recognition distribution parameters."

	# Initialize energy
	partial_α = 0.
	partial_θ = 0.

	# Iterate over edges connected to variable node
	for edge in edges(graph, graph[variable.id, :id])

		# Extract message from factor to variable
		message_f2v = get_prop(graph, edge, :message_factor2var)

		# Compute the gradient of expectation E_qx [log ν(x)] w.r.t. q(x) params
		partials = grad_expectation(variable.marginal, message_f2v)

		# Update partial parameters
		partial_α += partials[1]
		partial_θ += partials[2]
	end

	# Gradient of entropy
	partials = grad_entropy(variable)

	# Update partial parameters
	partial_α += partials[1]
	partial_θ += partials[2]

	return (partial_α, partial_θ)
end

function act!(graph::MetaGraph, variable::VarGamma)
    "Pass belief to edges."

    # Check each connected edge
    for edge in edges(graph, graph[variable.id, :id])

        # Check for incoming
		set_prop!(graph, edge, :message_var2factor, variable.marginal)
		set_prop!(graph, edge, :∇free_energy, norm(variable.grad_free_energy))
    end
end

function react!(graph::MetaGraph, variable::VarGamma)
    "React to incoming messages"

    # Check whether variable is blocked
    if !variable.block

        # Update variational distribution
        update!(graph, variable)

        # Compute free energy after update
        variable.free_energy =free_energy(graph, variable)

        # Compute gradient of free energy after update
        variable.grad_free_energy = grad_free_energy(graph, variable)

        # Message from edge to nodes
        act!(graph, variable)
    end
end

max_stars_repo_name: wmkouw/schedulefreeMP

max_issues_repo_name: wmkouw/schedulefreeMP

max_forks_repo_name: wmkouw/schedulefreeMP


text:
# Juliaを使った株価の取得と通知

using DataFrames, Dates, CSV
include("GetStockPrice.jl")

uniqueindex(x) = findfirst.(isequal.(unique(x)), [x]) # 重複しない行のindexを計算

# 株価データの初期化．新たに銘柄追加するときだけこれを実行する
function initdailystockprice(stock)
    range = "7d"    # 取得データの期間
    interval = "1m" # 取得データの間隔
    df = getstockprice(stock, range, interval)

    # CSV書き込み前に時刻のフォーマットを変更しておく
    dtfmt = "yyyy/mm/dd HH:MM:SS"
    df."timestamp" = Dates.format.(df."timestamp", dtfmt)

    # CSVへの書き込み
    println("write stock price: " * stock * ".csv")
    CSV.write("./data/" * stock * ".csv", df)
end

# 1日分の1分足データを取得し，CSVにマージする
function savedailystockprice(stock)
    range = "1d"    # 取得データの期間
    interval = "1m" # 取得データの間隔
    today = getstockprice(stock, range, interval)
    dtfmt = "yyyy/mm/dd HH:MM:SS"

    # CSV読み込み
    past = DataFrame(CSV.File("./data/" * stock * ".csv", header = true))
    past."timestamp" = DateTime.(past."timestamp", dtfmt)

    # データを結合
    all = vcat(past, today)

    # 重複データを削除
    all = all[uniqueindex(all.timestamp), :]

    # 日付順でソート
    all = sort(all)

    # CSV書き込み前に時刻のフォーマットを変更しておく
    all."timestamp" = Dates.format.(all."timestamp", dtfmt)

    # CSVへの書き込み
    println("write stock price: " * stock * ".csv")
    CSV.write("./data/" * stock * ".csv", all)
end

function main()
    # 取得対象の銘柄の指定．東証の株価は".T"をつける．
    stocks = ["1557.T" "SPY" "SPXL" "SPXS" "TECL" "CURE" "SOXL" "DRN" "TNA" "EDC"]
    try
        savedailystockprice.(stocks)
    catch err
        println(err)
    end
end

main()

ext: jl


text:
"""
    ART.jl

Description:
    Includes all of the unsupervised ART modules definitions.
"""

include("common.jl")    # train!, classify
include("FuzzyART.jl")  # FuzzyART
include("DDVFA.jl")     # DDVFA
include("DVFA.jl")      # DVFA


text:
"""
Module Mouse

Provides functions for controlling the mouse by voice.

# Functions

###### Button control
- [`click_left`](@ref)
- [`click_middle`](@ref)
- [`click_right`](@ref)
- [`press_left`](@ref)
- [`release_left`](@ref)
- [`click_double`](@ref)
- [`click_triple`](@ref)

To see a description of a function type `?<functionname>`.

See also: [`Keyboard`](@ref)
"""
module Mouse

using PyCall
import ..JustSayIt: pyimport_pip, controller, set_controller

## PYTHON MODULES

const Pynput = PyNULL()

function __init__()
    ENV["PYTHON"] = ""                                              # Force PyCall to use Conda.jl
    copy!(Pynput, pyimport_pip("pynput"))
    set_controller("mouse", Pynput.mouse.Controller())
end

## FUNCTIONS

"Click left mouse button."
click_left(;count::Integer=1) = controller("mouse").click(Pynput.mouse.Button.left, count)

"Click middle mouse button."
click_middle() = controller("mouse").click(Pynput.mouse.Button.middle)

"Click right mouse button."
click_right() = controller("mouse").click(Pynput.mouse.Button.right)

"Press and hold left mouse button."
press_left() = controller("mouse").press(Pynput.mouse.Button.left)

"Release left mouse button."
release_left() = controller("mouse").release(Pynput.mouse.Button.left)

"Doubleclick left mouse button."
click_double() = click_left(;count=2)

"Trippleclick left mouse button."
click_triple() = click_left(;count=3)

end # module Mouse


text:
# Load the module
using ROBDD

# Define your boolean variables
# (must be a vector of symbols)
ordered_variables = [:x1, :x2, :x3]

# Initiate a ROBDDTable. This datastructure
# encodes context information (the set of variables
# and their order); and *all* of the ROBDDs you
# construct over them. Each entry is a node
# in a ROBDD.
my_bdd_table = ROBDDTable(ordered_variables)

# At first it only has two entries.
# one for "false" and one for "true".
println(my_bdd_table)

# output:
#
# (0) false
# (1) true

# Construct the Boolean expression you're modeling
#my_expr = :(x1 | (x2 & !x2) | x3)
my_expr = :(x1 | x2 | x3)

# Build a ROBDD from the expression.
# The `build_robdd` function returns an
# integer index. The index points to a particular
# entry of the ROBDDTable.
idx = build_robdd(my_bdd_table, my_expr)

# During construction, several new
# entries were added to your table.
# `build_robdd` tells you which entry
# contains the root node of your BDD.
println(my_bdd_table)
println("INDEX: ", idx)

# output:
#
# (0) false
# (1) true
# (2) x1 0 1
# (3) x2 0 1
# (4) x2 1 0
# (5) x3 0 1
# (6) x1 5 1
# INDEX: 6

# We provide an interface for
# (1) getting an arbitrary SAT solution and
# (2) enumerating all SAT solutions
#     (modulo the unassigned variables--these
#      may take any value.)
println("ANYSAT: ", anysat(my_bdd_table, idx))
println("ALLSAT:")
for sat in allsat(my_bdd_table, idx)
    println("\t",sat)
end

# output:
#
# ANYSAT: Dict{Symbol,Bool}(:x2 => 0,:x3 => 1,:x1 => 0)
# ALLSAT:
# 	Dict{Symbol,Bool}(:x2 => 0,:x3 => 1,:x1 => 0)
# 	Dict{Symbol,Bool}(:x2 => 1,:x1 => 0)
# 	Dict{Symbol,Bool}(:x1 => 1)

# We can "restrict" a ROBDD by assigning
# concrete values to a subset of its variables.
constraints = Dict(:x3 => false)
restricted_idx = restrict(my_bdd_table, idx, constraints)
println("RESTRICTIONS: ", constraints)
println(my_bdd_table)
println("RESTRICTED IDX: ", restricted_idx)

# output:
#
# RESTRICTIONS: Dict{Symbol,Bool}(:x3 => 0)
# (0) false
# (1) true
# (2) x1 0 1
# (3) x2 0 1
# (4) x2 1 0
# (5) x3 0 1
# (6) x1 5 1
# RESTRICTED IDX: 2

# During construction the ROBDDTable accumulates
# many entries. However, some of them are not actually used by the
# ROBDD of interest. We provide a function for "cleaning"
# the table, leaving only the descendants of a specified BDD.
new_table, clean_idx = clean_table(my_bdd_table, restricted_idx)
println("CLEANED TABLE:")
println(new_table)
println("CLEAN IDX: ", clean_idx)

# output:
#
# CLEANED TABLE:
# (0) false
# (1) true
# (2) x1 0 1
# CLEAN IDX: 2

# You should only clean a ROBDDTable when you don't intend to
# construct any more ROBDDs. Otherwise, those extra ROBDDs
# may be useful for constructing other ROBDDs.

# ***MORE FEATURES TO COME***


text:
################################################################################
#
#	ABSTRACT TYPE
#
#   BondSpinHamiltonianHeisenberg <: AbstractBondHamiltonian{L,N}
#   --> L is the label type of bonds
#   --> N is the dimension the bond term matrix (NxN matrix)
#
#   FILE CONTAINS
#       - abstract type definition
#       - interface definition
#       - generator function
#
################################################################################

################################################################################
#
#   ABSTRACT TYPE DEFINITION
#
################################################################################
mutable struct BondSpinHamiltonianHeisenberg{N,L,NS} <: AbstractBondHamiltonian{L,NS}

    # couplings
    J :: Vector{Float64}

    # bond labels of these couplings
    J_bonds :: Vector{Vector{L}}

end

# export the type
export BondSpinHamiltonianHeisenberg

################################################################################
#
#	INTERFACING / ACCESSING BOND HAMILTONIANS
#	(functions have to be overwritten by concrete types)
#
################################################################################

# get bond term of Hamiltonian
function bondterm(
            h :: BondSpinHamiltonianHeisenberg{N,L,NS},
            b :: AbstractBond{L,NB}
        ) :: Matrix{Complex{Float64}} where {N,L,NB,NS}

    # get the bond label
    l = label(b)

    # check for all neighbors if the label is in one of the lists
    for i in 1:N
        if l in h.J_bonds[i]
            return zeros(NS,NS) + (h.J[i] * I)
        end
    end

    # if no bond found, return 0
    return zeros(NS,NS)
end

# export the function
export bondterm

# FUNCTIONS TO FIND BOND LABELS

# FALLBACK ANY LABEL TYPE
# if specific function for label type L not found, denote all bonds as nearest neighbors
function getBondLabelsHeisenberg(
        N :: Val{1},
        all_couplings :: Vector{L}
    ) :: Vector{L} where {L}

    # all couplings that either have 1 or no number in them
    string_couplings = string.(all_couplings)
    relevant = Bool[false for s in string_couplings]
    for i in 1:length(relevant)
        relevant[i] = occursin("1",string_couplings[i]) || sum([isdigit(c) for c in string_couplings[i]])::Int64==0
    end
    # return only the relevant couplings
    return L[all_couplings[i] for i in 1:length(all_couplings) if relevant[i]]
end
function getBondLabelsHeisenberg(
        N :: Val{NB},
        all_couplings :: Vector{L}
    ) :: Vector{L} where {NB,L}

    # all couplings that either have 1 or no number in them
    string_couplings = string.(all_couplings)
    relevant = Bool[false for s in string_couplings]
    for i in 1:length(relevant)
        relevant[i] = occursin(string(NB),string_couplings[i])
    end
    # return only the relevant couplings
    return L[all_couplings[i] for i in 1:length(all_couplings) if relevant[i]]
end

# CONSTRUCTION FUNCTION
function getSpinHamiltonianHeisenberg(
            unitcell :: U,
            N :: Int64 = 1,
            NS :: Int64 = 3
        ) :: BondSpinHamiltonianHeisenberg{N,L,NS} where {L,NB,S,B<:AbstractBond{L,NB},U<:AbstractUnitcell{S,B}}

    # obtain all couplings
    couplings = unique!(label.(bonds(unitcell)))

    # create lists and initialize default couplings
    J = [1 for i in 1:N]
    J_bonds = [getBondLabelsHeisenberg(Val(i),couplings) for i in 1:N]

    # create and return a new object
    return BondSpinHamiltonianHeisenberg{N,L,NS}(J,J_bonds)
end

# export the construction function
export getSpinHamiltonianHeisenberg

# CUSTOM SHOWING

# all labels
function show(io::IO, h::H) where {L,N,NS,H<:BondSpinHamiltonianHeisenberg{N,L,NS}}
    print(io, "Heisenberg (Bond) Hamiltonian of O("*string(NS)*") spins up to "*string(N)*" neighbors. Couplings strength and accepted labels are:")
    for i in 1:length(h.J)
        print(io, "\n--> J_"*string(i)*"="*string(h.J[i])*", labels="*string(h.J_bonds[i]))
    end
end


text:
#!/usr/bin/env julia

println("Loading libraries")
using Neurthino
using Unitful
using Plots
using ProgressMeter

function main()
    println("Initialising oscillation parameters")
    osc = OscillationParameters(3);
    mixingangle!(osc, 1 => 2, 0.58);
    mixingangle!(osc, 1 => 3, 0.15);
    mixingangle!(osc, 2 => 3, 0.737);
    cpphase!(osc, 1 => 3, 5.34);
    masssquareddiff!(osc, 2 => 3, -2.457e-3);
    masssquareddiff!(osc, 1 => 3, -2.457e-3-7.5e-5);
    masssquareddiff!(osc, 1 => 2, -7.5e-5);

    U = PMNSMatrix(osc)
    H = Hamiltonian(osc)

    coszenith = range(-1; stop=0, length=100)
    zenith = acos.(coszenith)
    energies = 10 .^(range(0;stop=2,length=100))

    data_points = zeros((3,3,length(zenith), length(energies)))

    @showprogress 1 "Calculating transition probabilities for different zenith angles" for (i,z) in enumerate(zenith)
        sec, dens = Neurthino.prempath(z, 2.5, samples=50)
        P = Neurthino.mattertransprob(osc, energies, dens, sec)
        for (k,l) in Iterators.product(fill(1:3, 2)...)
            data_points[k,l,i,:] = map(x->x[k,l], P)
        end
    end

    plot_arr = Array{Plots.Plot{Plots.GRBackend},1}()
    for (k,l) in Iterators.product(fill(1:3, 2)...)
        p = heatmap(data_points[k,l,:,:],  yticks=(1:20:100, coszenith[1:10:100]), xticks=(1:40:200, energies[1:10:50]), clim=(0,1))
        push!(plot_arr, p)
    end
    files = Vector{AbstractString}()
    @showprogress 1 "Creating plots" for (idx, p) in enumerate(plot_arr)
        f = "$idx.pdf"
        push!(files, f)
        savefig(p, f)
    end
end

main()


text:
function calc_discontinuity(p, v, d, t, annual_year, r)
    # setting the year for entry in lists.
    yr = annual_year - 2015 + 1 # + 1 because of 1-based indexing in Julia

    # global variables only have to be done once. [r goes from 1 to 8]
    if isequal(r, 1)

        v.idis_lossfromdisc_ann[yr] = max(0, p.rt_g_globaltemperature_ann[yr] - p.tdis_tolerabilitydisc) # global temperature minus 1.5 (default)

        if is_first(t)
            if v.idis_lossfromdisc_ann[yr] * (p.pdis_probability / 100) > p.rand_discontinuity # (temp - 1.5)*0.2 >? .5 (default)
                v.occurdis_occurrencedummy_ann[yr] = 1
            else
                v.occurdis_occurrencedummy_ann[yr] = 0
            end
            v.expfdis_discdecay_ann[yr] = exp(-(p.y_year[t] - p.y_year_0) / p.distau_discontinuityexponent) # left same as timestep version
        else
            if v.idis_lossfromdisc_ann[yr] * (p.pdis_probability / 100) > p.rand_discontinuity
                v.occurdis_occurrencedummy_ann[yr] = 1
            elseif v.occurdis_occurrencedummy_ann[yr - 1] == 1
                v.occurdis_occurrencedummy_ann[yr] = 1
            else
                v.occurdis_occurrencedummy_ann[yr] = 0
            end
            v.expfdis_discdecay_ann[yr] = exp(-(p.y_year[t] - p.y_year[t - 1]) / p.distau_discontinuityexponent) # left same as timestep version
        end
    end

    # in growth effects with feedback, `p.rgdp_per_cap_NonMarketRemainGDP_ann[yr,r]` goes below zero, leading to issues here.
    v.igdpeqdis_eqdiscimpact_ann[yr,r] = v.irefeqdis_eqdiscimpact[r] * (p.rgdp_per_cap_NonMarketRemainGDP_ann[yr,r] / p.GDP_per_cap_focus_0_FocusRegionEU)^p.ipow_incomeexponent

    if is_first(t)
        v.igdp_realizeddiscimpact_ann[yr,r] = v.occurdis_occurrencedummy_ann[yr] * (1 - v.expfdis_discdecay_ann[yr]) * v.igdpeqdis_eqdiscimpact_ann[yr,r]
    else
        v.igdp_realizeddiscimpact_ann[yr,r] = v.igdp_realizeddiscimpact_ann[yr - 1,r] + v.occurdis_occurrencedummy_ann[yr] * (1 - v.expfdis_discdecay_ann[yr]) * (v.igdpeqdis_eqdiscimpact_ann[yr,r] - v.igdp_realizeddiscimpact_ann[yr - 1,r])
    end

    if v.igdp_realizeddiscimpact_ann[yr,r] < p.isatg_saturationmodification    #### analyse this in dissertation?
        v.isat_satdiscimpact_ann[yr,r] = v.igdp_realizeddiscimpact_ann[yr,r]
    else
        v.isat_satdiscimpact_ann[yr,r] = p.isatg_saturationmodification + (100 - p.isatg_saturationmodification) * ((v.igdp_realizeddiscimpact_ann[yr,r] - p.isatg_saturationmodification) / ((100 - p.isatg_saturationmodification) + (v.igdp_realizeddiscimpact_ann[yr,r] - p.isatg_saturationmodification)))
    end

    v.isat_per_cap_DiscImpactperCapinclSaturation_ann[yr,r] = (v.isat_satdiscimpact_ann[yr,r] / 100) * p.rgdp_per_cap_NonMarketRemainGDP_ann[yr,r]
    v.rcons_per_cap_DiscRemainConsumption_ann[yr,r] = p.rcons_per_cap_NonMarketRemainConsumption_ann[yr,r] - v.isat_per_cap_DiscImpactperCapinclSaturation_ann[yr,r]
    if v.rcons_per_cap_DiscRemainConsumption_ann[yr,r] < 0
        v.rcons_per_cap_DiscRemainConsumption_ann[yr,r] = 0.
    end
end

@defcomp Discontinuity begin

    region = Index()
    year = Index()
    y_year = Parameter(index=[time], unit="year")
    y_year_0 = Parameter(unit="year")
    y_year_ann = Parameter(index=[year], unit="year")
    yagg_periodspan = Parameter(index=[time], unit="year") # for doing in-component aggregation

    rand_discontinuity = Parameter(unit="unitless", default=.5)

    irefeqdis_eqdiscimpact = Variable(index=[region], unit="%")
    wincf_weightsfactor_sea = Parameter(index=[region], unit="")
    wdis_gdplostdisc = Parameter(unit="%", default=3.)

    igdpeqdis_eqdiscimpact = Variable(index=[time,region], unit="%")
    igdpeqdis_eqdiscimpact_ann = Variable(index=[year,region], unit="%")
    rgdp_per_cap_NonMarketRemainGDP = Parameter(index=[time,region], unit="\$/person")
    rgdp_per_cap_NonMarketRemainGDP_ann = Parameter(index=[year,region], unit="\$/person")
    GDP_per_cap_focus_0_FocusRegionEU = Parameter(unit="\$/person", default=34298.93698672955)
    ipow_incomeexponent = Parameter(unit="unitless", default=-0.13333333333333333)

    igdp_realizeddiscimpact = Variable(index=[time,region], unit="%")
    igdp_realizeddiscimpact_ann = Variable(index=[year,region], unit="%")
    occurdis_occurrencedummy = Variable(index=[time], unit="unitless")
    occurdis_occurrencedummy_sum = Variable(unit="unitless") # for analysis
    occurdis_occurrencedummy_ann = Variable(index=[year], unit="unitless")
    occurdis_occurrencedummy_ann_sum = Variable(unit="unitless") # for analysis
    expfdis_discdecay = Variable(index=[time], unit="unitless")
    expfdis_discdecay_ann = Variable(index=[year], unit="unitless")

    distau_discontinuityexponent = Parameter(unit="unitless", default=20.)

    idis_lossfromdisc = Variable(index=[time], unit="degreeC")
    idis_lossfromdisc_ann = Variable(index=[year], unit="degreeC")
    tdis_tolerabilitydisc = Parameter(unit="degreeC", default=1.5)
    rt_g_globaltemperature = Parameter(index=[time], unit="degreeC")
    rt_g_globaltemperature_ann = Parameter(index=[year], unit="degreeC")
    pdis_probability = Parameter(unit="%/degreeC", default=20.)

    isatg_saturationmodification = Parameter(unit="unitless")
    isat_satdiscimpact = Variable(index=[time,region], unit="%")
    isat_satdiscimpact_ann = Variable(index=[year,region], unit="%")

    isat_per_cap_DiscImpactperCapinclSaturation = Variable(index=[time,region], unit="\$/person")
    isat_per_cap_DiscImpactperCapinclSaturation_sum = Variable(unit="\$/person") # for analysis
    isat_per_cap_DiscImpactperCapinclSaturation_ann = Variable(index=[year,region], unit="\$/person")
    isat_per_cap_DiscImpactperCapinclSaturation_ann_sum = Variable(unit="\$/person") # for analysis
    rcons_per_cap_DiscRemainConsumption = Variable(index=[time, region], unit="\$/person")
    rcons_per_cap_DiscRemainConsumption_sum =  Variable(unit="\$/person") # for analysis
    rcons_per_cap_DiscRemainConsumption_ann = Variable(index=[year, region], unit="\$/person")
    rcons_per_cap_DiscRemainConsumption_ann_sum = Variable(unit="\$/person") # for analysis
    rcons_per_cap_NonMarketRemainConsumption = Parameter(index=[time, region], unit="\$/person")
    rcons_per_cap_NonMarketRemainConsumption_ann = Parameter(index=[year, region], unit="\$/person")

    function run_timestep(p, v, d, t)

        v.idis_lossfromdisc[t] = max(0, p.rt_g_globaltemperature[t] - p.tdis_tolerabilitydisc)

        if is_first(t)
            if v.idis_lossfromdisc[t] * (p.pdis_probability / 100) > p.rand_discontinuity
                v.occurdis_occurrencedummy[t] = 1
            else
                v.occurdis_occurrencedummy[t] = 0
            end
            v.expfdis_discdecay[t] = exp(-(p.y_year[t] - p.y_year_0) / p.distau_discontinuityexponent)
        else
            if v.idis_lossfromdisc[t] * (p.pdis_probability / 100) > p.rand_discontinuity
                v.occurdis_occurrencedummy[t] = 1
            elseif v.occurdis_occurrencedummy[t - 1] == 1
                v.occurdis_occurrencedummy[t] = 1
            else
                v.occurdis_occurrencedummy[t] = 0
            end
            v.expfdis_discdecay[t] = exp(-(p.y_year[t] - p.y_year[t - 1]) / p.distau_discontinuityexponent)
        end

        for r in d.region
            v.irefeqdis_eqdiscimpact[r] = p.wincf_weightsfactor_sea[r] * p.wdis_gdplostdisc

            v.igdpeqdis_eqdiscimpact[t,r] = v.irefeqdis_eqdiscimpact[r] * (p.rgdp_per_cap_NonMarketRemainGDP[t,r] / p.GDP_per_cap_focus_0_FocusRegionEU)^p.ipow_incomeexponent

            if is_first(t)
                v.igdp_realizeddiscimpact[t,r] = v.occurdis_occurrencedummy[t] * (1 - v.expfdis_discdecay[t]) * v.igdpeqdis_eqdiscimpact[t,r]
            else
                v.igdp_realizeddiscimpact[t,r] = v.igdp_realizeddiscimpact[t - 1,r] + v.occurdis_occurrencedummy[t] * (1 - v.expfdis_discdecay[t]) * (v.igdpeqdis_eqdiscimpact[t,r] - v.igdp_realizeddiscimpact[t - 1,r])
            end

            if v.igdp_realizeddiscimpact[t,r] < p.isatg_saturationmodification
                v.isat_satdiscimpact[t,r] = v.igdp_realizeddiscimpact[t,r]
            else
                v.isat_satdiscimpact[t,r] = p.isatg_saturationmodification + (100 - p.isatg_saturationmodification) * ((v.igdp_realizeddiscimpact[t,r] - p.isatg_saturationmodification) / ((100 - p.isatg_saturationmodification) + (v.igdp_realizeddiscimpact[t,r] - p.isatg_saturationmodification)))
            end
            v.isat_per_cap_DiscImpactperCapinclSaturation[t,r] = (v.isat_satdiscimpact[t,r] / 100) * p.rgdp_per_cap_NonMarketRemainGDP[t,r]
            v.rcons_per_cap_DiscRemainConsumption[t,r] = p.rcons_per_cap_NonMarketRemainConsumption[t,r] - v.isat_per_cap_DiscImpactperCapinclSaturation[t,r]

            # calculate  for this specific year
            if is_first(t)
                # v.rcons_per_cap_DiscRemainConsumption_ann_sum = 0
                for annual_year = 2015:(gettime(t))
                    calc_discontinuity(p, v, d, t, annual_year, r)
                end
                if isequal(r, 8)
                    v.occurdis_occurrencedummy_sum = v.occurdis_occurrencedummy[t] * p.yagg_periodspan[t]
                    v.rcons_per_cap_DiscRemainConsumption_sum = sum(v.rcons_per_cap_DiscRemainConsumption[t,:]) * p.yagg_periodspan[t]
                    v.isat_per_cap_DiscImpactperCapinclSaturation_sum = sum(v.isat_per_cap_DiscImpactperCapinclSaturation[t,:]) * p.yagg_periodspan[t]
                end
            else
                for annual_year = (gettime(t - 1) + 1):(gettime(t))
                    calc_discontinuity(p, v, d, t, annual_year, r)

                    if isequal(annual_year, 2300)
                        if isequal(r, 8)
                            v.occurdis_occurrencedummy_ann_sum = sum(v.occurdis_occurrencedummy_ann[:])
                            v.rcons_per_cap_DiscRemainConsumption_ann_sum =  sum(v.rcons_per_cap_DiscRemainConsumption_ann[:,:])
                            v.isat_per_cap_DiscImpactperCapinclSaturation_ann_sum = sum(v.isat_per_cap_DiscImpactperCapinclSaturation_ann[:,:])
                        end
                    end
                end
                if isequal(r, 8)
                    # for analysis
                    v.occurdis_occurrencedummy_sum = v.occurdis_occurrencedummy_sum + v.occurdis_occurrencedummy[t] * p.yagg_periodspan[t]
                    v.rcons_per_cap_DiscRemainConsumption_sum = v.rcons_per_cap_DiscRemainConsumption_sum + sum(v.rcons_per_cap_DiscRemainConsumption[t,:]) * p.yagg_periodspan[t]
                    v.isat_per_cap_DiscImpactperCapinclSaturation_sum = v.isat_per_cap_DiscImpactperCapinclSaturation_sum + sum(v.isat_per_cap_DiscImpactperCapinclSaturation[t,:]) * p.yagg_periodspan[t]
                end

            end

        end
    end
end


text:
# Linear Discriminant Analysis

#### Type to represent a linear discriminant functional

abstract Discriminant

immutable LinearDiscriminant <: Discriminant
    w::Vector{Float64}
    b::Float64
end

length(f::LinearDiscriminant) = length(f.w)

evaluate(f::LinearDiscriminant, x::AbstractVector) = dot(f.w, x) + f.b

function evaluate(f::LinearDiscriminant, X::AbstractMatrix)
    R = At_mul_B(X, f.w)
    if f.b != 0
        broadcast!(+, R, R, f.b)
    end
    return R
end

predict(f::Discriminant, x::AbstractVector) = evaluate(f, x) > 0

predict(f::Discriminant, X::AbstractMatrix) = (Y = evaluate(f, X); Bool[y > 0 for y in Y])

#### function to solve linear discriminant

function ldacov(C::DenseMatrix{Float64},
                μp::DenseVector{Float64},
                μn::DenseVector{Float64})

    w = cholfact(C) \ (μp - μn)
    ap = dot(w, μp)
    an = dot(w, μn)
    c = 2 / (ap - an)
    LinearDiscriminant(scale!(w, c), 1 - c * ap)
end

ldacov(Cp::DenseMatrix{Float64},
       Cn::DenseMatrix{Float64},
       μp::DenseVector{Float64},
       μn::DenseVector{Float64}) = ldacov(Cp + Cn, μp, μn)

#### interface functions

function fit(::Type{LinearDiscriminant}, Xp::DenseMatrix{Float64}, Xn::DenseMatrix{Float64})
    μp = vec(mean(Xp, 2))
    μn = vec(mean(Xn, 2))
    Zp = Xp .- μp
    Zn = Xn .- μn
    Cp = A_mul_Bt(Zp, Zp)
    Cn = A_mul_Bt(Zn, Zn)
    ldacov(Cp, Cn, μp, μn)
end

#### Multiclass LDA Stats

type MulticlassLDAStats
    dim::Int                    # sample dimensions
    nclasses::Int               # number of classes
    cweights::Vector{Float64}   # class weights
    tweight::Float64            # total sample weight
    mean::Vector{Float64}       # overall sample mean
    cmeans::Matrix{Float64}     # class-specific means
    Sw::Matrix{Float64}         # within-class scatter matrix
    Sb::Matrix{Float64}         # between-class scatter matrix
end

Base.mean(S::MulticlassLDAStats) = S.mean
classweights(S::MulticlassLDAStats) = S.cweights
classmeans(S::MulticlassLDAStats) = S.cmeans

withclass_scatter(S::MulticlassLDAStats) = S.Sw
betweenclass_scatter(S::MulticlassLDAStats) = S.Sb

function MulticlassLDAStats(cweights::Vector{Float64},
                            mean::Vector{Float64},
                            cmeans::Matrix{Float64},
                            Sw::Matrix{Float64},
                            Sb::Matrix{Float64})
    d, nc = size(cmeans)
    length(mean) == d || throw(DimensionMismatch("Incorrect length of mean"))
    length(cweights) == nc || throw(DimensionMismatch("Incorrect length of cweights"))
    tw = sum(cweights)
    size(Sw) == (d, d) || throw(DimensionMismatch("Incorrect size of Sw"))
    size(Sb) == (d, d) || throw(DimensionMismatch("Incorrect size of Sb"))
    MulticlassLDAStats(d, nc, cweights, tw, mean, cmeans, Sw, Sb)
end

function multiclass_lda_stats(nc::Int, X::DenseMatrix{Float64}, y::AbstractVector{Int})
    # check sizes
    d = size(X, 1)
    n = size(X, 2)
    n >= nc || error("The number of samples is less than the number of classes")
    length(y) == n || throw(DimensionMismatch("Inconsistent array sizes."))

    # compute class-specific weights and means
    cweights = zeros(nc)
    cmeans = zeros(d, nc)
    for j = 1:n
        @inbounds c = y[j]
        1 <= c <= nc || error("class label must be in [1, nc].")
        cweights[c] += 1
        v = view(cmeans,:,c)
        x = view(X,:,j)
        for i = 1:d
            @inbounds v[i] += x[i]
        end
    end
    for j = 1:nc
        @inbounds cw = cweights[j]
        cw > 0 || error("The $(j)-th class has no sample.")
        scale!(view(cmeans,:,j), inv(cw))
    end

    # compute within-class scattering
    Z = Array(Float64, d, n)
    for j = 1:n
        z = view(Z,:,j)
        x = view(X,:,j)
        v = view(cmeans,:,y[j])
        for i = 1:d
            z[i] = x[i] - v[i]
        end
    end
    Sw = A_mul_Bt(Z, Z)

    # compute between-class scattering
    mean = cmeans * (cweights ./ n)
    U = scale!(cmeans .- mean, sqrt(cweights))
    Sb = A_mul_Bt(U, U)

    return MulticlassLDAStats(cweights, mean, cmeans, Sw, Sb)
end

#### Multiclass LDA

type MulticlassLDA
    proj::Matrix{Float64}
    pmeans::Matrix{Float64}
    stats::MulticlassLDAStats
end

indim(M::MulticlassLDA) = size(M.proj, 1)
outdim(M::MulticlassLDA) = size(M.proj, 2)

projection(M::MulticlassLDA) = M.proj

Base.mean(M::MulticlassLDA) = mean(M.stats)
classmeans(M::MulticlassLDA) = classmeans(M.stats)
classweights(M::MulticlassLDA) = classweights(M.stats)

withclass_scatter(M::MulticlassLDA) = withclass_scatter(M.stats)
betweenclass_scatter(M::MulticlassLDA) = betweenclass_scatter(M.stats)

transform{T<:FloatingPoint}(M::MulticlassLDA, x::AbstractVecOrMat{T}) = M.proj'x

function fit(::Type{MulticlassLDA}, nc::Int, X::DenseMatrix{Float64}, y::AbstractVector{Int};
             method::Symbol=:gevd,
             outdim::Int=min(size(X,1), nc-1),
             regcoef::Float64=1.0e-6)

    multiclass_lda(multiclass_lda_stats(nc, X, y);
                   method=method,
                   regcoef=regcoef,
                   outdim=outdim)
end

function multiclass_lda(S::MulticlassLDAStats;
                        method::Symbol=:gevd,
                        outdim::Int=min(size(X,1), S.nclasses-1),
                        regcoef::Float64=1.0e-6)

    P = mclda_solve(S.Sb, S.Sw, method, outdim, regcoef)
    MulticlassLDA(P, P'S.cmeans, S)
end

mclda_solve(Sb::DenseMatrix{Float64}, Sw::DenseMatrix{Float64}, method::Symbol, p::Int, regcoef::Float64) =
    mclda_solve!(copy(Sb), copy(Sw), method, p, regcoef)

function mclda_solve!(Sb::Matrix{Float64},
                      Sw::Matrix{Float64},
                      method::Symbol, p::Int, regcoef::Float64)

    p <= size(Sb, 1) || error("p cannot exceed sample dimension.")

    if method == :gevd
        regularize_symmat!(Sw, regcoef)
        E = eigfact!(Symmetric(Sb), Symmetric(Sw))
        ord = sortperm(E.values; rev=true)
        P = E.vectors[:, ord[1:p]]

    elseif method == :whiten
        W = _lda_whitening!(Sw, regcoef)
        wSb = At_mul_B(W, Sb * W)
        Eb = eigfact!(Symmetric(wSb))
        ord = sortperm(Eb.values; rev=true)
        P = W * Eb.vectors[:, ord[1:p]]

    else
        error("Invalid method name $(method)")
    end
    return P::Matrix{Float64}
end

function _lda_whitening!(C::Matrix{Float64}, regcoef::Float64)
    n = size(C,1)
    E = eigfact!(Symmetric(C))
    v = E.values
    a = regcoef * maximum(v)
    for i = 1:n
        @inbounds v[i] = 1.0 / sqrt(v[i] + a)
    end
    return scale!(E.vectors, v)
end


text:
# SPDX-License-Identifier: X11
# 2021-01-03
# Otoshidama (300pt)

function main()
  N, Y = parse.(Int, (split ∘ readline)())

  for i ∈ 0:N, j ∈ 0:N - i
    k = N - i - j
    k < 0 && continue
    i * 10000 + j * 5000 + k * 1000 == Y && (println("$i $j $k"); return)
  end

  println("-1 -1 -1")
end

main()

max_stars_repo_name: jthvai/competitive-programming

max_issues_repo_name: jthvai/competitive-programming

max_forks_repo_name: jthvai/competitive-programming


text:
using ArnoldiMethod
using ArnoldiMethodTransformations
using LinearAlgebra
using SparseArrays
using Test

# construct fixed eval matrix in random basis
D = diagm(0=>[0,1,2,3,4,5,6,7,8,9])

@testset "functionlity: " begin
    S = I + .1*randn(10,10) # random basis
    A = S\D*S
    S = I + .1*randn(10,10) # random basis
    B = S\D*S

    si = ArnoldiMethodTransformations.ShiftAndInvert(A;σ=5.001)
    decomp, hist = partialschur(si)
    λ, v = partialeigen(decomp,5.001)
    @test norm(A*v-v*diagm(0=>λ))/prod(size(v)) ≤ 5*sqrt(eps())
    decomp, hist = partialschur(A,5.001)
    λ, v = partialeigen(decomp,5.001)
    @test norm(A*v-v*diagm(0=>λ))/prod(size(v)) ≤ 5*sqrt(eps())

    As = sparse(A)
    si = ArnoldiMethodTransformations.ShiftAndInvert(As;σ=5.001)
    decomp, hist = partialschur(si)
    λ, v = partialeigen(decomp,5.001)
    @test norm(A*v-v*diagm(0=>λ))/prod(size(v)) ≤ 5*sqrt(eps())
    decomp, hist = partialschur(A,5.001)
    λ, v = partialeigen(decomp,5.001)
    @test norm(A*v-v*diagm(0=>λ))/prod(size(v)) ≤ 5*sqrt(eps())

    si = ArnoldiMethodTransformations.ShiftAndInvert(A,B;σ=5.001)
    decomp, hist = partialschur(si)
    λ, v = partialeigen(decomp,5.001)
    @test norm(A*v-B*v*diagm(0=>λ))/prod(size(v)) ≤ 5*sqrt(eps())
    decomp, hist = partialschur(A,B,5.001)
    λ, v = partialeigen(decomp,5.001)
    @test norm(A*v-B*v*diagm(0=>λ))/prod(size(v)) ≤ 5*sqrt(eps())

    si = ArnoldiMethodTransformations.ShiftAndInvert(A,B;σ=5.001)
    Bs = sparse(B)
    decomp, hist = partialschur(si)
    λ, v = partialeigen(decomp,5.001)
    @test norm(As*v-Bs*v*diagm(0=>λ))/prod(size(v)) ≤ 5*sqrt(eps())
    decomp, hist = partialschur(As,Bs,5.001)
    λ, v = partialeigen(decomp,5.001)
    @test norm(As*v-Bs*v*diagm(0=>λ))/prod(size(v)) ≤ 5*sqrt(eps())
end

# @code_warntype ArnoldiMethodTransformations.ShiftAndInvert(A;σ=5.001)
# @time si = ArnoldiMethodTransformations.ShiftAndInvert(A;σ=5.001)
#
# @code_warntype ArnoldiMethodTransformations.ShiftAndInvert(As;σ=5.001)
# @time si = ArnoldiMethodTransformations.ShiftAndInvert(As;σ=5.001)
#
# @code_warntype ArnoldiMethodTransformations.ShiftAndInvert(A,B;σ=5.001)
# @time si = ArnoldiMethodTransformations.ShiftAndInvert(A,B;σ=5.001)
#
# @code_warntype ArnoldiMethodTransformations.ShiftAndInvert(As,Bs;σ=5.001)
# @time si = ArnoldiMethodTransformations.ShiftAndInvert(As,Bs;σ=5.001)


text:
module FluxReconstruction

const FR = FluxReconstruction

using Base.Threads: @threads
using CUDA
using GSL
using LinearAlgebra
using NonlinearSolve
using OrdinaryDiffEq
using PyCall
using KitBase
using KitBase.FastGaussQuadrature
using KitBase.OffsetArrays
using KitBase.SpecialFunctions

export FR
export AbstractElementShape,
       Line,
       Quad,
       Tri,
       Hex,
       Wed,
       Pyr,
       Tet
export shock_detector,
       positive_limiter
export legendre_point,
       lagrange_point,
       ∂legendre,
       ∂radau,
       ∂lagrange,
       standard_lagrange,
       simplex_basis,
       ∂simplex_basis,
       vandermonde_matrix,
       ∂vandermonde_matrix,
       correction_field,
       JacobiP,
       ∂JacobiP,
       modal_filter!
export tri_quadrature, triface_quadrature
export FRPSpace1D,
       FRPSpace2D,
       UnstructFRPSpace,
       TriFRPSpace,
       global_sp,
       global_fp,
       rs_ab,
       xy_rs,
       rs_xy,
       rs_jacobi,
       neighbor_fpidx
export interp_face!
export poly_derivative!
export FREulerProblem

include("data.jl")
include("dissipation.jl")
include("Polynomial/polynomial.jl")
include("Transform/transform.jl")
include("quadrature.jl")
include("struct.jl")
include("Geometry/geometry.jl")
include("tools.jl")
include("interpolate.jl")
include("derivative.jl")
include("Equation/equation.jl")
include("integrator.jl")

end


text:
function analyze(code::String; input::IO = stdin, output::IO = stdout)
    raw_program = parse(code; input, output)
    if !isnothing(raw_program)
        macro_env = Environment(; input, output)
        program = define_macros!(macro_env, raw_program)
        expanded = expand_macros(program, macro_env)
        return analyze(expanded)
    end
end

function analyze(
    program::Program;
    existing_symbol_table::Union{SymbolTable,Nothing} = nothing,
    exisiting_env::Union{Environment,Nothing} = nothing,
)
    if !isnothing(existing_symbol_table)
        symbol_table = SymbolTable(existing_symbol_table)
    else
        symbol_table = SymbolTable()
        for (i, (name, _)) in enumerate(BUILTINS)
            define_builtin!(symbol_table, name, i - 1)
        end

        if !isnothing(exisiting_env)
            for key in keys(exisiting_env.store)
                define!(symbol_table, key)
            end
        end
    end

    for statement in program.statements
        result = analyze(statement, symbol_table)
        if isa(result, ErrorObj)
            return result
        end
    end
end

function analyze(bs::BlockStatement, symbol_table::SymbolTable)
    for statement in bs.statements
        result = analyze(statement, symbol_table)
        if isa(result, ErrorObj)
            return result
        end
    end
end

analyze(es::ExpressionStatement, symbol_table::SymbolTable) =
    analyze(es.expression, symbol_table)

function analyze(ls::LetStatement, symbol_table::SymbolTable)
    sym, _ = resolve(symbol_table, ls.name.value)

    if ls.reassign

        if isnothing(sym)
            return ErrorObj("identifier not found: $(ls.name.value)")
        end

        if sym.scope == FunctionScope ||
           (sym.scope == OuterScope && sym.ptr.scope == FunctionScope)
            return ErrorObj(
                "cannot reassign the current function being defined: $(ls.name.value)",
            )
        end
    else
        if !isnothing(sym) && (
            sym.scope == LocalScope ||
            (is_global(symbol_table) && sym.scope == GlobalScope)
        )
            return ErrorObj("$(ls.name.value) is already defined")
        end

        define!(symbol_table, ls.name.value)
    end

    return analyze(ls.value, symbol_table)
end

function analyze(ws::WhileStatement, symbol_table::SymbolTable)
    result = analyze(ws.condition, symbol_table)
    if isa(result, ErrorObj)
        return result
    end

    inner = SymbolTable(; outer = symbol_table, within_loop = true)
    return analyze(ws.body, inner)
end

analyze(rs::ReturnStatement, symbol_table::SymbolTable) =
    analyze(rs.return_value, symbol_table)

function analyze(::BreakStatement, symbol_table::SymbolTable)
    if !within_loop(symbol_table)
        return ErrorObj("syntax error: break outside loop")
    end
end

function analyze(::ContinueStatement, symbol_table::SymbolTable)
    if !within_loop(symbol_table)
        return ErrorObj("syntax error: continue outside loop")
    end
end

function analyze(ident::Identifier, symbol_table::SymbolTable)
    sym, _ = resolve(symbol_table, ident.value)
    if isnothing(sym)
        return ErrorObj("identifier not found: $(ident.value)")
    end
end

function analyze(al::ArrayLiteral, symbol_table::SymbolTable)
    for element in al.elements
        result = analyze(element, symbol_table)
        if isa(result, ErrorObj)
            return result
        end
    end
end

function analyze(hl::HashLiteral, symbol_table::SymbolTable)
    for (key, value) in hl.pairs
        result = analyze(key, symbol_table)
        if isa(result, ErrorObj)
            return result
        end

        result = analyze(value, symbol_table)
        if isa(result, ErrorObj)
            return result
        end
    end
end

function analyze(fl::FunctionLiteral, symbol_table::SymbolTable)
    inner = SymbolTable(; outer = symbol_table)

    if !isempty(fl.name)
        define_function!(inner, fl.name)
    end

    for param in fl.parameters
        define!(inner, param.value)
    end

    return analyze(fl.body, inner)
end

function analyze(ie::IfExpression, symbol_table::SymbolTable)
    result = analyze(ie.condition, symbol_table)
    if isa(result, ErrorObj)
        return result
    end

    result = analyze(ie.consequence, symbol_table)
    if isa(result, ErrorObj)
        return result
    end

    if !isnothing(ie.alternative)
        return analyze(ie.alternative, symbol_table)
    end
end

analyze(pe::PrefixExpression, symbol_table::SymbolTable) = analyze(pe.right, symbol_table)

function analyze(ie::InfixExpression, symbol_table::SymbolTable)
    result = analyze(ie.left, symbol_table)
    if isa(result, ErrorObj)
        return result
    end

    return analyze(ie.right, symbol_table)
end

function analyze(ie::IndexExpression, symbol_table::SymbolTable)
    result = analyze(ie.left, symbol_table)
    if isa(result, ErrorObj)
        return result
    end

    return analyze(ie.index, symbol_table)
end

function analyze(ce::CallExpression, symbol_table::SymbolTable)
    if (token_literal(ce.fn) == "quote")
        return nothing
    end

    result = analyze(ce.fn, symbol_table)
    if isa(result, ErrorObj)
        return result
    end

    for arg in ce.arguments
        result = analyze(arg, symbol_table)
        if isa(result, ErrorObj)
            return result
        end
    end
end

# The following functions do nothing in the current version.

analyze(il::IntegerLiteral, symbol_table::SymbolTable) = nothing

analyze(sl::StringLiteral, symbol_table::SymbolTable) = nothing

analyze(bl::BooleanLiteral, symbol_table::SymbolTable) = nothing

analyze(::NullLiteral, ::SymbolTable) = nothing


text:
"""
    SharedMATLABEngine

## Exports
- connect_matlab
- start_matlab
- find_matlab
- mat"

## Description
Run MATLAB code in an open MATLAB engine session. To get started, call the function
`connect_matlab(engine_name)` in the Julia REPL. To see a list of open engine names, call
`find_matlab()`. Once connected, the MATLAB command line can be accessedfrom the Julia REPL
by typing `>`. Julia variables can be interpolated into MATLAB commands via the `\$` operator.

## Examples
```julia
julia> using SharedMATLABEngine

julia> eng = connect_matlab(:default_engine)
REPL mode MATLAB initialized. Press > to enter and backspace to exit.
SharedMATLABEngine.Engine()

>> a = magic(3)

a =

     8     1     6
     3     5     7
     4     9     2

>> \$a21 = a(2,1)
3.0

julia> a21 .+ mat"a"
3×3 Array{Float64,2}:
 11.0   4.0   9.0
  6.0   8.0  10.0
  7.0  12.0   5.0
```
"""
module SharedMATLABEngine

using PyCall
using ReplMaker

include("utils.jl")
include("engine.jl")
include("matrepl_str.jl")

export connect_matlab, start_matlab, find_matlab
export @mat_str

# Note: Don't fix the indentation here, it breaks Julia syntax highlighting
function __init__()
py"""
import matlab.engine
import numpy as np
"""
end

end


text:
struct Cursor
    pos
end

struct Action
    getter
    permitted
end

Base.getindex(a::Action, n) = a.getter(n)

struct Mode
    signature::String
    actions::Action
end

function Mode(signature::String, d::Dict)
    function f(x)
        if in(x,keys(d))
            return d[x]
        else
            return y->push!(y.args, x)
        end
    end
    Mode(signature, Action(f, x->in(x,keys(d))))
end

struct Buffer
    text::Array{String,1}
    undo::Array{String,1} #Given time, this should be replaced with multi-level undo with diffs
    redo::Array{String,1}
    cursor::Cursor
    mode::Array{Mode,1}
    state::Dict{Symbol, String}
    registers::Dict{Char, String}
    args::Array{Char,1}
end

function settext(b::Buffer, text::Array{String,1})
    resize!(b.text, length(text))
    b.text .= text
end

function setundo(b::Buffer, text)
    b.state[:undo] = "undo"
    resize!(b.undo, length(text))
    b.undo .= text
end
setundo(b::Buffer) = setundo(b, b.text)

function setredo(b::Buffer, text)
    b.state[:undo] = "redo"
    resize!(b.redo, length(text))
    b.redo .= text
end
setredo(b::Buffer) = setredo(b, b.text)

mode(b::Buffer) = b.mode[1]
function setmode(b::Buffer, m::Mode)
    if m!=normal_mode
        setundo(b)
    end
    b.mode[1] = m
end

function Base.show(io::IO, b::Buffer)
    print(io, "Buffer: ")
    print(io, "$(length(b.text)) lines ")
    print(io, "$(mode(b).signature) ")
    print(io, "cursor position $(b.cursor.pos[1]),$(b.cursor.pos[2])")
end

function clear_screen()
    write(STDOUT, "\e[2J");
end

function move_sys_cursor(row,col)
    write(STDOUT, "\e[$(row);$(col)H");
end

function render_line(b::Buffer, l::String)
    s = Base.clamp(left(b), 1, max(1,length(l)))
    f = Base.clamp(right(b), 1, length(l))
    if length(l)>=1
        displ = l[chr2ind(l, s):chr2ind(l, f)]
    else
        displ = ""
    end
    if :syntax in keys(b.state) && b.state[:syntax] == "true"
        displ = highlight(displ)
    end
    write(STDOUT, "$displ\n")
end

function render(b::Buffer)
    resize(b)
    clear_screen()
    move_sys_cursor(1,1)
    #Change parsing to something else
    for l in b.text[max(1,top(b)):min(end, bottom(b))]
        render_line(b, l)
    end
    #Change hardcoded stuff
    tildes = max(0, (bottom(b)-top(b)+1)-length(b.text))
    write(STDOUT, string("~\n"^tildes, b.state[:console]))
end


text:
@recipe function fmultiple(o1::AbstractUncertainValue,
                            o2::AbstractUncertainValue;
                            mix = false,
                            n_samples = 10000, nbins = 50)
    size --> (300, 400)
    d1 = o1.distribution
    d2 = o1.distribution
    @series begin
        label --> "P1, $d1"
        seriestype := :bar
        fα --> 0.4
        fc --> :green
        fit(Histogram, resample(d1, n_samples), nbins = nbins)
    end
    @series begin
        label --> "P2, $d2"
        seriestype := :bar
        fc --> :blue
        fα --> 0.4
        fit(Histogram, resample(d2, n_samples), nbins = nbins)
    end

    if mix
        M = MixtureModel([d1, d2])

        @series begin
            label --> "MixtureModel with uniform priors"
            seriestype := :bar
            fα --> 0.6
            fc --> :black
            fit(Histogram, rand(M, n_samples), nbins = nbins)
        end
    end
end


text:
# Note: Test environment is setup by https://github.com/JuliaCI/PkgEvalHadoopEnv
include("hive_test.jl")


text:
####################
# value extraction #
####################

@inline extract_value!(::Type{T}, out::DiffResult, ydual) where {T} =
    DiffResults.value!(d -> value(T,d), out, ydual)
@inline extract_value!(::Type{T}, out, ydual) where {T} = out # ???

@inline function extract_value!(::Type{T}, out, y, ydual) where {T}
    map!(d -> value(T,d), y, ydual)
    copy_value!(out, y)
end

@inline copy_value!(out::DiffResult, y) = DiffResults.value!(out, y)
@inline copy_value!(out, y) = out

###################################
# vector mode function evaluation #
###################################

@generated function dualize(::Type{T}, x::SArray{S,V,D,N}) where {T,S,V,D,N}
    dx = Expr(:tuple, [:(Dual{T}(x[$i], chunk, Val{$i}())) for i in 1:N]...)
    return quote
        chunk = Chunk{N}()
        $(Expr(:meta, :inline))
        return SArray{S}($(dx))
    end
end

@inline static_dual_eval(::Type{T}, f, x::SArray) where {T} = f(dualize(T, x))

function vector_mode_dual_eval(f, x, cfg::Union{JacobianConfig,GradientConfig})
    xdual = cfg.duals
    seed!(xdual, x, cfg.seeds)
    return f(xdual)
end

function vector_mode_dual_eval(f!, y, x, cfg::JacobianConfig)
    ydual, xdual = cfg.duals
    seed!(xdual, x, cfg.seeds)
    seed!(ydual, y)
    f!(ydual, xdual)
    return ydual
end

##################################
# seed construction/manipulation #
##################################

@generated function construct_seeds(::Type{Partials{N,V}}) where {N,V}
    return Expr(:tuple, [:(single_seed(Partials{N,V}, Val{$i}())) for i in 1:N]...)
end

function seed!(duals::AbstractArray{Dual{T,V,N}}, x,
               seed::Partials{N,V} = zero(Partials{N,V})) where {T,V,N}
    for i in eachindex(duals)
        duals[i] = Dual{T,V,N}(x[i], seed)
    end
    return duals
end

function seed!(duals::AbstractArray{Dual{T,V,N}}, x,
               seeds::NTuple{N,Partials{N,V}}) where {T,V,N}
    for i in 1:N
        duals[i] = Dual{T,V,N}(x[i], seeds[i])
    end
    return duals
end

function seed!(duals::AbstractArray{Dual{T,V,N}}, x, index,
               seed::Partials{N,V} = zero(Partials{N,V})) where {T,V,N}
    offset = index - 1
    for i in 1:N
        j = i + offset
        duals[j] = Dual{T,V,N}(x[j], seed)
    end
    return duals
end

function seed!(duals::AbstractArray{Dual{T,V,N}}, x, index,
               seeds::NTuple{N,Partials{N,V}}, chunksize = N) where {T,V,N}
    offset = index - 1
    for i in 1:chunksize
        j = i + offset
        duals[j] = Dual{T,V,N}(x[j], seeds[i])
    end
    return duals
end


text:
# This file was generated, do not modify it. # hide
rpca = collect(values(report(spca).report_given_machine))[2]
cs = cumsum(rpca.principalvars ./ rpca.tvar)


text:
using LinearAlgebra, NLPModels, Stopping

include("backls.jl")
include("uncons.jl")

#https://juliasmoothoptimizers.github.io/SolverBenchmark.jl/latest/tutorial/
#In this tutorial we illustrate the main uses of SolverBenchmark.
using DataFrames, Printf, SolverBenchmark
#CUTEst is a collection of test problems
using CUTEst

problems_unconstrained = CUTEst.select(contype="unc")
n = length(problems_unconstrained) #240
#problems_boundconstrained = CUTEst.select(contype="bounds")
#n = length(problems_boundconstrained) #124

printstyled("Benchmark solvers: \n", color = :green)

n = min(n, 3)

#Names of 3 solvers:
names = [:armijo, :wolfe, :armijo_wolfe]
p1 = PrmUn(); p2 = PrmUn(ls_func = wolfe); p3 = PrmUn(ls_func = armijo_wolfe)
paramDict = Dict(:armijo => p1, :wolfe => p2, :armijo_wolfe => p3)
#Initialization of the DataFrame for n problems.
stats = Dict(name => DataFrame(:id => 1:n,
         :name => [@sprintf("prob%s", problems_unconstrained[i]) for i = 1:n],
         :nvar => zeros(Int64, n),
         :status => [:Unknown for i = 1:n],
         :f => NaN*ones(n),
         :t => NaN*ones(n),
         :iter => zeros(Int64, n),
         :eval_f => zeros(Int64, n),
         :eval_g => zeros(Int64, n),
         :eval_H => zeros(Int64, n),
         :score => NaN*ones(n)) for name in names)

for i=1:n
  nlp_cutest = CUTEst.CUTEstModel(problems_unconstrained[i])
  @show i, problems_unconstrained[i], nlp_cutest.meta.nvar
  #update the stopping with the new problem
  stop_nlp = NLPStopping(nlp_cutest,
                         NLPAtX(nlp_cutest.meta.x0),
                         max_iter = 20,
                         optimality_check = unconstrained_check)

  for name in names

    #solve the problem
    global_newton(stop_nlp, paramDict[name])

    #update the stats from the Stopping
    stats[name].nvar[i] = nlp_cutest.meta.nvar
    stats[name].status[i] = status(stop_nlp)
    stats[name].f[i] = stop_nlp.current_state.fx
    stats[name].t[i] = stop_nlp.current_state.current_time - stop_nlp.meta.start_time
    stats[name].iter[i] = stop_nlp.meta.nb_of_stop
    stats[name].score[i] = unconstrained_check(nlp_cutest, stop_nlp.current_state)
    stats[name].eval_f[i] = getfield(stop_nlp.current_state.evals, :neval_obj)
    stats[name].eval_g[i] = getfield(stop_nlp.current_state.evals, :neval_grad)
    stats[name].eval_H[i] = getfield(stop_nlp.current_state.evals, :neval_hess)

    #reinitialize the Stopping and the nlp
    reinit!(stop_nlp, rstate = true, x = nlp_cutest.meta.x0)
    reset!(stop_nlp.pb)
  end

  #finalize nlp
  finalize(nlp_cutest)
end #end of main loop

for name in names
@show stats[name]
end

#You can export the table in Latex
#latex_table(stdout, stats[:armijo])

#or run a performance profile:
#using Plots
#pyplot()
#cost(df) = (df.status != :Optimal) * Inf + df.t
#p = performance_profile(stats, cost)
#Plots.svg(p, "profile2")

#or a profile wall:
#solved(df) = (def.status .== :Optimal)
#costs = [df -> .!sovled(df) * Inf + df.t, df -> .!sovled(df) * Inf + df.iter]
#costnames = ["Time", "Iterations"]
#p = profile_solvers(stats, costs, costnames)
#Plots.svg(p, "profile3")

printstyled("The End.", color = :green)


text:
using Distributions
using Statistics
import LinearModel

function assert_same_rows(X, Y)
    if size(X,1) == size(Y,1)
        return true
    else
        throw(ArgumentError("X and Y must have same number of rows"))
    end
end

function assert_vector(Y)
    if size(Y,2) == 1
        return true
    else
        # more readable error message than struct type-checking error message
        throw(ArgumentError("Y must be a vector"))
    end
end

struct Sample
    X::Array
    Y::Vector
    Sample(X, Y) = (assert_same_rows(X, Y) & assert_vector(Y)
        ? new(X, Y)
        : throw(ArgumentError)
    )
end

# Create function that will return normal observations
function make_sample_generator_normal(; a, b, β_0, β, sd_e)
    k = size(β, 1)
    make_x_sample(n) = rand(Uniform(a, b), n, k)
    y_process(X) = β_0 .+ X * β
    noise_process(X) = rand(Normal(0, sd_e), size(X, 1), 1)
    function normal_sample(n)::Sample
        X = make_x_sample(n)
        Y = vec(y_process(X) + noise_process(X)) # vec reshapes Array{T, 2} to Vector
        return Sample(X, Y)
    end
    return normal_sample
end

function add_intercept(X::Array)
    nrows = size(X, 1)
    X0 = ones(nrows, 1)
    return [X0 X]
end

function add_intercept(sample::Sample)
    return Sample(add_intercept(sample.X), sample.Y)
end

normal_sampler = make_sample_generator_normal(a=0, b=20, β_0=5, β=[1,2], sd_e=1.1)
sam1 = normal_sampler(100)
lm1 = evaluate_ols(sam1, intercept=true)

println(show(lm1))
# TODO: add below this to show
println("R-squared: ", round(r2(lm1), digits=4))

# TODO: import data
# TODO: replicate https://itl.nist.gov/div898/strd/lls/data/LINKS/v-Norris.shtml

# the difference between two ols methods is really small
#@assert all(ols(obs1.X, obs1.Y)-ols_(obs1.X,obs1.Y) .< 1e-10) # was 14
#@assert add_intercept(collect(2:4)) == [1 1 1; 2 3 4]'
#@assert r2(lm1) > .8
using Test
@test_throws ArgumentError Sample([1, 2, 3], [1, 2])


text:
"""
	$(SIGNATURES)

This is used to list settings for nested `ModelObject`s.
"""
settings_table(x) = Matrix{String}();

"""
	$(SIGNATURES)

This is used to write literals from nested objects into preambles. Objects return `Vector{Vector{String}}` or `Vector{Tuple}`. Each entry contains (description, newcommand name, value).
"""
settings_list(x, st) = Vector{Vector{String}}();

"""
$(SIGNATURES)

Display current time and date as "2020-Jan-10 at 10:35am".
"""
function current_time()
    x = Dates.now();
    return Dates.format(x, "Y-u-d") * " at " * Dates.format(x, "HH:MM")
end

# Given a vector of things that can be converted to string, return a string of the form "abc_def_dhi" with connector character `linkChar`
function chain_strings(v :: AbstractVector{T}, linkStr = "_") where T
    s = string(v[1]);
    if length(v) > 1
        for j = 2 : length(v)
            s = s * string(linkStr) * string(v[j])
        end
    end
    return s
end

chain_strings(v :: AbstractString, linkStr = "_") = v;

"""
$(SIGNATURES)

Format a vector of numbers with a header. Output is a string.
"""
function format_vector(headerStr, numberV, nDigits :: Integer)
    roundV = round.(numberV, digits = nDigits);
    if isempty(headerStr)
        return "$roundV"
    else
        return headerStr * ":  $roundV"
    end
end

"""
$(SIGNATURES)

Show a text table.
"""
function show_text_table(tbM; io = stdout)
    pretty_table(io, tbM, noheader = true);
    return nothing
end

function show_text_table(tbM, headerV, fPath :: T1) where T1 <: AbstractString
    open(fPath, "w") do io
        pretty_table(io,  tbM,  headerV);
        println("Saved  $(fpath_to_show(fPath))")
    end
end

function show_matrix(io :: IO, m :: Matrix{F1},
ndigits :: Integer) where F1 <: AbstractFloat

    nr, nc = size(m);
    for ir = 1 : nr
        for ic = 1 : nc
            print(io, round(m[ir, ic], digits = ndigits), "\t");
        end
        println(io, " ");
    end
    println(io, " ");
end

## Format a single number in easy to read format.
function format_number(x :: Number; format :: Symbol = :default)
    if format == :dollars
        return format_dollars(x);
    elseif format == :default
        return format_number_default(x);
    elseif format == :none
        return string(x);
    else
        error("Invalid format: $format")
    end
end

## Format a number for display in a table
function format_number_default(x :: Number)
    absX = abs(x);
    if absX > 1e3
        y = format(x, precision = 1, commas = true);
    elseif absX > 1
        y = format(x, precision = 2, commas = false);
    elseif absX > 0.01
        y = format(x, precision = 4);
    else
        y = format(x, precision = 6);
    end
    return y
end

## Format a dollar number
function format_dollars(x; decimals = 0)
    format(x, precision = decimals, commas = true);
end

# Show a file name with `n` last directories
function fpath_to_show(fPath :: String; nDirs :: Integer = 3)
    fDir, fName = splitdir(fPath);
    s = right_dirs(fDir, nDirs) * " - " * fName;
    return s
end

# Display calibrated if true and fixed if false
function calibrated_string(isCal :: Bool; fixedValue = nothing)
    if isCal
        x = "calibrated"
    else
        if isnothing(fixedValue)
            valStr = "";
        else
            valStr = " with value $fixedValue";
        end
        x = "fixed" * valStr;
    end
    return x
end

# --------------


text:
using ApplicationBuilder

using Compat
using Compat.Pkg

# Allow this file to be called either as a standalone file to build the above
# example, or from runtests.jl using a provided builddir.
Compat.isdefined(:builddir) || (builddir="builddir")

libUIPkg = Pkg.dir("Libui")

@assert libUIPkg != nothing "LibUI is not installed!"

using Libui

ApplicationBuilder.build_app_bundle(joinpath(@__DIR__, "..", "..", "examples", "libui.jl");
    verbose = true,
    resources = [],
    libraries = [ Libui.libui ],
    appname="HelloLibui", builddir=builddir)


text:
# Functions for measuring model fit
include("DIC.jl")
include("LPML.jl")


text:
abstract type AbstractILMProblem{DT,ST} end

"""
$(TYPEDEF)

When defining a problem type with scalar data, make it a subtype of this.
"""
abstract type AbstractScalarILMProblem{DT,ST} <: AbstractILMProblem{DT,ST} end

"""
$(TYPEDEF)

When defining a problem type with vector data, make it a subtype of this.
"""
abstract type AbstractVectorILMProblem{DT,ST} <: AbstractILMProblem{DT,ST} end

"""
The `@ilmproblem` macro is used to automatically generate a type
particular to an immersed-layer problem, which can then be used for dispatch on
those types of problems. It takes two arguments: the name of the problem
(to which `Problem` will be automatically appended), and whether the problem
is of scalar or vector type. For example,

    @ilmproblem(StokesFlow,vector)

would generate a type `StokesFlowProblem` dealing with vector-valued data. The
resulting type then automatically has a constructor that allows one to pass
in the grid information and bodies, as well as optional choices for the
DDF function and the scaling type. For the example, this constructor would be

    StokesFlowProblem(grid,bodies[,ddftype=CartesianGrids.Yang3][,scaling=IndexScaling])

"""
macro ilmproblem(name,vector_or_scalar)
    vs_string = string(vector_or_scalar)
    abtype = Symbol("Abstract"*uppercasefirst(vs_string)*"ILMProblem")
    typename = Symbol(string(name)*"Problem")
    return esc(quote

          @doc """
              $($typename)

          ILM problem type dealing with $($vs_string)-type data.
          """
          struct $typename{DT,ST,PHT,BCF,FF} <: $abtype{DT,ST}
             g :: PhysicalGrid
             bodies :: BodyList
             phys_params :: PHT
             bc :: BCF
             f_funcs :: FF
             $typename(g::PT,bodies::BodyList;ddftype=CartesianGrids.Yang3,
                                              scaling=IndexScaling,
                                              phys_params=nothing,
                                              bc=nothing,
                                              f_funcs=nothing) where {PT} =
                    new{ddftype,scaling,typeof(phys_params),typeof(bc),typeof(f_funcs)}(g,bodies,phys_params,bc,f_funcs)
             $typename(g::PT,body::Body;ddftype=CartesianGrids.Yang3,
                                        scaling=IndexScaling,
                                        phys_params=nothing,
                                        bc=nothing,
                                        f_funcs=nothing) where {PT} =
                    new{ddftype,scaling,typeof(phys_params),typeof(bc),typeof(f_funcs)}(g,BodyList([body]),phys_params,bc,f_funcs)
          end

     end)
end

@ilmproblem BasicScalarILM scalar
@ilmproblem BasicVectorILM vector

# Extend this function for other problem types in order to create an extra cache
# of variables and operators for the problem
"""
    prob_cache(prob,base_cache::BasicILMCache)

This function is called by [`__init`](@ref) to generate a problem-specific extra
cache. Extend this function in order to generate an extra cache for a
user-defined problem type. The user must define the cache type itself.
""" prob_cache

prob_cache(prob::BasicScalarILMProblem,base_cache::BasicILMCache) = nothing
prob_cache(prob::BasicVectorILMProblem,base_cache::BasicILMCache) = nothing


text:
import Future
export Box

"Mutable wrapper for a value"
mutable struct Box{T}
  val::T
end

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
# Julia wrapper for header: /opt/ffmpeg/include/libavutil/timecode.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_timecode_adjust_ntsc_framenum2,
    av_timecode_get_smpte_from_framenum,
    av_timecode_make_string,
    av_timecode_make_smpte_tc_string,
    av_timecode_make_mpeg_tc_string,
    av_timecode_init,
    av_timecode_init_from_string,
    av_timecode_check_frame_rate

function av_timecode_adjust_ntsc_framenum2(framenum::Integer,fps::Integer)
    ccall((:av_timecode_adjust_ntsc_framenum2,libavutil),Cint,(Cint,Cint),framenum,fps)
end

function av_timecode_get_smpte_from_framenum(tc,framenum::Integer)
    ccall((:av_timecode_get_smpte_from_framenum,libavutil),UInt32,(Ptr{AVTimecode},Cint),tc,framenum)
end

function av_timecode_make_string(tc,buf,framenum::Integer)
    ccall((:av_timecode_make_string,libavutil),Ptr{UInt8},(Ptr{AVTimecode},Ptr{UInt8},Cint),tc,buf,framenum)
end

function av_timecode_make_smpte_tc_string(buf,tcsmpte::Integer,prevent_df::Integer)
    ccall((:av_timecode_make_smpte_tc_string,libavutil),Ptr{UInt8},(Ptr{UInt8},UInt32,Cint),buf,tcsmpte,prevent_df)
end

function av_timecode_make_mpeg_tc_string(buf,tc25bit::Integer)
    ccall((:av_timecode_make_mpeg_tc_string,libavutil),Ptr{UInt8},(Ptr{UInt8},UInt32),buf,tc25bit)
end

function av_timecode_init(tc,rate::AVRational,flags::Integer,frame_start::Integer,log_ctx)
    ccall((:av_timecode_init,libavutil),Cint,(Ptr{AVTimecode},AVRational,Cint,Cint,Ptr{Void}),tc,rate,flags,frame_start,log_ctx)
end

function av_timecode_init_from_string(tc,rate::AVRational,str,log_ctx)
    ccall((:av_timecode_init_from_string,libavutil),Cint,(Ptr{AVTimecode},AVRational,Ptr{UInt8},Ptr{Void}),tc,rate,str,log_ctx)
end

function av_timecode_check_frame_rate(rate::AVRational)
    ccall((:av_timecode_check_frame_rate,libavutil),Cint,(AVRational,),rate)
end

max_stars_repo_name: asafasban/AAVideoIO

max_issues_repo_name: asafasban/AAVideoIO

max_forks_repo_name: asafasban/AAVideoIO


text:
export RandomPolicy

using Random

"""
    RandomPolicy(action_space, rng)

Randomly return a valid action.
"""
struct RandomPolicy{S<:AbstractSpace,R<:AbstractRNG} <: AbstractPolicy
    action_space::S
    rng::R
end

Base.show(io::IO, p::RandomPolicy) = print(io, "RandomPolicy($(p.action_space))")

Random.seed!(p::RandomPolicy, seed) = Random.seed!(p.rng, seed)

"""
    RandomPolicy(action_space; seed=nothing)
"""
RandomPolicy(s; seed = nothing) = RandomPolicy(s, MersenneTwister(seed))

"""
    RandomPolicy(env::AbstractEnv; seed=nothing)
"""
RandomPolicy(env::AbstractEnv; seed = nothing) =
    RandomPolicy(get_action_space(env), MersenneTwister(seed))

function (p::RandomPolicy)(obs, ::FullActionSet)
    legal_actions = get_legal_actions(obs)
    length(legal_actions) == 0 ? get_invalid_action(obs) : rand(p.rng, legal_actions)
end

(p::RandomPolicy)(obs, ::MinimalActionSet) = rand(p.rng, p.action_space)
(p::RandomPolicy)(obs::BatchObs, ::MinimalActionSet) =
    [rand(p.rng, p.action_space) for _ in 1:length(obs)]

RLBase.update!(p::RandomPolicy, experience) = nothing

RLBase.get_prob(p::RandomPolicy, s) =
    fill(1 / length(p.action_space), length(p.action_space))
RLBase.get_prob(p::RandomPolicy, s, a) = 1 / length(p.action_space)


text:
# add GoogleCloud    # pkg mode

using GoogleCloud

creds = GoogleCredentials(expanduser("credentials.json"))
session = GoogleSession(creds, ["devstorage.full_control"])
set_session!(storage, session)    # storage is the API root, exported from GoogleCloud.jl
bkts = storage(:Bucket, :list)    # storage(:Bucket, :list; raw=true) returns addition information

for item in bkts
    display(item)
    println()
end

storage(:Bucket, :insert; data=Dict(:name => "a12345foo"))

# Verify the new bucket exists in the project
bkts = storage(:Bucket, :list)
for item in bkts
    display(item)
    println()
end

storage(:Object, :list, "a12345foo")

# String containing the contents of test_image.jpg. The semi-colon avoids an error caused by printing the returned value.
file_contents = readstring(open("test_image.jpg", "r"));

# Upload
storage(:Object, :insert, "a12345foo";     # Returns metadata about the object
    name="image.jpg",           # Object name is "image.jpg"
    data=file_contents,         # The data being stored on your project
    content_type="image/jpeg"   # The contents are specified to be in JPEG format
)

# Verify that the object is in the bucket
obs = storage(:Object, :list, "a12345foo")    # Ugly print
map(x -> x[:name], obs)                       # Pretty print

s = storage(:Object, :get, "a12345foo", "image.jpg");
s == file_contents    # Verify that the retrieved data is the same as that originally posted

storage(:Object, :delete, "a12345foo", "image.jpg")

# Verify that the bucket is now empty
storage(:Object, :list, "a12345foo")

storage(:Bucket, :delete, "a12345foo")

# Verify that the bucket has been deleted
bkts = storage(:Bucket, :list)
for item in bkts
    display(item)
    println()
end


text:
using MCMC
using Base.Test

# Define immutable class State
immutable State
  mu::Float64
  sig2::Float64
end

@testset "Posterior Means" begin
  srand(1)
  using Distributions

  # Generate data
  const n = 1000
  const (mu,sig2) = (5,2)
  const y = rand(Normal(mu,sqrt(sig2)),n)

  # precomputes
  const ybar = mean(y)
  const (sig2_a, sig2_b) = (2.0, 1.0)
  rig(shape::Float64,rate::Float64) = 1 / rand(Gamma(shape,1/rate))

  # define an update function for State
  function update(state::State)
    const newMu = rand(Normal(ybar,sqrt(state.sig2/n)))
    const newSig2 = rig(sig2_a+n/2, sig2_b + sum((y-newMu).^2)/2)
    State(newMu,newSig2)
  end

  # run gibbs
  @time out = gibbs(State(0,1),update,10000,1000);

  ## post processing
  const postMu = map(o -> o.mu, out)
  const postSig2 = map(o -> o.sig2, out)
  @test abs(mean(postMu)-mu) < .05
  @test abs(mean(postSig2)-sig2) < .05
end

#= Personal tests:
# to test: vulture test.jl julia --color=yes test.jl

include("runtests.jl")
using RCall
R"library(rcommon) # install_github('luairthur/rcommon')"
plotPosts = R"plotPosts"
plotPosts([postMu postSig2])
=#

max_stars_repo_name: luiarthur/mcmcJulia

max_issues_repo_name: luiarthur/mcmcJulia

max_forks_repo_name: luiarthur/mcmcJulia


text:
module PredictMDAPI

include("types.jl")
include("traits.jl")
include("functions.jl")

end


text:
#=
# 2D advection.
=#
if !@isdefined(Femshop)
    include("../Femshop.jl");
    using .Femshop
end
init_femshop("advection2d");

@useLog("advection2dlog")

@domain(2)
@solver(DG)
@functionSpace(LEGENDRE, 3)
@stepper(CRANK_NICHOLSON, 0.1)

@mesh(QUADMESH, 8)

@variable(u)
@testSymbol(v)

T = 2;
@timeInterval(T)
@initial(u, "0")

### Add boundary regions and set BCs ###############
@addBoundaryID(2, x <= 0 || y <= 0)
@addBoundaryID(3, x < 0.2 && y < 0.2)
@boundary(u, 1, NO_BC)
@boundary(u, 2, DIRICHLET, 0)
@boundary(u, 3, DIRICHLET, "sin(pi*(t - 0.7*(x+y)))^2")
####################################################

@coefficient(a, VECTOR, [0.5, 0.5])
@coefficient(beta, 0)
@weakForm(u, "Dt(u*v) + a*grad(u)*v + surface(-dot(a, normal())*u*v + dot(a, normal())*ave(u)*v - 0.5*(1-beta)*dot(a, jump(u))*v)")

solve(u);

# Uncomment to plot
# using Plots
# pyplot();
# display(plot(Femshop.grid_data.allnodes[1,:], Femshop.grid_data.allnodes[2,:], u.values[:], st=:surface, reuse=false))

@finalize()


text:
module Graphs

using StatsBase

export GraphVertex, NodeType, Person, Address,
       generate_random_graph, get_random_person, get_random_address, generate_random_nodes,
       convert_to_graph,
       bfs, check_euler, partition,
       graph_to_str, node_to_str,
       test_graph

 # Types of valid graph node's values.
 abstract NodeType

 type Person <: NodeType
   name :: AbstractString
 end

 type Address <: NodeType
   streetNumber :: Int64
 end

#= Single graph vertex type.
Holds node value and information about adjacent vertices =#

##################################################
############### opt. 1 - typy
##################################################
type GraphVertex{T<:NodeType}
  value :: T
  neighbors ::Vector{GraphVertex}
end

#= Generates random directed graph of size N with K edges
and returns its adjacency matrix.=#
function generate_random_graph(N::Int64, K::Int64)
    A = Array{Int64,2}(N, N)

    for i=1:N, j=1:N
      A[i,j] = 0
    end

    for i in sample(1:N*N, K, replace=false)
      A[i] = 1
    end
    A
end

# Generates random person object (with random name).
function get_random_person()
  Person(randstring())
end

# Generates random person object (with random name).
function get_random_address()
  Address(rand(1:100))
end

# Generates N random nodes (of random NodeType).
function generate_random_nodes(N :: Int64)
  nodes = Array{NodeType,1}()
  for i= 1:N
    push!(nodes, rand() > 0.5 ? get_random_person() : get_random_address())
  end
  nodes
end

#= Converts given adjacency matrix (NxN)
  into list of graph vertices (of type GraphVertex and length N). =#
###############################################
############### opt 1 - typy ##################
###############################################
function convert_to_graph(A::Array{Int64,2}, nodes:: Array{NodeType,1}, graph:: Array{GraphVertex,1})
  N = length(nodes)
  push!(graph, map(n -> GraphVertex(n, GraphVertex[]), nodes)...)

  for i = 1:N, j = 1:N
      if A[i,j] == 1
        push!(graph[i].neighbors, graph[j])
      end
  end
end

#= Groups graph nodes into connected parts. E.g. if entire graph is connected,
  result list will contain only one part with all nodes. =#
function partition(graph:: Array{GraphVertex,1})
  parts = []
  remaining = Set(graph)
  visited = bfs(graph, remaining=remaining)
  push!(parts, Set(visited))

  while !isempty(remaining)
    new_visited = bfs(graph, visited=visited, remaining=remaining)
    push!(parts, new_visited)
  end
  parts
end

#= Performs BFS traversal on the graph and returns list of visited nodes.
  Optionally, BFS can initialized with set of skipped and remaining nodes.
  Start nodes is taken from the set of remaining elements. =#
function bfs(graph:: Array{GraphVertex,1};visited=Set(), remaining=Set())
  if(remaining == Set())
    remaining = Set(graph)
  end
  first = next(remaining, start(remaining))[1]
  q = Array{GraphVertex,1}()
  push!(q, first)
  push!(visited, first)
  delete!(remaining, first)
  local_visited = Set(q)

  while !isempty(q)
    v = pop!(q)
    for n in v.neighbors
      if !(n in visited)
        push!(q, n)
        push!(visited, n)
        push!(local_visited, n)
        delete!(remaining, n)
      end
    end
  end
  local_visited
end

#= Checks if there's Euler cycle in the graph by investigating
   connectivity condition and evaluating if every vertex has even degree =#
function check_euler(graph:: Array{GraphVertex,1})
  if length(partition(graph)) == 1
    return all(map(v -> iseven(length(v.neighbors)), graph))
  end
    "Graph is not connected"
end

##################################################################
################# opt 3. - mniejsze funkcje bez isa ##############
##################################################################
function vertex_to_string(v::GraphVertex{Person})
  return "Person: $(v.value.name)\n"
end

function vertex_to_string(v::GraphVertex{Address})
  return "Street nr: $(v.value.streetNumber)\n"
end

#= Returns text representation of the graph consisiting of each node's value
   text and number of its neighbors. =#
###############################################################
############### opt. 2 - zmienne globalne #####################
###############################################################
function graph_to_str(graph:: Array{GraphVertex,1})
  graph_str = ""
  for v in graph
    graph_str *= "****\n"
    graph_str *= vertex_to_string(v)
    graph_str *= "Neighbors: $(length(v.neighbors))\n"
  end
  graph_str
end

#= Tests graph functions by creating 100 graphs, checking Euler cycle
  and creating text representation. =#
function test_graph()
  for i=1:100

    ###############################################################
    ############### opt. 2 - zmienne globalne #####################
    ###############################################################
    graph = GraphVertex[]
    # Number of graph nodes.
    N = 800

    # Number of graph edges.
    K = 10000

    A = generate_random_graph(N, K)
    nodes = generate_random_nodes(N)
    convert_to_graph(A, nodes, graph)

    str = graph_to_str(graph)
    # println(str)
    println(check_euler(graph))
  end
end

end


text:
"""
```
  genlasso(
      X    :: AbstractMatrix{T},
      y    :: AbstractVector{T};
      path :: Bool = true,
      ρ    :: Union{AbstractVector, Number} = zero(T),
      D    :: AbstractMatrix{T} = eye(size(X, 2)),
      solver = ECOSSolver(maxit=10e8, verbose=0)
    )
```
Solve generalized lasso problem by reformulating it as constrained lasso problem.
Note genralized lasso minimizes
    `0.5sumabs2(√obswt .* (y - X * β)) + ρ * sumabs(D * β)`.

### Arguments
- `X`       : predictor matrix.
- `y`       : response vector.

### Optional arguments
- `path`    : set `path=false` if user wishes to supply parameter value(s).
              Default is true.
- `ρ`       : tuning parameter value(s). Default is 0.
- `D`       : penalty matrix. Default is identity matrix.
- `solver`  : a solver Convex.jl supports. Default is ECOS.
              Note that Mosek and Gurobi are more robust than ECOS. Unlike ECOS or
              SCS, both Mosek and Gurobi require a license (free for academic
              use). <http://convexjl.readthedocs.io/en/latest/solvers.html>

### Returns
- `β`       : estimated coefficents.
- `objval`  : optimal objective value.
- `problem` : Convex.jl problem.
"""
function genlasso(
    X::AbstractMatrix{T},
    y::AbstractArray{T};
    path::Bool = true,
    ρ::Union{AbstractVector, Number} = zero(T),
    D::AbstractMatrix{T} = eye(size(X, 2)),
    solver = ECOSSolver(maxit=10e8, verbose=0)
    ) where T

  # singular value decomposition on D
  m = size(D, 1)
  F = svd!(D, full = true)
  # extract singular values and submatrices
  singvals = F.S
  rankD = count(!iszero,F.S .> abs(F.S[1]) * eps(F.S[1]) * maximum(size(D)))
  V1 = F.V[:, 1:rankD]
  V2 = F.V[:, rankD+1:end]
  U1 = F.U[:, 1:rankD]
  U2 = F.U[:, rankD+1:end]
  # calculate the MP-inverse of D
  Dplus = V1 * broadcast(*, U1', 1 ./ F.S)
  # transform the design matrix
  XDplus = X * Dplus
  XV2 = X * V2
  # projection matrix onto C(XV2)
  Pxv2 = (1 / dot(XV2, XV2)) * XV2 * transpose(XV2)
  # orthogonal projection matrix
  Mxv2 = Matrix{Float64}(I, size(XV2, 1),size(XV2, 1)) - Pxv2
  # obtain the new design matrix and response vector
  ỹ = vec(Mxv2 * y)
  X̃ = Mxv2 * XDplus
  # solve the constrained lasso problem
  if path
    α̂path, ρpath, = lsq_classopath(X̃, ỹ; solver=solver, Aeq = U2',
          beq = zeros(size(U2, 2)))
  else
    α̂path, = lsq_constrsparsereg(X̃, ỹ, ρ; solver=solver, Aeq = U2',
          beq = zeros(size(U, 2)))
    ρpath = ρ
  end
  # transform back to beta
  β̂path = Base.LinAlg.BLAS.ger!(1.0, vec(V2 * ((1 / dot(XV2, XV2)) *
    		At_mul_B(XV2, y))), ones(size(ρpath)), (Matrix{Float64}(I, size(V2, 1),size(V2, 1)) -
    		V2 * ((1 / dot(XV2, XV2)) * At_mul_B(XV2, X))) * Dplus * α̂path)

  return β̂path, ρpath

end # end of function


text:
#### Because charts are rendered in JavaScript, tests are checks for valid syntax
#### Visual correctness will need to be done manually, mostly via docs

using RDatasets
mtcars = dataset("datasets", "mtcars");
diamonds = dataset("ggplot2", "diamonds");

using ECharts, DataFrames, StatsBase, CSV, KernelDensity, Distributions
using Test
using JSON

#1: Homepage doc example
x = ["Monday","Tuesday","Wednesday","Thursday","Friday","Saturday","Sunday"]
y = [11, 11, 15, 13, 12, 13, 10]
ar = area(x, y, color = lineargradient("purple", "cyan"))
smooth!(ar)
xaxis!(ar, name = "Day of Week")
yaxis!(ar, name = "Daily High Temperature °C")
@test typeof(ar) == EChart

#2: area
x = ["Monday","Tuesday","Wednesday","Thursday","Friday","Saturday","Sunday"]
y = [11, 11, 15, 13, 12, 13, 10]
ar2 = area(x, y)
@test typeof(ar2) == EChart

df = DataFrame(x = x, y = y)
@test typeof(area(df, :x, :y)) == EChart
@test typeof(bar(df, :x, :y)) == EChart

x = ["Monday","Tuesday","Wednesday","Thursday","Friday","Saturday","Sunday"]
y = [11, 11, 15, 13, 12, 13, 10]
y2 = 3.7 .* y
as = area(x, hcat(y, y2))
@test typeof(as) == EChart

x = ["Monday","Tuesday","Wednesday","Thursday","Friday","Saturday","Sunday"]
y = [11, 11, 15, 13, 12, 13, 10]
y2 = 3.7 .* y
astep = area(x, hcat(y, y2), step = "middle")
@test typeof(astep) == EChart

x = [0,1,2,3,4,5,6,7,8,9,0,1,2,3,4,5,6,7,8,9]
y = [28, 43, 81, 19, 52, 24, 87, 17, 68, 49, 55, 91, 53, 87, 48, 49, 66, 27, 16, 15]
g = [0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1]
df_merged = DataFrame(hcat(x,y,g), [:x, :y, :g])
adfg = area(df_merged, :x, :y, :g)
@test typeof(adfg) == EChart

x = rand(Beta(3.0, 2.0), 10)
k = kde(x)
@test typeof(area(k)) == EChart

#3: bar
x = ["Monday","Tuesday","Wednesday","Thursday","Friday","Saturday","Sunday"]
y = [11, 11, 15, 13, 12, 13, 10]
b = bar(x, y)
@test typeof(b) == EChart

x = ["Monday","Tuesday","Wednesday","Thursday","Friday","Saturday","Sunday"]
y = [11, 11, 15, 13, 12, 13, 10]
bc = bar(x, y, color = "green")
@test typeof(bc) == EChart

x = ["Monday","Tuesday","Wednesday","Thursday","Friday","Saturday","Sunday"]
y = [11, 11, 15, 13, 12, 13, 10]
bch = bar(x, y, color = "lightgray", horizontal = true)
@test typeof(bch) == EChart

x = ["Monday","Tuesday","Wednesday","Thursday","Friday","Saturday","Sunday"]
y = [11, 11, 15, 13, 12, 13, 10]
bm = bar(x, hcat(0.95 .* y, 1.25 .* y, y))
@test typeof(bm) == EChart

x = ["Monday","Tuesday","Wednesday","Thursday","Friday","Saturday","Sunday"]
y = [11, 11, 15, 13, 12, 13, 10]
bcs = bar(x, hcat(0.95 .* y, 1.25 .* y, y), color = ["red", "gray", "blue"], stack = true)
@test typeof(bcs) == EChart

x = ["Monday","Tuesday","Wednesday","Thursday","Friday","Saturday","Sunday"]
y = [11, 11, 15, 13, 12, 13, 10]
bcsa = bar(x, hcat(0.95 .* y, 1.25 .* y, y), color = ["red", "gray", "blue"], stack = [1,1,2])
@test typeof(bcsa) == EChart

x = [0,1,2,3,4,5,6,7,8,9,0,1,2,3,4,5,6,7,8,9]
y = [28, 43, 81, 19, 52, 24, 87, 17, 68, 49, 55, 91, 53, 87, 48, 49, 66, 27, 16, 15]
g = [0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1]
df_merged = DataFrame(hcat(x,y,g), [:x, :y, :g])
bg = bar(df_merged, :x, :y, :g)
@test typeof(bg) == EChart

#4: box
data = [850, 740, 900, 1070, 930, 850, 950, 980, 980, 880, 1000, 980, 930, 650, 760, 810, 1000, 1000, 960, 960]
b = box(data)
@test typeof(b) == EChart

data = [
    [850, 740, 900, 1070, 930, 850, 950, 980, 980, 880, 1000, 980, 930, 650, 760, 810, 1000, 1000, 960, 960],
    [960, 940, 960, 940, 880, 800, 850, 880, 900, 840, 830, 790, 810, 880, 880, 830, 800, 790, 760, 800],
    [880, 880, 880, 860, 720, 720, 620, 860, 970, 950, 880, 910, 850, 870, 840, 840, 850, 840, 840, 840],
    [890, 810, 810, 820, 800, 770, 760, 740, 750, 760, 910, 920, 890, 860, 880, 720, 840, 850, 850, 780],
    [890, 840, 780, 810, 760, 810, 790, 810, 820, 850, 870, 870, 810, 740, 810, 940, 950, 800, 810, 870]
]

b2 = box(data, horizontal = true)
@test typeof(b2) == EChart

bdf = box(mtcars, :MPG, :Cyl)
@test typeof(bdf) == EChart

#5: candlestick
dt = ["2016/08/26", "2016/08/29", "2016/08/30", "2016/08/31", "2016/09/01", "2016/09/02",
"2016/09/06", "2016/09/07", "2016/09/08", "2016/09/09", "2016/09/12", "2016/09/13",
"2016/09/14", "2016/09/15", "2016/09/16", "2016/09/19", "2016/09/20", "2016/09/21",
"2016/09/22", "2016/09/23", "2016/09/26"]

close_ = [12.3800, 12.4700, 12.5500, 12.6000, 12.4400, 12.5000, 12.6700, 12.7000, 12.7300,
12.3800, 12.7000, 12.3800, 12.1400, 12.1100, 12.1100, 12.1100, 12.0000, 12.0900, 12.1800,
12.1700, 12.0100]

open_ = [12.4700, 12.3800, 12.4600, 12.4800, 12.6600, 12.5300, 12.4900, 12.6200, 12.5000,
12.6100, 12.3200, 12.5300, 12.1800, 12.1400, 12.0500, 12.1200, 12.0900, 12.0800, 12.1700,
12.1200, 12.1200]

high_ = [12.5500, 12.5000, 12.5600, 12.6100, 12.7200, 12.5700, 12.6700, 12.7500, 12.7500,
12.6850, 12.7700, 12.6800, 12.3100, 12.1800, 12.1300, 12.3300, 12.1900, 12.1100, 12.2900,
12.2200, 12.1200]

low_ = [12.3400, 12.3800, 12.4300, 12.4800, 12.3500, 12.4600, 12.4300, 12.6200, 12.5000,
12.3800, 12.2800, 12.3300,12.1100, 12.0600, 12.0100, 12.0586, 11.9600, 12.0100, 12.1600,
12.1200, 12.0000]

c = candlestick(dt, open_, close_, low_, high_)
@test typeof(c) == EChart

#6: funnel
v = [100,60,80,70,50]
n = ["A", "B", "C", "D", "E"]
fn = funnel(n, v)
@test typeof(fn) == EChart

#7: gauge
gg = gauge(27.64)
@test typeof(gg) == EChart

gg2 = gauge(27.64, colors = ["lightgreen", "lightblue","red"], breakpoints = [.4, .7, 1])
@test typeof(gg2) == EChart

#8: line
x = ["Monday","Tuesday","Wednesday","Thursday","Friday","Saturday","Sunday"]
y = [11, 11, 15, 13, 12, 13, 10]
l = line(x, y)
@test typeof(l) == EChart

x = ["Monday","Tuesday","Wednesday","Thursday","Friday","Saturday","Sunday"]
y = [11, 11, 15, 13, 12, 13, 10]
ls = line(x, y, step = "end")
@test typeof(ls) == EChart

x = ["Monday","Tuesday","Wednesday","Thursday","Friday","Saturday","Sunday"]
y = [11, 11, 15, 13, 12, 13, 10]
y2 = 0.6 .* y
l2s = line(x, hcat(y, y2))
@test typeof(l2s) == EChart

#9: pie
x = ["Monday","Tuesday","Wednesday","Thursday","Friday","Saturday","Sunday"]
y = [11, 11, 15, 13, 12, 13, 10]
p = pie(x, y)
@test typeof(p) == EChart

values = [335, 310, 234, 135, 1548]
names = ["direct traffic", "email marketing", "advertising alliance", "video ads", "search engine"]
explode = [true, false, false, false, false]
pe = pie(names, values, selected = explode)
@test typeof(pe) == EChart

x = ["Monday","Tuesday","Wednesday","Thursday","Friday","Saturday","Sunday"]
y = [11, 11, 15, 13, 12, 13, 10]
d = donut(x, y)
@test typeof(d) == EChart

values = [335, 310, 234, 135, 1548]
names = ["direct traffic", "email marketing", "advertising alliance", "video ads", "search engine"]
rose = pie(names, values, roseType = "angle")
@test typeof(rose) == EChart

values = [335, 310, 234, 135, 1548]
names = ["direct traffic", "email marketing", "advertising alliance", "video ads", "search engine"]
radiusd = donut(names, values, roseType = "radius")
@test typeof(radiusd) == EChart

#10: polar
t = 0:360
angle = [a / 180 * pi for a in t]
radius = [(sin(2t) * cos(2t)) for t in angle]
p = polar(angle, radius)
@test typeof(p) == EChart

t = 0:360
angle = [a / 180 * pi for a in t]
radius = [(sin(2t) * cos(2t)) for t in angle]
radius2 = [(sin(2.5t) * cos(2.5t)) for t in angle]
p2 = polar(angle, hcat(radius, radius2))
@test typeof(p2) == EChart

#11: radar
names = ["sales", "administration", "information technology", "customer support", "development",
"marketing"]
max = [6500, 16000, 30000, 38000, 52000, 25000]
values = [4300, 10000, 28000, 35000, 50000, 19000]
r = radar(names, values, max)
@test typeof(r) == EChart

names = ["sales", "administration", "information technology", "customer support", "development",
"marketing"]
max = [6500, 16000, 30000, 38000, 52000, 25000]
values = [4300, 10000, 28000, 35000, 50000, 19000]
values2 = 0.7 .* values
r2 = radar(names, hcat(values, values2), max)
@test typeof(r2) == EChart

#12: sankey
name = ["Agricultural waste", "Bio-conversion", "Liquid", "Losses", "Solid", "Gas", "Biofuel imports",
    "Biomass imports", "Coal imports", "Coal"]

source = [0, 1, 1, 1, 1, 6, 7, 8, 10]
target = [1, 2, 3, 4, 5, 2, 4, 9, 9]
value = [124.729, 0.597, 26.862, 280.322, 81.144, 35, 35, 11.606, 63.965]
s = sankey(name, source, target, value)
@test typeof(s) == EChart

#13: scatter
sc = scatter(rand(30), rand(30))
@test typeof(sc) == EChart

sc2 = scatter(rand(30), hcat(rand(30), rand(30)))
@test typeof(sc2) == EChart

sgrp = scatter(mtcars, :MPG, :HP, :Cyl)
@test typeof(sgrp) == EChart

#14: xy
x = ["Monday","Tuesday","Wednesday","Thursday","Friday","Saturday","Sunday"]
y = [11, 11, 15, 13, 12, 13, 10]
bmf = bar(x, hcat(0.95 .* y, 1.25 .* y, y), mark = ["line", "bar", "bar"])
@test typeof(bmf) == EChart

#15: waterfall
x = 1:5
y = [2900, -1200, -300, -200, -900]
w = waterfall(x, y)
@test typeof(w) == EChart

#16: histogram
h = fit(Histogram, diamonds[:Price], closed = :left)
hs = histogram(h)
@test typeof(hs) == EChart

h2 = fit(Histogram, (diamonds[:Price], diamonds[:Carat]), closed = :left)
hs2 = histogram(h2)
@test typeof(hs2) == EChart

#17: radialbar
x = ["Monday","Tuesday","Wednesday","Thursday","Friday","Saturday","Sunday"]
y = [11, 11, 15, 13, 12, 13, 10]
b = radialbar(x, y)
@test typeof(b) == EChart

x = ["Monday","Tuesday","Wednesday","Thursday","Friday","Saturday","Sunday"]
y = [11, 11, 15, 13, 12, 13, 10]
b = radialbar(x, hcat(0.95 .* y, 1.25 .* y,y))
@test typeof(b) == EChart

x = ["Monday","Tuesday","Wednesday","Thursday","Friday","Saturday","Sunday"]
y = [11, 11, 15, 13, 12, 13, 10]
b = radialbar(x, hcat(0.95 .* y, 1.25 .* y,y), stack = [1,1,2])
@test typeof(b) == EChart

#18: streamgraph
s_df = CSV.read(joinpath(dirname(@__FILE__), "..", "exampledata/streamdata.csv"), types=[String, Float64, String])
sg = streamgraph(s_df[:date], s_df[:value], s_df[:key], legend = true)
@test typeof(sg) == EChart

#19: bubble
sgrp = bubble(mtcars, :MPG, :HP, :Disp)
@test typeof(sgrp) == EChart

sgrp2 = bubble(mtcars[:MPG], mtcars[:HP], mtcars[:Disp])
@test typeof(sgrp2) == EChart

sgrp3 = bubble(mtcars, :MPG, :HP, :Disp, :Cyl)
@test typeof(sgrp3) == EChart

#19: corrplot
cplot = corrplot(mtcars)
@test typeof(cplot) == EChart

df_num = mtcars[[x <: Union{Number, Missing} for x in eltypes(mtcars)]]
c = cor(convert(Matrix, df_num))
ccc = corrplot(c)
@test typeof(ccc) == EChart


text:
module HonlyHelperTurretController
using ..Ahorn, Maple

@mapdef Trigger "HonlyHelper/TurretController" TurretController(x::Integer, y::Integer, width::Integer=Maple.defaultTriggerWidth, height::Integer=Maple.defaultTriggerHeight, turretID::String="TurretID", turretAction::String="Heli Fade In")

const placements = Ahorn.PlacementDict(
    "Turret Controller (HonlyHelper, cursed)" => Ahorn.EntityPlacement(
        TurretController,
        "rectangle",
    ),
)

const turretActions = String["HeliFadeIn", "HeliOn", "HeliLeave", "GunOnlyOn", "GunOnlyOff"]

Ahorn.editingOptions(trigger::TurretController) = Dict{String, Any}(
  "turretAction" => turretActions
)

end

max_stars_repo_name: HollyMagala/HonlyHelper

max_issues_repo_name: HollyMagala/HonlyHelper

max_forks_repo_name: HollyMagala/HonlyHelper


text:
#
# Functions pertaining to instantiated models and their components
#

"""
    modeldef(mi)

Return the `ModelDef` contained by ModelInstance `mi`.
"""
modeldef(mi::ModelInstance) = mi.md

compinstance(mi::ModelInstance, name::Symbol) = mi.components[name]

compdef(ci::ComponentInstance) = compdef(ci.comp_id)

"""
    name(ci::ComponentInstance)

Return the name of the component `ci`.
"""
name(ci::ComponentInstance) = ci.comp_name

"""
    components(mi::ModelInstance)

Return an iterator on the components in model instance `mi`.
"""
components(mi::ModelInstance) = values(mi.components)

"""
    add_comp!(mi::ModelInstance, ci::ComponentInstance)

Add the component `ci` to the `ModelInstance` `mi`'s list of components, and add
the `first` and `last` of `mi` to the ends of the `firsts` and `lasts` lists of
`mi`, respectively.
"""
function add_comp!(mi::ModelInstance, ci::ComponentInstance)
    mi.components[name(ci)] = ci

    push!(mi.firsts, ci.first)
    push!(mi.lasts, ci.last)
end

# Setting/getting parameter and variable values

# Get the object stored for the given variable, not the value of the variable.
# This is used in the model building process to connect internal parameters.
@inline function get_property_obj(obj::ComponentInstanceParameters{NT}, name::Symbol) where {NT}
    return getproperty(nt(obj), name)
end

@inline function get_property_obj(obj::ComponentInstanceVariables{NT}, name::Symbol) where {NT}
    return getproperty(nt(obj), name)
end

@inline function _get_prop(nt::NT, name::Symbol) where {NT <: NamedTuple}
    obj = getproperty(nt, name)
    return fieldtype(NT, name) <: ScalarModelParameter ? obj.value : obj
end

@inline function Base.getproperty(obj::ComponentInstanceParameters{NT}, name::Symbol) where {NT}
    return _get_prop(nt(obj), name)
end

@inline function Base.getproperty(obj::ComponentInstanceVariables{NT}, name::Symbol) where {NT}
    return _get_prop(nt(obj), name)
end

@inline function Base.setproperty!(obj::ComponentInstanceParameters{NT}, name::Symbol, value::VTYPE) where {NT, VTYPE}
    prop_obj = get_property_obj(obj, name)
    T = fieldtype(NT, name)

    if T <: ScalarModelParameter
        return setproperty!(prop_obj, :value, value)
    else
        error("You cannot override indexed parameter $name::$T.")
    end
end

@inline function Base.setproperty!(obj::ComponentInstanceVariables{NT}, name::Symbol, value::VTYPE) where {NT, VTYPE}
    prop_obj = get_property_obj(obj, name)
    T = fieldtype(NT, name)

    if T <: ScalarModelParameter
        return setproperty!(prop_obj, :value, value)
    else
        error("You cannot override indexed variable $name::$T.")
    end
end

"""
    get_param_value(ci::ComponentInstance, name::Symbol)

Return the value of parameter `name` in component `ci`.
"""
function get_param_value(ci::ComponentInstance, name::Symbol)
    try
        return getproperty(ci.parameters, name)
    catch err
        if isa(err, KeyError)
            error("Component $(ci.comp_id) has no parameter named $name")
        else
            rethrow(err)
        end
    end
end

"""
    get_var_value(ci::ComponentInstance, name::Symbol)

Return the value of variable `name` in component `ci`.
"""
function get_var_value(ci::ComponentInstance, name::Symbol)
    try
        # println("Getting $name from $(ci.variables)")
        return getproperty(ci.variables, name)
    catch err
        if isa(err, KeyError)
            error("Component $(ci.comp_id) has no variable named $name")
        else
            rethrow(err)
        end
    end
end

set_param_value(ci::ComponentInstance, name::Symbol, value) = setproperty!(ci.parameters, name, value)

set_var_value(ci::ComponentInstance, name::Symbol, value)  = setproperty!(ci.variables, name, value)

# Allow values to be obtained from either parameter type using one method name.
value(param::ScalarModelParameter) = param.value

value(param::ArrayModelParameter)  = param.values

dimensions(obj::ArrayModelParameter) = obj.dimensions

dimensions(obj::ScalarModelParameter) = []

"""
    variables(mi::ModelInstance, comp_name::Symbol)

Return the `ComponentInstanceVariables` for `comp_name` in ModelInstance 'mi'.
"""
variables(mi::ModelInstance, comp_name::Symbol) = variables(compinstance(mi, comp_name))

variables(ci::ComponentInstance) = ci.variables

"""
    parameters(mi::ModelInstance, comp_name::Symbol)

Return the `ComponentInstanceParameters` for `comp_name` in ModelInstance 'mi'.
"""
parameters(mi::ModelInstance, comp_name::Symbol) = parameters(compinstance(mi, comp_name))

"""
    parameters(ci::ComponentInstance)

Return an iterable over the parameters in `ci`.
"""
parameters(ci::ComponentInstance) = ci.parameters

function Base.getindex(mi::ModelInstance, comp_name::Symbol, datum_name::Symbol)
    if !(comp_name in keys(mi.components))
        error("Component :$comp_name does not exist in current model")
    end

    comp_inst = compinstance(mi, comp_name)
    vars = comp_inst.variables
    pars = comp_inst.parameters

    if datum_name in names(vars)
        which = vars
    elseif datum_name in names(pars)
        which = pars
    else
        error("$datum_name is not a parameter or a variable in component $comp_name.")
    end

    value = getproperty(which, datum_name)

    return value isa TimestepArray ? value.data : value
end

"""
    dim_count(mi::ModelInstance, dim_name::Symbol)

Return the size of index `dim_name`` in model instance `mi`.
"""
dim_count(mi::ModelInstance, dim_name::Symbol) = dim_count(mi.md, dim_name)

dim_key_dict(mi::ModelInstance) = dim_key_dict(mi.md)

dim_value_dict(mi::ModelInstance) = dim_value_dict(mi.md)

function make_clock(mi::ModelInstance, ntimesteps, time_keys::Vector{Int})
    last  = time_keys[min(length(time_keys), ntimesteps)]

    if isuniform(time_keys)
        first, stepsize = first_and_step(time_keys)
        return Clock{FixedTimestep}(first, stepsize, last)

    else
        last_index = findfirst(isequal(last), time_keys)
        times = (time_keys[1:last_index]...,)
        return Clock{VariableTimestep}(times)
    end
end

function reset_variables(ci::ComponentInstance)
    # println("reset_variables($(ci.comp_id))")
    vars = ci.variables

    for (name, T) in zip(names(vars), types(vars))
        value = getproperty(vars, name)

        if (T <: AbstractArray || T <: TimestepArray) && eltype(value) <: AbstractFloat
            fill!(value, NaN)

        elseif T <: AbstractFloat || (T <: ScalarModelParameter && T.parameters[1] <: AbstractFloat)
            setproperty!(vars, name, NaN)

        elseif (T <: ScalarModelParameter)    # integer or bool
            setproperty!(vars, name, 0)
        end
    end
end

function init(mi::ModelInstance)
    for ci in components(mi)
        init(ci)
    end
end

function init(ci::ComponentInstance)
    reset_variables(ci)
    if ci.init !== nothing
        ci.init(ci.parameters, ci.variables, DimDict(ci.dim_dict))
    end
end

function run_timestep(ci::ComponentInstance, clock::Clock)
    if ci.run_timestep === nothing
        return
    end

    pars = ci.parameters
    vars = ci.variables
    dims = ci.dim_dict
    t = clock.ts

    ci.run_timestep(pars, vars, DimDict(dims), t)
    advance(clock)
    nothing
end

function _run_components(mi::ModelInstance, clock::Clock,
                         firsts::Vector{Int}, lasts::Vector{Int}, comp_clocks::Vector{Clock{T}}) where {T <: AbstractTimestep}
    comp_instances = components(mi)
    tups = collect(zip(comp_instances, firsts, lasts, comp_clocks))

    while ! finished(clock)
        for (ci, first, last, comp_clock) in tups
            if first <= gettime(clock) <= last
                run_timestep(ci, comp_clock)
            end
        end
        advance(clock)
    end
    nothing
end

function Base.run(mi::ModelInstance, ntimesteps::Int=typemax(Int),
                  dimkeys::Union{Nothing, Dict{Symbol, Vector{T} where T <: DimensionKeyTypes}}=nothing)
    if length(mi.components) == 0
        error("Cannot run the model: no components have been created.")
    end

    t::Vector{Int} = dimkeys === nothing ? dim_keys(mi.md, :time) : dimkeys[:time]

    firsts = mi.firsts
    lasts = mi.lasts

    if isuniform(t)
        _, stepsize = first_and_step(t)
        comp_clocks = [Clock{FixedTimestep}(first, stepsize, last) for (first, last) in zip(firsts, lasts)]
    else
        comp_clocks = Array{Clock{VariableTimestep}}(undef, length(firsts))
        for i = 1:length(firsts)
            first_index = findfirst(isequal(firsts[i]), t)
            last_index  = findfirst(isequal(lasts[i]), t)
            times = (t[first_index:last_index]...,)
            comp_clocks[i] = Clock{VariableTimestep}(times)
        end
    end

    clock = make_clock(mi, ntimesteps, t)

    init(mi)    # call module's (or fallback) init function

    _run_components(mi, clock, firsts, lasts, comp_clocks)
    nothing
end