text:
module TestMechanics

println("TestMechanics: Tests how 3D mechanics could be combined with ModiaLang.")

using ModiaLang
using Unitful
@usingModiaPlot
using DifferentialEquations
using RuntimeGeneratedFunctions
RuntimeGeneratedFunctions.init(@__MODULE__)

# ModiaLang models
include("$(ModiaLang.path)/models/Blocks.jl")
include("$(ModiaLang.path)/models/Electric.jl")
include("$(ModiaLang.path)/models/Rotational.jl")

# Generic ModiaLang3D model stubs (always available)
RevoluteStub = Model(
    flange = Flange,

    equations = :[phi  = flange.phi
                  w    = der(phi)
                 der_w = der(w)
                 tau   = flange.tau]
)

# Generic Modia3D definitions (always available)
mutable struct SimulationModel3D
    # Internal memory of 3D model
    initial::Bool

    SimulationModel3D() = new(true)
end

# Function ..._f1 for 3D mechanics
# function Pendulum_f1(_m::SimulationModel3D, phi::Float64, w::Float64, tau::Float64,
#                       m::Float64, L::Float64, g::Float64)
function Pendulum_f1(_m::TestMechanics.SimulationModel3D, phi, w, tau, m, L, g)
#function Pendulum_f1(phi, w, tau, m, L, g)

#=
    if _m.initial
        println("... SimplePendulum_f1 called during initialization")
        _m.initial = false
    end
=#
    der_w = (tau - m*g*L*cos(phi))/m*L*L
    return der_w
end

# ModiaLang model for the system
Pendulum = Model(
    model3D = SimulationModel3D(),

    # Parameters
    m = 1.0u"kg",
	L = 0.5u"m",
	g = 9.81u"m/s^2",

    # 3D model stubs
    rev = RevoluteStub | Map(init = Map(phi=1.5u"rad", w=2u"rad/s")),

    # Standard ModiaLang models
    damper  = Damper | Map(d=0.4u"N*m*s/rad"),
	support = Fixed,

    connect = :[
        (rev.flange     , damper.flange_b),
        (damper.flange_a, support.flange)],

    equations = :(rev.der_w = Pendulum_f1(model3D, rev.phi, rev.w, rev.tau, m, L, g))
#    equations = :(rev.der_w = Pendulum_f1(rev.phi, rev.w, rev.tau, m, L, g))
)

pendulum = @instantiateModel(Pendulum, logCode=true, logExecution=false, unitless=true)

simulate!(pendulum, Tsit5(), stopTime=20.0, log=true)
plot(pendulum, ["rev.phi", "rev.w"])

end

max_stars_repo_name: ModiaSim/ModiaLang

max_issues_repo_name: ModiaSim/ModiaLang

max_forks_repo_name: ModiaSim/ModiaLang


text:
using UnicodeGraphics, OffsetArrays, Test

pac = [0 0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0;
       0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0;
       0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0;
       0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0;
       0 0 1 1 1 1 1 1 1 0 0 1 1 1 1 1 1 1 1 0;
       0 1 1 1 1 1 1 1 1 0 0 1 1 1 1 1 1 1 1 0;
       0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0;
       1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0;
       1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0;
       1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0;
       1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 1;
       1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0;
       1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0;
       0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0;
       0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0;
       0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0;
       0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0;
       0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0;
       0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0;
       0 0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0]

block_pac = """
     ▄▄██████▄▄
  ▄██████████████▄
 ▄███████  ████████
▄██████████████▀▀
███████████▀▀
███████████▄▄      ▀
▀██████████████▄▄
 ▀█████████████████
  ▀██████████████▀
     ▀▀██████▀▀     \0"""

test_pac = blockize(pac)
@test test_pac == block_pac
println(test_pac, "\n\n", block_pac, "\n\n")

brail_pac =
"""
⠀⣠⣴⣾⣿⣿⣷⣦⣄
⣰⣿⣿⣿⣧⣼⣿⡿⠟⠃
⣿⣿⣿⣿⣿⣏⡁⠀⠀⠠
⠹⣿⣿⣿⣿⣿⣿⣷⣦⡄
⠀⠙⠻⢿⣿⣿⡿⠟⠋⠀\0"""

test_pac = brailize(pac)
@test test_pac == brail_pac
println(test_pac, "\n\n", brail_pac, "\n\n")

ghost = [1.0 0.0 1.0 0.0 1.0 0.0 1.0 0.0 1.0 0.0 1.0 0.0 1.0 0.0 1.0 0.0 1.0 0.0 1.0;
         0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
         1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 1.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0;
         0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 0.0 0.0 0.0 0.0 0.0;
         1.0 0.0 0.0 0.0 0.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 0.0 0.0 0.0 1.0;
         0.0 0.0 0.0 0.0 1.0 0.0 0.0 1.0 1.0 1.0 1.0 0.0 0.0 1.0 1.0 1.0 0.0 0.0 0.0;
         1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 0.0 0.0 0.0 0.0 1.0 1.0 0.0 0.0 1.0;
         0.0 0.0 0.0 0.0 1.0 1.0 0.0 0.0 1.0 1.0 1.0 1.0 0.0 0.0 1.0 1.0 0.0 0.0 0.0;
         1.0 0.0 0.0 1.0 1.0 1.0 0.0 0.0 1.0 1.0 1.0 1.0 0.0 0.0 1.0 1.0 1.0 0.0 1.0;
         0.0 0.0 0.0 1.0 1.0 0.0 0.0 1.0 1.0 1.0 1.0 0.0 0.0 1.0 1.0 1.0 1.0 0.0 0.0;
         1.0 0.0 0.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 0.0 1.0;
         0.0 0.0 0.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 0.0 0.0;
         1.0 0.0 0.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 0.0 1.0;
         0.0 0.0 0.0 1.0 1.0 0.0 1.0 1.0 1.0 0.0 0.0 1.0 1.0 1.0 0.0 1.0 1.0 0.0 0.0;
         1.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0 1.0 0.0 0.0 1.0 1.0 0.0 0.0 0.0 1.0 0.0 1.0;
         0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
         1.0 0.0 1.0 0.0 1.0 0.0 1.0 0.0 1.0 0.0 1.0 0.0 1.0 0.0 1.0 0.0 1.0 0.0 1.0]

block_ghost =
"""
   ▄▄████▄▄
 ▄▀▀████▀▀██▄
 ▄▄  ██▄▄  ██
██▀ ▄███▀ ▄███
██████████████
██▀███▀▀███▀██
▀   ▀▀  ▀▀   ▀\0"""

test_ghost = blockize(OffsetArray(ghost[3:15, 4:17], 3:15, 4:17), 0.5)
@test test_ghost == block_ghost
println(test_ghost, "\n\n", block_ghost, "\n\n")

braile_ghost =
"""
⠀⣠⣴⣶⣦⣄
⣨⡄⢹⣯⡄⢹⣇
⣿⣶⣿⣿⣶⣿⣿
⠋⠈⠛⠀⠛⠁⠙\0"""

test_ghost = brailize(OffsetArray(ghost[2:15, 4:17], 2:15, 4:17), 0.5)
@test test_ghost == braile_ghost
println(test_ghost, "\n\n", braile_ghost)

ext: jl


text:
"""
    boostParticlesPosition!(beam::Particles, fp::SVector{3,T}, factor::T) where {T}

Boost the position of each particle contained in `beam` to the
inertial frame moving with the momentum `fp` by `factor`.
"""
function boostParticlesPosition!(beam::Particles, fp::SVector{3,T}, factor::T) where {T}
    pos = beam.positions
    @inbounds for k in 1:beam.npar
        x = pos[k]
        pos[k] = boostPosition(x, fp, factor)
    end
end

function boostPosition(x::SVector{3, T}, fp::SVector{3,T}, factor::T) where {T}
    fp_hat = fp / norm(fp)
    return x + (factor - 1.0) * dot(x, fp_hat) * fp_hat
end

"""
    transformParticlesMomentum!(beam::Particles, fp::SVector{3,T}) where {T}

Transform the momentum of each particle contained in `beam`
to the inertial frame moving with the momentum `fp`.
"""
function transformParticlesMomentum!(beam::Particles, fp::SVector{3,T}) where {T}
    fg = sqrt(1.0 + dot(fp, fp))
    mom = beam.momentums
    @inbounds for k in 1:beam.npar
       p = mom[k]
       mom[k] = transformMomentum(p, fg, fp)
    end
end
function transformMomentum(p::SVector{3,T}, fg::T, fp::SVector{3,T}) where {T}
    g = sqrt(1.0 + dot(p,p))
    return  p + 1.0/(fg + 1.0) * dot(fp, p)*fp - g*fp
end

"""
    transformParticlesField!(beam::Particles, fp::SVector{3,T}) where {T}

Transform the EM field experienced by each particle contained in `beam`
to the inertial frame moving with the momentum `fp`.
"""
function transformParticlesField!(beam::Particles, fp::SVector{3, T}) where {T}
    fg = sqrt(1.0 + dot(fp, fp))
    efields = beam.self_efields
    bfields = beam.self_bfields
    @inbounds for i in 1:beam.npar
      efield_f, bfield_f = transformEMField(efields[i], bfields[i], fg, fp)
      efields[i] = efield_f
      bfields[i] = bfield_f
    end
  end

function transformEMField(efield::SVector{3, T}, bfield::SVector{3, T}, fg::T, fp::SVector{3, T}) where {T}
  efield_f = fg*efield + cross(fp, bfield) - 1.0/(fg + 1.0)*dot(fp, efield)*fp
  bfield_f = fg*bfield - cross(fp, efield) - 1.0/(fg + 1.0)*dot(fp, bfield)*fp
  return efield_f, bfield_f
end


text:
using FFTW
using DSP
using LinearAlgebra

"make fourier filter size_padded is power of 2"
function make_filter(sz, filter_type="ramlak")
    # See Kak Chapter3 (61)

    tau = cat(Array(1:2:Int(sz/2-1)), Array(Int(sz/2-1):-2:0), dims=1);

    # note that the filter is already shifted
    spatial_filter = zeros(sz)
    spatial_filter[1] = 0.25
    spatial_filter[2:2:end] .= -1.0 ./ (pi * tau).^2

    fourier_filter = 2 * real(fft(spatial_filter))

    if filter_type == "shepplogan"
        omega = pi * FFTW.fftfreq(sz)[2:end]
        fourier_filter[2:end] .*= sin.(omega) ./ omega
    elseif filter_type == "hamming"
        fourier_filter .*= fftshift(hamming(sz))
    elseif filter_type == "hann"
        fourier_filter .*= fftshift(hanning(sz))
    elseif filter_type == "kaiser"
        fourier_filter .*= fftshift(kaiser(sz, 3.0))
    elseif filter_type == "consine"
        fourier_filter .*= fftshift(consine(sz))
    end

    return fourier_filter
end

@doc raw"""
    function filter_proj(p, filter_type)

Filter projections p by a given filter type (p: [detcount x nangles])

supported filter_type:
ramlak, shepplogan, hamming, hann
"""
function filter_proj(p::Array{T, 2}, filter_type="ramlak") where {T<:AbstractFloat}
    # See Kak ch3 (61)

    nangles, detcount = size(p)

    # compute next power of 2
    size_padded = Int(2^ceil(log(detcount)/log(2)))
    fourier_filter = make_filter(size_padded, filter_type)

    dpad = size_padded - detcount

    # pad data
    p_padded = zeros(nangles, size_padded)
    for i=1:nangles
        p_padded[i, 1:detcount] = p[i, :]
    end

    # filter projection in Fourier domain
    proj = fft(p_padded, 2) .* reshape(fourier_filter, 1, length(fourier_filter))
    pimg_ifft = real(ifft(proj, 2)[:, 1:detcount])
    return pimg_ifft / 2.0; # for matching with the original
end

"""
    function filter_proj(p::Array{T, 3}, filter_type="ramlak") where {T<:AbstractFloat}

Filter projections slice by slice where p: [nangles x H x W]
"""
function filter_proj(p::Array{T, 3}, filter_type="ramlak") where {T<:AbstractFloat}
    # See Kak ch3 (61)

    nangles, nz, detcount = size(p)
    out = similar(p)

    # compute next power of 2
    size_padded = Int(2^ceil(log(detcount)/log(2)))
    fourier_filter = make_filter(size_padded, filter_type)
    fourier_filter = reshape(fourier_filter, 1, length(fourier_filter))

    dpad = size_padded - detcount

    p_padded = zeros(nangles, size_padded)
    for z=1:nz
        # pad data
        fill!(p_padded, 0)
        for i=1:nangles
            p_padded[i, 1:detcount] = p[i, z, :]
        end

        # filter projection in Fourier domain
        proj = fft(p_padded, 2) .* fourier_filter

        pimg_fft = real(ifft(proj, 2)[:, 1:detcount])
        out[:, z, :] .= pimg_fft / 2.0
    end
    return out
end

"""
    bp_slices(p_, A, H, W, scaling=true)

Perform backprojection slice by slice
# Args
- p_ [nangles x nslices x detcount]
- A : forward operator
"""
function bp_slices(p_, A, H, W, scaling=true)
    nangles, nslice, detcount = size(p_)
    At = A'

    nslice = size(p_, 2)
    img = zeros(H*W, nslice)
    # temp = zeros(H*W)

    p_axWxH = permutedims(p_, [1, 3, 2])
    p = reshape(p_axWxH, :, nslice)

    Threads.@threads for slice=1:nslice
        img_slice = view(img, :, slice)
        pslice = view(p, :, slice)
        mul!(img_slice, At, pslice)
    end

    # scale for FBP
    if scaling
        img .*= (pi / nangles)
    end
    img = reshape(img, H, W, nslice)
    return img
end

# TODO
"""
    filter_proj_hilbert

Filter projections based on Hilbert transform.

(p: [detcount x nangles])
"""
# $q(s, \theta)=\frac{\mathrm{d} p(s, \theta)}{\mathrm{d} s} * \frac{-1}{2 \pi^{2} s}$
function filter_proj_hilbert(p)
    qimg = zeros(size(p))
    for i in 1:size(qimg, 1)
        qimg[i, 2:(end-1)] = pimg[i, 3:end] - pimg[i, 1:(end-2)]
    #     qimg[end,i] = 0
        qimg[i,:] *= -1.0 / (2*pi)
        qimg[i,:] = real(hilbert(qimg[i,:]))
    end
    return qimg
