text:
using KoalaTransforms
using Koala
using Test
using DataFrames

const EPS = eps(Float64)

## `UnivariateStandardizer`

t = UnivariateStandardizer()
tM = Machine(t, [0, 2, 4])

@test round.(Int, transform(tM, [0,4,8])) == [-1.0,1.0,3.0]
@test round.(Int, inverse_transform(tM, [-1, 1, 3])) == [0, 4, 8]

# create skewed non-negative vector with a zero value:
v = abs.(randn(1000))
v = v .- minimum(v)
KoalaTransforms.normality(v)

t = UnivariateBoxCoxTransformer(shift=true)
tM = Machine(t, v)
@test sum(abs.(v - inverse_transform(tM,transform(tM, v)))) <= 5000*EPS

const X, y = load_ames();
const train, test = partition(eachindex(y), 0.9);

# introduce a field of type `Char`:
X[:OverallQual] = map(Char, X[:OverallQual])

transformer = ToIntTransformer(sorted=true)
transformerM = Machine(transformer, X[:Neighborhood])
v = transform(transformerM, X[test,:Neighborhood])
@test X[test, :Neighborhood] == inverse_transform(transformerM, v)

transformer.map_unseen_to_minus_one = true
transformerM = Machine(transformer, [1,2,3,4])
@test transform(transformerM, 5) == -1
@test transform(transformerM, [5,1])[1] == -1

transformer = DataFrameToArrayTransformer(boxcox=true, standardize=true)
transformerM = Machine(transformer, X)
transform(transformerM, X);

transformer = RegressionTargetTransformer(standardize=false, boxcox=true)
transformerM = Machine(transformer, y)
transform(transformerM, y)

t = Standardizer()
tM = Machine(t, X)
transform(tM, X)

t = Standardizer(features=[:GrLivArea])
tM = Machine(t, X)
transform(tM, X)

t = OneHotEncoder(drop_last=true)
tM = Machine(t, X)
Xt = transform(tM, X)

t = BoxCoxTransformer(shift=true)
tM = Machine(t, X)
transform(tM, X)

t = BoxCoxTransformer(shift=true, features=[:GrLivArea])
tM = Machine(t, X)
Xt = transform(tM, X)

t = DataFrameToArrayTransformer()
tM = Machine(t, X[[:OverallQual]], verbosity=3)
Xt = transform(tM, X)
@test size(Xt, 2) == length(unique(X[:OverallQual]))

t = MakeCategoricalsIntTransformer(initial_label=0, sorted=true)
tM = Machine(t, X)
transform(tM, X)

# IntegerToInt64:
t = IntegerToInt64Transformer()
v = UInt8[4, 5, 2, 3, 1]
tM = Machine(t, v)
@test transform(tM, UInt8[3, 7, 8]) == [3, 7, 8]

# Univariate discretization:
v = randn(10000)
t = UnivariateDiscretizer(n_classes=100)
tM = Machine(t, v)
w = transform(tM, v)
bad_values = filter(v - inverse_transform(tM, w)) do x
    abs(x) > 0.05
end
@test length(bad_values)/length(v) < 0.06

# Discretization of DataFrame:
X = DataFrame(n1=["a", "b", "a"], n2=["g", "g", "g"], n3=[7, 8, 9],
              n4 =UInt8[3,5,10],  o1=[4.5, 3.6, 4.0], )

t = Discretizer(features=[:o1, :n3, :n2, :n1])
tM = Machine(t, X)
Xt = transform(tM, X)
@test Xt.features == [:o1, :n3, :n2, :n1]
@test Xt.is_ordinal == [true, false, false, false]
@test Xt.A == [512 1 1 1; 1 2 1 2; 256 3 1 1]

# testing transformer chains:

u = rand(100)
t = UnivariateBoxCoxTransformer()
tM = Machine(t, u)
v = transform(tM, u)
s = UnivariateStandardizer()
sM = Machine(s, v)
w = transform(sM, v)

st = ChainTransformer(s, t)
stM = Machine(st, u)
@test w ≈ transform(stM, u)

stM[1]
stM[2]


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiCoreV1ReplicationControllerSpec <: SwaggerModel
    minReadySeconds::Any # spec type: Union{ Nothing, Int32 } # spec name: minReadySeconds
    replicas::Any # spec type: Union{ Nothing, Int32 } # spec name: replicas
    selector::Any # spec type: Union{ Nothing, Dict{String, String} } # spec name: selector
    template::Any # spec type: Union{ Nothing, IoK8sApiCoreV1PodTemplateSpec } # spec name: template

    function IoK8sApiCoreV1ReplicationControllerSpec(;minReadySeconds=nothing, replicas=nothing, selector=nothing, template=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1ReplicationControllerSpec, Symbol("minReadySeconds"), minReadySeconds)
        setfield!(o, Symbol("minReadySeconds"), minReadySeconds)
        validate_property(IoK8sApiCoreV1ReplicationControllerSpec, Symbol("replicas"), replicas)
        setfield!(o, Symbol("replicas"), replicas)
        validate_property(IoK8sApiCoreV1ReplicationControllerSpec, Symbol("selector"), selector)
        setfield!(o, Symbol("selector"), selector)
        validate_property(IoK8sApiCoreV1ReplicationControllerSpec, Symbol("template"), template)
        setfield!(o, Symbol("template"), template)
        o
    end
end # type IoK8sApiCoreV1ReplicationControllerSpec

const _property_map_IoK8sApiCoreV1ReplicationControllerSpec = Dict{Symbol,Symbol}(Symbol("minReadySeconds")=>Symbol("minReadySeconds"), Symbol("replicas")=>Symbol("replicas"), Symbol("selector")=>Symbol("selector"), Symbol("template")=>Symbol("template"))
const _property_types_IoK8sApiCoreV1ReplicationControllerSpec = Dict{Symbol,String}(Symbol("minReadySeconds")=>"Int32", Symbol("replicas")=>"Int32", Symbol("selector")=>"Dict{String, String}", Symbol("template")=>"IoK8sApiCoreV1PodTemplateSpec")
Base.propertynames(::Type{ IoK8sApiCoreV1ReplicationControllerSpec }) = collect(keys(_property_map_IoK8sApiCoreV1ReplicationControllerSpec))
Swagger.property_type(::Type{ IoK8sApiCoreV1ReplicationControllerSpec }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiCoreV1ReplicationControllerSpec[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1ReplicationControllerSpec }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1ReplicationControllerSpec[property_name]

function check_required(o::IoK8sApiCoreV1ReplicationControllerSpec)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1ReplicationControllerSpec }, name::Symbol, val)
end

ext: jl


text:
using MccDaqHats
using Revise
includet(joinpath(@__DIR__, "scan_utils.jl"))

#Constants
CURSOR_BACK_2 = "\x1b[2D"
ERASE_TO_END_OF_LINE = "\x1b[0K"

"""
    function finite_scan()

MCC 172 Functions Demonstrated:

    mcc172_iepe_config_write
    mcc172_a_in_clock_config_write
    mcc172_a_in_clock_config_read
    mcc172_a_in_sensitivity_write
    mcc172_a_in_scan_start
    mcc172_a_in_scan_read
    mcc172_a_in_scan_stop

Purpose:
    Perform a finite acquisition on 1 or more channels.

Description:
    Acquires blocks of analog input data for a user-specified group
    of channels.  The RMS value for each channel is displayed for each
    block of data received from the device.  The acquisition is stopped
    when the specified number of samples is acquired for each channel.

Enahancements:
Deinterleave data using deinterleave()
Save data use one of the libraries in github.com/juliaIO, use similar format to matlab format
"""
function finite_scan()
    # Store the channels in a list and convert the list to a channel mask that
    # can be passed as a parameter to the MCC 172 functions.
    channels = [0, 1]
    num_channels = length(channels)

    sensitivity = 1000.0
    samples_per_channel = 10000
    scan_rate = 10240.0
    options = Set([:OPTS_DEFAULT])

    try
        # Select an MCC 172 HAT device to use.
        hat = select_hat_devices(:MCC_172, 1)
        address = hat[1].address

        println("\nSelected MCC 172 HAT device at address $address")

        # Turn on IEPE supply?
        iepe_enable = get_iepe()

        mcc172_open(address)
        for channel in channels
            mcc172_iepe_config_write(address, channel, iepe_enable)
            mcc172_a_in_sensitivity_write(address, channel, sensitivity)
        end

        # Configure the clock and wait for sync to complete.
        mcc172_a_in_clock_config_write(address, :SOURCE_LOCAL, scan_rate)

        synced = false
        actual_scan_rate = 0    # initial outside while loop
        while !synced
            (_source_type, actual_scan_rate, synced) = mcc172_a_in_clock_config_read(address)
            if !synced
                sleep(0.005)
            end
        end

        println("\nMCC 172 continuous scan example")
        println("    Functions demonstrated:")
        println("         mcc172_iepe_config_write")
        println("         mcc172_a_in_clock_config_write")
        println("         mcc172_a_in_clock_config_read")
        println("         mcc172_a_in_sensitivity_write")
        println("         mcc172_a_in_scan_start")
        println("         mcc172_a_in_scan_read")
        println("         mcc172_a_in_scan_stop")
        println("         mcc172_a_in_scan_cleanup")
        println("    IEPE power: $(iepe_enable ? "on" : "off")")
        println("    Channels: $(join(channels, ", "))")
        println("    Sensitivity: $sensitivity")
        println("    Requested scan rate: $scan_rate")
        println("    Actual scan rate: $(round(actual_scan_rate, digits=3))")
        println("    Samples per channel $(round(samples_per_channel, digits=3))")
        println("    Options: $(options)")

        try
            println("\nPress 'Enter' to continue")
            readline()
        catch
            error("^C to end program")
        end

        # Configure and start the scan.
        channel_mask = chan_list_to_mask(channels)
        mcc172_a_in_scan_start(address, channel_mask, UInt32(samples_per_channel), options)

        println("Starting scan ... Press Ctrl-C to stop\n")

        # Display the header row for the data table.
        print("Samples Read    Scan Count")
        for chan in channels
            print("      Ch , chan,  RMS")
        end
        println()

        try
            read_and_display_data(address, samples_per_channel, num_channels)

        catch e # KeyboardInterrupt
            if isa(e, InterruptException)
                # Clear the "^C" from the display.
                println("$CURSOR_BACK_2, $ERASE_TO_END_OF_LINE, Aborted\n")
            else
                println("\n $e")
            end
            mcc172_a_in_scan_stop(address)
        end

        mcc172_a_in_scan_cleanup(address)

    catch err  # (HatError, ValueError) as err
        println("\n $err")
    end
    return nothing
end

"""
    function read_and_display_data(address::Integer, samples_per_channel::Integer, num_channels::Integer)

Reads data from the specified channels on the specified DAQ HAT devices
and updates the data on the terminal display.  The reads are executed in a
loop that continues until either the scan completes or an overrun error
is detected.

Args:

    hat (mcc172): The mcc172 HAT device object.
    samples_per_channel: The number of samples to read for each channel.
    num_channels (int): The number of channels to display.

Returns:

    None
"""
function read_and_display_data(address::Integer, samples_per_channel::Integer, num_channels::Integer)
    total_samples_read = 0
    read_request_size = 1000
    timeout = 5.0

    # Since the read_request_size is set to a specific value, a_in_scan_read()
    # will block until that many samples are available or the timeout is
    # exceeded.

    # Continuously update the display value until Ctrl-C is
    # pressed or the number of samples requested has been read.
    while total_samples_read < samples_per_channel
        resultcode, statuscode, result, samples_read =
            mcc172_a_in_scan_read(address, Int32(read_request_size), num_channels, timeout)

        # Check for an overrun error
        status = mcc172_status_decode(statuscode)
        if status.hardwareoverrun
            println("\n\nHardware overrun\n")
            break
        elseif status.bufferoverrun
            println("\n\nBuffer overrun\n")
            break
        end

        samples_read_per_channel = Integer(length(result) / num_channels)
        total_samples_read += samples_read_per_channel

        print("\r$samples_read_per_channel,         $total_samples_read")

        # Display the RMS voltage for each channel.
        if samples_read_per_channel > 0
            for i in 0:num_channels -1
                # @show(length(result), i, num_channels, samples_read_per_channel)
                value = calc_rms(result, i, num_channels, samples_read_per_channel)
                print("             $i    $(round(value, digits=5))")
            end
            # stdout.flush()
        end
    end
    println("\n")
end

max_stars_repo_name: Spectrum-Tec/MccDaqHats

max_issues_repo_name: Spectrum-Tec/MccDaqHats

max_forks_repo_name: Spectrum-Tec/MccDaqHats


text:
"""
    RFrames

Ability to define frames in a relative fashion.

# Example
```julia
Background(1:100, ground)
Object(1:90, (args...)->circ("red"))
Object(RFrames(10), (args...)->circ("blue"))
Object((args...)->circ("red"))
```
is the same as
```julia
Background(1:100, ground)
Object(1:90, (args...)->circ("red"))
Object(91:100, (args...)->circ("blue"))
Object(91:100, (args...)->circ("red"))
```

# Fields
- frames::UnitRange defines the frames in a relative fashion.
"""
struct RFrames
    frames::UnitRange
