text:
module semiclassical

using QuantumOpticsBase
import Base: ==
import ..timeevolution: integrate, recast!, jump, integrate_mcwf, jump_callback, as_vector, QO_CHECKS
import LinearAlgebra: normalize, normalize!

using Random, LinearAlgebra
import OrdinaryDiffEq

# TODO: Remove imports
import DiffEqCallbacks, RecursiveArrayTools.copyat_or_push!
Base.@pure pure_inference(fout,T) = Core.Compiler.return_type(fout, T)

using ..timeevolution

const QuantumState{B} = Union{Ket{B}, Operator{B,B}}

"""
Semi-classical state.

It consists of a quantum part, which is either a `Ket` or a `DenseOperator` and
a classical part that is specified as a complex vector of arbitrary length.
"""
mutable struct State{B,T,C}
    quantum::T
    classical::C
    function State(quantum::T, classical::C) where {B,T<:QuantumState{B},C}
        new{B,T,C}(quantum, classical)
    end
end

Base.length(state::State) = length(state.quantum) + length(state.classical)
Base.copy(state::State) = State(copy(state.quantum), copy(state.classical))
Base.eltype(state::State) = promote_type(eltype(state.quantum),eltype(state.classical))
normalize!(state::State) where {B,T} = (normalize!(state.quantum); state)
normalize(state::State) = State(normalize(state.quantum),copy(state.classical))

function ==(a::State, b::State)
    QuantumOpticsBase.samebases(a.quantum, b.quantum) &&
    length(a.classical)==length(b.classical) &&
    (a.classical==b.classical) &&
    (a.quantum==b.quantum)
end

QuantumOpticsBase.expect(op, state::State) = expect(op, state.quantum)
QuantumOpticsBase.variance(op, state::State) = variance(op, state.quantum)
QuantumOpticsBase.ptrace(state::State, indices) = State(ptrace(state.quantum, indices), state.classical)

QuantumOpticsBase.dm(x::State) = State(dm(x.quantum), x.classical)

"""
    semiclassical.schroedinger_dynamic(tspan, state0, fquantum, fclassical[; fout, ...])

Integrate time-dependent Schrödinger equation coupled to a classical system.

# Arguments
* `tspan`: Vector specifying the points of time for which the output should
        be displayed.
* `psi0`: Initial semi-classical state [`semiclassical.State`](@ref).
* `fquantum`: Function `f(t, psi, u) -> H` returning the time and or state
        dependent Hamiltonian.
* `fclassical!`: Function `f!(du, u, psi, t)` calculating the possibly time and
        state dependent derivative of the classical equations and storing it
        in the vector `du`.
* `fout=nothing`: If given, this function `fout(t, state)` is called every time
        an output should be displayed. ATTENTION: The given state is neither
        normalized nor permanent!
* `kwargs...`: Further arguments are passed on to the ode solver.
"""
function schroedinger_dynamic(tspan, state0::State, fquantum, fclassical!;
                fout=nothing,
                kwargs...)
    dschroedinger_(t, state, dstate) = dschroedinger_dynamic!(dstate, fquantum, fclassical!, state, t)
    x0 = Vector{eltype(state0)}(undef, length(state0))
    recast!(x0,state0)
    state = copy(state0)
    dstate = copy(state0)
    integrate(tspan, dschroedinger_, x0, state, dstate, fout; kwargs...)
end

"""
    semiclassical.master_dynamic(tspan, state0, fquantum, fclassical!; <keyword arguments>)

Integrate time-dependent master equation coupled to a classical system.

# Arguments
* `tspan`: Vector specifying the points of time for which output should
        be displayed.
* `rho0`: Initial semi-classical state [`semiclassical.State`](@ref).
* `fquantum`: Function `f(t, rho, u) -> (H, J, Jdagger)` returning the time
        and/or state dependent Hamiltonian and Jump operators.
* `fclassical!`: Function `f!(du, u, rho, t)` calculating the possibly time and
        state dependent derivative of the classical equations and storing it
        in the complex vector `du`.
* `fout=nothing`: If given, this function `fout(t, state)` is called every time
        an output should be displayed. ATTENTION: The given state is not
        permanent!
* `kwargs...`: Further arguments are passed on to the ode solver.
"""
function master_dynamic(tspan, state0::State{B,T}, fquantum, fclassical!;
                rates=nothing,
                fout=nothing,
                tmp=copy(state0.quantum),
                kwargs...) where {B,T<:Operator}
    dmaster_(t, state, dstate) = dmaster_h_dynamic!(dstate, fquantum, fclassical!, rates, state, tmp, t)
    x0 = Vector{eltype(state0)}(undef, length(state0))
    recast!(x0,state0)
    state = copy(state0)
    dstate = copy(state0)
    integrate(tspan, dmaster_, x0, state, dstate, fout; kwargs...)
end

function master_dynamic(tspan, state0::State{B,T}, fquantum, fclassical!; kwargs...) where {B,T<:Ket{B}}
    master_dynamic(tspan, dm(state0), fquantum, fclassical!; kwargs...)
end

"""
    semiclassical.mcwf_dynamic(tspan, psi0, fquantum, fclassical!, fjump_classical!; <keyword arguments>)

Calculate MCWF trajectories coupled to a classical system. **NOTE**: The quantum
state with which `fquantum` and `fclassical!` are called is NOT NORMALIZED. Make
sure to take this into account when computing expectation values!

# Arguments
* `tspan`: Vector specifying the points of time for which output should
        be displayed.
* `psi0`: Initial semi-classical state [`semiclassical.State`](@ref) featuring
        a `Ket`(@ref).
* `fquantum`: Function `f(t, psi, u) -> (H, J, Jdagger)` returning the time
        and/or state dependent Hamiltonian and Jump operators.
* `fclassical!`: Function `f!(du, u, psi, t)` calculating the possibly time and
        state dependent derivative of the classical equations and storing it
        in the complex vector `du`.
* `fjump_classical!`: Function `f(t, psi, u, i)` performing a classical jump when a
        quantum jump of the i-th jump operator occurs.
* `fout=nothing`: If given, this function `fout(t, state)` is called every time
        an output should be displayed. ATTENTION: The given state is not
        permanent!
* `display_beforeevent`: Choose whether or not an additional point should be saved
        before a jump occurs. Default is false.
* `display_afterevent`: Choose whether or not an additional point should be saved
        after a jump occurs. Default is false.
* `display_jumps=false`: If set to true, an additional list of times and indices
        is returned. These correspond to the times at which a jump occured and
        the index of the jump operators with which the jump occured, respectively.
* `kwargs...`: Further arguments are passed on to the ode solver.
"""
function mcwf_dynamic(tspan, psi0::State{B,T}, fquantum, fclassical!, fjump_classical!;
                seed=rand(UInt),
                rates=nothing,
                fout=nothing,
                kwargs...) where {B,T<:Ket}
    tmp=copy(psi0.quantum)
    dmcwf_(t, psi, dpsi) = dmcwf_h_dynamic!(dpsi, fquantum, fclassical!, rates, psi, tmp, t)
    j_(rng, t, psi, psi_new) = jump_dynamic(rng, t, psi, fquantum, fclassical!, fjump_classical!, psi_new, rates)
    x0 = Vector{eltype(psi0)}(undef, length(psi0))
    recast!(x0,psi0)
    psi = copy(psi0)
    dpsi = copy(psi0)
    integrate_mcwf(dmcwf_, j_, tspan, psi, seed, fout; kwargs...)
end

function recast!(x::C,state::State{B,T,C}) where {B,T,C}
    N = length(state.quantum)
    copyto!(x, 1, state.quantum.data, 1, N)
    copyto!(x, N+1, state.classical, 1, length(state.classical))
    x
end

function recast!(state::State{B,T,C},x::C) where {B,T,C}
    N = length(state.quantum)
    copyto!(state.quantum.data, 1, x, 1, N)
    copyto!(state.classical, 1, x, N+1, length(state.classical))
end

"""
    dschroedinger_dynamic!(dstate, fquantum, fclassical!, state)

Update the semiclassical state `dstate` according to a time-dependent,
semiclassical Schrödinger equation.

See also: [`semiclassical.schroedinger_dynamic`](@ref)
"""
function dschroedinger_dynamic!(dstate, fquantum, fclassical!, state, t)
    fquantum_(t, psi) = fquantum(t, state.quantum, state.classical)
    timeevolution.dschroedinger_dynamic!(dstate.quantum, fquantum_, state.quantum, t)
    fclassical!(dstate.classical, state.classical, state.quantum, t)
    return dstate
end

"""
    dmaster_h_dynamic!(dstate, fquantum, fclassical!, rates, state, tmp, t)

Update the semiclassical state `dstate` according to a time-dependent,
semiclassical master eqaution.

See also: [`semiclassical.master_dynamic`](@ref)
"""
function dmaster_h_dynamic!(dstate, fquantum, fclassical!, rates, state, tmp, t)
    fquantum_(t, rho) = fquantum(t, state.quantum, state.classical)
    timeevolution.dmaster_h_dynamic!(dstate.quantum, fquantum_, rates, state.quantum, tmp, t)
    fclassical!(dstate.classical, state.classical, state.quantum, t)
    return dstate
end

"""
    dmcwf_h_dynamic!(dpsi, fquantum, fclassical!, rates, psi, tmp, t)

Update the semiclassical state `dpsi` according to a time-dependent, semiclassical
and non-Hermitian Schrödinger equation (MCWF).

See also: [`semiclassical.mcwf_dynamic`](@ref)
"""
function dmcwf_h_dynamic!(dpsi, fquantum, fclassical!, rates, psi, tmp, t)
    fquantum_(t, rho) = fquantum(t, psi.quantum, psi.classical)
    timeevolution.dmcwf_h_dynamic!(dpsi.quantum, fquantum_, rates, psi.quantum, tmp, t)
    fclassical!(dpsi.classical, psi.classical, psi.quantum, t)
    return dpsi
end

function jump_dynamic(rng, t, psi, fquantum, fclassical!, fjump_classical!, psi_new, rates)
    result = fquantum(t, psi.quantum, psi.classical)
    QO_CHECKS[] && @assert 3 <= length(result) <= 4
    J = result[2]
    if length(result) == 3
        rates_ = rates
    else
        rates_ = result[4]
    end
    i = jump(rng, t, psi.quantum, J, psi_new.quantum, rates_)
    fjump_classical!(psi.classical, psi_new.quantum, i, t)
    psi_new.classical .= psi.classical
    return i
end

function jump_callback(jumpfun, seed, scb, save_before!,
                        save_after!, save_t_index, psi0::State)
    tmp = copy(psi0)
    psi_tmp = copy(psi0)
    rng = MersenneTwister(convert(UInt, seed))
    jumpnorm = Ref(rand(rng))
    n = length(psi0.quantum)
    djumpnorm(x::Vector, t, integrator) = norm(x[1:n])^2 - (1-jumpnorm[])

    function dojump(integrator)
        x = integrator.u
        t = integrator.t

        affect! = scb.affect!
        save_before!(affect!,integrator)
        recast!(psi_tmp,x)
        i = jumpfun(rng, t, psi_tmp, tmp)
        recast!(x,tmp)
        save_after!(affect!,integrator)
        save_t_index(t,i)

        jumpnorm[] = rand(rng)
        return nothing
    end

    return OrdinaryDiffEq.ContinuousCallback(djumpnorm,dojump,
                     save_positions = (false,false))
end
as_vector(psi::State) = vcat(psi.quantum.data[:], psi.classical)

end # module


text:
# draw histogram of lidar's sensing data grouped by hour

module HistGrpByHour
    using Plots, DataFrames, CSV, Statistics
    pyplot()

    function main(is_test=false)
        data_path = joinpath(split(@__FILE__, "src")[1], "data/sensor_data_600.txt")
        df_org = CSV.read(data_path, DataFrame,
                      header=["date", "time", "ir", "lidar"],
                      delim=' ')

        hours = [Int64(floor(e/10000)) for e in df_org.time]

        # add hour array to df
        df_h = DataFrame(hour=hours)
        df_new = hcat(df_org, df_h)

        # sort
        df_sort = sort(df_new, "hour")

        # grouping
        df_grp = groupby(df_sort, "hour")

        # histogram
        df_grp_6 = df_grp[6]
        df_grp_14 = df_grp[14]
        histogram(df_grp_6.lidar, color=:blue, label="6")
        histogram!(df_grp_14.lidar, color=:orange, label="14")

        if is_test == false
            save_path = joinpath(split(@__FILE__, "src")[1], "src/prob_stats/complex_dist/hist_grp_by_hour/hist_grp_by_hour.png")
            savefig(save_path)
        end
    end
end

ext: jl

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
#
# Tests for the moular K-Profile-Parameterization module
#

function test_default_model_init(N=4, L=4.3)
    model = ModularKPP.Model(N=N, L=L)
    return model.grid.N == N && model.grid.L == L
end

function time_step_model(diffusivity, nonlocalflux, mixingdepth, kprofile)

    model = ModularKPP.Model(N=4, L=3,
                             diffusivity = diffusivity,
                            nonlocalflux = nonlocalflux,
                             mixingdepth = mixingdepth,
                                kprofile = kprofile)

    iterate!(model, 1e-16)
    return true
end

diffusivity_models = (
    ModularKPP.LMDDiffusivity(),
    ModularKPP.HoltslagDiffusivity()
)

nonlocal_flux_models = (
    ModularKPP.LMDCounterGradientFlux(),
    ModularKPP.DiagnosticPlumeModel()
)

mixing_depth_models = (
    ModularKPP.LMDMixingDepth(),
    ModularKPP.ROMSMixingDepth()
)

shape_functions = (
    ModularKPP.StandardCubicPolynomial(),
    ModularKPP.GeneralizedCubicPolynomial()
)

@testset "Modular KPP" begin
    @test test_default_model_init()

    # Modularity
    for diffusivity in diffusivity_models
        for nonlocalflux in nonlocal_flux_models
            for mixingdepth in mixing_depth_models
                for kprofile in shape_functions
                    @test time_step_model(diffusivity, nonlocalflux, mixingdepth, kprofile)
                end
            end
        end
    end
end


text:
using JuliaWorld
using Test

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


text:
########## Fish ##########

mutable struct Fish
    cycle::UInt8
end

cycle(f::Fish) = f.cycle

function next_day!(fish::Fish)
    (spawn = cycle(fish) == 0x00) ? fish.cycle = 0x06 : fish.cycle -= 0x01
    return spawn
end

function next_day!(fishes::Vector{Fish})
    n_spawns = sum(next_day!.(fishes))
    append!(fishes, Fish.(0x08 * ones(UInt8, n_spawns)))
    return n_spawns
end

########## Script ##########

function main()

    # open file and store content into an input string
    input = open("inputs//day6.txt", "r") do io
        read(io, String)
    end

    # split input into a data vector
    data = String.(split(input, ","))

    # parse data as a vector of fishes
    fishes = Fish.(parse.(UInt8, data))

    # part 1
    count::BigInt = length(fishes)
    for _ = 1:80
        count += next_day!(fishes)
    end
    println("part 1: ", count)

    # part 2
    for _ = 81:256
        count += next_day!(fishes)
    end
    println("part 2: ", count)

