text:
using Knet, Test, Base.Iterators, Printf, LinearAlgebra, Random, CUDA, IterTools
Knet.atype() = KnetArray{Float32}
include("tree.jl")
import Pkg; Pkg.add("PyCall")
using PyCall
@pyimport pickle
using JLD2

import Pkg
Pkg.add("PyCall")

@pyimport pickle

function mypickle(filename, obj)
    out = open(filename,"w")
    pickle.dump(obj, out)
    close(out)
 end

function myunpickle(filename)
    r = nothing
    @pywith pybuiltin("open")(filename,"rb") as f begin
        r = pickle.load(f)
    end
    return r
end

mutable struct SymbolsManager
    symbol2idx#::Dict{String,Int}
    idx2symbol#::Dict{Int,String}
    vocab_size#::Int
    whether_add_special_tags#::Bool
end

function SymbolsManager(whether_add_special_tags)
    symbol2idx = Dict()
    idx2symbol = Dict()
    vocab_size = 0
    symbols_manager = SymbolsManager(symbol2idx, idx2symbol, vocab_size, whether_add_special_tags)
    if whether_add_special_tags
        add_symbol(symbols_manager,"<S>")
        add_symbol(symbols_manager, "<E>")
        add_symbol(symbols_manager, "<U>")
        # non-terminal symbol = 4
        add_symbol(symbols_manager, "<N>")

    end
    #return SymbolsManager(symbol2idx, idx2symbol, vocab_size, whether_add_special_tags)
    return symbols_manager
end

function add_symbol(m::SymbolsManager, s)
    if s ∉ keys(m.symbol2idx)
        m.vocab_size = m.vocab_size + 1
        #s_man.symbol2idx[s] = s_man.vocab_size
        merge!(m.symbol2idx, Dict(s => m.vocab_size))
        #s_man.idx2symbol[s_man.vocab_size] = s
        merge!(m.idx2symbol, Dict(m.vocab_size=>s))

    end
    return m.symbol2idx[s]
end

function get_symbol_idx(s_man::SymbolsManager, s)
    if s ∉ keys(s_man.symbol2idx)
        if s_man.whether_add_special_tags
            return s_man.symbol2idx["<U>"]
        else
            print("this should never be reached (always add <U>")
            return 0
        end
    end
    return s_man.symbol2idx[s]
end

function get_idx_symbol(s_man::SymbolsManager, idx)
    if idx not in keys(s_man.idx2symbol)
        return "<U>"
    end
    return s_man.idx2symbol[idx]
end

function init_from_file(s_man::SymbolsManager, fn, min_freq, max_vocab_size)
    #print("loading vocabulary file: {}\n".format(fn))
    #f = open(fn, "r")
    open(fn) do file
    #for line in readlines(f)
        for line in eachline(file)
            l_list = split(strip(line), "\t")
            #c = Int(l_list[2])
	    c = parse.(Int,l_list[2] )
	    #c = l_list[1]
            if c >= min_freq
                add_symbol(s_man, l_list[1])
            end
            if s_man.vocab_size >= max_vocab_size
                break
            end
        end
    end
end

function get_symbol_idx_for_list(s_man::SymbolsManager,l)
    r = []
    for i in 1:(length(l))
        push!(r, get_symbol_idx(s_man, l[i]))
    end
    return r
end

struct MinibatchLoader
    enc_batch_list
    enc_len_batch_list
    dec_batch_list
    num_batch

end

function MinibatchLoader(opt, mode, using_gpu)
    #data = Pickle.load(open("test.pkl"))
    #data = myunpickle(opt["data_dir"] * "/" * string(mode) * ".pkl")
    @load  (opt["data_dir"] * "/" * string(mode) *".jld2") data
    #(opt["data_dir"] * "/" * string(mode) * ".pkl")
    if size(data,1) % opt["batch_size"] != 0
        n = length(data)
        for i in 1:(length(data)%opt["batch_size"])
            #data.insert(n-i+1, data[n-i+1])
            insert!(data, n-i+1, data[n-i+1])
        end
    end
    enc_batch_list = []
    enc_len_batch_list = []
    dec_batch_list = []
    p = 0
    while p + opt["batch_size"] <= size(data,1)
        max_len = length(data[p + opt["batch_size"] ][1])
        m_text = zeros(Float64, (opt["batch_size"], max_len + 2))
        if using_gpu
            m_text = Knet.KnetArray(m_text)
        end
        enc_len_list = []
        m_text[:,1] .= 1
        for i in 1:(opt["batch_size"] )
            w_list = data[p + i][1]
	    print(w_list)
	    println(m_text)
	    println(size(m_text))
	    println(size(w_list))
            for j in range(1, stop=length(w_list))
                m_text[i,j+1] = w_list[length(w_list) - j + 1]
            end
            for j in range(length(w_list)+2,stop=max_len+2)
                m_text[i, j] = 2
            end
            push!(enc_len_list, length(w_list)+2)
        end
        push!(enc_batch_list, m_text)
        push!(enc_len_batch_list, enc_len_list)
        tree_batch = []
        for i in range(1, stop=opt["batch_size"])
            push!(tree_batch, data[p+i][3])
        end
        push!(dec_batch_list, tree_batch)
        p += opt["batch_size"]
    end
    num_batch = size(enc_batch_list,1)
    return MinibatchLoader(enc_batch_list, enc_len_batch_list, dec_batch_list, num_batch)
end

function random_batch(mbl::MinibatchLoader)
    p = rand(1:mbl.num_batch,  1)
    return mbl.enc_batch_list[p], mbl.enc_len_batch_list[p], mbl.dec_batch_list[p]
end

function all_batch(mbl::MinibatchLoader)
    r = []
    for p in 1:(mbl.num_batch)
        push!(r, [mbl.enc_batch_list[p], mbl.enc_len_batch_list[p], mbl.dec_batch_list[p]])
    end
    return r
end

function convert_to_tree(r_list, i_left, i_right, form_manager)
    t = Tree()
    level = 0
    left = -1
    for i in range(i_left,stop=i_right)
        if r_list[i] ==  get_symbol_idx(form_manager,'(')
            if level == 0
                left = i
            end
            level = level + 1
        elseif r_list[i] == get_symbol_idx(form_manager ,')')
            level = level -1
            if level == 0
                if i == left+1
                    c = r_list[i]
                else
                    c = convert_to_tree(r_list, left + 1, i -1, form_manager)
                end
                add_child(t, c)
            end
        elseif level == 0
            add_child(t, r_list[i])

        end
    end
    return t
end

function norm_tree(r_list, form_manager)
    q = [convert_to_tree(r_list, 1, length(r_list), form_manager)]
    head = 1
    while head < length(q)
        t = q[head]
        if (t.children[1] == get_symbol_idx(form_manager, "and")) || (t.children[1] == get_symbol_idx(form_manager, "or"))
            # sort the following subchildren
            k = []
            for i in range(2, stop=length(t.children))
                if t.children[i] isa Tree
                    push!(k, (string.(t.children[i]), i))
                else
                    #k.append((str(t.children[i]), i))
                    push!(k, (string.(t.children[i]), i))
                end
            end
            sorted_t_dict = []
            k.sort(key=itemgetter(0))

            for key1 in k
                push!(sorted_t_dict, t.children[key1[1]])
            end
            for i in range(1, stop=t.num_children)
                #print('i ', i)
                t.children[i] = sorted_t_dict[i]
            end
        end
        # add children to q
        for i in range(1, stop=length(t.children))
            if t.children[i] isa Tree

                push!(q, t.children[i])
            end
        end
        head = head + 1
    end
    return q[0]
end

function is_all_same(c1, c2)
    if length(c1) == length(c2)
        all_same = true
        for j in 1:(length(c1))
            if c1[j] != c2[j]
                all_same = False
                break
            end
        end
        return all_same
    else
        return false
    end
end

function compute_accuracy(candidate_list, reference_list)
    if length(candidate_list) != length(reference_list)
        #print("candidate list has length {}, reference list has length {}\n".format(len(candidate_list), len(reference_list)))
    end
    len_min = min(length(candidate_list), length(reference_list))
    c = 0
    for i in 1:(len_min)
        if is_all_same(candidate_list[i], reference_list[i])
            print("above was all same")
            c = c+1
        end
    end
    return c/float.(len_min)
end

function compute_tree_accuracy(candidate_list_, reference_list_, form_manager)
    candidate_list = []
    for i in 1:(length(candidate_list_))
        #print("candidate\n\n")
        #push!(candidate_list, to_list(norm_tree(candidate_list_[i], form_manager), form_manager))
        push!(candidate_list, to_list(norm_tree(candidate_list_[i], form_manager), form_manager))
    end
    reference_list = []
    for i in 1:(length(reference_list_))
        #reference_list.append(norm_tree(reference_list_[i], form_manager).to_list(form_manager))
        push!(reference_list, to_list(norm_tree(reference_list_[i], form_manager), form_manager))
    end
    return compute_accuracy(candidate_list, reference_list)
end


text:
module SelectTP

using Serialization, LinearAlgebra
using DungBase
using AbstractPlotting, GLMakie, MakieLayout, Colors

export selectTP

_cmap(t) = range(HSL(colorant"green"), stop = HSL(colorant"red"), length = length(t))

function selectTP(file = "data")
    data = deserialize(file)
    scene, layout = layoutscene()
    video = Node("a")
    # layout[1, 1] = LText(scene, video)#, tellheight = true)
    ax = layout[1, 1] = LAxis(scene, aspect = DataAspect(), title = video)
    track = Node(rand(DungBase.Point, 10))
    rawtrack = Node(rand(DungBase.Point, 10))
    feeder = Node(rand(DungBase.Point, 1))
    lines!(ax, track, color = map(_cmap, track))
    scatter!(ax, rawtrack, markersize = 6px, color = map(_cmap, rawtrack))
    scatter!(ax, feeder, color = :blue, markersize = 8px)
    tp1 = Node(rand(2))
    mousestate = addmousestate!(ax.scene)
    onmouseover(mousestate) do state
        tp1[] = state.pos
    end
    tp = map(tp1) do xy
        _, i = findmin(norm.(track[] .- Ref(xy)))
        i
    end
    tp2 = map(tp) do i
        track[][i:i]
    end
    scatter!(ax, tp2, color = RGBA(0,0,0,0.5), markersize = 10px)
    onmouseleftclick(mousestate) do state
        println(string(video[], ",", tp[]))
    end
    c = Condition()
    h = layout[2,1] = LButton(scene, label = "Next!", tellwidth = false)
    on(h.clicks) do _
        notify(c)
    end
    for (k, v) in data
        @info "doing experiment: $k"
        for r in v.runs
            track[] = r.data.track.coords
            rawtrack[] = r.data.track.rawcoords.xy
            feeder[] = [r.data.feeder]
            video[] = r.metadata.comment
            display(scene)
            autolimits!(ax)
            wait(c)
        end
    end
    GLMakie.destroy!(display(scene))
    @info "Finished!!!"
end

end

ext: jl


text:
 addprocs(3)
@everywhere using Base.Test, Compat, IntervalArithmetic, ValidatedTransfer

# module Gauss
#   using IntervalArithmetic, ValidatedTransfer, Compat
#   # println(names(IntervalArithmetic))
#   # println(names(ValidatedTransfer))
#
#   λ = @biginterval 2
#   λ̌ = sqrt(λ)
#   C1 = @biginterval log(2)
#   δ = @biginterval acosh(16/5)
#   P = @biginterval acosh(2log2(1+2^1.1)-1)
#   H = @biginterval 2
#
#   sp = ChebyshevZeroOneBV(BigInterval)
#   entrybounds = AnalyticEntryBounds(sp,δ,P,H)
#   lybounds = FirstOrderLYBounds(sp,λ,C1)
#
#
#   # const biglog2 = (@biginterval(log(2)))::BigInterval
#   Compat.@compat immutable g end; Compat.@compat (::g){T}(x::T) = (1/(exp2(x)-1))::T
#   Compat.@compat immutable dg end; Compat.@compat (::dg){T}(x::T) = begin e2x = exp2(x); e2x1 = e2x-1; (-log(T(2))*e2x/(e2x1*e2x1))::T end
#   Compat.@compat immutable γa end; Compat.@compat (::γa){T}(x::T) = (log2(1+1/x))::T
#   Compat.@compat immutable h end; Compat.@compat (::h){T}(x::T) = (exp2(x)-1)::T
#   Compat.@compat immutable dh end; Compat.@compat (::dh){T}(x::T) = (log(T(2))*exp2(x))::T
#   α = IntervalArithmetic.@biginterval 1;
#   α1 = IntervalArithmetic.@biginterval -1;
#   G0 = IntervalArithmetic.@biginterval 1/log(2);
#   G1 = IntervalArithmetic.@biginterval 0.1
#   saferad = IntervalArithmetic.@biginterval 0.4;

# LANFORD TEST