end

"""
    RFrames(i::Int)

Shorthand for RFrames(1:i)
"""
RFrames(i::Int) = RFrames(1:i)


text:
function CQL_alpha_loss(π, 𝒫, 𝒟; info = Dict())
    ignore() do
        info["CQL alpha"] = exp(𝒫[:CQL_log_α][1])
    end
    -conservative_loss(π, 𝒫, 𝒟)
end

function importance_sampling(πsamp, π, obs, Nsamples)
    @assert ndims(obs) == 2 # does not support multidimensional observations yet
    @assert critic(π) isa DoubleNetwork # Assumes we have a double network

    rep_obs, flat_actions, logprobs = Zygote.ignore() do
        actions_and_logprobs = [exploration(πsamp, obs) for i=1:Nsamples]
        actions = cat([a for (a, _) in actions_and_logprobs]..., dims=3)
        logprobs = cat([lp for (_, lp) in actions_and_logprobs]..., dims=3)
        rep_obs = repeat(obs, 1, Nsamples)
        flat_actions = reshape(actions, size(actions)[1], :)
        rep_obs, flat_actions, logprobs
    end

    qvals = reshape(mean(value(π, rep_obs, flat_actions)), 1, :, Nsamples)

    return qvals .- logprobs
end

function conservative_loss(π, 𝒫, 𝒟; info=Dict())
    obs = 𝒟[:s]
    acts = 𝒟[:a]
    pol_values = importance_sampling(π, π, obs, 𝒫[:CQL_n_action_samples])
    unif_values = importance_sampling(𝒫[:CQL_is_distribution], π, obs, 𝒫[:CQL_n_action_samples])
    combined = cat(pol_values, unif_values, dims=3)
    lse = logsumexp(combined, dims=3)
    loss = mean(lse) - mean(mean(value(π, obs, acts)))

    β = clamp(exp(𝒫[:CQL_log_α][1]), 0f0, 1f6)
    β * (5f0*loss - 𝒫[:CQL_α_thresh])
end

function CQL_critic_loss(;kwargs...)
    Q2loss = double_Q_loss(;kwargs...)
    (π, 𝒫, 𝒟, y; info=Dict()) -> begin
        loss = Q2loss(π, 𝒫, 𝒟, y, info=info)
        c_loss = conservative_loss(π, 𝒫, 𝒟, info=info)
        loss + c_loss
    end
end

function CQL(;π::ActorCritic{T, DoubleNetwork{ContinuousNetwork, ContinuousNetwork}},
              solver_type=BatchSAC,
              CQL_α::Float32=1f0,
              CQL_is_distribution=DistributionPolicy(product_distribution([Uniform(-1,1) for i=1:dim(action_space(π))[1]])),
              CQL_α_thresh::Float32=10f0,
              CQL_n_action_samples::Int=10,
              CQL_α_opt::NamedTuple=(;),
              a_opt::NamedTuple=(;),
              c_opt::NamedTuple=(;),
              log::NamedTuple=(;),
              kwargs...) where T
    # Fill the parameters
    𝒫 = (CQL_log_α=[Base.log(CQL_α)],
          CQL_is_distribution=CQL_is_distribution,
          CQL_n_action_samples=CQL_n_action_samples,
          CQL_α_thresh=CQL_α_thresh)
    solver_type(;
        π=π,
        𝒫=𝒫,
        log=(;dir = "log/cql", log...),
        param_optimizers=Dict(Flux.params(𝒫[:CQL_log_α]) => TrainingParams(;loss=CQL_alpha_loss, name="CQL_alpha_", CQL_α_opt...)),
        a_opt=a_opt,
        c_opt=(loss=CQL_critic_loss(), name="critic_", c_opt...),
        kwargs...)
end


text:
using GoogleCloud
using Glob
using REPL
using REPL.TerminalMenus

const FIRESTORE_URL = "https://firestore.googleapis.com/v1"

function init(filepath)
    global CREDENTIALS
    CREDENTIALS = JSONCredentials(filepath)
    println("Project Admin SDK key is ", filepath, " now!!")
    println(CREDENTIALS)
end

function init()
    global CREDENTIALS
    try
        filenames = glob("*.json", "firebase")
        menu = RadioMenu(filenames, pagesize=4)
        file = request("Choose your admin sdk key for use:", menu)

        if file != -1
            println("Project Admin SDK key is ", filenames[file], " now!!")
            CREDENTIALS = JSONCredentials(filenames[file])
        else
            println("Menu canceled.")
        end
    catch e
        println("Either the folder firebase doesn't exist or it's empty")
    end
end

function projectid()
    global CREDENTIALS
    isnothing(CREDENTIALS) && throw("Run Firebase.init() to load credentials first")
    CREDENTIALS.project_id
end

function get_token()
    global CREDENTIALS
    isnothing(CREDENTIALS) && throw("Run Firebase.init() to load credentials first")
    session = GoogleSession(CREDENTIALS, ["datastore"])
    GoogleCloud.authorize(session)
end

function auth_header()
    println(CREDENTIALS)
    token = get_token()
    "Authorization" => "$(token[:token_type]) $(token[:access_token])"
end


text:
### A Pluto.jl notebook ###
# v0.17.1

using Markdown
using InteractiveUtils

# ╔═╡ 76f9e6bd-6c96-4441-8f80-dec874ea5a11
begin
	import Pkg
	Pkg.activate(".")
	#Pkg.add(; url="git@github.com:simeonschaub/YAPP.jl", rev="main")
	Pkg.develop(; path="/home/simeon/.julia/dev/YAPP")
	Pkg.develop(; path="..")
	Pkg.add(["AbstractAlgebra", "Transducers", "Dictionaries", "BenchmarkTools", "PlutoUI", "Plots"])
end

# ╔═╡ ee3c9bca-3fa0-11ec-3fc4-55bf8b43afb4
begin
	using SchurPolynomials, AbstractAlgebra, YAPP, Transducers, Dictionaries, LinearAlgebra
	using YAPP: Polynomial, Monomial
	using SchurPolynomials: schur
	using Transducers: Map
	using AbstractAlgebra: Partition
end

# ╔═╡ 2bf97bef-38c3-4140-b44d-b8bb9745673f
using Printf

# ╔═╡ 9e830e0c-2cab-4247-8487-ece4b6f82854
using BenchmarkTools

# ╔═╡ 555e1339-0d12-402b-957c-f0d8d70cbc42
using Plots

# ╔═╡ 480c3e86-ef27-437d-872d-612792f259f9
λ = Partition([2, 1, 1])

# ╔═╡ 5987e81e-b1d7-49cf-8bf1-83a43eece501
md"""
$$\begin{align}
s_\lambda(x_1, \dots, x_n) = \sum_{T \in \mathrm{SSYT}(\lambda)} x^T
= \sum_{T \in \mathrm{SSYT}(\lambda)} x_1^{T_1} \cdots x_n^{T_n}
\end{align}$$
"""

# ╔═╡ 53eee919-2c7b-4c0f-bd86-9e715d2b0129
n = 5

# ╔═╡ 97356e16-b186-4d7a-8dc4-0f062dc0c91f
S = Dict(YoungTableau(λ.part) => schur(n, λ) for λ in Generic.partitions(n))

# ╔═╡ d65b2b91-95cd-4bc5-a50c-60dfb994647a
md"""
$$\begin{align}
s_\lambda(x_1, \dots, x_n) = \sum_{\nu \vdash n} \frac{\chi_\nu^\lambda}{z_\nu} \cdot p_\nu(x_1, \dots, x_n)
\end{align}$$
"""

# ╔═╡ 5ce04038-e4a1-4ab2-bb67-b2e379747e56
powerfun(k::Int, x) = sum(xᵢ -> xᵢ^k, x)

# ╔═╡ b626a35f-6732-41db-a3d9-c53a906bc95a
powerfun(ν::Generic.Partition, x) = prod(k -> powerfun(k, x), ν)

# ╔═╡ c6903995-493e-40d6-bf40-4b1aef56f224
function z(ν::Generic.Partition)
	row_counts = foldl(right, GroupBy(identity, (x, _) -> x + 1, 0), ν)
	prod(pairs(row_counts)) do (k, rₖ)
		factorial(rₖ) * k^rₖ
	end
end

# ╔═╡ d407ded7-6613-4068-b55c-c1869e9178fe
x = Monomial{Vector{Int}}.(eachcol(I(n)))#@polyvar x[1:n]

# ╔═╡ b7c2ff8f-ff71-499b-84fe-d0cc77ee25d3
Base.copy(m::Monomial) = Monomial(copy(m.exponents))

# ╔═╡ d4b75503-c43a-4c58-880b-628f5e9e0972
S′ = Dict(
	YoungTableau(λ.part) => copy(sum(Generic.partitions(n)) do ν
		character(λ, ν) // z(ν) * powerfun(ν, x)
	end, Int)
	for λ in Generic.partitions(n)
)

# ╔═╡ ceebd22f-447c-4310-8a52-b22ec98ce876
S

# ╔═╡ 8a08ba51-f000-4771-854a-301133c7c405
let x = rand(n)
	(Ref(x) .|> Dictionary(S)) .- (Ref(x) .|> Dictionary(S′))
end

# ╔═╡ 9ee2519b-786a-428e-a517-7e599c11f3e4
let m=11, n=6
	x = Monomial{Vector{Int}}.(eachcol(I(n)))
	x_ = randn(n, m)
	p, s = (Matrix{Float64}(undef, length(Generic.partitions(n)), m) for _ in 1:2)
	for (i, λ) in enumerate(Generic.partitions(n))
		p[i, :] .= (eachcol(x_) .|> powerfun(λ, x[1:n])) ./ z(λ)
		s[i, :] .= (eachcol(x_) .|> schur(n, λ))
	end
	s / p .|> x -> round(Int, x)
	Text.(Printf.format.([Printf.format"%.2f"], s / p))
end

# ╔═╡ b3f710d2-5216-47a8-8427-3b0a5d0f2dd7
@benchmark schur(10, Partition([4, 3, 2, 1]))

# ╔═╡ f80625c9-0bc7-410c-a1aa-cbd80829520a
function next_partition!(p)
	i = rand(0:length(p))
	if i == 0
		push!(p, 1)
	elseif get(p, i-1, typemax(Int)) > p[i]
		p[i] += 1
	else
		next_partition!(p)
	end
	p
end

# ╔═╡ 3eef3ce3-c303-41ff-b608-6b5653fc1dc4
p = accumulate(1:12; init=Int[]) do p, _
	next_partition!(copy(p))
end

# ╔═╡ 766ca80a-5918-4d1b-bd6c-d0dba151208c
times = [@belapsed(schur($n, $(Partition(p)))) for (n, p) in enumerate(p)]

# ╔═╡ 5cbe2be9-1308-4bb0-8955-8666ad489110
plt1 = scatter(times; xscale=:log10, yscale=:log10, legend=:none,
	xlabel=raw"$n$", ylabel=raw"Runtime $t$ in s",
	title="Performance Scaling of Schur Polynomial Calculation")

# ╔═╡ 554ea451-6716-4824-adb3-0b95cd10036d
function _schur(n, λ)
	copy(sum(Generic.partitions(n)) do ν
		character(λ, ν) // z(ν) * powerfun(ν, x)
	end, Int)
end

# ╔═╡ c67d3c77-ba3e-4f1d-9410-9aa8e6675490
times2 = [@belapsed(_schur($n, $(Partition(p)))) for (n, p) in enumerate(p)]

# ╔═╡ bd628110-b25d-433d-938d-f993c5ea3d3f
plt2 = scatter([times times2]; xscale=:log10, yscale=:log10,# legend=:none,
	xlabel=raw"$n$", ylabel=raw"Runtime $t$ in s",
	title="Performance Scaling of Schur Polynomial Calculation",
	label=["sum over SSYT" "\$\\chi\$ times power polynomials"])

# ╔═╡ 9f2a664e-39e6-44de-be6d-b581ce5a8787
savefig(plt2, "performance.pdf")

# ╔═╡ 296f1da5-dae2-4e7a-ba19-bc5a683ee7ee
sprint() do io
	for (λ, s_λ) in (λ => schur(4, λ) for λ in Generic.partitions(4))
		print(io, "\\yng(")
		join(io, λ, ",")
		print(io, ") & \$ \\begin{autobreak}")
		print(io, replace(repr(MIME("text/latex"), s_λ), "\$\$" => ""))
		println(io, "\\end{autobreak} \$ \\\\")
	end
end |> Text