end

main()


text:
# Functions for STaylor1
@generated function exp(a::STaylor1{N,T}) where {N, T <: Number}
    ex_calc = quote end
    append!(ex_calc.args, Any[nothing for i in 1:N])
    syms = Symbol[Symbol("c$i") for i in 1:N]

    sym = syms[1]
    ex_line = :($(syms[1]) = exp(a[0]))
    ex_calc.args[1] = ex_line

    for k in 1:(N-1)
        kT = convert(T,k)
        sym = syms[k+1]
        ex_line = :($kT * a[$k] * $(syms[1]))
        @inbounds for i = 1:k-1
            ex_line = :($ex_line + $(kT-i) * a[$(k-i)] * $(syms[i+1]))
        end
        ex_line = :(($ex_line)/$kT)
        ex_line = :($sym = $ex_line)
        ex_calc.args[k+1] = ex_line
    end

    exout = :(($(syms[1]),))
    for i = 2:N
        push!(exout.args, syms[i])
    end
    return quote
               Base.@_inline_meta
               $ex_calc
               return STaylor1{N,T}($exout)
            end
end

@generated function log(a::STaylor1{N,T}) where {N, T <: Number}
    ex_calc = quote end
    append!(ex_calc.args, Any[nothing for i in 1:N])
    syms = Symbol[Symbol("c$i") for i in 1:N]

    (N >= 1) && (ex_calc.args[1] = :($(syms[1]) = log(constant_term(a))))
    (N >= 2) && (ex_calc.args[2] = :($(syms[2]) = a[1]/constant_term(a)))

    for k in 2:(N-1)
        ex_line = :($(k-1)*a[1]*$(syms[k]))
        @inbounds for i = 2:k-1
            ex_line = :($ex_line + $(k-i)*a[$i] * $(syms[k+1-i]))
        end
        ex_line = :((a[$k] - ($ex_line)/$(convert(T,k)))/constant_term(a))
        ex_line = :($(syms[k+1]) = $ex_line)
        ex_calc.args[k+1] = ex_line
    end

    exout = :(($(syms[1]),))
    for i = 2:N
        push!(exout.args, syms[i])
    end

    return quote
               Base.@_inline_meta
               iszero(constant_term(a)) && throw(ArgumentError("""
                       The 0-th order `STaylor1` coefficient must be non-zero
                       in order to expand `log` around 0.
                       """))
               $ex_calc
               return STaylor1{N,T}($exout)
            end
end

sin(a::STaylor1{N,T}) where {N, T <: Number} = sincos(a)[1]
cos(a::STaylor1{N,T}) where {N, T <: Number} = sincos(a)[2]
@generated function sincos(a::STaylor1{N,T}) where {N, T <: Number}

    ex_calc = quote end
    append!(ex_calc.args, Any[nothing for i in 1:(2*N)])

    syms_s = Symbol[Symbol("c$i") for i in 1:N]
    syms_c = Symbol[Symbol("c2$i") for i in 1:N]

    ex_line_s = :($(syms_s[1]) = sin(a[0]))
    ex_line_c = :($(syms_c[1]) = cos(a[0]))
    ex_calc.args[1] = ex_line_s
    ex_calc.args[2] = ex_line_c

    for k = 1:(N - 1)
        ex_line_s = :(a[1]*$(syms_c[k]))
        ex_line_c = :(-a[1]*$(syms_s[k]))
        for i = 2:k
            ex_line_s = :($ex_line_s + $i*a[$i]*$(syms_c[(k - i + 1)]))
            ex_line_c = :($ex_line_c - $i*a[$i]*$(syms_s[(k - i + 1)]))
        end
        ex_line_s = :($(syms_s[k + 1]) = ($ex_line_s)/$k)
        ex_line_c = :($(syms_c[k + 1]) = ($ex_line_c)/$k)
        ex_calc.args[2*k + 1] = ex_line_s
        ex_calc.args[2*k + 2] = ex_line_c
    end

    exout_s = :(($(syms_s[1]),))
    for i = 2:N
        push!(exout_s.args, syms_s[i])
    end

    exout_c = :(($(syms_c[1]),))
    for i = 2:N
        push!(exout_c.args, syms_c[i])
    end

    return quote
               Base.@_inline_meta
               $ex_calc
               return STaylor1{N,T}($exout_s), STaylor1{N,T}($exout_c)
            end
end

# Functions for STaylor1
#=
@generated function tan(a::STaylor1{N,T}) where {N, T <: Number}

    ex_calc = quote end
    append!(ex_calc.args, Any[nothing for i in 1:(4*N)])
    syms = Symbol[Symbol("c$i") for i in 1:N]
    syms2 = Symbol[Symbol("c2$i") for i in 1:N]

    for i = 1:N
        ex_calc.args[i] = :($(syms2[i]) = 0.0)
    end

    ex_line_c = :($(syms[1]) = tan(a[0]))
    ex_line_c2 = :($(syms2[1]) = ($(syms[1]))^2)
    ex_calc.args[N + 1] = ex_line_c
    ex_calc.args[N + 2] = ex_line_c2

    for k = 1:(N - 1)

        kodd = k%2
        kend = div(k - 2 + kodd, 2)
        kdiv2 = div(k, 2)

        ex_line_c = :($(k-1)*a[$(k-1)]*$(syms2[2]))
        for i = 1:(k - 1)
            q = k - i
            ex_line_c = :($ex_line_c + ($q)*a[$q]*$(syms2[i + 1]))
        end
        ex_line_c = :(a[$k] + ($ex_line_c)/$k)
        ex_line_c = :($(syms[k + 1]) = $ex_line_c)
        ex_calc.args[(3*k-2) + N + 2] = ex_line_c

        ex_line_c2 = :(a[0]*a[$k])
        for i = 1:kend
            ex_line_c2 = :($ex_line_c2 + a[$i]*a[$k - $i])
        end
        ex_line_c2 = :(2*$ex_line_c2)

        if kodd != 1
            ex_line_c2 = :($ex_line_c2 + a[$kdiv2]^2)
        end

        ex_line_c2 = :($(syms2[k + 1]) = $ex_line_c2)
        ex_calc.args[(3*k-1) + 2 + N] = ex_line_c2
        blar = k + 1
        blar_str = "c$(blar) "
        ex_calc.args[3*k + N + 2] = :(println($blar_str*string($(syms[k + 1]))))
    end

    exout = :(($(syms[1]),))
    for i = 2:N
        push!(exout.args, syms[i])
    end

    return quote
               Base.@_inline_meta
               $ex_calc
               return STaylor1{N,T}($exout)
            end
end
=#

# Functions for STaylor1
#=
@generated function asin(a::STaylor1{N,T}) where {N, T <: Number}
    ex_calc = quote end
    append!(ex_calc.args, Any[nothing for i in 1:N])
    syms = Symbol[Symbol("c$i") for i in 1:N]

    sym = syms[1]
    ex_line = :($(sym) = asin(a[0]))
    ex_calc.args[1] = ex_line

    #=
    for k in 1:(N-1)
        kT = convert(T,k)
        sym = syms[k+1]
        ex_line = :($kT * a[$k] * $(syms[1]))
        @inbounds for i = 1:k-1
            ex_line = :($ex_line + $(kT-i) * a[$(k-i)] * $(syms[i+1]))
        end
        ex_line = :(($ex_line)/$kT)
        ex_line = :($sym = $ex_line)
        ex_calc.args[k+1] = ex_line
    end
    =#

    exout = :(($(syms[1]),))
    for i = 2:N
        push!(exout.args, syms[i])
    end
    return quote
               Base.@_inline_meta
               $ex_calc
               return STaylor1{N,T}($exout)
            end
end
=#

# Functions for STaylor1
#=
@generated function acos(a::STaylor1{N,T}) where {N, T <: Number}
    ex_calc = quote end
    append!(ex_calc.args, Any[nothing for i in 1:N])
    syms = Symbol[Symbol("c$i") for i in 1:N]

    sym = syms[1]
    ex_line = :($(sym) = acos(a[0]))
    ex_calc.args[1] = ex_line

    #=
    for k in 1:(N-1)
        kT = convert(T,k)
        sym = syms[k+1]
        ex_line = :($kT * a[$k] * $(syms[1]))
        @inbounds for i = 1:k-1
            ex_line = :($ex_line + $(kT-i) * a[$(k-i)] * $(syms[i+1]))
        end
        ex_line = :(($ex_line)/$kT)
        ex_line = :($sym = $ex_line)
        ex_calc.args[k+1] = ex_line
    end
    =#

    exout = :(($(syms[1]),))
    for i = 2:N
        push!(exout.args, syms[i])
    end
    return quote
               Base.@_inline_meta
               $ex_calc
               return STaylor1{N,T}($exout)
            end
end
=#

# Functions for STaylor1
#=
@generated function atan(a::STaylor1{N,T}) where {N, T <: Number}
    ex_calc = quote end
    append!(ex_calc.args, Any[nothing for i in 1:N])
    syms = Symbol[Symbol("c$i") for i in 1:N]

    sym = syms[1]
    ex_line = :($(sym) = atan(a[0]))
    ex_calc.args[1] = ex_line

    #=
    for k in 1:(N-1)
        kT = convert(T,k)
        sym = syms[k+1]
        ex_line = :($kT * a[$k] * $(syms[1]))
        @inbounds for i = 1:k-1
            ex_line = :($ex_line + $(kT-i) * a[$(k-i)] * $(syms[i+1]))
        end
        ex_line = :(($ex_line)/$kT)
        ex_line = :($sym = $ex_line)
        ex_calc.args[k+1] = ex_line
    end
    =#

    exout = :(($(syms[1]),))
    for i = 2:N
        push!(exout.args, syms[i])
    end
    return quote
               Base.@_inline_meta
               $ex_calc
               return STaylor1{N,T}($exout)
            end
end
=#

sinh(a::STaylor1{N,T}) where {N, T <: Number} = sinhcosh(a)[1]
cosh(a::STaylor1{N,T}) where {N, T <: Number} = sinhcosh(a)[2]
@generated function sinhcosh(a::STaylor1{N,T}) where {N, T <: Number}

    ex_calc = quote end
    append!(ex_calc.args, Any[nothing for i in 1:(2*N)])

    syms_s = Symbol[Symbol("c$i") for i in 1:N]
    syms_c = Symbol[Symbol("c2$i") for i in 1:N]

    ex_line_s = :($(syms_s[1]) = sinh(a[0]))
    ex_line_c = :($(syms_c[1]) = cosh(a[0]))
    ex_calc.args[1] = ex_line_s
    ex_calc.args[2] = ex_line_c

    for k = 1:(N - 1)
        ex_line_s = :(a[1]*$(syms_c[k]))
        ex_line_c = :(a[1]*$(syms_s[k]))
        for i = 2:k
            ex_line_s = :($ex_line_s + $i*a[$i]*$(syms_c[(k - i + 1)]))
            ex_line_c = :($ex_line_c + $i*a[$i]*$(syms_s[(k - i + 1)]))
        end
        ex_line_s = :($(syms_s[k + 1]) = ($ex_line_s)/$k)
        ex_line_c = :($(syms_c[k + 1]) = ($ex_line_c)/$k)
        ex_calc.args[2*k + 1] = ex_line_s
        ex_calc.args[2*k + 2] = ex_line_c
    end

    exout_s = :(($(syms_s[1]),))
    for i = 2:N
        push!(exout_s.args, syms_s[i])
    end

    exout_c = :(($(syms_c[1]),))
    for i = 2:N
        push!(exout_c.args, syms_c[i])
    end

    return quote
               Base.@_inline_meta
               $ex_calc
               return STaylor1{N,T}($exout_s), STaylor1{N,T}($exout_c)
            end
end

# Functions for STaylor1
@generated function tanh(a::STaylor1{N,T}) where {N, T <: Number}

    ex_calc = quote end
    append!(ex_calc.args, Any[nothing for i in 1:(2*N)])
    syms = Symbol[Symbol("c$i") for i in 1:N]
    syms2 = Symbol[Symbol("c2$i") for i in 1:N]

    ex_line_c = :($(syms[1]) = tanh(a[0]))
    ex_line_c2 = :($(syms2[1]) = ($(syms[1]))^2)
    ex_calc.args[1] = ex_line_c
    ex_calc.args[2] = ex_line_c2

    for k = 1:(N - 1)

        kodd = k%2
        kend = div(k - 2 + kodd, 2)
        kdiv2 = div(k, 2)

        ex_line_c = :($k*a[$k]*$(syms2[1]))
        for i = 1:(k - 1)
            q = k - i
            ex_line_c = :($ex_line_c + ($q)*a[$q]*$(syms2[i+1]))
        end
        ex_line_c = :(a[$k] - $ex_line_c/$k)
        ex_line_c = :($(syms[k + 1]) = $ex_line_c)
        ex_calc.args[2*k + 1] = ex_line_c

        ex_line_c2 = :(a[1]*a[$k])
        @inbounds for i = 1:kend
            ex_line_c2 = :($ex_line_c2 + a[$i + 1]*a[$k - $i])
        end
        ex_line_c2 = :(2*$ex_line_c2)

        if kodd != 1
            ex_line_c2 = :($ex_line_c2 + a[$kdiv2 + 1]^2)
        end

        ex_line_c2 = :($(syms2[k + 1]) = $ex_line_c2)
        ex_calc.args[2*k + 2] = ex_line_c2
    end

    exout = :(($(syms[1]),))
    for i = 2:N
        push!(exout.args, syms[i])
    end

    return quote
               Base.@_inline_meta
               $ex_calc
               return STaylor1{N,T}($exout)
            end
end


text:
using MatrixProfile
# Download data from https://sites.google.com/site/snippetfinderinfo/home/MixedBag.zip
datapath = "/home/fredrikb/Downloads/MixedBag/"
cd(datapath)

## Some data
T = parse.(Float64, split(join(Char.(read("01911m_02019m_III_7680_200.txt"))), ','))
snips = snippets(T, 3, 200, m=50)
plot(snips)

## Power demand
T = parse.(Float64, readlines("Powerdemand_12_4500_200.txt"))
snips = snippets(T, 4, 24, m=24)
plot(snips)

profile = matrix_profile(T, 24)
mot = motifs(profile, 4, r=3, th=48)
plot(profile, mot, legend=false)
plot(mot)

## Walking styles
T = parse.(Float64, readlines("PAMAP_Subject4_NormalWalking_NordicWalking_17001_500.txt"))
T = T[1:2:end]
snips = snippets(T, 4, 200)
plot(snips)

profile = matrix_profile(T, 50)
mot = motifs(profile, 4, r=3)
plot(mot)
plot(profile, mot, legend=false)

