text:
# library and API version management

"""
Returns the CUDA driver version as a VersionNumber.
"""
function version()
    version_ref = Ref{Cint}()
    @apicall(:cuDriverGetVersion, (Ptr{Cint},), version_ref)

    major = version_ref[] ÷ 1000
    minor = mod(version_ref[], 100) ÷ 10

    return VersionNumber(major, minor)
end

#
# API versioning
#

const mapping = Dict{Symbol,Symbol}()
const minreq = Dict{Symbol,VersionNumber}()

function resolve(f::Symbol)
    global mapping, version_requirements
    versioned_f = get(mapping, f, f)
    if versioned_f == Symbol()
        throw(CuVersionError(f, minreq[f]))
    end
    return versioned_f
end

function __init_versioning__(driver_version=version())  # NOTE: performs unversioned API call
    ## Versioned API calls from cuda.h

    if driver_version >= v"3.2"
        mapping[:cuDeviceTotalMem]           = :cuDeviceTotalMem_v2
        mapping[:cuCtxCreate]                = :cuCtxCreate_v2
        mapping[:cuModuleGetGlobal]          = :cuModuleGetGlobal_v2
        mapping[:cuMemGetInfo]               = :cuMemGetInfo_v2
        mapping[:cuMemAlloc]                 = :cuMemAlloc_v2
        mapping[:cuMemAllocPitch]            = :cuMemAllocPitch_v2
        mapping[:cuMemFree]                  = :cuMemFree_v2
        mapping[:cuMemGetAddressRange]       = :cuMemGetAddressRange_v2
        mapping[:cuMemAllocHost]             = :cuMemAllocHost_v2
        mapping[:cuMemHostGetDevicePointer]  = :cuMemHostGetDevicePointer_v2
        mapping[:cuMemcpyHtoD]               = :cuMemcpyHtoD_v2
        mapping[:cuMemcpyDtoH]               = :cuMemcpyDtoH_v2
        mapping[:cuMemcpyDtoD]               = :cuMemcpyDtoD_v2
        mapping[:cuMemcpyDtoA]               = :cuMemcpyDtoA_v2
        mapping[:cuMemcpyAtoD]               = :cuMemcpyAtoD_v2
        mapping[:cuMemcpyHtoA]               = :cuMemcpyHtoA_v2
        mapping[:cuMemcpyAtoH]               = :cuMemcpyAtoH_v2
        mapping[:cuMemcpyAtoA]               = :cuMemcpyAtoA_v2
        mapping[:cuMemcpyHtoAAsync]          = :cuMemcpyHtoAAsync_v2
        mapping[:cuMemcpyAtoHAsync]          = :cuMemcpyAtoHAsync_v2
        mapping[:cuMemcpy2D]                 = :cuMemcpy2D_v2
        mapping[:cuMemcpy2DUnaligned]        = :cuMemcpy2DUnaligned_v2
        mapping[:cuMemcpy3D]                 = :cuMemcpy3D_v2
        mapping[:cuMemcpyHtoDAsync]          = :cuMemcpyHtoDAsync_v2
        mapping[:cuMemcpyDtoHAsync]          = :cuMemcpyDtoHAsync_v2
        mapping[:cuMemcpyDtoDAsync]          = :cuMemcpyDtoDAsync_v2
        mapping[:cuMemcpy2DAsync]            = :cuMemcpy2DAsync_v2
        mapping[:cuMemcpy3DAsync]            = :cuMemcpy3DAsync_v2
        mapping[:cuMemsetD8]                 = :cuMemsetD8_v2
        mapping[:cuMemsetD16]                = :cuMemsetD16_v2
        mapping[:cuMemsetD32]                = :cuMemsetD32_v2
        mapping[:cuMemsetD2D8]               = :cuMemsetD2D8_v2
        mapping[:cuMemsetD2D16]              = :cuMemsetD2D16_v2
        mapping[:cuMemsetD2D32]              = :cuMemsetD2D32_v2
        mapping[:cuArrayCreate]              = :cuArrayCreate_v2
        mapping[:cuArrayGetDescriptor]       = :cuArrayGetDescriptor_v2
        mapping[:cuArray3DCreate]            = :cuArray3DCreate_v2
        mapping[:cuArray3DGetDescriptor]     = :cuArray3DGetDescriptor_v2
        mapping[:cuTexRefSetAddress]         = :cuTexRefSetAddress_v2
        mapping[:cuTexRefGetAddress]         = :cuTexRefGetAddress_v2
        mapping[:cuGraphicsResourceGetMappedPointer] = :cuGraphicsResourceGetMappedPointer_v2
    end

    if driver_version >= v"4.0"
        mapping[:cuCtxDestroy]               = :cuCtxDestroy_v2
        mapping[:cuCtxPopCurrent]            = :cuCtxPopCurrent_v2
        mapping[:cuCtxPushCurrent]           = :cuCtxPushCurrent_v2
        mapping[:cuStreamDestroy]            = :cuStreamDestroy_v2
        mapping[:cuEventDestroy]             = :cuEventDestroy_v2
    end

    if driver_version >= v"4.1"
        mapping[:cuTexRefSetAddress2D]       = :cuTexRefSetAddress2D_v3
    end

    if driver_version >= v"6.5"
        mapping[:cuLinkCreate]              = :cuLinkCreate_v2
        mapping[:cuLinkAddData]             = :cuLinkAddData_v2
        mapping[:cuLinkAddFile]             = :cuLinkAddFile_v2
    end

    if driver_version >= v"6.5"
        mapping[:cuMemHostRegister]         = :cuMemHostRegister_v2
        mapping[:cuGraphicsResourceSetMapFlags] = :cuGraphicsResourceSetMapFlags_v2
    end

    if v"3.2" <= driver_version < v"4.1"
        mapping[:cuTexRefSetAddress2D]      = :cuTexRefSetAddress2D_v2
    end

    ## Version-dependent features

    minreq[:cuLinkCreate]       = v"5.5"
    minreq[:cuLinkDestroy]      = v"5.5"
    minreq[:cuLinkComplete]     = v"5.5"
    minreq[:cuLinkAddFile]      = v"5.5"
    minreq[:cuLinkAddData]      = v"5.5"

    minreq[:cuDummyAvailable]   = v"0"      # non-existing functions
    minreq[:cuDummyUnavailable] = v"999"    # for testing purposes

    # explicitly mark unavailable symbols, signaling `resolve` to error out
    for (api_function, minimum_version) in minreq
        if driver_version < minimum_version
            mapping[api_function]      = Symbol()
        end
    end
end


text:
replacements = split("""Al => ThF
Al => ThRnFAr
B => BCa
B => TiB
B => TiRnFAr
Ca => CaCa
Ca => PB
Ca => PRnFAr
Ca => SiRnFYFAr
Ca => SiRnMgAr
Ca => SiTh
F => CaF
F => PMg
F => SiAl
H => CRnAlAr
H => CRnFYFYFAr
H => CRnFYMgAr
H => CRnMgYFAr
H => HCa
H => NRnFYFAr
H => NRnMgAr
H => NTh
H => OB
H => ORnFAr
Mg => BF
Mg => TiMg
N => CRnFAr
N => HSi
O => CRnFYFAr
O => CRnMgAr
O => HP
O => NRnFAr
O => OTi
P => CaP
P => PTi
P => SiRnFAr
Si => CaSi
Th => ThCa
Ti => BP
Ti => TiTi
e => HF
e => NAl
e => OMg""", "\n")

inpt = """ORnPBPMgArCaCaCaSiThCaCaSiThCaCaPBSiRnFArRnFArCaCaSiThCaCaSiThCaCaCaCaCaCaSiRnFYFArSiRnMgArCaSiRnPTiTiBFYPBFArSiRnCaSiRnTiRnFArSiAlArPTiBPTiRnCaSiAlArCaPTiTiBPMgYFArPTiRnFArSiRnCaCaFArRnCaFArCaSiRnSiRnMgArFYCaSiRnMgArCaCaSiThPRnFArPBCaSiRnMgArCaCaSiThCaSiRnTiMgArFArSiThSiThCaCaSiRnMgArCaCaSiRnFArTiBPTiRnCaSiAlArCaPTiRnFArPBPBCaCaSiThCaPBSiThPRnFArSiThCaSiThCaSiThCaPTiBSiRnFYFArCaCaPRnFArPBCaCaPBSiRnTiRnFArCaPRnFArSiRnCaCaCaSiThCaRnCaFArYCaSiRnFArBCaCaCaSiThFArPBFArCaSiRnFArRnCaCaCaFArSiRnFArTiRnPMgArF"""

solutions = Set{ASCIIString}()
push!(solutions, inpt) # to cover the case of no matches. Will be subtracted at the end
subst = r"(\w+) => (\w+)"

#part 1
for line in replacements
	m = match(subst, line)
	for i in 1:length(inpt)
		sol = inpt[1:i-1]*replace(inpt[i:end], m.captures[1], m.captures[2], 1)
		push!(solutions, sol)
	end
end
println(length(solutions)-1) # one is the original molecule, in case of no substitutions

# part 2
# solution from here https://www.reddit.com/r/adventofcode/comments/3xflz8/day_19_solutions/cy4etju
# I can't figure this out ...

ext: jl

max_stars_repo_name: jstrube/adventofcode

max_issues_repo_name: jstrube/adventofcode

max_forks_repo_name: jstrube/adventofcode


text:
module Markov
include("parameter.jl")
include("grid.jl")
include("observable.jl")
include("utility/utility.jl")
include("diag/propagator.jl")
import .Propagator: green
using Random, StaticArrays, Printf, Dates
import .Vertex4, .Ver4Test, .Polar
const UpdateNum = 6
const INCREASE_ORDER, DECREASE_ORDER, CHANGE_EXTTAU, CHANGE_EXTK, CHANGE_TAU, CHANGE_K =
    1:UpdateNum
const Name = [
    "increase_order",
    "decrease_order",
    "change_ExtTau",
    "change_ExtK",
    "change_Tau",
    "change_K",
]
const Accepted = (@MArray zeros(Float, UpdateNum, Order + 1)) .+ 1.0e-10
const Proposed = (@MArray zeros(Float, UpdateNum, Order + 1)) .+ 1.0e-10

const curr = Main.Curr
const rng = Main.Curr.rng
const varK = Main.Curr.K
const varT = Main.Curr.T

const ver4 = Vector{Vertex4.Ver4}(undef, 0)
const oneBody = Observable.OneBody()
const polar = Vector{Polar.Polarization}(undef, 0)

# function init(_counter, _rng)
function init()
    #######  initialize MC variables  ################################
    # println(oneBody.norm)
    # println(typeof(oneBody))
    ###### initialized diagram trees #######################################

    curr.absWeight = Markov.eval(curr.order)
end

@fastmath function test() end

function eval(order)
    if order == 0
        return 1.0
    else
        Tau = varT[LastTidx] - varT[1]
        @assert 0.0 < Tau < Beta "t out of range"
        @assert abs(Grid.K.grid[curr.extKidx] - curr.K[1][1]) < 1.0e-15 "ExtK doesn't match!"
        @assert abs(Grid.tau.grid[curr.extTidx] - curr.T[LastTidx]) < 1.0e-15 "ExtT doesn't match! $(Grid.tau.grid[curr.extTidx]) vs $(curr.T[LastTidx])"
        @assert 0 < curr.extKidx <= KGridSize "K out of range"
        @assert 0 < curr.extTidx <= TauGridSize "Tau out of range"
        gWeight = green(Tau, varK[2]) * green(-Tau, varK[2] + varK[1])
        return -SPIN * gWeight * PhaseFactor
    end
    # return
    # order == 0 && return 1.0
    # if DiagType == POLAR
    #     return Polar.eval(polar[order])
    # elseif DiagType == GAMMA
    #     Vertex4.eval(ver4[order], varK[INL], varK[OUTL], varK[INR], varK[OUTR], 5, true)
    #     chanWeight = sum(ver4[order].weight)
    #     return chanWeight[DI] + chanWeight[EX] / SPIN
    # else
    #     throw("Not implemented!")
    # end
end

function measure()

    factor = 1.0 / curr.absWeight / ReWeight[curr.order+1]
    # println("$(curr.order), $(curr.absWeight), $(ReWeight[curr.order + 1])")

    if curr.order == 1
        Tau = curr.T[LastTidx]
        weight =
            -green(Tau, curr.K[2]) * green(-Tau, curr.K[2] + curr.K[1]) * SPIN * PhaseFactor
        @assert abs(weight - eval(curr.order)) < 1.0e-10 "weight wrong"
    end

    @assert isinf(factor) == false "factor is infinite at step $(curr.step)"
    @assert abs(curr.absWeight - abs(eval(curr.order))) < 1.0e-10 "weight is wrong!"

    if DiagType == POLAR || DiagType == SIGMA || DiagType == DELTA
        Observable.measure(oneBody, eval(curr.order), factor)
    end
end

function save()
    println("Saving data to disk ...")
    if DiagType == POLAR || DiagType == SIGMA || DiagType == DELTA
        Observable.save(oneBody)
    end
end

function reweight()
    return
end

@inline propose(name) = (Proposed[name, curr.order+1] += 1)
@inline accept(name) = (Accepted[name, curr.order+1] += 1)

function increaseOrder()
    curr.order == Order && return # already at the highest order
    newOrder = curr.order + 1
    prop = 1.0
    if curr.order == 0
        # create new external Tau
        # curr.extTidx, propT = createExtIdx(TauGridSize)
        # varT[LastTidx] = Grid.tau.grid[curr.extTidx]

        # curr.extKidx, propK = createExtIdx(KGridSize)
        # varK[1][1] = Grid.K.grid[curr.extKidx]
        # prop = propT * propK
    else
        # create new internal Tau
        varT[lastInnerTidx(newOrder)], prop = createTau()
    end
    # newOrder == 1 && println(lastInnerKidx(newOrder))
    # oldK = copy(varK[2])
    # prop *= createK!(varK[lastInnerKidx(newOrder)])
    prop *= createK!(varK[2])
    # @assert norm(oldK) != norm(varK[2]) "K remains the same"

    newAbsWeight = abs(eval(newOrder))
    # println(prop, ", ", newAbsWeight)
    R = prop * newAbsWeight * ReWeight[newOrder+1] / curr.absWeight / ReWeight[curr.order+1]
    propose(INCREASE_ORDER)
    if rand(rng) < R
        accept(INCREASE_ORDER)
        curr.order = newOrder
        curr.absWeight = newAbsWeight
    end
end

function decreaseOrder()
    curr.order == 0 && return
    newOrder = curr.order - 1
    prop = 1.0
    if newOrder == 0
        # remove external Tau
        # prop *= removeExtIdx(TauGridSize)
        # prop *= removeExtIdx(KGridSize)
    else
        # remove internal Tau
        prop *= removeTau()
    end
    # prop *= removeK(varK[lastInnerKidx(curr.order)])
    prop *= removeK(varK[2])
    newAbsWeight = abs(eval(newOrder))
    R = prop * newAbsWeight * ReWeight[newOrder+1] / curr.absWeight / ReWeight[curr.order+1]
    propose(DECREASE_ORDER)
    if rand(rng) < R
        accept(DECREASE_ORDER)
        curr.order = newOrder
        curr.absWeight = newAbsWeight
    end