# ╔═╡ Cell order:
# ╠═76f9e6bd-6c96-4441-8f80-dec874ea5a11
# ╠═ee3c9bca-3fa0-11ec-3fc4-55bf8b43afb4
# ╠═2bf97bef-38c3-4140-b44d-b8bb9745673f
# ╠═480c3e86-ef27-437d-872d-612792f259f9
# ╟─5987e81e-b1d7-49cf-8bf1-83a43eece501
# ╠═53eee919-2c7b-4c0f-bd86-9e715d2b0129
# ╠═97356e16-b186-4d7a-8dc4-0f062dc0c91f
# ╠═d65b2b91-95cd-4bc5-a50c-60dfb994647a
# ╠═5ce04038-e4a1-4ab2-bb67-b2e379747e56
# ╠═b626a35f-6732-41db-a3d9-c53a906bc95a
# ╠═c6903995-493e-40d6-bf40-4b1aef56f224
# ╠═d407ded7-6613-4068-b55c-c1869e9178fe
# ╠═b7c2ff8f-ff71-499b-84fe-d0cc77ee25d3
# ╠═d4b75503-c43a-4c58-880b-628f5e9e0972
# ╠═ceebd22f-447c-4310-8a52-b22ec98ce876
# ╠═8a08ba51-f000-4771-854a-301133c7c405
# ╠═9ee2519b-786a-428e-a517-7e599c11f3e4
# ╠═9e830e0c-2cab-4247-8487-ece4b6f82854
# ╠═b3f710d2-5216-47a8-8427-3b0a5d0f2dd7
# ╠═f80625c9-0bc7-410c-a1aa-cbd80829520a
# ╠═3eef3ce3-c303-41ff-b608-6b5653fc1dc4
# ╠═766ca80a-5918-4d1b-bd6c-d0dba151208c
# ╠═555e1339-0d12-402b-957c-f0d8d70cbc42
# ╠═5cbe2be9-1308-4bb0-8955-8666ad489110
# ╠═9f2a664e-39e6-44de-be6d-b581ce5a8787
# ╠═554ea451-6716-4824-adb3-0b95cd10036d
# ╠═c67d3c77-ba3e-4f1d-9410-9aa8e6675490
# ╠═bd628110-b25d-433d-938d-f993c5ea3d3f
# ╠═296f1da5-dae2-4e7a-ba19-bc5a683ee7ee


text:
mutable struct GridsContinuumParams
    epsilon::Float64
    num_generative_samples::Int64
    num_bridge_samples::Int64
    horizon_limit::Int64
    vals_along_axis::Int64 # Number of points along axis for interpolation
    move_coefficient::Float64
    time_coefficient::Float64
    beta::Float64
    move_amt::Float64
    move_std::Float64
    cworlds::Dict{Int64,CWorld} # The exemplar continuum worlds
end

function continuum_parse_params(filename::AbstractString)

    params_key = TOML.parsefile(filename)

    return GridsContinuumParams(params_key["EPSILON"],
                                params_key["NUM_GEN_SAMPLES"],
                                params_key["NUM_BRIDGE_SAMPLES"],
                                params_key["HORIZON_LIMIT"],
                                params_key["AXIS_VALS"],
                                params_key["MOVE_COEFFICIENT"],
                                params_key["TIME_COEFFICIENT"],
                                params_key["BETA"],
                                params_key["MOVE_AMT"],
                                params_key["MOVE_STD"],
                                Dict{Int64, CWorld}())

end

const GridsContinuumStateType = CMSSPState{Int64, Vec2}
const GridsContinuumActionType = CMSSPAction{Int64, Vec2}
const GridsContinuumContextType = Dict{Tuple{Int64, Int64}, Vec2}

mutable struct GridsContinuumContextSet
    curr_timestep::Int64
    curr_context::GridsContinuumContextType
    future_contexts::Vector{GridsContinuumContextType}
end

struct ContinuumDummyBookkeeping
    flag::Bool
end
ContinuumDummyBookkeeping() = ContinuumDummyBookkeeping(false)

struct ContinuumDummyMetadata
    flag::Bool
end
Base.zero(::Type{ContinuumDummyMetadata}) = ContinuumDummyMetadata(false)

GridsContinuumContextSet() = GridsContinuumContextSet(0, GridsContinuumContextType(), Vector{GridsContinuumContextType}(undef, 0))

const GridsContinuumCMSSPType = CMSSP{Int64, Vec2, Int64, Vec2, GridsContinuumParams, GridsContinuumContextSet}
const GridsContinuumMDPType = ModalFinHorMDP{Int64, Vec2, Vec2, GridsContinuumParams}
const GridsContinuumInfhorMDPType = ModalInfHorMDP{Int64, Vec2, Vec2, GridsContinuumParams}
const GridsContinuumBridgeSample = BridgeSample{Vec2, ContinuumDummyMetadata}
const GridsContinuumOLV = OpenLoopVertex{Int64, Vec2, Int64, ContinuumDummyMetadata}
const GridsContinuumGraphTracker = GraphTracker{Int64, Vec2, Int64, ContinuumDummyMetadata, ContinuumDummyBookkeeping, RNG} where {RNG <: AbstractRNG}
const GridsContinuumSolverType = HHPCSolver{Int64, Vec2, Int64, ContinuumDummyMetadata, ContinuumDummyBookkeeping, RNG} where {RNG <: AbstractRNG}

# Const mode values
const CONTINUUM_MODES = [1, 2, 3, 4]
const CONTINUUM_GOAL_MODE = 4

# Discrete MDP - 4 states, 2 actions
# 1 -- 2 -- 4
# |         |
# 3 --------|
function get_grids_continuum_mdp()

    T = zeros(4, 2, 4)
    T[1, 1, 2] = 1.0
    T[1, 2, 3] = 1.0
    T[2, 1, 4] = 1.0
    T[3, 1, 4] = 1.0

    # No cost for mode switches
    R = zeros(4, 2)

    return TabularMDP(T, R, 1.0)
end

function get_grids_continuum_actions(params::GridsContinuumParams)

    actions = Vector{GridsContinuumActionType}(undef, 0)

    # First enter the mode actions
    push!(actions, GridsContinuumActionType(1, 1))
    push!(actions, GridsContinuumActionType(2, 2))

    # Then enter the control actions
    # Use actions of CWorld 1
    cworld_actions = params.cworlds[1].actions

    for (i, cwa) in enumerate(cworld_actions)
        idx = i+2
        push!(actions, GridsContinuumActionType(cwa, idx))
    end

    return actions
end

function create_continuum_cmssp(params::GridsContinuumParams, rng::RNG=Random.GLOBAL_RNG) where {RNG <: AbstractRNG}

    actions = get_grids_continuum_actions(params)
    switch_mdp = get_grids_continuum_mdp()
    temp_goal_state = GridsContinuumStateType(CONTINUUM_GOAL_MODE, sample_continuum_state(rng))

    return GridsContinuumCMSSPType(actions, CONTINUUM_MODES, switch_mdp, temp_goal_state, params, GridsContinuumContextSet())
""
end

function POMDPs.isterminal(mdp::GridsContinuumMDPType, relative_state::Vec2)
    return norm(relative_state) < mdp.params.epsilon
end

function POMDPs.isterminal(cmssp::GridsContinuumCMSSPType, state::GridsContinuumStateType)
    return state.mode == CONTINUUM_GOAL_MODE
end

function HHPC.get_relative_state(mdp::Union{GridsContinuumMDPType, GridsContinuumInfhorMDPType},
                                 source::Vec2, target::Vec2)
    return (source - target)
end

function sample_continuum_state(rng::RNG=Random.GLOBAL_RNG) where {RNG <: AbstractRNG}
    d = Distributions.Uniform(0.0, 1.0)
    return Vec2(rand(rng, d), rand(rng, d))
end

function generate_start_state(cmssp::GridsContinuumCMSSPType, rng::RNG=Random.GLOBAL_RNG) where {RNG <: AbstractRNG}
    return GridsContinuumStateType(1, sample_continuum_state(rng))
end

function POMDPs.reward(mdp::GridsContinuumMDPType, state::Vec2, action::Vec2, statep::Vec2)

    # Copy over reward of underlying Continuum World
    delta_state = norm(statep - state)
    rew = -1.0*(mdp.params.move_coefficient*delta_state + mdp.params.time_coefficient)
    return rew

end

function POMDPs.convert_s(::Type{Vector{Float64}}, s::Vec2,
                          mdp::GridsContinuumMDPType) where V <: AbstractVector{Float64}
    return [s[1], s[2]]
end

function POMDPs.generate_s(mdp::GridsContinuumMDPType, state::Vec2, action::Vec2, rng::RNG=Random.GLOBAL_RNG) where {RNG <: AbstractRNG}

    params = mdp.params
    cworld = params.cworlds[mdp.mode]
    return generate_s(cworld, state, action, rng)
end

function POMDPs.generate_sr(mdp::Union{GridsContinuumMDPType, GridsContinuumInfhorMDPType},
                            state::Vec2, action::Vec2, rng::RNG=Random.GLOBAL_RNG) where {RNG <: AbstractRNG}

    statep = generate_s(mdp, state, action, rng)
    rew = reward(mdp, state, action, statep)
    return (statep, rew)
end

function HHPC.expected_reward(mdp::Union{GridsContinuumMDPType, GridsContinuumInfhorMDPType}, state::Vec2,
                              action::Vec2, rng::RNG=Random.GLOBAL_RNG) where {RNG <: AbstractRNG}

    params = mdp.params
    cworld = params.cworlds[mdp.mode]
    avg_reward = 0.0

    for i = 1:params.num_generative_samples
        statep = generate_s(cworld, state, action, rng)
        avg_reward += POMDPs.reward(mdp, state, action, statep)
    end

    avg_reward /= params.num_generative_samples
    return avg_reward
end

# Will never actually be called because any vertices in grid 4 satisfy the goal condition
# Just return the same vertex and the current goal idx
function HHPC.generate_goal_vertex_set!(cmssp::GridsContinuumCMSSPType, vertex::GridsContinuumOLV,
                                        graph_tracker::GridsContinuumGraphTracker, rng::RNG=Random.GLOBAL_RNG) where {RNG <: AbstractRNG}
    return [], [graph_tracker.curr_goal_idx]
end

# Just generate a new set of neighbours at each timestep
# So vertices_to_add has vertices and nbrs_to_add is empty
function HHPC.generate_bridge_vertex_set!(cmssp::GridsContinuumCMSSPType, vertex::GridsContinuumOLV,
                                          mode_pair::Tuple{Int64, Int64}, graph_tracker::GridsContinuumGraphTracker,
                                          action::Int64, rng::RNG=Random.GLOBAL_RNG) where {RNG <: AbstractRNG}

    context_set = cmssp.curr_context_set
    future_contexts = context_set.future_contexts # Vector of contexts
    curr_timestep = context_set.curr_timestep
    params = cmssp.params

    vertices_to_add = Vector{GridsContinuumOLV}(undef, 0)

    @assert !isempty(future_contexts)
    avg_samples = convert(Int64, ceil(cmssp.params.num_bridge_samples/length(future_contexts)))

    # @show avg_samples
    # @show length(future_contexts)

    for (hor, context) in enumerate(future_contexts)

        # Iterate over (M1,M2) -> point
        for (mpair, switch_point) in context

            # Generate bridge point
            if mpair == mode_pair

                d1 = Distributions.Uniform(switch_point[1] - params.epsilon/2.0, switch_point[1] + params.epsilon/2.0)
                d2 = Distributions.Uniform(switch_point[2] - params.epsilon/2.0, switch_point[2] + params.epsilon/2.0)

                for i = 1:avg_samples

                    bp1 = clamp(rand(rng, d1), 0.0, 1.0)
                    bp2 = clamp(rand(rng, d2), 0.0, 1.0)
                    cont_state = Vec2(bp1, bp2)

                    bridge_sample = GridsContinuumBridgeSample(cont_state, cont_state,
                                        TPDistribution([curr_timestep + hor], [1.0]), zero(ContinuumDummyMetadata))
                    new_bridge_vtx = GridsContinuumOLV(mode_pair[1], mode_pair[2], bridge_sample, action)

                    push!(vertices_to_add, new_bridge_vtx)
                end
            end
        end
    end

    # Return all new vertices and no existing neighbours
    return vertices_to_add, []
end

function HHPC.generate_next_valid_modes(cmssp::GridsContinuumCMSSPType, mode::Int64)

    context_set = cmssp.curr_context_set
    future_contexts = context_set.future_contexts

    next_actions_modes = Set{Tuple{Int64, Int64}}()

    for (hor, context) in enumerate(future_contexts) # Vector of dictionaries

        for (mpair, switch_point) in context # Iterate over current dictionary

            if mpair[1] == mode # Switch from current mode

                next_mode = mpair[2]

                for action = 1:2
                    if cmssp.modeswitch_mdp.T[mode, action, next_mode] > 0.0
                        push!(next_actions_modes, (action, next_mode))
                    end
                end
            end
        end
    end

    return next_actions_modes
end

# Dummy function as vertices will always be regenerated
function HHPC.update_vertices_with_context!(cmssp::GridsContinuumCMSSPType, graph_tracker::GridsContinuumGraphTracker,
                                            curr_timestep::Int64)
    return
end

function HHPC.update_next_target!(cmssp::GridsContinuumCMSSPType, solver::GridsContinuumSolverType, curr_timestep::Int64)
    return
end