## Robot dog
T = parse.(Float64, readlines("RoboticDogActivityY_64_4000_400.txt"))
snips = snippets(T, 4, 100)
plot(snips)

profile = matrix_profile(T, 50)
mot = motifs(profile, 4, r=2, th=200)
plot(profile, mot, legend=false)
plot(mot)

using DynamicAxisWarping

dist = DTW(3)
normalizer = ZNormalizer
profile2 = matrix_profile(T, 50, dist, normalizer=normalizer)

using Distances
struct ZWrapper{D} <: Distances.Metric
    dist::D
end
function Distances.evaluate(d::ZWrapper, x, y)
    evaluate(d.dist, znorm(x), znorm(y))
end
(d::ZWrapper)(x,y) = evaluate(d,x,y)

mot2 = motifs(profile2, 4, r=3, th=200, dist=ZWrapper(dist))
plot(profile2, mot2, legend=false)
plot(mot2)

snips = snippets(T, 4, 100, dist)
plot(snips)

@btime snippets($(T[1:2000]), 4, 100)
@profiler snippets((T[1:2000]), 4, 100)


text:
@testset "HTTP requests to mort.SOA.org" begin

    tbl = get_SOA_table(60029)

    @test tbl isa MortalityTable
    @test tbl[0] == 0.10139

    get_SOA_table(60029)

    @test MortalityTables.table("Australian Life Tables 1891-1900 Female") isa MortalityTable
    @test MortalityTables.table(60029) isa MortalityTable
    @test MortalityTables.table("Australian Life Tables 1891-1900 Female")[0] == 0.10139

    # issue #69
    tbl = get_SOA_table(887)
    @test tbl isa MortalityTable
    tbl = get_SOA_table(2585)
    @test tbl isa MortalityTable
end


text:
###
#   Set operations
###

# Envolpe/hull. Keeps it as a IntervalUnion
env(x :: IntervalU) = IntervalU([hull(x.v)])
hull(x :: IntervalU) = IntervalU([hull(x.v)])

# Computes the complement of a IntervalUnion
function complement(x :: IntervalU)

    v = sort(x.v)
    vLo = left.(v)
    vHi = right.(v)

    vLo[1] == -∞ ? popfirst!(vLo) : pushfirst!(vHi, -∞)
    vHi[end] == ∞ ? pop!(vHi) : push!(vLo, ∞)

    complements = interval.(vHi,vLo)

    return intervalU(complements)
end

complement(x :: Interval) = complement(intervalU(x))

# bisect x at α
function bisect( x :: IntervalU, α = 0)

    v = deepcopy(x.v)

    αIn = α .∈ v
    these = findall(αIn .== 1)

    if isempty(these); return x ; end

    this = popat!(v, these[1])

    β = (α - this.lo) /(this.hi - this.lo)

    bs = bisect(this, β)

    new = IntervalU(v ∪ bs[1] ∪ bs[2])
    new = remove_empties(new)
    sort!(new.v)
    return new

end

function intersect(x :: IntervalU, y :: IntervalU)
    intersects = [intersect(xv, yv) for xv in x.v, yv in y.v]
    return intervalU(intersects[:])
end

function intersect(x :: IntervalU, y :: Interval)
    intersects = [intersect(xv, y) for xv in x.v]
    return intervalU(intersects[:])
end

function intersect(x :: Interval, y :: IntervalU)
    intersects = [intersect(x, yv) for yv in y.v]
    return intervalU(intersects[:])
end

function setdiff(x :: IntervalU, y :: IntervalU)
    yc = complement(y)
    return intersect(x, yc)
end

\(x :: IntervalU, y :: IntervalU) = setdiff(x, y)
\(x :: Interval, y :: IntervalU) = setdiff(intervalU(x), y)
\(x :: IntervalU, y :: Interval) = setdiff(x, intervalU(y))

function ⊂(x :: IntervalU, y :: IntervalU)
    issubs = [xv .⊂ y.v for xv in x.v]
    return all(any.(issubs))
end

function ⊂(x :: Interval, y :: IntervalU)
    if length(y.v) == 1; return x ⊂ y; end
    issubs = x .⊆ y.v
    return any(issubs)
end

function ⊂( x::IntervalU, y :: Interval)
    issubs = x.v .⊂ y
    return all(issubs)
end

function ⊆(x :: IntervalU, y :: IntervalU)
    issubs = [xv .⊆ y.v for xv in x.v]
    return all(any.(issubs))
end

function ⊆(x :: Interval, y :: IntervalU)
    issubs = x .⊆ y.v
    return any(issubs)
end

function ⊆( x::IntervalU, y :: Interval)
    issubs = x.v .⊆ y
    return all(issubs)
end


text:
using DelimitedFiles

include("read_trajec.jl")
include("calcs.jl")

path = "csh2po4_444_standard-pos-1.xyz"
pbc_path = "cdp_444"

pbc = readdlm(pbc_path)
@time coord, atom = read_trajectory(path, true)
if coord == nothing || atom == nothing
    println("Nothing is fine.")
else
    println("Everything is fine.")
end

#@time time_msd, msd = msd_direct(coord, 0.5, atom, "H", 100, 0.3)
#println(msd)
#println(time_msd)
#open("msd_julia.tsv", "w") do outfile
#    writedlm(outfile, [time_msd msd .* 1E4])
#end


text:
"""
Julia program to reverse the bits of a number
Given an integer, reverse its bits in its binary equivalent and
print the new number obtained in its decimal form
"""

function reverse_bits(n)
    rev = 0

    while(n > 0)
        # Shift the bit of the reversed(answer) number to the right
        rev = rev << 1
        #  Stores the temporary lsb of the given number
        rem = n & 1
        # Set the lsb of the answer variable with the stored value
        rev = rem | rev
        # Drops the already processed lsb of the given number
        n = n >> 1
    end
    return rev
end

print("Enter the number? ")
num = readline()
num = parse(Int, num)
rev = reverse_bits(num)
print("The bits-reversed number is: $rev")

"""
Time Complexity: O(n)
Space Complexity: O(1)

SAMPLE INPUT AND OUTPUT

Enter the number? 39
The bits-reversed number is: 57
"""

max_issues_repo_name: adarshnjena/NeoAlgo

max_forks_repo_name: adarshnjena/NeoAlgo


text:
using DandelionWebSockets
using DandelionWebSockets: STATE_OPEN, STATE_CONNECTING, STATE_CLOSING, STATE_CLOSED
using DandelionWebSockets: SocketState, AbstractPonger, SendTextFrame, FrameFromServer
using DandelionWebSockets: masking!
import DandelionWebSockets: write, pong_received, ping_sent
import Base: write, close
using Random
import Random: rand

"InvalidPrecondition signals that a precondition to running the test was not met."
struct InvalidPrecondition <: Exception
    message::String
end

#
# WebSocketHandlerStub
#

"WebSocketHandlerStub acts as a handler for the tests, storing state and incoming messages."
mutable struct WebSocketHandlerStub <: WebSocketHandler
    state::SocketState
    texts::Vector{String}
    binaries::Vector{Vector{UInt8}}
    statesequence::Vector{SocketState}

    WebSocketHandlerStub() = new(STATE_CONNECTING, Vector{String}(), Vector{Vector{UInt8}}(), [])
end

function state_closed(h::WebSocketHandlerStub)
    h.state = STATE_CLOSED
    push!(h.statesequence, STATE_CLOSED)
end

function state_closing(h::WebSocketHandlerStub)
    h.state = STATE_CLOSING
    push!(h.statesequence, STATE_CLOSING)
end

function state_connecting(h::WebSocketHandlerStub)
    h.state = STATE_CONNECTING
    push!(h.statesequence, STATE_CONNECTING)
end

function state_open(h::WebSocketHandlerStub)
    h.state = STATE_OPEN
    push!(h.statesequence, STATE_OPEN)
end

on_text(h::WebSocketHandlerStub, text::String) = push!(h.texts, text)
on_binary(h::WebSocketHandlerStub, binary::AbstractVector{UInt8}) = push!(h.binaries, binary)

function getsingletext(h::WebSocketHandlerStub)
    if length(h.texts) == 0
        throw(InvalidPrecondition("exactly one text was expected, but none were received"))
    end
    if length(h.texts) > 1
        throw(InvalidPrecondition("exactly one text was expected, but more than one was received: $(h.texts)"))
    end

    h.texts[1]
end

function gettextat(h::WebSocketHandlerStub, i::Int)
    if length(h.texts) < i
        throw(InvalidPrecondition(
            "require text at index $i, but only $(length(h.texts)) messages received"))
    end
    h.texts[i]
end

function getbinaryat(h::WebSocketHandlerStub, i::Int)
    if length(h.binaries) < i
        throw(InvalidPrecondition(
            "require binary at index $i, but only $(length(h.binaries)) messages received"))
    end
    h.binaries[i]
end

#
# WriterStub
#

mutable struct FrameIOStub <: IO
    frames::Vector{Frame}
    isopen::Bool

    FrameIOStub() = new(Vector{Frame}(), true)
end

write(w::FrameIOStub, frame::Frame) = push!(w.frames, frame)
close(w::FrameIOStub) = w.isopen = false

function getframe(w::FrameIOStub, i::Int)
    if length(w.frames) < i
        throw(InvalidPrecondition("required frame at index $i, but only has $(length(w.frames))"))
    end
    w.frames[i]
end

function getframeunmasked(w::FrameIOStub, i::Int, mask::AbstractVector{UInt8})
    frame = getframe(w, i)
    masking!(frame.payload, mask)
    frame
end

clearframeswritten(w::FrameIOStub) = w.frames = []
get_no_of_frames_written(w::FrameIOStub) = length(w.frames)

#
# Ponger stub
#

mutable struct PongerStub <: AbstractPonger
    no_of_pongs::Int
    no_of_pings_sent::Int

    PongerStub() = new(0, 0)
end

ping_sent(p::PongerStub) = p.no_of_pings_sent += 1
pong_received(p::PongerStub) = p.no_of_pongs += 1

#
# A fake RNG allows us to deterministically test functions that would otherwise behave
# pseudo-randomly.
#

mutable struct FakeRNG{T} <: AbstractRNG
   values::AbstractArray{T, 1}

   FakeRNG{T}(v::AbstractArray{T, 1}) where {T} = new{T}(copy(v))
end

FakeRNG(::Type{T}) where T = FakeRNG{T}(AbstractArray{T, 1}())

function rand(rng::FakeRNG, ::Type{T}, n::Int) where T
    if isempty(rng.values)
        throw(InvalidPrecondition("FakeRNG requires more random data"))
    end
    splice!(rng.values, 1:n)
end


text:
export PBSManager, SGEManager, addprocs_pbs, addprocs_sge

immutable PBSManager <: ClusterManager
    queue::String
end

immutable SGEManager <: ClusterManager
    queue::String
end

function launch(manager::Union(PBSManager, SGEManager), np::Integer, config::Dict, instances_arr::Array, c::Condition)
    try
        exehome = config[:dir]
        exename = config[:exename]
        exeflags = config[:exeflags]
        home = ENV["HOME"]
        isPBS = isa(manager, PBSManager)
        if manager.queue == ""
            queue = ""
        else
            queue = "-q " * manager.queue
        end

        jobname = "julia-$(getpid())"
        qsub_cmd = `echo "cd $(pwd()) && $exehome/$exename $exeflags"` |> (isPBS ? `qsub -N $jobname $queue -j oe -k o -t 1-$np` : `qsub -N $jobname $queue -terse -j y -t 1-$np`)
        out,qsub_proc = open(qsub_cmd)
        if !success(qsub_proc)
            println("batch queue not available (could not run qsub)")
            return
        end
        id = chomp(split(readline(out),'.')[1])
        if endswith(id, "[]")
            id = id[1:end-2]
        end
        filename(i) = isPBS ? "$home/$jobname-$i.o$id" : "$home/$jobname.o$id.$i"
        print("job id is $id, waiting for job to start ")
        io_objs = cell(np)
        configs = cell(np)
        for i=1:np
            # wait for each output stream file to get created
            fname = filename(i)
            while !isfile(fname)
                print(".")
                sleep(1.0)
            end
            # Hack to get Base to get the host:port, the Julia process has already started.
            cmd = `tail -f $fname`
            cmd.detach = true
            io_objs[i],io_proc = open(cmd)
            io_objs[i].line_buffered = true
            configs[i] = merge(config, {:job => id, :task => i, :iofile => fname, :process => io_proc})
        end
        println("")

        push!(instances_arr, collect(zip(io_objs, configs)))
        notify(c)
    catch e
        println("Error launching qsub")
        println(e)
    end
end

function manage(manager::Union(PBSManager, SGEManager), id::Integer, config::Dict, op::Symbol)
    if op == :finalize
        kill(config[:process])
        if isfile(config[:iofile])
            rm(config[:iofile])
        end
#     elseif op == :interrupt
#         job = config[:job]
#         task = config[:task]
#         # this does not currently work
#         if !success(`qsig -s 2 -t $task $job`)
#             println("Error sending a Ctrl-C to julia worker $id (job: $job, task: $task)")
#         end
    end
end

addprocs_pbs(np::Integer, queue="") = addprocs(np, manager=PBSManager(queue))
addprocs_sge(np::Integer, queue="") = addprocs(np, manager=SGEManager(queue))


text:
module StippleUI

import Genie
import Stipple

using Stipple.Reexport

const assets_config = Genie.Assets.AssetsConfig(package = "StippleUI.jl")

#===#

function theme(; twbpatch::Bool = false) :: Vector{String}
  if ! Genie.Assets.external_assets(assets_config)
    Genie.Router.route(Genie.Assets.asset_path(assets_config, :css, file="quasar.min")) do
      Genie.Renderer.WebRenderable(
        Genie.Assets.embedded(Genie.Assets.asset_file(cwd=normpath(joinpath(@__DIR__, "..")), type="css", file="quasar.min")),
        :css) |> Genie.Renderer.respond
    end

    if twbpatch
      Genie.Router.route(Genie.Assets.asset_path(assets_config, :css, file="bootstrap-patch")) do
        Genie.Renderer.WebRenderable(
          Genie.Assets.embedded(Genie.Assets.asset_file(cwd=normpath(joinpath(@__DIR__, "..")), type="css", file="bootstrap-patch")),
          :css) |> Genie.Renderer.respond
      end
    end
  end

  [
    Stipple.Elements.stylesheet(Genie.Assets.asset_path(assets_config, :css, file="quasar.min")),
    (twbpatch ? Stipple.Elements.stylesheet(Genie.Assets.asset_path(assets_config, :css, file="bootstrap-patch")) : "")
  ]
end

#===#

function deps() :: Vector{String}
  if ! Genie.Assets.external_assets(assets_config)
    Genie.Router.route(Genie.Assets.asset_path(assets_config, :js, file="quasar.umd.min")) do
      Genie.Renderer.WebRenderable(
        Genie.Assets.embedded(Genie.Assets.asset_file(cwd=normpath(joinpath(@__DIR__, "..")), type="js", file="quasar.umd.min")),
        :javascript) |> Genie.Renderer.respond
    end
  end

  [
    Genie.Renderer.Html.script(src="$(Genie.Assets.asset_path(assets_config, :js, file="quasar.umd.min"))")
  ]