end

max_stars_repo_name: JuliaTomo/DynamicTomographyWithFlow

max_issues_repo_name: JuliaTomo/DynamicTomographyWithFlow

max_forks_repo_name: JuliaTomo/DynamicTomographyWithFlow


text:
using Survival: KaplanMeier
using Recombinase: Analysis
function _survival(times, events; axis = discrete_axis(times, npoints = npoints), kwargs...)
    km = fit(KaplanMeier, times, events)
    surv = zeros(length(axis))
    for (i, ax) in enumerate(axis)
        if ax < km.times[1]
            surv[i] = 1
        elseif ax > km.times[end]
            surv[i] = km.survival[end]
        else
            surv[i] = km.survival[searchsortedfirst(km.times, ax)]
        end
    end
    return zip(axis, surv)
end

const survival = Analysis((; discrete = _survival))


text:
using Ramnet.Encoders
using Ramnet.Models.AltDiscriminators:RegressionDiscriminator

using StaticArrays
using Random

using ..Approximators:DiscreteActionPolicy
using ..Buffers: MultiStepDynamicBuffer, add!

mutable struct DiscreteFunctionalAC{OS,AS,T <: Real,O <: StaticArray{Tuple{OS},T,1},E <: AbstractEncoder{T}} <: AbstractAgent{O,Int}
    steps::Int
    γ::Float64
    actor::DiscreteActionPolicy{OS,AS,T,O,E}
    critic::RegressionDiscriminator{1}
    buffer::MultiStepDynamicBuffer{O,Int}
    observation::Union{Nothing,O}
    done::Bool
    action::Union{Nothing,Int}
    rng::MersenneTwister # Is this even used?
end

function DiscreteFunctionalAC{AS}(::Type{O}, steps, n, start_learning_rate, end_learning_rate, learning_rate_decay, epochs, discount, forgetting_factor, encoder::E; seed::Union{Nothing,UInt}=nothing) where {OS,AS,T <: Real,O <: StaticArray{Tuple{OS},T,1},E <: AbstractEncoder{T}}
    !isnothing(seed) && seed < 0 && throw(DomainError(seed, "Seed must be non-negative"))
    rng = isnothing(seed) ? MersenneTwister() : MersenneTwister(seed)

    n < 1 && throw(DomainError(n, "Tuple size must be at least 1"))

    actor = DiscreteActionPolicy{AS}(O, n, encoder; start_learning_rate, end_learning_rate, learning_rate_decay, epochs, partitioner=:uniform_random, seed)
    critic = RegressionDiscriminator{1}(OS, 32, encoder; seed, γ=forgetting_factor)

    buffer = MultiStepDynamicBuffer{O,Int}()

    DiscreteFunctionalAC{OS,AS,T,O,E}(steps, discount, actor, critic, buffer, nothing, false, nothing, rng)
end

function DiscreteFunctionalAC(
  env, encoder::E; steps, tuple_size, start_learning_rate, end_learning_rate, learning_rate_decay, epochs=1, discount=1.0, forgetting_factor=1.0, seed=nothing) where {T <: Real,E <: AbstractEncoder{T}}
  DiscreteFunctionalAC{length(action_set(env))}(observation_type(env), steps, tuple_size, start_learning_rate, end_learning_rate, learning_rate_decay, epochs, discount, forgetting_factor, encoder; seed)
end

function observe!(agent::DiscreteFunctionalAC{OS,AS,T,O,E}, observation::O) where {OS,AS,T <: Real,O <: StaticArray{Tuple{OS},T,1},E <: AbstractEncoder{T}}
  agent.observation = observation
  agent.done = false
  agent.action = _select_action(agent, observation)

  nothing
end

# TODO: All observe! methods are the same for all agents. Generalize.
# TODO: This method does not need to take in the action
function observe!(agent::DiscreteFunctionalAC{OS,AS,T,O,E}, action::Int, reward::Float64, observation::O, done::Bool) where {OS,AS,T <: Real,O <: StaticArray{Tuple{OS},T,1},E <: AbstractEncoder{T}}
    add!(agent.buffer, agent.observation, agent.action, reward)

    if !done
        agent.action = _select_action(agent, observation)
    end

    agent.observation = observation
    agent.done = done
    nothing
end

function _select_action(agent::DiscreteFunctionalAC{OS,AS,T,O,E}, observation::O) where {OS,AS,T <: Real,O <: StaticArray{Tuple{OS},T,1},E <: AbstractEncoder{T}}
  return select_action(agent.actor, observation)
end

function select_action!(agent::DiscreteFunctionalAC{OS,AS,T,O,E}, observation::O) where {OS,AS,T,O,E}
  if !isnothing(agent.action)
      return agent.action
  end

  return _select_action(agent, observation)
end

function update!(agent::DiscreteFunctionalAC)
  if agent.done
      G = 0.0
      for transition in Iterators.reverse(agent.buffer)
          G = transition.reward + agent.γ * G

          δ = G - (predict(agent.critic, transition.observation) |> first)
          update!(agent.actor, transition.observation, transition.action, δ)
          train!(agent.critic, transition.observation, G)
      end

      reset!(agent.buffer)
  elseif length(agent.buffer) == agent.steps
      G = predict(agent.critic, agent.observation) |> first

      for transition in Iterators.reverse(agent.buffer)
          G = transition.reward + agent.γ * G
      end

      t = popfirst!(agent.buffer)

      δ = G - (predict(agent.critic, t.observation) |> first)
      update!(agent.actor, t.observation, t.action, δ)
      train!(agent.critic, t.observation, G)
  end

  nothing
end

function reset!(agent::DiscreteFunctionalAC; seed::Union{Nothing,UInt}=nothing)
  if !isnothing(seed)
      if seed ≥ 0
          Random.seed!(agent.rng, seed)
      else
          throw(DomainError(seed, "Seed must be non-negative"))
      end
  end

  reset!(agent.actor; seed)
  Ramnet.reset!(agent.critic)
  reset!(agent.buffer)

  nothing
end


text:
module TestTwoInertiasAndIdealGear

using Modia
@usingModiaPlot
using Test

TwoInertiasAndIdealGearTooManyInits = Model(
    J1 = 0.0025,
    J2 = 170,
    r  = 105,
    tau_max = 1,
    phi1 = Var(init = 0.0),
    w1   = Var(init = 1.0),
    phi2 = Var(init = 0.5),
    w2   = Var(init = 0.0),
    equations = :[
        tau = if time < 1u"s"; tau_max elseif time < 2u"s"; 0 elseif time < 3u"s"; -tau_max else 0 end,

        # inertia1
        w1 = der(phi1),
        J1*der(w1) = tau - tau1,

        # gear
        phi1   = r*phi2,
        r*tau1 = tau2,

        # inertia2]
        w2 = der(phi2),
        J2*der(w2) = tau2
    ]
)

TwoInertiasAndIdealGear = TwoInertiasAndIdealGearTooManyInits | Map(phi1 = Var(init=nothing), w1=Var(init=nothing))

twoInertiasAndIdealGearTooManyInits = @instantiateModel(TwoInertiasAndIdealGearTooManyInits, unitless=true)
twoInertiasAndIdealGear             = @instantiateModel(TwoInertiasAndIdealGear, unitless=true)

println("Next simulate! should result in an error:\n")
simulate!(twoInertiasAndIdealGearTooManyInits, Tsit5(), stopTime = 4.0, log=true)

simulate!(twoInertiasAndIdealGear, Tsit5(), stopTime = 4.0, log=true,
          logParameters=true, logStates=true,
          requiredFinalStates=[1.5628074713622309, -6.878080753044174e-5])

plot(twoInertiasAndIdealGear, ["phi2", "w2"])

simulate!(twoInertiasAndIdealGear, stopTime = 4.0,
          useRecursiveFactorizationUptoSize = 500,
          log=true, logParameters=true, logStates=true,
          requiredFinalStates=missing)
# Linearize
println("\n... Linearize at stopTime = 0 and 4")
(A_0, x_0) = linearize!(twoInertiasAndIdealGear, stopTime=0, analytic = true)
(A_4, x_4) = linearize!(twoInertiasAndIdealGear, stopTime=4, analytic = true)
(A_4_numeric, x_4_numeric) = linearize!(twoInertiasAndIdealGear, stopTime=4, analytic=false)

xNames = get_xNames(twoInertiasAndIdealGear)
@show xNames
@show A_0, x_0
@show A_4, x_4
@show A_4_numeric, x_4_numeric
@test isapprox(A_0,[0.0 1.0; 0.0 0.0])
@test isapprox(A_0, A_4)

end


text:
using Plots, Random, Statistics

N = 1000000
σ_η = 0.01  # variance
σ_ε = 0.05  # variance

# IID process
y = sqrt(σ_ε).*randn(N)
# Check data generation
abs(var(y) - σ_ε) < 0.001

# Plot discretizaiton vs simulated data
tauc = tauchen(19,0.,sqrt(σ_ε))
histogram(y, bins = range(tauc.state_values[1].-step(tauc.state_values)/2, step=step(tauc.state_values),
                    length=length(tauc.state_values)+1), normalize = :probability)
scatter!(tauc.state_values,tauc.p[1,:])

rouw = rouwenhorst(19,0.,sqrt(σ_ε))
histogram!(y, bins = range(rouw.state_values[1].-step(rouw.state_values)/2, step=step(rouw.state_values),
                    length=length(rouw.state_values)+1), normalize = :probability, alpha = 0.5)
compare = scatter!(rouw.state_values,rouw.p[1,:],alpha=0.5)
plot(compare)

# Sum of IID processes
y = sqrt(σ_ε).*randn(N) + sqrt(σ_η).*randn(N)
# Check data generation
abs(var(y) - (σ_ε + σ_η)) < 0.001

# Plot discretizaiton vs simulated data
tauc_ε = tauchen(19,0.,sqrt(σ_ε))
tauc_η = tauchen(19,0.,sqrt(σ_η))
plot(tauc_ε.state_values, tauc_η.state_values)
histogram(y, bins = range(tauc.state_values[1].-step(tauc.state_values)/2, step=step(tauc.state_values),
                    length=length(tauc.state_values)+1), normalize = :probability)
scatter!(tauc.state_values,tauc.p[1,:])

rouw = rouwenhorst(19,0.,sqrt(σ_ε))
histogram!(y, bins = range(rouw.state_values[1].-step(rouw.state_values)/2, step=step(rouw.state_values),
                    length=length(rouw.state_values)+1), normalize = :probability, alpha = 0.5)
compare = scatter!(rouw.state_values,rouw.p[1,:],alpha=0.5)
plot(compare)

# With persistent process
κ = zeros(Float64,N,T_ret-1) # *** Get from data
σ_η = 0.01
σ_z0 = 0.15
σ_ε = 0.05

function net_labor_income(N, T_ret, σ_ε, σ_η, σ_z0, κ)
    Y = zeros(Float64, N, T_ret-1)
    z_t1 = sqrt(σ_z0) .* randn(N)

    for t in 1:(T_ret-1)
        z_t = z_t1 + sqrt(σ_η) .* randn(N)
        Y[:,t] = exp.( κ[t] .+ z_t + sqrt(σ_ε) .* randn(N) )
        z_t1 = copy(z_t)
    end
    return Y
end

N = 100000
T_ret = 2
Y_l = net_labor_income(N, T_ret, σ_ε, σ_η, 0., κ)
histogram(Y_l, alpha = 0.5)
tauc_ε = tauchen(19,0.,σ_ε)
tauc_η1 = tauchen(39,0.,var(Y_l[:,1]))

# unconditional variance of z_t
T_ret = 40
ρ = 1.
σ

plot(σ)

using QuantEcon
rouwenhorst(19,0.,sqrt(0.2)).p == rouwenhorst(19,0.,sqrt(0.1)).p

import QuantEcon.rouwenhorst
function rouwenhorst(N::Integer, ρ::Real, σ::Real, μ::Real=0.0; stationary::Bool = true)
    if stationary
        σ_y = σ / sqrt(1-ρ^2)
        m = μ / (1 - ρ)
    else
        σ_y = σ
        m = 0.0  #μ not used if non-stationary
    end

    p  = (1+ρ)/2
    Θ = [p 1-p; 1-p p]
    ψ = sqrt(N-1) * σ_y
    state_values, p = _rouwenhorst(p, p, m, ψ, N)
    MarkovChain(p, state_values)
end

function _rouwenhorst(p::Real, q::Real, m::Real, Δ::Real, n::Integer)
    if n == 2
        return [m-Δ, m+Δ],  [p 1-p; 1-q q]
    else
        _, θ_nm1 = _rouwenhorst(p, q, m, Δ, n-1)
        θN = p    *[θ_nm1 zeros(n-1, 1); zeros(1, n)] +
             (1-p)*[zeros(n-1, 1) θ_nm1; zeros(1, n)] +
             q    *[zeros(1, n); zeros(n-1, 1) θ_nm1] +
             (1-q)*[zeros(1, n); θ_nm1 zeros(n-1, 1)]

        θN[2:end-1, :] ./= 2

        return range(m-Δ, stop=m+Δ, length=n), θN
    end
end

function rouwenhorst(N::Integer, ρ::Real, σ::Real, σ_0::Real, T::Int)
    σs = zeros(eltype(σ_0), T, 1)
    σs[1] = σ_0
    [ σs[t] = ρ^2*σs[t-1] + σ  for t in 2:T ] # from Fella et al. 2019
    θNs = rouwenhorst.(N, ρ, σs, stationary = false)
    return θNs
    # unnecessary to run this multiple times because transition matrix is same
end

T_ret = 40
σ_η, σ_z0
test = rouwenhorst(19, 1., σ_η, σ_z0, T_ret)
test[1].state_values
test[2].state_values
test[1].p == test[30].p

stop

plot(rouwenhorst(19,0.,sqrt(σ_ε),2.).state_values)
y = μ + ρ

#old