function set_start_context_set!(cmssp::GridsContinuumCMSSPType,  rng::RNG=Random.GLOBAL_RNG) where {RNG <: AbstractRNG}

    params = cmssp.params
    init_context = GridsContinuumContextType()

    # Generate initial switch point regions
    for m1 in cmssp.modes
        for m2 in cmssp.modes
            for a in cmssp.mode_actions
                if cmssp.modeswitch_mdp.T[m1, a ,m2] > 0.0
                    switch_point = sample_continuum_state(rng)
                    init_context[(m1, m2)] = switch_point
                end
            end
        end
    end

    future_contexts = Vector{GridsContinuumContextType}(undef, 0)
    temp_context = init_context

    # Sample within a radius of twice the epsilon value
    r = 2.0*params.epsilon

    for idx = 1:params.horizon_limit

        next_context = GridsContinuumContextType()

        # Iterate over current switch points and sample around them
        for (mpair, switch_point) in temp_context

            theta = rand(rng, Distributions.Uniform(0.0, 2.0*pi))
            new_pt1 = clamp(switch_point[1] + r*cos(theta), 0.0, 1.0)
            new_pt2 = clamp(switch_point[2] + r*sin(theta), 0.0, 1.0)
            new_point = Vec2(new_pt1, new_pt2)

            next_context[mpair] = new_point
        end

        push!(future_contexts, next_context)
        temp_context = next_context

    end

    cmssp.curr_context_set = GridsContinuumContextSet(0, init_context, future_contexts)
end

function HHPC.get_bridging_action(vertex::GridsContinuumOLV)
    return vertex.bridging_action
end

function HHPC.display_context_future(context_set::GridsContinuumContextSet, future_timestep::Int64)
    @assert future_timestep > context_set.curr_timestep
    println(context_set.future_contexts[future_timestep - context_set.curr_timestep])
end

function HHPC.update_context_set!(cmssp::GridsContinuumCMSSPType, rng::RNG=Random.GLOBAL_RNG) where {RNG <: AbstractRNG}

    curr_context_set = cmssp.curr_context_set
    params = cmssp.params

    # Update timestep
    curr_context_set.curr_timestep += 1

    # First update the current context based on the first future context
    curr_context = curr_context_set.curr_context
    next_context = curr_context_set.future_contexts[1]
    new_curr_context = GridsContinuumContextType()

    for (mpair, curr_switch_point) in curr_context

        next_switch_point = next_context[mpair]

        # Choose an interpolated point between curr and next chosen switch point
        interp_frac = rand(rng, Distributions.Uniform(0.5, 1.0))
        new_switch_point = curr_switch_point + interp_frac*(next_switch_point - curr_switch_point)
        new_curr_context[mpair] = new_switch_point
    end

    # Update new current context
    curr_context_set.curr_context = new_curr_context

    # Pop first future context
    popfirst!(curr_context_set.future_contexts)

    # Add a new future context
    pre_final_context = curr_context_set.future_contexts[end]
    new_final_context = GridsContinuumContextType()

    r = params.epsilon/2.0

    for (mpair, switch_point) in pre_final_context

        theta = rand(rng, Distributions.Uniform(0.0, 2.0*pi))
        new_pt1 = clamp(switch_point[1] + r*cos(theta), 0.0, 1.0)
        new_pt2 = clamp(switch_point[2] + r*sin(theta), 0.0, 1.0)
        new_point = Vec2(new_pt1, new_pt2)

        new_final_context[mpair] = new_point
    end

    push!(curr_context_set.future_contexts, new_final_context)

    @assert length(curr_context_set.future_contexts) == params.horizon_limit

end

function HHPC.simulate_cmssp!(cmssp::GridsContinuumCMSSPType, state::GridsContinuumStateType, a::Union{Nothing, GridsContinuumActionType},
                              curr_timestep::Int64, rng::RNG=Random.GLOBAL_RNG) where {RNG <: AbstractRNG}

    # Default - will delete when timeout removed
    timeout = false

    curr_context_set = cmssp.curr_context_set
    curr_context = curr_context_set.curr_context

    curr_mode = state.mode
    curr_cont_state = state.continuous
    params = cmssp.params

    if typeof(a.action) <: Int64

        # Action is mode switch
        next_mode = findfirst(x -> x > 0.0, cmssp.modeswitch_mdp.T[curr_mode, a.action, :])
        reward = cmssp.modeswitch_mdp.R[curr_mode, a.action]

        # Check if it can be done based on current context
        for (mpair, switch_point) in curr_context

            if mpair == (curr_mode, next_mode)

                if norm(switch_point - curr_cont_state) < 2.0*params.epsilon
                    @info "Successful mode switch to grid ",next_mode
                    next_state = GridsContinuumStateType(next_mode, switch_point)
                    return (next_state, reward, false, timeout)
                else
                    @info "Mode switch to grid ",next_mode, "failed "
                    next_state = state
                    return (next_state, reward, true, timeout)
                end
            end
        end

        @warn "Attempted mode switch to ", next_mode," unavailable from context"
        next_state = curr_state
        return (next_state, reward, true, timeout)
    else
        # control action within mode
        cworld = params.cworlds[curr_mode]
        new_cont_state = generate_s(cworld, state.continuous, a.action, rng)
        reward = POMDPs.reward(cworld, state.continuous, a.action, new_cont_state)
        next_state = GridsContinuumStateType(curr_mode, new_cont_state)
        return (next_state, reward, false, timeout)
    end
end


text:
function solve()
    n = parse(Int, readline())
    if n % 2 == 0
        "White"
    else
        "Black"
    end
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
module ClearSky

using Base: tail
using Base.Threads: @threads, @spawn, Task, fetch
using BasicInterpolators
using Cubature: hquadrature
using Faddeyeva985
using FastGaussQuadrature: gausslegendre, gausslobatto
using Memoize: @memoize
using ProgressMeter: Progress, next!
using ScalarRadau
using UnPack: @unpack

#order shouldn't matter for these
include("constants.jl")
include("util.jl")
include("radiation.jl")
include("orbits.jl")
include("insolation.jl")

#reading HITRAN data and molecular properties
include(joinpath("hitran", "par.jl"))
include(joinpath("hitran", "molparam.jl"))

#computing absorption cross-sections & coefficients
include(joinpath("absorption", "line_shapes.jl"))
include(joinpath("absorption", "gases.jl"))
include(joinpath("absorption", "collision_induced_absorption.jl"))
include(joinpath("absorption", "absorbers.jl"))

#adiabatic profiles, lapse rates, hydrostatics, etc.
include("atmospherics.jl")

#core modeling functions/numerics
include(joinpath("core", "shared.jl"))
include(joinpath("core", "radau.jl"))
include(joinpath("core", "discretized.jl"))

#wrappers of core functions for one shot calculations
include("fluxes.jl")
#highest-level modeling stuff
include("radiative_convective.jl")

end


text:
function colorbar!(fg::FigureGrid; position=:right,
                   vertical=default_isvertical(position), kwargs...)

    guide_pos = guides_position(fg.figure, position)
    return colorbar!(guide_pos, fg; vertical, kwargs...)
end

"""
    colorbar!(figpos, grid; kwargs...)

Compute colorbar for `grid` (which should be the output of [`draw!`](@ref)) and draw it in
position `figpos`. Attributes allowed in `kwargs` are the same as `MakieLayout.Colorbar`.
"""
function colorbar!(figpos, grid; kwargs...)
    colorbar = compute_colorbar(grid)
    return isnothing(colorbar) ? nothing : Colorbar(figpos; colorbar..., kwargs...)
end

function has_zcolor(entry::Entry)
    return entry.plottype <: Union{Heatmap, Contour, Contourf, Surface} &&
        !haskey(entry.primary, :color) &&
        !haskey(entry.named, :color) &&
        !haskey(entry.attributes, :color)
end

function getlabeledcolorrange(grid)
    zcolor = any(has_zcolor, entries(grid))
    key = zcolor ? 3 : :color
    continuous_color_entries = Iterators.filter(entries(grid)) do entry
        return get(entry, key, nothing) isa AbstractArray{<:Number}
    end
    colorrange = compute_extrema(continuous_color_entries, key)
    label = compute_label(continuous_color_entries, key)
    return isnothing(colorrange) ? nothing : (label, colorrange)
end

compute_colorbar(fg::FigureGrid) = compute_colorbar(fg.grid)
function compute_colorbar(grid::Matrix{AxisEntries})
    labeledcolorrange = getlabeledcolorrange(grid)
    isnothing(labeledcolorrange) && return
    label, limits = labeledcolorrange
    colormap = current_default_theme().colormap[]
    for entry in entries(grid)
        colormap = to_value(get(entry.attributes, :colormap, colormap))
    end
    return (; label, limits, colormap)
end


text:
@recipe function f(x::SmoothedFreqTab, fit::TableRegressionModel)
    label --> "Smoothed degree = $(length(coef(fit))-1)"
    x.tab.scale, predict(fit)
end

@recipe function f(x::SmoothedFreqTab)
    label --> "Smoothed observed probability"
    x.tab.scale, x.tab.freq
end

@recipe function f(x::KernelFreqTab)
    label --> "Kernel smoothed observed probability"
    x.tab.scale, x.tab.freq
end

# SG design
@recipe function f(x::SGFreqTab)
    label --> "Marginal Histogram under the SG design"
    seriestype := :heatmap
    # seriescolor --> cgrad([:white, :gray, :black])
    x.marginal
end

@recipe function f(x::NEATEquateResult)
    label --> string(x.method)
    xguide --> "Test X"
    yguide --> "Test Y (scaled)"
    legend --> :topleft
    x.table[!, r"Yx"][!, 1], x.table.scaleX
end


text:
using SIIPExamples
using PowerSystems
using TimeSeries
using Dates

PowerSystems.download(PowerSystems.TestData; branch = "master") # *note* add `force=true` to get a fresh copy
base_dir = dirname(dirname(pathof(PowerSystems)));

RTS_GMLC_DIR = joinpath(base_dir, "data", "RTS_GMLC");
rawsys = PowerSystems.PowerSystemTableData(
    RTS_GMLC_DIR,
    100.0,
    joinpath(RTS_GMLC_DIR, "user_descriptors.yaml"),
    timeseries_metadata_file = joinpath(RTS_GMLC_DIR, "timeseries_pointers.json"),
    generator_mapping_file = joinpath(RTS_GMLC_DIR, "generator_mapping.yaml"),
)

sys = System(rawsys; time_series_resolution = Dates.Hour(1));
horizon = 24;
interval = Dates.Hour(24);
transform_single_time_series!(sys, horizon, interval);
sys

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    precompile(Tuple{typeof(prepare_image),Vector{Vector{UInt}},Vector{UInt},Vector{Int},Dict{UInt,Vector{Base.StackTraces.StackFrame}},Vector{Vector{Base.StackTraces.StackFrame}},Bool,Bool,Bool,Vector{Any}})
    precompile(Tuple{typeof(hcat),Vector{TagData},Vector{TagData},Vector{TagData},Vararg{Vector{TagData},N} where N})
end


text:
module Rivers
"""
Module for the modelling of the dispersion of pollutant in rivers
"""

using ArchGDAL
using ArgParse
using Dates

include(".\\Utils\\Functions.jl")

export run_river

const agd = ArchGDAL

mutable struct River
  ma::Float64  # Pollutant mass
  t::Float64   # Time
  x::Float64
  dl::Float64
  v::Float64   # Speed
  w::Float64   # Hydraulic section
  k::Float64   # Decay coefficient

  concentration

  # element=river(args.concentration,args.time,args.distance,args.fickian,args.velocity)
  River( ma, t, x, dl, v, w, k ) = new(ma,t,x,dl,v,w,k)
end

function calc_concentration!( r::River )
  c1 = r.x - (r.v * r.t)
  c1_1 = -(c1^2)

  c2 = c1_1 / ( 4 * r.dl * r.t )
  c2_1 = exp( -r.k * r.t )
  c3 = exp(c2) * c2_1

  c4 = ( r.ma / r.w ) / ( √( 4π * r.dl * r.t ) )

  r.concentration = c4 * c3
  return r.concentration
end

         #                                                                                min         min_end   min_int