end

#===#

include("API.jl")
include("PopupProxies.jl")

include("Avatars.jl")
include("Badges.jl")
include("Banners.jl")
include("BigNumbers.jl")
include("Buttons.jl")
include("Cards.jl")
include("Checkboxes.jl")
include("Chips.jl")
include("DatePickers.jl")
include("Dialogs.jl")
include("Drawers.jl")
include("Editors.jl")
include("ExpansionItems.jl")
include("Forms.jl")
include("FormInputs.jl")
include("Headings.jl")
include("Icons.jl")
include("ImageViews.jl")
include("InnerLoaders.jl")
include("Intersections.jl")
include("Knobs.jl")
include("Layouts.jl")
include("Lists.jl")
include("Menus.jl")
include("Radios.jl")
include("Ranges.jl")
include("ScrollAreas.jl")
include("Selects.jl")
include("Separators.jl")
include("Spaces.jl")
include("Spinners.jl")
include("Tables.jl")
include("Toggles.jl")
include("Tooltips.jl")
include("Uploaders.jl")

#===#

import .API: quasar, quasar_pure, vue, vue_pure, xelem, xelem_pure, csscolors

export quasar, quasar_pure, vue, vue_pure, xelem, xelem_pure, @click, csscolors

@reexport using .Avatars
@reexport using .Badges
@reexport using .Banners
@reexport using .BigNumbers
@reexport using .Buttons
@reexport using .Cards
@reexport using .Checkboxes
@reexport using .Chips
@reexport using .DatePickers
@reexport using .Dialogs
@reexport using .Drawers
@reexport using .Editors
@reexport using .ExpansionItems
@reexport using .Forms
@reexport using .FormInputs
@reexport using .Headings
@reexport using .Icons
@reexport using .ImageViews
@reexport using .InnerLoaders
@reexport using .Intersections
          using .Layouts
@reexport using .Knobs
@reexport using .Lists
@reexport using .Menus
@reexport using .Radios
@reexport using .Ranges
@reexport using .PopupProxies
@reexport using .ScrollAreas
@reexport using .Selects
@reexport using .Separators
@reexport using .Spaces
@reexport using .Spinners
@reexport using .Tables
@reexport using .Toggles
@reexport using .Tooltips
@reexport using .Uploaders

export page_container # from Layouts

#===#

"""
    `@click(expr)`

Defines a js routine that is called by a click of the quasar component.
If a symbol argument is supplied, `@click` sets this value to true.

`@click("savefile = true")` or `@click("myjs_func();")` or `@click(:button)`

Modifers can be appended:
```
@click(:me, :native)
# "v-on:click.native='me = true'"
```
"""
macro click(expr, mode="")
  quote
    x = $(esc(expr))
    m = $(esc(mode))
    if x isa Symbol
      """v-on:click$(m == "" ? "" : ".$m")='$x = true'"""
    else
      "v-on:click='$(replace(x, "'" => raw"\'"))'"
    end
  end
end

#===#

function __init__()
  push!(Stipple.Layout.THEMES, theme)
  Stipple.DEPS[@__MODULE__] = deps
end

end


text:
# ---
# title: 959. Regions Cut By Slashes
# id: problem959
# author: Indigo
# date: 2021-01-25
# difficulty: Medium
# categories: Depth-first Search, Union Find, Graph
# link: <https://leetcode.com/problems/regions-cut-by-slashes/description/>
# hidden: true
# ---
#
# In a N x N `grid` composed of 1 x 1 squares, each 1 x 1 square consists of a
# `/`, `\`, or blank space.  These characters divide the square into contiguous
# regions.
#
# (Note that backslash characters are escaped, so a `\` is represented as
# `"\\"`.)
#
# Return the number of regions.
#
#
#
# **Example 1:**
#
#
#
#     Input: [
#       " /",
#       "/ "
#     ]
#     Output: 2
#     Explanation: The 2x2 grid is as follows:
#     ![](https://assets.leetcode.com/uploads/2018/12/15/1.png)
#
#
# **Example 2:**
#
#
#
#     Input: [
#       " /",
#       "  "
#     ]
#     Output: 1
#     Explanation: The 2x2 grid is as follows:
#     ![](https://assets.leetcode.com/uploads/2018/12/15/2.png)
#
#
# **Example 3:**
#
#
#
#     Input: [
#       "\\/",
#       "/\\"
#     ]
#     Output: 4
#     Explanation: (Recall that because \ characters are escaped, "\\/" refers to \/, and "/\\" refers to /\.)
#     The 2x2 grid is as follows:
#     ![](https://assets.leetcode.com/uploads/2018/12/15/3.png)
#
#
# **Example 4:**
#
#
#
#     Input: [
#       "/\\",
#       "\\/"
#     ]
#     Output: 5
#     Explanation: (Recall that because \ characters are escaped, "/\\" refers to /\, and "\\/" refers to \/.)
#     The 2x2 grid is as follows:
#     ![](https://assets.leetcode.com/uploads/2018/12/15/4.png)
#
#
# **Example 5:**
#
#
#
#     Input: [
#       "//",
#       "/ "
#     ]
#     Output: 3
#     Explanation: The 2x2 grid is as follows:
#     ![](https://assets.leetcode.com/uploads/2018/12/15/5.png)
#
#
#
#
# **Note:**
#
#   1. `1 <= grid.length == grid[0].length <= 30`
#   2. `grid[i][j]` is either `'/'`, `'\'`, or `' '`.
#
#
## @lc code=start
using LeetCode

function regions_by_slashes(grid::Vector{String})
    ## Like Problem 200
    function dfs(grid::Matrix{Int}, r, c)
        grid[r, c] = 1
        (r + 1 <= len * 3) && grid[r+1, c] == 0 && dfs(grid, r+1, c)
        (r - 1 > 0) && grid[r-1, c] == 0 && dfs(grid, r-1, c)
        (c + 1 <= len * 3) && grid[r, c+1] == 0 && dfs(grid, r, c+1)
        (c - 1 > 0) && grid[r, c-1] == 0 && dfs(grid, r, c-1)
        nothing
    end
    len = length(grid)
    new_grid = fill(0, len * 3, len * 3)
    res = 0
    for i in 1:len
        for j in 1:len
            if grid[i][j] == '\\'
                new_grid[3 * i - 2, 3 * j - 2] = new_grid[3 * i - 1, 3 * j - 1] = new_grid[3 * i, 3 * j] = 1
            elseif grid[i][j] == '/'
                new_grid[3 * i - 2, 3 * j] = new_grid[3 * i - 1, 3 * j - 1] = new_grid[3 * i, 3 * j - 2] = 1
            end
        end
    end
    for r in 1 : len * 3
        for c in 1 : len * 3
            if new_grid[r, c] == 0
                res += 1
                dfs(new_grid, r, c)
            end
        end
    end
    res
end
## @lc code=end


text:
# General disclaimer: This whole little mini-lib feels like a poor mans implementation of basic symbolic arithmetics and therefore it is hidden in here so I can just delete it silently when I finally realize how this should really be handled.
# It probably looks like it is capable of much more than it actually can do so if you stumble upon this, use with care. Any bug report are of course still extremely welcome, especially the kind which would motivate scrapping the whole thing in favour of a real solution.

"""
    ΔShape{N,M}

Abstract type describing how an operation changes the shape of the input.

Designed with convolutional and pooling layers in mind to deal with how parameters like kernel size, padding and stride affect the output size.

Typically has nothing to do with the channel or batch dimensions which are handled through other means.

`N` and `M` are integers and represent the number of dimensions of the input and output shape respectively.

For example, the ΔShapes describing a 2D convolutional layer has `N == M == 2` and only describes how the shape of the feature maps change.
"""
abstract type ΔShape{N,M} end

shapeΔ(s::ΔShape) = s.Δshape
ndimsin(::ΔShape{N,M}) where {N,M} = N
ndimsout(::ΔShape{N,M}) where {N,M} = M

ndimsin(s::Tuple{ΔShape, Vararg{ΔShape}}) = ndimsin(first(s))
ndimsout(s::Tuple{ΔShape, Vararg{ΔShape}}) = ndimsout(last(s))

"""
    fshape(s::ΔShape{N,M}, shape::Tuple{N,Integer})
    fshape(s::Tuple{ΔShape{N}, Vararg{ΔShape}}, shape::NTuple{N, Integer})

Return the output shape after passing input of shape `shape` through `s`.
"""
function fshape end

"""
    ShapeMul{N} <: ΔShape{N,N}
    ShapeMul(Δshape)
    ShapeMul(Δshape::Int...)

Size is multiplied by `Δshape`.
"""
struct ShapeMul{N} <: ΔShape{N,N}
    Δshape::NTuple{N, Int}
end
ShapeMul(Δshape::Int...) = ShapeMul(Δshape)
fshape(s::ShapeMul{N}, shape::NTuple{N, Integer}) where N = shape .* shapeΔ(s)

"""
    ShapeDiv{N} <: ΔShape{N,N}
    ShapeDiv(Δshape)
    ShapeDiv(Δshape::Int...)

Size is divided by `Δshape`, rounding up fractions.
"""
struct ShapeDiv{N} <: ΔShape{N,N}
    Δshape::NTuple{N, Int}
end
ShapeDiv(Δshape::Int...) = ShapeDiv(Δshape)
fshape(s::ShapeDiv{N}, shape::NTuple{N, T}) where {N, T<:Integer} = ceil.(T, shape ./ shapeΔ(s))

"""
    ShapeAdd{N} <: ΔShape{N,N}
    ShapeAdd(Δshape)
    ShapeAdd(Δshape::Int...)

Size is added to `Δshape`.
"""
struct ShapeAdd{N} <: ΔShape{N,N}
    Δshape::NTuple{N, Int}
end
ShapeAdd(Δshape::Int...) = ShapeAdd(Δshape)
fshape(s::ShapeAdd{N}, shape::NTuple{N, Integer}) where N = shape .+ shapeΔ(s)

"""
    ShapeFlatten{N} <: ΔShape{N, 0}

Flattens shape so it no longer exists.

Note that `N` can be `Any` since flatten operations typically don't care about the number of input dimensions.
"""
struct ShapeFlatten{N} <: ΔShape{N, 0} end
ShapeFlatten() = ShapeFlatten{Any}()
fshape(::ShapeFlatten, shape) = tuple()

"""
    RevertΔShape{N,M,T<:ΔShape{N,M}} <: ΔShape{N,M}
    RevertΔShape(Δshape)

Revert operation of `Δshape` typically used when the revert operation is undefined.
"""
struct RevertΔShape{N,M,T<:ΔShape{N,M}} <: ΔShape{N,M}
    Δshape::T
end

fshape(::Tuple{}, shape) = shape
fshape(s::Tuple{ΔShape{N}, Vararg{ΔShape}}, shape::NTuple{N, Integer}) where N = foldr(fshape, reverse(s); init=shape)

"""
    revert(s::ΔShape)
    revert(s::Tuple{Vararg{ΔShape}})

Return a `ΔShape` or tuple of `ΔShape`s which reverts the shape change of `s`, i.e `fshape((s..., revert(s)...), x) == x`.

Warning: Not guaranteed to produce exact results for all `ΔShape` types! For example, the operation performed by `ShapeDiv` is not generally invertible due to rounding.
"""
revert(s) = RevertΔShape(s)
revert(s::RevertΔShape) = shapeΔ(s)
revert(s::ShapeAdd) = ShapeAdd(.-shapeΔ(s))
revert(s::ShapeMul) = ShapeDiv(shapeΔ(s))
revert(s::ShapeDiv) = ShapeMul(shapeΔ(s))
revert(s::Tuple{Vararg{ΔShape}}) = reverse(revert.(s))

"""
    combine(s1::ΔShape,s2::ΔShape)

Return a `ΔShape` or tuple of `ΔShape`s which combines `s1` and `s2`, i.e `fshape((s1,s2), x) == fshape(combine(s1,s2), x)`
"""
combine(s1::ΔShape,s2::ΔShape) = s1,s2
combine(s1::Tuple{Vararg{ΔShape}}, s2::ΔShape) = (s1[1:end-1]..., combine(last(s1), s2)...)
combine(s1::ΔShape, s2::Tuple{Vararg{ΔShape}}) = (combine(s1, first(s2))..., s2[2:end]...)
combine(s1::ShapeAdd{N}, s2::ShapeAdd{N}) where N = tuple(ShapeAdd(shapeΔ(s1) .+ shapeΔ(s2)))
combine(s1::T, s2::T) where T <: Union{ShapeDiv{N}, ShapeMul{N}} where N = tuple(T(shapeΔ(s1) .* shapeΔ(s2)))
# Note: Combining ShapeDiv and ShapeMul not generally safe due to rounding when dividing
combine(s1::ShapeMul{N}, s2::ShapeDiv{N}) where N = isdiv(s1, s2) ? tuple(ShapeMul(shapeΔ(s1) .÷ shapeΔ(s2))) : (s1,s2)

isdiv(s1::ΔShape{N,N},s2::ΔShape{N,N}) where N = all(iszero, shapeΔ(s1) .% shapeΔ(s2))

"""
    swapΔshape(s1, s2)

Return a tuple of `ΔShape`s where `s1` and `s2` have swapped places without changing the shape mapping, i.e `fshape((s1,s2), x) == fshape(swapΔshape(s1,s2), x)`.

Basically a helper function for [`orderΔshapes`](@ref). Probably not useful on its own.
"""
swapΔshape(s1, s2) = s1,s2
swapΔshape(s1::T, s2::T) where T <: ΔShape = s2,s1

swapΔshape(s1::ShapeAdd{N}, s2::ShapeMul{N}) where N = s2, ShapeAdd(shapeΔ(s2) .* shapeΔ(s1))
swapΔshape(s1::ShapeMul{N}, s2::ShapeAdd{N}) where N = isdiv(s2, s1) ? (ShapeAdd(shapeΔ(s2) .÷ shapeΔ(s1)), s1) : (s1,s2)

swapΔshape(s1::ShapeAdd{N}, s2::ShapeDiv{N}) where N = isdiv(s1, s2) ? (s2, ShapeAdd(shapeΔ(s1) .÷ shapeΔ(s2))) : (s1,s2)
swapΔshape(s1::ShapeDiv{N}, s2::ShapeAdd{N}) where N = ShapeAdd(shapeΔ(s1) .* shapeΔ(s2)), s1

