text:
"""
---
Apply Universal Filtered Orthogonal Frequency Division Multiplexing (UF-OFDM) to the time frequency matrix qamMat and returns a time domain ufofdm signal [1]

ufofdm is parametrized by its FFT size, the filter length (in samples) and the allocatedSubcarriers vector. Optional parameters are carrier size in subcarrier (by default RB size which is 12)
# --- Syntax
qamRx	= genereSignalufofdm(sigRx,nFFT,nCp,allocatedSubcarriers;sizeRB=12)
# ---  Input parameters
- sigRx	  : Time UF-OFDM recevied signal [Array{Complex{Float64},nbEch}]
- nFFT	  : FFT size [Int]
- L		  : Dolph Chebyshev filter length [Int]
- allocatedSubcarriers : Vector of index of allocated subcarriers [Array{Int,nbSubcarriers}]
- sizeRB  : Carrier size in subcarriers (default : LTE RB size: 12) [Int]
# ---  Output parameters
- qamRx  : Time frequency matrix : [Array{Complex{Float64},nbSubcarriers,nbSymb}]
- nbSymb			: Number of OFDM symbol tro be transmitted
- nbSubcarriers	: Number of allocated subcarriers (shall be < nFFT)
# ---
# v 1.0 - Robin Gerzaguet.
"""

function ufofdmSigDecode(sigRx,nFFT,L,allocatedSubcarriers;sizeRB=12,window=0)
	# ----------------------------------------------------
	# --- UF-OFDM parameters
	# ----------------------------------------------------
	# --- Symbol size (classic convolution)
	sizeSymb	  = nFFT + L - 1;
	nbSymbDecode		  = Int( floor( length(sigRx)/sizeSymb  ));
	# --- Number of physical RB (i.e carriers)
	nbRB		  = Int(floor(length(allocatedSubcarriers) / sizeRB));
	nbDataSubcarrier = length(allocatedSubcarriers);
	if window == 1
		# TODO Create windowing on UFOFDM Rx
	end
	qamRx	= zeros(Complex{Float64},nbDataSubcarrier,nbSymbDecode);
	# --- Iterative decoding
	for iB = 1 : 1 :nbSymbDecode
		# --- Extract subsignal
		subSig	  = sigRx[ 1+(iB-1)*sizeSymb:iB*sizeSymb ];
		# --- Apodisation
		if window != 0
			subSig	= subSig .* filterTimeDomain;
		end
		# --- Zero Padding to get a 2nFFT
		sigZP	  = [subSig;zeros(Complex{Float64},nFFT-L+1)];
		# --- Apply FFT
		sigFFT	  = fft(sigZP);
		# --- Getting only odd subcarrier
		sigFFT	  = sigFFT[1:2:end];
		# --- Getting received QPSK symbols
		qamRx[ :,iB ]  = sigFFT[allocatedSubcarriers];
	end
	return qamRx;
end
function ufofdmSigDecode(sigChan,ufofdm::StrucUFOFDM;window=0)
	return ufofdmSigDecode(sigChan,ufofdm.nFFT,ufofdm.L,ufofdm.allocatedSubcarriers,sizeRB=ufofdm.sizeRB,window=window);
end


text:
# baiaiyu
# 2021,Jul,31
# Repetition of https://stone-zeng.github.io/2019-10-03-exact-diagonalization/
# 统摄不懈，只为破敌，岂可擅退

module BH_ED

using SparseArrays
using LinearAlgebra
using KrylovKit
using Memoize
# using Distributed

struct BH_Model
    states::Vector{String}
    # states_n::Matrix{Int}
    N::Int
    L::Int
    Lp::Int
    μ::Float64
    U::Float64
    J::Float64
    V::Float64
    ε::Float64
end

include("./bh_hamiltonian.jl")

export give_a_BH_Model, give_H, BH_solver, give_the_num_of_particles_of_eigvector

@doc raw"""
    give_a_BH_Model(; N::Int        = 4
                    , Lp::Int       = 2
                    , L::Int        = N*Lp
                    , μ::Float64    = 1.0
                    , U::Float64    = 1.0
                    , J::Float64    = 1.0
                    , V::Float64    = 1.0
                    , ε::Float64    = 1.0
                    )

定义一个Bose–Hubbard model

```math
\hat{H}=-J \sum_{i=1}^{N}\left(\hat{a}_{i}^{\dagger} \hat{a}_{i+1}+\hat{a}_{i+1}^{\dagger} \hat{a}_{i}\right)+\sum_{i=1}^{N}-\mu \hat{n}_{i}+\frac{U}{2} \hat{n}_{i}\left(\hat{n}_{i}-1\right)+V \sum_{i=1}^{N} \hat{n}_{i} \hat{n}_{i+1}-\frac{\varepsilon}{2} \sum_{i=1}^{N}\left(\hat{a}_{i}^{2}+\hat{a}_{i+1}^{+2}\right)
```
# Arguments
- `N::Int`: 一维链长度
- `Lp::Int`: 每个格点最大容纳粒子数
- `L::Int`: 最多粒子数考虑 所有 0 到 L 个粒子的模型；最大可设置为 Lp*N，

"""
function give_a_BH_Model(;N::Int        = 4
                        , Lp::Int       = 2
                        , L::Int        = N*Lp
                        , μ::Float64    = 1.0
                        , U::Float64    = 1.0
                        , J::Float64    = 1.0
                        , V::Float64    = 1.0
                        , ε::Float64    = 1.0
                        )
    states = String[]
    for i in 0:L
        append!( states, give_all_states(N, i, Lp) )
    end
    # states = Array{Int}( undef, 2, length(states) )
    # states_n[1, :] .= [i for i in 1:length(states)]
    # states_n[2, :] .= map(states) do x
    #     parse(Int, x, base = Lp + 1)
    # end
    BH_Model(states, N, L, Lp, μ, U, J, V, ε)
    # BH_Model(states, states_n, N, L, Lp, μ, J, V, ε)

end

```
给出一个BH模型的哈密顿矩阵
```
function give_H(the_model::BH_Model)
    nhil = length(the_model.states)
    H = spzeros(Float64, nhil, nhil)

    # if nhil > 10000
    #     H = spzeros(Float64, nhil, nhil)
    # else
    #     H = zeros(Float64, nhil, nhil)
    # end

    if the_model.μ != 0.0
        fill_H_μ_U!(H, the_model)
    end

    if the_model.J != 0.0
        fill_H_J!(H, the_model)
    end

    if the_model.V != 0.0
        fill_H_V!(H, the_model)
    end

    if the_model.ε != 0.0
        fill_H_ε!(H, the_model)
    end

    H
end

# ```
# 给出一个BH模型的哈密顿矩阵, 并行计算，最大可利用 4 核心
# ```
# function give_H_p(the_model::BH_Model)
#     nhil = length(the_model.states)
#     H = spzeros(Float64, nhil, nhil)
#     h1 = spzeros(Float64, nhil, nhil)
#     h2 = spzeros(Float64, nhil, nhil)
#     h3 = spzeros(Float64, nhil, nhil)
#     h4 = spzeros(Float64, nhil, nhil)
#     call_result =[]
#     # if nhil > 10000
#     #     H = spzeros(Float64, nhil, nhil)
#     # else
#     #     H = zeros(Float64, nhil, nhil)
#     # end

#     if the_model.μ != 0.0
#         # @spawn push!( call_state, fill_H_μ_U!(h1, the_model) )
#         push!( call_result, remotecall(BH_ED.fill_H_μ_U!, workers()[1], h1, the_model) )
#     end

#     if the_model.J != 0.0
#         # @spawn call_state[2] = fill_H_J!(h2, the_model)
#         push!( call_result, remotecall(BH_ED.fill_H_J!, workers()[2], h2, the_model) )
#     end

#     if the_model.V != 0.0
#         # @spawn call_state[3] = fill_H_V!(h3, the_model)
#         push!( call_result, remotecall(BH_ED.fill_H_V!, workers()[3], h3, the_model) )
#     end

#     if the_model.ε != 0.0
#         # @spawn call_state[4] = fill_H_ε!(h4, the_model)
#         push!( call_result, remotecall(BH_ED.fill_H_ε!, workers()[4], h4, the_model) )
#     end

#     for r in result
#         H = H + fetch(r)
#     end

#     H
# end

@doc raw"""
    BH_solver(; N::Int        = 4
              , Lp::Int       = 2
              , L::Int        = N*Lp
              , μ::Float64    = 1.0
              , U::Float64    = 1.0
              , J::Float64    = 1.0
              , V::Float64    = 1.0
              , ε::Float64    = 1.0
              , screen_p::Bool= true
              , bool_cf::Bool = true
                    )

定义并求解 Bose–Hubbard model, 返回最小特征值与对应的特征向量

```math
\hat{H}=-J \sum_{i=1}^{N}\left(\hat{a}_{i}^{\dagger} \hat{a}_{i+1}+\hat{a}_{i+1}^{\dagger} \hat{a}_{i}\right)+\sum_{i=1}^{N}-\mu \hat{n}_{i}+\frac{U}{2} \hat{n}_{i}\left(\hat{n}_{i}-1\right)+V \sum_{i=1}^{N} \hat{n}_{i} \hat{n}_{i+1}-\frac{\varepsilon}{2} \sum_{i=1}^{N}\left(\hat{a}_{i}^{2}+\hat{a}_{i+1}^{+2}\right)
```
# Arguments
- `N::Int`: 一维链长度
- `Lp::Int`: 每个格点最大容纳粒子数
- `L::Int`: 最多粒子数考虑 所有 0 到 L 个粒子的模型；最大可设置为 Lp*N
- `screen_p::Bool`: 是否屏幕打印信息
- `bool_cf::Bool`: 是否计算关联函数

"""
function BH_solver(;N::Int      = 4
                , Lp::Int       = 2
                , L::Int        = N*Lp
                , μ::Float64    = 1.0
                , U::Float64    = 1.0
                , J::Float64    = 1.0
                , V::Float64    = 1.0
                , ε::Float64    = 1.0
                , screen_p::Bool= true
                , bool_cf::Bool = true
    )::Tuple{ Vector{Float64}, Vector{Vector{Float64}}, Float64, Vector{Float64} }

    the_bhm = give_a_BH_Model(N = N, Lp = Lp, L = L, μ = μ, U = U, J = J, V = V, ε = ε)
    screen_p && println("- N: $(the_bhm.N) with Lp: $(the_bhm.Lp)")
    screen_p && println("- We got all required states.")
    screen_p && println("- Dimension of Hilbert space: $(length(the_bhm.states))")
    H = give_H(the_bhm)
    screen_p && println("- We get Hamiltonion Matrix.")
    screen_p && print("- Number of nonzero elements: ")
    screen_p && print(nnz(H))
    screen_p && print(" (",round((100.0*nnz(H))/(H.m*H.n), digits = 3))
    screen_p && println("%) of total.")
    if H == H'
        screen_p && println("- Hamiltonion Matrix is a Real Symmetric Matrix.")
    else
        throw("Hamiltonion Matrix is NOT a Real Symmetric Matrix!")
        # return H
    end
    eigsol = eigsolve(H, 1, :SR, eager = true)
    n_ev = give_the_num_of_particles_of_eigvector(eigsol[2][1], the_bhm)

    cf_r = fill(NaN64, the_bhm.N + 1)
    bool_cf && begin
        screen_p && println("- Strat to calculate the correlation function.")
        cf_r .= c_ff(the_bhm, eigsol[2][1])
    end

    screen_p && println("- Ground energy is $(eigsol[1][1])")
    screen_p && println("- The corresponding number of particles is $n_ev")
    screen_p && println("- Then We will return the eigenvalue, the corresponding eigenvector, the number of particles and the correlation function(as your setting).")

    eigsol[1], eigsol[2], n_ev, cf_r
    # H
end

end

ext: jl


text:
######################################################################
# Event System UTs
# -----
# Licensed under MIT License.
module TestEvents
using Test
using KirUtil
using KirUtil.Events

@testset "Events" begin
    @testset "meta" begin
        foo_success = false
        @event function foo{:test_event1}(arg1, arg2; kwarg = false)
            foo_success = arg1 == 24 && arg2 == 25 && kwarg == 42
        end

        @test emit(:test_event1, 24, 25, kwarg=42) === nothing
        @test foo_success
    end

    @testset "callback" begin
        success = false
        listen(:test_event2) do arg1, arg2
            @assert arg1 == 42
            @assert arg2 == 420
            success = true
        end
        @test emit(:test_event2, 42, 420) === nothing
        @test success

        count = false
        listen_once(:test_event3) do arg1, arg2
            @assert arg1 == 69
            @assert arg2 == .69
            count = true
        end
        @test emit(:test_event3, 69, .69) === nothing
        @test emit(:test_event3, 42, 420) === nothing
        @test count == 1
    end
end

end # module TestEvents


text:
const weber_portaudio = "/Users/davidlittle/.julia/v0.6/Homebrew/deps/usr/lib/libportaudio.2.dylib"
const weber_sound = "/Users/davidlittle/MEGA/Research/software/TimedSound/deps/usr/lib/libweber-sound.4.dylib"