#
# σ_η = 0.01
# σ_ε = 0.05
#
# # Set up check
# N = 1000000
# y = sqrt(σ_ε).*randn(N)
# abs(var(y) - σ_ε) < 0.001
# tauc = tauchen(19,0.,sqrt(σ_ε))
# rouw = rouwenhorst(19,0.,sqrt(σ_ε))
# histogram(y, bins = range(tauc.state_values[1].-step(tauc.state_values)/2, step=step(tauc.state_values),
#                     length=length(tauc.state_values)+1), normalize = :probability)
# scatter!(tauc.state_values,tauc.p[1,:])
# histogram(y, bins = range(rouw.state_values[1].-step(rouw.state_values)/2, step=step(rouw.state_values),
#                     length=length(rouw.state_values)+1), normalize = :probability, alpha = 0.5)
# compare = scatter!(rouw.state_values,rouw.p[1,:],alpha=0.5)
# plot(compare)
#
#
#
#
#
# π_εr = rouwenhorst(19,0.,sqrt(σ_ε)).p
# distt = π_εt[10,:]
# distr = π_εr[10,:]
# plot()
# # histogram!(y, bins = range(-3*sqrt(σ_ε)-3*sqrt(σ_ε)/19, stop = 3*sqrt(σ_ε)+3*sqrt(σ_ε)/19, length = 20), normalize = :probability)
# histogram!(y, bins = range(-3*sqrt(var(y))-3*sqrt(var(y))/19, stop = 3*sqrt(var(y))+3*sqrt(var(y))/19, length = 20), normalize = :probability)
# scatter!(range(-3*sqrt(var(y)), stop = 3*sqrt(var(y)), length = 19),distt)
# sqrt(var(y))*sqrt(18)
# scatter!(range(-3*sqrt(var(y)), stop = 3*sqrt(var(y)), length = 19),distr)
# ylims!(0.,0.2)
# count( x -> (x < sqrt(var(y))/6 && x > -sqrt(var(y))/6), y)/N
# π_ε[10,10]
#
#
#
# # test sum of random vars
# σ_η = 0.2
# # σ_ε
# N = 1000000
# y = sqrt(σ_ε).*randn(N) .+ sqrt(σ_η).*randn(N)
# abs(var(y) - (σ_ε + σ_η)) < 0.001
# π_ε = tauchen(19,0.,sqrt(σ_ε)).p
# π_η = tauchen(19,0.,sqrt(σ_η),5).p
# dist = (π_ε * π_η)[10,:]
# # count( x -> (x < sqrt(σ_ε)/6 && x > -sqrt(σ_ε)/6), y)/N
# # count( x -> (x < sqrt(var(y))/6 && x > -sqrt(var(y))/6), y)/N
# (π_ε * π_η)[10,10]
# plot()
# histogram!(y, bins = range(-3*sqrt(var(y))-3*sqrt(var(y))/19, stop = 3*sqrt(var(y))+3*sqrt(var(y))/19, length = 20), normalize = :probability)
# scatter!(range(-3*sqrt(var(y)), stop = 3*sqrt(var(y)), length = 19),dist)
# ylims!(0.,0.2)
#
#
#
#
#
#
#
# y = sqrt(σ_ε).*randn(N) .+ sqrt(σ_η).*randn(N)
# histogram(y, bins = 100, normalize = true)
# count( x -> (x < 0 && x > -0.05), y)/N
# π_ε = rouwenhorst(19,0.,σ_ε)
# π_η = rouwenhorst(19,0.,σ_η)
#
#
# plot(π_η[1])
# π_η.p[1,:] .== π_η.p[2]
# π_η.p[10,10]


text:
# Length

using JSON
import Base: convert, *, /, +, -

export em, cm, mm, inch, pt, pc, px, ex, vw, vh, vmin, cent

@doc """
Type for representing CSS length measures.
""" ->
immutable Length{unit}
    value::Float64
end

@doc " 1mm " ->
const mm = Length{:mm}(1.0)
@doc " 1em " ->
const em = Length{:em}(1.0)
@doc " 1cm " ->
const cm = Length{:cm}(1.0)
@doc " 1inch " ->
const inch = Length{:in}(1.0)
@doc " 1pt " ->
const pt = Length{:pt}(1.0)
@doc " 1px " ->
const px = Length{:px}(1.0)
@doc " 1vw " ->
const vw = Length{:vw}(1.0)
@doc " 1vh " ->
const vh = Length{:vh}(1.0)
@doc " 1vmin " ->
const vmin = Length{:vmin}(1.0)
@doc " 1 percent " ->
const cent = Length{:cent}(1.0)

*{T <: Length}(l::T, n::Real) = T(n * l.value)
*(n::Real, l::Length) = l * n
/{T <: Length}(l::T, n::Real) = T(l.value / n)
+{T <: Length}(l::T, r::T) = T(l.value + r.value)
-{T <: Length}(l::T, r::T) = T(l.value - r.value)

convert(::Type{Length}, x::Real) =
    x * 100.0cent

Base.string{unit}(x::Length{unit}) = string(x.value, unit)
Base.string(x::Length{:cent}) = string(x.value, :%)
JSON.lower(x::Length) = string(x)


text:
using ImagePhaseCongruency
using Test

include("test_frequencyfilt.jl")
include("test_phasecongruency.jl")
include("test_syntheticimages.jl")
include("test_utilities.jl")

max_stars_repo_name: peterkovesi/ImagePhaseCongruency

max_issues_repo_name: peterkovesi/ImagePhaseCongruency

max_forks_repo_name: peterkovesi/ImagePhaseCongruency


text:
# transmcp.jl
# translated by Changhyun Kwon from transmcp.gms
# http://www.gams.com/modlib/libhtml/transmcp.htm
#
# Transportation model as equilibrium problem (TRANSMCP,SEQ=126)
#
#    Dantzig's original transportation model (TRNSPORT) is
#    reformulated as a linear complementarity problem.  We first
#    solve the model with fixed demand and supply quantities, and
#    then we incorporate price-responsiveness on both sides of the
#    market.
#
#
# Dantzig, G B, Chapter 3.3. In Linear Programming and Extensions.
# Princeton University Press, Princeton, New Jersey, 1963.

using Complementarity
using Test

@testset "transmcp_numeric.jl 1" begin

    # plants = ["seattle", "san-diego"]
    plants = 1:2
    # markets = ["new-york", "chicago", "topeka"]
    markets = 1:3

    capacity = [350, 600]
    a = Dict(zip(plants, capacity))

    demand = [325, 300, 275]
    b = Dict(zip(markets, demand))

    elasticity = [1.5, 1.2, 2.0]
    esub = Dict(zip(markets, elasticity))

    distance = [ 2.5 1.7 1.8 ;
                 2.5 1.8 1.4  ]
    d = Dict()
    for i in 1:length(plants), j in 1:length(markets)
        d[plants[i], markets[j]] = distance[i,j]
    end

    f = 90

    m = MCPModel()
    @variable(m, w[i in plants] >= 0)
    @variable(m, p[j in markets] >= 0)
    @variable(m, x[i in plants, j in markets] >= 0)

    @NLexpression(m, c[i in plants, j in markets], f * d[i,j] / 1000)
    @mapping(m, profit[i in plants, j in markets],    w[i] + c[i,j] - p[j])
    @mapping(m, supply[i in plants],                  a[i] - sum(x[i,j] for j in markets))
    @mapping(m, fxdemand[j in markets],               sum(x[i,j] for i in plants) - b[j])

    @complementarity(m, profit, x)
    @complementarity(m, supply, w)
    @complementarity(m, fxdemand, p)

    status = solveMCP(m; convergence_tolerance=1e-8, output="yes", time_limit=3600)

    @show result_value.(x)
    @show result_value.(w)
    @show result_value.(p)

    @show status
    @test status == :Solved
    @test isapprox(result_value(x[1,2]), 300.0)
    @test isapprox(result_value(p[3]),  0.126)

end

@testset "transmcp_numeric.jl 2" begin

    # plants = ["seattle", "san-diego"]
    plants = 1:2
    # markets = ["new-york", "chicago", "topeka"]
    markets = 1:3

    capacity = [350, 600]
    a = Dict(zip(plants, capacity))

    demand = [325, 300, 275]
    b = Dict(zip(markets, demand))

    elasticity = [1.5, 1.2, 2.0]
    esub = Dict(zip(markets, elasticity))

    distance = [ 2.5 1.7 1.8 ;
                 2.5 1.8 1.4  ]
    d = Dict()
    for i in 1:length(plants), j in 1:length(markets)
        d[plants[i], markets[j]] = distance[i,j]
    end

    f = 90

    m = MCPModel()
    @variable(m, w[1:2] >= 0)
    @variable(m, p[1:3] >= 0)
    @variable(m, x[1:2, 1:3] >= 0)

    @NLexpression(m, c[i in 1:2, j in 1:3], f * d[i,j] / 1000)
    @mapping(m, profit[i in 1:2, j in 1:3], w[i] + c[i,j] - p[j])
    @mapping(m, supply[i in 1:2], a[i] - sum(x[i,j] for j in markets))
    @mapping(m, fxdemand[j in 1:3], sum(x[i,j] for i in plants) - b[j])

    @complementarity(m, profit, x)
    @complementarity(m, supply, w)
    @complementarity(m, fxdemand, p)

    status = solveMCP(m; convergence_tolerance=1e-8, output="yes", time_limit=3600)

    @show result_value.(x)
    @show result_value.(w)
    @show result_value.(p)

    @show status
    @test status == :Solved
    @test isapprox(result_value(x[1,2]), 300.0)
    @test isapprox(result_value(p[3]),  0.126)

end


text:
using ContinuousPOMDPTreeSearchExperiments
using ProgressMeter
using POMCPOW
using POMDPs
@everywhere using POMDPToolbox

m = LQG1D(1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 3)
filter = Kalman1D(m)
policy = LinearFeedback(1/2)

# @show N = 1_000_000_000
@show N = 1_000_000

# XXX must be on polynomial branch of POMCPOW for this to work

ps = POMCPOWSolver(estimate_value=FORollout(policy),
                   k_action=1.0,
                   k_observation=1.0,
                   alpha_action=NaN,
                   alpha_observation=NaN,
                   check_repeat_act=true,
                   check_repeat_obs=false,
                   criterion=MaxPoly(2.0),
                   tree_queries=N,
                   max_depth=2
                  )
# pps = ParallelPOMCPOWSolver(ps, 500)
planner = solve(ps, m)
b = initial_state_distribution(m)

rng = MersenneTwister(7)

# a, info = action_info(planner, b)
# @show maximum(info[:qs])

tree = POMCPOW.make_tree(planner, b)
@showprogress for i in 1:N
    POMCPOW.simulate(planner, POWTreeObsNode(tree, 1), rand(rng, b), 2)
end

best_node = POMCPOW.select_best(MaxTries(), POWTreeObsNode(tree, 1), Base.GLOBAL_RNG)
@show tree.v[best_node]
@show [tree.v[i] for i in tree.tried[1]]
# @show [tree.a_labels[i] for i in tree.tried[1]]

#=
if !isdefined(:results)
    results = DataFrame()
end

if isempty(results)
    sims = []
    @showprogress for i in 1:N
        rng = MersenneTwister(i)
        sim = Sim(m, policy, filter, rng=rng)
        push!(sims, sim)
    end
    results = run_parallel(sims) do sim, hist
        sh = state_hist(hist)
        ah = action_hist(hist)
        oh = observation_hist(hist)
        @assert length(sh) == 3
        @assert length(ah) == 2
        return [:x1=>sh[1].x,
                :x2=>sh[2].x,
                :x3=>sh[3].x,
                :a2=>ah[2],
                :z2=>oh[1],
                :reward=>discounted_reward(hist)]
    end
end

@show var(results[:x1])
@show var(results[:x2])
@show var(results[:x3])
@show var(results[:a2])
@show mean(results[:reward])
# @show results[:a2]./results[:z2]
=#

# @time rsum = @parallel (+) for i in 1:N
#     rng = MersenneTwister(i+5_000_000)
#     ro = RolloutSimulator(rng=rng)
#     simulate(ro, m, policy, filter)
# end
# @show rsum/N


text:
### transformations for dynamic datums.  Can Proj4 do this?

#=TODO: make this
function convert{T,Y}(::Type{StaticLLA{GDA94}}, X::StaticLLA{T, Y})

    // ITRF transformation parameters from
    // J. Dawson and A. Woods, "ITRF to GDA94 coordinate transforms",
    // Journal of Applied Geodesy, 4, p. 189 (2010):
    //
    // Format of the following table is:
    //   tx ty tz sc rx ry rz (value at 1994 reference epoch)
    //   tx ty tz sc rx ry rz (derivatives)
    //
    // Units are:
    //   mm    mm    mm    ppb    mas    mas    mas
    //   mm/yr mm/yr mm/yr ppb/yr mas/yr mas/yr mas/yr
    //
    // where ppb = parts per billion; mas = milli arc seconds; yr = year
    //

    const double* params = 0;
    const double* rates = 0;
    const double itrf2008_params[7] = {-84.68, -19.42, 32.01, 9.710, -0.4254, 2.2578, 2.4015};
    const double itrf2008_rates [7] = {1.42, 1.34, 0.90, 0.109, 1.5461, 1.1820, 1.1551};
    const double itrf2005_params[7] = {-79.73, -6.86, 38.03, 6.636, -0.0351, 2.1211, 2.1411};
    const double itrf2005_rates [7] = {2.25, -0.62, -0.56, 0.294, 1.4707, 1.1443, 1.1701};
    const double itrf2000_params[7] = {-45.91, -29.85, -20.37, 7.070, -1.6705, 0.4594, 1.9356};
    const double itrf2000_rates [7] = {-4.66, 3.55, 11.24, 0.249, 1.7454, 1.4868, 1.2240};
    const double itrf1997_params[7] = {-14.63, -27.62, -25.32, 6.695, -1.7893, -0.6047, 0.9962};
    const double itrf1997_rates [7] = {-8.60, 0.36, 11.25, 0.007, 1.6394, 1.5198, 1.3801};
    const double itrf1996_params[7] = {24.54, -36.43, -68.12, 6.901, -2.7359, -2.0431, 0.3731};
    const double itrf1996_rates [7] = {-21.80, 4.71, 26.27, 0.388, 2.0203, 2.1735, 1.6290};
    if      (inputDatumName == "ITRF2008") { params = itrf2008_params; rates = itrf2008_rates; }
    else if (inputDatumName == "ITRF2005") { params = itrf2005_params; rates = itrf2005_rates; }
    else if (inputDatumName == "ITRF2000") { params = itrf2000_params; rates = itrf2000_rates; }
    else if (inputDatumName == "ITRF1997") { params = itrf1997_params; rates = itrf1997_rates; }
    else if (inputDatumName == "ITRF1996") { params = itrf1996_params; rates = itrf1996_rates; }
    else
    {
        throw RoamesError("Transform from input datum %s not supported", inputDatum);
    }
    if (!inputDatumHadTime)
    {
        throw RoamesError("Date must be provided for %s, "
                          "since the datum is time dependent", inputDatum);
    }

    // Unit conversion factors to get into meters & radians
    const double mas2rad = 1e-3/(60*60) * M_PI/180;
    const double unitConv[7] = {1e-3, 1e-3, 1e-3, 1e-9, mas2rad, mas2rad, mas2rad};
    // time is in fractional years, relative to 1994.
    double t = inputDatumTime - 1994;
    double coeffs[7] = {0};
    for (int i = 0; i < 7; ++i)
        coeffs[i] = unitConv[i] * (params[i] + rates[i]*t);
    m_Tx = coeffs[0];
    m_Ty = coeffs[1];
    m_Tz = coeffs[2];
    m_Sc = coeffs[3];
    m_Rx = coeffs[4];
    m_Ry = coeffs[5];
    m_Rz = coeffs[6];
    LOG_DEBUG("%s -> %s tansformation parameters: "
              "T = (%.5f %.5f %.5f), Sc = %.5e, R = (%.5e, %.5e, %.5e)",
              inputDatum, outputDatum, m_Tx, m_Ty, m_Tz, m_Sc, m_Rx, m_Ry, m_Rz);