"""
    filter_noops(s::ΔShape)
    filter_noops(s::ΔShape...)
    filter_noops(s::Tuple{Vararg{ΔShape}})
    filter_noops(s::ΔShape)

Return a tuple of `ΔShape`s where all identity mappings (e.g things like `ShapeAdd(0)`) are removed.

If called with a single identity mapping and empty tuple is returned.
"""
filter_noops(s::ΔShape) = tuple(s)
filter_noops(s::ΔShape...) = filter_noops(s)
filter_noops(s::Tuple{Vararg{ΔShape}}) = mapreduce(filter_noops, (s1,s2) -> (s1...,s2...), s; init=tuple())
filter_noops(s::Union{ShapeMul, ShapeDiv}) = all(x -> x == 1, shapeΔ(s)) ? tuple() : tuple(s)
filter_noops(s::ShapeAdd) = all(x -> x == 0, shapeΔ(s)) ? tuple() : tuple(s)

"""
    orderΔshapes(s::Tuple{Vararg{ΔShape}}; order=allΔshapetypes(s))

Return a tuple of `ΔShape`s which has the same shape mapping as `s` (i.e `fshape(s, x) == fshape(orderΔshapes(s), x)`) but where `ΔShape`s to the extent possible are ordered according to `order`.

Useful to determine whether two arbitrary sequences of `ΔShape`s result in the same shape mapping for all shapes.

Warning: Sort is not stable due to lazy implementation, i.e `orderΔshapes(orderΔshapes(s; order=someorder);order=someorder)` is not guaranteed to return the same thing as `orderΔshapes(s; order=someorder)`.
"""
function orderΔshapes(s::Tuple{Vararg{ΔShape}}; order=allΔshapetypes(s))
    # Yeah, this is bubble sort :/ Given the constraint that ΔShapes can't always swap places along with the fact that swapping generally changes the swapped elements I couldn't think up and other sorting algo that works
    sprev = tuple()
    snew = s
    nlook = length(s)
    # I'm a little worried that lack of stability guarantee will cause this to loop forever, but I have not been able to trigger it despite trying.
    while sprev != snew
        sprev = snew
        nlook -= 1
        for i in 1:nlook
            s1,s2 = snew[i], snew[i+1]
            # Check if s1 shall be before s2 in the ordering
            # Note: We need to "bubble" on equality because s1 might be prevented from bubbling up while s2 isn't.
            if findfirst(st -> s1 isa st, order) >= findfirst(st -> s2 isa st, order)
                snew = (snew[1:i-1]..., swapΔshape(s1,s2)..., snew[i+2:end]...)
            end
        end
    end
    return snew
end

allΔshapetypes(s::T) where T <: ΔShape = T
allΔshapetypes(s::Tuple{Vararg{ΔShape}}) = unique(allΔshapetypes.(s))

"""
    squashshapes(s::ΔShape...; order=allΔshapetypes(s))
    squashshapes(s::Tuple{Vararg{ΔShape}}; order=allΔshapetypes(s))

Return a tuple of `ΔShape`s with the same shape mapping as `s` (i.e `fshape(s, x) == fshape(squashshapes(s), x)`) with as few `ΔShape`s as possible for the given `order`.

Useful to determine whether two arbitrary sequences of `ΔShape`s result in the same shape mapping for all shapes.
"""
squashshapes(s::Tuple{}; order=nothing) = s
squashshapes(s::ΔShape; order=nothing) = tuple(s)
squashshapes(s::ΔShape...; order=allΔshapetypes(s)) = squashshapes(s; order = order)
squashshapes(s::Tuple{Vararg{ΔShape}}; order=allΔshapetypes(s)) = _squashshapes(orderΔshapes(s; order=order))

_squashshapes(s::ΔShape) = tuple(s)
_squashshapes(s::Tuple{ΔShape}) = s
function _squashshapes(s::Tuple{Vararg{ΔShape}})
    squashed = filter_noops(foldr(combine, s)...)
    squashed == s && return s
    isempty(squashed) && return squashed
    return _squashshapes(squashed)
end

Δshapediff(s1,s2) = filter_noops(squashshapes(_Δshapediff(s1,s2)))
_Δshapediff(s1::ΔShape{N}, s2::ΔShape{M}) where {N,M} = N == M ? (revert(s2), s1) : (s1,s2)
_Δshapediff(s1::ΔShape{N}, s2::ΔShape{N}) where N = s1 == s2 ? tuple() : (revert(s2), s1)
function _Δshapediff(s1::Tuple{Vararg{ΔShape}}, s2::Tuple{Vararg{ΔShape}})
    # Pretty crappy heurisic tbh, but I couldn't think of anything better:
    # Step 1: Remove all identical ΔShapes
    # Step 2: Squash shapes and try again
    # Step 3: revert s2 and concat s1

    firstdiff = findfirst(((ss1,ss2)::Tuple) -> ss1 != ss2, collect(zip(s1,s2)))
    firstdiff = isnothing(firstdiff) ? min(length(s1), length(s2))+1 : firstdiff
    sd1 = s1[firstdiff:end]
    sd2 = s2[firstdiff:end]

    isempty(sd1) && isempty(sd2) && return tuple()

    ts1 = allΔshapetypes(sd1)
    ts2 = allΔshapetypes(sd2)
    front = intersect(ts1,ts2)
    back = symdiff(ts1, ts2)

    so1 = squashshapes(sd1; order=vcat(front, back))
    so2 = squashshapes(sd2; order=vcat(front, back))

    (so1 != s1 || so2 != s2) && return _Δshapediff(so1, so2)
    return (revert(so2)..., so1...)
end

"""
    AbstractShapeTrace

Abstract type for tracing shape properties of a graph.

A `MutationVertex v` called with an `AbstractShapeTrace tr` as input will call `visitvertex(tr, v)` or `merge(v, trs...)` if called with more than one `AbstractShapeTrace`.
"""
abstract type AbstractShapeTrace end

"""
    ShapeTrace{T,V1,V2} <: AbstractShapeTrace

Records `ΔShape`s of any visited vertices as well as the origin and destination vertex.
"""
struct ShapeTrace{T,V1,V2} <: AbstractShapeTrace
    origin::V1
    dest::V2
    trace::T
end
ShapeTrace(v) = ShapeTrace(v, v, Δshapes(v))

allΔshapetypes(t::ShapeTrace) = allΔshapetypes(t.trace)
allΔshapetypes(t::Tuple) = unique(mapreduce(allΔshapetypes, vcat, t))

squashshapes(t::ShapeTrace; order=allΔshapetypes(t)) = squashshapes(t.trace;order=order)
# TODO: One can probably do better here when parallel paths can't be squashed
# Now the first instance of such a path will basically prevent squashing of any subsequent paths, even if they are not parallel
squashshapes(t::Tuple; order=allΔshapetypes(t)) = mapfoldr(tt -> squashshapes(tt;order=order), (t1,t2) -> squashshapes(t1,t2; order=order), t)
function squashshapes(t::Tuple{Vararg{ShapeTrace}}; order=allΔshapetypes(t))
     squashed = unique(map(tt -> squashshapes(tt;order=order), t))
     length(squashed) == 1 && return first(squashed)
     return Tuple(squashed) # Danger danger! Graph probably only works for one single input shape
end
# This is the reason for "TODO: One can probably do better here when parallel paths can't be squashed" above
squashshapes(s1, s2; order=missing) = s1, s2
squashshapes(s1::Tuple{Vararg{ΔShape}}, s2::Tuple{Vararg{ΔShape}}; order=allΔshapetypes((s1,s2))) = squashshapes((s1...,s2...); order=order)

visitvertex(tr::ShapeTrace, v) = ShapeTrace(tr.origin, v, (tr.trace..., Δshapes(v)...))

Base.merge(::AbstractVertex, tr::ShapeTrace) = tr
Base.merge(v::AbstractVertex, trs::ShapeTrace...) = ShapeTrace(v, v, tuple(tuple((ShapeTrace(t.origin, v, (t.trace..., Δshapes(v)...)) for t in trs)...)))

"""
    shapetrace(v::AbstractVertex, vs::AbstractVertex...; trfun = v -> ShapeTrace(v))

Return a `AbstractShapeTrace` (default `ShapeTrace`) between `v` and `vs` where `vs` must be input ascendants to `v`. If `vs` is omitted then all input vertices will be used.
"""
function shapetrace(v::AbstractVertex; trfun = v -> ShapeTrace(v))
    ins = filter(v -> isempty(inputs(v)), ancestors(v))
    memo = Dict{AbstractVertex, Any}(ins .=> trfun.(ins))
    return output!(memo, v)
end

function shapetrace(v::AbstractVertex, vs::AbstractVertex...; trfun = v -> ShapeTrace(v))
    memo = Dict{AbstractVertex, Any}(vs .=> trfun.(vs))
    return output!(memo, v)
end

(v::NaiveNASlib.MutationVertex)(trs::AbstractShapeTrace...) = merge(v, trs...)
(v::NaiveNASlib.MutationVertex)(tr::AbstractShapeTrace) = visitvertex(tr, v)

"""
    Δshapes(v::AbstractVertex)

Return a tuple of `ΔShape`s describing the shape mapping of `v`.

More concretely, if `xs = size(x)[sdims]` then `size(v(x))[sdims] == fshape(Δshapes(v), xs)` where `sdims` are the shape dimensions of `x`, e.g. the height and width in case of 2D convolutions.
"""
Δshapes(v::AbstractVertex) = Δshapes(base(v))
Δshapes(::InputVertex) = tuple()
Δshapes(v::MutationVertex) = Δshapes(trait(v), v)
Δshapes(t::DecoratingTrait, v) = Δshapes(base(t), v)
Δshapes(::MutationSizeTrait, v) = _Δshapes(layertype(v), v)

_Δshapes(t::Any, v) = tuple()

function _Δshapes(::FluxConv{N}, v) where N
    c = layer(v)
    ks = size(NaiveNASflux.weights(c))[1:N]
    Δwindow = Δshape_from_window(c, ks, c.dilation, c.pad)
    Δstride = ShapeDiv(c.stride)
    return Δwindow, Δstride
end

_Δshapes(::FluxNoParLayer, v) = _Δshapes(layer(v), v)
_Δshapes(p::Union{MeanPool, MaxPool}, v) = Δshape_from_window(p, p.k, 1, p.pad), ShapeDiv(p.stride)
_Δshapes(::Type{typeof(Flux.flatten)}, v) = _Δshapes_gp(NaiveNASflux.actdim(v) .- 2)
_Δshapes(::GlobalPool, v) = _Δshapes_gp(NaiveNASflux.actdim(v) .- 2)

Δshape_from_window(l::LT, ws::NTuple{N}, dilation::Integer, pad) where {N, LT} = Δshape_from_window(l, ws, ntuple(i -> dilation, N), pad)
function Δshape_from_window(::LT, ws::NTuple{N}, dilation, pad) where {LT, N}
    padact = length(pad) == N ? 2 .* pad : ntuple(i -> pad[2(i-1)+1] + pad[2(i-1)+2], N)
    padref = ntuple(i -> sum(Flux.calc_padding(LT, SamePad(), tuple(ws[i]), dilation[i], 1)), N)
    return ShapeAdd(padact .- padref)
end

function _Δshapes_gp(shapedims)
    udims = unique(shapedims)
    # Graphs with different shapedims are probably not correct, but there is always the chance that some op in between which fixes this does not have a defined ΔShape.
    length(udims) != 1 && return tuple(ShapeFlatten())
    return tuple(ShapeFlatten{first(udims)}())
end


text:
let
    d = SparseCat([:a, :b, :d], [0.4, 0.5, 0.1])
    c = collect(weighted_iterator(d))
    @test c == [:a=>0.4, :b=>0.5, :d=>0.1]
    @test pdf(d, :c) == 0.0
    @test pdf(d, :a) == 0.4
    @test mode(d) == :b
    @test sampletype(d) == Symbol
    @test sampletype(typeof(d)) == Symbol
    @inferred rand(Random.GLOBAL_RNG, d)

    dt = SparseCat((:a, :b, :d), (0.4, 0.5, 0.1))
    c = collect(weighted_iterator(dt))
    @test c == [:a=>0.4, :b=>0.5, :d=>0.1]
    @test pdf(dt, :c) == 0.0
    @test pdf(dt, :a) == 0.4
    @test mode(dt) == :b
    @test sampletype(dt) == Symbol
    @test sampletype(typeof(dt)) == Symbol
    @inferred rand(Random.GLOBAL_RNG, dt)

    rng = MersenneTwister(14)
    samples = Symbol[]
    N = 100_000
    @time for i in 1:N
        push!(samples, rand(rng, d))
    end
    @test isapprox(count(samples.==:a)/N, pdf(d,:a), atol=0.005)
    @test isapprox(count(samples.==:b)/N, pdf(d,:b), atol=0.005)
    @test isapprox(count(samples.==:c)/N, pdf(d,:c), atol=0.005)
    @test isapprox(count(samples.==:d)/N, pdf(d,:d), atol=0.005)

    @test_throws ErrorException rand(Random.GLOBAL_RNG, SparseCat([1], [0.0]))
end


text:
using InternedStrings
using Strided
using ..Pickle: store, save_global, save_get, save_object, memoize, OpCodes, hasref
import ..Pickle: save

const THTensorElType = Union{Float64, Float32, Float16, UInt8,
                             Int8, Int16, Int32, Int16, Bool}

"""
  THsave(file::AbstractString, x)

save data that can be load by `torch.load`. `Array` and
`Strided.StridedView` will be treated as `torch.tensor`.
"""
THsave(file::AbstractString, x) = open(file, "w+") do io
  THsave(TorchPickler(), io, x)
end

function THsave(tp::TorchPickler, io, x)
  global MAGIC, TORCH_PROTOCOL
  store(io, MAGIC; proto=protocol(tp))
  store(io, TORCH_PROTOCOL; proto=protocol(tp))
  store(io,
        Dict(
          "little_endian" => Base.ENDIAN_BOM == 0x04030201,
          "protocol_version" => 1001,
          "type_sizes" => Dict(
            "int"=>4,"long"=>4,"short"=>2
          )
        ); proto=protocol(tp))

  store(tp, io, x)
  store(io, collect(keys(tp.storage.data)); proto=protocol(tp))
  write_storage(io, tp.storage)
end

Strided.StridedView(x::PermutedDimsArray{T,N,perm}) where {T,N,perm} = permutedims(StridedView(x.parent), perm)

save(p::TorchPickler, io::IO, a::A) where {T <: THTensorElType, A <: AbstractArray{T}} = save(p, io, StridedView(a))
function save(p::TorchPickler, io::IO, x::S) where {T <: THTensorElType, S <: StridedView{T}}
  save_global(p, io, i"__torch__.rebuild_tensor")

  write(io, OpCodes.MARK)
  save_storage(p, io, x.parent)
  save(p, io, x.offset)
  save_object(p, io, x.size)
  save_object(p, io, x.strides)
  save(p, io, false)
  save_global(p, io, i"OrderedDict")
  write(io, OpCodes.EMPTY_TUPLE)
  write(io, OpCodes.REDUCE)
  write(io, OpCodes.TUPLE)
  write(io, OpCodes.REDUCE)
end

"dummy type for storage type. This is use to inform julia that some arrays are actually the same
 since `objectid` is not consistent on array with same data."
struct StorageTensor{P}
  ptr::P