@everywhere @eval begin
  N = 100
  N_interp = nextpow2(N)

  prec = 8round(Int,(N+512)*0.2/8)
  @compat setprecision(BigFloat,prec)

  import Base.sinpi
  Base.sinpi(x::IntervalArithmetic.Interval) = sin(pi*x)

  ## MAP CONSTANTS
  @compat lan_v(x,b) = (5 - sqrt.(25-8(x+b)))/2
  @compat lan_dv(x,b) = 2./sqrt.(25-8(x+b))
  λ = @biginterval 3/2
  λ̌ = sqrt(λ)
  C1 = @biginterval 4/9
  δ = @biginterval acosh( 7/4)
  P = @biginterval acosh(4 - sqrt(6))
  H = @biginterval sqrt((7+sqrt(33))/2) #5/2sqrt(39-24sqrt(2)) + 3/2sqrt(40sqrt(2)-55)

  # Pointwise aliasing error
  t(j) = (1+(j!=0))
  aliaserror(j,k,N) = (H * t(j-1) * exp(P*(k-1)) * 2cosh(δ*(j-1)) /(exp(2δ*N)-1)).hi
  entrybound(j,k) = (H*t(j)*exp((k-1)*P-(j-1)*δ)).hi

  @compat @assert precision(BigFloat) == prec

  function acos_lan_v_norm{T}(A::Vector{T},b)
    B = similar(A)
    for i in eachindex(B)
      B[i] = acos(2lan_v((A[i]+1)/2,b)-1)
    end
    B
  end

  function lan_dv_norm{T}(A::Vector{T},b)
    B = similar(A)
    for i in eachindex(B)
      B[i] = lan_dv((A[i]+1)/2,b)
    end
    B
  end

  cheb_pts = ValidatedTransfer.chebyshevpoints(N_interp,BigInterval)
  acos_v0_cheb_pts = acos_lan_v_norm(cheb_pts,0); acos_v1_cheb_pts = acos_lan_v_norm(cheb_pts,1);
  dv0_cheb_pts = lan_dv_norm(cheb_pts,0); dv1_cheb_pts = lan_dv_norm(cheb_pts,1);
  helperdata = [acos_v0_cheb_pts,dv0_cheb_pts,acos_v1_cheb_pts,dv1_cheb_pts]

  println(size(acos_v0_cheb_pts))
  lanford_transferfn(i,k,helperdata) = cos((k-1)*helperdata[1][i])*helperdata[2][i] + cos((k-1)*helperdata[3][i])*helperdata[4][i]
end

@time M = maketransfer(N,lanford_transferfn,helperdata,aliaserror,entrybound,workers(),N_interp)
@time S = getsolution(M,workers())
@time savezip(tempname(),S)


text:
using Pkg

# Activate test environment on older Julia versions
@static if VERSION < v"1.2"
    Pkg.activate(@__DIR__)
    Pkg.develop(PackageSpec(; path=dirname(@__DIR__)))
    Pkg.instantiate()
end

using MCMCDiagnosticTools
using FFTW

using Random
using Statistics
using Test

Random.seed!(1)

@testset "MCMCDiagnosticTools.jl" begin
    @testset "discrete diagnostic" begin
        include("discretediag.jl")
    end
    @testset "ESS" begin
        include("ess.jl")
    end
    @testset "Gelman, Rubin and Brooks diagnostic" begin
        include("gelmandiag.jl")
    end
    @testset "Geweke diagnostic" begin
        include("gewekediag.jl")
    end
    @testset "Heidelberger and Welch diagnostic" begin
        include("heideldiag.jl")
    end
    @testset "Monte Carlo standard error" begin
        include("mcse.jl")
    end
    @testset "Raftery and Lewis diagnostic" begin
        include("rafterydiag.jl")
    end
    @testset "R⋆ diagnostic" begin
        # MLJXGBoostInterface requires Julia >= 1.3
        # XGBoost errors on 32bit systems: https://github.com/dmlc/XGBoost.jl/issues/92
        if VERSION >= v"1.3" && Sys.WORD_SIZE == 64
            # run tests related to rstar statistic
            Pkg.activate("rstar")
            Pkg.develop(; path=dirname(dirname(pathof(MCMCDiagnosticTools))))
            Pkg.instantiate()
            include(joinpath("rstar", "runtests.jl"))
        else
            @info "R⋆ not tested: requires Julia >= 1.3 and a 64bit architecture"
        end
    end
end


text:
using Flux: train!, params
using IterTools:ncycle
using Statistics:mean
using LinearAlgebra:norm

"""
    ESADDetector(encoder = Chain(),
                decoder = Chain(),
                batchsize = 32,
                epochs = 1,
                shuffle = false,
                partial = true,
                opt = ADAM(),
                λ1 = 1,
                λ2 = 1,
                noise = identity)

End-to-End semi-supervised anomaly detection algorithm similar to DeepSAD, but without the pretraining phase. The
algorithm was published by Huang et al., see [1].

Parameters
----------
$AE_LOSS

    λ1::Real
Weighting parameter of the norm loss, which minimizes the empirical variance and thus minimizes entropy.

    λ2::Real
Weighting parameter of the assistent loss function to define the consistency between the two encoders.

    noise::Function (AbstractArray{T} -> AbstractArray{T})
A function to be applied to a batch of input data to add noise, see [1] for an explanation.

Examples
--------
$(SCORE_SUPERVISED("ESADDetector"))

References
----------
[1] Huang, Chaoqin; Ye, Fei; Zhang, Ya; Wang, Yan-Feng; Tian, Qi (2020): ESAD: End-to-end Deep Semi-supervised Anomaly
Detection.
"""
OD.@detector mutable struct ESADDetector <: SupervisedDetector
    encoder::Chain = Chain()
    decoder::Chain = Chain()
    batchsize::Integer = 32::(_ > 0)
    epochs::Integer = 1::(_ > 0)
    shuffle::Bool = false
    partial::Bool = true
    opt::Any = ADAM()
    λ1::Number = 1
    λ2::Number = 1
    noise::Function = identity
end

struct ESADModel <: DetectorModel
    chain::Chain
end

function OD.fit(detector::ESADDetector, X::Data, y::Labels; verbosity)::Fit
    loader = DataLoader((X, y), batchsize=detector.batchsize, shuffle=detector.shuffle, partial=detector.partial)

    # Create the autoencoder // TODO: deepcopy the encoder/decoder?
    model = Chain(detector.encoder, detector.decoder, detector.encoder)

    # Precalculate dimensions
    dims = ndims(X)

    # Calculate loss as described in the paper
    train!((x, y) -> _esadloss(model[1:2](x), x, model(x), model[1](x), y, detector.λ1, detector.λ2, detector.noise,
    dims), params(model), ncycle(loader, detector.epochs), detector.opt)

    # Score as described in the paper
    scores = _esadscore(model[1:2](X), X, model(X), dims)
    return ESADModel(model), scores
end

function OD.transform(_::ESADDetector, model::ESADModel, X::Data)::Scores
    _esadscore(model.chain[1:2](X), X, model.chain(X), ndims(X))
end

function _esadloss(x̂, x, ẑ, z, y, λ1, λ2, noise, dims)
    # The esad loss function is based on the distance to the hypersphere center. The inverse distance is used if an
    # example is an outlier and labeled samples are weighted using the hyperparameter eta
    rec_loss(y, rec_outlier, rec_normal) = ifelse(y == CLASS_OUTLIER, rec_outlier, rec_normal)
    norm_loss(y, origin_dist) = ifelse(y == CLASS_OUTLIER, origin_dist ^ -1, origin_dist)

    rec_outlier = dropdims(sum((x̂ .- noise(x)).^2, dims=1:dims - 1), dims=1)
    rec_normal = dropdims(sum((x̂ .- x).^2, dims=1:dims - 1), dims=1)
    origin_dist = map(norm, eachslice(ẑ; dims=ndims(ẑ)))

    # reconstruction error (optimize mutual information of latent)
    l_rec_semi = rec_loss.(y, rec_outlier, rec_normal) |> mean

    # minimize distance to origin (thus minimize entropy for normal and maximize entropy for outliers)
    l_norm_semi = norm_loss.(y, origin_dist) |> mean

    # similarity constraint of latent embeddings
    l_ass_semi = mse(ẑ, z)

    # combine the (weighted) loss terms (scalars)
    l_rec_semi + λ1 * l_norm_semi + λ2 * l_ass_semi
end

function _esadscore(x̂, x, ẑ, dims)
    # Combines reconstruction loss and distance to the origin
    origin_dist = map(norm, eachslice(ẑ; dims=ndims(ẑ)))
    rec_loss = dropdims(mean((x̂ .- x).^2, dims=1:dims - 1), dims=1)
    rec_loss .+ origin_dist
end

max_stars_repo_name: davnn/OutlierDetectionNetworks

max_issues_repo_name: davnn/OutlierDetectionNetworks

max_forks_repo_name: davnn/OutlierDetectionNetworks


text:
### gls.jl ---  Perform Generalised Lomb–Scargle periodogram
#
# Copyright (C) 2017 Mosè Giordano.
#
# Maintainer: Mosè Giordano <mose AT gnu DOT org>
# Keywords: periodogram, lomb scargle
#
# This file is a part of LombScargle.jl.
#
# License is BSD 3-clause "New" or "Revised".
#
### Code:

struct GLSPlan{T,A,B<:AbstractVector{T},C,D,E,F,G} <: PeriodogramPlan
    times::A
    signal::B
    freq::C
    sumw::D
    w::E
    y::B
    YY::T
    noise::F
    norm::Symbol
    P::G
end

struct GLSPlan_fit_mean{T,A,B<:AbstractVector{T},C,D,E,F,G,H} <: PeriodogramPlan
    times::A
    signal::B
    freq::C
    sumw::D
    w::E
    y::B
    Y::F
    YY::T
    noise::G
    norm::Symbol
    P::H
end

# Generalised Lomb–Scargle algorithm: this takes into account uncertainties and
# fit the mean of the signal.  This is implemented following the recipe by
# * Zechmeister, M., Kürster, M. 2009, A&A, 496, 577  (URL:
#   http://dx.doi.org/10.1051/0004-6361:200811296,
#   Bibcode: http://adsabs.harvard.edu/abs/2009A%26A...496..577Z)
# In addition, some tricks suggested by
# * Press, W. H., Rybicki, G. B. 1989, ApJ, 338, 277 (URL:
#   http://dx.doi.org/10.1086/167197,
#   Bibcode: http://adsabs.harvard.edu/abs/1989ApJ...338..277P)
# to make computation faster are adopted.

function _generalised_lombscargle!(P, freqs, times, y, w, Y, YY, nil)
    @inbounds Threads.@threads for n in eachindex(freqs)
        ω = freqs[n] * 2 * pi
        # Find τ for current angular frequency
        C = S = CS = CC = nil
        @inbounds for i in eachindex(times)
            W    = w[i]
            s, c = sincos(ω * times[i])
            CS  += W*c*s
            CC  += W*c*c
            C   += W*c
            S   += W*s
        end
        CS -= C*S
        SS  = 1 - CC - S*S
        CC -= C*C
        ωτ   = atan(2CS, CC - SS) / 2
        # Now we can compute the power
        YC_τ = YS_τ = CC_τ = nil
        @inbounds for i in eachindex(times)
            W     = w[i]
            s, c  = sincos(ω*times[i] - ωτ)
            YC_τ += W*y[i]*c
            YS_τ += W*y[i]*s
            CC_τ += W*c*c
        end
        # "C_τ" and "S_τ" are computed following equation (7) of Press &
        # Rybicki, this formula simply comes from angle difference trigonometric
        # identities.
        sin_ωτ, cos_ωτ = sincos(ωτ)
        C_τ    = C*cos_ωτ + S*sin_ωτ
        S_τ    = S*cos_ωτ - C*sin_ωτ
        YC_τ  -= Y*C_τ
        YS_τ  -= Y*S_τ
        SS_τ   = 1 - CC_τ - S_τ*S_τ
        CC_τ  -= C_τ*C_τ
        P[n] = (abs2(YC_τ)/CC_τ + abs2(YS_τ)/SS_τ)/YY
    end
    return P
end

function _generalised_lombscargle!(P, freqs, times, y, w, YY, nil)
    @inbounds Threads.@threads for n in eachindex(freqs)
        ω = freqs[n] * 2 * pi
        # Find τ for current angular frequency
        C = S = CS = CC = nil
        @inbounds for i in eachindex(times)
            W    = w[i]
            s, c = sincos(ω * times[i])
            CS  += W*c*s
            CC  += W*c*c
        end
        ωτ   = atan(2CS, 2CC - 1) / 2
        # Now we can compute the power
        YC_τ = YS_τ = CC_τ = nil
        @inbounds for i in eachindex(times)
            W     = w[i]
            s, c  = sincos(ω*times[i] - ωτ)
            YC_τ += W*y[i]*c
            YS_τ += W*y[i]*s
            CC_τ += W*c*c
        end
        # P[n] should be (abs2(YC_τ)/CC_τ + abs2(YS_τ)/SS_τ)/YY, but we guard the values of
        # CC_τ and SS_τ in case they're zeros (note that SS_τ = 1 - CC_τ).
        frac_C = ifelse(iszero(CC_τ), nil, abs2(YC_τ) / CC_τ)
        frac_S = ifelse(CC_τ == 1,    nil, abs2(YS_τ) / (1 - CC_τ))
        P[n] = (frac_C + frac_S)/YY
    end
    return P