end

function changeTau()
    # Proposed[CHANGE_TAU, curr.order + 1] += 1
    # Accepted[CHANGE_TAU, curr.order + 1] += 1
    return
end

function changeK()
    # return
    curr.order == 0 && return
    loopidx = rand(rng, firstInnerKidx():lastInnerKidx(curr.order))
    oldK = copy(varK[loopidx])
    prop = shiftK!(oldK, varK[loopidx])

    propose(CHANGE_K)
    newAbsWeight = abs(eval(curr.order))
    if rand(rng) < prop * newAbsWeight / curr.absWeight
        accept(CHANGE_K)
        curr.absWeight = newAbsWeight
    else
        varK[loopidx] .= oldK
    end
end

function changeExtTau()
    # curr.order == 0 && return
    oldTidx = curr.extTidx
    curr.extTidx, prop = shiftExtIdx(TauGridSize)
    varT[LastTidx] = Grid.tau.grid[curr.extTidx]
    propose(CHANGE_EXTTAU)
    newAbsWeight = abs(eval(curr.order))
    if rand(rng) < prop * newAbsWeight / curr.absWeight
        accept(CHANGE_EXTTAU)
        curr.absWeight = newAbsWeight
    else
        curr.extTidx = oldTidx
        varT[LastTidx] = Grid.tau.grid[curr.extTidx]
    end
end

function changeExtK()
    # curr.order == 0 && return
    oldKidx = curr.extKidx
    prop = 1.0
    if DiagType == POLAR
        curr.extKidx, prop = shiftExtIdx(KGridSize)
        varK[1][1] = Grid.K.grid[curr.extKidx]
    else
        return
    end
    propose(CHANGE_EXTK)
    newAbsWeight = abs(eval(curr.order))
    if rand(rng) < prop * newAbsWeight / curr.absWeight
        accept(CHANGE_EXTK)
        curr.absWeight = newAbsWeight
    else
        curr.extKidx = oldKidx
        varK[1][1] = Grid.K.grid[curr.extKidx]
    end
end

const barbar = "====================================================================================="
const bar = "-------------------------------------------------------------------------------------"

function printStatus()
    # Var.counter += 1
    println(barbar)
    printstyled(Dates.now(), color = :green)
    println("\nStep:", curr.step)
    println(bar)
    for i = 1:UpdateNum
        @printf("%-14s %12s %12s %12s\n", Name[i], "Proposed", "Accepted", "Ratio  ")
        for o = 0:Order
            @printf(
                "  Order%2d:     %12.0f %12.0f %12.6f\n",
                o,
                Proposed[i, o+1],
                Accepted[i, o+1],
                Accepted[i, o+1] / Proposed[i, o+1]
            )
        end
        println(bar)
    end
    println(progressBar(round(curr.step / 1000_000, digits = 2), TotalBlock))
    println()
end

end

max_stars_repo_name: kunyuan/DiagramCalculator

max_issues_repo_name: kunyuan/DiagramCalculator

max_forks_repo_name: kunyuan/DiagramCalculator


text:
using Distributions
using Random
using LinearAlgebra
using PDMats
using ProgressMeter
E = LinearAlgebra.I
eye(n) = Matrix(1.0E,n,n) #identity matrix

#import Fontconfig, Cairo
using Gadfly

using SpecialFunctions

# fundamaental setting
## calc. RLCT
SIM_ITERS = 20
LOG_FILE_PATH = joinpath("log","result.csv")
## MCMC
MCMC_K = 1000
BURNIN = 20000
THIN = 20
MCMC_ITER = BURNIN + THIN*MCMC_K
## true and model
#Random.seed!(1959)
eps = 1e-5
U_0 = [1; 2; 3; 1-eps]
V_0 = [1 9 9 1+eps]
#const A_0 = ones(10,1).*[0.1; 0.05; 0.02; 0.03; 0.15; 0.25; 0.1; 0.1; 0.15; 0.05] #[0.4 0.3; 0.1 0.3; 0.4 0.2; 0.1 0.2;]#np.array([[0.5,0.3],[0.1,0.3],[0.4,0.4]])
#const B_0 = ones(Float64,1,9) #[0.1 0.1 0.4; 0.9 0.9 0.6;] #np.array([[0.1,0.1,0.4],[0.9,0.9,0.6]])
M = 4 #size(A_0)[1] #vocab. size
N = 4 #size(B_0)[2] #num of docs
H = 2 #num of topics in learner
H_0 = 1 #size(A_0)[2] #true num of topics
## hyper parameter
PHI_U = 0.25
THETA_U = 1.0
PHI_V = 0.75
THETA_V = 1.0
## data
SEED = 1994
#Random.seed!(SEED);
L = 500
SAMPLESIZE = L
TESTSIZE = 100*SAMPLESIZE #15000 #100*SAMPLESIZE

println(U_0)
println(V_0)
println(U_0 * V_0)

function generate_matrix_from_prior(nR, nC; phi=PHI_V, theta=THETA_V, seed=10000)
    W = zeros(Float64, nR, nC);
    for i in 1:nR
        for j in 1:nC
            W[i,j] = rand(Gamma(phi, 1/theta))
        end
    end
    return W
end