function run_river( dem_file::AbstractString, slope_file::AbstractString, river_file::AbstractString, source_file::AbstractString, start_time, end_time, time_interval,
                  #                    C / conc                radius                                                   w / sez                         k
                    resolution::Int64, concentration::Float64, mean_hydraulic_radius::Float64, fickian_x::Float64=0.05, hydraulic_section::Float64=1.0, decay_coeff::Float64=0.0,
                  # manning
                    manning_coeff::Float64=0.05, output_path::AbstractString )

 """ CONTROLLO SUI RASTER `dem` E `slope`
    if not self.slope.isValid():
        QMessageBox.warning(self,"Warning", u"La mappa delle pendenza è mancante o non valida" )
        return False

    if not self.dem.isValid():
        QMessageBox.warning(self,"Warning", u"Il DEM è mancante o non valido" )
        return False
 """

    src_geom = agd.getgeom(collect(agd.getlayer(agd.read(source_file), 0))[1])

    if agd.geomdim(src_geom) != 0
        throw(DomainError(source_file, "The source shapefile must contain a point."))
    end

   river_layer = agd.getlayer(agd.read(river_file), 0)
   river_geoms = agd.getgeom.(collect(river_layer))

   if any( geom -> agd.geomdim(geom != 1), river_geoms)
      throw(DomainError(river_file, "The river shapefile geometry in not valid."))
   end

   refsys = agd.getspatialref(src_geom)

   dem = agd.read(dem_file)
   slope = agd.read(slope_file)

   if agd.getspatialref(river_layer) != refsys || agd.importWKT(agd.getproj(slope)) != refsys || agd.importWKT(agd.getproj(dem)) != refsys
      throw(DomainError("The reference systems are not uniform. Aborting analysis."))
   end

   if start_time >= end_time
      throw(DomainError("`end_time` must be greater than `start_time`"))
   end

 # messaggio+='ALGORITMO UTILIZZATO: Fickian Mixing Process (Hemond, Harold F., and Elizabeth J. Fechner. Chemical fate and transport in the environment. Elsevier, 2014.)\n\n'

    start_sec, end_sec, int_sec = 60 .* ( time_start, time_end, time_interval )
    #calcolo ciclo intervallo temporale di analisi
    cicles = ( (end_time - start_time) / time_interval ) + 1

 # CREDO SIA L'EQUIVALENTE DI:
    #   feature = next(self.river.getFeatures())
    #   geomfeature = feature.geometry()
    #   features = self.river.getFeatures()
    features = agd.getgeom.(collect(layer))
    x_source = agd.getx(src_geom, 0)
    y_source = agd.gety(src_geom, 0)
    r_source, c_source = toIndexes(dem, x_source, y_source)

 # DEVO TROVARE L'EQUIVALENTE DI INTERPOLATE
    firstpoint = features[1].interpolate(0)

    x_first = agd.getx(firstpoint, 0)
    y_first = agd.gety(firstpoint, 0)
    r_first ,c_first = toIndexes(dem, x_first, y_first)

    slopeband = agd.getband(slope, 1)

    demband = agd.getband(dem, 1)
    demfirstpoint = demband[r_first, c_first]
    demsource = demband[r_source, c_source]

    trend = 0
    old_x = x_first
    old_r = r_first
    old_y = y_first
    old_c = c_first
    if demfirstpoint >= demsource
        trend = 1
        old_x = x_source
        old_r = r_source
        old_y = y_source
        old_c = c_source
    end

    for geom in features
        length = agd.geomlength(geom)
        avanzamento = 1
        realdistance = 0
        feats = []
        featlines=[]
        list_result=[]

        #   start_time = time.time()

        if outputname == ""
            outputname = "concentrazione"
        end

     """ CREA UN LAYER
        vl = QgsVectorLayer("Point?crs=EPSG:"+self.refsys,self.outputname, "memory")
        pr = vl.dataProvider()
        prfield = pr.addAttributes( [ QgsField("distance", QVariant.Int) ] )
        prfield2 = pr.addAttributes( [ QgsField("vmedia", QVariant.Double) ] )
     """
     # NON FUNZIONA
        vl = agd.createlayer( outputname, agd.wkbPoint, refsys )

        list_vmedia = []

     """ CREA UN'ALTRO LAYER
        vline = QgsVectorLayer("LineString?crs=EPSG:"+self.refsys, self.outputname, "memory")
        prline = vline.dataProvider()
        prlfield=prline.addAttributes( [ QgsField("distance", QVariant.Int) ] )
        prlfield2=prline.addAttributes( [ QgsField("vmedia", QVariant.Double) ] )
     """
     # NON FUNZIONA
        vline = agd.createlayer( outputname, agd.wkbLineString, refsys )

        sec_cicli = start_sec
        checknumcampo = 0
        for sec_cicli in start_sec:int_sec:end_sec
            checknumcampo += 1
            fieldname = "conc $(sec_cicli/60)"

         """ AGGIUNGE CAMPI AI LAYER CREATI
            prfield1 = pr.addAttributes( [ QgsField(nomecampo, QVariant.Double) ] )
            prlfield1 = prline.addAttributes( [ QgsField(nomecampo, QVariant.Double) ] )
         """
            agd.addfielddefn!(vl, fieldname, agd.OFTReal )
            agd.addfielddefn!(vline, fieldname, agd.OFTReal )
        end

     """ CREA UNA NUOVA FEATURE E LA AGGIUNGE A vline
        fetline = QgsFeature()
        fetline.setGeometry( QgsGeometry.fromPolyline( [s_geom] ))
        prline.addFeatures( [ fetline ] )
        geomline = fetline.geometry()
     """
        fetline = agd.createfeature( x -> x, vline )
        agd.setgeom!(fetline, agd.wkbMultiLineString)
        agd.addfeature!(vline, fetline)
        geomline = agd.getgeom(fetline)

        controllo = 1
        if trend == 1
            controllo = 0
        end

        for currentdistance in 1:length
            #   point = geom.interpolate(currentdistance)
            point = agd.getpoint(geom, currentdistance)
            x = agd.getx(point, 0)
            y = agd.gety(point, 0)

            dist = √( (x - x_source)^2 + (y - y_source)^2 )
            if trend == 1 && dist <= 1
                controllo = 1
            end
            if trend == 0 && dist <= 1
                controllo = 0
            end

            if controllo == 1
                if avanzamento == resolution
                    realdistance = realdistance + resolution
                    count_index += 1
                    z = slopeband[x, y]
                    v_inst = ( mean_hydraulic_radius^(2/3) * √(z/100) ) * manning_coeff
                    push!( list_vmedia, v_inst )
                    mean_v = sum(list_vmedia) / count_index

                 """ STA INIZIALIZZANDO DELLE FEATURES
                    fet = QgsFeature()
                    fet.initAttributes(2+cicli)

                    fetline = QgsFeature()
                    fetline.initAttributes(2+cicli)

                    fet.setAttribute(0,realdistance)
                    fet.setAttribute(1,self.vmedia)

                    fetline.setAttribute(0,realdistance)
                    fetline.setAttribute(1,self.vmedia)
                 """
                    fet = agd.createfeature(x->x, vl)
                    agd.fillunsetwithdefault!(fet)
                    agd.setfield!( Ref(fet), [0, 1], [realdistance, vmedia])

                    fetline = agd.createfeature(x->x, vline)
                    agd.fillunsetwithdefault!(fetline)
                    agd.setfield!( Ref(fetline), [0, 1], [realdistance, vmedia])

                    for (cicle, t) in enumerate(start_sec:int_sec:end_sec)
                        element = River( concentration, t, realdistance, fickian_x, mean_v, hydraulic_section, decay_coeff )
                        finalC = calc_concentration!(element)
                        if cicle == 1
                            push!( list_result, Cfinal )
                        end

                     """ AGGIORNA I VALORI DELLE FEATURES
                        fet.setAttribute(1+ciclo,Cfinal)
                        fetline.setAttribute(1+ciclo,Cfinal)
                     """
                        agd.setfield!(fet, cicle+1, finalC )
                        agd.setfield!(fetline, cicle+1, finalC )
                    end

                 """ AGGIORNA I VALORI DELLA FEATURE(?) E NE CAMBIA LA GEOMETRIA
                    vl.updateFeature(fet)
                    fet.setGeometry(point)
                 """
                    agd.setgeom!(fet, point)

                 """ COME SOPRA
                    fetline.setGeometry( QgsGeometry.fromPolyline( [QgsPoint(old_x,old_y),QgsPoint(x,y)] ))
                    vline.updateFeature(fetline)
                 """
                    line = agd.createlinestring( Flloat64.([old_x, old_y]), Flloat64.([x, y]) )
                    agd.setgeom!(fetline, line)

                    push!(feats, fet)
                    push!(featlines, fetline)

                    old_x = x
                    old_y = y
                    avanzamento = 0
                end
                avanzamento += 1
            end
        end

        agd.addfeature!.(Ref(vl), feats)
        agd.addfeature!.(Ref(vline), featlines)

     """ AGGIUNGE I VETTORI APPENA CREATI
        QgsProject.instance().addMapLayer(vl)
        QgsProject.instance().addMapLayer(vline)
     """
     """ PIRNT DI COSE
        tempoanalisi=time.time() - start_time
        tempostimato=time.strftime("%H:%M:%S", time.gmtime(tempoanalisi))
        messaggio="---------------------------------\n"
        messaggio+="Fine modellazione\n"
        messaggio+="\nTempo di analisi: "+tempostimato+"\n"
        messaggio+="---------------------------------\n\n"
        self.console.appendPlainText(messaggio)
     """
     """ NON SO COSA FACCIA QUESTA PARTE
        self.tab_2.setEnabled(True)

        self.run_temp_river()
     """
    end
end

end # module

max_stars_repo_name: DavideFavaro/Tirocinio

max_issues_repo_name: DavideFavaro/Tirocinio

max_forks_repo_name: DavideFavaro/Tirocinio


text:
@testset "$(_gen("Optimizer: Poisson Loss [Small][Forgiving][AutoGrad]"))"  begin
    @timeit to "Optimizer: Poisson Loss [Small][Forgiving][AutoGrad]"  begin
        @test unit_test_train_subloss(gradient_eval      = provide(Loss{AbstractPoisson}()),
                                      input_distribution = Distributions.Poisson(10),
                                      input_size         = 1000,
                                      ρ                  = 0,
                                      step_size          = 0.1,
                                      max_iter           = 100) > 0.9
        @test unit_test_train_subloss(gradient_eval      = provide(Loss{AbstractPoisson}()),
                                      input_distribution = Distributions.Poisson(10),
                                      input_size         = 3000,
                                      ρ                  = 0,
                                      step_size          = 0.1,
                                      max_iter          = 100) > 0.9
        @test unit_test_train_subloss(gradient_eval      = provide(Loss{AbstractPoisson}()),
                                      input_distribution = Distributions.Poisson(10),
                                      input_size         = 5000,
                                      ρ                  = 0,
                                      step_size          = 0.1,
                                      max_iter           = 100) > 0.9
        @test unit_test_train_subloss(gradient_eval      = provide(Loss{AbstractPoisson}()),
                                      input_distribution = Distributions.Poisson(10),
                                      input_size         = 1000,
                                      ρ                  = 0.2,
                                      step_size          = 0.1,
                                      max_iter           = 100) > 0.9
        @test unit_test_train_subloss(gradient_eval      = provide(Loss{AbstractPoisson}()),
                                      input_distribution = Distributions.Poisson(10),
                                      input_size         = 3000,
                                      ρ                  = 0.2,
                                      step_size          = 0.1,max_iter = 100) > 0.9
        @test unit_test_train_subloss(gradient_eval      = provide(Loss{AbstractPoisson}()),
                                      input_distribution = Distributions.Poisson(10),
                                      input_size         = 5000,
                                      ρ                  = 0.2,
                                      step_size          = 0.1,
                                      max_iter           = 100) > 0.9
    end
end

@testset "$(_gen("Optimizer: Poisson Loss [Small][Forgiving][Native]"))" begin
    @timeit to "Optimizer: Poisson Loss [Small][Forgiving][Native]" begin
        @test unit_test_train_subloss(gradient_eval      = Loss{AbstractPoisson}(),
                                      input_distribution = Distributions.Poisson(10),
                                      input_size         = 1000,
                                      ρ                  = 0,
                                      step_size          = 0.1,
                                      max_iter           = 100) > 0.9
        @test unit_test_train_subloss(gradient_eval      = Loss{AbstractPoisson}(),
                                      input_distribution = Distributions.Poisson(10),
                                      input_size         = 3000,
                                      ρ                  = 0,
                                      step_size          = 0.1,
                                      max_iter          = 100) > 0.9
        @test unit_test_train_subloss(gradient_eval      = Loss{AbstractPoisson}(),
                                      input_distribution = Distributions.Poisson(10),
                                      input_size         = 5000,
                                      ρ                  = 0,
                                      step_size          = 0.1,
                                      max_iter           = 100) > 0.9
        @test unit_test_train_subloss(gradient_eval      = Loss{AbstractPoisson}(),
                                      input_distribution = Distributions.Poisson(10),
                                      input_size         = 1000,
                                      ρ                  = 0.2,
                                      step_size          = 0.1,
                                      max_iter           = 100) > 0.9
        @test unit_test_train_subloss(gradient_eval      = Loss{AbstractPoisson}(),
                                      input_distribution = Distributions.Poisson(10),
                                      input_size         = 3000,
                                      ρ                  = 0.2,
                                      step_size          = 0.1,max_iter = 100) > 0.9
        @test unit_test_train_subloss(gradient_eval      = Loss{AbstractPoisson}(),
                                      input_distribution = Distributions.Poisson(10),
                                      input_size         = 5000,
                                      ρ                  = 0.2,
                                      step_size          = 0.1,
                                      max_iter           = 100) > 0.9
    end
end


text:
outliar = RGB(0.941176, 0.901961, 0.54902)

writeRGB(opt::Cmd,name::AbstractString) = run(`convert originalRGB.jpg $opt $assets/$name.jpg`)

function updateRGB(opts...)
	_, opt = build_opt(value(SZ),opts...)
	name = string(abs(rand(Int)))
	writeRGB(opt,name)
	name
end

function cleanassets()
	files = readdir(assets)
	filter!(r".*\.jpg",files)
	map(x -> rm(joinpath(assets,x)), files)
end

function copyassets()
	files = readdir(assets)
	filter!(r".*\.jpg",files)
	cp(joinpath(assets,files[1]),"RGB.jpg",remove_destination=true)
end

function recolor_I(img,t)
	N = length(img.cmap)
	n = round(Int,N*t/100)
	cmap = [linspace(RGB(0,0,0),RGB(1,1,1),n);fill(RGB(1,1,1),N - n)]
	img.cmap = cmap
	return convert(Image, img)