end

_periodogram!(P::AbstractVector, times, p::GLSPlan_fit_mean) =
    _generalised_lombscargle!(P, p.freq, times, p.y, p.w, p.Y, p.YY, zero(eltype(p.P)))
_periodogram!(p::GLSPlan_fit_mean) = _periodogram!(p.P, p.times, p)
_periodogram!(P::AbstractVector, times, p::GLSPlan) =
    _generalised_lombscargle!(P, p.freq, times, p.y, p.w, p.YY, zero(eltype(p.P)))
_periodogram!(p::GLSPlan) = _periodogram!(p.P, p.times, p)


text:
# `_vscalef` for architectures without `vscalef`.
# magic rounding constant: 1.5*2^52 Adding, then subtracting it from a float rounds it to an Int.
MAGIC_ROUND_CONST(::Type{Float64}) = 6.755399441055744e15
MAGIC_ROUND_CONST(::Type{Float32}) = 1.2582912f7
@inline function vscalef(x::Union{T,AbstractSIMD{<:Any,T}}, y::Union{T,AbstractSIMD{<:Any,T}}, ::False) where {T<:Union{Float32,Float64}}
    _vscalef(x, floor(y))
end
@inline signif_bits(::Type{Float32}) = 0x00000017 # 23
@inline signif_bits(::Type{Float64}) = 0x0000000000000034 # 52
@inline function _vscalef(x::Union{T,AbstractSIMD{<:Any,T}}, y::Union{T,AbstractSIMD{<:Any,T}}) where {T<:Union{Float32,Float64}}
    N = reinterpret(Base.uinttype(T), y + MAGIC_ROUND_CONST(T))
    k = N# >>> 0x00000008

    small_part = reinterpret(Base.uinttype(T), x)
    twopk = (k % Base.uinttype(T)) << signif_bits(T)
    reinterpret(T, twopk + small_part)
end
@inline vscalef(m::AbstractMask, v1::AbstractSIMD, v2::AbstractSIMD, v3::AbstractSIMD, ::False) = vifelse(m, vscalef(v1, v2, False()), v3)
@inline vscalef(v1::T,v2::T) where {T<:AbstractSIMD} = vscalef(v1,v2,has_feature(Val(:x86_64_avx512f)))
@inline vscalef(m::AbstractMask, v1::T, v2::T, v3::T) where {T<:AbstractSIMD} = vscalef(m,v1,v2,v3,has_feature(Val(:x86_64_avx512f)))
@inline vscalef(v1::T,v2::T) where {T<:Union{Float32,Float64}} = vscalef(v1,v2,False())
@inline vscalef(b::Bool,v1::T,v2::T,v3::T) where {T<:Union{Float32,Float64}} = b ? vscalef(v1,v2,False()) : zero(T)
@inline vscalef(v1,v2) = ((v3,v4) = promote(v1,v2); vscalef(v3,v4))
@generated function vscalef(v1::Vec{W,T}, v2::Vec{W,T}, ::True) where {W,T<:Union{Float32,Float64}}
    bits = 8W*sizeof(T)
    bits ∈ (128,256,512) || throw(ArgumentError("Vectors are $bits bits, but only 128, 256, and 512 bits are supported."))
    ltyp = LLVM_TYPES[T]
    vtyp = "<$W x $ltyp>"
    dors = T === Float64 ? 'd' : 's'
    instr = "$vtyp @llvm.x86.avx512.mask.scalef.p$(dors).$bits"
    mtyp = W == 16 ? "i16" : "i8"
    if bits == 512
        decl = "declare $instr($vtyp, $vtyp, $vtyp, $mtyp, i32)"
        # instrs = "%res = call $instr($vtyp %0, $vtyp %1, $vtyp undef, $mtyp -1, i32 11)\nret $vtyp %res"
        instrs = "%res = call $instr($vtyp %0, $vtyp %1, $vtyp undef, $mtyp -1, i32 8)\nret $vtyp %res"
    else
        decl = "declare $instr($vtyp, $vtyp, $vtyp, $mtyp)"
        instrs = "%res = call $instr($vtyp %0, $vtyp %1, $vtyp undef, $mtyp -1)\nret $vtyp %res"
    end
    arg_syms = [:(data(v1)), :(data(v2))]
    llvmcall_expr(decl, instrs, :(_Vec{$W,$T}), :(Tuple{_Vec{$W,$T},_Vec{$W,$T}}), vtyp, fill(vtyp, 2), arg_syms)
end
@generated function vscalef(m::AbstractMask{W}, v1::Vec{W,T}, v2::Vec{W,T}, v3::Vec{W,T}, ::True) where {W,T<:Union{Float32,Float64}}
    bits = 8W*sizeof(T)
    bits ∈ (128,256,512) || throw(ArgumentError("Vectors are $bits bits, but only 128, 256, and 512 bits are supported."))
    ltyp = LLVM_TYPES[T]
    vtyp = "<$W x $ltyp>"
    dors = T === Float64 ? 'd' : 's'
    mtyp = W == 16 ? "i16" : "i8"
    mtypj = W == 16 ? :UInt16 : :UInt8
    instr = "$vtyp @llvm.x86.avx512.mask.scalef.p$(dors).$bits"
    if bits == 512
        decl = "declare $instr($vtyp, $vtyp, $vtyp, $mtyp, i32)"
        instrs = "%res = call $instr($vtyp %0, $vtyp %1, $vtyp %2, $mtyp %3, i32 11)\nret $vtyp %res"
        instrs = "%res = call $instr($vtyp %0, $vtyp %1, $vtyp %2, $mtyp %3, i32 8)\nret $vtyp %res"
    else
        decl = "declare $instr($vtyp, $vtyp, $vtyp, $mtyp)"
        instrs = "%res = call $instr($vtyp %0, $vtyp %1, $vtyp %2, $mtyp %3)\nret $vtyp %res"
    end
    arg_syms = [:(data(v1)), :(data(v2)), :(data(v3)), :(data(m))]
    llvmcall_expr(decl, instrs, :(_Vec{$W,$T}), :(Tuple{_Vec{$W,$T},_Vec{$W,$T},_Vec{$W,$T},$mtypj}), vtyp, [vtyp,vtyp,vtyp,mtyp], arg_syms)
end
@generated function vsreduce(v::Vec{W,T}, ::Val{M}) where {W,T<:Union{Float32,Float64}, M}
    bits = 8W*sizeof(T)
    bits ∈ (128,256,512) || throw(ArgumentError("Vectors are $bits bits, but only 128, 256, and 512 bits are supported."))
    M isa Integer || throw(ArgumentError("M must be an integer, but received $M of type $(typeof(M))."))
    ltyp = LLVM_TYPES[T]
    vtyp = "<$W x $ltyp>"
    dors = T === Float64 ? 'd' : 's'
    instr = "$vtyp @llvm.x86.avx512.mask.reduce.p$(dors).$bits"
    mtyp = W == 16 ? "i16" : "i8"
    if bits == 512
        decl = "declare $instr($vtyp, i32, $vtyp, $mtyp, i32)"
        instrs = "%res = call $instr($vtyp %0, i32 $M, $vtyp undef, $mtyp -1, i32 8)\nret $vtyp %res"
    else
        decl = "declare $instr($vtyp, i32, $vtyp, $mtyp)"
        instrs = "%res = call $instr($vtyp %0, i32 $M, $vtyp undef, $mtyp -1)\nret $vtyp %res"
    end
    arg_syms = [:(data(v))]
    llvmcall_expr(decl, instrs, :(_Vec{$W,$T}), :(Tuple{_Vec{$W,$T}}), vtyp, [vtyp], arg_syms)
end

@generated function vpermi2pd(c::Vec{8,UInt64}, v1::Vec{8,Float64}, v2::Vec{8,Float64}) #where {W,T<:Union{Float32,Float64}, M}
    W = 8; T = Float64
    bits = 8W*sizeof(T)
    # bits ∈ (128,256,512) || throw(ArgumentError("Vectors are $bits bits, but only 128, 256, and 512 bits are supported."))
    ityp = "i$(8sizeof(T))"
    vityp = "<$W x $ityp>"
    ltyp = LLVM_TYPES[T]
    vtyp = "<$W x $ltyp>"
    dors = T === Float64 ? 'd' : 's'
    instr = "$vtyp @llvm.x86.avx512.vpermi2var.p$(dors).$bits"
    decl = "declare $instr($vtyp, $vityp, $vtyp)"
    instrs = "%res = call $instr($vtyp %0, $vityp %1, $vtyp %2)\nret $vtyp %res"
    arg_syms = [:(data(v1)), :(data(c)), :(data(v2))]
    jityp = T === Float64 ? :UInt64 : :UInt32
    llvmcall_expr(decl, instrs, :(_Vec{$W,$T}), :(Tuple{_Vec{$W,$T},_Vec{$W,$jityp},_Vec{$W,$T}}), vtyp, [vtyp,vityp,vtyp], arg_syms)
end

@inline vscalef(v1::VecUnroll, v2::VecUnroll) = VecUnroll(fmap(vscalef, getfield(v1, :data), getfield(v2, :data)))
@inline vscalef(m::VecUnroll, v1::VecUnroll, v2::VecUnroll, v3::VecUnroll) = VecUnroll(fmap(vscalef, getfield(m, :data), getfield(v1, :data), getfield(v2, :data), getfield(v3, :data)))
@inline vsreduce(v::VecUnroll, ::Val{M}) where {M} = VecUnroll(fmap(vsreduce, getfield(v, :data), Val{M}()))
@inline vpermi2pd(v1::VecUnroll, v2::VecUnroll, v3::VecUnroll) = VecUnroll(fmap(vpermi2pd, getfield(v1, :data), getfield(v2, :data), getfield(v3, :data)))
@inline vpermi2pd(v1::VecUnroll, v2::Vec, v3::Vec) = VecUnroll(fmap(vpermi2pd, getfield(v1, :data), v2, v3))

# magic rounding constant: 1.5*2^52 Adding, then subtracting it from a float rounds it to an Int.

# min and max arguments by base and type
MAX_EXP(::Val{2},  ::Type{Float64}) =  1024.0                   # log2 2^1023*(2-2^-52)
MIN_EXP(::Val{2},  ::Type{Float64}) = -1022.0                   # log2(big(2)^-1023*(2-2^-52))
MAX_EXP(::Val{2},  ::Type{Float32}) =  128f0                    # log2 2^127*(2-2^-52)
MIN_EXP(::Val{2},  ::Type{Float32}) = -126f0                    # log2 2^-1075
MAX_EXP(::Val{ℯ},  ::Type{Float64}) =  709.782712893383996732   # log 2^1023*(2-2^-52)
MIN_EXP(::Val{ℯ},  ::Type{Float64}) = -708.396418532264106335  # log 2^-1075
MAX_EXP(::Val{ℯ},  ::Type{Float32}) =  88.72283905206835f0      # log 2^127 *(2-2^-23)
MIN_EXP(::Val{ℯ},  ::Type{Float32}) = -87.3365448101577555f0#-103.97207708f0           # log 2^-150
MAX_EXP(::Val{10}, ::Type{Float64}) =  308.25471555991675       # log10 2^1023*(2-2^-52)
MIN_EXP(::Val{10}, ::Type{Float64}) = -307.65260000       # log10 2^-1075
MAX_EXP(::Val{10}, ::Type{Float32}) =  38.531839419103626f0     # log10 2^127 *(2-2^-23)
MIN_EXP(::Val{10}, ::Type{Float32}) = -37.9297794795476f0      # log10 2^-127 *(2-2^-23)

# 256/log(base, 2) (For Float64 reductions)
LogBo256INV(::Val{2}, ::Type{Float64})    = 256.
LogBo256INV(::Val{ℯ}, ::Type{Float64})    = 369.3299304675746
LogBo256INV(::Val{10}, ::Type{Float64})   = 850.4135922911647
LogBo16INV(::Val{2}, ::Type{Float64})    = 16.
LogBo16INV(::Val{ℯ}, ::Type{Float64})    = 23.083120654223414
LogBo16INV(::Val{10}, ::Type{Float64})   = 53.150849518197795
# -log(base, 2)/256 in upper and lower bits
LogBo256U(::Val{2}, ::Type{Float64})      = -0.00390625
LogBo256U(::Val{ℯ}, ::Type{Float64})      = -0.0027076061740622863
LogBo256U(::Val{10}, ::Type{Float64})     = -0.0011758984205624266
LogBo256L(base::Val{2}, ::Type{Float64})  =  0.
LogBo256L(base::Val{ℯ}, ::Type{Float64})  = -9.058776616587108e-20
LogBo256L(base::Val{10}, ::Type{Float64}) = 1.0952062999160822e-20