function generate_data(U, V; n::Int=SAMPLESIZE, seed=1)
    #===
    for l=1:n
        generate words from Cat(x|AB[:,j]) for j=1,...,N
        as two format: docs array (N, n//N, M) and words array (n, M)
    end

    :param U: stochastic matrix A= (the probability that the word is i when the topic is k)
    :type U: shape is (M,H)
    :param V: stochastic matrix B= (the probability that the topic is k when the doc is j)
    :type V: shape is (H,N)
    :param n: the number of all words in all document. sample size
    :type n: Int
    :return: generated data matrices whose [i,j,l] entry is the (i,j)-element in the l-th matrix.
    :rtype: Array{Int,3}
    ===#
    Random.seed!(seed);
    C = U*V
    #M = size(C)[1]
    #N = size(C)[2]
    Xs_arr = zeros(Int,M,N,n)
    for i in 1:M
        for j in 1:N
            X_ij = rand(Poisson(C[i,j]), n)
            Xs_arr[i,j,:] = X_ij;
        end
    end
    return Xs_arr
end

# each sampling step for Gibbs sampler
function latent_step(Xs_arr::Array{Int64,3}, U, V; n=SAMPLESIZE)
    #===
    latent variable sampling step.

    :param Xs_arr: data
    :param U: param matrix
    :param V: param matrix
    :return: a latent variable s whose shape is (M,H,N,n)
    ===#
    #H = size(A)[2]
    #N = size(B)[2]
    #L = div(n, N)
    W = U*V #W_ij = sum_k U_ik V_kj
    π = zeros(Float64, M, H, N)
    s = zeros(Float64, M, H, N, n)
    for i in 1:M
        U_i = H>1 ? U[i,:] : U
        for j in 1:N
            V_j = H>1 ? V[:,j] : V
            π[i,:,j] = U_i.*V_j ./ W[i,j]
            for l in 1:n
                s[i,:,j,l] = rand(Multinomial(Xs_arr[i,j,l], π[i,:,j]))
            end
        end
    end
    return s
end

function U_step(s, V, phi_U, theta_U; n=SAMPLESIZE)
    #===
    matrix U sampling step.

    :param s: latent variable tensor; shape is (M, H, N, n)
    :param phi_U: hyperparameter of Gamma prior for U (RLCT depends on it)
    :param theta_U: hyperparameter of Gamma prior for U (RLCT does not depend on it)
    :return: param matrix U
    ===#
    hat_phi_U = sum(s, dims=(3,4)) .+ phi_U
    hat_theta_U = n .* sum(V, dims=2) .+ theta_U
    ## Note: Gamma dist in Julia is x^{phi-1} e^{-x/theta_julia}.
    ## Thus: our theta must be inversed for plugin: Gamma(phi, theta_julia = 1/theta ).
    U = zeros(Float64, M, H)
    for i in 1:M
        for k in 1:H
            U[i,k] = rand(Gamma(hat_phi_U[i,k,1,1], 1 ./ hat_theta_U[k,1]))
        end
    end
    return U
end

function V_step(s, U, phi_V, theta_V; n=SAMPLESIZE)
    #===
    matrix V sampling step.

    :param s: latent variable tensor; shape is (M, H, N, n)
    :param phi_V: hyperparameter of Gamma prior for V (RLCT depends on it)
    :param theta_V: hyperparameter of Gamma prior for V (RLCT does not depend on it)
    :return: param matrix U
    ===#
    hat_phi_V = sum(s, dims=(1,4)) .+ phi_V
    hat_theta_V = n .* sum(U, dims=1) .+ theta_V
    ## Note: Gamma dist in Julia is x^{phi-1} e^{-x/theta_julia}.
    ## Thus: our theta must be inversed for plugin: Gamma(phi, theta_julia = 1/theta ).
    V = zeros(Float64, H, N)
    for k in 1:H
        for j in 1:N
            V[k,j] = rand(Gamma(hat_phi_V[1,k,j,1], 1 ./ hat_theta_V[1,k]))
        end
    end
    return V
end

function model_pmf(X::Array{Int64,2}, U, V)
    #===
    probability mass function of model p(X | A,B).
    See also the above markdown cell.

    :param X: an (M,N) non-negative integer matrix which means a data.
    :param U: an (M,H) non-negative matrix.
    :param V: an (H,N) non-negative matrix.
    ===#
    #mass = 1.0
    W = U*V
    #Brute Forse
    #===for i in 1:M
        for j in 1:N
            mass *= pdf(Poisson(W[i,j]), X[i,j])
        end
    end===#
    log_mass = sum(X.*log.(W) - W - loggamma.(X.+1))
    mass = exp(log_mass)
    return mass
end

# Gibbs sampler
function run_all_sampling!(
        Xs_arr::Array{Int64,3}, init_S,
        #allSs::Array{Int64,5},
        allUs::Array{Float64,3}, allVs::Array{Float64,3})
    #===
    run Gibbs sampling before doing burn-in and thining.

    :param Xs_arr: dataset; shape is (M, N, n)
    :param init_S: initial value of the latent variable; shape is (M, H, N, n)
    :param allSs: (to be mutated!) all sample of latent variables; shape is (M, H, N, n, sampling_iters)
    :param allUs: (to be mutated!) all samle of parameter U; shape is (M, H, sampling_iters)
    :param allVs: (to be mutated!) all samle of parameter V; shape is (H, N, sampling_iters)
    :return: void (this function is mutator of allSs, allUs, and allVs.)
    ===#
    ## get the number of iterations and check array size
    iters = size(allUs)[3];
    @assert size(allVs)[3]==iters
    ## get hyperparameter
    ϕ_U = PHI_U;
    θ_U = THETA_U;
    ϕ_V = PHI_V;
    θ_V = THETA_V;
    ## define Progress object for checking iteration progress
    sampling_progress = Progress(iters);
    ## set S to the initial value
    S = init_S;
    ## set V to the initial value from the prior
    V = generate_matrix_from_prior(H, N)
    ## start all sampling
    @time for k in 1:iters
        ## parameter matrix sampling and saving
        U = U_step(S, V, ϕ_U, θ_U);
        V = V_step(S, U, ϕ_V, θ_V);
        allUs[:,:,k] = U;
        allVs[:,:,k] = V;
        ## latent variable sampling and saving
        S = latent_step(Xs_arr, U, V);
        #allSs[:,:,:,:,k] = S;
        ## Progress update
        next!(sampling_progress);
    end
    #return (allYs, allAs, allBs)
end

function run_process_posterior!(
        Xs_arr::Array{Int64,3},
        #allSs::Array{Int64,5},
        allUs::Array{Float64,3}, allVs::Array{Float64,3},
        #gsSs::Array{Int64,5},
        gsUs::Array{Float64,3}, gsVs::Array{Float64,3}, likelihoodMat::Array{Float64,2};
        K::Int=MCMC_K, burn::Int=BURNIN, th::Int=THIN, n::Int=SAMPLESIZE)
    #===
    For sampling result, this function do burn-in, thining, and calculating the likelihood matrix (for WAIC).
    ===#
    # burnin and thining
    ## check MCMC samplesize
    #@assert size(gsUs)[3]==size(gsSs)[5]
    @assert size(gsUs)[3]==size(gsVs)[3]
    @assert size(gsUs)[3]==K
    ## define Progress object.
    thining_progress = Progress(K);
    ## run thining and burn-in
    @time for k in 1:K
        U = allUs[:,:,burn+k*th];
        gsUs[:,:,k] = U;
        V = allVs[:,:,burn+k*th];
        gsVs[:,:,k] = V;
        #gsSs[:,:,:,:,k] = allSs[:,:,:,:,burn+k*th];
        ## generated quantity: likelihood matrix
        #likelihoodMat = calc_likelihoodMat!(k, words_arr, doc_onehot, A, B, likelihoodMat)
        for l in 1:n
            X = Xs_arr[:,:,l];
            likelihoodMat[l,k] = model_pmf(X,U,V);
        end
        next!(thining_progress);
    end
    #return (gsYs, gsAs, gsBs, likelihoodMat)
end

function run_Gibbs_sampler_cored(
    Xs_arr::Array{Int64,3}; K::Int=MCMC_K, burn::Int=BURNIN, th::Int=THIN, n::Int=SAMPLESIZE, seed_MCMC=2)
    #===
    Version: separating core functions
    ===#
    Random.seed!(seed_MCMC);
    ## initial value of A and B is sampled from the prior distributuion
    #### Assume that prior is SYMMETRIC Gamma distribution. Hyparam is scalar.
    init_U = generate_matrix_from_prior(M, H; phi=PHI_U, theta=THETA_U);
    init_V = generate_matrix_from_prior(H, N; phi=PHI_V, theta=THETA_V);
    init_S = latent_step(Xs_arr, init_U, init_V);
    ## sampling iteration
    #### allocate tensors for MCMC sample
    iters = burn + K*th;
    allUs = zeros(Float64, M, H, iters);
    gsUs = zeros(Float64, M, H, K);
    allVs = zeros(Float64, H, N, iters);
    gsVs = zeros(Float64, H, N, K);
    #allSs = zeros(Int64, M, H, N, n, iters);
    #gsSs = zeros(Int64, M, H, N, n, K);
    #### all sampling
    println("Start $iters iteration for GS")
    run_all_sampling!(Xs_arr, init_S, allUs, allVs)
    ## burnin and thining
    likelihoodMat = zeros(Float64, n, K)
    println("Start burn-in and thining from $iters to $K")
    run_process_posterior!(Xs_arr, allUs, allVs, gsUs, gsVs, likelihoodMat)
    return (gsUs, gsVs, likelihoodMat)
end

function run_prior_sampling(Xs_arr; K=MCMC_K, burn=BURNIN, th=THIN, n=SAMPLESIZE, seed_MCMC=2)
    #===
    Run prior sampling K times.

    :return: sample of Y(topic indicator variable),A and B from GS and likelihood matrix
    :rtype: tuple(Array{Bool,3};shape(H,n,K), Array{Float64,3};shape(M,H,K), Array{Float64,3};shape(H,N,K), Array{Float64,2};shape(n,K))
    ===#
    Random.seed!(seed_MCMC);
    ## sampling iteration
    #### allocate tensors for MCMC sample
    gsUs = zeros(Float64, M, H, K);
    gsVs = zeros(Float64, H, N, K);
    #gsSs = zeros(Int64, M, H, N, n, K);
    #### all sampling and calc. likelihood Mat
    println("Start $K iteration for prior sampling (no thining)")
    sampling_progress = Progress(K);
    likelihoodMat = zeros(Float64, n, K)
    @time for k in 1:K
        ## sampling and saving
        U = generate_matrix_from_prior(M, H; phi=PHI_U, theta=THETA_U);
        V = generate_matrix_from_prior(H, N; phi=PHI_V, theta=THETA_V);
        S = latent_step(Xs_arr, U, V);
        gsUs[:,:,k] = U;
        gsVs[:,:,k] = V;
        #gsSs[:,:,:,:,k] = S;
        ## calc likelihood matrix
        for l in 1:n
            X = Xs_arr[:,:,l];
            likelihoodMat[l,k] = model_pmf(X,U,V);
        end
        ## Progress update
        next!(sampling_progress);
    end
    return (gsUs, gsVs, likelihoodMat)
end

# define functions for calculating RLCT
function calc_functional_var(loglikeMat::Array{Float64,2})
    #===
    Calculate functional variance from loglike matrix.

    :param loglikeMat: a matrix whose (l,k) element is log p(x_l|z_l,A_k,B_k)
    ===#
    n = size(loglikeMat)[1];
    K = size(loglikeMat)[2];
    first_term = reshape(mean(loglikeMat.^2, dims=2),n);
    second_term = reshape(mean(loglikeMat, dims=2).^2, n);
    func_var = sum(first_term - second_term);
    return func_var
end

function calc_predict_dist(X, Us, Vs)
    #===
    Calculate pmf of predictive distribution.

    :param X: an (M,N) matrix which means a new data.
    :param Us: an Array whose [:,:,k] element means an MCMC sample of U.
    :param Vs: an Array whose [:,:,k] element means an MCMC sample of V.
    ===#
    @assert size(Us)[3]==size(Vs)[3];
    K = size(Vs)[3];
    mass = 0.0;
    for k in 1:K
        mass += model_pmf(X, Us[:,:,k], Vs[:,:,k]);
    end
    return mass/K
end

function calc_normalized_WAIC(Xs_arr, true_U, true_V, likelihoodMat)
    #===
    Calculating normalized WAIC.
    ===#
    n = SAMPLESIZE;
    emp_loss = -mean(log.(mean(likelihoodMat, dims=2)));
    emp_entropy = -mean([log(model_pmf(Xs_arr[:,:,l], true_U, true_V)) for l in 1:n])
    func_var = calc_functional_var(log.(likelihoodMat))
    #println(typeof(emp_entropy))
    normalized_WAIC = emp_loss + func_var/n - emp_entropy
    return normalized_WAIC
end

function calc_generalization_error(true_U, true_V, Us, Vs; nT=TESTSIZE, seed_T=3)
    #===
    Calculating generalization error using test data generated by true distribution.
    ===#
    test_Xs_arr = generate_data(true_U, true_V, n=nT, seed=seed_T)
    ge = 0.0
    calc_gerr_progress = Progress(nT)
    println("Test sample size nT = $nT")
    @time for t in 1:nT
        q = model_pmf(test_Xs_arr[:,:,t], true_U, true_V)
        pred = calc_predict_dist(test_Xs_arr[:,:,t], Us, Vs)
        ge += log(q) - log(pred) #log(q/pred), but q and pred is small number.
        next!(calc_gerr_progress)
    end
    ge /= nT
    return ge
end

function calc_hamiltonian_unif(likelihoodMat)
    #===
    Calculate hamiltonian if the prior is uniform.
    I.e. this function calclutates negative logarithm likelihood at each MCMC sample.
    ===#
    K = size(likelihoodMat)[2]
    #@assert K==size(As)[3]
    #@assert K==size(Bs)[3]
    loglikeMat = log.(likelihoodMat)
    loglike = sum(loglikeMat, dims=1)[1,:]
    @assert length(loglike) == K
    return -loglike
end

function run_single_inference(true_U, true_V, seed)
    #===
    Run single inference for debug mode.
    ===#
    Random.seed!(seed);
    train_X = generate_data(true_U, true_V, n=SAMPLESIZE, seed=seed+1);
    println("Gibbs Sampling")
    param_Us, param_Vs, likelihoodMat = run_Gibbs_sampler_cored(train_X, seed_MCMC=seed+2);
    println("Calculation Hamiltonian trace")
    hamils = calc_hamiltonian_unif(likelihoodMat)
    println("Calculation Normalized WAIC")
    normalized_WAIC = calc_normalized_WAIC(train_X, true_U, true_V, likelihoodMat);
    println("Calculaton Generalization Error")
    ge = calc_generalization_error(true_U, true_V, param_Us, param_Vs, seed_T=seed+3);
    return (hamils, ge, normalized_WAIC)
end

function run_single_inference_debug(true_U, true_V, seed)
    #===
    Run single inference for debug mode.
    ===#
    Random.seed!(seed);
    train_X = generate_data(true_U, true_V, n=SAMPLESIZE, seed=seed+1);
    println("Gibbs Sampling")
    param_Us, param_Vs, likelihoodMat = run_Gibbs_sampler_cored(train_X, seed_MCMC=seed+2);
    println("Calculation Normalized WAIC")
    normalized_WAIC = calc_normalized_WAIC(train_X, true_U, true_V, likelihoodMat);
    println("Calculaton Generalization Error")
    ge = calc_generalization_error(true_U, true_V, param_Us, param_Vs, seed_T=seed+3);
    return (param_Us, param_Vs, likelihoodMat, ge, normalized_WAIC)
end

function run_single_inference_PS(true_U, true_V, seed)
    #===
    Run single inference using prior sampling (worst benchmark).
    ===#
    Random.seed!(seed);
    train_X = generate_data(true_U, true_V, n=SAMPLESIZE, seed=seed+1);
    println("Prior Sampling")
    latent_Ss, param_Us, param_Vs, likelihoodMat = run_prior_sampling(train_X, seed_MCMC=seed+2)
    println("Calculation Normalized WAIC")
    normalized_WAIC = calc_normalized_WAIC(train_X, true_U, true_V, likelihoodMat);
    if normalized_WAIC>1e10
        normalized_WAIC = 1e10
    elseif isnan(normalized_WAIC)
        normalized_WAIC = 1e12 + 2e11
    elseif normalized_WAIC<0
        println("Warning: normalized_WAIC became negative!")
    end
    println("Calculaton Generalization Error")
    ge = calc_generalization_error(true_U, true_V, param_Us, param_Vs, seed_T=seed+3);
    if ge>1e10
        ge = 1e10
    elseif isnan(ge)
        ge = 1e12 + 2e11
    elseif ge<0
        println("Warning: Generalization Error became negative!")
    end
    return (latent_Ss, param_Us, param_Vs, likelihoodMat, ge, normalized_WAIC)
end

function run_multi_inference(true_U, true_V, sim_iters, log_file_path)
    #===
    for it in 1:sim_iters
        run_single_inference
    end
    and calculate RLCT from above $sim_iters simulation results.

    :param true_U: true parameter non-negative matrix U
    :param true_V: true parameter non-negative matrix V
    :param sim_iters: number of simulations
    :param log_file_path: path of the file to write the experimental log
    :return: generalization errors and normalized WAICs in the simulations
    :rtype: tuple(Array{Float64,1}, Array{Float64,1})
    ===#

    seeds = SEED .+ SEED * (1:sim_iters)
    gerrors = zeros(Float64, sim_iters)
    normalized_WAICs = zeros(Float64, sim_iters)
    hamiltonian_vecs = zeros(Float64, MCMC_K, sim_iters)
    open(log_file_path, "a") do fp
        println(fp, "## Simulation Setting")
        println(fp, "M,N,H,H_0,PHI_U,THETA_U,PHI_V,THETA_V,SAMPLESIZE,MCMC_K,BURNIN,THIN,TESTSIZE,SIM_ITERS,SEED")
        println(fp, "$M,$N,$H,$H_0,$PHI_U,$THETA_U,$PHI_V,$THETA_V,$SAMPLESIZE,$MCMC_K,$BURNIN,$THIN,$TESTSIZE,$sim_iters,$SEED")
        println(fp, "## Simulation Log")
        println(fp, "iter,gerror,normalized_WAIC,RLCT,RLCT_SEM")
    end
    simulation_progress = Progress(sim_iters)
    for it in 1:sim_iters
        println("# start $it th simulation")
        hamils, ge, norm_W = run_single_inference(true_U, true_V, seeds[it])
        #plot(x=1:MCMC_K, y=hamils, Geom.line, Geom.point)
        hamiltonian_vecs[:,it] = hamils
        gerrors[it] = ge
        normalized_WAICs[it] = norm_W
        now_ges = gerrors[1:it]
        now_nwaics = normalized_WAICs[1:it]
        rlct = (SAMPLESIZE/2)*mean(now_ges + now_nwaics)
        rlct_sem = (SAMPLESIZE/2)*std(now_ges + now_nwaics)/sqrt(it)
        open(log_file_path, "a") do fp
            println(fp,"$it,$ge,$norm_W,$rlct,$rlct_sem")
            println("$it,$ge,$norm_W,$rlct,$rlct_sem")
        end
        next!(simulation_progress)
        #@assert ge>=0
        #@assert norm_W>=0
    end
    return (hamiltonian_vecs, gerrors, normalized_WAICs, seeds)
end

# main process
#===#
hamivecs, gerrors, normalized_WAICs, seeds = @time run_multi_inference(U_0, V_0, SIM_ITERS, LOG_FILE_PATH)
ge_and_nW = gerrors + normalized_WAICs
each_lams = (SAMPLESIZE/2) .* ge_and_nW
lam = (SAMPLESIZE/2)*mean(ge_and_nW)
lam_sd = (SAMPLESIZE/2)*std(ge_and_nW)/sqrt(SIM_ITERS)
println("The RLCT is $lam ± $lam_sd")
#===#

each_lams_mean = mean(each_lams)
each_lams_std = std(each_lams)
upperline = each_lams_mean + each_lams_std
middleline = each_lams_mean
lowerline = each_lams_mean - each_lams_std
plot(y=each_lams, yintercept = [upperline, middleline, lowerline], Geom.violin, Geom.hline(color=["orange","red","orange"], style=[:dash,:solid,:dash]))

lam2 = SAMPLESIZE*mean(gerrors)

lam3 = SAMPLESIZE*mean(normalized_WAICs)

println("$M, $N, $H, $H_0, $PHI_U, $PHI_V")
lam_ub = (H_0*(M+N-1) + (H-H_0)*min(PHI_U*M, PHI_V*N))/2

#Phase Transition Line in VBNMF
phase = M*PHI_U + N*PHI_V < (M+N)/2

lam_reg = (M+N)*H/2

lam_reg_rrr = (M*H + H*N - H*H)/2

#exact lambda_vb
if phase
    lam_vb = (H-H_0)*(M*PHI_U + N*PHI_V) + (M+N)*H_0/2
else
    lam_vb = lam_reg
end
lam_vb

#memo: script for convergence diagnox from visual.

#===
for i in 1:M
    for k in 1:H
        p = plot(x=debug_As[i,k,:], Geom.histogram(bincount=12),
                Guide.title("posthist20191226_A[$i.$k]"))
        img = PNG(joinpath("img_MSI","20191226-d1","posthist_A_$i-$k.png"), 16cm, 16cm)
        draw(img,p)
        p = plot(x=1:MCMC_K, y=debug_As[i,k,:], Geom.point, Geom.line,
                Guide.title("posttrace20191226_A[$i.$k]"))
        img = PNG(joinpath("img_MSI","20191226-d1","posttrace_A_$i-$k.png"), 16cm, 16cm)
        draw(img,p)
    end
end

for k in 1:H
    for j in 1:N
        p = plot(x=debug_Bs[k,j,:], Geom.histogram(bincount=12),
                Guide.title("posthist20191226_B[$k.$j]"))
        img = PNG(joinpath("img_MSI","20191226-d1","posthist_B_$k-$j.png"), 16cm, 16cm)
        draw(img,p)
        p = plot(x=1:MCMC_K, y=debug_Bs[k,j,:], Geom.point, Geom.line,
                Guide.title("posttrace20191226_B[$k.$j]"))
        img = PNG(joinpath("img_MSI","20191226-d1","posttrace_B_$k-$j.png"), 16cm, 16cm)
        draw(img,p)
    end
end
===#


text:
@testset "Forcings" begin
    fx = zeros(5,5)
    fy = zeros(5,5)
    f1 = zeros(30)
    @testset "Slippage" begin
        # No velocities
        Swalbe.slippage!(fx, fy, ones(5,5), zeros(5,5), zeros(5,5), 1.0, 1/6)
        @test all(fx .== 0.0)
        @test all(fy .== 0.0)
        # Velocity in x
        Swalbe.slippage!(fx, fy, ones(5,5), fill(0.1,5,5), zeros(5,5), 1.0, 1/6)
        @test all(fx .== 0.1/11)
        @test all(fy .== 0.0)
        # Velocity in y
        Swalbe.slippage!(fx, fy, ones(5,5), zeros(5,5), fill(0.1,5,5), 1.0, 1/6)
        @test all(fx .== 0.0)
        @test all(fy .== 0.1/11)
        # Velocity
        Swalbe.slippage!(fx, fy, ones(5,5), fill(-0.1,5,5), fill(0.1,5,5), 1.0, 1/6)
        @test all(fx .== -0.1/11)
        @test all(fy .== 0.1/11)
        # No slip
        Swalbe.slippage!(fx, fy, ones(5,5), fill(-0.1,5,5), fill(0.1,5,5), 0.0, 1/6)
        @test all(fx .== -0.1/2)
        @test all(fy .== 0.1/2)
    end
    @testset "Slippage 1D" begin
        # No velocities
        Swalbe.slippage!(f1, ones(30), zeros(30), 1.0, 1/6)
        @test all(f1 .== 0.0)
        # Velocity in x
        Swalbe.slippage!(f1, ones(30), fill(0.1,30), 1.0, 1/6)
        @test all(f1 .== 0.1/11)
        # No slip
        Swalbe.slippage!(f1, ones(30), fill(-0.1,30), 0.0, 1/6)
        @test all(f1 .== -0.1/2)
    end

    @testset "Thermal" begin
        f1 = ones(50,50)
        f2 = ones(50,50)
        vartest = 2*0.01/11
        Swalbe.thermal!(f1, f2, ones(50,50), 0.01, 1/6, 1.0)
        @test mean(f1) ≈ 0.0 atol=1e-2
        @test mean(f2) ≈ 0.0 atol=1e-2
        @test var(f1) ≈ vartest atol=vartest/10
        @test var(f2) ≈ vartest atol=vartest/10
        vartest = 0.2/11
        Swalbe.thermal!(f1, f2, ones(50,50), 0.1, 1/6, 1.0)
        @test mean(f1) ≈ 0.0 atol=1e-2
        @test mean(f2) ≈ 0.0 atol=1e-2
        @test var(f1) ≈ vartest atol=vartest/10
        @test var(f2) ≈ vartest atol=vartest/10
    end
    @testset "Thermal 1D" begin
        f1D = ones(100000)
        vartest = 2*0.01/11
        Swalbe.thermal!(f1D, ones(100000), 0.01, 1/6, 1.0)
        @test mean(f1D) ≈ 0.0 atol=1e-2
        @test var(f1D) ≈ vartest atol=vartest/10
        vartest = 0.2/11
        Swalbe.thermal!(f1D, ones(100000), 0.1, 1/6, 1.0)
        @test mean(f1D) ≈ 0.0 atol=1e-2
        @test var(f1D) ≈ vartest atol=vartest/10
    end

    @testset "rho update" begin
        @testset "Constant fields" begin
            rho = ones(25)
            height = ones(25)
            output = zeros(25)
            Swalbe.update_rho!(rho, output, height, zeros(25,2), zeros(25,4))
            @test all(rho .== 1)
        end
    end

    @testset "view four" begin
        dummy = reshape(collect(1:20),5,4)
        d1, d2, d3, d4 = Swalbe.view_four(dummy)
        @test all(d1 .== dummy[:,1])
        @test all(d2 .== dummy[:,2])
        @test all(d3 .== dummy[:,3])
        @test all(d4 .== dummy[:,4])
    end
end


text:
struct QuerySynthesisKNNOracle <: Oracle
    data
    labels::Vector{Symbol}
    k::Int
    dist_func::Function
end

function QuerySynthesisKNNOracle(data::Array{T, 2}, labels::Vector{Symbol}, params::Dict{Symbol, Any}=Dict{Symbol, Any}()) where T <: Real
    k = get(params, :k, 1)
    k % 2 == 0 && throw(ArgumentError("Even value k = $k not allowed in QuerySynthesisKNNOracle to prevent decision ties."))
    dist_func = eval(get(params, :dist_func, :euclidean))
    return QuerySynthesisKNNOracle(data, labels, k, dist_func)
end

function ask_oracle(oracle::QuerySynthesisKNNOracle, query_objects::Array{T, 2})::Vector{Symbol} where T <: Real
    function ask_oracle(q)
        sorted_distances = sort([(i, oracle.dist_func(q, oracle.data[:, i]), oracle.labels[i]) for i in 1:size(oracle.data, 2)], by = x -> x[2])[1:oracle.k]
        return count([x[3] for x in sorted_distances] .== :inlier) > length(sorted_distances) * 0.5 ? :inlier : :outlier
    end
    return [ask_oracle(query_objects[:, i]) for i in 1:size(query_objects, 2)]
end


text:
adj = [0 1 0 1;
       1 0 1 0;
       0 1 0 1;
       1 0 1 0]
nf = rand(3, 4)

@testset "cuda" begin
    fg = FeaturedGraph(cu(adj))
    @test has_graph(fg)
    @test has_node_feature(fg) == false
    @test has_edge_feature(fg) == false
    @test has_global_feature(fg) == false
    @test typeof(graph(fg)) == CuMatrix{Int64}
    @test typeof(node_feature(fg)) == Fill{Int64,2,Tuple{Base.OneTo{Int64},Base.OneTo{Int64}}}
    @test typeof(edge_feature(fg)) == Fill{Int64,2,Tuple{Base.OneTo{Int64},Base.OneTo{Int64}}}
    @test typeof(global_feature(fg)) == Fill{Int64,1,Tuple{Base.OneTo{Int64}}}

    fg = FeaturedGraph(adj, cu(nf))
    @test has_graph(fg)
    @test has_node_feature(fg)
    @test has_edge_feature(fg) == false
    @test has_global_feature(fg) == false
    @test typeof(graph(fg)) == CuMatrix{Float32}
    @test typeof(node_feature(fg)) == CuMatrix{Float32}
    @test typeof(edge_feature(fg)) == Fill{Float32,2,Tuple{Base.OneTo{Int64},Base.OneTo{Int64}}}
    @test typeof(global_feature(fg)) == Fill{Float32,1,Tuple{Base.OneTo{Int64}}}
end


text:
using MinimaxCalibratedEBayes
using Test

include("run_steinminimax.jl")
include("run_hermite.jl")


text:
@doc raw"""
    Vector3

An object `v` of the type [`Vector3`](@ref) is a vector in the three-dimensional
vector space ``ℝ^3``.

Geometrically spoken, a vector is the difference of two points ``p_2, p_1`` and
and denotes the direction and the distance from ``p_1`` to ``p_2``.

CGAL defines a symbolic constant [`NULL_VECTOR`](@ref). We will explicitly
state where you can pass this constant as an argument instead of a vector
initialized with zeros.
"""
Vector3

"""
    Vector3(a::Point3, b::Point3)

Introduces the vector `b - a`.
"""
Vector3(a::Point3, b::Point3)

"""
    Vector3(s::Segment3)

Introduces the vector `target(s) - source(s)`.
"""
Vector3(s::Segment3)

"""
    Vector3(::NullVector = NULL_VECTOR)

Introduces a null vector `v`.
"""
Vector3(::NullVector = NULL_VECTOR)

"""
    Vector3(x::Real, y::Real, z::Real)

Introduces a vector `v` initialized to ``(x, y, z)``.
"""
Vector3(x::Real, y::Real, z::Real) =
    Vector3(convert(FT, x), convert(FT, y), convert(FT, z))

@doc raw"""
    Vector3(hx::Real, hy::Real, hz::Real, hw::Real = 1)

Introduces a vector `v` initialized to `(hx/hw,hy/hw,hz/hw)`.

!!! info "Precondition"

    `hw` ``≠`` `0`
"""
Vector3(hx::Real, hy::Real, hz::Real, hw::Real) =
    Vector3(convert(RT, x), convert(RT, y), convert(RT, z), convert(RT, hw))

"""
    hx(p::Vector3)

Returns the homogeneous ``x`` coordinate.
"""
hx(p::Vector3)

"""
    hy(p::Vector3)

Returns the homogeneous ``y`` coordinate.
"""
hy(p::Vector3)

"""
    hz(p::Vector3)

Returns the homogeneous ``z`` coordinate.
"""
hz(p::Vector3)

"""
    hw(p::Vector3)

Returns the homogenizing coordinate.
"""
hw(p::Vector3)

"""
    x(p::Vector3)

Returns the Cartesian ``x`` coordinate, that is [`hx()`](@ref)/[`hw()`](@ref).
"""
x(p::Vector3)

"""
    y(p::Vector3)

Returns the Cartesian ``y`` coordinate, that is [`hy()`](@ref)/[`hw()`](@ref).
"""
y(p::Vector3)

"""
    z(p::Vector3)

Returns the Cartesian ``z`` coordinate, that is [`hz()`](@ref)/[`hw()`](@ref).
"""
z(p::Vector3)

@doc raw"""
    homogeneous(v::Vector3, i::Integer)

Returns the `i`ᵗʰ homogeneous coordinate of `v`.

!!! info "Precondition"

    ``0 ≤ i ≤ 3``.
"""
homogeneous(v::Vector3, i::Integer)

@doc raw"""
    cartesian(v::Vector3, i::Integer)

Returns the `i`ᵗʰ Cartesian coordinate of `v`.

!!! info "Precondition"

    ``0 ≤ i ≤ 2``.
"""
cartesian(v::Vector3, i::Integer)

"""
    dimension(v::Vector3)

Returns the dimension (the constant 3).
"""
dimension(v::Vector3)

"""
    squared_length(v::Vector3)

Returns the squared length of `v`.
"""
squared_length(v::Vector3)

"""
    ==(v::Vector3, w::Vector3)

Test for equality: two vectors are equal, iff their ``x``, ``y`` and ``z``
coordinates are equal.

You can compare a vector with the `NULL_VECTOR`.
"""
==(v::Vector3, w::Vector3)

"""
    +(v::Vector3, w::Vector3)

Addition.
"""
+(v::Vector3, w::Vector3)

"""
    -(v::Vector3, w::Vector3)

Subtraction.
"""
-(v::Vector3, w::Vector3)

"""
    -(v::Vector3)

Returns the opposite vector.
"""
-(v::Vector3)

"""
    *(v::Vector3, w::Vector3)

Returns the scalar product (= inner product) of the two vectors.
"""
*(v::Vector3, w::Vector3)

"""
    /(v::Vector3, s::Real)

Division by a scalar.
"""
/(v::Vector3, s::Real)

@cxxdereference Base.:/(v::Vector3, s::Real) = v / convert(RT, s)


text:
module MiniBee

using LinearAlgebra, RobotDynamics, TrajectoryOptimization, StaticArrays, Rotations, Parameters, Altro, Sockets, JSON
import RobotDynamics: mass, inertia, forces, moments

include("inertia.jl")
include("spacecraft.jl")
include("waypoints.jl")
include("costs.jl")
include("problem.jl")
include("postprocessing.jl")

"""
Build a problem, solve it and return a `SpacecraftTrajectory` to be used in postprocessing.
"""
function main(model, opts, N, tf, dt, xf_pos, u0; waypoints=[])

    prob = build_problem(model, N, tf, dt, xf_pos, u0; waypoints)

    solver = ALTROSolver(prob, opts)
    solve!(solver)
    x = states(solver)
    u = controls(solver)
    @show size(x)
    @show size(u)
    x_mat, u_mat = hcat(Array.(x)...), hcat(Array.(u)...)
    postprocess(solver, x_mat, u_mat)
    println(size(u_mat))
    time = vcat(0, cumsum(dt))
    SpacecraftTrajectory(x_mat, u_mat, time)
end

export main,
       Spacecraft,
       postprocess,
       look_for_port,
       Waypoint

end # module


text:
"""
Independently Recurrent Neural Networks, i.e. ⤦\n
    hᵗ = f(w*xᵗ + u .* hᵗ⁻¹ + b), where u is a vector
"""
mutable struct IndRNN <: Block
    w::VarOrNil # input to hidden weights
    b::VarOrNil # bias of hidden units
    u::VarOrNil # recurrent weights
    f::Function # activation function
    h::Any      # hidden variable
    function IndRNN(inputSize::Int, hiddenSize::Int, fn::Function=relu; type::Type=Array{Float32})
        T = eltype(type)
        w = randn(T, hiddenSize, inputSize) .* sqrt( T(2/inputSize) )
        b = zeros(T, hiddenSize, 1)
        u = uniform(T, (hiddenSize, 1); from=-0.1990, to=0.1997)
        new(Variable{type}(w,true,true,true),
            Variable{type}(b,true,true,true),
            Variable{type}(u,true,true,true), fn, nothing)
    end
    function IndRNN(fn::Function)
        new(nothing, nothing, nothing, fn, nothing)
    end
end

function clone(this::IndRNN; type::Type=Array{Float32})
    cloned = IndRNN(this.f)
    cloned.w = clone(this.w, type=type)
    cloned.b = clone(this.b, type=type)
    cloned.u = clone(this.u, type=type)
    return cloned
end

mutable struct IndRNNs <: Block
    layers::Vector{IndRNN}
    function IndRNNs(topology::Vector{Int}, fn::Array{F}; type::Type=Array{Float32}) where F
        n = length(topology) - 1
        layers = Vector{IndRNN}(undef, n)
        for i = 1:n
            layers[i] = IndRNN(topology[i], topology[i+1], fn[i]; type=type)
        end
        new(layers)
    end
end

Base.getindex(m::IndRNNs,     k...) =  m.layers[k...]
Base.setindex!(m::IndRNNs, v, k...) = (m.layers[k...] = v)
Base.length(m::IndRNNs)       = length(m.layers)
Base.lastindex(m::IndRNNs)    = length(m.layers)
Base.firstindex(m::IndRNNs)   = 1
Base.iterate(m::IndRNNs, i=firstindex(m)) = i>length(m) ? nothing : (m[i], i+1)

function Base.show(io::IO, m::IndRNN)
    SIZE = size(m.w)
    TYPE = typeof(m.w.value)
    print(io, "IndRNN($(SIZE[2]), $(SIZE[1]), $(m.f); type=$TYPE)")
end

function Base.show(io::IO, m::IndRNNs)
    print(io, "IndRNNs\n      (\n          ")
    join(io,      m.layers, ",\n          ")
    print(io,                   "\n      )")
end

function resethidden(m::IndRNN)
    m.h = nothing
end

function resethidden(model::IndRNNs)
    for m in model
        resethidden(m)
    end
end

function forward(m::IndRNN, x::Variable{T}) where T
    f = m.f  # activition function
    w = m.w  # input's weights
    b = m.b  # input's bias
    u = m.u  # memory's weights
    h = m.h ≠ nothing ? m.h : Variable{T}(Zeros(T, size(w,1), size(x,2)))
    x = f(matAddVec(matMulVec(h, u) + w*x, b))
    m.h = x
    return x
end

function forward(model::IndRNNs, x::Variable)
    for m in model
        x = forward(m, x)
    end
    return x
end

function predict(m::IndRNN, x::T) where T
    f = m.f        # activition function
    w = m.w.value  # input's weights
    b = m.b.value  # input's bias
    u = m.u.value  # memory's weights
    h = m.h ≠ nothing ? m.h : Zeros(T, size(w,1), size(x,2))
    x = f(w*x + h .* u .+ b)
    m.h = x
    return x
end

function predict(model::IndRNNs, x)
    for m in model
        x = predict(m, x)
    end
    return x
end

"""
    unbiasedof(m::IndRNN)

unbiased weights of IndRNN block
"""
function unbiasedof(m::IndRNN)
    weights = Vector(undef, 1)
    weights[1] = m.w.value
    return weights
end

function weightsof(m::IndRNN)
    weights = Vector(undef,3)
    weights[1] = m.w.value
    weights[2] = m.b.value
    weights[3] = m.u.value
    return weights
end

"""
    unbiasedof(model::IndRNNs)

unbiased weights of IndRNNs block
"""
function unbiasedof(model::IndRNNs)
    weights = Vector(undef, 0)
    for m in model
        append!(weights, unbiasedof(m))
    end
    return weights
end

function weightsof(m::IndRNNs)
    weights = Vector(undef,0)
    for i = 1:length(m)
        append!(weights, weightsof(m[i]))
    end
    return weights
end

function gradsof(m::IndRNN)
    grads = Vector(undef,3)
    grads[1] = m.w.delta
    grads[2] = m.b.delta
    grads[3] = m.u.delta
    return grads
end

function gradsof(m::IndRNNs)
    grads = Vector(undef,0)
    for i = 1:length(m)
        append!(grads, gradsof(m[i]))
    end
    return grads
end

function zerograds!(m::IndRNN)
    for v in gradsof(m)
        v .= zero(v)
    end
end

function zerograds!(m::IndRNNs)
    for v in gradsof(m)
        v .= zero(v)
    end
end

function paramsof(m::IndRNN)
    params = Vector{Variable}(undef,3)
    params[1] = m.w
    params[2] = m.b
    params[3] = m.u
    return params
end

function xparamsof(m::IndRNN)
    xparams = Vector{XVariable}(undef,3)
    xparams[1] = ('w', m.w)
    xparams[2] = ('b', m.b)
    xparams[3] = ('u', m.u)
    return xparams
end

function paramsof(m::IndRNNs)
    params = Vector{Variable}(undef,0)
    for i = 1:length(m)
        append!(params, paramsof(m[i]))
    end
    return params
end

function xparamsof(m::IndRNNs)
    xparams = Vector{XVariable}(undef,0)
    for i = 1:length(m)
        append!(xparams, xparamsof(m[i]))
    end
    return xparams
end

function nparamsof(m::IndRNN)
    lw = length(m.w)
    lb = length(m.b)
    lu = length(m.u)
    return (lw + lb + lu)
end

function bytesof(model::IndRNN, unit::String="MB")
    n = nparamsof(model) * elsizeof(model.w)
    return blocksize(n, uppercase(unit))
end

function nparamsof(m::IndRNNs)
    num = 0
    for i = 1:length(m)
        num += nparamsof(m[i])
    end
    return num
end

function bytesof(model::IndRNNs, unit::String="MB")
    n = nparamsof(model) * elsizeof(model[1].w)
    return blocksize(n, uppercase(unit))
end

function to(type::Type, m::IndRNN)
    m.w = to(type, m.w)
    m.b = to(type, m.b)
    m.u = to(type, m.u)
    return m
end

function to!(type::Type, m::IndRNN)
    m = to(type, m)
    return nothing
end

function to(type::Type, m::IndRNNs)
    for layer in m
        layer = to(type, layer)
    end
    return m
end

function to!(type::Type, m::IndRNNs)
    for layer in m
        to!(type, layer)
    end
end

max_stars_repo_name: sonosole/Delta

max_issues_repo_name: sonosole/Delta

max_forks_repo_name: sonosole/Delta


text:
using ReachabilityAnalysis, Plots

@taylorize function laubloomis!(dx, x, params, t)
    dx[1] = 1.4*x[3] - 0.9*x[1]
    dx[2] = 2.5*x[5] - 1.5*x[2]
    dx[3] = 0.6*x[7] - 0.8*(x[2]*x[3])
    dx[4] = 2 - 1.3*(x[3]*x[4])
    dx[5] = 0.7*x[1] - (x[4]*x[5])
    dx[6] = 0.3*x[1] - 3.1*x[6]
    dx[7] = 1.8*x[6] - 1.6*(x[2]*x[7])
    return dx
end

function laubloomis(; W=0.01)
    # initial states
    X0c = [1.2, 1.05, 1.5, 2.4, 1.0, 0.1, 0.45]
    X0 = Hyperrectangle(X0c, fill(W, 7))

    # initial-value problem
    prob = @ivp(x' = laubloomis!(x), dim: 7, x(0) ∈ X0)

    return prob
end

prob = laubloomis(W=0.01)
alg = TMJets(abstol=1e-11, orderT=7, orderQ=1, adaptive=true);

sol_1 = solve(prob, T=20.0, alg=alg);

sol_1z = overapproximate(sol_1, Zonotope);

# canonical direction along x₄
const e4 = [0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0];

ρ(e4, sol_1z)

ρ(e4, sol_1z) < 4.5

@show ρ(e4, sol_1z[end]) + ρ(-e4, sol_1z[end])

prob = laubloomis(W=0.05)
alg = TMJets(abstol=1e-12, orderT=7, orderQ=1, adaptive=false);

sol_2 = solve(prob, T=20.0, alg=alg)
sol_2z = overapproximate(sol_2, Zonotope);

ρ(e4, sol_2z)

ρ(e4, sol_2z) < 4.5

# width of final box
@show ρ(e4, sol_2z[end]) + ρ(-e4, sol_2z[end])

prob = laubloomis(W=0.1)
alg = TMJets(abstol=1e-12, orderT=7, orderQ=1, adaptive=false);

sol_3 = solve(prob, T=20.0, alg=alg)
sol_3z = overapproximate(sol_3, Zonotope);

ρ(e4, sol_3z)

ρ(e4, sol_3z) < 5.0

# width of final box
@show ρ(e4, sol_3z[end]) + ρ(-e4, sol_3z[end])

fig = plot()
plot!(fig, sol_3z, vars=(0, 4), linecolor="green", color=:green, alpha=0.8)
plot!(fig, sol_2z, vars=(0, 4), linecolor="blue",  color=:blue, alpha=0.8)
plot!(fig, sol_1z, vars=(0, 4), linecolor="yellow", color=:yellow, alpha=0.8,
      tickfont=font(10, "Times"), guidefontsize=15,
      xlab=L"t",
      ylab=L"x_4",
      xtick=[0., 5., 10., 15., 20.], ytick=[1.5, 2., 2.5, 3., 3.5, 4., 4.5, 5.],
      xlims=(0., 20.), ylims=(1.5, 5.02),
      bottom_margin=6mm, left_margin=2mm, right_margin=4mm, top_margin=3mm,
      size=(600, 600))

plot!(fig, x->x, x->4.5, 0., 20., line=2, color="red", linestyle=:dash, legend=nothing)
plot!(fig, x->x, x->5., 0., 20., line=2, color="red", linestyle=:dash, legend=nothing)


text:
module X86_64
  using ..Registers
  abstract RegisterSet <: Registers.RegisterSet
  using MachO
  import Base: copy

  immutable X86_64Arch <: Registers.Architecture
  end
  Registers.intptr(::X86_64Arch) = UInt64

  # See http://www.x86-64.org/documentation/abi-0.99.7.pdf
  const dwarf_numbering = Dict{Int, Symbol}(
  0  => :rax,   1  => :rdx,   2  => :rcx,
  3  => :rbx,   4  => :rsi,   5  => :rdi,
  6  => :rbp,   7  => :rsp,
  (i => Symbol("r$i") for i = 8:15)...,
  16 => :rip,
  (17+i => Symbol("xmm$i") for i = 0:15)...,
  (33+i => Symbol("st$i") for i = 0:7)...,
  (41+i => Symbol("mm$i") for i = 0:7)...,
  49 => :rflags,
  50 => :es, 51 => :cs, 52 => :ss,
  53 => :ds, 54 => :fs, 55 => :gs,
  58 => :fs_base, 59 => :gs_base,
  62 => :tr, 63 => :ldtr, 64 => :mxcsr, 65 => :fcw,
  66 => :fsw)
  const inverse_dwarf = map(p->p[2]=>p[1], dwarf_numbering)
  const basic_regs = 0:16
  const extended_registers = 17:32

  const gdb_numbering = Dict{Int, Symbol}(
     0 => :rax, 1 => :rbx, 2 => :rcx, 3 => :rdx,
     4 => :rsi, Dict(i=>dwarf_numbering[i] for i in 5:16)...,
    17 => :rflags, 18 => :cs, 19 => :ds, 20 => :es, 21 => :fs, 22 => :gs,
    #Dict(23+i => Symbol("st$i") for i = 0:7)...,
    31 => :fctrl, 32 => :fstat, 33 => :ftag, 34 => :fiseg, 35 => :fioff,
    36 => :foseg, 37 => :fooff, 38 => :fop,
    Dict(39+i => Symbol("xmm$i") for i = 0:15)...
  )
  const inverse_gdb = map(p->p[2]=>p[1], gdb_numbering)

  const seh_numbering = Dict{Int, Symbol}(
     0 => :rax, 1 => :rcx, 2 => :rdx, 3 => :rbx, 4 => :rsp, 5 => :rbp,
     6 => :rsi, 7 => :rdi, (i => Symbol("r$i") for i = 8:15)...
  )
  const inverse_seh = map(p->p[2]=>p[1], seh_numbering)

  const kernel_order = [
    [Symbol("r$r") for r in 15:-1:12]; :rbp; :rbx;
    [Symbol("r$r") for r in 11:-1:8]; :rax; :rcx;
    :rdx; :rsi; :rdi; :orig_rax; :rip; :cs; :eflags;
    :rsp; :ss; :fs_base; :gs_base; :ds; :es; :fs; :gs
  ]
  const kernel_numbering = Dict(sym=>idx for (idx,sym) in enumerate(kernel_order))

  # This operation can be performance critical, precompute it.
  const dwarf2gdbmap = [inverse_gdb[dwarf_numbering[regno]] for regno in basic_regs]
  dwarf2gdb(regno) = regno in basic_regs ? dwarf2gdbmap[regno+1] :
    inverse_gdb[dwarf_numbering[regno]]

  # Basic Register Set
  const RegT = RegisterValue{UInt64}
  const ExtendedRegT = RegisterValue{UInt128}
  @eval type BasicRegs <: RegisterSet
      $(Expr(:block, (
        :($(dwarf_numbering[i]) :: RegT) for i in basic_regs
      )...))
      BasicRegs() = new()
  end
  type ExtendedRegs <: RegisterSet
      basic::BasicRegs
      xsave_state::NTuple{832, UInt8}
  end
  Registers.getarch(::Union{BasicRegs,ExtendedRegs}) = X86_64Arch()

  import ..Registers: ip, set_ip!, set_sp!, invalidate_regs!

  ip(regs::BasicRegs) = regs.rip
  set_ip!(regs::BasicRegs, ip) = regs.rip = RegisterValue{UInt64}(ip)
  set_sp!(regs::BasicRegs, sp) = regs.rsp = RegisterValue{UInt64}(sp)
  function invalidate_regs!(regs::BasicRegs)
      for fieldi = 1:nfields(regs)
          setfield!(regs, fieldi, Registers.invalidated(getfield(regs, fieldi)))
      end
  end

  for f in (:ip, :set_ip!, :set_sp!, :invalidate_regs!)
    @eval $(f)(regs::ExtendedRegs, args...) = $(f)(regs.basic, args...)
  end

  copy(regs::ExtendedRegs) = ExtendedRegs(copy(regs.basic),regs.xsave_state)
  function copy(regs::BasicRegs)
      ret = BasicRegs()
      for i = 1:nfields(regs)
          setfield!(ret, i, getfield(regs, i))
      end
      ret
  end

  function Base.show(io::IO, regs::BasicRegs)
      for (i,reg) in enumerate(fieldnames(typeof(regs)))
          println(io," "^(3-length(string(reg))),reg," ",getfield(regs, i))
      end
  end

  function Registers.set_dwarf!(regs::BasicRegs, reg::Integer, value)
      (reg <= last(basic_regs)) && setfield!(regs, reg+1, RegisterValue{UInt64}(value))
  end

  function Registers.get_dwarf(regs::BasicRegs, reg::Integer)
      (reg <= last(basic_regs)) ? getfield(regs, reg+1) :
        RegisterValue{UInt64}(0, 0)
  end

  function Registers.get_dwarf(regs::ExtendedRegs, reg::Integer)
      if reg <= last(basic_regs)
        get_dwarf(regs.basic, reg)
      elseif reg <= last(extended_registers)
        # XMM registers
        startoff = 160+16*(reg - first(extended_registers))
        endoff = startoff + 16
        reinterpret(UInt128,
          reinterpret(UInt8,[regs.xsave_state])[startoff:endoff])[]
      else
        error("Extraction not yet implemented for this register")
      end
  end

  function Registers.set_dwarf!(regs::ExtendedRegs, reg::Integer, value)
    if reg <= last(basic_regs)
      set_dwarf!(regs.basic, reg, value)
    else
      error("Extraction not yet implemented for this register")
    end
  end

  function Registers.get_dwarf(::X86_64Arch, regs, sym)
    get_dwarf(regs, inverse_dwarf[sym])
  end

  function Registers.set_dwarf!(::X86_64Arch, RC, sym, val)
    set_dwarf!(RC, inverse_dwarf[sym], val)
  end

  function Registers.get_syscallarg(::X86_64Arch, regs, idx)
    @assert 1 <= idx <= 6
    get_dwarf(regs, [:rdi, :rsi, :rdx, :r10, :r8, :r9][idx])
  end

  const state_64_regs = [:rax, :rbx, :rcx, :rdx, :rdi, :rsi, :rbp, :rsp,
    (Symbol("r$i") for i = 8:15)..., :rip, #= :rflags, :cs, :fs, :gs =#]
  function BasicRegs(thread::MachO.thread_command)
    RC = BasicRegs()
    if thread.flavor == MachO.x86_THREAD_STATE64
      for (i,reg) in enumerate(state_64_regs)
        set_dwarf!(RC, inverse_dwarf[reg], thread.data[i])
      end
    else
      error("Unknown flavor")
    end
    return RC
  end

  include("instemulate.jl")
end


text:
export map

import Base: map
import Base: show

"""
    map(::Type{R}, mappingFn::F) where { F <: Function }

Creates a map operator, which applies a given `mappingFn` to each value emmited by the source
Observable, and emits the resulting values as an Observable. You have to specify output R type after
`mappingFn` projection.

# Arguments
- `::Type{R}`: the type of data of transformed value, may be or may not be the same as source type
- `mappingFn::Function`: transformation function with `(data::L) -> R` signature, where L is type of data in input source

# Producing

Stream of type `<: Subscribable{R}`

# Examples
```jldoctest
using Rocket

source = from([ 1, 2, 3 ])
subscribe!(source |> map(Int, (d) -> d ^ 2), logger())
;

# output

[LogActor] Data: 1
[LogActor] Data: 4
[LogActor] Data: 9
[LogActor] Completed
```

See also: [`AbstractOperator`](@ref), [`RightTypedOperator`](@ref), [`ProxyObservable`](@ref), [`logger`](@ref)
"""
map(::Type{R}, mappingFn::F) where { R, F <: Function } = MapOperator{R, F}(mappingFn)

struct MapOperator{R, F} <: RightTypedOperator{R}
    mappingFn::F
end

function on_call!(::Type{L}, ::Type{R}, operator::MapOperator{R, F}, source) where { L, R, F }
    return proxy(R, source, MapProxy{L, F}(operator.mappingFn))
end

struct MapProxy{L, F} <: ActorProxy
    mappingFn::F
end

actor_proxy!(::Type, proxy::MapProxy{L, F}, actor::A) where { L, A, F } = MapActor{L, A, F}(proxy.mappingFn, actor)

struct MapActor{L, A, F} <: Actor{L}
    mappingFn  :: F
    actor      :: A
end

on_next!(actor::MapActor{L},  data::L) where L = next!(actor.actor, actor.mappingFn(data))
on_error!(actor::MapActor, err)                = error!(actor.actor, err)
on_complete!(actor::MapActor)                  = complete!(actor.actor)

Base.show(io::IO, ::MapOperator{R}) where R   = print(io, "MapOperator( -> $R)")
Base.show(io::IO, ::MapProxy{L})    where L   = print(io, "MapProxy($L)")
Base.show(io::IO, ::MapActor{L})    where L   = print(io, "MapActor($L)")


text:
let
  # testing constructor and pdf
  d = BoolDistribution(0.3)
  @test pdf(d, true) == 0.3
  @test pdf(d, false) == 0.7

  # testing support
  @test support(d) == [true, false]

  # testing ==
  d2 = BoolDistribution(0.3)
  @test d == d2

  # testing hash
  @test hash(d) == hash(d.p)
end


text:
function get_Kσs(sp::LinearElasticityProblem{dim, TT}, dofs, cellvalues) where {dim, TT}
    E = YoungsModulus(sp)
    ν = PoissonRatio(sp)
    dh = sp.ch.dh
    n = ndofs_per_cell(dh)
    global_dofs = zeros(Int, n)
    Kσs = [zeros(TT, n, n) for i in 1:getncells(dh.grid)]
    Kσ_e = zeros(TT, n, n)
    ψ_e = zeros(TT, dim*n, dim*n)
    G = zeros(TT, dim*n, n)
    δ = eye(TT, dim)
    ϵ = zeros(TT, dim, dim)
    σ = zeros(TT, dim, dim)
    u = zeros(TT, dim, dim)
    n_basefuncs = getnbasefunctions(cellvalues)
    for (cellidx, cell) in enumerate(CellIterator(dh))
        Kσ_e .= 0
        reinit!(cellvalues, cell)
        celldofs!(global_dofs, dh, cellidx)
        for q_point in 1:getnquadpoints(cellvalues)
            dΩ = getdetJdV(cellvalues, q_point)
            for d in 1:dim
                ψ_e[(d-1)*dim+1:d*dim, (d-1)*dim+1:d*dim] .= 0
            end
            for a in 1:n_basefuncs
                ∇ϕ = shape_gradient(cellvalues, q_point, a)
                _u = @view dofs[(@view global_dofs[dim*(a-1) + (1:dim)])]
                @einsum u[i,j] = _u[i]*∇ϕ[j]
                @einsum ϵ[i,j] = 1/2*(u[i,j] + u[j,i])
                @einsum σ[i,j] = E*ν/(1-ν^2)*δ[i,j]*ϵ[k,k] + E*ν*(1+ν)*ϵ[i,j]
                for d in 1:dim
                    ψ_e[(d-1)*dim+1:d*dim, (d-1)*dim+1:d*dim] .+= σ
                    G[(dim*(d-1)+1):(dim*d), (a-1)*dim+d] .= ∇ϕ
                end
            end
            Kσ_e .+= G'*ψ_e*G*dΩ
        end
        Kσs[cellidx] .= Kσ_e
    end

    return Kσs
end

function buckling(problem::LinearElasticityProblem{dim, T}, ginfo, einfo) where {dim, T}
    dh = problem.ch.dh

    u = ginfo.K \ ginfo.f
    Kσs = get_Kσs(problem, u, einfo.cellvalues)
    Kσ = deepcopy(ginfo.K)

    if Kσ isa Symmetric
        Kσ.data.nzval .= 0
        assembler = JuAFEM.AssemblerSparsityPattern(Kσ.data, T[], Int[], Int[])
    else
        Kσ.nzval .= 0
        assembler = JuAFEM.AssemblerSparsityPattern(Kσ, T[], Int[], Int[])
    end

    global_dofs = zeros(Int, ndofs_per_cell(dh))
    Kσ_e = zeros(T, size(Kσs[1]))
    celliteratortype = CellIterator{typeof(dh).parameters...}
    _celliterator::celliteratortype = CellIterator(dh)
    TK = eltype(Kσs)
    for (i,cell) in enumerate(_celliterator)
        celldofs!(global_dofs, dh, i)
        if TK <: Symmetric
            JuAFEM.assemble!(assembler, global_dofs, Kσs[i].data)
        else
            JuAFEM.assemble!(assembler, global_dofs, Kσs[i])
        end
    end

    return ginfo.K, Kσ
end


text:
@safetestset "getpoints, getpointspec and getpointnum" begin
    using PoincareInvariants

    @test getpointspec(1, FirstFinDiffPlan) == 2
    @test getpointspec(50, FirstFinDiffPlan) == 50
    @test getpointspec(345, FirstFinDiffPlan) == 346

    @inferred Matrix{Float32} getpoints(x -> (x, 2x), Float32, 8, FirstFinDiffPlan)
    @inferred Vector{Float64} getpoints(x -> sin(x), Float64, 124, FirstFinDiffPlan)

    @test getpoints(identity, Float64, 4, FirstFinDiffPlan) ≈ [0, 0.25, 0.5, 0.75]
    @test getpoints(identity, Float64, 6, FirstFinDiffPlan) ≈ [n / 6 for n in 0:5]
    @test getpoints(cos, Float64, 124, FirstFinDiffPlan) ≈ [cos(n / 124) for n in 0:123]

    N = 6342
    testpnts = getpoints(identity, Float64, N, FirstFinDiffPlan)
    pnts3 = getpoints(Float32, N, FirstFinDiffPlan) do x
        x, x^2, sin(x)
    end
    @test pnts3 isa Matrix{Float32}
    @test size(pnts3) == (N, 3)

    @test pnts3[:, 1] ≈ testpnts
    @test pnts3[:, 2] ≈ testpnts .* testpnts
    @test pnts3[:, 3] ≈ sin.(testpnts)
end

@safetestset "compute! with FirstFinDiffPlan" begin
    using PoincareInvariants
    using StaticArrays: SVector
    using LinearAlgebra: dot

    T = Float64
    D = 3

    θ(z, t, p) = SVector{3}(p * z[2], -t * z[1], z[3])
    f(x) = ((s, c) = sincospi(2x); SVector{3}(2c,  5s, c + s))

    T = Float64
    D = 3
    N = 6
    fvals = getpoints(f, Float64, N, FirstFinDiffPlan)
    Δx = 1 / N

    testI = sum(1:6) do i
        dfi = if i == 1
            (fvals[2, :] .- fvals[6, :]) ./ (2 * Δx)
        elseif i == 6
            (fvals[1, :] .- fvals[5, :]) ./ (2 * Δx)
        else
            (fvals[i+1, :] .- fvals[i-1, :]) ./ (2 * Δx)
        end

        fi = fvals[i, :]
        return Δx * dot(θ(fi, 2, 3), dfi)
    end

    pinv = FirstPoincareInvariant{T, D}(θ, N, FirstFinDiffPlan{T, D}())
    @test compute!(pinv, fvals, 2, 3) ≈ testI atol=500eps()
end


text:
using LinearAlgebra
using SparseArrays
using MultivariatePolynomials
using DynamicPolynomials
using JuMP
using SumOfSquares
using CSDP
using Test

solver = with_optimizer(CSDP.Optimizer, printlevel=0)

    @polyvar x[1:8]

    # The matrix under consideration
    alpha = 3 + sqrt(3);
    beta = sqrt(3) - 1;
    a = sqrt(2/alpha);
    b = 1/sqrt(alpha);
    c = b;
    d = -sqrt(beta/alpha);
    f = (1 + im)*sqrt(1/(alpha*beta));
    U = [a 0; b b; c im*c; d f];
    V = [0 a; b -b; c -im*c; -im*f -d];
    M = U*V';

    for (gam, expected) in [(0.8723, :Infeasible), (0.8724, :Optimal)]
        @show gam
        @show expected
        Z = monomials(x, 1)

        function build_A(i)
            H = M[i,:]*M[i,:]' - (gam^2)*sparse([i],[i],[1],4,4)
            H = [real(H) -imag(H); imag(H) real(H)]
            dot(Z, H*Z)
        end
        A = build_A.(1:4)

        m = SOSModel(solver)

        # -- Q(x)'s -- : sums of squares
        # Monomial vector: [x1; ... x8]
        @variable m Q[1:4] SOSPoly(Z)

        # -- r's -- : constant sum of squares
        Z = monomials(x, 0)
        @variable m r[i=1:4,j=(i+1):4] >= 0

        # Constraint : -sum(Qi(x)*Ai(x)) - sum(rij*Ai(x)*Aj(x)) + I(x) >= 0
        expr = 0
        # Adding term
        for i = 1:4
            expr -= A[i]*Q[i]
        end
        for i = 1:4
            for j = (i+1):4
                expr -= A[i]*A[j]*r[i,j]
            end
        end
        # Constant term: I(x) = -(x1^4 + ... + x8^4)
        I = -sum(x.^4)
        expr = expr + I

        @constraint m expr >= 0

        println("solve")
        optimize!(m)
        println("solved")

        # Program is feasible, thus 0.8724 is an upper bound for mu.
        @test termination_status(m) == expected
    end


text:
norecurselist = ["+", "-", "*", "/"]

function add_norecurse(f::String)
    global norecurselist
    if findfirst(isequal(f), norecurselist) == nothing
        push!(norecurselist, f)
    end
end

add_norecurse(f::Method) = add_norecurse(string(f.name))

function add_norecurse(f::Function)
    if !isa(f, Core.Builtin)
        add_norecurse(first(methods(f)))
    end
end

function rm_norecurse(f::String)
    global norecurselist
    index = findfirst(isequal(f), norecurselist)
    if findfirst(isequal(f), norecurselist) == nothing
        deleteat!(norecurselist, index)
    end
end

rm_norecurse(f::Method) = rm_norecurse(string(f.name))

function rm_norecurse(f::Function)
    if !isa(f, Core.Builtin)
        rm_norecurse(first(methods(f)))
    end
end

function show_norecurse()
    global norecurselist
    for s in norecurselist
        println(s)
    end
end


text:
using RoME
using Test
using Statistics

## MWE Pose2Point2 from #388
@testset "basic Pose2Point2 test" begin

    fg = initfg()

    addVariable!(fg, :x1, Pose2)
    addVariable!(fg, :l1, Point2)

    addFactor!(fg, [:x1], PriorPose2(MvNormal([0.,0, 0], [0.01, 0.01, 0.01])))

    addFactor!(fg, [:x1; :l1], Pose2Point2(MvNormal([0.0,-1], [0.1,0.1])))

    ensureAllInitialized!(fg)

    tree = solveTree!(fg)

    M = getManifold(Pose2)
    @test isapprox(M, mean(M, getVal(fg, :x1)), ProductRepr([0,0], [1 0; 0 1]), atol=0.05)
    @test isapprox(mean(getVal(fg, :l1)), [0,-1], atol = 0.05)
end


text:
using FFTW

function FourierCoeffsOdd(f,ta,tb,M)
  N = 2*M+1
  h = (tb-ta)/N
  tt = range(ta,step=h,length=N)
  ff = f.(tt)
  hf = 1/N*rfft(ff)
  return (tt,ff,hf)
end

function FourierEvalOdd(ta,tb,hf,N)
    M = length(hf)-1
    if (N<2M+1)
        error("$M=length(hf)-1 should be less than (N+1)/2=$div(N+1,2)")
    end
    h = (tb-ta)/N
    TT = range(ta,step=h,length=N)
    MM = div(N,2)
    hF = zeros(Complex{Float64},MM+1)
    hF[1:(M+1)] = hf[1:(M+1)]
    PP = N*real.(irfft(hF,N))
    return (TT, PP)
end

function FourierCoeffsEven(f,ta,tb,M)
  N = 2*M
  h = (tb-ta)/N
  tt = range(ta,step=h,length=N)
  ff = f.(tt)
  hf = 1/N*rfft(ff)
  return (tt,ff,hf)
end

function FourierEvalEven(ta,tb,hf,N)
    M = length(hf)-1
    if (N<2M)
        error("$M=length(hf)-1 should be less than or equal to N/2=$div(N,2)")
    end
    h = (tb-ta)/N
    TT = range(ta,step=h,length=N)
    MM = div(N,2)
    hF = zeros(Complex{Float64},MM+1)
    hF[1:M] = hf[1:M]
    hF[M+1] = hf[M+1]/2
    PP = N*real.(irfft(hF,N))
    return (TT, PP)
end

function FourierCoeffs(tt,ff)
    N = length(ff)
    hf = 1/N*rfft(ff)
    return (tt,ff,hf)
end

function FourierEval(ta,tb,hf,N)
    if mod(N, 2) == 0
        return FourierEvalEven(ta,tb,hf,N)
    else
        return FourierEvalOdd(ta,tb,hf,N)
    end
end

function DiffMatrix(N)
    D = zeros(N,N)
    isodd(N) ? auxfcn = csc : auxfcn = cot
    for j in 1:N
        for k in 1:j-1
            D[j,k] = 0.5*(-1)^(j-k)*auxfcn((j-k)*pi/N)
            D[k,j] = -D[j,k]
        end
    end
    return D
end


text:
function make_count_matrix(counts)
    idx = findall(counts .> 0)
    countmat = Array{Int}(undef, length(idx), 2)
    countmat[:, 1] = idx
    countmat[:, 2] = counts[idx]
    return countmat
end

function format_counts_lda(countsdf::DataFrame)
    counts = Matrix{Int}[]
    for lib in names(countsdf)
        if lib == :term
            continue
        end

        #
        colcounts = convert(Array, countsdf[lib])
        countmat = make_count_matrix(colcounts)

        push!(counts, countmat)
    end
    return counts
end

function format_counts_ctm(countsdf::DataFrame)
    return format_counts_mmctm(countsdf)
end

function format_counts_mmctm(countsdfs::DataFrame...)
    counts = Vector{Matrix{Int}}[]

    if length(countsdfs) == 0
        return counts
    end

    for lib in names(countsdfs[1])
        if lib == :term
            continue
        end

        libcounts = Matrix{Int}[]
        for countsdf in countsdfs
            modality_counts = convert(Array, countsdf[lib])
            modality_countmat = make_count_matrix(modality_counts)
            push!(libcounts, modality_countmat)
        end

        push!(counts, libcounts)
    end

    return counts
end


text:
using MatrixNetworks: scomponents, bfs
using LinearAlgebra: det

export orientation
"""
    orientation(s::Simplex{N, K}) where {N, K}

Compute the relative orientation of the simplex with respect to the embedding space. Returns
`1` or `-1` if the simplex is non-degenerate and `0` otherwise. Requires `K == N+1`.
"""
function orientation(s::Simplex{N, K}) where {N, K}
    @assert K == N+1
    return sign(det(hcat(Wedge(s).vectors...)))
end

export change_orientation!
"""
    change_orientation!(c::Cell)

Reverse the orientation of a cell and return the cell.
"""
function change_orientation!(c::Cell)
    if length(c.points) >= 2
        a = c.points[1]
        c.points[1] = c.points[2]
        c.points[2] = a
    end
    for parent in keys(c.parents)
        c.parents[parent] = !c.parents[parent]
    end
    for child in c.children
        child.parents[c] = !child.parents[c]
    end
    return c
end

export orient_component!
"""
    orient_component!(cells::AbstractVector{Cell{N}}, adj::AbstractMatrix{<:Real},
        faces::Dict{Tuple{Cell{N}, Cell{N}}, Cell{N}}, root::Int) where N

Use breadth first search to orient the connected component of `cells` containing
`cells[root]` in the graph with adjacency matrix `adj`. `(adj, faces)` should be the output
of `adjacency(cells)`. If the cells are orientable, a consistent orientation will be
produced. Return `cells`.
"""
function orient_component!(cells::AbstractVector{Cell{N}},
    adj::AbstractMatrix{<:Real},
    faces::Dict{Tuple{Cell{N}, Cell{N}}, Cell{N}}, root::Int) where N
    orient!(cells[root])
    dists, _, predecessors = bfs(adj, root)
    for j in sortperm(dists)
        if dists[j] > 0 # exclude unreachables and the root
            c1, c2 = cells[predecessors[j]], cells[j]
            f = faces[(c1, c2)]
            if f.parents[c1] == f.parents[c2]
                change_orientation!(c2)
            end
        end
    end
    return cells
end

export orient!
"""
    orient!(cells::AbstractVector{Cell{N}}) where N

For each connected component of `cells`, use `orient_component!` to produce a consistent
orientation, if one exists. Return `cells`.
"""
function orient!(cells::AbstractVector{Cell{N}}) where N
    adj, faces = adjacency(cells)
    cc = scomponents(adj)
    roots = [findfirst(isequal(i), cc.map) for i in 1:cc.number]
    # start with higest dimensional component
    for i in reverse(sortperm([cells[r].K for r in roots]))
        orient_component!(cells, adj, faces, roots[i])
    end
    return cells
end

"""
    orient!(cell::Cell{N}) where N

If `cell` is of maximum dimension and is simplicial, orient it according to the orientation
of the embedding space. If `cell` is on the boundary, orient it accoring to its parent.
Otherwise leave its orientation as it is. Return `cell`.
"""
function orient!(cell::Cell{N}) where N
    num_parents = length(cell.parents)
    if (num_parents == 0) && (cell.K == length(cell.points) == N + 1)
        if orientation(Simplex(cell)) < 0
            change_orientation!(cell)
        end
    elseif (num_parents == 1) && (!(collect(values(cell.parents))[1]))
        change_orientation!(cell)
    end
    return cell
end

"""
    orient!(comp::CellComplex{N, K}) where {N, K}

If `K == N+1`, orient each highest dimensional cell according to the embedding space.
Otherwise orient the highest dimensional cells consistently with eachother, if such an
orientation exists. Return `comp`.
"""
function orient!(comp::CellComplex{N, K}) where {N, K}
    if K == N + 1
        map(orient!, comp.cells[end])
    else
        orient!(comp.cells[end])
    end
    return comp
end


text:
__precompile__(true)

module Geometry

using DocStringExtensions
using GeometryTypes: HomogenousMesh
using Rotations
using LinearAlgebra
using StaticArrays
using DataStructures
using NearestNeighbors
using ..Binary_BB_Trees
using ..MathKernel

include("mesh.jl")
include("mesh_create_rot_sym.jl")
include("mesh_create_swept.jl")
include("extensions.jl")
include("blob_types.jl")
include("top_down.jl")
include("util.jl")

export
    # mesh.jl
    Tri,
    Tet,
    eMesh,
    as_tet_eMesh,
    as_tri_eMesh,
    vertex_pos_for_tri_ind,
    vertex_pos_for_tet_ind,
    n_point,
    n_tri,
    n_tet,
    get_tri,
    get_tet,
    get_point,
    verify_mesh,
    transform!,
    invert!,
    mesh_inplace_rekey!,
    mesh_remove_unused_points!,
    delete_triangles!,
    sub_div_mesh,
    mesh_repair!,
    eMesh_half_plane,
    eMesh_sphere,
    eMesh_box,
    crop_mesh,
    extrude_mesh,
    eMesh_cylinder,

    # mesh_create_rot_sym.jl
    obj_from_point_sequence,

    # mesh_create_swept.jl
    create_swept_mesh,
    f_swept_circle,
    f_swept_triv,

    # blob_types.jl
    blob,
    eMesh_to_tree,

    # util.jl
    calc_min_max,
    calc_obb

    # extensions.jl

end


text:
export flex_tnep

"TNEP with flexible loads and storage, for transmission networks"
function flex_tnep(data::Dict{String,Any}, model_type::Type, optimizer; kwargs...)
    require_dim(data, :hour, :year)
    return _PM.run_model(
        data, model_type, optimizer, post_flex_tnep;
        ref_extensions = [ref_add_gen!, _PMACDC.add_ref_dcgrid!, _PMACDC.add_candidate_dcgrid!, ref_add_storage!, ref_add_ne_storage!, ref_add_flex_load!, _PM.ref_add_on_off_va_bounds!, _PM.ref_add_ne_branch!],
        multinetwork = true,
        kwargs...
    )
end

"TNEP with flexible loads and storage, for distribution networks"
function flex_tnep(data::Dict{String,Any}, model_type::Type{<:_PM.AbstractBFModel}, optimizer; kwargs...)
    require_dim(data, :hour, :year)
    return _PM.run_model(
        data, model_type, optimizer, post_flex_tnep;
        ref_extensions = [ref_add_gen!, ref_add_storage!, ref_add_ne_storage!, ref_add_flex_load!, _PM.ref_add_on_off_va_bounds!, ref_add_ne_branch_allbranches!, ref_add_frb_branch!, ref_add_oltc_branch!],
        solution_processors = [_PM.sol_data_model!],
        multinetwork = true,
        kwargs...
    )
end

"TNEP with flexible loads and storage, combines transmission and distribution networks"
function flex_tnep(t_data::Dict{String,Any}, d_data::Dict{String,Any}, t_model_type::Type, d_model_type::Type{<:_PM.AbstractBFModel}, optimizer; kwargs...)
    require_dim(t_data, :hour, :year)
    require_dim(d_data, :hour, :year)
    return run_model(
        t_data, d_data, t_model_type, d_model_type, optimizer, post_flex_tnep;
        t_ref_extensions = [ref_add_gen!, _PM.ref_add_on_off_va_bounds!, _PM.ref_add_ne_branch!, _PMACDC.add_ref_dcgrid!, _PMACDC.add_candidate_dcgrid!, ref_add_storage!, ref_add_ne_storage!, ref_add_flex_load!],
        d_ref_extensions = [ref_add_gen!, _PM.ref_add_on_off_va_bounds!, ref_add_ne_branch_allbranches!, ref_add_frb_branch!, ref_add_oltc_branch!, ref_add_storage!, ref_add_ne_storage!, ref_add_flex_load!],
        t_solution_processors = [_PM.sol_data_model!],
        d_solution_processors = [_PM.sol_data_model!, sol_td_coupling!],
        kwargs...
    )
end

# Here the problem is defined, which is then sent to the solver.
# It is basically a declaration of variables and constraints of the problem

"Builds transmission model."
function post_flex_tnep(pm::_PM.AbstractPowerModel; objective::Bool=true)
    # VARIABLES: defined within PowerModels(ACDC) can directly be used, other variables need to be defined in the according sections of the code
    for n in nw_ids(pm)

        # AC Bus
        _PM.variable_bus_voltage(pm; nw = n)

        # AC branch
        _PM.variable_branch_power(pm; nw = n)

        # DC bus
        _PMACDC.variable_dcgrid_voltage_magnitude(pm; nw = n)

        # DC branch
        _PMACDC.variable_active_dcbranch_flow(pm; nw = n)
        _PMACDC.variable_dcbranch_current(pm; nw = n)

        # AC-DC converter
        _PMACDC.variable_dc_converter(pm; nw = n)

        # Generator
        _PM.variable_gen_power(pm; nw = n)
        expression_gen_curtailment(pm; nw = n)

        # Storage
        _PM.variable_storage_power(pm; nw = n)
        variable_absorbed_energy(pm; nw = n)

        # Candidate AC branch
        variable_ne_branch_investment(pm; nw = n)
        variable_ne_branch_indicator(pm; nw = n, relax=true) # FlexPlan version: replaces _PM.variable_ne_branch_indicator().
        _PM.variable_ne_branch_power(pm; nw = n)
        _PM.variable_ne_branch_voltage(pm; nw = n)

        # Candidate DC bus
        _PMACDC.variable_dcgrid_voltage_magnitude_ne(pm; nw = n)

        # Candidate DC branch
        variable_ne_branchdc_investment(pm; nw = n)
        variable_ne_branchdc_indicator(pm; nw = n, relax=true) # FlexPlan version: replaces _PMACDC.variable_branch_ne().
        _PMACDC.variable_active_dcbranch_flow_ne(pm; nw = n)
        _PMACDC.variable_dcbranch_current_ne(pm; nw = n)

        # Candidate AC-DC converter
        variable_dc_converter_ne(pm; nw = n) # FlexPlan version: replaces _PMACDC.variable_dc_converter_ne().
        _PMACDC.variable_voltage_slack(pm; nw = n)

        # Candidate storage
        variable_storage_power_ne(pm; nw = n)
        variable_absorbed_energy_ne(pm; nw = n)

        # Flexible demand
        variable_flexible_demand(pm; nw = n)
        variable_energy_not_consumed(pm; nw = n)
        variable_total_demand_shifting_upwards(pm; nw = n)
        variable_total_demand_shifting_downwards(pm; nw = n)
    end

    # OBJECTIVE: see objective.jl
    if objective
        objective_min_cost_flex(pm)
    end

    # CONSTRAINTS: defined within PowerModels(ACDC) can directly be used, other constraints need to be defined in the according sections of the code
    for n in nw_ids(pm)
        _PM.constraint_model_voltage(pm; nw = n)
        _PM.constraint_ne_model_voltage(pm; nw = n)
        _PMACDC.constraint_voltage_dc(pm; nw = n)
        _PMACDC.constraint_voltage_dc_ne(pm; nw = n)
        for i in _PM.ids(pm, n, :ref_buses)
            _PM.constraint_theta_ref(pm, i, nw = n)
        end

        for i in _PM.ids(pm, n, :bus)
            constraint_power_balance_acne_dcne_flex(pm, i; nw = n)
        end
        if haskey(pm.setting, "allow_line_replacement") && pm.setting["allow_line_replacement"] == true
            for i in _PM.ids(pm, n, :branch)
                constraint_ohms_yt_from_repl(pm, i; nw = n)
                constraint_ohms_yt_to_repl(pm, i; nw = n)
                constraint_voltage_angle_difference_repl(pm, i; nw = n)
                constraint_thermal_limit_from_repl(pm, i; nw = n)
                constraint_thermal_limit_to_repl(pm, i; nw = n)
            end
        else
            for i in _PM.ids(pm, n, :branch)
                _PM.constraint_ohms_yt_from(pm, i; nw = n)
                _PM.constraint_ohms_yt_to(pm, i; nw = n)
                _PM.constraint_voltage_angle_difference(pm, i; nw = n)
                _PM.constraint_thermal_limit_from(pm, i; nw = n)
                _PM.constraint_thermal_limit_to(pm, i; nw = n)
            end
        end
        for i in _PM.ids(pm, n, :ne_branch)
            _PM.constraint_ne_ohms_yt_from(pm, i; nw = n)
            _PM.constraint_ne_ohms_yt_to(pm, i; nw = n)
            _PM.constraint_ne_voltage_angle_difference(pm, i; nw = n)
            _PM.constraint_ne_thermal_limit_from(pm, i; nw = n)
            _PM.constraint_ne_thermal_limit_to(pm, i; nw = n)
        end

        for i in _PM.ids(pm, n, :busdc)
            _PMACDC.constraint_power_balance_dc_dcne(pm, i; nw = n)
        end
        for i in _PM.ids(pm, n, :busdc_ne)
            _PMACDC.constraint_power_balance_dcne_dcne(pm, i; nw = n)
        end

        for i in _PM.ids(pm, n, :branchdc)
            _PMACDC.constraint_ohms_dc_branch(pm, i; nw = n)
        end
        for i in _PM.ids(pm, n, :branchdc_ne)
            _PMACDC.constraint_ohms_dc_branch_ne(pm, i; nw = n)
            _PMACDC.constraint_branch_limit_on_off(pm, i; nw = n)
        end

        for i in _PM.ids(pm, n, :convdc)
            _PMACDC.constraint_converter_losses(pm, i; nw = n)
            _PMACDC.constraint_converter_current(pm, i; nw = n)
            _PMACDC.constraint_conv_transformer(pm, i; nw = n)
            _PMACDC.constraint_conv_reactor(pm, i; nw = n)
            _PMACDC.constraint_conv_filter(pm, i; nw = n)
            if _PM.ref(pm,n,:convdc,i,"islcc") == 1
                _PMACDC.constraint_conv_firing_angle(pm, i; nw = n)
            end
        end
        for i in _PM.ids(pm, n, :convdc_ne)
            _PMACDC.constraint_converter_losses_ne(pm, i; nw = n)
            _PMACDC.constraint_converter_current_ne(pm, i; nw = n)
            _PMACDC.constraint_converter_limit_on_off(pm, i; nw = n)
            _PMACDC.constraint_conv_transformer_ne(pm, i; nw = n)
            _PMACDC.constraint_conv_reactor_ne(pm, i; nw = n)
            _PMACDC.constraint_conv_filter_ne(pm, i; nw = n)
            if _PM.ref(pm,n,:convdc_ne,i,"islcc") == 1
                _PMACDC.constraint_conv_firing_angle_ne(pm, i; nw = n)
            end
        end

        for i in _PM.ids(pm, n, :flex_load)
            constraint_total_flexible_demand(pm, i; nw = n)
            constraint_flex_bounds_ne(pm, i; nw = n)
        end
        for i in _PM.ids(pm, n, :fixed_load)
            constraint_total_fixed_demand(pm, i; nw = n)
        end

        for i in _PM.ids(pm, :storage, nw=n)
            constraint_storage_excl_slack(pm, i, nw = n)
            _PM.constraint_storage_thermal_limit(pm, i, nw = n)
            _PM.constraint_storage_losses(pm, i, nw = n)
        end
        for i in _PM.ids(pm, :ne_storage, nw=n)
            constraint_storage_excl_slack_ne(pm, i, nw = n)
            constraint_storage_thermal_limit_ne(pm, i, nw = n)
            constraint_storage_losses_ne(pm, i, nw = n)
            constraint_storage_bounds_ne(pm, i, nw = n)
        end

        if is_first_id(pm, n, :hour)
            for i in _PM.ids(pm, :storage, nw = n)
                constraint_storage_state(pm, i, nw = n)
            end
            for i in _PM.ids(pm, :storage_bounded_absorption, nw = n)
                constraint_maximum_absorption(pm, i, nw = n)
            end

            for i in _PM.ids(pm, :ne_storage, nw = n)
                constraint_storage_state_ne(pm, i, nw = n)
            end
            for i in _PM.ids(pm, :ne_storage_bounded_absorption, nw = n)
                constraint_maximum_absorption_ne(pm, i, nw = n)
            end

            for i in _PM.ids(pm, :flex_load, nw = n)
                constraint_red_state(pm, i, nw = n)
                constraint_shift_up_state(pm, i, nw = n)
                constraint_shift_down_state(pm, i, nw = n)
            end
        else
            if is_last_id(pm, n, :hour)
                for i in _PM.ids(pm, :storage, nw = n)
                    constraint_storage_state_final(pm, i, nw = n)
                end

                for i in _PM.ids(pm, :ne_storage, nw = n)
                    constraint_storage_state_final_ne(pm, i, nw = n)
                end

                for i in _PM.ids(pm, :flex_load, nw = n)
                    constraint_shift_state_final(pm, i, nw = n)
                end
            end

            # From second hour to last hour
            prev_n = prev_id(pm, n, :hour)
            first_n = first_id(pm, n, :hour)
            for i in _PM.ids(pm, :storage, nw = n)
                constraint_storage_state(pm, i, prev_n, n)
            end
            for i in _PM.ids(pm, :storage_bounded_absorption, nw = n)
                constraint_maximum_absorption(pm, i, prev_n, n)
            end
            for i in _PM.ids(pm, :ne_storage, nw = n)
                constraint_storage_state_ne(pm, i, prev_n, n)
            end
            for i in _PM.ids(pm, :ne_storage_bounded_absorption, nw = n)
                constraint_maximum_absorption_ne(pm, i, prev_n, n)
            end
            for i in _PM.ids(pm, :flex_load, nw = n)
                constraint_red_state(pm, i, prev_n, n)
                constraint_shift_up_state(pm, i, prev_n, n)
                constraint_shift_down_state(pm, i, prev_n, n)
                constraint_shift_duration(pm, i, first_n, n)
            end
        end

        # Constraints on investments
        if is_first_id(pm, n, :hour)
            # Inter-year constraints
            prev_nws = prev_ids(pm, n, :year)
            for i in _PM.ids(pm, :ne_branch; nw = n)
                constraint_ne_branch_activation(pm, i, prev_nws, n)
            end
            for i in _PM.ids(pm, :branchdc_ne; nw = n)
                constraint_ne_branchdc_activation(pm, i, prev_nws, n)
            end
            for i in _PM.ids(pm, :convdc_ne; nw = n)
                constraint_ne_converter_activation(pm, i, prev_nws, n)
            end
            for i in _PM.ids(pm, :ne_storage; nw = n)
                constraint_ne_storage_activation(pm, i, prev_nws, n)
            end
            for i in _PM.ids(pm, :flex_load; nw = n)
                constraint_flexible_demand_activation(pm, i, prev_nws, n)
            end
        end
    end
end

"Builds distribution model."
function post_flex_tnep(pm::_PM.AbstractBFModel; objective::Bool=true, intertemporal_constraints::Bool=true)

    for n in nw_ids(pm)

        # AC Bus
        _PM.variable_bus_voltage(pm; nw = n)

        # AC branch
        _PM.variable_branch_power(pm; nw = n)
        _PM.variable_branch_current(pm; nw = n)
        variable_oltc_branch_transform(pm; nw = n)

        # Generator
        _PM.variable_gen_power(pm; nw = n)
        expression_gen_curtailment(pm; nw = n)

        # Storage
        _PM.variable_storage_power(pm; nw = n)
        variable_absorbed_energy(pm; nw = n)

        # Candidate AC branch
        variable_ne_branch_investment(pm; nw = n)
        variable_ne_branch_indicator(pm; nw = n, relax=true) # FlexPlan version: replaces _PM.variable_ne_branch_indicator().
        _PM.variable_ne_branch_power(pm; nw = n, bounded = false) # Bounds computed here would be too limiting in the case of ne_branches added in parallel
        variable_ne_branch_current(pm; nw = n)
        variable_oltc_ne_branch_transform(pm; nw = n)

        # Candidate storage
        variable_storage_power_ne(pm; nw = n)
        variable_absorbed_energy_ne(pm; nw = n)

        # Flexible demand
        variable_flexible_demand(pm; nw = n)
        variable_energy_not_consumed(pm; nw = n)
        variable_total_demand_shifting_upwards(pm; nw = n)
        variable_total_demand_shifting_downwards(pm; nw = n)
    end

    if objective
        objective_min_cost_flex(pm)
    end

    for n in nw_ids(pm)
        _PM.constraint_model_current(pm; nw = n)
        constraint_ne_model_current(pm; nw = n)

        if haskey(dim_prop(pm), :sub_nw)
            constraint_td_coupling_power_reactive_bounds(pm; nw = n)
        end

        for i in _PM.ids(pm, n, :ref_buses)
            _PM.constraint_theta_ref(pm, i, nw = n)
        end

        for i in _PM.ids(pm, n, :bus)
            constraint_power_balance_acne_flex(pm, i; nw = n)
        end

        for i in _PM.ids(pm, n, :branch)
            constraint_dist_branch_tnep(pm, i; nw = n)
        end

        for i in _PM.ids(pm, n, :ne_branch)
            constraint_dist_ne_branch_tnep(pm, i; nw = n)
        end

        for i in _PM.ids(pm, n, :flex_load)
            constraint_total_flexible_demand(pm, i; nw = n)
            constraint_flex_bounds_ne(pm, i; nw = n)
        end
        for i in _PM.ids(pm, n, :fixed_load)
            constraint_total_fixed_demand(pm, i; nw = n)
        end

        for i in _PM.ids(pm, :storage, nw=n)
            constraint_storage_excl_slack(pm, i, nw = n)
            _PM.constraint_storage_thermal_limit(pm, i, nw = n)
            _PM.constraint_storage_losses(pm, i, nw = n)
        end
        for i in _PM.ids(pm, :ne_storage, nw=n)
            constraint_storage_excl_slack_ne(pm, i, nw = n)
            constraint_storage_thermal_limit_ne(pm, i, nw = n)
            constraint_storage_losses_ne(pm, i, nw = n)
            constraint_storage_bounds_ne(pm, i, nw = n)
        end

        if intertemporal_constraints
            if is_first_id(pm, n, :hour)
                for i in _PM.ids(pm, :storage, nw = n)
                    constraint_storage_state(pm, i, nw = n)
                end
                for i in _PM.ids(pm, :storage_bounded_absorption, nw = n)
                    constraint_maximum_absorption(pm, i, nw = n)
                end

                for i in _PM.ids(pm, :ne_storage, nw = n)
                    constraint_storage_state_ne(pm, i, nw = n)
                end
                for i in _PM.ids(pm, :ne_storage_bounded_absorption, nw = n)
                    constraint_maximum_absorption_ne(pm, i, nw = n)
                end

                for i in _PM.ids(pm, :flex_load, nw = n)
                    constraint_red_state(pm, i, nw = n)
                    constraint_shift_up_state(pm, i, nw = n)
                    constraint_shift_down_state(pm, i, nw = n)
                end
            else
                if is_last_id(pm, n, :hour)
                    for i in _PM.ids(pm, :storage, nw = n)
                        constraint_storage_state_final(pm, i, nw = n)
                    end

                    for i in _PM.ids(pm, :ne_storage, nw = n)
                        constraint_storage_state_final_ne(pm, i, nw = n)
                    end

                    for i in _PM.ids(pm, :flex_load, nw = n)
                        constraint_shift_state_final(pm, i, nw = n)
                    end
                end

                # From second hour to last hour
                prev_n = prev_id(pm, n, :hour)
                first_n = first_id(pm, n, :hour)
                for i in _PM.ids(pm, :storage, nw = n)
                    constraint_storage_state(pm, i, prev_n, n)
                end
                for i in _PM.ids(pm, :storage_bounded_absorption, nw = n)
                    constraint_maximum_absorption(pm, i, prev_n, n)
                end
                for i in _PM.ids(pm, :ne_storage, nw = n)
                    constraint_storage_state_ne(pm, i, prev_n, n)
                end
                for i in _PM.ids(pm, :ne_storage_bounded_absorption, nw = n)
                    constraint_maximum_absorption_ne(pm, i, prev_n, n)
                end
                for i in _PM.ids(pm, :flex_load, nw = n)
                    constraint_red_state(pm, i, prev_n, n)
                    constraint_shift_up_state(pm, i, prev_n, n)
                    constraint_shift_down_state(pm, i, prev_n, n)
                    constraint_shift_duration(pm, i, first_n, n)
                end
            end
        end

        # Constraints on investments
        if is_first_id(pm, n, :hour)
            for i in _PM.ids(pm, n, :branch)
                if !isempty(ne_branch_ids(pm, i; nw = n))
                    constraint_branch_complementarity(pm, i; nw = n)
                end
            end

            # Inter-year constraints
            prev_nws = prev_ids(pm, n, :year)
            for i in _PM.ids(pm, :ne_branch; nw = n)
                constraint_ne_branch_activation(pm, i, prev_nws, n)
            end
            for i in _PM.ids(pm, :ne_storage; nw = n)
                constraint_ne_storage_activation(pm, i, prev_nws, n)
            end
            for i in _PM.ids(pm, :flex_load; nw = n)
                constraint_flexible_demand_activation(pm, i, prev_nws, n)
            end
        end
    end
end

"Builds combined transmission and distribution model."
function post_flex_tnep(t_pm::_PM.AbstractPowerModel, d_pm::_PM.AbstractBFModel)

    # Transmission variables and constraints
    post_flex_tnep(t_pm; objective = false)

    # Distribution variables and constraints
    post_flex_tnep(d_pm; objective = false)

    # Variables related to the combined model
    # (No new variables are needed here.)

    # Constraints related to the combined model
    constraint_td_coupling(t_pm, d_pm)

    # Objective function of the combined model
    objective_min_cost_flex(t_pm, d_pm)

end


text:
struct SBitstreamHandler end

is_hardware_primitive(::Type{typeof(SBitstream)}, ::Type{<:Real}) = true
is_hardware_primitive(::Type{typeof(Base.broadcasted)},
                      ::Type{typeof(SBitstream)},
                      ::Type{<:Union{Real, AbstractArray{<:Real}}}) = true
gethandler(::Bool, ::Type{typeof(SBitstream)}, ::Type{<:Union{Real, AbstractArray{<:Real}}}) = SBitstreamHandler()
init_state(::SBitstreamHandler) = (id = 0,)

function (handler::SBitstreamHandler)(buffer, netlist, state, inputs, outputs)
    # set output as signed
    setsigned!(netlist, outputs[1], true)

    num_elements = join(netsize(outputs[1]), "*")

    # value
    if value(inputs[1]) isa AbstractArray
        val = lstrip.(name(inputs[1]), '-')
        isneg = Int.(value(inputs[1]) .< 0)
    else
        val = [lstrip(name(inputs[1]), '-')]
        isneg = [Int(value(inputs[1]) < 0)]
    end

    write(buffer, """
        // BEGIN bitstream_rng$(state.id)
        """)

    # create parameter if necessary
    if isparameter(inputs[1])
        input_string = name(inputs[1])
        input_size = bitwidth(inputs[1])
    else
        input_string = "bitstream_rng$(state.id)_value"
        input_size = bitwidth(inputs[1])

        write(buffer, "localparam $input_string = {")
        write(buffer, join(val, ", "))
        write(buffer, "};\n")
    end

    write(buffer, "localparam $(input_string)_isneg = $(length(isneg))'b")
    write(buffer, join(isneg, ""))
    write(buffer, ";\n")
    write_bcast_instantiation(buffer, "bitstream_rng$(state.id)", netsize(outputs[1]), """
        bitstream_rng #(
                .BITWIDTH($input_size),
                .VALUE($input_string[bitstream_rng$(state.id)_i*$input_size +: $input_size]),
                .IS_NEGATIVE($(input_string)_isneg[bitstream_rng$(state.id)_i])
            ) bitstream_rng$(state.id) (
                .CLK(CLK),
                .nRST(nRST),
                .out_p($(name(outputs[1]))_p[bitstream_rng$(state.id)_i]),
                .out_m($(name(outputs[1]))_m[bitstream_rng$(state.id)_i])
            );""")
    write(buffer, """
        // END bitstream_rng$(state.id)
        \n""")

    return buffer, (id = state.id + 1,)
end


text:
import OpenCL: cl
import OpenCL.cl: CLArray, to_host
using CLBLAS
import CLBLAS: gemm

const single_kernel = """
__kernel void single(
                     const __global float* x,
                     __global float* y)
 {
   const int i = get_global_id(0);
   y[i] = 1.0 / (1.0 + exp(-x[i]));
 }
"""

const exp_kernel = """
__kernel void exp_k(
                     const __global float* x,
                     __global float* y)
 {
   const int i = get_global_id(0);
   y[i] = exp(-x[i]);
 }
"""

const log_kernel = """
__kernel void log_k(
                     const __global float* x,
                     __global float* y)
 {
   const int i = get_global_id(0);
   y[i] = 1.0 + x[i];
 }
"""

const div_kernel = """
__kernel void div_k(
                     const __global float* x,
                     __global float* y)
 {
   const int i = get_global_id(0);
   y[i] = 1.0 / x[i];
 }
"""

function main2()
    x = rand(Float32, 50_000)
    y = similar(x)

    device, ctx, queue = cl.create_compute_context()
    x_buff = cl.Buffer(Float32, ctx, (:r, :copy), hostbuf=x)
    y_buff = cl.Buffer(Float32, ctx, :w, length(x))

    p = cl.Program(ctx, source=single_kernel) |> cl.build!
    k = cl.Kernel(p, "single")
    ev = queue(k, size(x), nothing, x_buff, y_buff)
    @time begin
        for i=1:1_000
            ev = queue(k, size(x), nothing, x_buff, y_buff; wait_on=[ev])
        end
        cl.wait(ev)
    end

    exp_p = cl.Program(ctx, source=exp_kernel) |> cl.build!
    exp_k = cl.Kernel(exp_p, "exp_k")
    log_p = cl.Program(ctx, source=log_kernel) |> cl.build!
    log_k = cl.Kernel(log_p, "log_k")
    div_p = cl.Program(ctx, source=div_kernel) |> cl.build!
    div_k = cl.Kernel(div_p, "div_k")
    @time begin
        for i=1:1_000
            ev = queue(exp_k, size(x), nothing, x_buff, y_buff; wait_on=[ev])
            ev = queue(log_k, size(x), nothing, y_buff, y_buff; wait_on=[ev])
            ev = queue(div_k, size(x), nothing, y_buff, y_buff; wait_on=[ev])
        end
        cl.wait(ev)
    end

end

function main()
    device, ctx, queue = cl.create_compute_context()
    p = cl.Program(ctx, source=template_kernel) |> cl.build!
    W = rand(Float32, 128*16, 64*16)
    x = rand(Float32, 64*16, 1)
    # b = rand(Float32, 128, 1)
    y = zeros(Float32, 128*16, 1)
    d_W = CLArray(queue, W)
    d_x = CLArray(queue, x)
    # d_b = CLArray(queue, b)
    d_y = CLArray(queue, y)

    @time begin
        ev = nothing
        for i=1:1000
            ev = gemm!('N', 'N', Float32(1.0), d_W, d_x, Float32(0.0), d_y)
        end
        cl.wait(ev)
    end

    @time begin
        for i=1:1000
            gemm!('N', 'N', Float32(1.0), W, x, Float32(0.0), y)
        end
    end
end