text:
@testset verbose = true "Tests for Predefined Conversions" begin

    using JavaCall.Conversions
    using JavaCall.JNI

    @testset "Primitive Conversions" begin
        @testset "Boolean conversions" begin
            @test convert_to_julia(Bool, JNI_TRUE)
            @test_false convert_to_julia(Bool, JNI_FALSE)

            @test convert_to_jni(jboolean, true) == JNI_TRUE
            @test convert_to_jni(jboolean, false) == JNI_FALSE
        end

        @testset "Byte conversions" begin
            @test convert_to_julia(Int8, jbyte(1)) == Int8(1)
            @test convert_to_jni(jbyte, Int8(1)) == jbyte(1)
        end

        @testset "Char conversions" begin
            @test convert_to_julia(Char, jchar('a')) == 'a'
            @test Int(convert_to_julia(Char, jchar('a'))) == 97
            @test convert_to_jni(jchar, 'a') == jchar('a')
            @test Int(convert_to_jni(jchar, 'a')) == 97
        end

        @testset "Short conversions" begin
            @test convert_to_julia(Int16, jshort(1)) == Int16(1)
            @test convert_to_jni(jshort, Int16(1)) == jshort(1)
        end

        @testset "Integer conversions" begin
            @test convert_to_julia(Int32, jint(1)) == Int32(1)
            @test convert_to_jni(jint, Int32(1)) == jint(1)
        end

        @testset "Long conversions" begin
            @test convert_to_julia(Int64, jlong(1)) == 1
            @test convert_to_jni(jlong, 1) == jlong(1)
        end

        @testset "Float conversions" begin
            @test convert_to_julia(Float32, jfloat(1.1)) == Float32(1.1)
            @test convert_to_jni(jfloat, Float32(1.1)) == jfloat(1.1)
        end

        @testset "Double conversions" begin
            @test convert_to_julia(Float64, jdouble(1.1)) == 1.1
            @test convert_to_jni(jdouble, 1.1) == jdouble(1.1)
        end

        @testset "Void conversions" begin
            @test convert_to_julia(Nothing, nothing) === nothing
            @test convert_to_jni(jvoid, nothing) === nothing
        end
    end

    @testset "Primitive Arrays Conversions" begin

        iscopy = Ref(JNI_TRUE)
        iscopyptr = convert(Ptr{jboolean}, pointer_from_objref(iscopy))

        @testset "Boolean Array" begin
            jni_array = JNI.new_boolean_array(5)
            JNI.set_boolean_array_region(jni_array, 0, 5, jboolean[
                JNI_TRUE, JNI_FALSE, JNI_FALSE, JNI_TRUE, JNI_TRUE
            ])
            julia_array = convert_to_julia(Vector{Bool}, jni_array)
            @test julia_array == [true, false, false, true, true]

            julia_array = [true, false, false, true, true]
            jni_array = convert_to_jni(jbooleanArray, julia_array)

            jni_ptr = JNI.get_boolean_array_elements(jni_array, iscopyptr)
            @test map(x -> convert(jboolean, x), unsafe_wrap(Vector{jboolean}, jni_ptr, (5,))) ==
                jboolean[JNI_TRUE, JNI_FALSE, JNI_FALSE, JNI_TRUE, JNI_TRUE]
        end

        @testset "Char Array" begin
            jni_array = JNI.new_char_array(5)
            JNI.set_char_array_region(jni_array, 0, 5, jchar[
                jchar('h'), jchar('e'), jchar('l'), jchar('l'), jchar('o')
            ])
            julia_array = convert_to_julia(Vector{Char}, jni_array)
            @test julia_array == ['h', 'e', 'l', 'l', 'o']

            julia_array = ['h', 'e', 'l', 'l', 'o']
            jni_array = convert_to_jni(jcharArray, julia_array)

            jni_ptr = JNI.get_char_array_elements(jni_array, iscopyptr)
            @test map(x -> convert(jchar, x), unsafe_wrap(Vector{jchar}, jni_ptr, (5,))) ==
                jchar[jchar('h'), jchar('e'), jchar('l'), jchar('l'), jchar('o')]
        end

        @testset "Short Array" begin
            jni_array = JNI.new_short_array(5)
            JNI.set_short_array_region(jni_array, 0, 5, jshort[
                jshort(1), jshort(2), jshort(3), jshort(4), jshort(5)
            ])
            julia_array = convert_to_julia(Vector{Int16}, jni_array)
            @test julia_array == [Int16(1), Int16(2), Int16(3), Int16(4), Int16(5)]

            julia_array = [Int16(1), Int16(2), Int16(3), Int16(4), Int16(5)]
            jni_array = convert_to_jni(jshortArray, julia_array)

            jni_ptr = JNI.get_short_array_elements(jni_array, iscopyptr)
            @test map(x -> convert(jshort, x), unsafe_wrap(Vector{jshort}, jni_ptr, (5,))) ==
                jshort[jshort(1), jshort(2), jshort(3), jshort(4), jshort(5)]
        end

        @testset "Int Array" begin
            jni_array = JNI.new_int_array(5)
            JNI.set_int_array_region(jni_array, 0, 5, jint[
                jint(1), jint(2), jint(3), jint(4), jint(5)
            ])
            julia_array = convert_to_julia(Vector{Int32}, jni_array)
            @test julia_array == [Int32(1), Int32(2), Int32(3), Int32(4), Int32(5)]

            julia_array = [Int32(1), Int32(2), Int32(3), Int32(4), Int32(5)]
            jni_array = convert_to_jni(jintArray, julia_array)

            jni_ptr = JNI.get_int_array_elements(jni_array, iscopyptr)
            @test map(x -> convert(jint, x), unsafe_wrap(Vector{jint}, jni_ptr, (5,))) ==
                jint[jint(1), jint(2), jint(3), jint(4), jint(5)]
        end

        @testset "Long Array" begin
            jni_array = JNI.new_long_array(5)
            JNI.set_long_array_region(jni_array, 0, 5, jlong[
                jlong(1), jlong(2), jlong(3), jlong(4), jlong(5)
            ])
            julia_array = convert_to_julia(Vector{Int64}, jni_array)
            @test julia_array == [1, 2, 3, 4, 5]

            julia_array = [1, 2, 3, 4, 5]
            jni_array = convert_to_jni(jlongArray, julia_array)

            jni_ptr = JNI.get_long_array_elements(jni_array, iscopyptr)
            @test map(x -> convert(jlong, x), unsafe_wrap(Vector{jlong}, jni_ptr, (5,))) ==
                jlong[jlong(1), jlong(2), jlong(3), jlong(4), jlong(5)]
        end

        @testset "Float Array" begin
            jni_array = JNI.new_float_array(5)
            JNI.set_float_array_region(jni_array, 0, 5, jfloat[
                jfloat(1.1), jfloat(2.2), jfloat(3.3), jfloat(4.4), jfloat(5.5)
            ])
            julia_array = convert_to_julia(Vector{Float32}, jni_array)
            @test julia_array == [
                Float32(1.1), Float32(2.2), Float32(3.3), Float32(4.4), Float32(5.5)
            ]

            julia_array = [
                Float32(1.1), Float32(2.2), Float32(3.3), Float32(4.4), Float32(5.5)
            ]
            jni_array = convert_to_jni(jfloatArray, julia_array)

            jni_ptr = JNI.get_float_array_elements(jni_array, iscopyptr)
            @test map(x -> convert(jfloat, x), unsafe_wrap(Vector{jfloat}, jni_ptr, (5,))) ==
                jfloat[jfloat(1.1), jfloat(2.2), jfloat(3.3), jfloat(4.4), jfloat(5.5)]
        end

        @testset "Double Array" begin
            jni_array = JNI.new_double_array(5)
            JNI.set_double_array_region(jni_array, 0, 5, jdouble[
                jdouble(1.1), jdouble(2.2), jdouble(3.3), jdouble(4.4), jdouble(5.5)
            ])
            julia_array = convert_to_julia(Vector{Float64}, jni_array)
            @test julia_array == [1.1, 2.2, 3.3, 4.4, 5.5]

            julia_array = [1.1, 2.2, 3.3, 4.4, 5.5]
            jni_array = convert_to_jni(jdoubleArray, julia_array)

            jni_ptr = JNI.get_double_array_elements(jni_array, iscopyptr)
            @test map(x -> convert(jdouble, x), unsafe_wrap(Vector{jdouble}, jni_ptr, (5,))) ==
                jdouble[jdouble(1.1), jdouble(2.2), jdouble(3.3), jdouble(4.4), jdouble(5.5)]
        end
    end
end


text:
using MultiTankMaterialBalance
using FiniteDiff
using Test

const Float = Float64

# Function to calculate the value and gradient of the objective function
function fun(fset::FittingSet, solver::NewtonSolver, targ::TargetFunction, adjoint::AdjointSolver)
    return (x, grad) -> begin
        setparams!(fset, x)
        solve!(solver)
        update_targ!(targ)
        if length(grad) > 0
            solve!(adjoint)
            copyto!(grad, adjoint.g)
        end
        return getvalue(targ)
    end
end

@testset "MultiTankMaterialBalance.jl" begin
    # Path to csv files with source data
    workdir = "./data"

    # Initial data for calculation
    opts_csv = Dict("dateformat" => "dd.mm.yyyy", "delim" => ";")
    df_rates = read_rates(joinpath(workdir, "tank_prod.csv"), opts_csv)
    df_params = read_params(joinpath(workdir, "tank_params.csv"), opts_csv)

    # Additional parameter preprocessing
    process_params!(df_params, df_rates)

    # Description of the forward problem
    prob = NonlinearProblem{Float}(df_rates, df_params)

    # Linear equation solver
    linalg = DenseLinearSolver{Float}(prob)

    # Algorithm to solve forward problem
    opts_nsol = Dict("maxiters" => 10, "P_tol" => 1e-7, "r_tol" => 1e-5)
    solver = NewtonSolver{Float}(prob, linalg, opts_nsol)

    # Solve forward problem
    solve!(solver)
    @test solver.success

    # Parameter scaling method
    scale = SigmoidScaling{Float}(df_params)

    # Fitting parameter list
    fset = FittingSet{Float}(df_params, prob, scale)

    # Objective function
    opts_targ = Dict("alpha_resp" => 1, "alpha_bhp" => 0.01, "alpha_inj" => 0.01, "alpha_lb" => 10, "alpha_ub" => 10, "alpha_l2" => 1)
    targ = TargetFunction{Float}(df_rates, df_params, prob, fset, opts_targ)

    # Algorithm to calculate objective function gradient
    adjoint_ = AdjointSolver{Float}(prob, targ, linalg, fset)

    # Function to calculate the value and gradient of the objective function
    optim_fun = fun(fset, solver, targ, adjoint_)

    # Initial parameter values
    initial_x = copy(getparams!(fset))
    # Preallocate gradient buffer
    grad = similar(initial_x)

    # Computing the gradient by the finite difference method
    grad_fd = FiniteDiff.finite_difference_gradient(x -> optim_fun(x, Float[]), initial_x, Val(:central); relstep=1e-5)

    # Computing the gradient by the adjoint equation method
    _ = optim_fun(initial_x, grad)

    # Comparing gradients computed in different ways
    @test all(isapprox.(grad, grad_fd; rtol=0.0001))
end


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

se33 = speye(3)
do33 = ones(3)
@test isequal(se33 \ do33, do33)

# based on deps/Suitesparse-4.0.2/UMFPACK/Demo/umfpack_di_demo.c

using Base.SparseArrays.UMFPACK.increment!

A0 = sparse(increment!([0,4,1,1,2,2,0,1,2,3,4,4]),
           increment!([0,4,0,2,1,2,1,4,3,2,1,2]),
           [2.,1.,3.,4.,-1.,-3.,3.,6.,2.,1.,4.,2.], 5, 5)