end

function recolor_dolp(img,t)
	N = length(img.cmap)
	n = round(Int,N*t/100)
	cmap = [linspace(RGB(0,0,0),RGB(1,1,1),n);fill(outliar,N - n)]
	img.cmap = cmap
	return convert(Image, img)
end

function recolor_docp(img,b,t)
	N2 = round(Int,length(img.cmap)/2)
	nb = round(Int,N2*b/100)
	nt = round(Int,N2*t/100)
	cmap = [fill(outliar,N2 - nb); linspace(RGB(1,0,0),RGB(0,0,0),nb); linspace(RGB(0,0,0),RGB(0,1,0),nt); fill(outliar,N2 - nt)]
	img.cmap = cmap
	return convert(Image, img)
end

write_originalRGB(fname::AbstractString) = run(pipeline(`dcraw -w -H 0 -o 0 -h -c $fname`,`convert - -auto-level originalRGB.jpg`))

function calculatepolar(s)
	mats = stokes2mats(s)
	imgs = mats2images(mats...)
end

function recolor_profile(img,b,t)
	N = length(img.cmap)
	nb = round(Int,N*b/100)
	nt = round(Int,N*t/100)
	cmap = [fill(outliar,nb); linspace(RGB(0,0,0),RGB(1,1,1),nt - nb); fill(outliar,N - nt)]
	img.cmap = cmap
	return img
end

function latexellipse(s,window,maxpol,minpol,eccentricity)
	copyassets()
	stokes2ellipse(s,window_in_cm = window, maxpol = maxpol, minpol = minpol, eccentricity = eccentricity)
	open("ellipse.tex","w") do o
		println(o,"""\\documentclass{standalone}
		\\usepackage{tikz}
		\\usetikzlibrary{arrows.meta,positioning,calc}
		\\usepackage{pgfplots}
		\\pgfplotsset{compat=1.12}
		\\definecolor{dolpcolor}{rgb}{0,1,0}
		\\definecolor{dolcpcolor}{rgb}{1,0,0}
		\\definecolor{dorcpcolor}{rgb}{0,0,1}
		\\begin{document}
				\\input{ellipse.tikz}
		\\end{document}""")
	end
	run(`pdflatex ellipse.tex`)
	name = string(abs(rand(Int)))
	name = "assets/$name.png"
	run(`convert -density 96 -quality 100 ellipse.pdf $name`)
	#run(`convert -density 96 -quality 100 ellipse.pdf -size $(value(SZ)[1])x $name`)
	return name
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "ATKmm"
version = v"2.36.1"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://download.gnome.org/sources/atkmm/$(version.major).$(version.minor)/atkmm-$(version).tar.xz",
                  "e11324bfed1b6e330a02db25cecc145dca03fb0dff47f0710c85e317687da458")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/atkmm*/
mkdir output && cd output/
meson --cross-file=${MESON_TARGET_TOOLCHAIN} ..
ninja -j${nproc}
ninja install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(; experimental=true)