LogBo16U(::Val{2}, ::Type{Float64})      = -0.0625
LogBo16U(::Val{ℯ}, ::Type{Float64})      = -0.04332169878499658183857700759113603550471875839751595338254250059333710137310597
LogBo16U(::Val{10}, ::Type{Float64})     = -0.01881437472899882470085868092028081417301186759138178383190171632044426182965149
LogBo16L(base::Val{2}, ::Type{Float64})  = Zero()
LogBo16L(base::Val{ℯ}, ::Type{Float64})  = -1.4494042586539372e-18
LogBo16L(base::Val{10}, ::Type{Float64}) = 1.7523300798657315e-19

# 1/log(base, 2) (For Float32 reductions)
LogBINV(::Val{2}, ::Type{Float32})    =  1f0
LogBINV(::Val{ℯ}, ::Type{Float32})    =  1.442695f0
LogBINV(::Val{10}, ::Type{Float32})   =  3.321928f0
# -log(base, 2) in upper and lower bits
LogBU(::Val{2}, ::Type{Float32})      = -1f0
LogBU(::Val{ℯ}, ::Type{Float32})      = -0.6931472f0
LogBU(::Val{10}, ::Type{Float32})     = -0.30103f0
LogBL(base::Val{2}, ::Type{Float32})  =  0f0
LogBL(base::Val{ℯ}, ::Type{Float32})  =  1.9046542f-9
LogBL(base::Val{10}, ::Type{Float32}) =  1.4320989f-8

const FloatType64 = Union{Float64,AbstractSIMD{<:Any,Float64}}
const FloatType32 = Union{Float32,AbstractSIMD{<:Any,Float32}}
# Range reduced kernels
@inline function expm1b_kernel(::Val{2}, x::FloatType64)
    # c6 = 0.6931471807284470571335252997834339128744539291358546258980326560263434831636494
    # c5 = 0.2402265119815758621794630410361025063296309075509484445159065872903725193960909
    # c4 = 0.05550410353447979823044149277158612685395896745775325243210607075766620053156177
    # c3 = 0.009618027253668450057706478612143223628979891379942570690446533010539871321541621
    # c2 = 0.001333392256353875413926876917141786686018234585146983223440727245459444740967253
    # c1 = 0.0001546929114168849728971327603158937595919966441732209337930866845915899223829891
    # c0 = 1.520192159457321441849564286267892034534060236471603225598783028117591315796835e-05
    # c5 = 0.6931472067096466099497350107329038640311532915014328403152023514862215769240471
    # c4 = 0.2402265150505520926831534797602254284855354178135282005410994764797061783074115
    # c3 = 0.05550327215766594452554739168596479012109775178907146059734348050730310383358696
    # c2 = 0.00961799451416147891836707565892019722415069604010430702590438969636015825337285
    # c1 = 0.001340043166700788064581996335332499076913713747844545061461125917537819216640726
    # c0 = 0.0001547802227945780278842074081393459334029013917349238368250485753166922523500416
    # x * muladd(muladd(muladd(muladd(muladd(muladd(muladd(muladd(c0,x,c1),x,c2),x,c3),x,c4),x,c5),x,c6),x,c7),x,c8)
    # x * muladd(muladd(muladd(muladd(muladd(muladd(c0,x,c1),x,c2),x,c3),x,c4),x,c5),x,c6)
    # x * muladd(muladd(muladd(muladd(muladd(c0,x,c1),x,c2),x,c3),x,c4),x,c5)
  x * muladd(muladd(muladd(0.009618130135925114, x, 0.055504115022757844), x, 0.2402265069590989), x, 0.6931471805599393)
end
@inline function expm1b_kernel(::Val{ℯ}, x::FloatType64)
  x * muladd(muladd(muladd(0.04166666762124105, x, 0.1666666704849642), x, 0.49999999999999983), x, 0.9999999999999998)
end
@inline function expm1b_kernel(::Val{10}, x::FloatType64)
  x * muladd(muladd(muladd(muladd(0.5393833837413015, x, 1.1712561359457612), x, 2.0346785922926713), x, 2.6509490552382577), x, 2.302585092994046)
end
@inline function expb_kernel(::Val{2}, x::FloatType32)
  muladd(muladd(muladd(muladd(muladd(muladd(muladd(1.5316464f-5, x, 0.00015478022f0), x, 0.0013400431f0), x, 0.009617995f0), x, 0.05550327f0), x, 0.24022652f0), x, 0.6931472f0), x, 1.0f0)
end
@inline function expb_kernel(::Val{ℯ}, x::FloatType32)
  muladd(muladd(muladd(muladd(muladd(muladd(muladd(0.00019924171f0, x, 0.0013956056f0), x, 0.008375129f0), x, 0.041666083f0), x, 0.16666415f0), x, 0.5f0), x, 1.0f0), x, 1.0f0)
end
@inline function expb_kernel(::Val{10}, x::FloatType32)
  muladd(muladd(muladd(muladd(muladd(muladd(muladd(0.06837386f0, x, 0.20799689f0), x, 0.54208815f0), x, 1.1712388f0), x, 2.034648f0), x, 2.6509492f0), x, 2.3025851f0), x, 1.0f0)
end

const J_TABLE= Float64[2.0^(big(j-1)/256) for j in 1:256];

@inline fast_fma(a, b, c, ::True) = fma(a, b, c)
@inline function fast_fma(a, b, c, ::False)
    d = dadd(dmul(Double(a),Double(b),False()),Double(c))
    add_ieee(d.hi, d.lo)
end

@static if (Sys.ARCH === :x86_64) | (Sys.ARCH === :i686)
const TABLE_EXP_64_0= Vec(ntuple(j -> Core.VecElement(Float64(2.0^(big(j-1)/16))), Val(8)))
const TABLE_EXP_64_1= Vec(ntuple(j -> Core.VecElement(Float64(2.0^(big(j+7)/16))), Val(8)))

@inline target_trunc(v, ::VectorizationBase.True) = v
@inline target_trunc(v, ::VectorizationBase.False) = v % UInt32
@inline target_trunc(v) = target_trunc(v, VectorizationBase.has_feature(Val(:x86_64_avx512dq)))

# @inline function vexp2_v1(x::AbstractSIMD{8,Float64})
#     x16 = x
#     # x16 = 16x
#     r =  vsreduce(x16, Val(4))
#     m = x16 - r
#     mfrac = m
#     inds = (reinterpret(UInt64, mfrac) >> 0x000000000000002d) & 0x000000000000000f
#     # @show r m mfrac reinterpret(UInt64, m) reinterpret(UInt64, mfrac)
#     # js = vpermi2pd(inds, TABLE_EXP_64_0, TABLE_EXP_64_1)
#     # @show m mfrac r
#     small_part = expm1b_kernel(Val(2), r) + 1.0
#     # js = 1.0
#     # small_part = vfmadd(js, expm1b_kernel(Val(2), r), js)
#     vscalef(small_part, mfrac)
# end

@inline function expm1b_kernel_16(::Val{2}, x)
    c5 = 0.6931471805599457533351827593325319924753473772859614915719486459595837313933663
    c4 = 0.2402265069591009431089489060897837825648676480621950809556237945205562267112511
    c3 = 0.05550410865663929372911461843767669974316894963735870580154522796380984673567634
    c2 = 0.00961812910613182376367867689426991348318504321185094738294343470767871628697879
    c1 = 0.001333378157683735211078403326604752238340853209789619792858391909299167771871147
    c0 = 0.0001540378851029625623114826398060979330719673345637296642237670082377277446583639

    x * vmuladd_fast(vmuladd_fast(vmuladd_fast(vmuladd_fast(vmuladd_fast(c0,x,c1),x,c2),x,c3),x,c4),x,c5)
end
@inline function expm1b_kernel_16(::Val{ℯ}, x)
    c5 = 1.000000000000000640438225946852982701258391604480638275588427888399057119915227
    c4 = 0.5000000000000004803287542332217715766116071510522583538834274474935075866898906
    c3 = 0.1666666666420970554527207281834431226735741940161719645985080160507073865025253
    c2 = 0.04166666666018301921665935823120024420659933010915524936059730214858712998209511
    c1 = 0.008333472974984405879148046293292753978131598285368755170712233472964279745777974
    c0 = 0.001388912162496018623851591184688043066476532389093553363939521815388727152058955
    x * vmuladd_fast(vmuladd_fast(vmuladd_fast(vmuladd_fast(vmuladd_fast(c0,x,c1),x,c2),x,c3),x,c4),x,c5)
end
@inline function expm1b_kernel_16(::Val{10}, x)
    c5 = 2.302585092994047158681503503460480873860999793973827515869365761962430703319599
    c4 = 2.650949055239201551934947671858339219424413336755573194223955910456821842421188
    c3 = 2.034678591993528625083054018110717593492391962926107057414972651492780932119609
    c2 = 1.171255148730010832148986815739840479496404930758611797062000621828706728750671
    c1 = 0.5393919676343165962473794696403862709895176169091768718271240564783020989456929
    c0 = 0.2069993173257113377172910397724414085027323868592924170210484263853229417141011
    x * vmuladd_fast(vmuladd_fast(vmuladd_fast(vmuladd_fast(vmuladd_fast(c0,x,c1),x,c2),x,c3),x,c4),x,c5)
end

# @inline function expm1b_kernel_6(::Val{2}, x)
#     c6 = 0.6931471805599453094157857128856867906777808839530204388709017060018142940776171
#     c5 = 0.2402265069591008469523412357777710806331204676211882505524452303635329805655358
#     c4 = 0.05550410866482161698117481145997657675888637052903928922860089114193615876343007
#     c3 = 0.009618129106525681209446910869436511674612664593938506025921861852732097374219765
#     c2 = 0.001333355814485111028814491629997422287837080138262113640339762815224264878390602
#     c1 = 0.0001540375624549508734984308915404920724081959405205139643458055193131059150331185
#     c0 = 1.525295744513115862409484203574886089068776710280414076942740530839031792795809e-05
#     x * vmuladd_fast(vmuladd_fast(vmuladd_fast(vmuladd_fast(vmuladd_fast(vmuladd_fast(c0,x,c1),x,c2),x,c3),x,c4),x,c5),x,c6)
# end
# # using Remez
# # N,D,E,X = ratfn_minimax(x -> (exp2(x) - big"1")/x, [big(nextfloat(-0.03125)),big(0.03125)], 4, 0); @show(E); N
# @inline function expm1b_kernel_4(::Val{2}, x)
#     c4 = 0.6931471805599461972549081995383434692316977327912755704234013405443109498729026
#     c3 = 0.2402265069131940842333497738928958607054740795078596615709864445611497846077303
#     c2 = 0.05550410865300270379171299778517151376504051870524903806295523325435530249981495
#     c1 = 0.009618317140648284298097106744730186251913149278152053357630395863210686828434175
#     c0 = 0.001333381881551676348461495248002642715422207072457864472267417920610122672570108
#     x * vmuladd_fast(vmuladd_fast(vmuladd_fast(vmuladd_fast(c0,x,c1),x,c2),x,c3),x,c4)
# end

# @inline function vexp2_v4(x::AbstractSIMD{W,Float64}) where {W}
#     # r - VectorizationBase.vroundscale(r, Val(16*(4)))
#     r = VectorizationBase.vsreduce(x,Val(0))
#     rscale = VectorizationBase.vroundscale(r, Val(64))
#     rs = r - rscale
#     inds = convert(UInt, vsreduce(rscale, Val(1))*16.0)
#     expr = expm1b_kernel_5(Val(2), rs)
#     N_float = x - rs
#     # @show inds rs N_float

#     js = vpermi2pd(inds, TABLE_EXP_64_0, TABLE_EXP_64_1)
#     small_part = vfmadd(js, expr, js)
#     res = vscalef(small_part, N_float)
#     return res

# end

####################################################################################################
################################## Non-AVX512 implementation #######################################
####################################################################################################

# With AVX512, we use a tiny look-up table of just 16 numbers for `Float64`,
# because we can perform the lookup using `vpermi2pd`, which is much faster than gather.
# To compensate, we need a larger polynomial.

# Because of the larger polynomial, this implementation works better on systems with 2 FMA units.