=#


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export
    sssr, sssrmembership,
    fusedrings, fusedringmembership,
    sssrsizes, sssrcount,
    isringatom, isringbond,
    atomsymbol, charge, multiplicity, bondorder,
    nodedegree, valence, lonepair,
    heavyatomconnected, explicithconnected, implicithconnected,
    hydrogenconnected, connectivity,
    ishdonor, hacceptorcount, ishacceptor, hdonorcount,
    isrotatable, rotatablecount,
    atomcounter, heavyatomcount, molecularformula, empiricalformula,
    pielectron, hybridization,
    isaromaticring, isaromatic, isaromaticbond,
    precalculate!

const LONEPAIR_COUNT = Dict(
    :H => 0, :B => -1, :C => 0, :N => 1, :O => 2, :F => 3,
    :Si => 0, :P => 1, :S => 2, :Cl => 3,
    :As => 1, :Se => 2, :Br => 3, :I => 3
)

# Molecular graph topology

"""
    sssr(mol::UndirectedGraph) -> Vector{Vector{Int}}

Calculate Small set of smallest rings (SSSR).

This returns a vector of rings represented as a vector of atom node index along with the cycle path.
"""
sssr = Graph.mincycles

"""
    sssrmembership(mol::UndirectedGraph) -> Vector{Set{Int}}

Return size n vector of SSSR membership set where n is the number of atom nodes.

The numbers in the set correspond to the index of [`sssr`](@ref) that the node belongs.
"""
sssrmembership = Graph.mincyclemembership

"""
    fusedrings(mol::UndirectedGraph) -> Vector{Set{Int}}

Return the size ``n`` vector of ring labels within the molecule that have ``n``
fused rings (2-edge connected components in the context of graph theory).
"""
@cachefirst function fusedrings(mol::UndirectedGraph)
    cobr = setdiff(edgeset(mol), bridges(mol))
    subg = plaingraph(edgesubgraph(mol, cobr))
    return Graph.connectedcomponents(subg)
end

"""
    fusedringmembership(mol::UndirectedGraph) -> Vector{Int}

Return the membership of fused rings. See [`fusedrings`](@ref).
"""
@cachefirst function fusedringmembership(mol::UndirectedGraph)
    arr = [Set{Int}() for i in 1:nodecount(mol)]
    for (i, conn) in enumerate(fusedrings(mol))
        for n in conn
            push!(arr[n], i)
        end
    end
    return arr
end

@cachefirst function sssrsizes(mol::UndirectedGraph)
    ssrs_ = sssr(mol)
    arr = Set{Int}[]
    for ks in sssrmembership(mol)
        push!(arr, Set{Int}([length(ssrs_[k]) for k in ks]))
    end
    return arr
end

@cachefirst sssrcount(mol::UndirectedGraph) = length.(sssrmembership(mol))

@cachefirst isringatom(mol::UndirectedGraph) = .!isempty.(sssrmembership(mol))
isringatom(view::SubgraphView) = isringatom(view.graph)

@cachefirst isringbond(mol::UndirectedGraph
    ) = .!isempty.(edgemincyclemembership(mol))
isringbond(view::SubgraphView) = isringbond(view.graph)

# Elemental properties

"""
    atomsymbol(mol::GraphMol) -> Vector{Int}

Return the size ``n`` vector of atom symbols within the molecule that have ``n``
atoms.
"""
@cachefirst atomsymbol(mol::GraphMol) = getproperty.(nodeattrs(mol), :symbol)
atomsymbol(view::SubgraphView) = atomsymbol(view.graph)

"""
    charge(mol::GraphMol) -> Vector{Int}

Return the size ``n`` vector of atom charges within the molecule that have ``n``
atoms.
"""
@cachefirst charge(mol::GraphMol) = getproperty.(nodeattrs(mol), :charge)
charge(view::SubgraphView) = charge(view.graph)

"""
    multiplicity(mol::GraphMol) -> Vector{Int}

Return the size ``n`` vector of atom multiplicities within the molecule that
have ``n`` atoms (1: non-radical, 2: radical, 3: biradical).
"""
@cachefirst multiplicity(mol::GraphMol
    ) = getproperty.(nodeattrs(mol), :multiplicity)
multiplicity(view::SubgraphView) = multiplicity(view.graph)

"""
    bondorder(mol::GraphMol) -> Vector{Int}

Return the size ``n`` vector of bond orders within the molecule that have ``n``
bonds.
"""
@cachefirst bondorder(mol::GraphMol) = getproperty.(edgeattrs(mol), :order)
bondorder(view::SubgraphView) = bondorder(view.graph)

"""
    nodedegree(mol::GraphMol) -> Vector{Int}

Return the size ``n`` vector of node degrees within the molecular graph that
have ``n`` atom nodes.

`nodedegree` values correspond to SMARTS `D` property.
"""
@cachefirst nodedegree(mol::GraphMol
    ) = [degree(mol, n) for n in 1:nodecount(mol)]
nodedegree(view::SubgraphView) = nodedegree(view.graph)

# Valence

"""
    lonepair(mol::GraphMol) -> Vector{Union{Int,Nothing}}

Return the size ``n`` vector of the number of lone pairs within the molecule that
have ``n`` atoms.

Note that implicit hydrogens are available for only organic atoms. The lonepair
value of inorganic atoms would be `nothing`. The result can take negative value if the atom has empty valence shells.
"""
@cachefirst function lonepair(mol::GraphMol)
    vec = Union{Int,Nothing}[]
    atomsymbol_ = atomsymbol(mol)
    charge_ = charge(mol)
    for i in 1:nodecount(mol)
        num = get(LONEPAIR_COUNT, atomsymbol_[i], nothing)
        v = num === nothing ? nothing : num - charge_[i]
        push!(vec, v)
    end
    return vec
end

lonepair(view::SubgraphView) = lonepair(view.graph)

@cachefirst function apparentvalence(mol::GraphMol)
    vec = zeros(Int, nodecount(mol))
    bondorder_ = bondorder(mol)
    for i in 1:nodecount(mol)
        for (inc, adj) in neighbors(mol, i)
            vec[i] += bondorder_[inc]
        end
    end
    return vec
end

apparentvalence(view::SubgraphView) = apparentvalence(view.graph)

"""
    valence(mol::GraphMol) -> Vector{Union{Int,Nothing}}

Return the size ``n`` vector of intrinsic valences (with considering implicit
hydrogens) within the molecule that have ``n`` atoms.

`valence` values correspond to SMARTS `v` property. `valence` values of inorganic atoms would be `nothing`.
"""
@cachefirst function valence(mol::GraphMol)
    vec = Union{Int,Nothing}[]
    atomsymbol_ = atomsymbol(mol)
    lonepair_ = lonepair(mol)
    for i in 1:nodecount(mol)
        if lonepair_[i] === nothing
            push!(vec, nothing)
        elseif atomsymbol_[i] === :H
            push!(vec, 1)
        else
            push!(vec, 4 - abs(lonepair_[i]))
        end
    end
    return vec
end

valence(view::SubgraphView) = valence(view.graph)

"""
    explicithconnected(mol::GraphMol) -> Vector{Int}

Return the size ``n`` vector of the number of adjacent explicit hydrogen nodes within the molecule that have ``n`` atoms.
"""
@cachefirst function explicithconnected(mol::GraphMol)
    vec = zeros(Int, nodecount(mol))
    atomsymbol_ = atomsymbol(mol)
    for i in 1:nodecount(mol)
        for (inc, adj) in neighbors(mol, i)
            if atomsymbol_[adj] === :H
                vec[i] += 1
            end
        end
    end
    return vec
end

explicithconnected(view::SubgraphView) = explicithconnected(view.graph)

"""
    implicithconnected(mol::GraphMol) -> Vector{Int}

Return the size ``n`` vector of the number of implicit hydrogens within
the molecule that have ``n`` atoms.
"""
@cachefirst function implicithconnected(mol::GraphMol)
    hcnt = (v, av) -> v === nothing ? 0 : max(0, v - av)
    return hcnt.(valence(mol), apparentvalence(mol))
end

implicithconnected(view::SubgraphView) = implicithconnected(view.graph)

"""
    heavyatomconnected(mol::GraphMol) -> Vector{Int}

Return the size ``n`` vector of the number of adjacent non-hydrogen atoms within
the molecule that have ``n`` atoms.
"""
@cachefirst heavyatomconnected(mol::GraphMol
    ) = nodedegree(mol) - explicithconnected(mol)
heavyatomconnected(view::SubgraphView) = heavyatomconnected(view.graph)

"""
    hydrogenconnected(mol::GraphMol) -> Vector{Int}

Return the size ``n`` vector of the total number of hydrogens attached to the atom within the molecule that have ``n`` atoms.

`hydrogenconnected` values correspond to SMARTS `H` property.
"""
@cachefirst hydrogenconnected(mol::GraphMol
    ) = explicithconnected(mol) + implicithconnected(mol)
hydrogenconnected(view::SubgraphView) = hydrogenconnected(view.graph)

"""
    connectivity(mol::GraphMol) -> Vector{Int}

Return the size ``n`` vector of the total number of adjacent atoms (including implicit hydrogens) within the molecule that have ``n`` atoms.

`connectivity` values correspond to SMARTS `X` property.
"""
@cachefirst connectivity(mol::GraphMol
    ) = nodedegree(mol) + implicithconnected(mol)
connectivity(view::SubgraphView) = connectivity(view.graph)

# Hydrogen bond donor/acceptor

@cachefirst function ishacceptor(mol::GraphMol)
    ac = (sym, lp) -> lp === nothing ? false : sym in (:N, :O, :F) && lp > 0
    return ac.(atomsymbol(mol), lonepair(mol))
end

ishacceptor(view::SubgraphView) = ishacceptor(view.graph)

"""
    hacceptorcount(mol::GraphMol) -> Int

Return the number of hydrogen bond acceptors (N, O and F).
"""
hacceptorcount(mol::GraphMol) = reduce(+, ishacceptor(mol); init=0)

@cachefirst function ishdonor(mol::GraphMol)
    dc = (sym, h) -> sym in (:N, :O) && h > 0
    return dc.(atomsymbol(mol), hydrogenconnected(mol))
end

ishdonor(view::SubgraphView) = ishdonor(view.graph)

"""
    hdonorcount(mol::GraphMol) -> Int

Return the number of hydrogen bond donors (O and N attached to hydrogens).
"""
hdonorcount(mol::GraphMol) = reduce(+, ishdonor(mol); init=0)

# Rotatable bonds

"""
    isrotatable(mol::GraphMol)

Return whether the bonds are rotatable or not.
"""
@cachefirst function isrotatable(mol::GraphMol)
    nodedegree_ = nodedegree(mol)
    isringbond_ = isringbond(mol)
    bondorder_ = bondorder(mol)
    vec = Bool[]
    for (i, (u, v)) in enumerate(edgesiter(mol))
        rot = (!isringbond_[i] && bondorder_[i] == 1
            && nodedegree_[u] != 1 && nodedegree_[v] != 1)
        push!(vec, rot)
    end
    return vec
end

isrotatable(view::SubgraphView) = isrotatable(view.graph)

"""
    rotatablecount(mol::GraphMol) -> Int

Return the number of rotatable bonds.
"""
rotatablecount(mol::GraphMol) = reduce(+, isrotatable(mol); init=0)

# Composition

"""
    atomcounter(mol::GraphMol) -> Dict{Symbol,Int}

Count the number of atoms and return symbol => count dict.
"""
@cachefirst function atomcounter(mol::GraphMol)
    counter = Dict{Symbol,Int}()
    for sym in atomsymbol(mol)
        if !haskey(counter, sym)
            counter[sym] = 1
        else
            counter[sym] += 1
        end
    end
    hcnt = reduce(+, implicithconnected(mol); init=0)
    if hcnt > 0
        if !haskey(counter, :H)
            counter[:H] = hcnt
        else
            counter[:H] += hcnt
        end
    end
    return counter
end

"""
    heavyatomcount(mol::GraphMol) -> Int

Return the total number of non-hydrogen atoms.
"""
function heavyatomcount(mol::GraphMol)
    counter = atomcounter(mol)
    delete!(counter, :H)
    return reduce(+, values(counter); init=0)
end

function writeformula(counter::Dict{Symbol,Int})
    strs = []
    if haskey(counter, :C)
        push!(strs, "C")
        c = pop!(counter, :C)
        c > 1 && push!(strs, string(c))
        if haskey(counter, :H)
            push!(strs, "H")
            h = pop!(counter, :H)
            h > 1 && push!(strs, string(h))
        end
    end
    for sym in sort(collect(keys(counter)))
        push!(strs, string(sym))
        cnt = counter[sym]
        cnt > 1 && push!(strs, string(cnt))
    end
    return join(strs)
end

"""
    molecularformula(mol::GraphMol) -> String

Return the molecular formula in Hill system.
"""
function molecularformula(mol::GraphMol)
    counter = atomcounter(mol)
    return writeformula(counter)
end