# The products that we will ensure are always built
products = [
    LibraryProduct(["libatkmm-$(version.major).$(version.minor)", "libatkmm-$(version.major)"], :libatkmm)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="ATK_jll", uuid="7b86fcea-f67b-53e1-809c-8f1719c154e8"); compat=string(version))
    Dependency(PackageSpec(name="Glibmm_jll", uuid="5d85a9da-21f7-5855-afec-cdc5039c46e8"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6", preferred_gcc_version = v"7.1.0")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
module ConnectFourSolver

using CxxWrap
using ConnectFourSolver_jll
using Pkg.Artifacts

export Solver, Position, load_book, get_book

@wrapmodule libconnect4jl

function get_book(book="7x6.book")
    book = joinpath(artifact"books", "books", book)
    !isfile(book) && error("book not found: $book")
    book
end

Base.copy(p::Position) = clone(p)

"move starts with 0"
(p::Position)(move::Int) = play(p, move)

"the returned score is negmax"
(solver::Solver)(p::Position, is_weak=false) = solve(solver, p, is_weak)

function __init__()
    @initcxx
end

end # module


text:
"""
$(SIGNATURES)

Fill empty sparse matrix A with random nonzero elements from interval [1,2]
using incremental assembly.

"""
function sprand!(A::AbstractSparseMatrix{Tv,Ti},xnnz::Int) where {Tv,Ti}
    m,n=size(A)
    for i=1:xnnz
        i=rand((1:m))
        j=rand((1:n))
        a=1.0+rand(Tv)
        A[i,j]+=a
    end
end


text:
@doc raw"""
    Torus{N} <: AbstractPowerManifold

The n-dimensional torus is the $n$-dimensional product of the [`Circle`](@ref).

The [`Circle`](@ref) is stored internally within `M.manifold`, such that all functions of
[`AbstractPowerManifold`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/manifolds.html#ManifoldsBase.AbstractPowerManifold)  can be used directly.
"""
struct Torus{N} <: AbstractPowerManifold{ℝ,Circle{ℝ},ArrayPowerRepresentation}
    manifold::Circle{ℝ}
end

Torus(n::Int) = Torus{n}(Circle())

Base.:^(M::Circle, n::Int) = Torus{n}(M)

@doc raw"""
    check_point(M::Torus{n},p)

Checks whether `p` is a valid point on the [`Torus`](@ref) `M`, i.e. each of
its entries is a valid point on the [`Circle`](@ref) and the length of `x` is `n`.
"""
check_point(::Torus, ::Any)
function check_point(M::Torus{N}, p; kwargs...) where {N}
    return check_point(PowerManifold(M.manifold, N), p; kwargs...)
end
@doc raw"""
    check_vector(M::Torus{n}, p, X; kwargs...)

Checks whether `X` is a valid tangent vector to `p` on the [`Torus`](@ref) `M`.
This means, that `p` is valid, that `X` is of correct dimension and elementwise
a tangent vector to the elements of `p` on the [`Circle`](@ref).
"""
function check_vector(M::Torus{N}, p, X; kwargs...) where {N}
    return check_vector(PowerManifold(M.manifold, N), p, X; kwargs...)
end

get_iterator(::Torus{N}) where {N} = 1:N

manifold_dimension(::Torus{N}) where {N} = N

power_dimensions(::Torus{N}) where {N} = (N,)

representation_size(::Torus{N}) where {N} = (N,)

Base.show(io::IO, ::Torus{N}) where {N} = print(io, "Torus($(N))")


text:
using Distributed

# Add this file to following folder:
#-------------------------------------------------------------------------------
# "./.julia/config"
#-------------------------------------------------------------------------------

#-------------------------------------------------------------------------------
addprocs(6)  #ToDo: Modify for the amount of processors you can provide!!!
#-------------------------------------------------------------------------------

#Insert the directory of this project here!
#-------------------------------------------------------------------------------
@everywhere yourdir = "/Desktop/REINFORCEMENT/Master_TRPO_Niko"
#-------------------------------------------------------------------------------

# Include Model Directories
@everywhere begin
	modeldir = homedir()*yourdir
	alldirs = readdir(modeldir)
	for i = 1:length(alldirs)
		alldirs[i] = joinpath(modeldir, alldirs[i], "Julia")

		if isdir(alldirs[i])
			push!(LOAD_PATH, alldirs[i])
		end
	end

	#Include Module Directories
	moduledir = homedir()*yourdir
	alldirs = readdir(moduledir)
	for i = 1:length(alldirs)
		alldirs[i] = joinpath(moduledir, alldirs[i])

		if isdir(alldirs[i])
			push!(LOAD_PATH, alldirs[i])
		end
	end
end

# Set this to have PyPlot as default backend in "Plots.jl"
ENV["PLOTS_DEFAULT_BACKEND"]="PyPlot"

max_stars_repo_name: NikonPic/REINFORCEMENT-JULIA

max_issues_repo_name: NikonPic/REINFORCEMENT-JULIA

max_forks_repo_name: NikonPic/REINFORCEMENT-JULIA


text:
using ..Geometry

import .Rendering.load_font!

export
    VectorFont,
    load_font!, get_glyph

mutable struct VectorGlyph
    vectors::Array{Geometry.Point{Float64},1}

    function VectorGlyph()
        new([])
    end
end

function add_vector!(glyph::VectorGlyph, x1::Float64, y1::Float64, x2::Float64, y2::Float64)
    push!(glyph.vectors, Point{Float64}(x1, y1))
    push!(glyph.vectors, Point{Float64}(x2, y2))
end

mutable struct VectorFont
    vectors::Array{VectorGlyph,1}
    glyphs::Dict{Char,Int64}

    horizontal_offset::Float64
    vertical_offset::Float64
    scale::Float64

    function VectorFont()
        new([],
        Dict{Char,Int64}(),
        1.2,
        1.2,
        3.0)
    end
end

function load_font!(font::VectorFont, file::String)
    try
        idx = 1
        lines = readlines(file)
        vectors = VectorGlyph()

        for line in lines
            if length(line) == 1
                # println("[", line, "]", idx)
                # Add character to glyph dictionary
                font.glyphs[line[1]] = idx
                idx += 1
                vectors = VectorGlyph()
                continue
            end

            if line == "||"
                push!(font.vectors, vectors)
                continue
            end

            # readlines until end of pixel marker: "||"
            if line ≠ "||"
                ele = split(line, " ")
                # println("ele: ", ele)
                # Add data to pixel array
                v1 = parse(Float64, ele[1])
                v2 = parse(Float64, ele[2])
                v3 = parse(Float64, ele[3])
                v4 = parse(Float64, ele[4])
                add_vector!(vectors, v1, v2, v3, v4)
            end
        end
    catch ex
        println("******************************************************")
        println(ex)
        println("Exception loading: ", file)
        println("******************************************************")
        return false
    end

    true

end

function get_glyph(font::VectorFont, char::Char) # returns a VectorGlyph
    glyph = font.glyphs[char]
    font.vectors[glyph]
end

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


text:
using BlockArrays
using ControlSystems
using LinearAlgebra
using PredictiveControl
using Test

# Setup the problem using a test case from Richter's paper
n = 4
m = 2
N = 5

A = [ 0.7 -0.1  0.0  0.0;
      0.2 -0.5  0.1  0.0;
      0.0  0.1  0.1  0.0;
      0.5  0.0  0.5  0.5 ]

B = [ 0.0  0.1;
      0.1  1.0;
      0.1  0.0;
      0.0  0.0 ]

C = Matrix{Float64}( I, 4, 4 )

D = fill( 0.0, (4, 2) )

sys = StateSpace( A, B, C, D, 0.01 )

Gₛ = [ 1.0  0.0;
      -1.0  0.0;
       0.0  1.0;
       0.0 -1.0 ]
gₛ = [ 1.0;
      -1.0;
       1.0;
      -1.0 ]

clqr = ConstrainedTimeInvariantLQR( sys, N, Matrix( 1.0*I(n) ), Matrix( 1.0*I(m) ), :Q )

H = PredictiveControl.hessian( clqr, PrimalFullCondensing )
G = kron( I(N), Gₛ )
g = reshape( kron( ones( Float64, N, 1 ), gₛ ), (:) )

@test isapprox( FGM.coldstartdelta( G, g, H = Matrix( H ) ), 19.026675, atol = 1e-4 )
@test isapprox( FGM.coldstartdelta( G, g, L = eigmax( Matrix( H ) ) ), 19.026675, atol = 1e-4 )


text:
mutable struct MWSingleLayerField3D{K}
    wavenumber::K
end

"""
    MWSingleLayerField3D(wavenumber = error())

Create the single layer near field operator, for use with `potential`.
"""
MWSingleLayerField3D(;wavenumber = error("Missing arg: `wavenumber`")) = MWSingleLayerField3D(wavenumber)

quaddata(op::MWSingleLayerField3D,rs,els) = quadpoints(rs,els,(2,))
quadrule(op::MWSingleLayerField3D,refspace,p,y,q,el,qdata) = qdata[1,q]

function kernelvals(op::MWSingleLayerField3D,y,p)

        k = op.wavenumber
        r = y - cartesian(p)
        R = norm(r)

        kr = k * R
        expn = exp(-im * k * R)
        green = expn / (4pi*R)
        gradgreen = -(im*k + 1/R) * green / R * r

        KernelValsMaxwell3D(im*k, r, R, green, gradgreen)
end

function integrand(op::MWSingleLayerField3D, krn, y, fp, p)

    j = fp[1]
    ρ = fp[2]

    #κ = krn.
    γ = krn.gamma
    G = krn.green
    ∇G = krn.gradgreen

    -γ*G*j + ∇G*ρ/γ
end


text:
module TestMacroTrap
using Debug
import Debug.Node, Debug.@instrument, Debug.Scope, Debug.BPNode

trap(::Node, ::Scope) = nothing
function trap(::BPNode, scope::Scope)
    @assert scope[:z] == 3
    scope[:z] = 4
end

macro m()
    esc(quote
        @bp
    end)
end

@instrument trap let
    z = 3
    @m
    @assert z == 4
end

end # module


text:
# Autogenerated wrapper script for MiniFB_jll for x86_64-linux-gnu
export libminifb

using Xorg_libX11_jll
using xkbcommon_jll
using Libglvnd_jll
JLLWrappers.@generate_wrapper_header("MiniFB")
JLLWrappers.@declare_library_product(libminifb, "libminifb.so")
function __init__()
    JLLWrappers.@generate_init_header(Xorg_libX11_jll, xkbcommon_jll, Libglvnd_jll)
    JLLWrappers.@init_library_product(
        libminifb,
        "lib/libminifb.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# Maximum adjacency visit / traversal

#################################################
#
#  Maximum adjacency visit
#
#################################################

mutable struct MaximumAdjacency <: AbstractGraphVisitAlgorithm
end

abstract type AbstractMASVisitor <: AbstractGraphVisitor end
# @static if VERSION > v"0.6-"
#   abstract type AbstractMASVisitor <: AbstractGraphVisitor end
# else
#   abstract AbstractMASVisitor <: AbstractGraphVisitor
#   import Base.Collections.PriorityQueue
# end

function maximum_adjacency_visit_impl!(
  graph::AbstractGraph{V,E},	                      # the graph
  pq::PriorityQueue{V,W},                           # priority queue
  visitor::AbstractMASVisitor,                      # the visitor
  colormap::Vector{Int}) where {V,E,W}              # traversal status

  while !isempty(pq)
    u = DataStructures.dequeue!(pq)
    # u = VERSION > v"0.6-" ? DataStructures.dequeue!(pq) : Collections.dequeue!(pq)
    discover_vertex!(visitor, u)
    for e in out_edges(u, graph)
      examine_edge!(visitor, e, 0)
      v = e.target

      if haskey(pq,v)
        pq[v] += edge_property(visitor.edge_weights, e, graph)
      end
    end
    close_vertex!(visitor, u)
  end

end

function traverse_graph(
  graph::AbstractGraph{V,E},
  alg::MaximumAdjacency,
  s::V,
  visitor::AbstractMASVisitor,
  colormap::Vector{Int},
  ::Type{W}) where {V,E,W}

  # @show Base.Order.Reverse
  # @show Base.Order.ReverseOrdering
  # PriorityQueue{K, V}(::Type{K}, ::Type{V}, o::Ordering) is deprecated, use PriorityQueue{K, V, typeof(o)}(o)
  # pq = DataStructures.PriorityQueue(V,W,Base.Order.Reverse)
  # pq = DataStructures.PriorityQueue(V,W,Base.Order.Reverse)
  # @show pq
  pq = DataStructures.PriorityQueue{V,W}(Base.Order.Reverse)
  # @show pq
  # pq = DataStructures.PriorityQueue{V,W,Base.Order.Reverse}()

  # Set number of visited neighbours for all vertices to 0
  for v in vertices(graph)
    pq[v] = zero(W)
  end

  @graph_requires graph incidence_list vertex_list
  @assert haskey(pq,s)
  @assert num_vertices(graph) >= 2

  #Give the starting vertex high priority
  pq[s] = one(W)

  #start traversing the graph
  maximum_adjacency_visit_impl!(graph, pq, visitor, colormap)
end

#################################################
#
#  Visitors
#
#################################################

#################################################
#
#  Minimum Cut Visitor
#
#################################################

mutable struct MinCutVisitor{G<:AbstractGraph,V,W} <: AbstractMASVisitor
  graph::G
  parities::Vector{Bool}
  colormap::Vector{Int}
  bestweight::W
  cutweight::W
  visited::Integer
  edge_weights::AbstractEdgePropertyInspector{W}
  vertices::Vector{V}
end

function MinCutVisitor(graph::AbstractGraph{V,E}, edge_weights::AbstractEdgePropertyInspector{W}) where {V,E,W}
  n = num_vertices(graph)
  parities = falses(n)
  MinCutVisitor{typeof(graph),V,W}(graph, parities, zeros(n), Inf, 0, 0, edge_weights, V[])
end

function discover_vertex!(vis::MinCutVisitor, v)
  vi = vertex_index(v,vis.graph)
  vis.parities[vi] = false
  vis.colormap[vi] = 1
  push!(vis.vertices,v)
  true
end

function examine_edge!(vis::MinCutVisitor, e, color::Int)
  vi = vertex_index(e.target,vis.graph)
  ew = edge_property(vis.edge_weights, e, vis.graph)

  # if the target of e is already marked then decrease cutweight
  # otherwise, increase it

  if vis.colormap[vi] != color # here color is 0
    vis.cutweight -= ew
  else
    vis.cutweight += ew
  end
end

function close_vertex!(vis::MinCutVisitor, v)
  vi = vertex_index(v,vis.graph)
  vis.colormap[vi] = 2
  vis.visited += 1

  if vis.cutweight < vis.bestweight && vis.visited < num_vertices(vis.graph)
    vis.bestweight = vis.cutweight
    for u in vertices(vis.graph)
      ui = vertex_index(u, vis.graph)
      vis.parities[ui] = ( vis.colormap[ui] == 2)
    end
  end
end

#################################################
#
#  MAS Visitor
#
#################################################

mutable struct MASVisitor{V,W} <: AbstractMASVisitor
  io::IO
  vertices::Vector{V}
  edge_weights::AbstractEdgePropertyInspector{W}
  log::Bool
end

function discover_vertex!(visitor::MASVisitor, v)
  push!(visitor.vertices,v)
  visitor.log ? println(visitor.io, "discover vertex: $v") : nothing
end

function examine_edge!(visitor::MASVisitor, e, color::Int)
  visitor.log ? println(visitor.io, " -- examine edge: $e") : nothing
end

function close_vertex!(visitor::MASVisitor, v)
  visitor.log ? println(visitor.io, "close vertex: $v") : nothing
end

#################################################
#
#  Minimum Cut
#
#################################################

function min_cut(
  graph::AbstractGraph{V,E},
  edge_weights::AbstractEdgePropertyInspector{W}) where {V,E,W}

  @graph_requires graph incidence_list vertex_list
  visitor = MinCutVisitor(graph, edge_weights)
  colormap = zeros(Int, num_vertices(graph))

  traverse_graph(graph, MaximumAdjacency(), first( vertices(graph) ), visitor, colormap, W)

  return( visitor.parities, visitor.bestweight)
end

function min_cut(
  graph::AbstractGraph{V,E},
  edge_weight_vec::Vector{W}) where {V,E,W}

  @graph_requires graph incidence_list vertex_list

  edge_weights = VectorEdgePropertyInspector(edge_weight_vec)
  visitor = MinCutVisitor(graph, edge_weights)
  colormap = zeros(Int, num_vertices(graph))

  traverse_graph(graph, MaximumAdjacency(), first( vertices(graph) ), visitor, colormap, W)

  return( visitor.parities, visitor.bestweight)
end

function min_cut(graph::AbstractGraph{V,E}) where {V,E}
  m = num_edges(graph)
  min_cut(graph,ones(m))
end

function maximum_adjacency_visit(graph::AbstractGraph{V,E}, edge_weights::AbstractEdgePropertyInspector{W}; log::Bool=false, io::IO=stdout) where {V,E,W}
  visitor = MASVisitor(io, V[],edge_weights,log)
  traverse_graph(graph, MaximumAdjacency(), first(vertices(graph)), visitor, zeros(Int, num_vertices(graph)), W)
  visitor.vertices
end

function maximum_adjacency_visit(graph::AbstractGraph{V,E}, edge_weight_vec::Vector{W}; log::Bool=false, io::IO=stdout) where {V,E,W}
  edge_weights = VectorEdgePropertyInspector(edge_weight_vec)
  visitor = MASVisitor(io, V[],edge_weights,log)
  traverse_graph(graph, MaximumAdjacency(), first(vertices(graph)), visitor, zeros(Int, num_vertices(graph)), W)
  visitor.vertices
end

function maximum_adjacency_visit(graph::AbstractGraph{V,E}; log::Bool=false, io::IO=stdout) where {V,E}
  m = num_edges(graph)
  maximum_adjacency_visit(graph,ones(m); log=log, io=io)
end


text:
mutable struct MMCTM
    K::Vector{Int}          # topics
    D::Int                  # documents
    N::Vector{Vector{Int}}  # observations per document modality
    M::Int                  # modalities
    V::Vector{Int}          # vocab items per modality

    μ::Vector{Float64}                      # doc-topic mean
    Σ::Matrix{Float64}                      # doc-topic covariance
    invΣ::Matrix{Float64}                   # inverse Σ
    props::Vector{Vector{Vector{Float64}}}  # doc-topic
    α::Vector{Float64}                      # topic-word hyperparameter
    ϕ::Vector{Vector{Vector{Float64}}}      # topic-word

    # variational parameters
    ζ::Vector{Vector{Float64}}              # doc-topic normalizer
    θ::Vector{Vector{Matrix{Float64}}}      # Z
    λ::Vector{Vector{Float64}}              # doc-topic mean
    ν::Vector{Vector{Float64}}              # doc-topic variance
    γ::Vector{Vector{Vector{Float64}}}      # topic-word
    Elnϕ::Vector{Vector{Vector{Float64}}}   # expectation of ln(ϕ)

    X::Vector{Vector{Matrix{Int}}}

    converged::Bool
    elbo::Float64
    ll::Vector{Float64}

    function MMCTM(k::Vector{Int}, α::Vector{Float64}, V::Vector{Int},
            X::Vector{Vector{Matrix{Int}}}; init=:random)
        model = new()

        model.K = copy(k)
        model.α = copy(α)
        model.X = X
        model.D = length(X)
        model.M = length(k)
        model.N = [[sum(X[d][m][:, 2]) for m in 1:model.M] for d in 1:model.D]

        model.V = copy(V)

        MK = sum(model.K)

        model.μ = zeros(MK)
        model.Σ = Matrix{Float64}(I, MK, MK)
        model.invΣ = Matrix{Float64}(I, MK, MK)
        model.props = [
            [Array{Float64}(undef, model.K[m]) for m in 1:model.M]
            for d in 1:model.D
        ]

        model.θ = [
            [
                fill(1.0 / model.K[m], model.K[m], size(model.X[d][m])[1])
                for m in 1:model.M
            ] for d in 1:model.D
        ]

        if init == :random
            model.γ = [
                [rand(1:100, model.V[m]) for k in 1:model.K[m]]
                for m in 1:model.M
            ]
        elseif init == :document
            model.γ = [
                [ones(model.V[m]) for k in 1:model.K[m]] for m in 1:model.M
            ]
            for m in 1:model.M
                seed_docs = sample(1:model.D, model.K[m], replace=false)
                for k in model.K[m]
                    d = seed_docs[k]
                    model.γ[m][k][X[d][m][:, 1]] .+= X[d][m][:, 2]
                end
            end
        else
            error("init must be either :random or :document")
        end
        model.Elnϕ = deepcopy(model.γ)
        update_Elnϕ!(model)
        model.ϕ = deepcopy(model.γ)

        model.λ = [zeros(MK) for d in 1:model.D]
        model.ν = [ones(MK) for d in 1:model.D]

        model.ζ = [Array{Float64}(undef, model.M) for d in 1:model.D]
        for d in 1:model.D update_ζ!(model, d) end

        model.converged = false

        return model
    end
end

function MMCTM(k::Vector{Int}, α::Vector{Float64},
        X::Vector{Vector{Matrix{Int}}})
    D = length(X)
    M = length(k)
    V = zeros(Int, M)
    for d in 1:D
        for m in 1:M
            if size(X[d][m], 1) > 0
                V[m] = max(V[m], maximum(X[d][m][:, 1]))
            end
        end
    end

    return MMCTM(k, α, V, X)
end

function calculate_sumθ(model::MMCTM, d::Int)
    return vcat(
        [
            vec(sum(model.θ[d][m] .* model.X[d][m][:, 2]', dims=2))
            for m in 1:model.M
        ]...
    )
end

function calculate_Ndivζ(model::MMCTM, d::Int)
    return vcat(
        [
            fill(model.N[d][m] / model.ζ[d][m], model.K[m]) for m in 1:model.M
        ]...
    )
end

function update_λ!(model::MMCTM, d::Int)
    opt = Opt(:LD_MMA, sum(model.K))
    xtol_rel!(opt, 1e-4)
    xtol_abs!(opt, 1e-4)

    Ndivζ = calculate_Ndivζ(model, d)
    sumθ = calculate_sumθ(model, d)

    max_objective!(
        opt,
        (λ, ∇λ) -> λ_objective(
            λ, ∇λ, model.ν[d], Ndivζ, sumθ, model.μ, model.invΣ
        )
    )
    (optobj, optλ, ret) = optimize(opt, model.λ[d])
    model.λ[d] .= optλ
end

function update_props!(model::MMCTM)
    for d in 1:model.D
        offset = 1
        for m in 1:model.M
            η = model.λ[d][offset:(offset + model.K[m] - 1)]
            model.props[d][m] .= exp.(η) ./ sum(exp.(η))
            offset += model.K[m]
        end
    end
end

function update_ν!(model::MMCTM, d::Int)
    opt = Opt(:LD_MMA, sum(model.K))
    lower_bounds!(opt, 1e-7)
    xtol_rel!(opt, 1e-4)
    xtol_abs!(opt, 1e-4)

    Ndivζ = calculate_Ndivζ(model, d)

    max_objective!(
        opt,
        (ν, ∇ν) -> ν_objective(ν, ∇ν, model.λ[d], Ndivζ, model.μ, model.invΣ)
    )
    (optobj, optν, ret) = optimize(opt, model.ν[d])
    model.ν[d] .= optν
end

function update_ζ!(model::MMCTM, d::Int)
    start = 1
    for m in 1:model.M
        stop = start + model.K[m] - 1
        model.ζ[d][m] = sum(
            exp.(model.λ[d][start:stop] .+ 0.5 * model.ν[d][start:stop])
        )
        start += model.K[m]
    end
end

function update_θ!(model::MMCTM, d::Int)
    offset = 0
    for m in 1:model.M
        for w in 1:size(model.X[d][m], 1)
            v = model.X[d][m][w, 1]

            for k in 1:model.K[m]
                model.θ[d][m][k, w] = exp(
                    model.λ[d][offset + k] + model.Elnϕ[m][k][v]
                )
            end
        end
        model.θ[d][m] ./= sum(model.θ[d][m], dims=1)
        offset += model.K[m]
    end
end

function update_μ!(model::MMCTM)
    model.μ .= mean(model.λ)
end

function update_Σ!(model::MMCTM)
    model.Σ .= sum(diagm.(0 .=> model.ν))
    for d in 1:model.D
        diff = model.λ[d] .- model.μ
        model.Σ .+= diff * diff'
    end
    model.Σ ./= model.D
    model.invΣ .= inv(model.Σ)
end

function update_Elnϕ!(model::MMCTM)
    for m in 1:model.M
        for k in 1:model.K[m]
            model.Elnϕ[m][k] .= (
                digamma.(model.γ[m][k]) .- digamma(sum(model.γ[m][k]))
            )
        end
    end
end

function update_γ!(model::MMCTM)
    for m in 1:model.M
        for k in 1:model.K[m]
            model.γ[m][k] .= model.α[m]
        end
    end
    for d in 1:model.D
        for m in 1:model.M
            Nθ = model.θ[d][m] .* model.X[d][m][:, 2]'
            for w in 1:size(model.X[d][m])[1]
                v = model.X[d][m][w, 1]
                for k in 1:model.K[m]
                    model.γ[m][k][v] += Nθ[k, w]
                end
            end
        end
    end
    update_Elnϕ!(model)
end

function update_ϕ!(model::MMCTM)
    for m in 1:model.M
        for k in 1:model.K[m]
            model.ϕ[m][k] .= model.γ[m][k] ./ sum(model.γ[m][k])
        end
    end
end

function update_α!(model::MMCTM)
    opt = Opt(:LD_MMA, 1)
    lower_bounds!(opt, 1e-7)
    xtol_rel!(opt, 1e-5)
    xtol_abs!(opt, 1e-5)

    for m in 1:model.M
        sum_Elnϕ = sum(sum(model.Elnϕ[m][k] for k in 1:model.K[m]))

        max_objective!(
            opt,
            (α, ∇α) -> α_objective(α, ∇α, sum_Elnϕ, model.K[m], model.V[m])
        )

        (optobj, optα, ret) = optimize(opt, model.α[m:m])
        model.α[m] = optα[1]
    end
end

function calculate_ElnPϕ(model::MMCTM)
    lnp = 0.0

    for m in 1:model.M
        for k in 1:model.K[m]
            lnp -= logmvbeta(fill(model.α[m], model.V[m]))
            for v in 1:model.V[m]
                lnp += (model.α[m] - 1) * model.Elnϕ[m][k][v]
            end
        end
    end

    return lnp
end

function calculate_ElnPη(model::MMCTM)
    lnp = 0.0

    for d in 1:model.D
        diff = model.λ[d] .- model.μ
        lnp += 0.5 * (
            logdet(model.invΣ) -
            sum(model.K) * log(2π) -
            tr(diagm(0 => model.ν[d]) * model.invΣ) -
            (diff' * model.invΣ * diff)[1]
        )
    end

    return lnp
end

function calculate_ElnPZ(model::MMCTM)
    lnp = 0.0

    for d in 1:model.D
        Eeη = exp.(model.λ[d] .+ 0.5model.ν[d])
        sumθ = calculate_sumθ(model, d)
        Ndivζ = calculate_Ndivζ(model, d)

        lnp += sum(model.λ[d] .* sumθ)
        lnp -= sum(Ndivζ .* Eeη) - sum(model.N[d])
        lnp -= sum(model.N[d] .* log.(model.ζ[d]))
    end

    return lnp
end

function calculate_ElnPX(model::MMCTM)
    lnp = 0.0

    for d in 1:model.D
        for m in 1:model.M
            for w in 1:size(model.X[d][m])[1]
                v = model.X[d][m][w, 1]
                for k in 1:model.K[m]
                    lnp += (
                        model.X[d][m][w, 2] * model.θ[d][m][k, w] *
                        model.Elnϕ[m][k][v]
                    )
                end
            end
        end
    end

    return lnp
end

function calculate_ElnQϕ(model::MMCTM)
    lnq = 0.0

    for m in 1:model.M
        for k in 1:model.K[m]
            lnq += -logmvbeta(model.γ[m][k])
            for v in 1:model.V[m]
                lnq += (model.γ[m][k][v] - 1) * model.Elnϕ[m][k][v]
            end
        end
    end
    return lnq
end

function calculate_ElnQη(model::MMCTM)
    lnq = 0.0
    for d in 1:model.D
        lnq += -0.5 * (sum(log.(model.ν[d])) + sum(model.K) * (log(2π) + 1))
    end
    return lnq
end

function calculate_ElnQZ(model::MMCTM)
    lnq = 0.0
    for d in 1:model.D
        for m in 1:model.M
            lnq += sum(
                model.X[d][m][:, 2]' .* log.(model.θ[d][m] .^ model.θ[d][m])
            )
        end
    end
    return lnq
end

function calculate_elbo(model::MMCTM)
    elbo = 0.0
    elbo += calculate_ElnPϕ(model)
    elbo += calculate_ElnPη(model)
    elbo += calculate_ElnPZ(model)
    elbo += calculate_ElnPX(model)
    elbo -= calculate_ElnQϕ(model)
    elbo -= calculate_ElnQη(model)
    elbo -= calculate_ElnQZ(model)
    return elbo
end

function calculate_docmodality_loglikelihood(X::Matrix{Int},
        props::Vector{Float64}, ϕ::Vector{Vector{Float64}})

    K = length(ϕ)

    ll = 0.0
    for w in 1:size(X, 1)
        v = X[w, 1]
        pw = 0.0
        for k in 1:K
            pw += props[k] * ϕ[k][v]
        end
        ll += X[w, 2] * log(pw)
    end

    return ll / sum(X[:, 2])
end

function calculate_modality_loglikelihood(X::Vector{Matrix{Int}},
        props::Vector{Vector{Float64}}, ϕ::Vector{Vector{Float64}})
    D = length(X)

    ll = 0.0
    N = 0
    for d in 1:D
        doc_N = sum(X[d][:, 2])
        if doc_N > 0
            doc_ll = calculate_docmodality_loglikelihood(X[d], props[d], ϕ)
            ll += doc_ll * doc_N
            N += doc_N
        end
    end

    return ll / N
end

function calculate_loglikelihoods(X::Vector{Vector{Matrix{Int}}},
                                  props::Vector{Vector{Vector{Float64}}},
                                  ϕ::Vector{Vector{Vector{Float64}}})
    D = length(X)
    M = length(ϕ)
    K = [length(ϕ[m]) for m in 1:M]

    ll = Array{Float64}(undef, M)

    offset = 1
    for m in 1:M
        Xm = [X[d][m] for d in 1:D]
        propsm = [props[d][m] for d in 1:D]

        ll[m] = calculate_modality_loglikelihood(Xm, propsm, ϕ[m])

        offset += K[m]
    end

    return ll
end

function calculate_loglikelihoods(X::Vector{Vector{Matrix{Int}}}, model::MMCTM)
    return calculate_loglikelihoods(X, model.props, model.ϕ)
end

function calculate_loglikelihoods(model::MMCTM)
    return calculate_loglikelihoods(model.X, model.props, model.ϕ)
end

function fitdoc!(model::MMCTM, d::Int)
    update_ζ!(model, d)
    update_θ!(model, d)
    update_ν!(model, d)
    update_λ!(model, d)
end

function fit!(model::MMCTM; maxiter=100, tol=1e-4, verbose=true, autoα=false,
        updateΣ=true)

    ll = Vector{Float64}[]

    for iter in 1:maxiter
        for d in 1:model.D
            fitdoc!(model, d)
        end

        update_μ!(model)
        if updateΣ
            update_Σ!(model)
        end
        update_γ!(model)
        if autoα
            update_α!(model)
        end

        update_props!(model)
        update_ϕ!(model)

        push!(ll, calculate_loglikelihoods(model))

        if verbose
            println("$iter\tLog-likelihoods: ", join(ll[end], ", "))
        end

        if length(ll) > 10 && check_convergence(ll, tol=tol)
            model.converged = true
            break
        end
    end
    model.elbo = calculate_elbo(model)
    model.ll = ll[end]

    return ll
end

function unsmoothed_update_θ!(model::MMCTM, d::Int)
    offset = 0
    for m in 1:model.M
        for w in 1:size(model.X[d][m], 1)
            v = model.X[d][m][w, 1]

            for k in 1:model.K[m]
                model.θ[d][m][k, w] = exp(model.λ[d][offset + k]) * model.ϕ[m][k][v]
            end
        end
        model.θ[d][m] ./= sum(model.θ[d][m], dims=1)
        offset += model.K[m]
    end
end

function transform(model::MMCTM, X::Vector{Vector{Matrix{Int}}};
                   maxiter=1000, tol=1e4, verbose=false)

    newmodel = MMCTM(model.K, model.α, X)
    newmodel.ϕ = deepcopy(model.ϕ)

    ll = Vector{Float64}[]
    for iter in 1:maxiter
        for d in 1:newmodel.D
            update_ζ!(newmodel, d)
            unsmoothed_update_θ!(newmodel, d)
            update_ν!(newmodel, d)
            update_λ!(newmodel, d)
        end

        update_μ!(newmodel)
        update_Σ!(newmodel)

        update_props!(newmodel)

        push!(ll, calculate_loglikelihoods(newmodel))

        if verbose
            println("$iter\tLog-likelihoods: ", join(ll[end], ", "))
        end

        if length(ll) > 10 && check_convergence(ll, tol=tol)
            newmodel.converged = true
            break
        end
    end

    return newmodel
end

function fit_heldout(Xheldout::Vector{Vector{Matrix{Int}}}, model::MMCTM;
        maxiter=100, verbose=false)

    heldout_model = MMCTM(model.K, model.α, Xheldout)
    heldout_model.μ .= model.μ
    heldout_model.Σ .= model.Σ
    heldout_model.invΣ .= model.invΣ
    heldout_model.γ = deepcopy(model.γ)
    heldout_model.Elnϕ = deepcopy(model.Elnϕ)
    heldout_model.ϕ = deepcopy(model.ϕ)

    ll = Vector{Float64}[]
    for iter in 1:maxiter
        for d in 1:heldout_model.D
            fitdoc!(heldout_model, d)
        end

        update_props!(heldout_model)

        push!(ll, calculate_loglikelihoods(heldout_model))

        if verbose
            println("$iter\tLog-likelihoods: ", join(ll[end], ", "))
        end

        if length(ll) > 10 && check_convergence(ll)
            heldout_model.converged = true
            break
        end
    end

    return heldout_model
end

function predict_modality_η(Xobs::Vector{Vector{Matrix{Int}}}, m::Int,
        model::MMCTM; maxiter=100, verbose=false)
    obsM = setdiff(1:model.M, m)

    moffset = sum(model.K[1:(m - 1)])
    unobsMK = (moffset + 1):(moffset + model.K[m])
    obsMK = setdiff(1:sum(model.K), unobsMK)

    obsmodel = MMCTM(model.K[obsM], model.α[obsM], Xobs)
    obsmodel.μ .= model.μ[obsMK]
    obsmodel.Σ .= model.Σ[obsMK, obsMK]
    obsmodel.invΣ .= model.invΣ[obsMK, obsMK]
    obsmodel.γ = deepcopy(model.γ[obsM])
    obsmodel.Elnϕ = deepcopy(model.Elnϕ[obsM])

    ll = Vector{Float64}[]
    for iter in 1:maxiter
        for d in 1:obsmodel.D
            fitdoc!(obsmodel, d)
        end

        push!(ll, calculate_loglikelihoods(Xobs, obsmodel))

        if verbose
            println("$iter\tLog-likelihoods: ", join(ll[end], ", "))
        end

        if length(ll) > 10 && check_convergence(ll)
            obsmodel.converged = true
            break
        end
    end

    if !obsmodel.converged
        warn("model not converged.")
    end

    η = [
        (
            model.μ[unobsMK] .+ model.Σ[unobsMK, obsMK] *
            model.invΣ[obsMK, obsMK] * (obsmodel.λ[d] .- model.μ[obsMK])
        )
        for d in 1:obsmodel.D
    ]

    return η
end