@inline function vexp2(x::AbstractSIMD{8,Float64}, ::True)
    r = vsreduce(16.0x, Val(0)) * 0.0625
    N_float = x - r
    expr = expm1b_kernel_16(Val(2), r)
    inds = convert(UInt64, vsreduce(N_float, Val(1))*16.0)
    # inds = ((trunc(Int64, 16.0*N_float)%UInt64)) & 0x000000000000000f
    js = vpermi2pd(inds, TABLE_EXP_64_0, TABLE_EXP_64_1)
    small_part = vfmadd(js, expr, js)
    res = vscalef(small_part, N_float)
    return res
end
# @inline function vexp2_v2(x::AbstractSIMD{8,Float64}, ::True)#, ::Val{N}) where {N}
#     r1 = vsreduce(x, Val(0))
#     m = x - r1
#     r = vfmsub(vsreduce(r1 * 16.0, Val(1)), 0.0625, 0.5)
#     j = r1 - r
#     js = vpermi2pd(convert(UInt, j), TABLE_EXP_64_0, TABLE_EXP_64_1)

#     expr = expm1b_kernel_5(Val(2), r) # 2^r - 1

#     small_part = vfmadd(js, expr, js)
#     # @show r m j
#     vscalef(small_part, m)
# end
# @inline function vexp_v3(x::AbstractSIMD{8,Float64}, ::True)#, ::Val{N}) where {N}
#     xl2e = mul_ieee(1.4426950408889634, x)
#     r1 = vsreduce(xl2e, Val(0))
#     m = xl2e - r1
#     r = vfmsub(vsreduce(r1 * 16.0, Val(1)), 0.0625, 0.5)
#     j = r1 - r
#     js = vpermi2pd(convert(UInt, j), TABLE_EXP_64_0, TABLE_EXP_64_1)
#     rs = vfnmadd(0.6931471805599453094172321214581765680755001343602552541206800094933936219696955, m+j, x)
#     expr = expm1b_kernel_5(Val(ℯ), r) # 2^r - 1

#     small_part = vfmadd(js, expr, js)
#     @show r m j
#     vscalef(small_part, m)
# end

# _log2(::Val{ℯ}) = 1.4426950408889634
# invlog2hi(::Val{ℯ}) = 0.6931471805599453094172321214581765680755001343602552541206800094933936219696955
# invlog2lo(::Val{ℯ}) = -2.319046813846299615494855463875478650412068000949339362196969553467383712860567e-17
# _log2(::Val{10}) = 3.321928094887362347870319429489390175864831393024580612054756395815934776608624
# invlog2hi(::Val{10}) = 0.3010299956639811952137388947244930267681898814621085413104274611271081892744238
# invlog2lo(::Val{10}) = 2.803728127785170339013117338996875833689572538872891810725576172209659522828247e-18

# Requires two more floating point μops, but 8 less loading μops than the default version.
# This thus microbenchmarks a little worse, but the theory is that using less cache than the
# 256 Float64 * 8 bytes/Float64 = 2 KiB table may improve real world performance / reduce
# random latency.
@inline function vexp_avx512(x::AbstractSIMD{8,Float64}, ::Val{B}) where {B}
    N_float = round(x*LogBo16INV(Val(B), Float64))
    r = muladd(N_float, LogBo16U(Val(B), Float64), x)
    r = muladd(N_float, LogBo16L(Val(B), Float64), r)
    inds = ((trunc(Int64, N_float)%UInt64)) & 0x000000000000000f
    expr = expm1b_kernel_16(Val(B), r)
    js = vpermi2pd(inds, TABLE_EXP_64_0, TABLE_EXP_64_1)
    small_part = vfmadd(js, expr, js)
    res = vscalef(small_part, 0.0625*N_float)
    return res
end

@inline function vexp_avx512(x::AbstractSIMD{W,Float64}, ::Val{B}) where {W,B}
    N_float = muladd(x, LogBo256INV(Val{B}(), Float64), MAGIC_ROUND_CONST(Float64))
    N = target_trunc(reinterpret(UInt64, N_float))
    N_float = N_float - MAGIC_ROUND_CONST(Float64)
    r = fma(N_float, LogBo256U(Val{B}(), Float64), x)
    r = fma(N_float, LogBo256L(Val{B}(), Float64), r)
    # @show (N & 0x000000ff) % Int
    # @show N N & 0x000000ff
    js = vload(VectorizationBase.zero_offsets(stridedpointer(J_TABLE)), (N & 0x000000ff,))
    # k = N >>> 0x00000008
    # small_part = reinterpret(UInt64, vfmadd(js, expm1b_kernel(Val{B}(), r), js))
    small_part = vfmadd(js, expm1b_kernel(Val{B}(), r), js)
    # return reinterpret(Float64, small_part), r, k, N_float, js
    res = vscalef(small_part, 0.00390625*N_float)
    # twopk = (k % UInt64) << 0x0000000000000034
    # res = reinterpret(Float64, twopk + small_part)
    return res
end
@inline function vexp_avx512(x::Union{Float32,AbstractSIMD{<:Any,Float32}}, ::Val{B}) where {B}
    N_float = vfmadd(x, LogBINV(Val{B}(), Float32), MAGIC_ROUND_CONST(Float32))
    N = reinterpret(UInt32, N_float)
    N_float = (N_float - MAGIC_ROUND_CONST(Float32))

    r = fast_fma(N_float, LogBU(Val{B}(), Float32), x, fma_fast())
    r = fast_fma(N_float, LogBL(Val{B}(), Float32), r, fma_fast())

    small_part = expb_kernel(Val{B}(), r)
    res = vscalef(small_part, N_float)
    # twopk = N << 0x00000017
    # res = reinterpret(Float32, twopk + small_part)
    return res
end

else# if !((Sys.ARCH === :x86_64) | (Sys.ARCH === :i686))
  const target_trunc = identity
end

# @inline function vexp_avx512(vu::VecUnroll{1,8,Float64,Vec{8,Float64}}, ::Val{B}) where {B}
#   x, y = data(vu)
#   N_float₁ = round(x*LogBo16INV(Val(B), Float64))
#   N_float₂ = muladd(y, LogBo256INV(Val{B}(), Float64), MAGIC_ROUND_CONST(Float64))
#   r₁ = muladd(N_float₁, LogBo16U(Val(B), Float64), x)
#   N₂ = target_trunc(reinterpret(UInt64, N_float₂))
#   r₁ = muladd(N_float₁, LogBo16L(Val(B), Float64), r₁)
#   js₂ = vload(VectorizationBase.zero_offsets(stridedpointer(J_TABLE)), (N₂ & 0x000000ff,))
#   N_float₂ = N_float₂ - MAGIC_ROUND_CONST(Float64)
#   inds₁ = ((trunc(Int64, N_float₁)%UInt64)) & 0x000000000000000f
#   r₂ = fma(N_float₂, LogBo256U(Val{B}(), Float64), y)
#   expr₁ = expm1b_kernel_16(Val(B), r₁)
#   r₂ = fma(N_float₂, LogBo256L(Val{B}(), Float64), r₂)
#   js₁ = vpermi2pd(inds₁, TABLE_EXP_64_0, TABLE_EXP_64_1)
#   small_part₁ = vfmadd(js₁, expr₁, js₁)
#   small_part₂ = vfmadd(js₂, expm1b_kernel(Val{B}(), r₂), js₂)
#   res₁ = vscalef(small_part₁, 0.0625*N_float₁)
#   res₂ = vscalef(small_part₂, 0.00390625*N_float₂)
#   return VecUnroll((res₁, res₂))
# end

# @inline _vexp(x, ::True) = vexp2( 1.4426950408889634 * x, True() )
# @inline _vexp(x, ::True) = vexp2( mul_ieee(1.4426950408889634, x), True() )
# @inline _vexp10(x, ::True) = vexp2( 3.321928094887362 * x, True() )
# @inline _vexp(x) = _vexp(x, has_feature(Val(:x86_64_avx512f)))
# @inline _vexp10(x) = _vexp10(x, has_feature(Val(:x86_64_avx512f)))

@inline vexp(x::AbstractSIMD, ::True) = vexp_avx512(x, Val(ℯ))
@inline vexp2(x::AbstractSIMD, ::True) = vexp_avx512(x, Val(2))
@inline vexp10(x::AbstractSIMD, ::True) = vexp_avx512(x, Val(10))
# @inline vexp(x::AbstractSIMD{W,Float32}, ::True) where {W} = vexp_generic(x, Val(ℯ))
# @inline vexp2(x::AbstractSIMD{W,Float32}, ::True) where {W} = vexp_generic(x, Val(2))
# @inline vexp10(x::AbstractSIMD{W,Float32}, ::True) where {W} = vexp_generic(x, Val(10))

@inline Base.exp(v::AbstractSIMD{W}) where {W} = vexp(float(v))
@inline Base.exp2(v::AbstractSIMD{W}) where {W} = vexp2(float(v))
@inline Base.exp10(v::AbstractSIMD{W}) where {W} = vexp10(float(v))
@static if (Base.libllvm_version ≥ v"11") & ((Sys.ARCH === :x86_64) | (Sys.ARCH === :i686))
    @inline vexp(v::AbstractSIMD) = vexp(v, has_feature(Val(:x86_64_avx512f)))
    @inline vexp2(v::AbstractSIMD) = vexp2(v, has_feature(Val(:x86_64_avx512f)))
    @inline vexp10(v::AbstractSIMD) = vexp10(v, has_feature(Val(:x86_64_avx512f)))
else
    @inline vexp(v::AbstractSIMD) = vexp(v, False())
    @inline vexp2(v::AbstractSIMD) = vexp2(v, False())
    @inline vexp10(v::AbstractSIMD) = vexp10(v, False())
end
@inline vexp(v::Union{Float32,Float64}) = vexp(v, False())
@inline vexp2(v::Union{Float32,Float64}) = vexp2(v, False())
@inline vexp10(v::Union{Float32,Float64}) = vexp10(v, False())
@inline vexp(v::AbstractSIMD{2,Float32}) = vexp(v, False())
@inline vexp2(v::AbstractSIMD{2,Float32}) = vexp2(v, False())
@inline vexp10(v::AbstractSIMD{2,Float32}) = vexp10(v, False())

####################################################################################################
################################## Non-AVX512 implementation #######################################
####################################################################################################

@inline function vexp_generic_core(x::Union{Float64,AbstractSIMD{<:Any,Float64}}, ::Val{B}) where {B}
    N_float = muladd(x, LogBo256INV(Val{B}(), Float64), MAGIC_ROUND_CONST(Float64))
    N = target_trunc(reinterpret(UInt64, N_float))
    N_float = N_float - MAGIC_ROUND_CONST(Float64)
    r = fast_fma(N_float, LogBo256U(Val{B}(), Float64), x, fma_fast())
    r = fast_fma(N_float, LogBo256L(Val{B}(), Float64), r, fma_fast())
    # @show (N & 0x000000ff) % Int
    js = vload(VectorizationBase.zero_offsets(stridedpointer(J_TABLE)), (N & 0x000000ff,))
    k = N >>> 0x00000008
    small_part = reinterpret(UInt64, vfmadd(js, expm1b_kernel(Val{B}(), r), js))
    # return reinterpret(Float64, small_part), r, k, N_float, js
    twopk = (k % UInt64) << 0x0000000000000034
    res = reinterpret(Float64, twopk + small_part)
    return res
end
@inline function vexp_generic(x::Union{Float64,AbstractSIMD{<:Any,Float64}}, ::Val{B}) where {B}
    res = vexp_generic_core(x, Val{B}())
    res = ifelse(x >= MAX_EXP(Val{B}(), Float64), Inf, res)
    res = ifelse(x <= MIN_EXP(Val{B}(), Float64), 0.0, res)
    res = ifelse(isnan(x), x, res)
    return res
end

@inline function vexp_generic_core(x::Union{Float32,AbstractSIMD{<:Any,Float32}}, ::Val{B}) where {B}
    N_float = vfmadd(x, LogBINV(Val{B}(), Float32), MAGIC_ROUND_CONST(Float32))
    N = reinterpret(UInt32, N_float)
    N_float = (N_float - MAGIC_ROUND_CONST(Float32))

    r = fast_fma(N_float, LogBU(Val{B}(), Float32), x, fma_fast())
    r = fast_fma(N_float, LogBL(Val{B}(), Float32), r, fma_fast())

    small_part = reinterpret(UInt32, expb_kernel(Val{B}(), r))
    twopk = N << 0x00000017
    res = reinterpret(Float32, twopk + small_part)
    return res
end
@inline function vexp_generic(x::Union{Float32,AbstractSIMD{<:Any,Float32}}, ::Val{B}) where {B}
    res = vexp_generic_core(x, Val{B}())
    res = ifelse(x >= MAX_EXP(Val{B}(), Float32), Inf32, res)
    res = ifelse(x <= MIN_EXP(Val{B}(), Float32), 0.0f0, res)
    res = ifelse(isnan(x), x, res)
    return res
end
for (func, base) in (:vexp2=>Val(2), :vexp=>Val(ℯ), :vexp10=>Val(10))
    @eval @inline $func(x, ::False) = vexp_generic(x, $base)