end

StorageTensor(x::A) where {A<:AbstractArray} = StorageTensor(pointer(x))

function save_storage(p::TorchPickler, io::IO, x::AbstractArray{T}) where T
  s = StorageTensor(x)
  hasref(p.memo, s) && return save_get(p, io, s)

  write(io, OpCodes.MARK)

  save_object(p, io, i"storage")
  save_global(p, io, i"__torch__.StorageType.$T")
  save_object(p, io, string(objectid(x)))
  save_object(p, io, i"cpu")
  save(p, io, length(x))
  save(p, io, nothing)

  write(io, OpCodes.TUPLE)
  write(io, OpCodes.BINPERSID)

  memoize(p, io, s)
  p.storage[string(objectid(x))] = (jltype2dtype(T), length(x), i"cpu", x)

end

function write_storage(io::IO, st::StorageManager)
  for (k, v) in pairs(st.data)
    _, numel, _, storage = v
    write(io, Int64(numel))
    write(io, reinterpret(UInt8, reshape(storage, numel)))
  end
end

function save(p::TorchPickler, io, x::OrderedDict)
  save_global(p, io, i"OrderedDict")
  write(io, OpCodes.EMPTY_TUPLE)
  write(io, OpCodes.REDUCE)

  write(io, OpCodes.MARK)
  for k in keys(x)
    save_object(p, io, k)
    save_object(p, io, x[k])
  end
  write(io, OpCodes.SETITEMS)

  save_object(p, io, Dict(
    "_metadata" => Dict("version"=>1)
  ))

  write(io, OpCodes.BUILD)
  memoize(p, io, x)
end


text:
"""
`Maker.PhonyTarget`

The type created by `phony()`. Fields expected to be accessed publicly
include:

- `name::UTF8String`
- `dependencies::Vector{UTF8String}`
- `description::UTF8String`
"""
type PhonyTarget <: AbstractTarget
    name::UTF8String
    dependencies::Vector{UTF8String}
    description::UTF8String
    action::Function
    timestamp::DateTime
    funhash::UInt64
    isstale::Bool
end

timestamp(t::PhonyTarget) = DateTime()

function isstale(t::PhonyTarget)
    if t.isstale
        return true
    end
    ds = dependencies(t)
    isempty(ds) ? true : t.timestamp < maximum([timestamp(d) for d in ds])
end

"""
```julia
phony(action::Function, name::AbstractString, dependencies=[])
phony(name::AbstractString, dependencies=[])
```
Define and register a `phony` target for Maker.jl.

- `action` is the function that operates when the target is
  executed.

- `name` refers to the name of the task.

- `dependencies` refers to the name (`AbstractString`) or names
  (`Vector{AbstractString}`) of targets that need to be satisfied
  for this target before running the `action`. These are also referred
  to as prerequisites.

Targets are registered globally.

`phony` targets are generic targets that follow the same rules as a `task`, but
`phony` always has an "ancient" timestamp meaning it won't trigger upstream
changes. It does store a timestamp to  determine whether to execute its own
action based on changes to dependencies.

If the `action` or `dependencies` of a target are redefined, the target will be
marked as stale, and the action will be updated at the next target check.

See also `make`, `file`, `task`, and `variable`. `phony` registers a
`PhonyTarget` type.
"""
phony(action::Function, name::AbstractString, dependencies=UTF8String[]) =
    target(PhonyTarget, name, action, dependencies)
phony(name::AbstractString, dependencies=UTF8String[]) =
    target(PhonyTarget, name, ()->nothing, dependencies)


text:
using LinearAlgebra

const GAS_CONSTANT_SI = 8.31446261815324

"""
    solve_cubic(a::Real, b::Real, c::Real, d::Real)

Finds real roots of a cubic equation
```math
a x^3 + b x^2 + c x + d = 0
```

Return Tuple{Float64, Float64, Float64} where first roots are real,
complex roots are defined to `NaN`

Reference: J. F. Blinn, "How to Solve a Cubic Equation, Part 5: Back to Numerics,"
in IEEE Computer Graphics and Applications, vol. 27, no. 3, pp. 78-89, May-June 2007.
doi: 10.1109/MCG.2007.60
"""
function solve_cubic(a::Real, b::Real, c::Real, d::Real)::NTuple{3,Float64}
    # convert to Ax³ + 3Bx² + 3Cx + D = 0
    A, B, C, D = Float64(a), b / 3.0, c / 3.0, Float64(d)
    δ₁ = A * C - B * B
    δ₂ = A * D - B * C
    lostbits2 = -log2(abs(A*D/(B*C) - 1))
    δ₃ = B * D - C * C
    d13 = δ₁ * δ₃
    d22 = δ₂ * δ₂
    Δ = 4 * d13 - d22
    Δ = 52 + log2(abs(Δ/d22)) < lostbits2+1 ? 0.0 : Δ

    if Δ < 0
        At, Cb, Db = 0., 0., 0. # A-tilde, C-bar, D-bar
        if B^3 * D >= A * C^3
            At, Cb, Db = A, δ₁, -2 * B * δ₁ + A * δ₂
        else
            At, Cb, Db = D, δ₃, -D * δ₂ + 2 * C * δ₃
        end
        T₀ = -copysign(At, Db) * sqrt(-Δ)
        T₁ = -Db + T₀
        p = cbrt(0.5 * T₁)
        q = T₁ == T₀ ? -p : -Cb / p
        x₁ = Cb <= 0 ? p + q : -Db / (p^2 + q^2 + Cb)
        x, w = B^3 * D >= A * C^3 ? (x₁ - B, A) : (-D, x₁ + C)
        return (x/w, NaN, NaN)
    else
        δ₁ == δ₂ == δ₃ == 0 && return (-B/A, -B/A, -B/A)
        sΔ = sqrt(Δ)
        θA, θD = (atan(A*sΔ, 2*B*δ₁ - A*δ₂), atan(D*sΔ, D*δ₂ - 2*C*δ₃)) ./ 3 .|> abs
        sCA, sCD = sqrt(-δ₁), sqrt(-δ₃)
        x₁A, x₁D = 2*sCA*cos(θA), 2*sCD*cos(θD)
        x₃A, x₃D = -sCA*(cos(θA)+sqrt(3)*sin(θA)), -sCD*(cos(θD)+sqrt(3)*sin(θD))
        xlt = (x₁A+x₃A > 2*B) ? x₁A : x₃A
        xst = (x₁D+x₃D < 2*C) ? x₁D : x₃D
        xl, wl = xlt - B, A
        xs, ws = -D, xst + C
        Δ == 0 && return (xs/ws, xl/wl, NaN)
        E = wl * ws
        F = -xl * ws - wl * xs
        G = xl * xs
        xm, wm = C * F - B * G, C * E - B * F
        return (xs/ws, xm/wm, xl/wl)
    end
end

struct BrusilovskyEoSComponent
    # meta information
    name::String

    # physical parameters
    Pc::Float64  # critical pressure
    Tc::Float64  # critical temperature
    acentric_factor::Float64
    RTc::Float64   # R * critical temperature
    molar_mass::Float64  # [kg mol⁻¹] molar mass
    carbon_number::Int64  # [dimless] number of carbons

    # eos parameters
    ac::Float64     # explicit coefficient of the eos a_c
    b::Float64      # explicit coefficient of the eos b
    c::Float64      # explicit coefficient of the eos c
    d::Float64      # explicit coefficient of the eos d
    Psi::Float64    # primary coefficient of the eos - \Psi

    function BrusilovskyEoSComponent(
        ;
        name::String="no name",
        critical_pressure::Real=NaN,
        critical_temperature::Real=NaN,
        acentric_factor::Real=NaN,
        Omegac::Real=NaN,
        Zc::Real=NaN,
        Psi::Real=NaN,
        molar_mass::Real=NaN,
        carbon_number::Int64=0
    )
        alpha = Omegac^3
        beta = Zc + Omegac - 1.0
        ds = sqrt(Omegac - 0.75)
        sigma = -Zc + Omegac * (0.5 + ds)
        delta = -Zc + Omegac * (0.5 - ds)
        RTc = GAS_CONSTANT_SI * critical_temperature
        rtp = RTc / critical_pressure
        ac = alpha * RTc * rtp
        b = beta * rtp
        c = sigma * rtp
        d = delta * rtp
        new(
            name,
            critical_pressure,
            critical_temperature,
            acentric_factor,
            RTc,
            molar_mass,
            carbon_number,
            ac, b, c, d, Psi
        )
    end
end

"""
    brusilovsky_a_coef(substance::BrusilovskyEoSComponent, RT::Real)

Returns EoS coefficient ``a(T, Ψ) = a_c ϕ(T, Ψ)`` of `substance`.
```math
ϕ(T, psi) = [ 1 + Ψ (1 - T_r^0.5) ]^2
```
Reference: Brusylovsky2002[section: 5.5.2 (algorithm step 3), eq: (4.34), see p.142 and p.164]
"""
function brusilovsky_a_coef(substance::BrusilovskyEoSComponent, RT::Real)
    psi = substance.Psi
    phi = ( 1 + psi * (1 - sqrt(RT / substance.RTc)) )^2
    return substance.ac * phi
end

function brusilovsky_compressibility(substance::BrusilovskyEoSComponent,
                                     pressure::Real,
                                     RT::Real,
                                     phase::AbstractChar = 'g')

    prt = pressure / RT
    am = brusilovsky_a_coef(substance, RT) * prt / RT
    bm = substance.b * prt
    cm = substance.c * prt
    dm = substance.d * prt
    zf = solve_cubic(
        1.0,
        cm + dm - bm - 1.0,
        am - bm * cm + cm * dm - bm * dm - dm - cm,
        -bm * cm * dm - cm * dm - am * bm
    )
    if phase == 'g'
        return maximum(z for z in zf if z > 0.0) # NaNs are filtered too
    else
        return minimum(z for z in zf if z > 0.0) # NaNs are filtered too
    end
end

function logΦ(subst::BrusilovskyEoSComponent, V, RT, nmoles)
    b, c, d = subst.b, subst.c, subst.d
    a = brusilovsky_a_coef(subst, RT)
    B, C, D = nmoles .* (b, c, d)
    A = a * nmoles / RT
    return log(1 - B / V) - B / (V - B) +
           A / (C - D) * log((V + C) / (V + D)) +
           A * V / ((V + C) * (V + D))
end

function wilson_psat(subst, RT)
    return subst.Pc * exp(5.37 * (1 + subst.acentric_factor) * (1 - subst.RTc / RT))
end

"""
    brusilovsky_pressure(substance::BrusilovskyEoSComponent, mvol::Real, RT::Real)

Computes pressure of fluid component `substance` at given molar volume `mvol`
and thermal energy `RT`.
"""
function brusilovsky_pressure(substance::BrusilovskyEoSComponent, mvol::Real, RT::Real)
    acoeff = brusilovsky_a_coef(substance, RT)
    b = substance.b
    c = substance.c
    d = substance.d
    return RT / (mvol - b) - acoeff / ((mvol + c) * (mvol + d))
end

const methane = BrusilovskyEoSComponent(
                    name = "CH4",
                    critical_pressure = 4.5992e6,
                    critical_temperature = 190.56,
                    acentric_factor = 0.01142,
                    Omegac = 0.7563,
                    Zc = 0.33294,
                    Psi = 0.37447,
                    molar_mass = 0.016043)

function vt_stability(subst, RT, molar_dens;
                      optmethod::T = DescentMethods.BFGS
                     ) where T<:Type{<:DescentMethods.DescentMethod}
    chempot(rho) = log(rho) - logΦ(subst, 1.0, RT, rho)
    μ0 = chempot(molar_dens)
    p_orig = brusilovsky_pressure(subst, 1/molar_dens, RT) / RT
    function Dfunc!(rho, dro)
        Δμ = chempot(rho[1]) - μ0
        d = Δμ * rho[1] - brusilovsky_pressure(subst, 1/rho[1], RT) / RT
        dro[1] = Δμ
        return d, dro
    end

    P_init = wilson_psat(subst, RT)
    z1 = brusilovsky_compressibility(subst, P_init, RT, 'g')
    z2 = brusilovsky_compressibility(subst, P_init, RT, 'l')
    rho1, rho2 = P_init ./ (RT .* (z1, z2))

    rhovec = [rho1]
    dd = [0.0]

    function maxstep(x, d)
        to_zero = -x[1] / d[1]
        to_b = (1 / subst.b - x[1]) / d[1]
        return minimum(a for a in (to_zero, to_b) if a > 0)
    end
    opt = optmethod(rhovec)
    #opt.α0 = 1e-8
    optresult = DescentMethods.optimize!(opt,
                                         Dfunc!,
                                         rhovec,
                                         gtol = 1e-9,
                                         maxiter = 1000,
                                         maxcalls = 2000,
                                         constrain_step = maxstep)
    p_gas = brusilovsky_pressure(subst, 1/optresult.argument[1], RT) / RT
    if p_gas >= p_orig &&
        norm(x - y for(x, y) in zip(rhovec, optresult.argument)) > 1e-5 * norm(rhovec) #+ sqrt(eps(one(p_orig))) * abs(p_orig)
        return optresult, p_gas, p_orig
    end
    rhovec[1] = rho2
    optresult = DescentMethods.optimize!(opt,
                                         Dfunc!,
                                         rhovec,
                                         gtol = 1e-9,
                                         maxiter = 1000,
                                         maxcalls = 2000,
                                         constrain_step = maxstep)
    p_liq = brusilovsky_pressure(subst, 1/optresult.argument[1], RT) / RT
    return optresult, p_liq, p_orig
end