for Tv in (Float64, Complex128)
    for Ti in Base.uniontypes(Base.SparseArrays.UMFPACK.UMFITypes)
        A = convert(SparseMatrixCSC{Tv,Ti}, A0)
        lua = lufact(A)
        @test nnz(lua) == 18
        @test_throws KeyError lua[:Z]
        L,U,p,q,Rs = lua[:(:)]
        @test (Diagonal(Rs) * A)[p,q] ≈ L * U

        det(lua) ≈ det(Array(A))

        b = [8., 45., -3., 3., 19.]
        x = lua\b
        @test x ≈ float([1:5;])

        @test norm(A*x-b,1) < eps(1e4)
        z = complex.(b,zeros(b))
        x = Base.SparseArrays.A_ldiv_B!(lua, z)
        @test x ≈ float([1:5;])
        @test z === x
        y = similar(z)
        A_ldiv_B!(y, lua, complex.(b,zeros(b)))
        @test y ≈ x

        @test norm(A*x-b,1) < eps(1e4)

        b = [8., 20., 13., 6., 17.]
        x = lua'\b
        @test x ≈ float([1:5;])

        @test norm(A'*x-b,1) < eps(1e4)
        z = complex.(b,zeros(b))
        x = Base.SparseArrays.Ac_ldiv_B!(lua, z)
        @test x ≈ float([1:5;])
        @test x === z
        y = similar(x)
        Base.SparseArrays.Ac_ldiv_B!(y, lua, complex.(b,zeros(b)))
        @test y ≈ x

        @test norm(A'*x-b,1) < eps(1e4)
        x = lua.'\b
        @test x ≈ float([1:5;])

        @test norm(A.'*x-b,1) < eps(1e4)
        x = Base.SparseArrays.At_ldiv_B!(lua,complex.(b,zeros(b)))
        @test x ≈ float([1:5;])
        y = similar(x)
        Base.SparseArrays.At_ldiv_B!(y, lua,complex.(b,zeros(b)))
        @test y ≈ x

        @test norm(A.'*x-b,1) < eps(1e4)

        # Element promotion and type inference
        @inferred lua\ones(Int, size(A, 2))
    end
end

Ac0 = complex.(A0,A0)
for Ti in Base.uniontypes(Base.SparseArrays.UMFPACK.UMFITypes)
    Ac = convert(SparseMatrixCSC{Complex128,Ti}, Ac0)
    lua = lufact(Ac)
    L,U,p,q,Rs = lua[:(:)]
    @test (Diagonal(Rs) * Ac)[p,q] ≈ L * U
end

for elty in (Float64, Complex128)
    for (m, n) in ((10,5), (5, 10))
        A = sparse([1:min(m,n); rand(1:m, 10)], [1:min(m,n); rand(1:n, 10)], elty == Float64 ? randn(min(m, n) + 10) : complex.(randn(min(m, n) + 10), randn(min(m, n) + 10)))
        F = lufact(A)
        L, U, p, q, Rs = F[:(:)]
        @test (Diagonal(Rs) * A)[p,q] ≈ L * U
    end
end

#4523 - complex sparse \
x = speye(2) + im * speye(2)
@test (x*(lufact(x) \ ones(2))) ≈ ones(2)

@test det(sparse([1,3,3,1], [1,1,3,3], [1,1,1,1])) == 0

# UMFPACK_ERROR_n_nonpositive
@test_throws ArgumentError lufact(sparse(Int[], Int[], Float64[], 5, 0))

#15099
for (Tin, Tout) in (
        (Complex32, Complex128),
        (Complex64, Complex128),
        (Complex128, Complex128),
        (Float16, Float64),
        (Float32, Float64),
        (Float64, Float64),
        (Int, Float64),
    )

    F = lufact(sparse(ones(Tin, 1, 1)))
    L = sparse(ones(Tout, 1, 1))
    @test F[:p] == F[:q] == [1]
    @test F[:Rs] == [1.0]
    @test F[:L] == F[:U] == L
    @test F[:(:)] == (L, L, [1], [1], [1.0])
end

for T in (BigFloat, Complex{BigFloat})
    @test_throws ArgumentError lufact(sparse(ones(T, 1, 1)))
end

#size(::UmfpackLU)
let
    m = n = 1
    F = lufact(sparse(ones(m, n)))
    @test size(F) == (m, n)
    @test size(F, 1) == m
    @test size(F, 2) == n
    @test size(F, 3) == 1
    @test_throws ArgumentError size(F,-1)
end

let
    a = rand(5)
    @test_throws ArgumentError Base.SparseArrays.UMFPACK.solve!(a, lufact(speye(5,5)), a, Base.SparseArrays.UMFPACK.UMFPACK_A)
    aa = complex(a)
    @test_throws ArgumentError Base.SparseArrays.UMFPACK.solve!(aa, lufact(complex(speye(5,5))), aa, Base.SparseArrays.UMFPACK.UMFPACK_A)
end

#18246,18244-lufact sparse pivot
let A = speye(4)
    A[1:2,1:2] = [-.01 -200; 200 .001]
    F = lufact(A)
    @test F[:p] == [3 ; 4 ; 2 ; 1]
end

# Test that A[c|t]_ldiv_B!{T<:Complex}(X::StridedMatrix{T}, lu::UmfpackLU{Float64},
# B::StridedMatrix{T}) works as expected.
let N = 10, p = 0.5
    A = N*speye(N) + sprand(N, N, p)
    X = zeros(Complex{Float64}, N, N)
    B = complex.(rand(N, N), rand(N, N))
    luA, lufA = lufact(A), lufact(Array(A))
    @test A_ldiv_B!(copy(X), luA, B) ≈ A_ldiv_B!(copy(X), lufA, B)
    @test At_ldiv_B!(copy(X), luA, B) ≈ At_ldiv_B!(copy(X), lufA, B)
    @test Ac_ldiv_B!(copy(X), luA, B) ≈ Ac_ldiv_B!(copy(X), lufA, B)
end


text:
using GNSSTools
using Statistics
using PyPlot
pygui(true)
include("ch4_data_plot_functions.jl")

directory = "/home/sjbilardi/Dropbox/Apps/Overleaf/MastersThesis/"
folder = "figures/l1ca_skydel/"

file = "/home/sjbilardi/Documents/SeNSeLab/GNSSTools.jl/test/data/skydel/20200301_000000_L1_1582.0_60.0.sc16"
data = loaddata(:sc16, file, 60e6, 1582e6, L1_freq, 2, skip_to=2)

replica = definereplica(define_l1ca_code_type(), 60e6, 2e-3)
detectable_prns = []
for i in 1:32
    fd_est, n0_est, SNR_est,
    P_d, above_threshold = courseacquisition(data, replica, i; M=20,
                                             Δfd=1/(1*replica.t_length))
    if above_threshold && (P_d >= 0.9)
        detectable = true
        push!(detectable_prns, i)
    else
        detectable = false
    end
    f_code_d, f_code_dd = GNSSTools.calc_doppler_code_rate(l1ca_chipping_rate,
                                                           replica.signal_type.sig_freq,
                                                           fd_est, 0)
    n0 = calcinitcodephase(l1ca_code_length, f_code_d, f_code_dd, data.f_s,
                           n0_est)
    println("PRN $i\t$SNR_est\t$fd_est\t$n0_est\t$n0\t$(Float64(round(P_d, digits=2)))\t$detectable")
end
println("$(length(detectable_prns)) Detected PRNs: $detectable_prns")

fig = figure(figsize=(7.5, 5))
acquisition_T = 2e-3
fine_acq_T = 10e-3
fine_acq_sub_T = 1e-3
tracking_T = 1e-3
M = 20
linewidth = 0.75
tracking_σᵩ = Array{Float64}(undef, length(detectable_prns))
tracking_snr = Array{Float64}(undef, length(detectable_prns))
tracking_code_phase_err = Array{Float64}(undef, length(detectable_prns))
for i in 1:length(detectable_prns)
    prn = detectable_prns[i]
    ax = fig.add_subplot(2, 4, i)
    tick_params(left=false, bottom=false, labelleft=false, labelbottom=false)
    title("PRN $prn")
    acqresults, trackresults = process(data, define_l1ca_code_type(), prn;
                                       show_plot=false,
                                       fine_acq_method=:fft,
                                       state_num=2,
                                       dynamickf=true,
                                    #    h₀=h_parms[3], h₋₂=h_parms[1],
                                       dll_b=0.01,
                                       acquisition_T=acquisition_T,
                                       fine_acq_T=fine_acq_T,
                                       fine_acq_sub_T=fine_acq_sub_T,
                                       tracking_T=tracking_T,
                                       use_fine_acq=true,
                                       σω=1,
                                       M=M)
    ax.plot(trackresults.t, real.(trackresults.ZP), linewidth=linewidth)
    ax.plot(trackresults.t, imag.(trackresults.ZP), linewidth=linewidth)
    N = floor(Int, length(trackresults.t)/2)
    tracking_σᵩ[i] = std(trackresults.dphi_meas[N:end]) * 180 / π
    tracking_snr[i] = mean(trackresults.SNR[N:end])
    tracking_code_phase_err[i] = std(trackresults.code_err_meas[N:end])
end
subplots_adjust(hspace=0.4, wspace=0.4,  top=0.88, left=0.02, right=0.98,
                bottom=0.02)
plot_title = string("Course Acq. T=$(floor(Int, acquisition_T*1000))ms    ",
                    "Fine Acq. T=$(floor(Int, fine_acq_T*1000))ms    ",
                    "Tracking T=$(floor(Int, tracking_T*1000))ms    ",
                    "M=$(M)")
suptitle(plot_title)
savefig(string(directory, folder, "ch4_l1ca_acquired_prns.svg"), dpi=300)

for prn in detectable_prns
    acqresults, trackresults, corr_result, SNR_est,
    P_d, above_threshold = process(data,
                                define_l1ca_code_type(),
                                prn;
                                fine_acq_method=:fft,
                                return_corrresult=true,
                                return_Pd=true,
                                state_num=3,
                                #    h₀=h_parms[3],
                                #    h₋₂=h_parms[1],
                                q_a=1,
                                dll_b=1,
                                acquisition_T=2e-3,
                                fine_acq_T=10e-3,
                                fine_acq_sub_T=1e-3,
                                tracking_T=1e-3,
                                show_plot=false,
                                σω=1,
                                M=10,
                                use_fine_acq=true);

    name = "pdfs/ch4_l1ca_prn_$(trackresults.prn).pdf"
    file_name = string(directory, folder, name)
    plot_example_track_results(trackresults; f_if=0, code_length=l1ca_code_length,
                            save_to_file=file_name,
                            bottom=0.075, top=0.95, left=0.125, right=0.875,
                            hspace=0.75, wspace=0.35, marker_size=1,
                            line_width=0.8)
end


text:
const IDType = UInt64

"""
An `Index` represents a single tensor index with fixed dimension `dim`. Copies of an Index compare equal unless their
`tags` are different.

An Index carries a `TagSet`, a set of tags which are small strings that specify properties of the `Index` to help
distinguish it from other Indices. There is a special tag which is referred to as the integer tag or prime
level which can be incremented or decremented with special priming functions.

Internally, an `Index` has a fixed `id` number, which is how the ITensor library knows two indices are copies of a
single original `Index`. `Index` objects must have the same `id`, as well as the `tags` to compare equal.
"""
struct Index{T}
  id::IDType
  space::T
  dir::Arrow
  tags::TagSet
  plev::Int
  function Index{T}(id, space::T, dir, tags, plev) where {T}
    return new{T}(id, space, dir, tags, plev)
  end
end

# Used in NDTensors for generic code,
# mostly for internal usage
Index{T}(dim::T) where {T} = Index(dim)

function Index(id, space::T, dir, tags, plev) where {T}
  return Index{T}(id, space, dir, tags, plev)
end

"""
    Index(dim::Int; tags::Union{AbstractString, TagSet} = "",
                    plev::Int = 0)

Create an `Index` with a unique `id`, a TagSet given by `tags`,
and a prime level `plev`.

# Examples
```jldoctest; filter=r"id=[0-9]{1,3}"
julia> i = Index(2; tags = "l", plev = 1)
(dim=2|id=818|"l")'

julia> dim(i)
2

julia> plev(i)
1

julia> tags(i)
"l"
```
"""
function Index(dim::Int; tags="", plev=0)
  return Index(rand(IDType), dim, Neither, tags, plev)
end

"""
    Index(dim::Integer, tags::Union{AbstractString, TagSet}; plev::Int = 0)

Create an `Index` with a unique `id` and a tagset given by `tags`.

# Examples
```jldoctest; filter=r"id=[0-9]{1,3}"
julia> i = Index(2, "l,tag")
(dim=2|id=58|"l,tag")

julia> dim(i)
2

julia> plev(i)
0

julia> tags(i)
"l,tag"
```
"""
Index(dim::Int,
      tags::Union{AbstractString, TagSet};
      plev::Int = 0) = Index(dim; tags = tags, plev = plev)

"""
    id(i::Index)

Obtain the id of an Index, which is a unique 64 digit integer.
"""
id(i::Index) = i.id

"""
    dim(i::Index)

Obtain the dimension of an Index.

For a QN Index, this is the sum of the block dimensions.
"""
NDTensors.dim(i::Index) = i.space

space(i::Index) = i.space

"""
    dir(i::Index)

Obtain the direction of an Index (In, Out, or Neither).
"""
dir(i::Index) = i.dir

# Used for generic code in NDTensors
NDTensors.dir(i::Index) = dir(i)

"""
    setdir(i::Index, dir::Arrow)

Create a copy of Index i with the specified direction.
"""
function setdir(i::Index, dir::Arrow)
  return Index(id(i),copy(space(i)),dir,copy(tags(i)),plev(i))
end

"""
    tags(i::Index)

Obtain the TagSet of an Index.
"""
tags(i::Index) = i.tags

"""
    plev(i::Index)

Obtain the prime level of an Index.
"""
plev(i::Index) = i.plev

"""
    ==(i1::Index, i1::Index)

Compare indices for equality. First the id's are compared,
then the prime levels are compared, and finally the
tags are compared.
"""
function Base.:(==)(i1::Index, i2::Index)
  return id(i1) == id(i2) &&
         plev(i1) == plev(i2) &&
         tags(i1) == tags(i2)
end

# This is so that when IndexSets are converted
# to Julia Base Sets, the hashing is done correctly
function Base.hash(i::Index, h::UInt)
  return hash((id(i), tags(i), plev(i)), h)
end

"""
    copy(i::Index)

Create a copy of index `i` with identical `id`, `dim`, `dir` and `tags`.
"""
Base.copy(i::Index) = Index(id(i),
                            copy(space(i)),
                            dir(i),
                            tags(i),
                            plev(i))

"""
    sim(i::Index; tags = tags(i), plev = plev(i), dir = dir(i))

Produces an `Index` with the same properties (dimension or QN structure)
but with a new `id`.
"""
sim(i::Index;
    tags=copy(tags(i)),
    plev=plev(i),
    dir=dir(i)) = Index(rand(IDType),
                        copy(space(i)),
                        dir,
                        tags,
                        plev)

# Used for internal use in NDTensors
NDTensors.sim(i::Index) = sim(i)

"""
    dag(i::Index)

Copy an index `i` and reverse its direction.
"""
dag(i::Index) = Index(id(i),
                      copy(space(i)),
                      -dir(i),
                      tags(i),
                      plev(i))

# For internal use in NDTensors
NDTensors.dag(i::Index) = dag(i)

"""
    hastags(i::Index, ts::Union{AbstractString,TagSet})

Check if an `Index` `i` has the provided tags,
which can be a string of comma-separated tags or
a TagSet object.

# Examples
```jldoctest; filter=r"id=[0-9]{1,3}"
julia> i = Index(2, "SpinHalf,Site,n=3")
(dim=2|id=861|"Site,SpinHalf,n=3")

julia> hastags(i, "SpinHalf,Site")
true

julia> hastags(i, "Link")
false
```
"""
hastags(i::Index,
        ts::Union{AbstractString,TagSet}) = hastags(tags(i),ts)

hastags(ts::Union{AbstractString,TagSet}) = x->hastags(x,ts)

"""
    hasplev(i::Index, plev::Int)

Check if an `Index` `i` has the provided prime level.

# Examples
```jldoctest; filter=r"id=[0-9]{1,3}"
julia> i = Index(2; plev=2)
(dim=2|id=543)''

julia> hasplev(i, 2)
true

julia> hasplev(i, 1)
false
```
"""
hasplev(i::Index, pl::Int) = plev(i) == pl

hasplev(pl::Int) = x -> hasplev(x, pl)

"""
    hasid(i::Index, id::ITensors.IDType)

Check if an `Index` `i` has the provided id.

# Examples
```jldoctest; filter=r"id=[0-9]{1,3}"
julia> i = Index(2)
(dim=2|id=321)

julia> hasid(i, id(i))
true

julia> j = Index(2)
(dim=2|id=17)

julia> hasid(i, id(j))
false
```
"""
hasid(ind::Index, i::IDType) = id(ind) == i

hasid(i::IDType) = x -> hasid(x, i)

"""
    settags(i::Index, ts)

Return a copy of Index `i` with
tags replaced by the ones given
The `ts` argument can be a comma-separated
string of tags or a TagSet.

# Examples
```jldoctest; filter=r"id=[0-9]{1,3}"
julia> i = Index(2, "SpinHalf,Site,n=3")
(dim=2|id=543|"Site,SpinHalf,n=3")

julia> hastags(i, "Link")
false

julia> j = settags(i,"Link,n=4")
(dim=2|id=543|"Link,n=4")

julia> hastags(j, "Link")
true

julia> hastags(j, "n=4,Link")
true
```
"""
settags(i::Index, ts) = Index(id(i),
                              copy(space(i)),
                              dir(i),
                              ts,
                              plev(i))

"""
    addtags(i::Index,ts)

Return a copy of Index `i` with the
specified tags added to the existing ones.
The `ts` argument can be a comma-separated
string of tags or a TagSet.
"""
addtags(i::Index, ts) = settags(i, addtags(tags(i), ts))

"""
    removetags(i::Index, ts)

Return a copy of Index `i` with the
specified tags removed. The `ts` argument
can be a comma-separated string of tags or a TagSet.
"""
removetags(i::Index, ts) = settags(i, removetags(tags(i), ts))

"""
    replacetags(i::Index, tsold, tsnew)

If the tag set of `i` contains the tags specified by `tsold`,
replaces these with the tags specified by `tsnew`, preserving
any other tags. The arguments `tsold` and `tsnew` can be
comma-separated strings of tags, or TagSet objects.
"""
replacetags(i::Index,
            tsold,
            tsnew) = settags(i,
                             replacetags(tags(i), tsold, tsnew))

"""
    prime(i::Index, plinc::Int = 1)

Return a copy of Index `i` with its
prime level incremented by the amount `plinc`
"""
prime(i::Index, plinc::Int = 1) = setprime(i, plev(i) + plinc)

"""
    setprime(i::Index, plev::Int)

Return a copy of Index `i` with its
prime level set to `plev`
"""
setprime(i::Index, plev::Int) = Index(id(i),
                                      copy(space(i)),
                                      dir(i),
                                      tags(i),
                                      plev)

"""
    noprime(i::Index)

Return a copy of Index `i` with its
prime level set to zero.
"""
noprime(i::Index) = setprime(i, 0)

"""
    adjoint(i::Index)

Prime an Index using the notation `i'`.
"""
Base.adjoint(i::Index) = prime(i)

"""
    ^(i::Index, pl::Int)

Prime an Index using the notation `i^3`.
"""
Base.:^(i::Index, pl::Int) = prime(i, pl)

"""
Iterating over Index `I` gives the IndexVals `I(1)` through `I(dim(I))`.
"""
function Base.iterate(i::Index, state::Int = 1)
  (state > dim(i)) && return nothing
  return (i(state), state+1)
end

# This is a trivial definition for use in NDTensors
NDTensors.outer(i::Index;
                dir = dir(i),
                tags = "",
                plev::Int = 0) =
  sim(i; tags = tags,
         plev = plev,
         dir = dir)

# This is for use in NDTensors
function NDTensors.outer(i1::Index, i2::Index; tags = "")
  return Index(dim(i1) * dim(i2), tags)
end

#
# QN related functions
#

"""
    hasqns(::Index)

Checks of the Index has QNs or not.
"""
hasqns(::Index) = false

"""
    removeqns(::Index)

Removes the QNs from the Index, if it has any.
"""
removeqns(i::Index) = i

#
# IndexVal
#

"""
An IndexVal represents an Index object set to a certain value.
"""
struct IndexVal{IndexT<:Index}
  ind::IndexT
  val::Int
  function IndexVal(i::IndexT,
                    n::Int) where {IndexT <: Index}
    n>dim(i) && throw(ErrorException("Value $n greater than size of Index $i"))
    dim(i)>0 && n<1 && throw(ErrorException("Index value must be >= 1 (was $n)"))
    dim(i)==0 && n<0 && throw(ErrorException("Index value must be >= 1 (was $n)"))
    return new{IndexT}(i, n)
  end
end

"""
    IndexVal(i::Index, n::Int)

    IndexVal(iv::Pair{<:Index, Int})

    (i::Index)(n::Int)

Create an `IndexVal` from a pair of `Index` and `Int`.

Alternatively, you can use the syntax `i(n)`.
"""
IndexVal(iv::Pair{<:Index,Int}) = IndexVal(iv.first,iv.second)

# Help treat a Pair{IndexT, Int} like an IndexVal{IndexT}
const IndexValOrPairIndexInt{IndexT} = Union{IndexVal{IndexT},
                                             Pair{IndexT, Int}}

Base.convert(::Type{IndexVal}, iv::Pair{<:Index,Int}) = IndexVal(iv)

Base.convert(::Type{IndexVal{IndexT}},
             iv::Pair{IndexT, Int}) where {IndexT<:Index} = IndexVal(iv)

Base.getindex(i::Index, j::Int) = IndexVal(i, j)

(i::Index)(n::Int) = IndexVal(i, n)

"""
    ind(iv::IndexVal)

Return the Index of the IndexVal.
"""
NDTensors.ind(iv::IndexVal) = iv.ind

NDTensors.ind(iv::Pair{<:Index,Int}) = iv.first

"""
    val(iv::IndexVal)

Return the value of the IndexVal.
"""
val(iv::IndexVal) = iv.val

val(iv::Pair{<:Index,Int}) = iv.second

"""
    isindequal(i::Index, iv::IndexVal)

    isindequal(i::IndexVal, iv::Index)

    isindequal(i::IndexVal, iv::IndexVal)

Check if the Index and IndexVal have the same indices.
"""
isindequal(i::Index,
           iv::IndexValOrPairIndexInt) = i == ind(iv)

isindequal(iv::IndexValOrPairIndexInt,
           i::Index) = isindequal(i, iv)

isindequal(iv1::IndexValOrPairIndexInt,
           iv2::IndexValOrPairIndexInt) = ind(iv1) == ind(iv2)

plev(iv::IndexValOrPairIndexInt) = plev(ind(iv))

prime(iv::IndexVal,
      inc::Integer = 1) = IndexVal(prime(ind(iv), inc), val(iv))

dag(iv::IndexVal) = IndexVal(dag(ind(iv)), val(iv))

Base.adjoint(iv::IndexVal) = IndexVal(prime(ind(iv)), val(iv))

#
# Printing, reading, and writing
#

function primestring(plev)
  if plev<0
    return " (warning: prime level $plev is less than 0)"
  end
  if plev==0
    return ""
  elseif plev > 3
    return "'$plev"
  else
    return "'"^plev
  end
end

function Base.show(io::IO,
                   i::Index)
  idstr = "$(id(i) % 1000)"
  if length(tags(i)) > 0
    print(io,"(dim=$(space(i))|id=$(idstr)|\"$(tagstring(tags(i)))\")$(primestring(plev(i)))")
  else
    print(io,"(dim=$(space(i))|id=$(idstr))$(primestring(plev(i)))")
  end
end

Base.show(io::IO, iv::IndexVal) = print(io, ind(iv), "=>$(val(iv))")

function readcpp(io::IO, ::Type{Index}; kwargs...)
  format = get(kwargs, :format, "v3")
  if format != "v3"
    throw(ArgumentError("read Index: format=$format not supported"))
  end
  tags = readcpp(io, TagSet; kwargs...)
  id = read(io, IDType)
  dim = convert(Int64, read(io, Int32))
  dir_int = read(io, Int32)
  dir = dir_int < 0 ? In : Out
  read(io, 8) # Read default IQIndexDat size, 8 bytes
  return Index(id, dim, dir, tags)
end

function HDF5.write(parent::Union{HDF5File, HDF5Group},
                    name::AbstractString,
                    I::Index)
  g = g_create(parent, name)
  attrs(g)["type"] = "Index"
  attrs(g)["version"] = 1
  write(g, "id", id(I))
  write(g, "dim", dim(I))
  write(g, "dir", Int(dir(I)))
  write(g, "tags", tags(I))
  write(g, "plev", plev(I))
  if typeof(space(I)) == Int
    attrs(g)["space_type"] = "Int"
  elseif typeof(space(I)) == QNBlocks
    attrs(g)["space_type"] = "QNBlocks"
    write(g,"space",space(I))
  else
    error("Index space type not recognized")
  end
end

function HDF5.read(parent::Union{HDF5File,HDF5Group},
                   name::AbstractString,
                   ::Type{Index})
  g = g_open(parent,name)
  if read(attrs(g)["type"]) != "Index"
    error("HDF5 group or file does not contain Index data")
  end
  id = read(g,"id")
  dim = read(g,"dim")
  dir = Arrow(read(g,"dir"))
  tags = read(g,"tags",TagSet)
  plev = read(g,"plev")
  space_type = "Int"
  if exists(attrs(g),"space_type")
    space_type = read(attrs(g)["space_type"])
  end
  if space_type == "Int"
    space = dim
  elseif space_type == "QNBlocks"
    space = read(g,"space",QNBlocks)
  end
  return Index(id,space,dir,tags,plev)
end


text:
pcount = 4
psize = 2

srand(1)

context("mutation.jl") do
  landscapes = [
    NK.NKLandscape(4, 1),
    NK.NKqLandscape(4, 1, 2),
    NK.NKpLandscape(4, 1, 0.90),
  ]

  for ls = landscapes

    context("$(ls)") do
      context("bwmutate(...)") do
        p = rand(NK.Population, ls, psize)
        np = NK.bwmutate(p, 1.0)
        @fact NK.popsize(np) --> psize
        @fact np --> not(p)
        for i = 1:psize
          @fact np.genotypes[i] --> not(p.genotypes[i])
        end

        mp = rand(NK.MetaPopulation, ls, psize, pcount)
        nmp = NK.bwmutate(mp, 1.0)
        @fact NK.popct(nmp) --> pcount
        @fact NK.popsize(nmp) --> psize
        @fact nmp --> not(mp)
        for i = 1:pcount
          for j = 1:psize
            @fact nmp.populations[i].genotypes[j] --> not(mp.populations[i].genotypes[j])
          end
        end
      end

      context("bwmutate(::Genotype)") do
        g = rand(NK.Genotype, ls)
        ng = NK.bwmutate(g, 1.0)
        @fact g.alleles --> not(ng.alleles)
        g = rand(NK.Genotype, ls)
        ng = NK.bwmutate(g, 0.0)
        @fact g.alleles --> ng.alleles
      end

      context("bsmutate(...)") do
        p = rand(NK.Population, ls, psize)
        np = NK.bsmutate(p, 1.0)
        @fact NK.popsize(np) --> psize
        @fact np --> not(p)
        for i = 1:psize
          @fact np.genotypes[i] --> not(p.genotypes[i])
        end

        mp = rand(NK.MetaPopulation, ls, psize, pcount)
        nmp = NK.bsmutate(mp, 1.0)
        @fact NK.popct(nmp) --> pcount
        @fact NK.popsize(nmp) --> psize
        @fact nmp --> not(mp)
        for i = 1:pcount
          for j = 1:psize
            @fact nmp.populations[i].genotypes[j] --> not(mp.populations[i].genotypes[j])
          end
        end
      end

      context("bsmutate(::Genotype)") do
        g = rand(NK.Genotype, ls)
        ng = NK.bsmutate(g, 1.0)
        @fact g.alleles --> not(ng.alleles)
        g = rand(NK.Genotype, ls)
        ng = NK.bsmutate(g, 0.0)
        @fact g.alleles --> ng.alleles
      end

      context("bsmutate!(::Genotype)") do
        g = rand(NK.Genotype, ls)
        og = deepcopy(g)
        NK.bsmutate!(g)
        @fact g.alleles --> not(og.alleles)
      end
    end
  end
end


text:
using FunctionalCollections: PersistentVector, push, assoc

####################
# markov generator #
####################

struct Markov{T,U} <: Generator{PersistentVector{T},VectorTrace{T,U}}
    kernel::Generator{T,U}
end

function markov(kernel::Generator{T,U}) where {T,U}
    kernel_arg_types = get_static_argument_types(kernel)
    if length(kernel_arg_types) < 3 || kernel_arg_types[1] != Int
        error("markov requires a kernel with arguments (t::Int, state, params...)")
    end
    Markov{T,U}(kernel)
end

function get_static_argument_types(markov::Markov)
    kernel_arg_types = get_static_argument_types(markov.kernel)
    state_type = kernel_arg_types[2]
    params_types = kernel_arg_types[3:end]
    # 1 total number of time steps
    # 2 initial state (this must also be the return type of the kernel, not currently checked)
    # 3 parameters (shared across all time steps)
    [Int, state_type, params_types...]
end

function unpack_args(args::Tuple)
    len = args[1]
    init_state = args[2]
    params = args[3:end]
    (len, init_state, params)
end

function check_length(len::Int)
    #if len < 1
    if len < 0
        error("markov got length of $len < 0")
    end
end

accepts_output_grad(gen::Markov) = false # TODO add support for this

# does not have grad for len
# does not have grad for initial state (TODO add support for this)
# may or may not has grad for parameters, depending the kernel
has_argument_grads(gen::Markov) = (false, false, has_argument_grads(gen.kernel)[3:end]...)

############
# generate #
############

function generate(gen::Markov{T,U}, args, constraints) where {T,U}
    # NOTE: could be strict and check there are no extra constraints
    # probably we want to have this be an option that can be turned on or off?
    (len, init_state, params) = unpack_args(args)
    check_length(len)
    states = Vector{T}(undef, len)
    subtraces = Vector{U}(undef, len)
    weight = 0.
    score = 0.
    state::T = init_state
    is_empty = false
    for key=1:len
        if has_internal_node(constraints, key)
            node = get_internal_node(constraints, key)
        else
            node = EmptyAssignment()
        end
        kernel_args = (key, state, params...)
        (subtrace::U, w) = generate(gen.kernel, kernel_args, node)
        subtraces[key] = subtrace
        weight += w
        call = get_call_record(subtrace)
        states[key] = call.retval
        state = call.retval
        score += call.score
        is_empty = is_empty && !has_choices(subtrace)
    end
    call = CallRecord(score, PersistentVector{T}(states), args)
    trace = VectorTrace{T,U}(PersistentVector{U}(subtraces), call, is_empty)
    (trace, weight)
end

function simulate(gen::Markov{T,U}, args) where {T,U}
    (trace, weight) = generate(gen, args, EmptyAssignment())
    trace
end

##################################
# update, fix_update, and extend #
##################################

# TODO fix_update

struct MarkovChange
    len_changed::Bool
    init_changed::Bool
    params_changed::Bool
end

function extend(gen::Markov{T,U}, args, change::Nothing, trace::VectorTrace{T,U},
                constraints) where {T,U}
    change = MarkovChange(true, true, true)
    extend(gen, args, change, trace, constraints)
end

function extend(gen::Markov{T,U}, args, change::MarkovChange, trace::VectorTrace{T,U},
                constraints) where {T,U}
    (len, init_state, params) = unpack_args(args)
    check_length(len)
    prev_call = get_call_record(trace)
    prev_args = prev_call.args
    prev_len = prev_args[1]
    if len < prev_len
        error("Extend cannot remove addresses or namespaces")
    end
    if change.params_changed
        to_visit = Set{Int}(1:len)
    else
        to_visit = Set{Int}(prev_len+1:len)
    end
    if change.init_changed
        push!(to_visit, 1)
    end
    for (key::Int, _) in get_internal_nodes(constraints)
        push!(to_visit, key)
    end
    subtraces::PersistentVector{U} = trace.subtraces
    states::PersistentVector{T} = trace.call.retval
    weight = 0.
    score = prev_call.score
    is_empty = !has_choices(trace)
    for key in sort(collect(to_visit))
        state = key > 1 ? states[key-1] : init_state
        kernel_args = (key, state, params...)
        if has_internal_node(constraints, key)
            node = get_internal_node(constraints, key)
        else
            node = EmptyAssignment()
        end
        local call::CallRecord{T}
        if key > prev_len
            (subtrace::U, w) = generate(gen.kernel, kernel_args, node)
            call = get_call_record(subtrace)
            score += call.score
            states = push(states, call.retval)
            @assert length(states) == key
            subtraces = push(subtraces, subtrace)
            @assert length(subtraces) == key
        else
            prev_subtrace::U = subtraces[key]
            prev_score = get_call_record(prev_subtrace).score
            kernel_args_change = nothing # TODO permit user to pass through change info to kernel
            (subtrace, w, retchange) = extend(
                gen.kernel, kernel_args, kernel_args_change, prev_subtrace, node)
            call = get_call_record(subtrace)
            score += call.score - prev_score
            @assert length(states) == key
            subtraces = assoc(subtraces, key, subtrace)
            states = assoc(states, key, call.retval)
        end
        weight += w
        is_empty = is_empty && !has_choices(subtrace)
    end
    call = CallRecord(score, states, args)
    trace = VectorTrace(subtraces, call, is_empty)
    (trace, weight, nothing)
end

function update(gen::Markov{T,U}, args, change::Nothing, trace::VectorTrace{T,U},
                constraints) where {T,U}
    change = MarkovChange(true, true, true)
    update(gen, args, change, trace, constraints)
end

function update(gen::Markov{T,U}, args, change::MarkovChange,
                trace::VectorTrace{T,U}, constraints) where {T,U}
    (len, init_state, params) = unpack_args(args)
    check_length(len)
    prev_call = get_call_record(trace)
    prev_args = prev_call.args
    prev_len = prev_args[1]

    subtraces::PersistentVector{U} = trace.subtraces
    states::PersistentVector{T} = trace.call.retval

    # discard deleted applications
    discard = DynamicAssignment()
    if prev_len > len
        for key=len+1:prev_len
            set_internal_node!(discard, key, get_assignment(get_subtrace(trace, key)))
        end
        n_delete = prev_len - len
        for i=1:n_delete
            subtraces = pop(subtraces)
            states = pop(states)
        end
    end
    @assert length(subtraces) == min(prev_len, len)
    @assert length(states) == min(prev_len, len)

    # which retained (not deleted or new) applications to visit
    if change.params_changed
        to_visit = Set{Int}(1:min(prev_len, len))
    else
        to_visit = Set{Int}()
    end
    if change.init_changed
        push!(to_visit, 1)
    end
    for (key::Int, _) in get_internal_nodes(constraints)
        if key <= min(prev_len, len)
            push!(to_visit, key)
        end
    end

    # handle retained applications
    weight = 0.
    score = prev_call.score
    is_empty = !has_choices(trace)
    for key in sort(collect(to_visit))
        state = key > 1 ? states[key-1] : init_state
        kernel_args = (key, state, params...)
        if has_internal_node(constraints, key)
            node = get_internal_node(constraints, key)
        else
            node = EmptyAssignment()
        end
        prev_subtrace::U = subtraces[key]
        prev_score = get_call_record(prev_subtrace).score
        args_change = nothing # NOTE we could propagate detailed change information
        (subtrace, w, kern_discard, retchange) = update(
            gen.kernel, kernel_args, args_change, prev_subtrace, node)
        set_internal_node!(discard, key, kern_discard)
        call = get_call_record(subtrace)
        score += call.score - prev_score
        subtraces = assoc(subtraces, key, subtrace)
        states = assoc(states, key, call.retval)
        weight += w
        is_empty = is_empty && !has_choices(subtrace)
    end

    # handle new applications
    for key=prev_len+1:len
        state = states[key-1]
        kernel_args = (key, state, params...)
        if has_internal_node(constraints, key)
            node = get_internal_node(constraints, key)
        else
            node = EmptyAssignment()
        end
        subtrace::U = assess(gen.kernel, kernel_args, node)
        call = get_call_record(subtrace)
        score += call.score
        weight += call.score
        states = push(states, call.retval)
        subtraces = push(subtraces, subtrace)
        @assert length(states) == key
        @assert length(subtraces) == key
        is_empty = is_empty && !has_choices(subtrace)
    end

    call = CallRecord(score, states, args)
    new_trace = VectorTrace(subtraces, call, is_empty)
    retchange = nothing # NOTE we could provide some information
    (new_trace, weight, discard, retchange)
end

##################
# backprop_trace #
##################

function backprop_trace(gen::Markov{T,U}, trace::VectorTrace{T,U},
                        selection::AddressSet, retval_grad) where {T,U}
    args = get_call_record(trace).args
    (len, init_state, params) = unpack_args(args)
    has_grads = has_argument_grads(gen.kernel)
    if has_grads[1]
        error("Cannot compute gradients for length of markov module")
    end
    if has_grads[2]
        # we ignore this, since the module must be absorbing, so we don't need to compute the grads
    end
    param_has_grad = has_grads[3:end]

    # initialize parameter gradient accumulators
    param_grads = [
        has_grad ? zero(param) : nothing for (param, has_grad) in zip(params, param_has_grad)]

    value_trie = DynamicAssignment()
    gradient_trie = DynamicAssignment()

    # NOTE: order does not matter
    for (key, sub_selection) in get_internal_nodes(selection)
        subtrace = trace.subtraces[key]
        kernel_retval_grad = nothing
        (kernel_arg_grad::Tuple, kernel_value_trie, kernel_gradient_trie) = backprop_trace(
            gen.kernel, subtrace, sub_selection, kernel_retval_grad)
        @assert length(kernel_arg_grad) >= 3
        @assert kernel_arg_grad[1] === nothing
        set_internal_node!(value_trie, key, kernel_value_trie)
        set_internal_node!(gradient_trie, key, kernel_gradient_trie)
        for (i, kernel_param_grad) in enumerate(kernel_arg_grad[3:end])
            if param_has_grad[i]
                @assert kernel_param_grad !== nothing
                param_grads[i] += kernel_param_grad
            end
        end
    end
    arg_grad = (nothing, nothing, param_grads...)
    ((arg_grad...,), value_trie, gradient_trie)
end

export markov
export MarkovChange

max_stars_repo_name: BoltzmannBrain/Gen

max_issues_repo_name: BoltzmannBrain/Gen

max_forks_repo_name: BoltzmannBrain/Gen


text:
# canonical data version 1.1.0
using Test

include("perfect-numbers.jl")

@testset "Perfect numbers" begin

    @testset "Smallest perfect number is classified correctly" begin
        @test classify(6) == "perfect"
    end

    @testset "Medium perfect number is classified correctly" begin
        @test classify(28) == "perfect"
    end

    @testset "Large perfect number is classified correctly" begin
        @test classify(33550336) == "perfect"
    end
end

@testset "Abundant numbers" begin

    @testset "Smallest abundant number is classified correctly" begin
        @test classify(12) == "abundant"
    end

    @testset "Medium abundant number is classified correctly" begin
        @test classify(30) == "abundant"
    end

    @testset "Large abundant number is classified correctly" begin
        @test classify(33550335) == "abundant"
    end
end

@testset "Deficient numbers" begin

    @testset "Smallest prime deficient number is classified correctly" begin
        @test classify(2) == "deficient"
    end

    @testset "Smallest non-prime deficient number is classified correctly" begin
        @test classify(4) == "deficient"
    end

    @testset "Medium deficient number is classified correctly" begin
        @test classify(32) == "deficient"
    end

    @testset "Large deficient number is classified correctly" begin
        @test classify(33550337) == "deficient"
    end

    @testset "Edge case (no factors other than itself) is classified correctly" begin
        @test classify(1) == "deficient"
    end
end

@testset "Invalid inputs" begin

    @testset "Zero is rejected (not a natural number)" begin
        @test_throws DomainError classify(0)
    end

    @testset "Negative integer is rejected (not a natural number)" begin
        @test_throws DomainError classify(-1)
    end
end

max_stars_repo_name: Akshat-mehrotra/julia

max_issues_repo_name: Akshat-mehrotra/julia

max_forks_repo_name: Akshat-mehrotra/julia


text:
@testset "Dynamics" begin
    T = 3
    num_state = 2
    num_action = 1
    num_parameter = 0
    w_dim = [num_parameter for t = 1:T]

    function pendulum(z, u, w)
        mass = 1.0
        lc = 1.0
        gravity = 9.81
        damping = 0.1
        [z[2], (u[1] / ((mass * lc * lc)) - gravity * sin(z[1]) / lc - damping * z[2] / (mass * lc * lc))]
    end

    function euler_explicit(x, u, w)
        h = 0.1
        x + h * pendulum(x, u, w)
    end

    dynamics = Dynamics(euler_explicit, num_state, num_action, num_parameter=num_parameter)
    model = [dynamics for t = 1:T-1]

    x1 = ones(num_state)
    u1 = ones(num_action)
    w1 = zeros(num_parameter)
    X = [x1 for t = 1:T]
    U = [u1 for t = 1:T]
    W = [w1 for t = 1:T]

    dynamics.evaluate(dynamics.evaluate_cache, x1, u1, w1)
    @test norm(dynamics.evaluate_cache - euler_explicit(x1, u1, w1)) < 1.0e-8
    dynamics.jacobian_state(dynamics.jacobian_state_cache, x1, u1, w1)
    dynamics.jacobian_action(dynamics.jacobian_action_cache, x1, u1, w1)

    jac_fd = ForwardDiff.jacobian(a -> euler_explicit(a[1:num_state], a[num_state .+ (1:num_action)], w1), [x1; u1])
    @test norm([dynamics.jacobian_state_cache dynamics.jacobian_action_cache] - jac_fd) < 1.0e-8

    d = [zeros(num_state) for t = 1:T-1]
    for (t, dynamics) in enumerate(model)
        d[t] .= dynamics!(dynamics, X[t], U[t], W[t])
    end
    @test norm(vcat(d...) - vcat([euler_explicit(X[t], U[t], W[t]) for t = 1:T-1]...)) < 1.0e-8

    jx = [zeros(num_state, num_state) for t = 1:T-1]
    ju = [zeros(num_state, num_action) for t = 1:T-1]
    IterativeLQR.jacobian!(jx, ju, model, X, U, W)
    jac_dense = [jx[1] ju[1] zeros(num_state, num_state + num_action);
                 zeros(num_state, num_state + num_action) jx[2] ju[2]]
    @test norm(jac_dense - [jac_fd zeros(model[2].num_state, model[2].num_action + model[2].num_next_state); zeros(model[2].num_next_state, model[1].num_state + model[1].num_action) jac_fd]) < 1.0e-8
end


text:
using Distributions
using Lora

# using Gadfly

nrats, ndays = 30, 5
nvars = 9

a_μ = BasicContUnvParameter(:a_μ, logtarget=(p::Float64, v::Vector) -> logpdf(Normal(0, 100), p), nkeys=nvars)

a_σ2 = BasicContUnvParameter(:a_σ2, logtarget=(p::Float64, v::Vector) -> logpdf(InverseGamma(0.01, 0.01), p), nkeys=nvars)

b_μ = BasicContUnvParameter(:b_μ, logtarget=(p::Float64, v::Vector) -> logpdf(Normal(0, 100), p), nkeys=nvars)

b_σ2 = BasicContUnvParameter(:b_σ2, logtarget=(p::Float64, v::Vector) -> logpdf(InverseGamma(0.01, 0.01), p), nkeys=nvars)

y_σ2 = BasicContUnvParameter(:y_σ2, logtarget=(p::Float64, v::Vector) -> logpdf(InverseGamma(0.01, 0.01), p), nkeys=nvars)

x = Data(:x)

a = BasicContMuvParameter(
  :a,
  logtarget=(p::Vector{Float64}, v::Vector) -> logpdf(MvNormal(v[1]*ones(nrats), sqrt(v[2])), p),
  nkeys=nvars
)

b = BasicContMuvParameter(
  :b,
  logtarget=(p::Vector{Float64}, v::Vector) -> logpdf(MvNormal(v[3]*ones(nrats), sqrt(v[4])), p),
  nkeys=nvars
)

ylogtarget(p::Vector{Float64}, v::Vector) =
  logpdf(MvNormal(repeat(v[7], inner=[ndays])+repeat(v[8], inner=[ndays]).*repeat(v[6], outer=[nrats]), sqrt(v[5])), p)

y = BasicContMuvParameter(:y, logtarget=ylogtarget, nkeys=nvars)

model = GenericModel([a_μ, a_σ2, b_μ, b_σ2, y_σ2, x, a, b, y], isindexed=false)

mcrange = BasicMCRange(nsteps=120000, burnin=100000)

X, = dataset("rats", "age")
X = vec(X)-mean(X)
Y, = dataset("rats", "weight")
v0 = Dict{Symbol, Any}(
  :a_μ=>150.,
  :a_σ2=>1.,
  :b_μ=>10.,
  :b_σ2=>1.,
  :y_σ2=>1.,
  :x=>X,
  :a=>fill(250., nrats),
  :b=>fill(6., nrats),
  :y=>vec(hcat([Y[i, :] for i in 1:30]...))
)

samplers = [
  [
    MH(1.), MH(x::Float64 -> Truncated(Normal(x, 1.), 0, Inf), symmetric=false),
    MH(1.), MH(x::Float64 -> Truncated(Normal(x, 1.), 0, Inf), symmetric=false),
    MH(x::Float64 -> Truncated(Normal(x, 1.), 0, Inf), symmetric=false)
  ];
  fill(MH(1*ones(nrats)), 2); MH(1*ones(nrats*ndays))
]

kk = [:a_μ, :a_σ2, :b_μ, :b_σ2, :y_σ2, :a, :b, :y]

samplers = Dict(zip(kk, samplers))

basicmcjobs =
  Dict([(k, BasicMCJob(
    model,
    samplers[k],
    BasicMCRange(nsteps=1, burnin=0),
    v0,
    pindex=model[k].index,
    jointtarget=true,
    resetpstate=false
  )) for k in kk])

tuners = Dict(zip(kk, fill(RobertsRosenthalMCTuner(0.44), 8)))

gibbsjob = GibbsJob(model, basicmcjobs, mcrange, v0, jointtarget=true, tuner=tuners, verbose=true)

@time run(gibbsjob)

output(gibbsjob)

chains = Dict(gibbsjob)

mean(chains[:b_μ].value)

# plot(x=1:length(chains[:b_μ].value), y=chains[:b_μ].value, Geom.line)


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

export
    SubgraphView, DiSubgraphView,
    supergraph, nodesubgraph, edgesubgraph

struct SubgraphView{T<:UndirectedGraph} <: UndirectedGraph
    graph::T
    nodes::Set{Int}
    edges::Set{Int}
end

struct DiSubgraphView{T<:DirectedGraph} <: DirectedGraph
    graph::T
    nodes::Set{Int}
    edges::Set{Int}
end

# TODO: use trait
AllSubgraphView = Union{SubgraphView,DiSubgraphView}

supergraph(view::AllSubgraphView) = view.graph

"""
    nodesubgraph(graph::UndirectedGraph, nodes::Set{Int}) -> SubgraphView
    nodesubgraph(graph::DirectedGraph, nodes::Set{Int}) -> DiSubgraphView

Generate node-induced subgraph view.
"""
function nodesubgraph(graph::AbstractGraph, nodes::Set{Int})
    # TODO: test for updating in nodes
    edges = Set{Int}()
    for n in nodes
        for (ninc, nadj) in neighbors(graph, n)
            if nadj in nodes
                push!(edges, ninc)
            end
        end
    end
    if graph isa UndirectedGraph
        return SubgraphView(graph, Set(nodes), edges)
    elseif graph isa DirectedGraph
        return DiSubgraphView(graph, Set(nodes), edges)
    end
end

"""
    edgesubgraph(graph::UndirectedGraph, edges::Set{Int}) -> SubgraphView
    edgesubgraph(graph::DirectedGraph, edges::Set{Int}) -> DiSubgraphView

Generate edge-induced subgraph view.
"""
function edgesubgraph(graph::AbstractGraph, edges::Set{Int})
    # TODO: test for updating in edges
    nodes = Set{Int}()
    for e in edges
        union!(nodes, getedge(graph, e))
    end
    if graph isa UndirectedGraph
        return SubgraphView(graph, nodes, Set(edges))
    elseif graph isa DirectedGraph
        return DiSubgraphView(graph, nodes, Set(edges))
    end
end

function neighbors(view::AllSubgraphView, idx::Int)
    (idx in view.nodes) || throw(KeyError(idx))
    return Dict(
        n for n in neighbors(view.graph, idx)
        if n.first in view.edges && n.second in view.nodes
    )
end

function outneighbors(view::DiSubgraphView, idx::Int)
    (idx in view.nodes) || throw(KeyError(idx))
    return Dict(
        n for n in outneighbors(view.graph, idx)
        if n.first in view.edges && n.second in view.nodes
    )
end

function inneighbors(view::DiSubgraphView, idx::Int)
    (idx in view.nodes) || throw(KeyError(idx))
    return Dict(
        n for n in inneighbors(view.graph, idx)
        if n.first in view.edges && n.second in view.nodes
    )
end

function getedge(view::AllSubgraphView, idx::Int)
    idx in view.edges || throw(KeyError(idx))
    return getedge(view.graph, idx)
end

function hasedge(view::AllSubgraphView, u::Int, v::Int)
    u in view.nodes || return false
    v in view.nodes || return false
    return findedgekey(view.graph, u, v) !== nothing
end

function nodeattr(view::AllSubgraphView, idx::Int)
    idx in view.nodes || throw(KeyError(idx))
    return nodeattr(view.graph, idx)
end

function edgeattr(view::AllSubgraphView, idx::Int)
    idx in view.edges || throw(KeyError(idx))
    return edgeattr(view.graph, idx)
end

function edgeattr(view::AllSubgraphView, u::Int, v::Int)
    u in view.nodes || return nothing
    v in view.nodes || return nothing
    k = findedgekey(view.graph, u, v)
    return k === nothing ? nothing : edgeattr(view.graph, k)
end

nodeset(view::AllSubgraphView) = copy(view.nodes)
edgeset(view::AllSubgraphView) = copy(view.edges)

nodecount(view::Union{SubgraphView,DiSubgraphView}) = length(view.nodes)
edgecount(view::Union{SubgraphView,DiSubgraphView}) = length(view.edges)

nodeattrs(view::AllSubgraphView) = nodeattrs(view.graph)
edgeattrs(view::AllSubgraphView) = edgeattrs(view.graph)

nodeattrtype(view::AllSubgraphView) = nodeattrtype(view.graph)
edgeattrtype(view::AllSubgraphView) = edgeattrtype(view.graph)


text:
using RLESUtils, RunUtils, IFTTTUtils

const COMP = gethostname()
const DATA = "nvn_dascfilt"
const CONFIG = "normal"
const I_START = 1
const I_END = 20
const NOTIFY = false

template_gp(i) =
"""
using GrammarExpts,ACASX_GP,RLESUtils,IFTTTUtils
config=configure(ACASX_GP,"$CONFIG","$DATA")
acasx_gp(;seed=$i, outdir=Pkg.dir("GrammarExpts/results/ACASX_Compare/ACASX_GP/seed$i"), config...)
$NOTIFY && sendifttt(;value1="gp,$i", value2="$COMP", value3="$DATA")
"""

template_ge(i) =
"""
using GrammarExpts,ACASX_GE,RLESUtils,IFTTTUtils
config=configure(ACASX_GE,"$CONFIG","$DATA")
acasx_ge(;seed=$i, outdir=Pkg.dir("GrammarExpts/results/ACASX_Compare/ACASX_GE/seed$i"), config...)
$NOTIFY && sendifttt(;value1="ge,$i", value2="$COMP", value3="$DATA")
"""

template_mc(i) =
"""
using GrammarExpts,ACASX_MC,RLESUtils,IFTTTUtils
config=configure(ACASX_MC,"$CONFIG","$DATA")
acasx_mc1(;seed=$i, outdir=Pkg.dir("GrammarExpts/results/ACASX_Compare/ACASX_MC/seed$i"), config...)
$NOTIFY && sendifttt(;value1="mc,$i", value2="$COMP", value3="$DATA")
"""

template_mcts(i) =
"""
using GrammarExpts,ACASX_MCTS,RLESUtils,IFTTTUtils
config=configure(ACASX_MCTS,"$CONFIG","$DATA")
acasx_mcts(;seed=$i, outdir=Pkg.dir("GrammarExpts/results/ACASX_Compare/ACASX_MCTS/seed$i"), config...)
$NOTIFY && sendifttt(;value1="mcts,$i", value2="$COMP", value3="$DATA")
"""

template_ce(i) =
"""
using GrammarExpts,ACASX_CE,RLESUtils,IFTTTUtils
config=configure(ACASX_CE,"$CONFIG","$DATA")
acasx_ce(;seed=$i, outdir=Pkg.dir("GrammarExpts/results/ACASX_Compare/ACASX_CE/seed$i"), config...)
$NOTIFY && sendifttt(;value1="ce,$i", value2="$COMP", value3="$DATA")
"""

notifydone() = sendifttt(;value1="$COMP done")

A = JuliaSource[]
append!(A, [JuliaSource(template_ce(i)) for i=I_START:I_END])
append!(A, [JuliaSource(template_ge(i)) for i=I_START:I_END])
append!(A, [JuliaSource(template_gp(i)) for i=I_START:I_END])
append!(A, [JuliaSource(template_mc(i)) for i=I_START:I_END])
#append!(A, [JuliaSource(template_mcts(i)) for i=I_START:I_END])

#include("script.jl")
#fs=julia_process(A, np); notifydone()


text:
module InvariantEnvironmentPredictionModule

using Distributed
using Hwloc
if length(workers()) == 1
    addprocs(Hwloc.num_physical_cores())  # get ready for parallelism
end
@everywhere begin
    using Distributed, MLJ, DataFrames
    using StatsBase: sample
    using HypothesisTests:  UnequalVarianceZTest, pvalue, EqualVarianceZTest, FisherExactTest, MannWhitneyUTest
end
include("ICPlibrary.jl")
@everywhere forest_model = MLJ.@load RandomForestClassifier pkg = DecisionTree

export InvariantEnvironmentPrediction

@everywhere function fitAndPredict(
    tmRandomForestClassifier,
    X::DataFrame,
    y::CategoricalArray,
)
    mtm = machine(tmRandomForestClassifier, X, y)
    Base.invokelatest( fit!, mtm)
    ŷ = MLJ.predict(mtm, X)
    ŷMode = [mode(ŷ[i]) for i in 1:length(ŷ)]
	return ((ŷMode = ŷMode, mtm = mtm))
end

@everywhere function InvariantEnvironmentPrediction(
    X::DataFrame,   # X must have scientific types
    y::Vector{Float64},
    ExpInd,
    numberOfTargets::Int64,
    numberOfEnvironments::Int64,
    alpha = 0.05,
    verbose = false,
)
    println("InvariantEnvironmentPrediction: ", names(X))
    Xy = copy(X)
    yPermutedIndex = sample(1:length(y), length(y), replace = false, ordered = false)
    Xy.Y = y[yPermutedIndex,1]

    @everywhere pipeRandomForestClassifier = @pipeline FeatureSelector RandomForestClassifier
    cases = [[Symbol(names(Xy)[j]) for j in 1:i] for i in 1:ncol(Xy)]
    r1 = range(pipeRandomForestClassifier, :(feature_selector.features), values=cases)
    tmRandomForestClassifier = TunedModel(
        model = pipeRandomForestClassifier,
        range = r1,
        measures = [cross_entropy, BrierScore()],
        resampling = CV(nfolds = 10)
    )
    # train model with X and permuted Y
    #BrefOnlyX = @spawnat :any fitAndPredict(tmRandomForestClassifier, Xy, E)

    BrefOnlyX = fitAndPredict( tmRandomForestClassifier, Xy, ExpInd)

    # train model with X and Y without permuting
    Xy = copy(X)
    Xy.Y = y[:,1]
    #BrefXY = @spawnat :any fitAndPredict(tmRandomForestClassifier, Xy, E)
    BrefXY = fitAndPredict( tmRandomForestClassifier, Xy, ExpInd)

    #fetchXY = fetch(BrefXY)
    #fetchOnlyX = fetch(BrefOnlyX)

    pval = MannWhitneyUTestTargetEnv(DataFrame(ExpInd),
                                        convert(Vector{Float64}, BrefOnlyX.ŷMode),
                                        convert(Vector{Float64}, BrefXY.ŷMode))

    return ((pvalue = pval, testSet = names(X)))
end #  InvariantEnvironmentPrediction

@everywhere function MannWhitneyUTestTargetEnv(
    ExpInd::DataFrame,  # remember that E was the target
    predictedOnlyX::Vector{Float64},
    predictedXY::Vector{Float64},
)
    residOnlyX = abs.(ExpInd[:,1] - predictedOnlyX)
    residXY = abs.(ExpInd[:,1] - predictedXY)
    pval = pvalue(MannWhitneyUTest(residOnlyX, residXY))
    return (pval)
end # MannWhitneyUTestTargetEnv

end # module


text:
function visualize(io::IO, image::AbstractMatrix{Bool})
    height, width = size(image)
    for i in 1:height
        for j in 1:width
            if image[i, j]
                print(io, "██")
            else
                if iseven(i + j)
                    print(io, "░░")
                else
                    print(io, "▒▒")
                end
            end
        end
        println(io)
    end

    return nothing
end

visualize(image::AbstractMatrix) = visualize(stdout, image)


text:
"""
Creates a Cell line object from JSON. Creates Gene (or Protein) objects on the fly as needed.
"""
function create_cell_line(experiment_id::String, data::Dict{String, Any}; for_prediction::Bool=false)
	experiment = get_experiment(experiment_id)
	log_message("fetching/creating cell line...")
	cl = get_cell_line(experiment_id, data, for_prediction = for_prediction)
	if haskey(data, "views")
		# construct the views
		views = data["views"]
		# should be a dict of views
		for (k,v) in views
			view_type = eval(Meta.parse(k))
			if view_type == RPPA
				key_type = Protein
			else
				key_type = Gene
			end
			data_view = get_dataview!(cl, view_type, DataView{key_type, view_type}(cl.id))
			populate_data_view!(experiment_id, data_view, v)
			add_view!(experiment, view_type)
		end

	end
end

function create_cell_lines(experiment_id::String, data::Dict{String, Any}; for_prediction::Bool=false)
    log_message("importing cell lines")
    if haskey(data, "cell_lines")
        cell_lines = data["cell_lines"]
        for cl in cell_lines
            create_cell_line(experiment_id, cl, for_prediction = for_prediction)
        end
    end
end

function create_genes(experiment_id::String, data::Dict{String, Any})
    log_message("importing genes")
    if haskey(data, "genes")
        for gene in data["genes"]
            add_gene(experiment_id, gene)
        end
    end
end

function create_proteins(experiment_id::String, data::Dict{String, Any})
    log_message("importing proteins")
    if haskey(data, "proteins")
        for protein in data["proteins"]
            add_protein(experiment_id, protein)
        end
    end
end

function create_drugs(experiment_id::String, data::Dict{String, Any})
    log_message("importing drugs")
    if haskey(data, "drugs")
        for drug in data["drugs"]
            add_drug(experiment_id, drug)
        end
    end
end

function create_outcomes(experiment_id::String, data::Dict{String, Any})
    log_message("importing outcomes")
    if !haskey(data, "outcome_type")
        throw(ArgumentError("No outcome type specified."))
    end
    outcome_type = data["outcome_type"]
    # TODO: change to "outcomes"
    if !haskey(data, "outcomes")
        throw(ArgumentError("No outcomes specified, provide a dictionary of drug => outcome under the identifier 'outcomes'"))
    end
    outcomes = data["outcomes"]
	experiment = get_experiment(experiment_id)
    for (drug_name, outcome_data) in outcomes
        log_message("importing outcome for drug $drug_name")
        drug = get_drug!(experiment, drug_name)
        add_outcome(experiment_id, drug, outcome_data, outcome_type = outcome_type)
    end
end

function create_pathways(experiment_id::String, data::Dict{String, Any})
    log_message("importing pathway information")
    if haskey(data, "pathways")
        pathways = data["pathways"]
        for pathway_data in pathways
            add_pathway(experiment_id, pathway_data)
            log_message("done for $(pathway_data["name"])")
        end
    end
end

"""
Populates a data view with ExomeSeq data from JSON data.
If no protein change is provided, we assume the empty string.
"""
function populate_data_view!(experiment_id::String, d::DataView{Gene,ExomeSeq}, data::Vector{Any})
	for entry in data
		gene = add_gene(experiment_id, entry)

		# mandatory values
        if !haskey(entry, "protein_change")
            throw(ArgumentError("You have to specify at least a protein change for ExomeSeq data."))
        end

		# if !(haskey(entry,"reference_mismatch_avg") && haskey(entry,"variant_mismatch_avg")
		# 		&& haskey(entry,"reference_mismatch_sum") && haskey(entry, "variant_mismatch_sum")
		# 		&& haskey(entry, "reference_dist3effective_avg") && haskey(entry, "variant_dist3effective_avg"))
		# 	throw(ArgumentError("For ExomeSeq data, you have to specify at least values for reference_mismatch_avg, variant_mismatch_avg,
		# 							reference_mismatch_sum, variant_mismatch_sum, reference_dist3effective_avg and variant_dist3effective_avg!"))
		# end
        protein_change = entry["protein_change"]

		# optional values
        reference_mismatch_avg = get(entry, "reference_mismatch_avg", 0.)
        reference_mismatch_sum = get(entry, "reference_mismatch_sum", 0.)
        reference_dist3effective_avg = get(entry, "reference_dist3effective_avg", 0.)
        variant_mismatch_avg = get(entry, "variant_mismatch_avg", 0.)
        variant_mismatch_sum = get(entry, "variant_mismatch_sum", 0.)
        variant_dist3effective_avg = get(entry, "variant_dist3effective_avg", 0.)

		num_cosmic = get(entry, "num_cosmic", 0)
		is_cancer_gene = get(entry, "is_cancer_gene", false)
		variant_effect = get(entry, "variant_effect", "")
		nucleotid_change = get(entry, "nucleotid_change", "")
		variant_confidence = get(entry, "variant_confidence", 1.)
		norm_zygosity = get(entry, "norm_zygosity", "")
		norm_reference_count = get(entry, "norm_reference_count", 0.)
		norm_variant_count = get(entry, "norm_variant_count", 0.)
		tumor_zygosity = get(entry, "tumor_zygosity", "")
		tumor_reference_count = get(entry, "tumor_reference_count", 0.)
		tumor_variant_count = get(entry, "tumor_variant_count", 0.)
		details = get(entry, "details", "")

		exome_data = ExomeSeq(protein_change, reference_mismatch_sum=reference_mismatch_sum, reference_mismatch_avg=reference_mismatch_avg,
                                reference_dist3effective_avg=reference_dist3effective_avg, variant_mismatch_sum=variant_mismatch_sum,
                                variant_mismatch_avg=variant_mismatch_avg, variant_dist3effective_avg=variant_dist3effective_avg,
								num_cosmic=num_cosmic, variant_effect=variant_effect,
								nucleotid_change=nucleotid_change, variant_confidence=variant_confidence,
								norm_zygosity=norm_zygosity, norm_reference_count=norm_reference_count, norm_variant_count=norm_variant_count,
								tumor_zygosity=tumor_zygosity, tumor_reference_count=tumor_reference_count, tumor_variant_count=tumor_variant_count,
								details=details)
		add_measurement!(d, gene, exome_data)
	end
end

"""
Populates a data view with Methylation data from JSON.
"""
function populate_data_view!(experiment_id::String, d::DataView{Gene, Methylation}, data::Vector{Any})
	for entry in data
		gene = add_gene(experiment_id, entry)

        # TODO: decide for key name
		# beta_value = entry["illumina_beta_value"]
        beta_value = entry["beta_value"]
		#optionals

        cgct1_value = get(entry, "cgct1", 1)
        cct1_value = get(entry, "cct1", 1)
		methylated_threshold = get(entry, "methylated_threshold", .2)
		illumina_id = get(entry, "illumina_id", "")

		methylation_data = Methylation(beta_value, cgct1_value, cct1_value, methylated_threshold, illumina_id=illumina_id)
		add_measurement!(d, gene, methylation_data)
	end
end

"""
Populate a data view with gene expression data from JSON.
"""
function populate_data_view!(experiment_id::String, d::DataView{Gene, GeneExpression}, data::Vector{Any})
	for entry in data
		gene = add_gene(experiment_id, entry)

		expression_value = Float64(entry["expression_value"])

		add_measurement!(d, gene, GeneExpression(expression_value))
	end
end

"""
Populate a data view with RNASeq data from JSON.
"""
function populate_data_view!(experiment_id::String, d::DataView{Gene, RNASeq}, data::Vector{Any})
	for entry in data
		gene = add_gene(experiment_id, entry)

		expression_value = entry["expression_value"]

		add_measurement!(d, gene, RNASeq(expression_value, expression_status))
	end
end

"""
Populate a data view with RNASeqCall data from JSON.
"""
function populate_data_view!(experiment_id::String, d::DataView{Gene, RNASeqCall}, data::Vector{Any})
	for entry in data
		gene = add_gene(experiment_id, entry)

		expression_status = entry["expression_status"]

		add_measurement!(d, gene, RNASeqCall(expression_value, expression_status))
	end
end

"""
Populate a data view with RPPA data from JSON.
"""
function populate_data_view!(experiment_id::String, d::DataView{Protein, RPPA}, data::Vector{Any})
	for entry in data
		protein = add_protein(experiment_id, entry)

		abundance = entry["protein_abundance"]

		add_measurement!(d, protein, RPPA(abundance))
	end
end

"""
Populate a data view with CNV data form JSON.
"""
function populate_data_view!(experiment_id::String, d::DataView{Gene, CNV}, data::Vector{Any})
	for entry in data
		gene = add_gene(experiment_id, entry)
		cnv_value = Float64(entry["gene_level_cnv"])
		add_measurement!(d, gene, CNV(cnv_value))
	end
end


text:
"""
DR_BlockModel

Block model struture contrains a set of `JuMP.Model` objects, each of which
represents a sub-block model with the information of how to couple these block
models.
"""

mutable struct DR_BlockModel <: DD.AbstractBlockModel
    model::Dict{Int,JuMP.Model} # Dictionary of block models
    coupling_variables::Vector{DD.CouplingVariableRef} # array of variables that couple block models
    variables_by_couple::Dict{Any,Vector{DD.CouplingVariableKey}} # maps `couple_id` to `CouplingVariableKey`

    dual_bound::Float64
    dual_solution::Vector{Float64}

    P_solution::Dict{Int,Float64}

    # TODO: These may be available with heuristics.
    primal_bound::Float64
    primal_solution::Dict{Any, Float64} #coupling_id : value
    combined_weights::Dict{Int, Float64} # block_id : value
    record::Dict{Any, Any}

    function DR_BlockModel()
        return new(
            Dict(),
            [],
            Dict(),
            0.0,
            [],
            Dict(),
            +Inf,
            Dict(),
            Dict(),
            Dict())
    end
end


text:
using OpenHPLJulia
using Test
@testset "OpenHPLJulia.jl" begin
    # Write your own tests here.
    #trunc.(Int,penstock(200, 20; design_by=["height_diff",100], Ld="default" )) ==[400 20 100 44 1 2 101]
end


text:
abstract type DynamicInverterComponent <: DynamicComponent end
abstract type Converter <: DynamicInverterComponent end
abstract type DCSource <: DynamicInverterComponent end
abstract type Filter <: DynamicInverterComponent end
abstract type FrequencyEstimator <: DynamicInverterComponent end
abstract type InnerControl <: DynamicInverterComponent end

abstract type ActivePowerControl <: DeviceParameter end
abstract type ReactivePowerControl <: DeviceParameter end

"""
Parameters of a Outer-Loop controller using a virtual inertia with VSM for active power controller
and a reactive power droop controller.
# Conmutable structor
```julia
VirtualInertiaQDroop(A, R)
```
# Arguments
*  `A`::Float64 : Active power controller using virtual inertia with VSM
*  `R`::Float64 : Reactive power controller using reactive power droop
"""
mutable struct OuterControl{A <: ActivePowerControl, R <: ReactivePowerControl} <:
               DynamicInverterComponent
    active_power::A
    reactive_power::R
    ext::Dict{String, Any}
    states::Vector{Symbol}
    n_states::Int64
end

function OuterControl(
    active_power::A,
    reactive_power::R,
    ext = Dict{String, Any}(),
) where {A <: ActivePowerControl, R <: ReactivePowerControl}
    OuterControl(
        active_power,
        reactive_power,
        ext,
        vcat(active_power.states, reactive_power.states),
        active_power.n_states + reactive_power.n_states,
    )
end

function OuterControl(; active_power, reactive_power, ext = Dict{String, Any}())
    OuterControl(active_power, reactive_power, ext)
end

get_active_power(value::OuterControl) = value.active_power
get_reactive_power(value::OuterControl) = value.reactive_power
get_ext(value::OuterControl) = value.ext
get_states(value::OuterControl) = value.states
get_n_states(value::OuterControl) = value.n_states


text:
module MaxHelpingHandMultiNodeFlagSwitchGate

using ..Ahorn, Maple

@mapdef Entity "MaxHelpingHand/MultiNodeFlagSwitchGate" MultiNodeFlagSwitchGate(x::Integer, y::Integer, width::Integer=16, height::Integer=16,
    inactiveColor::String="5FCDE4", activeColor::String="FFFFFF", finishColor::String="F141DF", flags::String="switch1,switch2", shakeTime::Number=0.5, moveTime::Number=2.0, easing::String="CubeOut",
    sprite::String="block", icon::String="vanilla", resetFlags::Bool=true, canReturn::Bool=true, progressionMode::Bool=false, persistent::Bool=true, pauseTimes::String="", doNotSkipNodes::Bool=false, smoke::Bool=true)

function gateFinalizer(entity)
    x, y = Ahorn.position(entity)

    width = Int(get(entity.data, "width", 8))
    height = Int(get(entity.data, "height", 8))

    entity.data["nodes"] = [(x + width, y)]
end

const textures = ["block", "mirror", "temple", "stars"]
const iconTypes = String["vanilla", "tall", "triangle", "circle", "diamond", "double", "heart", "square", "wide", "winged"]

const easeTypes = String["Linear", "SineIn", "SineOut", "SineInOut", "QuadIn", "QuadOut", "QuadInOut", "CubeIn", "CubeOut", "CubeInOut", "QuintIn", "QuintOut", "QuintInOut", "BackIn", "BackOut", "BackInOut", "ExpoIn", "ExpoOut", "ExpoInOut", "BigBackIn", "BigBackOut", "BigBackInOut", "ElasticIn", "ElasticOut", "ElasticInOut", "BounceIn", "BounceOut", "BounceInOut"]

const placements = Ahorn.PlacementDict(
    "Flag Switch Gate ($(uppercasefirst(texture)), Multi-Node) (max480's Helping Hand)" => Ahorn.EntityPlacement(
        MultiNodeFlagSwitchGate,
        "rectangle",
        Dict{String, Any}(
            "sprite" => texture
        ),
        gateFinalizer
    ) for texture in textures
)

Ahorn.editingOptions(entity::MultiNodeFlagSwitchGate) = Dict{String, Any}(
    "sprite" => textures,
    "icon" => iconTypes,
    "easing" => easeTypes
)

Ahorn.nodeLimits(entity::MultiNodeFlagSwitchGate) = 1, -1

Ahorn.minimumSize(entity::MultiNodeFlagSwitchGate) = 16, 16
Ahorn.resizable(entity::MultiNodeFlagSwitchGate) = true, true

function Ahorn.selection(entity::MultiNodeFlagSwitchGate)
    x, y = Ahorn.position(entity)
    stopX, stopY = Int.(entity.data["nodes"][1])

    width = Int(get(entity.data, "width", 8))
    height = Int(get(entity.data, "height", 8))

    res = [Ahorn.Rectangle(x, y, width, height)]

    for node in get(entity.data, "nodes", ())
        stopX, stopY = Int.(node)
        push!(res, Ahorn.Rectangle(stopX, stopY, width, height))
    end

    return res
end

iconResource = "objects/switchgate/icon00"

function renderGateSwitch(ctx::Ahorn.Cairo.CairoContext, x::Number, y::Number, width::Number, height::Number, sprite::String, icon::String)
    iconSprite = Ahorn.getSprite(icon, "Gameplay")

    tilesWidth = div(width, 8)
    tilesHeight = div(height, 8)

    frame = "objects/switchgate/$sprite"

    for i in 2:tilesWidth - 1
        Ahorn.drawImage(ctx, frame, x + (i - 1) * 8, y, 8, 0, 8, 8)
        Ahorn.drawImage(ctx, frame, x + (i - 1) * 8, y + height - 8, 8, 16, 8, 8)
    end

    for i in 2:tilesHeight - 1
        Ahorn.drawImage(ctx, frame, x, y + (i - 1) * 8, 0, 8, 8, 8)
        Ahorn.drawImage(ctx, frame, x + width - 8, y + (i - 1) * 8, 16, 8, 8, 8)
    end

    for i in 2:tilesWidth - 1, j in 2:tilesHeight - 1
        Ahorn.drawImage(ctx, frame, x + (i - 1) * 8, y + (j - 1) * 8, 8, 8, 8, 8)
    end

    Ahorn.drawImage(ctx, frame, x, y, 0, 0, 8, 8)
    Ahorn.drawImage(ctx, frame, x + width - 8, y, 16, 0, 8, 8)
    Ahorn.drawImage(ctx, frame, x, y + height - 8, 0, 16, 8, 8)
    Ahorn.drawImage(ctx, frame, x + width - 8, y + height - 8, 16, 16, 8, 8)

    Ahorn.drawImage(ctx, iconSprite, x + div(width - iconSprite.width, 2), y + div(height - iconSprite.height, 2))
end

function Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::MultiNodeFlagSwitchGate, room::Maple.Room)
    sprite = get(entity.data, "sprite", "block")
    startX, startY = Int(entity.data["x"]), Int(entity.data["y"])

    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    iconName = get(entity.data, "icon", "vanilla")
    icon = "objects/switchgate/icon00"
    if iconName != "vanilla"
        icon = "objects/MaxHelpingHand/flagSwitchGate/"*iconName*"/icon00"
    end

    for node in get(entity.data, "nodes", ())
        stopX, stopY = Int.(node)

        renderGateSwitch(ctx, stopX, stopY, width, height, sprite, icon)
        Ahorn.drawArrow(ctx, startX + width / 2, startY + height / 2, stopX + width / 2, stopY + height / 2, Ahorn.colors.selection_selected_fc, headLength=6)
        startX, startY = stopX, stopY
    end
end

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::MultiNodeFlagSwitchGate, room::Maple.Room)
    sprite = get(entity.data, "sprite", "block")

    iconName = get(entity.data, "icon", "vanilla")
    icon = "objects/switchgate/icon00"
    if iconName != "vanilla"
        icon = "objects/MaxHelpingHand/flagSwitchGate/"*iconName*"/icon00"
    end

    x = Int(get(entity.data, "x", 0))
    y = Int(get(entity.data, "y", 0))

    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    renderGateSwitch(ctx, x, y, width, height, sprite, icon)
end

end


text:
using Test

@testset "Implementation of WikstromTererlius prover and verifier using Haines paper" begin
    include("proof.jl")
end

@testset "Some Verificatum proof parsing tests" begin
    include("verificatum/runtests.jl")
end

@testset "Verificatum proof verification API" begin
    include("verificatum.jl")
end

@testset "Verifiactum verifier compatible proof generation" begin
    include("vproof.jl")
end

@testset "Elliptic curve tests (under development)" begin
    include("ecverificatum/ecparser.jl")
end


text:
####################################
####################################
type PSF
    nu::Array{Float64}
    nu0::Float64
    mypsf::Array{Float64}
    mypsfadj::Array{Float64}
end

type SKY
    sky::Array{Float64}
    alpha::Array{Float64}
    sig::Float64
    var::Float64
    noise::Array{Float64}
    mydata::Array{Float64}
    sumsky2::Array{Float64}
end
####################################
####################################
type PSF_dirty
    nu::Array{Float64}
    nu0::Float64
    mypsf::Array{Float64}
    mypsfadj::Array{Float64}
end

type SKY_dirty
    sky::Array{Float64}
    mydata::Array{Float64}
    sumsky2::Array{Float64}
end
####################################
####################################
type Algo_param
    nspat::Int64
    nfreq::Int64
    nspec::Int64
    nxy::Int64
    niter::Int64
    lastiter::Int64
    nitermax::Int64
end

type Admm_array
    s::Array{Float64}
    taus::Array{Float64}
    rhos::Float64
    sh::Array{Float64}

    taup::Array{Float64}
    p::Array{Float64}
    t::Array{Float64}
    taut::Array{Float64}
    wlt::Array{Float64}
    x::Array{Float64}

    rhop::Float64
    tauv::Array{Float64}
    v::Array{Float64}
    rhov::Float64
    rhot::Float64
    xmm::Array{Float64}

    spectralwlt::Array{Float64}

    fty::Array{Float64}

    μt::Float64
    μv::Float64
    mueps::Float64
    tt::Float64
    mu::Float64
end

type Admm_array_parallel
    s::Array{Float64}
    taus::Array{Float64}
    rhos::Float64
    sh::Array{Float64}

    taup::SharedArray{Float64}
    p::SharedArray{Float64}
    t::SharedArray{Float64}
    taut::SharedArray{Float64}
    wlt::SharedArray{Float64}
    x::SharedArray{Float64}

    rhop::Float64
    tauv::Array{Float64}
    v::Array{Float64}
    rhov::Float64
    rhot::Float64
    xmm::Array{Float64}

    spectralwlt::Array{Float64}

    fty::Array{Float64}

    μt::Float64
    μv::Float64
    mueps::Float64
    tt::Float64
    mu::Float64
end

type TOOLS
    snr::Array{Float64}
    tol1::Array{Float64}
    tol2::Array{Float64}
    tol3::Array{Float64}
    tol4::Array{Float64}
    tol5::Array{Float64}
    err::Array{Float64}
    errorrec::Array{Float64}
    errorest::Array{Float64}
    errorraw::Array{Float64}
end
####################################
####################################
function init_PSF()
    return PSF([],0.,[],[])
end
function init_SKY()
    return SKY([],[],0.,0.,[],[],[])
end
####################################
####################################
function init_PSF_dirty()
    return PSF_dirty([],0.,[],[])
end
function init_SKY_dirty()
    return SKY_dirty([],[],[])
end
####################################
####################################
function init_Algoparam()
    return Algo_param(0,0,0,0,0,0,0)
end
# function init_Admmarray()
#     return Admm_array([],[],0.,[],[],[],0.,[],[],0.,[],[],0.,[],[],[],[],[],0.,0.,0.,0.,0.)
# end

function init_Admmarray(parallel)
    if parallel == "true"
        return Admm_array_parallel([],[],0.,[],[],[],[],[],[],[],0.,[],[],0.,0.,[],[],[],0.,0.,0.,0.,0.)
    else
        return Admm_array([],[],0.,[],[],[],[],[],[],[],0.,[],[],0.,0.,[],[],[],0.,0.,0.,0.,0.)
    end
end

function init_TOOLS()
    return TOOLS([],[],[],[],[],[],[],[],[],[])
end


text:
@testset "Transform" begin
    include("fourier_transform.jl")
end