end

####################################################################################################
#################################### LOG HELPERS ###################################################
####################################################################################################

# TODO: move these back to log.jl when when the log implementations there are good & well tested enough to use.

@generated function vgetexp(v::Vec{W,T}) where {W,T<:Union{Float32,Float64}}
    bits = 8W*sizeof(T)
    bits ∈ (128,256,512) || throw(ArgumentError("Vectors are $bits bits, but only 128, 256, and 512 bits are supported."))
    ltyp = LLVM_TYPES[T]
    vtyp = "<$W x $ltyp>"
    dors = T === Float64 ? 'd' : 's'
    instr = "$vtyp @llvm.x86.avx512.mask.getexp.p$(dors).$bits"
    mtyp = W == 16 ? "i16" : "i8"
    if bits == 512
        decl = "declare $instr($vtyp, $vtyp, $mtyp, i32)"
        instrs = "%res = call $instr($vtyp %0, $vtyp undef, $mtyp -1, i32 12)\nret $vtyp %res"
    else
        decl = "declare $instr($vtyp, $vtyp, $mtyp)"
        instrs = "%res = call $instr($vtyp %0, $vtyp undef, $mtyp -1)\nret $vtyp %res"
    end
    arg_syms = [:(data(v))]
    llvmcall_expr(decl, instrs, :(_Vec{$W,$T}), :(Tuple{_Vec{$W,$T}}), vtyp, [vtyp], arg_syms)
end
@generated function vgetmant(v::Vec{W,T}) where {W,T<:Union{Float32,Float64}}
    bits = 8W*sizeof(T)
    bits ∈ (128,256,512) || throw(ArgumentError("Vectors are $bits bits, but only 128, 256, and 512 bits are supported."))
    ltyp = LLVM_TYPES[T]
    vtyp = "<$W x $ltyp>"
    dors = T === Float64 ? 'd' : 's'
    instr = "$vtyp @llvm.x86.avx512.mask.getmant.p$(dors).$bits"
    mtyp = W == 16 ? "i16" : "i8"
    if bits == 512
        decl = "declare $instr($vtyp, i32, $vtyp, $mtyp, i32)"
        instrs = "%res = call $instr($vtyp %0, i32 11, $vtyp undef, $mtyp -1, i32 12)\nret $vtyp %res"
    else
        decl = "declare $instr($vtyp, i32, $vtyp, $mtyp)"
        instrs = "%res = call $instr($vtyp %0, i32 11, $vtyp undef, $mtyp -1)\nret $vtyp %res"
    end
    arg_syms = [:(data(v))]
    llvmcall_expr(decl, instrs, :(_Vec{$W,$T}), :(Tuple{_Vec{$W,$T}}), vtyp, [vtyp], arg_syms)
end
@generated function vgetmant(v::Vec{W,T}, ::Union{StaticInt{N},Val{N}}) where {W,T<:Union{Float32,Float64},N}
    bits = 8W*sizeof(T)
    bits ∈ (128,256,512) || throw(ArgumentError("Vectors are $bits bits, but only 128, 256, and 512 bits are supported."))
    ltyp = LLVM_TYPES[T]
    vtyp = "<$W x $ltyp>"
    dors = T === Float64 ? 'd' : 's'
    instr = "$vtyp @llvm.x86.avx512.mask.getmant.p$(dors).$bits"
    mtyp = W == 16 ? "i16" : "i8"
    if bits == 512
        decl = "declare $instr($vtyp, i32, $vtyp, $mtyp, i32)"
        instrs = "%res = call $instr($vtyp %0, i32 $N, $vtyp undef, $mtyp -1, i32 12)\nret $vtyp %res"
    else
        decl = "declare $instr($vtyp, i32, $vtyp, $mtyp)"
        instrs = "%res = call $instr($vtyp %0, i32 $N, $vtyp undef, $mtyp -1)\nret $vtyp %res"
    end
    arg_syms = [:(data(v))]
    llvmcall_expr(decl, instrs, :(_Vec{$W,$T}), :(Tuple{_Vec{$W,$T}}), vtyp, [vtyp], arg_syms)
end
@generated function vgetmant12(v::Vec{W,T}) where {W,T<:Union{Float32,Float64}}
    bits = 8W*sizeof(T)
    bits ∈ (128,256,512) || throw(ArgumentError("Vectors are $bits bits, but only 128, 256, and 512 bits are supported."))
    ltyp = LLVM_TYPES[T]
    vtyp = "<$W x $ltyp>"
    dors = T === Float64 ? 'd' : 's'
    instr = "$vtyp @llvm.x86.avx512.mask.getmant.p$(dors).$bits"
    mtyp = W == 16 ? "i16" : "i8"
    if bits == 512
        decl = "declare $instr($vtyp, i32, $vtyp, $mtyp, i32)"
        instrs = "%res = call $instr($vtyp %0, i32 12, $vtyp undef, $mtyp -1, i32 8)\nret $vtyp %res"
    else
        decl = "declare $instr($vtyp, i32, $vtyp, $mtyp)"
        instrs = "%res = call $instr($vtyp %0, i32 12, $vtyp undef, $mtyp -1)\nret $vtyp %res"
    end
    arg_syms = [:(data(v))]
    llvmcall_expr(decl, instrs, :(_Vec{$W,$T}), :(Tuple{_Vec{$W,$T}}), vtyp, [vtyp], arg_syms)
end

@generated function vroundscale(v::Vec{W,T}, ::Union{Val{N},StaticInt{N}}) where {W,T,N}
    bits = 8W*sizeof(T)
    bits ∈ (128,256,512) || throw(ArgumentError("Vectors are $bits bits, but only 128, 256, and 512 bits are supported."))
    ltyp = LLVM_TYPES[T]
    vtyp = "<$W x $ltyp>"
    dors = T === Float64 ? 'd' : 's'
    instr = "$vtyp @llvm.x86.avx512.mask.rndscale.p$(dors).$bits"
    mtyp = W == 16 ? "i16" : "i8"
    if bits == 512
        decl = "declare $instr($vtyp, i32, $vtyp, $mtyp, i32)"
        instrs = "%res = call $instr($vtyp %0, i32 $N, $vtyp undef, $mtyp -1, i32 4)\nret $vtyp %res"
    else
        decl = "declare $instr($vtyp, i32, $vtyp, $mtyp)"
        instrs = "%res = call $instr($vtyp %0, i32 $N, $vtyp undef, $mtyp -1)\nret $vtyp %res"
    end
    arg_syms = [:(data(v))]
    llvmcall_expr(decl, instrs, :(_Vec{$W,$T}), :(Tuple{_Vec{$W,$T}}), vtyp, [vtyp], arg_syms)
end

@inline vgetexp(v::VecUnroll) = VecUnroll(fmap(vgetexp, getfield(v, :data)))
@inline vgetmant(v::VecUnroll) = VecUnroll(fmap(vgetmant, getfield(v, :data)))
@inline vgetmant(v::VecUnroll,::Union{StaticInt{N},Val{N}}) where {N} = VecUnroll(fmap(vgetmant, getfield(v, :data), StaticInt{N}()))
@inline Base.significand(v::VecUnroll, ::True) = VecUnroll(fmap(vgetmant12, getfield(v, :data)))
@inline vroundscale(v::VecUnroll,::Union{StaticInt{N},Val{N}}) where {N} = VecUnroll(fmap(vroundscale, getfield(v, :data), StaticInt{N}()))
@inline Base.significand(v::Vec, ::True) = vgetmant12(v)
@inline Base.significand(v::AbstractSIMDVector, ::True) = vgetmant12(Vec(v))

mask_exponent(::Val{Float64}) = 0x000f_ffff_ffff_ffff
set_exponent(::Val{Float64}) = 0x3ff0_0000_0000_0000

mask_exponent(::Val{Float32}) = 0x007fffff
set_exponent(::Val{Float32}) = 0x3f800000

mask_exponent(::Val{Float16}) = 0x07ff
set_exponent(::Val{Float16}) = 0x3c00

@inline function Base.significand(v::AbstractSIMD{W,T}, ::False) where {W,T}
  reinterpret(T, (reinterpret(Base.uinttype(T), v) & mask_exponent(Val(T))) | set_exponent(Val(T)))
end
@inline Base.exponent(v::Vec, ::True)= vgetexp(v)
@inline Base.exponent(v::AbstractSIMDVector, ::True)= vgetexp(Vec(v))
@inline Base.exponent(v::VecUnroll, ::True)= VecUnroll(fmap(vgetexp, getfield(v, :data)))
@static if VERSION ≥ v"1.7.0-beta"
  using Base: inttype
else
  @inline inttype(::Type{T}) where {T} = signed(Base.uinttype(T))
end
@inline function Base.exponent(v::AbstractSIMD{W,T}, ::False) where {W,T}
  I = inttype(T)
  vshift = reinterpret(I, v) >> (Base.significand_bits(T)%I)
  e = ((vshift % Int) & Base.exponent_raw_max(T)) - I(Base.exponent_bias(T))
  convert(T, e % Int32)
end

@inline Base.significand(v::AbstractSIMD{W,T}) where {W,T<:Union{Float32,Float64}} = significand(v, has_feature(Val(:x86_64_avx512f)))
@inline Base.exponent(v::AbstractSIMD{W,T}) where {W,T<:Union{Float32,Float64}} = exponent(v, has_feature(Val(:x86_64_avx512f)))
@inline Base.significand(v::AbstractSIMD{W,Float16}) where {W} = significand(v, False())
@inline Base.exponent(v::AbstractSIMD{W,Float16}) where {W} = exponent(v, False())
@inline Base.ldexp(v::AbstractSIMD, e::AbstractSIMD) = vscalef(v, e, has_feature(Val(:x86_64_avx512f)))

# @inline function vexp2_v2(x::AbstractSIMD{8,Float64})
#     x16 = 16.0*x
#     # x8 = 8x
#     r =  vsreduce(x16, Val(0)) * 0.0625
#     @fastmath  begin
#     m = x - r
#     # m + r = x16, r ∈ (-0.5,0.5]
#     # m/16 + r/16 = x, r ∈ (-1/32, 1/32]
#     # we must now vreduce `mfrac`
#     # return m
#     end
#     # expr = expm1b_kernel_4(Val(2), r)
#     expr = expm1b_kernel_5(Val(2), r)
#     inds = convert(UInt64, vsreduce(m, Val(1)) * 16.0)
#     # inds = (reinterpret(UInt64, mfrac) >> 0x000000000000002d) & 0x000000000000000f
#     # inds = (reinterpret(UInt64, mfrac) >> 0x0000000000000035) & 0x000000000000000f
#     # @show r m mfrac reinterpret(UInt64, m) reinterpret(UInt64, mfrac)
#     js = vpermi2pd(inds, TABLE_EXP_64_0, TABLE_EXP_64_1)
#     # return r, mfrac, js
#     # @show js m 16mfrac r mfrac inds%Int ((inds>>1)%Int)
#     # js = 1.0
#     small_part = vfmadd(js, expr, js)
#     vscalef(small_part, m)#, r, mfrac, js, inds
# end

# for (func, base) in (:vexp2=>Val(2), :vexp=>Val(ℯ), :vexp10=>Val(10))
#     @eval begin
#         @inline function $func(x::AbstractSIMD{W,Float64}) where {W}
#             N_float = muladd(x, LogBo256INV($base, Float64), MAGIC_ROUND_CONST(Float64))
#             N = target_trunc(reinterpret(UInt64, N_float))
#             N_float = N_float - MAGIC_ROUND_CONST(Float64)
#             r = fast_fma(N_float, LogBo256U($base, Float64), x, fma_fast())
#             r = fast_fma(N_float, LogBo256L($base, Float64), r, fma_fast())
#             js = vload(VectorizationBase.zero_offsets(stridedpointer(J_TABLE)), (N & 0x000000ff,))
#             k = N >>> 0x00000008
#             small_part = vfmadd(js, expm1b_kernel($base, r), js)
#             twopk = (k % UInt64) << 0x0000000000000034
#             @show N_float k twopk small_part
#             return small_part
#             # res = reinterpret(Float64, twopk + small_part)
#             # return res
#             vscalef(small_part, N_float)
#         end
#     end
# end


text:
"""
Abstract equation.
"""
abstract type Equation end

"""
    AdvectionEquation(domain, filter = IdentityFilter())

Filtered advection equation.
"""
Base.@kwdef struct AdvectionEquation{D<:Domain,F<:Filter} <: Equation
    domain::D
    filter::F = IdentityFilter()
end

"""
    DiffusionEquation(
        domain,
        filter = IdentityFilter(),
        f = nothing,
        g_a = nothing,
        g_b = nothing
    )

Filtered diffusion equation.
"""
Base.@kwdef struct DiffusionEquation{D<:Domain,F<:Filter} <: Equation
    domain::D
    filter::F = IdentityFilter()
    f::Union{Nothing,Function} = nothing
    g_a::Union{Nothing,Function} = nothing
    g_b::Union{Nothing,Function} = nothing