"""
    empiricalformula(mol::GraphMol) -> String

Return the empirical formula in Hill system.
"""
function empiricalformula(mol::GraphMol)
    counter = atomcounter(mol)
    if length(counter) > 1
        dv = reduce(gcd, values(counter))
        for k in keys(counter)
            counter[k] = div(counter[k], dv)
        end
    end
    return writeformula(counter)
end

# Hybridization

"""
    pielectron(mol::GraphMol) -> Vector{Int}

Return the size ``n`` vector of the number of ``\\pi`` electrons (including implicit hydrogens) within the molecule that have ``n`` atoms.

The ``\\pi`` electron enumelation is based on the following simple rules.

- Any atom incident to a double bond -> +1
- Any atom incident to two double bond -> +2
- Any atom incident to a triple bond -> +1
- Uncharged N and O not incident to any multiple bonds and adjacent to atoms described above -> +2

These rules are applied for only typical organic atoms. The values for inorganic atoms will be 0.
"""
@cachefirst function pielectron(mol::GraphMol)
    atomsymbol_ = atomsymbol(mol)
    charge_ = charge(mol)
    pie_ = apparentvalence(mol) - nodedegree(mol)
    vec = zeros(Int, nodecount(mol))
    for i in 1:nodecount(mol)
        vec[i] = pie_[i]
        atomsymbol_[i] in (:N, :O) || continue
        for (inc, adj) in neighbors(mol, i)
            if pie_[i] == 0 && pie_[adj] > 0 && charge_[i] == 0
                vec[i] = 2
                break
            end
        end
    end
    return vec
end

pielectron(view::SubgraphView) = pielectron(view.graph)

"""
    hybridization(mol::GraphMol) -> Vector{Int}

Return the size ``n`` vector of orbital hybridization symbols (`:sp3`, `:sp2`, `:sp`, `:none`) within the molecule that have ``n`` atoms.

These hybridizations are for only typical organic atoms. Inorganic atoms and other orbitals like s, sp3d and sp3d2 will be `:none`.
"""
@cachefirst function hybridization(mol::GraphMol)
    atomsymbol_ = atomsymbol(mol)
    pielectron_ = pielectron(mol)
    connectivity_ = connectivity(mol)
    lonepair_ = lonepair(mol)
    vec = Symbol[]
    for i in 1:nodecount(mol)
        if lonepair_[i] === nothing || atomsymbol_[i] === :H
            push!(vec, :none)
            continue
        end
        orbitals = connectivity_[i] + lonepair_[i]
        if orbitals == 4
            if (atomsymbol_[i] in [:N, :O] && pielectron_[i] == 2)
                push!(vec, :sp2)  # adjacent to conjugated bonds
            else
                push!(vec, :sp3)
            end
        elseif orbitals == 3
            push!(vec, :sp2)
        elseif orbitals == 2
            push!(vec, :sp)
        else
            push!(vec, :none)
        end
    end
    return vec
end

hybridization(view::SubgraphView) = hybridization(view.graph)

# Aromatic rings

@cachefirst function isaromaticring(mol::GraphMol)
    atomsymbol_ = atomsymbol(mol)
    nodedegree_ = nodedegree(mol)
    pie_ = apparentvalence(mol) - nodedegree(mol)
    lonepair_ = lonepair(mol)
    carbonylO = findall(
        (atomsymbol_ .=== :O) .* (nodedegree_ .== 1) .* (pie_ .== 1))
    carbonylC = falses(nodecount(mol))
    for o in carbonylO
        c = pop!(adjacencies(mol, o))
        if atomsymbol_[c] === :C
            carbonylC[c] = true
        end
    end
    arr = falses(circuitrank(mol))
    for (i, ring) in enumerate(sssr(mol))
        cnt = 0
        for r in ring
            if carbonylC[r]
                continue
            elseif pie_[r] == 1
                cnt += 1
            elseif lonepair_[r] === nothing
                cnt = 0
                break
            elseif lonepair_[r] > 0
                cnt += 2
            elseif lonepair_[r] < 0
                continue
            else
                cnt = 0
                break
            end
        end
        if cnt % 4 == 2
            # Huckel rule check
            arr[i] = true
            continue
        end
        if nodeattrtype(mol) === SmilesAtom
            # SMILES aromatic atom
            sub = nodesubgraph(mol, Set(ring))
            if all(nodeattr(mol, n).isaromatic for n in nodeset(sub))
                arr[i] = true
            end
        end
    end
    return arr
end

"""
    isaromatic(mol::GraphMol)

Return the vector whether the atom belongs to an aromatic ring.

Note that aromaticity described here means simplified binary descriptor
(aromatic or not) based on classical Huckel's rule. This is intended for use in
some kind of pharmaceutical researches. Non-classical aromaticity such as
Moebius aromaticity is not considered.
"""
@cachefirst function isaromatic(mol::GraphMol)
    aromatic = falses(nodecount(mol))
    for ring in sssr(mol)[isaromaticring(mol)]
        sub = nodesubgraph(mol, Set(ring))
        for n in nodeset(sub)
            aromatic[n] = true
        end
    end
    return aromatic
end

isaromatic(view::SubgraphView) = isaromatic(view.graph)

@cachefirst function isaromaticbond(mol::GraphMol)
    aromaticbond = falses(edgecount(mol))
    for ring in sssr(mol)[isaromaticring(mol)]
        sub = nodesubgraph(mol, Set(ring))
        for e in edgeset(sub)
            aromaticbond[e] = true
        end
    end
    return aromaticbond
end

isaromaticbond(view::SubgraphView) = isaromaticbond(view.graph)

"""
    precalculate!(mol::GraphMol)

Convenient method to pre-calculate and cache performance bottleneck descriptors.
"""
function precalculate!(mol)
    @cache edgemincycles(mol)
    @cache mincycles(mol)
    @cache lonepair(mol)
    @cache apparentvalence(mol)
    @cache valence(mol)
    nodeattrtype(mol) === SmilesAtom && @cache coordgen(mol)
end


text: Dict{String, Any}("kilo" => Dict{String, Any}("value" => "1000", "type" => "integer"), "x" => Dict{String, Any}("value" => "1111", "type" => "integer"))


text:
"""
    MathOptLMO{OT <: MOI.Optimizer} <: LinearMinimizationOracle

Linear minimization oracle with feasible space defined through a MathOptInterface.Optimizer.
The oracle call sets the direction and reruns the optimizer.

The `direction` vector has to be set in the same order of variables as the `MOI.ListOfVariableIndices()` getter.
"""
struct MathOptLMO{OT<:MOI.AbstractOptimizer} <: LinearMinimizationOracle
    o::OT
end

function compute_extreme_point(lmo::MathOptLMO{OT}, direction::AbstractVector{T}) where {OT,T<:Real}
    variables = MOI.get(lmo.o, MOI.ListOfVariableIndices())
    terms = [MOI.ScalarAffineTerm(d, v) for (d, v) in zip(direction, variables)]
    obj = MOI.ScalarAffineFunction(terms, zero(T))
    MOI.set(lmo.o, MOI.ObjectiveFunction{typeof(obj)}(), obj)
    MOI.set(lmo.o, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    return _optimize_and_return(lmo, variables)
end

function compute_extreme_point(lmo::MathOptLMO{OT}, direction::AbstractMatrix{T}) where {OT,T<:Real}
    n = size(direction, 1)
    v = compute_extreme_point(lmo, vec(direction))
    return reshape(v, n, n)
end

function Base.copy(lmo::MathOptLMO{OT}; ensure_identity=true) where {OT}
    opt = OT() # creates the empty optimizer
    index_map = MOI.copy_to(opt, lmo.o)
    if ensure_identity
        for (src_idx, des_idx) in index_map.var_map
            if src_idx != des_idx
                error("Mapping of variables is not identity")
            end
        end
    end
    return MathOptLMO(opt)
end

function Base.copy(lmo::MathOptLMO{OT}; ensure_identity=true) where {OTI, OT <: MOIU.CachingOptimizer{OTI}}
    opt = MOIU.CachingOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
        OTI(),
    )
    index_map = MOI.copy_to(opt, lmo.o)
    if ensure_identity
        for (src_idx, des_idx) in index_map.var_map
            if src_idx != des_idx
                error("Mapping of variables is not identity")
            end
        end
    end
    return MathOptLMO(opt)
end