function vt_flash(subst, RT, molar_dens;
                  optmethod::T = DescentMethods.BFGS
                 ) where T<:Type{<:DescentMethods.DescentMethod}
    chempot(rho) = log(rho) - logΦ(subst, 1.0, RT, rho)
    μ0 = chempot(molar_dens)
    p_orig = brusilovsky_pressure(subst, 1/molar_dens, RT) / RT

    A0 = μ0 * molar_dens - p_orig
    function twophaseA!(rhov, grad)
        mol1, V1 = rhov[1], rhov[2]
        mol2, V2 = molar_dens - mol1, 1 - V1
        p1 = brusilovsky_pressure(subst, V1/mol1, RT) / RT
        p2 = brusilovsky_pressure(subst, V2/mol2, RT) / RT
        μ1, μ2 = chempot(mol1 / V1), chempot(mol2 / V2)
        dA = μ1 * mol1 + μ2 * mol2 - (p1 * V1 + p2 * V2)
        grad[1] = μ1 - μ2
        grad[2] = p2 - p1
        return dA, grad
    end

    stabtest = vt_stability(subst, RT, molar_dens, optmethod = optmethod)

    if stabtest[2] < p_orig #+ sqrt(eps(one(p_orig))) * abs(p_orig)
        return
    end

    znew = stabtest[1].argument
    rhov = [float(molar_dens), one(molar_dens)]
    dd = [0.0, 0.0]

    function maxstep(x, d)
        to_zero = minimum(zip(x, d)) do (xi, di)
            a = -xi / di
            a > 0 ? a : convert(typeof(a), Inf)
        end
        to_one = minimum(zip(x, d, rhov)) do (xi, di, init)
            a = (init - xi) / di
            a > 0 ? a : convert(typeof(a), Inf)
        end
        to_b =
            let Δ1 = subst.b * x[1] - x[2],
                Δ2 = subst.b * (rhov[1] - x[1]) - (1 - x[2]),
                Δ = subst.b * d[1] - d[2]

                minimum((-Δ1 / Δ, Δ2 / Δ)) do x
                    x > 0 ? x : convert(typeof(x), Inf)
                end
            end
        return min(to_zero, to_one, to_b)
    end

    opt = optmethod(rhov)
    #opt.α0 = 1e-8

    arg = DescentMethods.argumentvec(opt)
    gradA = DescentMethods.gradientvec(opt)

    let a = maxstep(rhov, (-znew[1], -1)) / 2
        while a > 1e-20
            arg[1] = rhov[1] - a * znew[1]
            arg[2] = 1 - a
            y, g = twophaseA!(arg, gradA)
            d = g ⋅ (znew[1], one(znew[1]))
            if y < A0 && d > 0
                if opt isa DescentMethods.BFGS
                    opt.xdiff .= arg .- rhov
                    DescentMethods.reset!(opt, arg, abs(dot(opt.xdiff, g)) / dot(g, g))
                end
                break
            end
            a /= 2
        end
    end

    track_file = open("track.vtflash.$(optmethod).$(RT / GAS_CONSTANT_SI)K.$(round(molar_dens, digits=2))mol_m3.txt", "w")
    optresult = DescentMethods.optimize!(opt,
                                         twophaseA!,
                                         arg,
                                         gtol = 1e-9,
                                         maxiter = 1000,
                                         maxcalls = 20000,
                                         constrain_step = maxstep,
                                         reset = false,
                                         track_io = track_file)

    close(track_file)
    xfinal = optresult.argument
    d1 = xfinal[1] / xfinal[2]
    d2 = (rhov[1] - xfinal[1]) / (1 - xfinal[2])
    p1 = brusilovsky_pressure(subst, 1/d1, RT)
    p2 = brusilovsky_pressure(subst, 1/d2, RT)
    if d1 < d2
        return (ρgas = d1, ρliq = d2, s = xfinal[2], pgas = p1, pliq = p2, opt = optresult)
    else
        return (ρgas = d2, ρliq = d1, s = 1 - xfinal[2], pgas = p2, pliq = p1, opt = optresult)
    end
end

let phase_split_cg = vt_flash(methane,
                              GAS_CONSTANT_SI * 110,
                              1000,
                              optmethod = DescentMethods.CGDescent),
    phase_split_bfgs = vt_flash(methane,
                                GAS_CONSTANT_SI * 110,
                                1000,
                                optmethod = DescentMethods.BFGS)

    println("""
    CG: $phase_split_cg
    BFGS: $phase_split_bfgs
    """)
end


text:
module GeneralizedSylvesterSolver
###
# solving a x + b x (c ⊗ c ⊗ ... ⊗ c) = d
# using (I + s^T ⊗ s^T ⊗ ... \otimes s^T \otimes t)x = d
###

using QuasiTriangular
using FastLapackInterface
using KroneckerTools
using LinearAlgebra
export GeneralizedSylvesterWs, generalized_sylvester_solver!

struct GeneralizedSylvesterWs
    ma::Int64
    mb::Int64
    b1::Matrix{Float64}
    c1::Matrix{Float64}
    vs_b::Matrix{Float64}
    vs_c::Matrix{Float64}
    s2::QuasiUpperTriangular{Float64,Matrix{Float64}}
    t2::QuasiUpperTriangular{Float64,Matrix{Float64}}
    work1::Vector{Float64}
    work2::Vector{Float64}
    work3::Vector{Float64}
    result::Matrix{Float64}
    linsolve_ws::LinSolveWs
    dgees_ws_b::DgeesWs
    dgees_ws_c::DgeesWs
    function GeneralizedSylvesterWs(ma::Int64, mb::Int64, mc::Int64, order::Int64)
        if mb != ma
            DimensionMismatch("a has $ma rows but b has $mb rows")
        end
        b1 = Matrix{Float64}(undef, mb,mb)
        c1 = Matrix{Float64}(undef, mc,mc)
        vs_b = Matrix{Float64}(undef, mb,mb)
        vs_c = Matrix{Float64}(undef, mc,mc)
        s2 = QuasiUpperTriangular(Matrix{Float64}(undef, mc,mc))
        t2 = QuasiUpperTriangular(Matrix{Float64}(undef, mb,mb))
        linsolve_ws = LinSolveWs(ma)
        dgees_ws_b = DgeesWs(mb)
        dgees_ws_c = DgeesWs(mc)
        work1 = Vector{Float64}(undef, ma*mc^order)
        work2 = Vector{Float64}(undef, ma*mc^order)
        work3 = Vector{Float64}(undef, ma*mc^order)
        result = Matrix{Float64}(undef, ma,mc^order)
        new(ma, mb, b1, c1, vs_b, vs_c, s2, t2, work1, work2, work3, result, linsolve_ws, dgees_ws_b, dgees_ws_c)
    end
end

function generalized_sylvester_solver!(a::AbstractMatrix,b::AbstractMatrix,c::AbstractMatrix,
                                   d::AbstractMatrix,order::Int64,ws::GeneralizedSylvesterWs)
    copy!(ws.b1,b)
    copy!(ws.c1,c)
    linsolve_core!(a, ws.b1, ws.linsolve_ws)
    linsolve_core_no_lu!(a, d, ws.linsolve_ws)
    dgees!(ws.dgees_ws_b,ws.b1)
    dgees!(ws.dgees_ws_c,ws.c1)
    t = QuasiUpperTriangular(ws.b1)
    A_mul_B!(ws.t2,t,t)
    s = QuasiUpperTriangular(ws.c1)
    A_mul_B!(ws.s2,s,s)
    at_mul_b_kron_c!(ws.result, ws.dgees_ws_b.vs, d, ws.dgees_ws_c.vs, order, ws.work2, ws.work3)
    copy!(d,ws.result)
    solve1!(1.0, order, t, ws.t2, s, ws.s2, vec(d), ws)
    a_mul_b_kron_ct!(ws.result, ws.dgees_ws_b.vs, d, ws.dgees_ws_c.vs, order, ws.work2, ws.work3)
    copy!(d,reshape(ws.result, size(a, 1), size(c, 2)^order))
end

function solver!(t::QuasiUpperTriangular,s::QuasiUpperTriangular,d::AbstractVector,order::Int64,ws::GeneralizedSylvesterWs)
    s2 = QuasiUpperTriangular(s*s)
    t2 = QuasiUpperTriangular(t*t)
    solve1!(1.0,order,t,t2,s,s2,d,ws)
    d
end

function solve1!(r::Float64, depth::Int64, t::AbstractArray{Float64,2}, t2::AbstractArray{Float64,2}, s::AbstractArray{Float64,2}, s2::AbstractArray{Float64,2}, d::AbstractVector{Float64}, ws::GeneralizedSylvesterWs)
    m = size(t,2)
    n = size(s,1)
    if depth == 0
        I_plus_rA_ldiv_B!(r,t,d)
    else
        nd = m*n^(depth-1)
        nd2 = 2*nd
        drange1 = 1:nd
        drange2 = 1:nd2
        i = 1
        while i <= n
            if i == n || s[i+1,i] == 0
                dv = view(d,drange1)
                solve1!(r*s[i,i],depth-1,t,t2,s,s2,dv,ws)
                if i < n
                    solvi_real_eliminate!(i,n,nd,drange1,depth-1,r,t,s,d,ws)
                end
                drange1 = drange1 .+ nd
                drange2 = drange2 .+ nd
                i += 1
            else
                dv = view(d,drange2)
                solvii(r*s[i,i],r*s[i+1,i],r*s[i,i+1],depth-1,t,t2,s,s2,dv,ws)
                if i < n - 1
                    solvi_complex_eliminate!(i,n,nd,drange1,depth-1,r,t,s,d,ws)
                end
                drange1 = drange1 .+ nd2
                drange2 = drange2 .+ nd2
                i += 2
            end
        end
    end
end

"""
solvi_real_eliminate!(i::Int64, n::Int64, nd::Int64, drange::UnitRange{Int64},
                      depth::Int64, r::Float64, t::QuasiUpperTriangular, s::QuasiUpperTriangular,
                      d::AbstractVector, ws::GeneralizedSylvesterWs)
    updates d[k] with d[k] - ∑_{j=i+1}^n r*s[i, j]*(s ⊗ s ⊗ ... ⊗ s ⊗ t)*d with
    depth occurences of s
"""
function solvi_real_eliminate!(i::Int64, n::Int64, nd::Int64, drange::UnitRange{Int64},
                               depth::Int64, r::Float64, t::QuasiUpperTriangular, s::QuasiUpperTriangular,
                               d::AbstractVector, ws::GeneralizedSylvesterWs)
    work1 = ws.work1
    work2 = ws.work2
    work3 = ws.work3
    kron_at_kron_b_mul_c!(work1,1,s,depth,t,d,drange[1],work2,work3,1)
    k1 = drange[1] + nd
    @inbounds for j = i+1:n
        m = r*s[i,j]
        @simd for k2 = 1:nd
            d[k1] -= m*work1[k2]
            k1 += 1
        end
    end
end

function solvi_complex_eliminate!(i::Int64,n::Int64,nd::Int64,drange::UnitRange{Int64},
                                  depth::Int64,r::Float64,t::QuasiUpperTriangular,s::QuasiUpperTriangular,
                                  d::AbstractVector,ws::GeneralizedSylvesterWs)
    work1 = ws.work1
    work2 = ws.work2
    work3 = ws.work3
    kron_at_kron_b_mul_c!(work1, 1, s, depth, t, d, drange[1], work2, work3, 1)
    drange = drange .+ nd
    kron_at_kron_b_mul_c!(work1, nd+1, s, depth, t, d, drange[1], work2, work3, 1)
    k1 = drange[1] + nd
    @inbounds for j = i + 2 : n
        m1 = r*s[i,j]
        m2 = r*s[i+1,j]
        @simd for k2 = 1:nd
            d[k1] -= m1*work1[k2] + m2*work1[k2 + nd]
            k1 += 1
        end
    end
end

function solvii(alpha::Float64,beta1::Float64,beta2::Float64,depth::Int64,
                t::QuasiUpperTriangular,t2::QuasiUpperTriangular,s::QuasiUpperTriangular,
                s2::QuasiUpperTriangular,d::AbstractVector,ws::GeneralizedSylvesterWs)
    m = size(t,2)
    n = size(s,1)
    nd = m*n^depth
    transformation1(alpha,beta1,beta2,depth,t,s,d,ws)
    dv = view(d,1:nd)
    solviip(alpha,sqrt(-beta1*beta2),depth,t,t2,s,s2,dv,ws)
    dv = view(d,nd+1:2*nd)
    solviip(alpha,sqrt(-beta1*beta2),depth,t,t2,s,s2,dv,ws)
end

function transformation1(a::Float64,b1::Float64,b2::Float64,depth::Int64,
                         t::QuasiUpperTriangular,s::QuasiUpperTriangular,
                         d::AbstractVector,ws::GeneralizedSylvesterWs)
    m = size(t, 2)
    n = size(s, 1)
    nd = m*n^depth
    copyto!(ws.work3, d)
    drange = 1:nd
    d1 = view(ws.work3, drange)
    d2 = view(ws.work3, drange .+ nd)
    work = view(ws.work2,drange)
    kron_at_kron_b_mul_c!(s,depth,t,d1,work)
    kron_at_kron_b_mul_c!(s,depth,t,d2,work)
    @inbounds @simd for i = drange
        d[i] += a*d1[i] - b1*d2[i]
        d[i+nd] += -b2*d1[i] + a*d2[i]
    end
end

diag_zero_sq = 1e-30

function solviip(alpha::Float64,beta::Float64,depth::Int64,t::QuasiUpperTriangular,t2::QuasiUpperTriangular,
                 s::QuasiUpperTriangular,s2::QuasiUpperTriangular,d::AbstractVector,ws::GeneralizedSylvesterWs)
    m = size(t,2)
    n = size(s,1)
    if beta*beta < diag_zero_sq
        solve1!(alpha,depth,t,t2,s,s2,d,ws)
        solve1!(alpha,depth,t,t2,s,s2,d,ws)
        return
    end

    if depth == 0
        I_plus_rA_plus_sB_ldiv_C!(2*alpha,alpha*alpha+beta*beta,t,t2,d)
    else
        nd = m*n^(depth-1)
        nd2 = 2*nd
        drange1 = 1:nd
        drange2 = 1:nd2
        i = 1
        while i <= n
            if i == n || s[i+1,i] == 0
                dv = view(d,drange1)
                if s[i,i]*s[i,i]*(alpha*alpha+beta*beta) > diag_zero_sq
                    solviip(s[i,i]*alpha,s[i,i]*beta,depth-1,t,t2,s,s2,dv,ws)
                end
                if i < n
                    solviip_real_eliminate!(i,n,nd,drange1,depth-1,alpha,beta,t,t2,s,s2,d,ws)
                end
                drange1 = drange1 .+ nd
                drange2 = drange2 .+ nd
                i += 1
            else
                dv = view(d,drange2)
                solviip2(alpha,beta,s[i,i],s[i+1,i],s[i,i+1],depth,t,t2,s,s2,dv,ws)
                if i < n - 1
                    solviip_complex_eliminate!(i,n,nd,drange1,depth-1,alpha,beta,t,t2,s,s2,d,ws)
                end
                drange1 = drange1 .+ nd2
                drange2 = drange2 .+ nd2
                i += 2
            end
        end
    end
end

function solviip_real_eliminate!(i::Int64,n::Int64,nd::Int64,drange::UnitRange{Int64},
                                 depth::Int64,alpha::Float64,beta::Float64,t::QuasiUpperTriangular,
                                 t2::QuasiUpperTriangular,s::QuasiUpperTriangular,
                                 s2::QuasiUpperTriangular,d::AbstractVector,ws::GeneralizedSylvesterWs)
    y1 = view(ws.work1,drange)
    y2 = view(ws.work1,drange .+ nd)
    copyto!(y1,1,d,drange[1],nd)
    copyto!(y2,1,d,drange[1],nd)
    work = view(ws.work2,1:length(drange))
    kron_at_kron_b_mul_c!(s,depth,t,y1,work)
    kron_at_kron_b_mul_c!(s2,depth,t2,y2,work)
    k1 = drange[1] + nd
    @inbounds for j = i+1:n
        m1 = 2*alpha*s[i,j]
        m2 = (alpha*alpha+beta*beta)*s2[i,j]
        @simd for k2 = 1:nd
            d[k1] -= m1*y1[k2] + m2*y2[k2]
            k1 += 1
        end
    end