end

"""
    BurgersEquation(domain, filter = IdentityFilter())

Filtered Burgers equation.
"""
Base.@kwdef struct BurgersEquation{D<:Domain,F<:Filter} <: Equation
    domain::D
    filter::F = IdentityFilter()
end


text:
module TwoDGridWorldUtils

using SparseArrays
using Distributions

import ..TMazeCumulantSchedules
import ..ContGridWorldParams
import ..ContGridWorld
import ..Learner
import ..check_goal
import ..range_check
import ..get_action_probs
import ..GVFHordes
import ..update
import ..Curiosity
import ..GVFSRHordes
import ..SRCreationUtils
import ..FeatureCreator
using ...StatsBase

const TMCS = TMazeCumulantSchedules
const CGWP = ContGridWorldParams
const SRCU = Curiosity.SRCreationUtils

struct GoalTermination <: GVFHordes.GVFParamFuncs.AbstractDiscount
    γ::Float64
end

function Base.get(gt::GoalTermination; state_tp1, kwargs...)
    any(state_tp1[3:end] .!= 0) ? 0.0 : gt.γ
end

struct GoalPolicy <: GVFHordes.GVFParamFuncs.AbstractPolicy
    goal::Int
    normalized::Bool
end

(π::GoalPolicy)(s) = sample(Weights([get(π, state_t=s, action_t=a) for a ∈ 1:4]))

function Base.get(π::GoalPolicy; state_t, action_t, kwargs...)
    cur_y = state_t[1]
    cur_x = state_t[2]

    if π.goal == 1
        boundry_y = π.normalized ? 0.9 : 9.0
        boundry_x = π.normalized ? 0.1 : 1.0
        if cur_x > boundry_x && cur_y < boundry_y
            if CGWP.UP == action_t || CGWP.LEFT == action_t
                0.5
            else
                0.0
            end
        elseif cur_x > boundry_x
            if CGWP.LEFT == action_t
                1.0
            else
                0.0
            end
        elseif cur_y < boundry_y
            if CGWP.UP == action_t
                1.0
            else
                0.0
            end
        else
            if CGWP.UP == action_t
                1.0
            else
                0.0
            end
            # @show action_t, state_t
            # throw("What?")
        end
    elseif π.goal == 2
        boundry_y = π.normalized ? 0.1 : 1.0
        boundry_x = π.normalized ? 0.1 : 1.0
        if cur_x > boundry_x && cur_y > boundry_y
            if CGWP.DOWN == action_t || CGWP.LEFT == action_t
                0.5
            else
                0.0
            end
        elseif cur_x > boundry_x
            if CGWP.LEFT == action_t
                1.0
            else
                0.0
            end
        elseif cur_y > boundry_y
            if CGWP.DOWN == action_t
                1.0
            else
                0.0
            end
        else
            if CGWP.DOWN == action_t
                1.0
            else
                0.0
            end
            # @show action_t, state_t
            # throw("What?")
        end
    elseif π.goal == 3
        boundry_y = π.normalized ? 0.9 : 9.0
        boundry_x = π.normalized ? 0.9 : 9.0
        if cur_x < boundry_x && cur_y < boundry_y
            if CGWP.UP == action_t || CGWP.RIGHT == action_t
                0.5
            else
                0.0
            end
        elseif cur_x < boundry_x
            if CGWP.RIGHT == action_t
                1.0
            else
                0.0
            end
        elseif cur_y < boundry_y
            if CGWP.UP == action_t
                1.0
            else
                0.0
            end
        else
            if CGWP.UP == action_t
                1.0
            else
                0.0
            end
            # @show action_t, state_t
            # throw("What?")
        end
    elseif π.goal == 4
        boundry_y = π.normalized ? 0.1 : 1.0
        boundry_x = π.normalized ? 0.9 : 9.0
        if cur_x < boundry_x && cur_y > boundry_y
            if CGWP.DOWN == action_t || CGWP.RIGHT == action_t
                0.5
            else
                0.0
            end
        elseif cur_x < boundry_x
            if CGWP.RIGHT == action_t
                1.0
            else
                0.0
            end
        elseif cur_y > boundry_y
            if CGWP.DOWN == action_t
                1.0
            else
                0.0
            end
        else
            if CGWP.DOWN == action_t
                1.0
            else
                0.0
            end
            # @show action_t, state_t
            # throw("What?")
        end
    end
end

struct NaiveGoalPolicy <: GVFHordes.GVFParamFuncs.AbstractPolicy
    goal::Int
end

(π::NaiveGoalPolicy)(s) = sample(Weights([Base.get(π, state_t=s, action_t=a) for a ∈ 1:4]))

function Base.get(π::NaiveGoalPolicy; state_t, action_t, kwargs...)
    if π.goal == 1
        if CGWP.UP == action_t || CGWP.LEFT == action_t
            0.5
        else
            0.0
        end
    elseif π.goal == 2
        if CGWP.DOWN == action_t || CGWP.LEFT == action_t
            0.5
        else
            0.0
        end
    elseif π.goal == 3
        if CGWP.UP == action_t || CGWP.RIGHT == action_t
            0.5
        else
            0.0
        end
    elseif π.goal == 4
        if CGWP.DOWN == action_t || CGWP.RIGHT == action_t
            0.5
        else
            0.0
        end
    end
end

function create_demons(parsed, demon_projected_fc = nothing)
    action_space = 4
    demons = if parsed["demon_learner"] != "SR"
        GVFHordes.Horde(
            [GVFHordes.GVF(GVFHordes.GVFParamFuncs.FeatureCumulant(i+2),
                           GoalTermination(parsed["demon_gamma"]),
                           NaiveGoalPolicy(i)) for i in 1:4])
    elseif parsed["demon_learner"] == "SR"
        @assert demon_projected_fc != nothing
        pred_horde =  GVFHordes.Horde(
                [GVFHordes.GVF(GVFHordes.GVFParamFuncs.FeatureCumulant(i+2),
                     GVFHordes.GVFParamFuncs.ConstantDiscount(0.0),
                     NaiveGoalPolicy(i)) for i in 1:4])

        SF_policies = [NaiveGoalPolicy(i) for i in 1:4]
        SF_discounts = [GoalTermination(parsed["demon_gamma"]) for i in 1:4]
        num_SFs = length(SF_policies)
        SF_horde = SRCU.create_SF_horde(SF_policies, SF_discounts, demon_projected_fc,1:action_space)
        GVFSRHordes.SRHorde(pred_horde, SF_horde, num_SFs, demon_projected_fc)
    else
        throw(ArgumentError("Cannot create demons"))
    end
    return demons
end

function b_π(state_constructor, learner, exploration_strategy; kwargs...)
    s = state_constructor(kwargs[:state_t])
    preds = learner(s)
    return exploration_strategy(preds)[kwargs[:action_t]]
end

function make_behaviour_gvf(learner, γ, state_constructor, expl_strat) #discount, state_constructor_func, learner, exploration_strategy)
    GVF_policy = GVFHordes.GVFParamFuncs.FunctionalPolicy((;kwargs...) -> b_π(state_constructor, learner, expl_strat; kwargs...))
    BehaviourGVF = GVFHordes.GVF(GVFHordes.GVFParamFuncs.RewardCumulant(),
                       GoalTermination(γ),
                       GVF_policy)
end

function check_goal(goal, state, epsilon=0.0)
    cur_y = state[1]
    cur_x = state[2]
    if goal == 1
        boundry_y = 0.9 - epsilon
        boundry_x = 0.1 + epsilon
        cur_y > boundry_y && cur_x < boundry_x
    elseif goal == 2
        boundry_y = 0.1 + epsilon
        boundry_x = 0.1 + epsilon
        cur_y < boundry_y && cur_x < boundry_x
    elseif goal == 3
        boundry_y = 0.9 - epsilon
        boundry_x = 0.9 - epsilon
        cur_y > boundry_y && cur_x > boundry_x
    elseif goal == 4
        boundry_y = 0.1 + epsilon
        boundry_x = 0.9 - epsilon
        cur_y < boundry_y && cur_x > boundry_x
    else
        false
    end
end

####
# Ideal Feature Creator
####
struct IdealDemonFeatures <: FeatureCreator
end

function project_features(fc::IdealDemonFeatures, state, action, next_state)
    new_state = sparsevec(convert(Array{Int,1}, [check_goal(i, next_state) for i in 1:4]))
    return new_state
end

(FP::IdealDemonFeatures)(s_t,a_t,s_tp1) = project_features(FP, s_t, a_t, s_tp1)
Base.size(FP::IdealDemonFeatures) = 4

####
# Ideal State Action Feature Creator
####
struct IdealStateActionDemonFeatures <: FeatureCreator
    num_actions::Int
end

function project_features(fc::IdealStateActionDemonFeatures, s_t, a_t, s_tp1)
    goal_ind = findfirst([check_goal(i, s_tp1) for i in 1:4])
    reward_feature = zeros(Int(fc.num_actions*4))
    if !(goal_ind isa Nothing)
        reward_feature[Int((goal_ind-1)*fc.num_actions + a_t)] = 1
    end
    return sparsevec(reward_feature)
end

(FP::IdealStateActionDemonFeatures)(s_t,a_t,s_tp1) = project_features(FP, s_t, a_t, s_tp1)
Base.size(FP::IdealStateActionDemonFeatures) = 4 * FP.num_actions

struct MarthaIdealDemonFeatures <: FeatureCreator
end

function project_features(fc::MarthaIdealDemonFeatures, state)
    new_state = sparsevec(convert(Array{Int,1}, [check_goal(i, state, 0.05) for i in 1:4]))
    return new_state
end

(FP::MarthaIdealDemonFeatures)(state) = project_features(FP, state)
Base.size(FP::MarthaIdealDemonFeatures) = 4

struct StateAggregation <: FeatureCreator
    bins_per_dim::Int
end

function project_features(sa::StateAggregation, state)
    bpd = sa.bins_per_dim
    new_state = spzeros(Int, bpd*bpd)
    y = state[1]
    x = state[2]

    idx_y = bpd
    idx_x = bpd
    for i in 1:(bpd)
        if y < i//bpd
            idx_y = i
        end
        if x < i//bpd
            idx_x = i
        end
    end
    idx = (idx_x - 1) * bpd + idx_y
    new_state[idx] = 1
    new_state
end

(FP::StateAggregation)(state) = project_features(FP, state)
Base.size(FP::StateAggregation) = begin; bpd = FP.bins_per_dim; bpd*bpd; end

struct SmallStateAggregation <: FeatureCreator end

function project_features(::SmallStateAggregation, state)
    new_state = spzeros(Int, 9)

    y = state[1]
    x = state[2]
    idx_y = if y < 1//3
        1
    elseif y < 2//3
        2
    else
        3
    end

    idx_x = if x < 1//3
        1
    elseif x < 2//3
        2
    else
        3
    end

    idx = (idx_x - 1) * 3 + idx_y
    new_state[idx] = 1
    new_state
end

(fp::SmallStateAggregation)(state) = project_features(fp, state)
Base.size(::SmallStateAggregation) = 9

DrifterDistractor(parsed) = begin
    c_dist = Uniform(parsed["constant_target"][1],parsed["constant_target"][2])
    c1,c2 = rand(c_dist,2)
    if "drifter" ∈ keys(parsed)
        TMCS.DrifterDistractor(
            c1,
            c2,
            parsed["drifter"][1],
            parsed["drifter"][2],
            parsed["distractor"][1],
            parsed["distractor"][2])
    else
        TMCS.DrifterDistractor(
            c1,
            c2,
            parsed["drifter_init"],
            parsed["drifter_std"],
            parsed["distractor_mean"],
            parsed["distractor_std"])
    end
end

DrifterDrifterDistractor(parsed) = begin
    c_dist = Uniform(parsed["constant_target"][1],parsed["constant_target"][2])
    c1 = rand(c_dist)
    if "drifter_1" ∈ keys(parsed)
        TMCS.DrifterDrifterDistractor(
            c1,
            parsed["drifter_1"][1],
            parsed["drifter_1"][2],
            parsed["drifter_2"][1],
            parsed["drifter_2"][2],
            parsed["distractor"][1],
            parsed["distractor"][2])
    else
        TMCS.DrifterDrifterDistractor(
            c1,
            parsed["drifter_1_init"],
            parsed["drifter_1_std"],
            parsed["drifter_2_init"],
            parsed["drifter_2_std"],
            parsed["distractor_mean"],
            parsed["distractor_std"])
    end
end