function compute_extreme_point(
    lmo::MathOptLMO{OT},
    direction::AbstractVector{MOI.ScalarAffineTerm{T}},
) where {OT,T}
    variables = [term.variable for term in direction]
    obj = MOI.ScalarAffineFunction(direction, zero(T))
    MOI.set(lmo.o, MOI.ObjectiveFunction{typeof(obj)}(), obj)
    MOI.set(lmo.o, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    return _optimize_and_return(lmo, variables)
end

function _optimize_and_return(lmo, variables)
    MOI.optimize!(lmo.o)
    term_st = MOI.get(lmo.o, MOI.TerminationStatus())
    if term_st ∉ (MOI.OPTIMAL, MOI.ALMOST_OPTIMAL, MOI.SLOW_PROGRESS)
        @error "Unexpected termination: $term_st"
        return MOI.get.(lmo.o, MOI.VariablePrimal(), variables)
    end
    return MOI.get.(lmo.o, MOI.VariablePrimal(), variables)
end

"""
    convert_mathopt(lmo::LMO, optimizer::OT; kwargs...) -> MathOptLMO{OT}

Converts the given LMO to its equivalent MathOptInterface representation using `optimizer`.
Must be implemented by LMOs.
"""
function convert_mathopt end


text:
using GeometricFlux
using GraphSignals
using Flux
using Flux: onehotbatch, onecold, logitcrossentropy, throttle
using Flux: @epochs
using JLD2
using Statistics
using SparseArrays
using Graphs.SimpleGraphs
using CUDA

@load "data/cora_features.jld2" features
@load "data/cora_labels.jld2" labels
@load "data/cora_graph.jld2" g

num_nodes = 2708
num_features = 1433
hidden = 16
target_catg = 7
epochs = 200
λ = 5e-4

## Preprocessing data
train_X = Matrix{Float32}(features) |> gpu  # dim: num_features * num_nodes
train_y = Matrix{Float32}(labels) |> gpu  # dim: target_catg * num_nodes
fg = FeaturedGraph(g)  # pass to gpu together in model layers

## Model
model = Chain(GCNConv(fg, num_features=>hidden, relu),
              Dropout(0.5),
              GCNConv(fg, hidden=>target_catg),
              ) |> gpu;
# do not show model architecture, showing CuSparseMatrix will trigger errors

## Loss
l2norm(x) = sum(abs2, x)
# cross entropy with first layer L2 regularization
loss(x, y) = logitcrossentropy(model(x), y) + λ*sum(l2norm, Flux.params(model[1]))
accuracy(x, y) = mean(onecold(softmax(cpu(model(x)))) .== onecold(cpu(y)))

## Training
ps = Flux.params(model)
train_data = [(train_X, train_y)]
opt = ADAM(0.01)
evalcb() = @show(accuracy(train_X, train_y))

@epochs epochs Flux.train!(loss, ps, train_data, opt, cb=throttle(evalcb, 10))


text:
using Metatheory.Rewriters

function EGraphs.preprocess(t::Symbolic)
    toterm(unflatten(t))
end

function symbolicegraph(ex)
    g = EGraph(ex)
    analyze!(g, SymbolicUtils.SymtypeAnalysis)
    settermtype!(g, Term{symtype(ex), Any})
    return g
end

"""
Equational rewrite rules for optimizing expressions
"""
opt_theory = @theory a b x y  begin
    a + b == b + a
    a * b == b * a
    a * x + a * y == a*(x+y)
    -1 * a == -a
    a + (-1 * b) == a - b
    x^-1 == 1/x
    1/x * a == a/x
    # fraction rules
    # (a/b) + (c/b) => (a+c)/b
    # trig functions
    sin(x)/cos(x) == tan(x)
    cos(x)/sin(x) == cot(x)
    sin(x)^2 + cos(x)^2 --> 1
    sin(2a) == 2sin(a)cos(a)
end

"""
Approximation of costs of operators in number
of CPU cycles required for the numerical computation

See
 * https://latkin.org/blog/2014/11/09/a-simple-benchmark-of-various-math-operations/
 * https://streamhpc.com/blog/2012-07-16/how-expensive-is-an-operation-on-a-cpu/
 * https://github.com/triscale-innov/GFlops.jl
"""
const op_costs = Dict(
    (+)     => 1,
    (-)     => 1,
    abs     => 2,
    (*)     => 3,
    exp     => 18,
    (/)     => 24,
    (^)     => 100,
    log1p   => 124,
    deg2rad => 125,
    rad2deg => 125,
    acos    => 127,
    asind   => 128,
    acsch   => 133,
    sin     => 134,
    cos     => 134,
    atan    => 135,
    tan     => 156,
)
# TODO some operator costs are in FLOP and not in cycles!!

function costfun(n::ENodeTerm, g::EGraph, an)
    op = operation(n)
    cost = 0
    cost += get(op_costs, op, 1)

    for id ∈ n.args
        eclass = g[id]
        !hasdata(eclass, an) && (cost += Inf; break)
        cost += last(getdata(eclass, an))
    end
    cost
end

costfun(n::ENodeLiteral, g::EGraph, an) = 0

egraph_simterm(x, head, args, symtype=nothing; metadata=nothing, exprhead=exprhead(x)) =
TermInterface.similarterm(typeof(x), head, args, symtype; metadata=metadata, exprhead=exprhead)

# Custom similarterm to use in EGraphs on <:Symbolic types that treats everything as a Term
function egraph_simterm(x::Type{<:Term}, f, args, symtype=nothing; metadata=nothing, exprhead=:call)
    T = symtype
    if T === nothing
        T = _promote_symtype(f, args)
    end
    res = Term{T}(f isa Symbol ? eval(f) : f, args; metadata=metadata);
    return res
end

function optimize(ex; params=SaturationParams(timeout=20))
    # @show ex
    g = symbolicegraph(ex)
    params.simterm = egraph_simterm
    saturate!(g, opt_theory, params)
    return extract!(g, costfun)
end


text:
# CUDA Execution control

get_dim_x(g::Int) = g
get_dim_x(g::(Int, Int)) = g[1]
get_dim_x(g::(Int, Int, Int)) = g[1]

get_dim_y(g::Int) = 1
get_dim_y(g::(Int, Int)) = g[2]
get_dim_y(g::(Int, Int, Int)) = g[2]

get_dim_z(g::Int) = 1
get_dim_z(g::(Int, Int)) = 1
get_dim_z(g::(Int, Int, Int)) = g[3]

typealias CuDim Union(Int, (Int, Int), (Int, Int, Int))

# Stream management

function launch(f::CuFunction, grid::CuDim, block::CuDim, args::Tuple; shmem_bytes::Int=4, stream::CuStream=null_stream())
	gx = get_dim_x(grid)
	gy = get_dim_y(grid)
	gz = get_dim_z(grid)

	tx = get_dim_x(block)
	ty = get_dim_y(block)
	tz = get_dim_z(block)

	kernel_args = [cubox(arg) for arg in args]

	@cucall(:cuLaunchKernel, (
		Ptr{Void},  # function
		Cuint,  # grid dim x
		Cuint,  # grid dim y
		Cuint,  # grid dim z
		Cuint,  # block dim x
		Cuint,  # block dim y
		Cuint,  # block dim z
		Cuint,  # shared memory bytes,
		Ptr{Void}, # stream
		Ptr{Ptr{Void}}, # kernel parameters,
		Ptr{Ptr{Void}}), # extra parameters
		f.handle, gx, gy, gz, tx, ty, tz, shmem_bytes, stream.handle, kernel_args, 0)
end


text:
module Editors

using Genie, Stipple, StippleUI, StippleUI.API
import Genie.Renderer.Html: HTMLString, normal_element, select, template

export editor

function __init__()
  Genie.Renderer.Html.register_normal_element("q__editor", context = Genie.Renderer.Html)
end

"""
  editor(fieldname, args...; wrap, kwargs...)

Creates a Quasar editor view editing the content stored in the fieldname property
of the model.

Some optional arguments are:
    definitions : A dict or list of toolbar command descriptions
    toolbar: A dict or list describing the arrangement of the toolbar elements
    autocorrect: "on/off", whever to activate autocorrect
    autocapitalize: "on/off"
    autocomplete: "on/off
    spellcheck: "true/false"
"""
function editor(fieldname::Symbol,
                args...;
                wrap::Function = StippleUI.DEFAULT_WRAPPER,
                kwargs...)
  wrap() do
    Genie.Renderer.Html.q__editor(args...; attributes([:fieldname => fieldname, kwargs...], StippleUI.API.ATTRIBUTES_MAPPINGS)...)
  end
end

end


text:
function string_zeros(s::AbstractString)
x = Array(s=="Int64"?Int64:Float64, 1_000_000)
for i in 1:length(x)
x[i] = 0
end
return x
end

max_stars_repo_path: Module 2/Chapter03/stringzero.jl

max_issues_repo_path: Module 2/Chapter03/stringzero.jl

max_forks_repo_path: Module 2/Chapter03/stringzero.jl


text:
export BrickletAnalogInVoltageCallbackThreshold
struct BrickletAnalogInVoltageCallbackThreshold
    option::Char
    min::Integer
    max::Integer
end

export BrickletAnalogInAnalogValueCallbackThreshold
struct BrickletAnalogInAnalogValueCallbackThreshold
    option::Char
    min::Integer
    max::Integer
end

export BrickletAnalogInIdentity
struct BrickletAnalogInIdentity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletAnalogIn
"""
Measures DC voltage between 0V and 45V
"""
mutable struct BrickletAnalogIn <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletAnalogIn(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_analog_in")
        deviceInternal = package.BrickletAnalogIn(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export get_voltage
"""
    $(SIGNATURES)

Returns the voltage of the sensor. The resolution between 0 and 6V is about 2mV.
Between 6 and 45V the resolution is about 10mV.

If you want to get the voltage periodically, it is recommended to use the
:cb:`Voltage` callback and set the period with
:func:`Set Voltage Callback Period`.
"""
function get_voltage(device::BrickletAnalogIn)
    return device.deviceInternal.get_voltage()
end

export get_analog_value
"""
    $(SIGNATURES)

Returns the value as read by a 12-bit analog-to-digital converter.

.. note::
 The value returned by :func:`Get Voltage` is averaged over several samples
 to yield less noise, while :func:`Get Analog Value` gives back raw
 unfiltered analog values. The only reason to use :func:`Get Analog Value` is,
 if you need the full resolution of the analog-to-digital converter.

If you want the analog value periodically, it is recommended to use the
:cb:`Analog Value` callback and set the period with
:func:`Set Analog Value Callback Period`.
"""
function get_analog_value(device::BrickletAnalogIn)
    return device.deviceInternal.get_analog_value()
end

export set_voltage_callback_period
"""
    $(SIGNATURES)

Sets the period with which the :cb:`Voltage` callback is triggered
periodically. A value of 0 turns the callback off.

The :cb:`Voltage` callback is only triggered if the voltage has changed since
the last triggering.
"""
function set_voltage_callback_period(device::BrickletAnalogIn, period)
    device.deviceInternal.set_voltage_callback_period(period)
end

export get_voltage_callback_period
"""
    $(SIGNATURES)

Returns the period as set by :func:`Set Voltage Callback Period`.
"""
function get_voltage_callback_period(device::BrickletAnalogIn)
    return device.deviceInternal.get_voltage_callback_period()
end

export set_analog_value_callback_period
"""
    $(SIGNATURES)

Sets the period with which the :cb:`Analog Value` callback is triggered
periodically. A value of 0 turns the callback off.

The :cb:`Analog Value` callback is only triggered if the analog value has
changed since the last triggering.
"""
function set_analog_value_callback_period(device::BrickletAnalogIn, period)
    device.deviceInternal.set_analog_value_callback_period(period)
end

export get_analog_value_callback_period
"""
    $(SIGNATURES)

Returns the period as set by :func:`Set Analog Value Callback Period`.
"""
function get_analog_value_callback_period(device::BrickletAnalogIn)
    return device.deviceInternal.get_analog_value_callback_period()
end

export set_voltage_callback_threshold
"""
    $(SIGNATURES)

Sets the thresholds for the :cb:`Voltage Reached` callback.

The following options are possible:

.. csv-table::
 :header: "Option", "Description"
 :widths: 10, 100

 "'x'",    "Callback is turned off"
 "'o'",    "Callback is triggered when the voltage is *outside* the min and max values"
 "'i'",    "Callback is triggered when the voltage is *inside* the min and max values"
 "'<'",    "Callback is triggered when the voltage is smaller than the min value (max is ignored)"
 "'>'",    "Callback is triggered when the voltage is greater than the min value (max is ignored)"
"""
function set_voltage_callback_threshold(device::BrickletAnalogIn, option, min, max)
    device.deviceInternal.set_voltage_callback_threshold(option, min, max)
end

export get_voltage_callback_threshold
"""
    $(SIGNATURES)

Returns the threshold as set by :func:`Set Voltage Callback Threshold`.
"""
function get_voltage_callback_threshold(device::BrickletAnalogIn)
    return device.deviceInternal.get_voltage_callback_threshold()
end

export set_analog_value_callback_threshold
"""
    $(SIGNATURES)

Sets the thresholds for the :cb:`Analog Value Reached` callback.

The following options are possible:

.. csv-table::
 :header: "Option", "Description"
 :widths: 10, 100

 "'x'",    "Callback is turned off"
 "'o'",    "Callback is triggered when the analog value is *outside* the min and max values"
 "'i'",    "Callback is triggered when the analog value is *inside* the min and max values"
 "'<'",    "Callback is triggered when the analog value is smaller than the min value (max is ignored)"
 "'>'",    "Callback is triggered when the analog value is greater than the min value (max is ignored)"
"""
function set_analog_value_callback_threshold(device::BrickletAnalogIn, option, min, max)
    device.deviceInternal.set_analog_value_callback_threshold(option, min, max)
end

export get_analog_value_callback_threshold
"""
    $(SIGNATURES)

Returns the threshold as set by :func:`Set Analog Value Callback Threshold`.
"""
function get_analog_value_callback_threshold(device::BrickletAnalogIn)
    return device.deviceInternal.get_analog_value_callback_threshold()
end

export set_debounce_period
"""
    $(SIGNATURES)

Sets the period with which the threshold callbacks

* :cb:`Voltage Reached`,
* :cb:`Analog Value Reached`

are triggered, if the thresholds

* :func:`Set Voltage Callback Threshold`,
* :func:`Set Analog Value Callback Threshold`

keep being reached.
"""
function set_debounce_period(device::BrickletAnalogIn, debounce)
    device.deviceInternal.set_debounce_period(debounce)
end

export get_debounce_period
"""
    $(SIGNATURES)

Returns the debounce period as set by :func:`Set Debounce Period`.
"""
function get_debounce_period(device::BrickletAnalogIn)
    return device.deviceInternal.get_debounce_period()
end

export set_range
"""
    $(SIGNATURES)

Sets the measurement range. Possible ranges:

* 0: Automatically switched
* 1: 0V - 6.05V, ~1.48mV resolution
* 2: 0V - 10.32V, ~2.52mV resolution
* 3: 0V - 36.30V, ~8.86mV resolution
* 4: 0V - 45.00V, ~11.25mV resolution
* 5: 0V - 3.3V, ~0.81mV resolution, new in version 2.0.3\$nbsp;(Plugin)

.. versionadded:: 2.0.1\$nbsp;(Plugin)
"""
function set_range(device::BrickletAnalogIn, range)
    device.deviceInternal.set_range(range)
end

export get_range
"""
    $(SIGNATURES)

Returns the measurement range as set by :func:`Set Range`.

.. versionadded:: 2.0.1\$nbsp;(Plugin)
"""
function get_range(device::BrickletAnalogIn)
    return device.deviceInternal.get_range()
end

export set_averaging
"""
    $(SIGNATURES)

Set the length of a averaging for the voltage value.

Setting the length to 0 will turn the averaging completely off. If the
averaging is off, there is more noise on the data, but the data is without
delay.

.. versionadded:: 2.0.3\$nbsp;(Plugin)
"""
function set_averaging(device::BrickletAnalogIn, average)
    device.deviceInternal.set_averaging(average)
end

export get_averaging
"""
    $(SIGNATURES)

Returns the averaging configuration as set by :func:`Set Averaging`.

.. versionadded:: 2.0.3\$nbsp;(Plugin)
"""
function get_averaging(device::BrickletAnalogIn)
    return device.deviceInternal.get_averaging()
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletAnalogIn)
    return device.deviceInternal.get_identity()
end

export register_callback
"""
Registers the given *function* with the given *callback_id*.
"""
function register_callback(device::BrickletAnalogIn, callback_id, function_)
    if isnothing(function_)
        device.registered_callbacks.pop(callback_id, None)
    else
        device.registered_callbacks[callback_id] = function_
    end
end


text:
function FMM_mainMVP_pre!(output, beta, scatteringMatrices, φs::Vector{Float64}, ids::Vector{Int}, P, mFMM, pre_agg, translated_sum)
    #@simd does not have a positive effect in 0.6.0
    #calculate matrix-vector product - devectorized with pre-preconditioning

    @inbounds begin
    mul!(output, mFMM.Znear, beta)
    G = length(mFMM.groups)
    fill!(pre_agg,0.0)
    #first preagg all
    for ig2 = 1:G
        for is = 1:mFMM.groups[ig2].size
            indices = (mFMM.groups[ig2].point_ids[is]-1)*(2*P+1)
            for ii = 1:2*P+1
                for iq = 1:mFMM.Q
                    pre_agg[iq,ig2] += mFMM.Agg[ig2][iq,(is-1)*(2*P+1) + ii]*beta[indices + ii]
                end
            end
        end
    end

    for ig1 = 1:G
        #translate plane waves from ig2 to ig1
        fill!(translated_sum,0.0)
        for ig2 = 1:G
            if isempty(mFMM.Trans[(ig1-1)*G + ig2])
                continue
            else
                for iQ = 1:mFMM.Q
                    translated_sum[iQ] -= mFMM.Trans[(ig1-1)*G + ig2][iQ]*pre_agg[iQ,ig2]
                end
                #minus sign because the real equation is (I-ST)x=b
            end
        end
        #disaggregate from ig1 center to ig1's scatterers, producing -Tx
        for is = 1:mFMM.groups[ig1].size
            for ip = 1:2*P+1
                disagged = 0.0im
                for iq = 1:mFMM.Q
                    disagged += conj(mFMM.Agg[ig1][iq,(is-1)*(2*P+1) + ip])*translated_sum[iq]
                end
                output[(mFMM.groups[ig1].point_ids[is]-1)*(2*P+1) + ip] += disagged
            end
        end
    end
    #multiply by S to produce -S
    #temp can be moved outward, but for now this preallocation prevents most
    #dynamic mem alloc by this MVP
    temp = Array{Complex{Float64}}(undef, 2*P+1)
    for ic = 1:length(ids)
        rng = (ic-1)*(2*P+1) .+ (1:2*P+1)
        #copy values to avoid dynamic allocation - copyto!(temp,1:2*P+1,output,rng) doesn't work because output is a subarray
        for ip = 1:2*P+1
            temp[ip] = output[(ic-1)*(2*P+1) + ip]
        end
        v = view(output,rng)
        if φs[ic] == 0.0
            mul!(v, scatteringMatrices[ids[ic]], temp)
        else
            #rotate without matrix
            rotateMultipole!(temp,-φs[ic],P)
            mul!(v, scatteringMatrices[ids[ic]], temp)
            rotateMultipole!(v,φs[ic],P)
        end
    end
    #add identity matrix (Ix)
    output .+= beta
    end #inbounds
    return output
end

function FMM_mainMVP_transpose!(output, beta, scatteringMatrices, φs::Vector{Float64}, ids::Vector{Int}, P, mFMM, pre_agg, translated_sum)
    #calculate matrix^T-vector product - devectorized with pre-preconditioning

    @inbounds begin
    #here we first multiply by Xᵀ

    #Xβ storage can be moved outward, but for now this preallocation prevents most
    #dynamic mem alloc by this MVP
    Xβ = Array{Complex{Float64}}(undef, length(beta))
    temp = Array{Complex{Float64}}(undef, 2*P+1)
    for ic = 1:length(ids)
        rng = (ic-1)*(2*P+1) .+ (1:2*P+1)
        v = view(Xβ, rng)
        if φs[ic] == 0.0
            mul!(v, transpose(scatteringMatrices[ids[ic]]), view(beta, rng))
        else
            #rotate without matrix - transposed
            rotateMultipole!(temp, view(beta, rng), φs[ic], P)
            mul!(v, transpose(scatteringMatrices[ids[ic]]), temp)
            rotateMultipole!(v, -φs[ic], P)
        end
    end

    mul!(output, transpose(mFMM.Znear), Xβ)
    G = length(mFMM.groups)
    fill!(pre_agg,0.0)
    #first preagg all
    for ig2 = 1:G
        for is = 1:mFMM.groups[ig2].size
            indices = (mFMM.groups[ig2].point_ids[is]-1)*(2*P+1)
            for ii = 1:2*P+1
                for iq = 1:mFMM.Q
                    pre_agg[iq,ig2] += conj(mFMM.Agg[ig2][iq,(is-1)*(2*P+1) + ii])*Xβ[indices + ii]
                end
            end
        end
    end

    for ig1 = 1:G
        #translate plane waves from ig2 to ig1
        fill!(translated_sum,0.0)
        for ig2 = 1:G
            if isempty(mFMM.Trans[(ig2-1)*G + ig1]) #transposed
                continue
            else
                for iQ = 1:mFMM.Q
                    translated_sum[iQ] -= mFMM.Trans[(ig2-1)*G + ig1][iQ]*pre_agg[iQ,ig2] #transposed
                end
                #minus sign because the real equation is (I-TᵀXᵀ)β=b
            end
        end
        #disaggregate from ig1 center to ig1's scatterers
        for is = 1:mFMM.groups[ig1].size
            for ip = 1:2*P+1
                disagged = 0.0im
                for iq = 1:mFMM.Q
                    disagged += mFMM.Agg[ig1][iq,(is-1)*(2*P+1) + ip]*translated_sum[iq]
                end
                output[(mFMM.groups[ig1].point_ids[is]-1)*(2*P+1) + ip] += disagged
            end
        end
    end

    #add identity matrix (Iβ)
    output .+= beta
    end #inbounds
    return output
end


text:
using Coverage

coverage_src = process_folder("../src")
coverage_test = process_folder("../test")
coverage = merge_coverage_counts(coverage_src, coverage_test)
LCOV.writefile("htmlcov/coverage.info", coverage)


text:
module analyticalmaterials

using ..utilities

export θ_plasmone, n_drude, θ_grating_coupling, ϵ, n, hs_Palm

function ϵ(n::Number)::Number
    real(n)^2 - imag(n)^2 + 1im * 2 * real(n) * imag(n)
end

function n(ϵ::Number)::Number
    √ϵ
end

function n_drude(σ_0::Real, τ::Real, λ::Real)::Number
    ω = 2 * π * c_0 / λ
    √ϵ_c_drude(σ_0, τ, ω)
end

function ϵ_c_drude(σ_0::Real, τ::Real, ω::Real)::Number
    # See Saleh & Teich 3.ed eq.8.2-11

    ω_p = √(σ_0 / (ϵ_0 * τ))
    ϵ_0 * (1 + ω_p^2 / (-(ω^2) + 1im * ω / τ))
end

function β_plasmone(n1::Number, n2::Number, k_0::Number)::Number
    # See Maier eq. 2.14 and eq. 1.11a and b

    ϵ_1 = ϵ(n1)
    ϵ_2 = ϵ(n2)
    k_0 * √(ϵ_1 * ϵ_2 / (ϵ_1 + ϵ_2))
end

function θ_plasmone(ni::Number, n1::Number, n2::Number, k_0::Number)::Number
    β = real(β_plasmone(n1, n2, k_0))
    ki = real(ni) * k_0

    if ((β / ki) > 1)
        return NaN
    end
    asin(β / ki)
end

function θ_grating_coupling(λ, neff_core, n_cladding, Λ)
    asin(((λ * 2 / Λ) - neff_core) / n_cladding)
end

function θ_grating_coupling(λ, n_core)
    # A default for the grating in question
    θ_grating_coupling(λ, 1.4676, n_core, 1073.81e-9)
end

module hs_Palm
# An approximation of the nonlinear function h(c) such that ϵ(c) = h(c)ϵ(0) for palladium and gold-palladium alloys.
# Estimate based on linear approximation from Palm2019 and theory form Tobiska2001.
# This module should only be used as a first approximation, not for spesific optimization.

using Interpolations
using DelimitedFiles
using Match

import ..ϵ

export h

rawdatafile034 = readdlm("materials/palm_PdAu-H/HM_Pd034.csv", ',', skipstart=2)
rawdatafile042 = readdlm("materials/palm_PdAu-H/HM_Pd042.csv", ',', skipstart=2)
rawdatafile052 = readdlm("materials/palm_PdAu-H/HM_Pd052.csv", ',', skipstart=2)
rawdatafile073 = readdlm("materials/palm_PdAu-H/HM_Pd073.csv", ',', skipstart=2)
rawdatafile100 = readdlm("materials/palm_PdAu-H/HM_Pd100.csv", ',', skipstart=2)

function loadLine(filerawdata, λ::Float64, Pd_c::Float64)::Function
    point = @match λ begin
        300e-9 => 1
        600e-9 => 2
        900e-9 => 3
        1200e-9 => 4
        1500e-9 => 5
        _ => 0
    end
    offsett = 10

    raw_n = filerawdata[:, (point*2 - 1):(point*2)]
    raw_k = filerawdata[:, (point*2 - 1 + offsett):(point*2 + offsett)]

    raw_n = raw_n[(raw_n[:, 1] .!= ""), :]
    raw_k = raw_k[(raw_k[:, 1] .!= ""), :]

    sorted_n = raw_n[sortperm(raw_n[:, 1]), :]
    sorted_k = raw_k[sortperm(raw_k[:, 1]), :]

    n_estimator = LinearInterpolation(sorted_n[:, 1], sorted_n[:, 2])
    k_estimator = LinearInterpolation(sorted_k[:, 1], sorted_k[:, 2])

    h_total(c) = ϵ(n_estimator(c_to_HM(c, Pd_c)) + k_estimator(c_to_HM(c, Pd_c)) * 1im) / ϵ(sorted_n[1, 2] + sorted_k[1, 2] * 1im)

    h_total
end

function c_to_HM(H_c::Float64, Pd_c::Float64)::Float64
    # Based on Palm2019 fig. S7
    # Assumes concentration under 1 atmosphere pressure
    # Assumes the effect of a low pressure of pure hydrogen is the same as the equivalent partial pressure hydrogen in an otherwize inert atmosphere totalling 1 atmosphere
    # Assumes the relation between HM and H_c is proportional to the square root of H_c, as shown in separate note
    # H_c is the partial pressure of hydrogen in the atmosphere given in atmospheres, the approximation is only valid for H_c ⪅ 0.25
    # Pd_c is the atomic alloying ratio between Pd and Au, and should be between 0.4 and 1
    # The following data was extracted from the figure:
    # 41.50974025974026	0.0011235955056179137
    # 99.9512987012987	0.6089887640449437

    HM_Pd_c025 = 1.04012484394507 * Pd_c - 0.43062952559301

    sqrt(H_c * 1.01325) * HM_Pd_c025 / sqrt(0.25)
end

function h(λ::Float64, H_c::Float64, Pd_c::Float64)::Number
    # H_c and Pd_c is the hydrogen concentration and the Pd atomic concentration in the Pd-Au alloy, both between 0 and 1

    filerawdata = @match Pd_c begin
        0.34 => rawdatafile034
        0.42 => rawdatafile042
        0.52 => rawdatafile052
        0.73 => rawdatafile073
        1.0 => rawdatafile100
        _ => "Illegal Pd_c"
    end

    λ_lower = 300e-9
    λ_upper = 1500e-9
    if 300e-9 <= λ <= 600e-9
        λ_lower = 300e-9
        λ_upper = 600e-9
    elseif 600e-9 <= λ <= 900e-9
        λ_lower = 600e-9
        λ_upper = 900e-9
    elseif 900e-9 <= λ <= 1200e-9
        λ_lower = 900e-9
        λ_upper = 1200e-9
    elseif 1200e-9 <= λ <= 1500e-9
        λ_lower = 1200e-9
        λ_upper = 1500e-9
    end

    h_lower = loadLine(filerawdata, λ_lower, Pd_c)(H_c)
    h_upper = loadLine(filerawdata, λ_upper, Pd_c)(H_c)

    λ_partial = (λ - λ_lower) / (λ_upper - λ_lower)
    h_real = (h_upper - h_lower) * λ_partial + h_lower

    h_real
end

end

end # analyticalmaterials


text:
using DelimitedFiles

inputs = readdlm("input.txt", '\n', Int)
bigger_count = 0
for i in 2:size(inputs)[1]
    if inputs[i] > inputs[i-1]
        global bigger_count += 1
    end
end

print("Number of increases ",bigger_count)


text:
include("../src/TrueSkill.jl")
using .TrueSkill
global const ttt = TrueSkill
using Test
using CSV
using DataFrames

@testset "Examples" begin

@testset "same_strength" begin
    events = [ [["aj"],["bj"]],[["bj"],["cj"]], [["cj"],["aj"]] ,[["aj"],["bj"]],[["bj"],["cj"]], [["cj"],["aj"]]]
    results = [[1.,0.],[1.,0.],[1.,0.],[1.,0.],[1.,0.],[1.,0.]]
    priors = Dict{String,ttt.Player}()
    for k in ["aj", "bj", "cj"]
        priors[k] = ttt.Player(ttt.Gaussian(0., 3.), 0.5, 0.0)
    end
    h = ttt.History(events, results, [1,2,3,4,5,6], priors)
    lc_ts = ttt.learning_curves(h)
    ttt.convergence(h)
    lc_ttt = ttt.learning_curves(h)

    df = DataFrame(mu_a_ts = [ N.mu for (k,N) in lc_ts["aj"]]
                    ,sigma_a_ts = [ N.sigma for (k,N) in lc_ts["aj"]]
                    ,mu_a_ttt = [ N.mu for (k,N) in lc_ttt["aj"]]
                    ,sigma_a_ttt = [ N.sigma for (k,N) in  lc_ttt["aj"]]
                    ,time_a = [ k for (k,N) in lc_ts["aj"]]
                    ,mu_b_ts = [ N.mu for (k,N) in lc_ts["bj"]]
                    ,sigma_b_ts = [ N.sigma for (k,N) in lc_ts["bj"]]
                    ,mu_b_ttt = [ N.mu for (k,N) in lc_ttt["bj"]]
                    ,sigma_b_ttt = [ N.sigma for (k,N) in lc_ttt["bj"]]
                    ,time_b = [ k for (k,N) in lc_ts["bj"]]
                    ,mu_c_ts = [ N.mu for (k,N) in lc_ts["cj"]]
                    ,sigma_c_ts = [ N.sigma for (k,N) in lc_ts["cj"]]
                    ,mu_c_ttt = [ N.mu for (k,N) in lc_ttt["cj"]]
                    ,sigma_c_ttt = [ N.sigma for (k,N) in lc_ttt["cj"]]
                    ,time_c = [ k for (k,N) in lc_ts["cj"]]
                    )

    CSV.write("output/same_strength.csv", df; header=true)
    @test true
end
@testset "same_strength_two_groups" begin
    predicciones_mle = Float64[]
    predicciones_aiaj = Float64[]
    predicciones_bibj = Float64[]
    predicciones_aibi = Float64[]
    predicciones_bici = Float64[]
    predicciones_ajbj = Float64[]
    predicciones_bjcj = Float64[]
    composition = [ [["aj"],["bj"]], [["bj"],["cj"]], [["cj"],["aj"]],[["aj"],["bj"]], [["bj"],["cj"]], [["cj"],["aj"]],
                [["ai"],["bi"]], [["bi"],["ci"]], [["ci"],["ai"]], [["ai"],["bi"]], [["bi"],["ci"]], [["ci"],["ai"]],
                [["aj"],["ai"]] ]
    times = [1,1,1,1,1,1,1,1,1,1,1,1,1]
    Nbeta = ttt.Gaussian(0.0,0.5)

    for i in 1:150
        push!(composition, [["ai"],["aj"]])
        push!(times, 1)
        priors = Dict{String,ttt.Player}()
        for k in ["aj", "bj", "cj", "ai", "bi", "ci"]
            priors[k] = ttt.Player(ttt.Gaussian(0., 5.0), 0.5, 0.0)
        end

        h = ttt.History(composition=composition, times=times, priors=priors, iter=40)
        ttt.convergence(h)

        lc = ttt.learning_curves(h)

        push!(predicciones_mle, i/(1.0+i))

        push!(predicciones_aiaj, 1-ttt.cdf(lc["ai"][1][2]+Nbeta - lc["aj"][1][2]+Nbeta, 0.))
        push!(predicciones_bibj, 1-ttt.cdf(lc["bi"][1][2]+Nbeta - lc["bj"][1][2]+Nbeta, 0.))
        push!(predicciones_aibi, 1-ttt.cdf(lc["ai"][1][2]+Nbeta - lc["bi"][1][2]+Nbeta, 0.))
        push!(predicciones_bici, 1-ttt.cdf(lc["bi"][1][2]+Nbeta - lc["ci"][1][2]+Nbeta, 0.))
        push!(predicciones_ajbj, 1-ttt.cdf(lc["aj"][1][2]+Nbeta - lc["bj"][1][2]+Nbeta, 0.))
        push!(predicciones_bjcj, 1-ttt.cdf(lc["bj"][1][2]+Nbeta - lc["cj"][1][2]+Nbeta, 0.))
    end

    # S - Shape
    using Plots
    plot(predicciones_mle,predicciones_mle,legend=(-10.0, 10.0))
    plot!(predicciones_mle,predicciones_aiaj)
    #plot!(predicciones_mle,predicciones_bibj)

    df = DataFrame(mle = predicciones_mle
                    ,aiaj = predicciones_aiaj
                    ,bibj = predicciones_bibj
                    ,aibi = predicciones_aibi
                    ,bici = predicciones_bici
                    )

    CSV.write("output/same_strength_two_groups.csv", df; header=true)
    @test true
end
@testset "smoothing" begin
    events = [ [["a"],["b"]], [["a"],["b"]]]
    results = [[1.,0.],[0.,1.]]
    times = [1,2]
    priors = Dict{String,ttt.Player}()
    Nbeta = ttt.Gaussian(0.,0.5)
    for k in ["a", "b"]
        priors[k] = ttt.Player(ttt.Gaussian(0., 3.), 0.5, 0.0 )
    end
    h = ttt.History(events, results, times, priors)
    fp_a_1 = Vector{ttt.Gaussian}()
    bp_a_1 = Vector{ttt.Gaussian}()
    lh_a_1 = Vector{ttt.Gaussian}()
    wp_a_1 = Vector{ttt.Gaussian}()

    fp_a_2 = Vector{ttt.Gaussian}()
    bp_a_2 = Vector{ttt.Gaussian}()
    lh_a_2 = Vector{ttt.Gaussian}()
    wp_a_2 = Vector{ttt.Gaussian}()

    lh_b_1 = Vector{ttt.Gaussian}()
    lh_b_2 = Vector{ttt.Gaussian}()

    wp_b_1 = Vector{ttt.Gaussian}()
    wp_b_2 = Vector{ttt.Gaussian}()

    p_a_1 =  Vector{ttt.Gaussian}()
    p_a_2 =  Vector{ttt.Gaussian}()
    p_b_1 =  Vector{ttt.Gaussian}()
    p_b_2 =  Vector{ttt.Gaussian}()

    e_1 = Float64[]
    e_2 = Float64[]
    d_div_1 = Vector{ttt.Gaussian}()
    d_div_2 = Vector{ttt.Gaussian}()

    push!(fp_a_1, h.batches[1].skills["a"].forward)
    push!(bp_a_1, h.batches[1].skills["a"].backward)
    push!(lh_a_1, h.batches[1].skills["a"].likelihood)
    push!(wp_a_1, fp_a_1[end]*bp_a_1[end])
    push!(wp_b_1, h.batches[1].skills["b"].forward*h.batches[1].skills["b"].backward)
    push!(fp_a_2, h.batches[2].skills["a"].forward)
    push!(bp_a_2, h.batches[2].skills["a"].backward)
    push!(lh_a_2, h.batches[2].skills["a"].likelihood)
    push!(wp_a_2, fp_a_2[end]*bp_a_2[end])
    push!(wp_b_2, h.batches[2].skills["b"].forward*h.batches[2].skills["b"].backward)
    push!(e_1, h.batches[1].events[1].evidence)
    push!(e_2, h.batches[2].events[1].evidence)
    push!(lh_b_1, h.batches[1].skills["b"].likelihood)
    push!(lh_b_2, h.batches[2].skills["b"].likelihood)

    push!(p_a_1, wp_a_1[end]*lh_a_1[end])
    push!(p_a_2, wp_a_2[end]*lh_a_2[end])
    push!(p_b_1, wp_b_1[end]*lh_b_1[end])
    push!(p_b_2, wp_b_2[end]*lh_b_2[end])

    d_1 = wp_a_1[end]+Nbeta - wp_b_1[end]+Nbeta
    push!(d_div_1 , ttt.approx(d_1,0.,false)/d_1)
    d_2 = wp_a_1[end]+Nbeta - wp_b_1[end]+Nbeta
    push!(d_div_2 , ttt.approx(d_2,0.,false)/d_2)

    for _ in 1:10
        ttt.iteration(h)
        push!(fp_a_1, h.batches[1].skills["a"].forward)
        push!(bp_a_1, h.batches[1].skills["a"].backward)
        push!(lh_a_1, h.batches[1].skills["a"].likelihood)
        push!(wp_a_1, fp_a_1[end]*bp_a_1[end])
        push!(wp_b_1, h.batches[1].skills["b"].forward*h.batches[1].skills["b"].backward)
        push!(fp_a_2, h.batches[2].skills["a"].forward)
        push!(bp_a_2, h.batches[2].skills["a"].backward)
        push!(lh_a_2, h.batches[2].skills["a"].likelihood)
        push!(wp_a_2, fp_a_2[end]*bp_a_2[end])
        push!(wp_b_2, h.batches[2].skills["b"].forward*h.batches[2].skills["b"].backward)
        push!(e_1, h.batches[1].events[1].evidence)
        push!(e_2, h.batches[2].events[1].evidence)
        push!(lh_b_1, h.batches[1].skills["b"].likelihood)
        push!(lh_b_2, h.batches[2].skills["b"].likelihood)

        d_1 = wp_a_1[end]+Nbeta - wp_b_1[end]+Nbeta
        push!(d_div_1 , ttt.approx(d_1,0.,false)/d_1)
        d_2 = wp_a_1[end]+Nbeta - wp_b_1[end]+Nbeta
        push!(d_div_2 , ttt.approx(d_2,0.,false)/d_2)
        push!(p_a_1, wp_a_1[end]*lh_a_1[end])
        push!(p_a_2, wp_a_2[end]*lh_a_2[end])
        push!(p_b_1, wp_b_1[end]*lh_b_1[end])
        push!(p_b_2, wp_b_2[end]*lh_b_2[end])
    end

    df = DataFrame(fp_a_1_mu = [N.mu for N in fp_a_1]
                    ,fp_a_1_sigma = [N.sigma for N in fp_a_1]
                    ,bp_a_1_mu = [N.mu for N in bp_a_1]
                    ,bp_a_1_sigma = [N.sigma for N in bp_a_1]
                    ,lh_a_1_mu = [N.mu for N in lh_a_1]
                    ,lh_a_1_sigma = [N.sigma for N in lh_a_1]
                    ,wp_a_1_mu = [N.mu for N in wp_a_1]
                    ,wp_a_1_sigma = [N.sigma for N in wp_a_1]
                    ,p_a_1_mu = [N.mu for N in p_a_1]
                    ,p_a_1_sigma = [N.sigma for N in p_a_1]
                    ,p_a_2_mu = [N.mu for N in p_a_2]
                    ,p_a_2_sigma = [N.sigma for N in p_a_2]
                    ,p_b_1_mu = [N.mu for N in p_b_1]
                    ,p_b_1_sigma = [N.sigma for N in p_b_1]
                    ,p_b_2_mu = [N.mu for N in p_b_2]
                    ,p_b_2_sigma = [N.sigma for N in p_b_2]
                    ,wp_b_1_mu = [N.mu for N in wp_b_1]
                    ,wp_b_1_sigma = [N.sigma for N in wp_b_1]
                    ,fp_a_2_mu = [N.mu for N in fp_a_2]
                    ,fp_a_2_sigma = [N.sigma for N in fp_a_2]
                    ,bp_a_2_mu = [N.mu for N in bp_a_2]
                    ,bp_a_2_sigma = [N.sigma for N in bp_a_2]
                    ,lh_a_2_mu = [N.mu for N in lh_a_2]
                    ,lh_a_2_sigma = [N.sigma for N in lh_a_2]
                    ,wp_a_2_mu = [N.mu for N in wp_a_2]
                    ,wp_a_2_sigma = [N.sigma for N in wp_a_2]
                    ,wp_b_2_mu = [N.mu for N in wp_b_2]
                    ,wp_b_2_sigma = [N.sigma for N in wp_b_2]
                    ,lh_b_1_mu = [N.mu for N in lh_b_1]
                    ,lh_b_2_mu = [N.mu for N in lh_b_2]
                    ,lh_b_1_sigma = [N.sigma for N in lh_b_1]
                    ,lh_b_2_sigma = [N.sigma for N in lh_b_2]
                    ,e_1 = e_1
                    ,e_2 = e_2
                    ,d_div_1_mu = [N.mu for N in d_div_1]
                    ,d_div_1_sigma = [N.sigma for N in d_div_1]
                    ,d_div_2_mu = [N.mu for N in d_div_2]
                    ,d_div_2_sigma = [N.sigma for N in d_div_2]
                    )
    CSV.write("output/smoothing.csv", df; header=true)
    @test true
end
@testset "ttt_vs_ts" begin

    events = Array{Array{Array{String,1},1},1}()
    results =  Array{Array{Float64,1},1}()
    times = Int64[]
    ts_log_evidence = Float64[]
    ttt_log_evidence = Float64[]

    ts_last = Float64[]; ts_midle = Float64[]; ts_second = Float64[]
    ttt_last = Float64[]; ttt_midle = Float64[]; ttt_second = Float64[]
    for i in 0:127
        push!(events ,  [["a"],["b"]], [["a"],["b"]])
        push!(results ,  [1.,0.],[0.,1.] )
        push!(times,  i*2+1, i*2+2 )
        priors = Dict{String,ttt.Player}()
        for k in ["a", "b"]
            priors[k] = ttt.Player(ttt.Gaussian(0., 3.0), 0.5, 0.0)
        end
        h = ttt.History(events, results, times, priors)
        push!(ts_log_evidence , ttt.log_evidence(h))
        push!(ts_last, h.batches[end].events[1].evidence); push!(ts_midle, h.batches[i*2+1].events[1].evidence);
        push!(ts_second, h.batches[2].events[1].evidence);
        ttt.convergence(h)
        push!(ttt_last, h.batches[end].events[1].evidence); push!(ttt_midle, h.batches[i*2+1].events[1].evidence);
        push!(ttt_log_evidence, ttt.log_evidence(h))
        push!(ttt_second, h.batches[2].events[1].evidence)
    end
    df = DataFrame(ts_log_evidence = ts_log_evidence
                    ,ttt_log_evidence = ttt_log_evidence
                    ,ts_last = ts_last
                    ,ts_midle = ts_midle
                    ,ts_second = ts_second
                    ,ttt_second = ttt_second
                    ,ttt_last = ttt_last
                    ,ttt_midle = ttt_midle )

    CSV.write("output/ttt_vs_ts.csv", df; header=true)
    @test true
end
@testset "Best gamma" begin

    function skill(exp::Int64, alpha::Float64=0.133)
        return exp^alpha - 1
    end
    mean_agent = [skill(i) for i in 1:1000]
    beta = 0.5

    using Random
    Random.seed!(1)

    mean_target = [(Random.randn(1)[1]*beta + skill(i)) for i in 1:1000]
    perf_target = [(Random.randn(1)[1]*beta + mean_target[i]) for i in 1:1000]
    perf_agent = [(Random.randn(1)[1]*beta + mean_agent[i]) for i in 1:1000]
    events = [ [["a"], [string(i)] ] for i in 1:1000]
    results = [ perf_agent[i] > perf_target[i] ? [1.,0.] : [0.,1.] for i in 1:1000 ]
    batches= [i for i in 1:1000 ]

    selected_gammas = [0.005,0.01,0.015,0.02,0.025]

    gammas = [gamma for gamma in 0.001:0.001:0.04]
    evidencias_ts = Float64[]
    evidencias_ttt = Float64[]
    learning_curves = Vector{Vector{Float64}}()
    for gamma in gammas#gamma=0.015

        priors = Dict{String,ttt.Player}()
        priors["a"] = ttt.Player(ttt.Gaussian(0., 3.0), beta, gamma)
        for k in 1:1000
            priors[string(k)] = ttt.Player(ttt.Gaussian(mean_target[k], 0.5), beta, 0.0)
        end
        h = ttt.History(events, results, batches, priors)
        push!(evidencias_ts, ttt.log_evidence(h))
        ttt.convergence(h)
        push!(evidencias_ttt, ttt.log_evidence(h))
        if gamma in selected_gammas
            push!(learning_curves, [r.mu for (t, r) in ttt.learning_curves(h)["a"]] )
        end
    end
    #using Plots
    #plot([g for g in gammas ],evidencias_ts)
    #plot!([g for g in gammas ] ,evidencias_ttt)

    #plot(mean_agent )
    #plot!(learning_curves[3])

    @test 0.015 == gammas[argmax(evidencias_ttt)]

    df = DataFrame(gammas = gammas
                    ,evidencias_ts = evidencias_ts
                    ,evidencias_ttt = evidencias_ttt)
    CSV.write("output/best_gamma-evidences.csv", df; header=true)

    df = DataFrame(mean_agent = mean_agent
                    ,lc_005 = learning_curves[1]
                    ,lc_01  = learning_curves[2]
                    ,lc_015 = learning_curves[3]
                    ,lc_020 = learning_curves[4]
                    ,lc_025 = learning_curves[5])
    CSV.write("output/best_gamma-learning_curves.csv", df; header=true)

    #plot(mean_agent)
    #plot!(learning_curves[3])
end
@testset "Chain of betas" begin
    using Random
    Random.seed!(1)
    beta = 1.0
    agents = [string(i) for i in -5:5]
    true_skills = [i for i in -5.0:5.0]
    mu_prior = [ abs(i)==5.0 ? i : 0.0 for i in -5.0:5.0]
    sigma_prior = [ abs(i)==5.0 ? 1e-7 : 6.0 for i in -5.0:5.0]
    priors = Dict{String,ttt.Player}()
    for i in -5:5
        priors[string(i)] = ttt.Player(ttt.Gaussian(mu_prior[i+6], sigma_prior[i+6]), beta, 0.0)
    end

    events = [ [[string(a-1)],[string(a)]]  for e in 1:100 for a in -4:5]
    results = [ (Random.randn(1)[1]*beta + a-1) > (Random.randn(1)[1]*beta + a) ? [0.,1.] : [1.,0.] for i in 1:100 for a in -4.0:5.0]
    times = [ e  for e in 1:100 for a in -4:5]
    h = ttt.History(events , results, times, priors, iter=100)
    ttt.convergence(h)
    diffs = [(ttt.posterior(h.batches[1],string(a))-ttt.posterior(h.batches[1],string(a-1))).mu for a in -3:4]
    probs = [ 1-ttt.cdf(ttt.Gaussian(d,sqrt(2)),0.0)  for d in diffs]
    for p in probs
        @test abs(1-ttt.cdf(ttt.Gaussian(1.0,sqrt(2)),0.0)-p) < 0.04
    end
end
@testset "Online predictions" begin
    priors = Dict{String,ttt.Player}()
    priors["a"] = ttt.Player(ttt.Gaussian(0.0, 3.0), 1.0, 0.0)
    priors["b"] = ttt.Player(ttt.Gaussian(0.0, 3.0), 1.0, 0.0)
    priors["c"] = ttt.Player(ttt.Gaussian(2.0, 0.5), 1.0, 0.0)
#
    composition = [[["a"],["b"]], [["a"],["c"]], [["a"],["b"]]]
    results = [ [0.,1.], [1.,0.], [1.,0.]]
    times = [1,2,3]
    h = ttt.History(composition , results, times, priors, iter=100, gamma=0.0,online=true)

    h.batches[3].skills["b"].forward
    h.batches[3].skills["b"].online
end
end