end

function solviip2(alpha::Float64,beta::Float64,gamma::Float64,delta1::Float64,delta2::Float64,
                  depth::Int64,t::QuasiUpperTriangular,t2::QuasiUpperTriangular,
                  s::QuasiUpperTriangular,s2::QuasiUpperTriangular,d::AbstractVector,ws::GeneralizedSylvesterWs)
    m = size(t,2)
    n = size(s,1)
    aspds = alpha*alpha + beta*beta
    gspds = gamma*gamma - delta1*delta2
    nd = m*n^(depth-1)
    dv1 = view(d,1:nd)
    dv2 = view(d,nd+1:2*nd)
    if aspds*gspds > diag_zero_sq
        transform2(alpha, beta, gamma, -delta1, -delta2, nd, depth, t, t2, s, s2, d,ws)

        delta = sqrt(-delta1*delta2)
	a1 = alpha*gamma - beta*delta
	b1 = alpha*delta + gamma*beta
	a2 = alpha*gamma + beta*delta
	b2 = alpha*delta - gamma*beta
	solviip(a2, b2, depth-1, t, t2, s, s2, dv1, ws);
	solviip(a1, b1, depth-1, t, t2, s, s2, dv1, ws);
        solviip(a2, b2, depth-1, t, t2, s, s2, dv2, ws);
        solviip(a1, b1, depth-1, t, t2, s, s2, dv2, ws);
    end
end

function transform2(alpha::Float64, beta::Float64, gamma::Float64, delta1::Float64, delta2::Float64,
                    nd::Int64, depth::Int64, t::QuasiUpperTriangular, t2::QuasiUpperTriangular,
                    s::QuasiUpperTriangular, s2::QuasiUpperTriangular, d::AbstractVector, ws::GeneralizedSylvesterWs)
    d1 = ws.work1
    kron_at_kron_b_mul_c!(d1,1,s,depth-1,t,d,1,ws.work2,ws.work3,1)
    kron_at_kron_b_mul_c!(d1,nd+1,s,depth-1,t,d,nd+1,ws.work2,ws.work3,1)

    m1 = 2*alpha*gamma
    m2 = 2*alpha*delta1
    @inbounds @simd for i = 1:nd
        dtmp = d1[i]
        d1[i] = d[i] + m1*dtmp + m2*d1[i+nd]
        d1[i+nd] = d[i+nd] + 2*alpha*(delta2*dtmp + gamma*d1[i+nd])
    end

    kron_at_kron_b_mul_c!(d,1,s2,depth-1,t2,d,1,ws.work2,ws.work3,1)
    kron_at_kron_b_mul_c!(d,nd+1,s2,depth-1,t2,d,nd+1,ws.work2,ws.work3,1)

    aspds = alpha*alpha + beta*beta;
    gspds = gamma*gamma + delta1*delta2;
    m1 = aspds*gspds
    m2 = 2*aspds*gamma*delta1
    m3 = 2*aspds*gamma*delta2
    @inbounds @simd for i = 1:nd
        dtmp = d[i]
        d[i] = d1[i] + m1*dtmp + m2*d[i+nd]
        d[i+nd] = d1[i+nd] + m3*dtmp + m1*d[i+nd]
    end
end

"""
    solviip_complex_eliminate!(i,n,nd,drange,depth,alpha,beta,t,t2,s,s2,d)

perfoms elimination after solving for a complex diagonal block of size 2*n^depth

d n^(depth+2) x 1

The solution is stored in d[drange; drange + nd]

The function updates d[i*nd+1:n*nd]
"""
function solviip_complex_eliminate!(i::Int64,n::Int64,nd::Int64,drange::UnitRange{Int64},depth::Int64,
                                    alpha::Float64,beta::Float64,t::QuasiUpperTriangular,t2::QuasiUpperTriangular,
                                    s::QuasiUpperTriangular,s2::QuasiUpperTriangular,d::AbstractVector,ws::GeneralizedSylvesterWs)
    y11 = view(ws.work1,drange)
    y12 = view(ws.work1, drange .+ nd)
    copyto!(y11,1,d,drange[1],nd)
    copyto!(y12,1,d,drange[1],nd)
    drange = drange .+ nd
    y21 = view(ws.work2,drange)
    y22 = view(ws.work2, drange .+ nd)
    copyto!(y21,1,d,drange[1],nd)
    copyto!(y22,1,d,drange[1],nd)
    work = view(ws.work3,drange)

    kron_at_kron_b_mul_c!(s,depth,t,y11,work)
    kron_at_kron_b_mul_c!(s2,depth,t2,y12,work)
    kron_at_kron_b_mul_c!(s,depth,t,y21,work)
    kron_at_kron_b_mul_c!(s2,depth,t2,y22,work)

    alpha2beta2 = alpha*alpha + beta*beta
    k1 = drange[1] + nd
    @inbounds for j = i+2:n
        m1 = 2*alpha*s[i,j]
        m2 = alpha2beta2*s2[i,j]
        m3 = 2*alpha*s[i+1,j]
        m4 = alpha2beta2*s2[i+1,j]
        @simd for k2 = 1:nd
            d[k1] -= m1*y11[k2] + m2*y12[k2] + m3*y21[k2] + m4*y22[k2]
            k1 += 1
        end
    end
end

end


text:
include("PadTest.jl")
include("../src/PositionalEncoder.jl")

Enc = Encoder(TimeSteps, EmbDims)

SeqLens = SequenceLengths(X)
mask = Flux.param(SequenceMask(output, SeqLens))
@time output = Enc(output, mask) .+ output


text:
"""
Convenience code to store info about what should be (re)calculated, saved to disk and/or plotted,
and what can be skipped.

Author: Eric Ford
Created: Sept 2020
"""

"""
The PipelinePlan stores what work needs to be done, what plots should be made, and what data/plots should be saved to disk.
Users will querty it via  need_to(plan, symbol), make_plot(plan, symbol), save_plot(plan, symbol), and save_data(plan, symbol).
Users can update the plan with make_plot!(plan, symbol),  dont_make_plot!(plan, symbol), etc.

In principle, we could cache data here, too.  But I'm not sure if we'll actulaly use that.  So consider the cache as experimental.
"""
module Pipeline

export PipelinePlan
export make_plot, save_plot, save_data, need_to, has_cache

export make_plot!, dont_make_plot!, make_all_plots!,  make_no_plots!
export save_plot!, dont_save_plot!, save_all_plots!,  save_no_plots!
export save_data!, dont_save_data!, save_all_data!,   save_no_data!
export need_to!,   dont_need_to!,   reset_all_needs!, reset_no_needs!
export has_cache,  read_cache, set_cache!, reset_all_cache!

# Default steps for the PipelinePlan to track
default_pipeline_input_collect(keys)  = [ :read_spectra, :read_line_list ]
default_pipeline_work_collect(keys)   = [ :extract_orders, :clean_line_list_tellurics, :fit_lines, :clean_line_list_blends  ]
default_pipeline_output_collect(keys) = [ :ccf_pass0, :ccf_total, :rvs_ccf_total, :ccf_orders, :rvs_ccf_orders, :scalpels, :template, :dcpca ]
default_pipeline_need_to_collect(keys) = vcat( default_pipeline_input_collect(keys),  default_pipeline_work_collect(keys), default_pipeline_output_collect(keys) )

DictSB = Dict{Symbol,Bool}
CacheT = Dict{Symbol,Any}

struct PipelinePlan
    make_plot::DictSB
    save_plot::DictSB
    save_data::DictSB
    need_to::DictSB
    cache::CacheT
end

make_all_plots_dict = DictSB( zip(default_pipeline_output_collect(keys),  trues(length(default_pipeline_output_collect(keys)))  ) )
save_all_plots_dict = DictSB( zip(default_pipeline_output_collect(keys),  trues(length(default_pipeline_output_collect(keys)))  ) )
save_all_files_dict = DictSB( zip(default_pipeline_output_collect(keys),  trues(length(default_pipeline_output_collect(keys)))  ) )
need_to_do_everything_dict = DictSB( zip(default_pipeline_need_to_collect(keys), trues(length(default_pipeline_need_to_collect(keys)))  ) )

make_no_plots_dict = DictSB( zip(default_pipeline_output_collect(keys),  falses(length(default_pipeline_output_collect(keys)))  ) )
save_no_plots_dict = DictSB( zip(default_pipeline_output_collect(keys),  falses(length(default_pipeline_output_collect(keys)))  ) )
save_no_files_dict = DictSB( zip(default_pipeline_output_collect(keys),  falses(length(default_pipeline_output_collect(keys)))  ) )
need_to_do_nothing_dict = DictSB( zip(default_pipeline_need_to_collect(keys), falses(length(default_pipeline_need_to_collect(keys)))  ) )

default_pipeline_make_plot = make_all_plots_dict
default_pipeline_save_plot = save_no_plots_dict
default_pipeline_save_file = save_no_files_dict
default_pipeline_need_to   = need_to_do_everything_dict
default_pipeline_cached_data = CacheT()

function PipelinePlan(; make_plot::DictSB = default_pipeline_make_plot,
                        save_plot::DictSB = default_pipeline_save_plot,
                        save_file::DictSB = default_pipeline_save_file,
                        need_to::DictSB   = default_pipeline_need_to,
                        cached_data::Dict{Symbol,Any} = default_pipeline_cached_data )
    PipelinePlan(make_plot, save_plot, save_file, need_to, cached_data )
end

make_plot(p::PipelinePlan, s::Symbol) = haskey(p.make_plot,s) && p.make_plot[s]
make_plot!(p::PipelinePlan, s::Symbol) = p.make_plot[s] = true
dont_make_plot!(p::PipelinePlan, s::Symbol) = p.make_plot[s] = false
make_all_plots!(p::PipelinePlan ) = map(k -> p.make_plot[k]=true, collect(keys(p.make_plot)) )
make_no_plots!(p::PipelinePlan ) = map(k -> p.make_plot[k]=false, collect(keys(p.make_plot)) )

save_plot(p::PipelinePlan, s::Symbol) = haskey(p.save_plot,s) && p.save_plot[s]
save_plot!(p::PipelinePlan, s::Symbol) = p.save_plot[s] = true
dont_save_plot!(p::PipelinePlan, s::Symbol) = p.save_plot[s] = false
save_all_plots!(p::PipelinePlan ) = map(k -> p.save_plot[k]=true, collect(keys(p.save_plot)) )
save_no_plots!(p::PipelinePlan ) = map(k -> p.save_plot[k]=false, collect(keys(p.save_plot)) )

save_data(p::PipelinePlan, s::Symbol) = haskey(p.save_data,s) && p.save_data[s]
save_data!(p::PipelinePlan, s::Symbol) = p.save_data[s] = true
dont_save_data!(p::PipelinePlan, s::Symbol) = p.save_data[s] = false
save_all_data!(p::PipelinePlan ) = map(k -> p.save_data[k]=true, collect(keys(p.save_data)) )
save_no_data!(p::PipelinePlan ) = map(k -> p.save_data[k]=false, collect(keys(p.save_data)) )

need_to(p::PipelinePlan, s::Symbol) = haskey(p.need_to,s) && p.need_to[s]
function need_to!(p::PipelinePlan, s::Symbol)
    p.need_to[s] = true
    invalidate!(p,s)
end
dont_need_to!(p::PipelinePlan, s::Symbol) = p.need_to[s] = false
reset_all_needs!(p::PipelinePlan) = map(k -> p.need_to[k]=true, collect(keys(p.need_to)) )
reset_no_needs!(p::PipelinePlan) = map(k -> p.need_to[k]=false, collect(keys(p.need_to)) )

has_cache(p::PipelinePlan, s::Symbol ) = haskey(p.cache,s)
read_cache(p::PipelinePlan, s::Symbol ) = p.cache[s]
set_cache!(p::PipelinePlan, s::Symbol, data ) = p.cache[s] = data
reset_all_cache!(p::PipelinePlan ) = p.cache = needs_to_do_everything

# TODO: Eventually make invalidate! more intelligent, so it automatically invalidates other steps/cached data that depend on this step.
invalidate!(p::PipelinePlan, s::Symbol ) = delete!(p.cache, s)

end # module Pipeline

max_stars_repo_name: RvSpectML/RvSpectMLTemplate

max_issues_repo_name: RvSpectML/RvSpectMLTemplate

max_forks_repo_name: RvSpectML/RvSpectMLTemplate


text:
using Distributions, Calculus, SpecialFunctions, Plots; pyplot()

xGrid = -5:0.01:5

PhiA(x) = 0.5*(1+erf(x/sqrt(2)))
PhiB(x) = cdf(Normal(),x)

println("Maximum difference between two CDF implementations: ",
            maximum(PhiA.(xGrid) - PhiB.(xGrid)))

normalDensity(z) = pdf(Normal(),z)

d0 = normalDensity.(xGrid)
d1 = derivative.(normalDensity,xGrid)
d2 = second_derivative.(normalDensity, xGrid)

plot(xGrid, [d0 d1 d2], c=[:blue :red :green],label=[L"f(x)" L"f'(x)" L"f''(x)"])
plot!([-5,5],[0,0],  color=:black, lw=0.5, xlabel="x", xlims=(-5,5), label="")

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
using FractionalDiffEq, Plots
function chua!(du, x, p, t)
    a = 10.725; b = 10.593
    c = 0.268
    m0 = -1.1726
    m1 = -0.7872
    du[1] = a*(x[2]-x[1]-(m1*x[1]+0.5*(m0-m1)*(abs(x[1]+1)-abs(x[1]-1))))
    du[2] = x[1]-x[2]+x[3]
    du[3] = -b*x[2]-c*x[3]
end
α = [0.93, 0.99, 0.92];
x0 = [0.2; -0.1; 0.1];
h = 0.01; tspan = (0, 100);
prob = FODESystem(chua!, α, x0, tspan)
result = solve(prob, h, NonLinearAlg())

plot(result[:, 1], result[:, 2], title="Chua System", legend=:bottomright)


text:
@kwdef mutable struct FxpAddHandler <: AbstractHandler
    id = 0
end

@register(FxpAddHandler, +, [Number, Number] => [Number])

function (handler::FxpAddHandler)(netlist::Netlist,
                                  inputs::Vector{Variable},
                                  outputs::Vector{Variable})
    # check size
    if getsize(netlist, getname(inputs[1])) != (1, 1) && getsize(netlist, getname(inputs[1])) != (1, 1)
        error("FXP add does not support matrix hardware generation yet.")
    end

    outstring = """
        $stdcomment
        // BEGIN f_add$(handler.id)
        fxp_add #(
                .BIT_WIDTH(BIT_WIDTH)
            ) f_add$(handler.id) (
                .a($(inputs[1].name)),
                .b($(inputs[2].name)),
                .y($(outputs[1].name))
            );
        // END f_add$(handler.id)
        \n"""

    handler.id += 1

    return outstring
end