function get_cumulant_schedule(parsed)
    sched = parsed["cumulant_schedule"]
    if parsed["cumulant_schedule"] == "DrifterDistractor"
        DrifterDistractor(parsed)
    elseif sched == "DrifterDrifterDistractor"
        DrifterDrifterDistractor(parsed)
    elseif parsed["cumulant_schedule"] == "Constant"
        if parsed["cumulant"] isa Number
            TMCS.Constant(parsed["cumulant"])
        else
            TMCS.Constant(parsed["cumulant"]...)
        end

    else
        throw("$(sched) Not Implemented")
    end
end

function get_true_values(env::Curiosity.ContGridWorld, eval_set)
    copy_eval_est = deepcopy(eval_set)
    num_gvfs = 4
    goal_cumulants = TMCS.get_cumulant_eval_values(env.cumulant_schedule)
    for i in 1:num_gvfs
        copy_eval_est[i, :] .*= goal_cumulants[i]
    end
    return copy_eval_est
end

function get_true_values(env::Curiosity.ContGridWorld, eval_set, gvf_idx)
    copy_eval_est = deepcopy(eval_set)
    goal_cumulants = TMCS.get_cumulant_eval_values(env.cumulant_schedule)
    copy_eval_est .*= goal_cumulants[gvf_idx]
    return copy_eval_est
end

Base.@kwdef mutable struct RoundRobinPolicy{GT} <: Learner
    cur_goal::Int = 1
    num_goals::Int = 4
    goal_type::GT = NaiveGoalPolicy
    update = nothing
end

Curiosity.update!(rrp::RoundRobinPolicy, args...) = if args[end-2]
    rrp.cur_goal += 1
    if rrp.cur_goal > rrp.num_goals
        rrp.cur_goal = 1
    end
end

Base.get(rrp::RoundRobinPolicy; state_t, action_t, kwargs...) =
    get(rrp.goal_type(rrp.cur_goal); state_t = state_t, action_t = action_t)

function get_action_probs(π::RoundRobinPolicy, features, state)
    [get(π; state_t=state, action_t=a) for a ∈ 1:4]
end

end

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
# Invertible conditional HINT network from Kruse et. al (2020)
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: February 2020

export NetworkConditionalHINT

"""
    CH = NetworkConditionalHINT(nx, ny, n_in, batchsize, n_hidden, depth; k1=3, k2=3, p1=1, p2=1, s1=1, s2=1)

 Create a conditional HINT network for data-driven generative modeling based
 on the change of variables formula.

 *Input*:

 - `nx`, `ny`, `n_in`, `batchsize`: spatial dimensions, number of channels and batchsize of input tensors `X` and `Y`

 - `n_hidden`: number of hidden units in residual blocks

 - `depth`: number network layers

 - `k1`, `k2`: kernel size for first and third residual layer (`k1`) and second layer (`k2`)

 - `p1`, `p2`: respective padding sizes for residual block layers

 - `s1`, `s2`: respective strides for residual block layers

 *Output*:

 - `CH`: conditioinal HINT network

 *Usage:*

 - Forward mode: `Zx, Zy, logdet = CH.forward(X, Y)`

 - Inverse mode: `X, Y = CH.inverse(Zx, Zy)`

 - Backward mode: `ΔX, X = CH.backward(ΔZx, ΔZy, Zx, Zy)`

 *Trainable parameters:*

 - None in `CH` itself

 - Trainable parameters in activation normalizations `CH.AN_X[i]` and `CH.AN_Y[i]`,
 and in coupling layers `CH.CL[i]`, where `i` ranges from `1` to `depth`.

 See also: [`ActNorm`](@ref), [`ConditionalLayerHINT!`](@ref), [`get_params`](@ref), [`clear_grad!`](@ref)
"""
mutable struct NetworkConditionalHINT <: InvertibleNetwork
    AN_X::AbstractArray{ActNorm, 1}
    AN_Y::AbstractArray{ActNorm, 1}
    CL::AbstractArray{ConditionalLayerHINT, 1}
    logdet::Bool
    is_reversed::Bool
end

@Flux.functor NetworkConditionalHINT

# Constructor
function NetworkConditionalHINT(nx, ny, n_in, batchsize, n_hidden, depth; k1=3, k2=3, p1=1, p2=1, s1=1, s2=1, logdet=true)

    AN_X = Array{ActNorm}(undef, depth)
    AN_Y = Array{ActNorm}(undef, depth)
    CL = Array{ConditionalLayerHINT}(undef, depth)

    # Create layers
    for j=1:depth
        AN_X[j] = ActNorm(n_in; logdet=logdet)
        AN_Y[j] = ActNorm(n_in; logdet=logdet)
        CL[j] = ConditionalLayerHINT(nx, ny, n_in, n_hidden, batchsize;
                                     permute=true, k1=k1, k2=k2, p1=p1, p2=p2, s1=s1, s2=s2, logdet=logdet)
    end

    return NetworkConditionalHINT(AN_X, AN_Y, CL, logdet, false)
end

# Forward pass and compute logdet
function forward(X, Y, CH::NetworkConditionalHINT; logdet=nothing, x_lane=false)
    isnothing(logdet) ? logdet = (CH.logdet && ~CH.is_reversed) : logdet = logdet

    depth = length(CH.CL)
    logdet_ = 0f0
    logdet_x = 0f0
    for j=1:depth
        logdet ? (X_, logdet1) = CH.AN_X[j].forward(X) : X_ = CH.AN_X[j].forward(X)
        logdet ? (Y_, logdet2) = CH.AN_Y[j].forward(Y) : Y_ = CH.AN_Y[j].forward(Y)
        logdet ? (X, Y, logdet3) = CH.CL[j].forward(X_, Y_; x_lane=x_lane) : (X, Y) = CH.CL[j].forward(X_, Y_)
        logdet && (logdet_ += (logdet1 + logdet2 + logdet3))
        logdet && (logdet_x += (logdet1 + logdet3))
    end
    if x_lane == false
        logdet ? (return X, Y, logdet_) : (return X, Y)
    else
        logdet ? (return X, Y, logdet_x) : (return X, Y)
    end
end

# Inverse pass and compute gradients
function inverse(Zx, Zy, CH::NetworkConditionalHINT; logdet=nothing, x_lane=false)
    isnothing(logdet) ? logdet = (CH.logdet && CH.is_reversed) : logdet = logdet

    depth = length(CH.CL)
    logdet_ = 0f0
    logdet_x = 0f0
    for j=depth:-1:1
        logdet ? (Zx_, Zy_, logdet1) = CH.CL[j].inverse(Zx, Zy; logdet=true, x_lane=x_lane) : (Zx_, Zy_) = CH.CL[j].inverse(Zx, Zy; logdet=false)
        logdet ? (Zy, logdet2) = CH.AN_Y[j].inverse(Zy_; logdet=true) : Zy = CH.AN_Y[j].inverse(Zy_; logdet=false)
        logdet ? (Zx, logdet3) = CH.AN_X[j].inverse(Zx_; logdet=true) : Zx = CH.AN_X[j].inverse(Zx_; logdet=false)
        logdet && (logdet_ += (logdet1 + logdet2 + logdet3))
        logdet && (logdet_x += (logdet1 + logdet3))
    end
    if x_lane == false
        logdet ? (return Zx, Zy, logdet_) : (return Zx, Zy)
    else
        logdet ? (return Zx, Zy, logdet_x) : (return Zx, Zy)
    end
end

# Backward pass and compute gradients
function backward(ΔZx, ΔZy, Zx, Zy, CH::NetworkConditionalHINT)
    depth = length(CH.CL)
    for j=depth:-1:1
        ΔZx_, ΔZy_, Zx_, Zy_ = CH.CL[j].backward(ΔZx, ΔZy, Zx, Zy)
        ΔZx, Zx = CH.AN_X[j].backward(ΔZx_, Zx_)
        ΔZy, Zy = CH.AN_Y[j].backward(ΔZy_, Zy_)
    end
    return ΔZx, ΔZy, Zx, Zy
end

# Backward reverse pass and compute gradients
function backward_inv(ΔX, ΔY, X, Y, CH::NetworkConditionalHINT)
    depth = length(CH.CL)
    for j=1:depth
        ΔX_, X_ = backward_inv(ΔX, X, CH.AN_X[j])
        ΔY_, Y_ = backward_inv(ΔY, Y, CH.AN_Y[j])
        ΔX, ΔY, X, Y = backward_inv(ΔX_, ΔY_, X_, Y_, CH.CL[j])
    end
    return ΔX, ΔY, X, Y
end

# Forward pass and compute logdet
function forward_Y(Y, CH::NetworkConditionalHINT)
    depth = length(CH.CL)
    for j=1:depth
        Y_ = CH.AN_Y[j].forward(Y; logdet=false)
        Y = CH.CL[j].forward_Y(Y_)
    end
    return Y
end

# Inverse pass and compute gradients
function inverse_Y(Zy, CH::NetworkConditionalHINT)
    depth = length(CH.CL)
    for j=depth:-1:1
        Zy_ = CH.CL[j].inverse_Y(Zy)
        Zy = CH.AN_Y[j].inverse(Zy_; logdet=false)
    end
    return Zy
end

# Clear gradients
function clear_grad!(CH::NetworkConditionalHINT)
    depth = length(CH.CL)
    for j=1:depth
        clear_grad!(CH.AN_X[j])
        clear_grad!(CH.AN_Y[j])
        clear_grad!(CH.CL[j])
    end
end

# Get parameters
function get_params(CH::NetworkConditionalHINT)
    depth = length(CH.CL)
    p = []
    for j=1:depth
        p = cat(p, get_params(CH.AN_X[j]); dims=1)
        p = cat(p, get_params(CH.AN_Y[j]); dims=1)
        p = cat(p, get_params(CH.CL[j]); dims=1)
    end
    return p
end

# Put parameters
function put_params!(CH::NetworkConditionalHINT, Params::Array{Any,1})
    depth = length(CH.CL)
    idx_s = 1
    counter = 0

    for j = 1:depth
        idx_s += counter
        counter = 2
        put_params!(CH.AN_X[j], Params[idx_s:idx_s+counter-1])

        idx_s += counter
        counter = 2
        put_params!(CH.AN_Y[j], Params[idx_s:idx_s+counter-1])

        idx_s += counter
        counter = 0
        nlayers_x = length(CH.CL[j].CL_X.CL)
        isnothing(CH.CL[j].CL_X.C) ? counter += 5*nlayers_x : counter += 5*nlayers_x + 3
        nlayers_y = length(CH.CL[j].CL_Y.CL)
        isnothing(CH.CL[j].CL_Y.C) ? counter += 5*nlayers_y : counter += 5*nlayers_y + 3
        counter += 5
        ~isnothing(CH.CL[j].C_X) && (counter += 3)
        ~isnothing(CH.CL[j].C_Y) && (counter += 3)

        put_params!(CH.CL[j], Params[idx_s:idx_s+counter-1])
    end
end

# Set is_reversed flag in full network tree
function tag_as_reversed!(CH::NetworkConditionalHINT, tag::Bool)
    depth = length(CH.CL)
    CH.is_reversed = tag
    for j=1:depth
        tag_as_reversed!(CH.AN_X[j], tag)
        tag_as_reversed!(CH.AN_Y[j], tag)
        tag_as_reversed!(CH.CL[j], tag)
    end
    return CH
end


text:
using VLConstraintBasedModelGenerationUtilities
using DataFrames

# setup path to sequence file -
path_to_gene_sequence_file = "/Users/jeffreyvarner/Desktop/julia_work/VLConstraintBasedModelGenerationUtilities.jl/test/data/G-MZ373340.fasta"

# load -
gene_table = build_gene_table(path_to_gene_sequence_file) |> check

# get the sequence -
gene_seq = gene_table[!,:gene_sequence][1]
gene_seq_name = gene_table[!,:id][1]

# transcription table -
transcription_table = build_transcription_reaction_table(gene_seq_name, gene_seq; polymeraseSymbol=:RdRp) |> check

# setup path to protein sequence file -
path_to_protein_sequence_file = "/Users/jeffreyvarner/Desktop/julia_work/VLConstraintBasedModelGenerationUtilities.jl/test/data/P-MZ373340.fasta"

# load -
protein_table = build_protein_table(path_to_protein_sequence_file) |> check

# build the translation table -
translation_dictionary = Dict{String,DataFrame}()
translation_array = Array{DataFrame,1}()
(number_of_proteins, _) = size(protein_table)
for protein_index = 1:number_of_proteins

    local_seq = protein_table[protein_index,:protein_sequence]
    local_seq_name = protein_table[protein_index,:id]
    translation_table = build_translation_reaction_table(local_seq_name, local_seq) |> check

    # grab -
    # translation_dictionary[local_seq_name] = translation_table
    push!(translation_array, translation_table)
end

# add together -
master_translation_table = reduce(vcat, translation_array)

# lastly - let's build the transport reactions -
transport_reaction_table = build_transport_reaction_table() |> check

# append reactions into master reaction table -
master_reaction_table = reduce(vcat, [transcription_table, master_translation_table, transport_reaction_table])