text:
using JudiLing
using CSV
using Test
using DataFrames

@testset "eval_SC" begin
    A = [
        1 0 0
        0 1 0
        0 0 1
        1 0 0
        0 1 0
        0 0 1
        1 0 0
        0 1 0
        0 0 1
        0 0 1
    ]

    B = [
        1 0 0
        0 1 0
        0 0 1
        1 0 0
        1 0 0
        0 0 1
        1 0 0
        0 1 0
        0 0 1
        0 0 1
    ]

    C = [
        1 0 0
        0 1 0
        0 0 1
        1 0 0
        1 0 0
        1 0 0
        1 0 0
        1 0 0
        1 0 0
        0 0 1
    ]

    @test JudiLing.eval_SC(A, B) == 0.9
    @test JudiLing.eval_SC(A, C) == 0.6
    @test JudiLing.eval_SC(C, B) == 0.7
    @test JudiLing.eval_SC(A, B, 2) == 0.9
    @test JudiLing.eval_SC(A, B, 3) == 0.9
    @test JudiLing.eval_SC(A, B, 4) == 0.9
    @test JudiLing.eval_SC(A, B, 4) == 0.9
    @test JudiLing.eval_SC(A, C, 2) == 0.6
    @test JudiLing.eval_SC(C, B, 2) == 0.7
end

@testset "eval_SC homophone" begin
    data = DataFrame(
        A = ["aba", "aba", "abc", "abd"],
        B = ["M", "F", "F", "M"],
        C = ["S", "L", "S", "L"],
    )

    cue_obj = JudiLing.make_cue_matrix(
        data,
        grams = 3,
        target_col = :A,
        tokenized = false,
        keep_sep = false,
    )

    n_features = size(cue_obj.C, 2)

    S = JudiLing.make_S_matrix(data, ["B"], ["C"], add_noise = true)

    G = JudiLing.make_transform_matrix(S, cue_obj.C)

    Chat = S * G

    @test JudiLing.eval_SC(Chat, cue_obj.C, data, :A) == 1
    @test JudiLing.eval_SC(Chat, cue_obj.C, data, :A, 2) == 1
    @test JudiLing.eval_SC(Chat, cue_obj.C, data, :A, 3) == 1
end

@testset "eval_SC_loose" begin
    latin = DataFrame(
        Word = ["ABC", "BCD", "CDE", "BCD"],
        Lexeme = ["A", "B", "C", "B"],
        Person = ["B", "C", "D", "D"],
    )

    cue_obj = JudiLing.make_cue_matrix(
        latin,
        grams = 3,
        target_col = :Word,
        tokenized = false,
        keep_sep = false,
    )

    n_features = size(cue_obj.C, 2)
    S = JudiLing.make_S_matrix(
        latin,
        [:Lexeme],
        [:Person],
        ncol = n_features,
    )

    G = JudiLing.make_transform_matrix(S, cue_obj.C)
    Chat = S * G
    F = JudiLing.make_transform_matrix(cue_obj.C, S)
    Shat = cue_obj.C * F

    @test JudiLing.eval_SC_loose(Chat, cue_obj.C, 1) == 0.75
    @test JudiLing.eval_SC_loose(Shat, S, 1) == 0.75
    @test JudiLing.eval_SC_loose(Chat, cue_obj.C, 1, latin, :Word) == 1
    @test JudiLing.eval_SC_loose(Shat, S, 1, latin, :Word) == 1

    for k = 2:4
        @test JudiLing.eval_SC_loose(Chat, cue_obj.C, k) == 1
        @test JudiLing.eval_SC_loose(Shat, S, k) == 1
        @test JudiLing.eval_SC_loose(Chat, cue_obj.C, k, latin, :Word) == 1
        @test JudiLing.eval_SC_loose(Shat, S, k, latin, :Word) == 1
    end
end

@testset "accuracy_comprehension" begin
    latin = DataFrame(
        Word = ["ABC", "BCD", "CDE", "BCD"],
        Lexeme = ["A", "B", "C", "B"],
        Person = ["B", "C", "D", "D"],
    )

    cue_obj = JudiLing.make_cue_matrix(
        latin,
        grams = 3,
        target_col = :Word,
        tokenized = false,
        keep_sep = false,
    )

    n_features = size(cue_obj.C, 2)
    S = JudiLing.make_S_matrix(
        latin,
        [:Lexeme],
        [:Person],
        ncol = n_features,
    )

    G = JudiLing.make_transform_matrix(S, cue_obj.C)
    Chat = S * G
    F = JudiLing.make_transform_matrix(cue_obj.C, S)
    Shat = cue_obj.C * F

    acc_comp = JudiLing.accuracy_comprehension(
        S,
        Shat,
        latin;
        target_col = :Word,
        base = [:Lexeme],
        inflections = [:Person],
    )
end

@testset "eval_acc tests" begin
    latin_train =
        DataFrame(CSV.File(joinpath("data", "latin_mini.csv")))
    cue_obj_train = JudiLing.make_cue_matrix(
        latin_train,
        grams = 3,
        target_col = :Word,
        tokenized = false,
        keep_sep = false,
    )

    latin_val = latin_train[101:150, :]
    cue_obj_val = JudiLing.make_cue_matrix(
        latin_val,
        cue_obj_train,
        grams = 3,
        target_col = :Word,
        tokenized = false,
        keep_sep = false,
    )

    n_features = size(cue_obj_train.C, 2)

    S_train, S_val = JudiLing.make_S_matrix(
        latin_train,
        latin_val,
        ["Lexeme"],
        ["Person", "Number", "Tense", "Voice", "Mood"],
        ncol = n_features,
    )

    G_train = JudiLing.make_transform_matrix(S_train, cue_obj_train.C)

    Chat_train = S_train * G_train
    Chat_val = S_val * G_train

    F_train = JudiLing.make_transform_matrix(cue_obj_train.C, S_train)

    Shat_train = cue_obj_train.C * F_train
    Shat_val = cue_obj_val.C * F_train

    A = cue_obj_train.A

    max_t = JudiLing.cal_max_timestep(latin_train, latin_val, :Word)

    res_train, gpi_train = JudiLing.learn_paths(
        latin_train,
        latin_train,
        cue_obj_train.C,
        S_train,
        F_train,
        Chat_train,
        A,
        cue_obj_train.i2f,
        cue_obj_train.f2i,
        gold_ind = cue_obj_train.gold_ind,
        Shat_val = Shat_train,
        check_gold_path = true,
        max_t = max_t,
        max_can = 10,
        grams = 3,
        threshold = 0.1,
        tokenized = false,
        sep_token = "_",
        keep_sep = false,
        target_col = :Word,
        issparse = :dense,
        verbose = false,
    )

    res_val, gpi_val = JudiLing.learn_paths(
        latin_train,
        latin_val,
        cue_obj_train.C,
        S_val,
        F_train,
        Chat_val,
        A,
        cue_obj_train.i2f,
        cue_obj_train.f2i,
        gold_ind = cue_obj_val.gold_ind,
        Shat_val = Shat_val,
        check_gold_path = true,
        max_t = max_t,
        max_can = 10,
        grams = 3,
        threshold = 0.1,
        is_tolerant = true,
        tolerance = 0.1,
        max_tolerance = 0,
        tokenized = false,
        sep_token = "-",
        keep_sep = false,
        target_col = :Word,
        issparse = :dense,
        verbose = false,
    )

    JudiLing.eval_acc(res_train, cue_obj_train.gold_ind)
    JudiLing.eval_acc(res_train, cue_obj_train)

    JudiLing.eval_acc(res_val, cue_obj_val.gold_ind)
    JudiLing.eval_acc(res_val, cue_obj_val)

    res_train = JudiLing.learn_paths(
        latin_train,
        cue_obj_train,
        S_train,
        F_train,
        Chat_train,
        verbose = false,
    )

    JudiLing.eval_acc(res_train, cue_obj_train.gold_ind)
    JudiLing.eval_acc(res_train, cue_obj_train)
end

@testset "eval_SC train and val tests" begin
    Chat_train = [
    1 0 0
    0 1 0
    0 0 1
    ]

    C_train = [
    1 0 0
    0 1 0
    0 0 1
    ]

    Chat_val = [
    1 1 0
    0 1 1
    1 0 1
    ]

    C_val = [
    1 1 0
    0 1 1
    1 0 1
    ]

    @test JudiLing.eval_SC(Chat_train, C_train) == 1
    @test JudiLing.eval_SC(Chat_val, C_val) == 1
    @test JudiLing.eval_SC(Chat_train, C_train, C_val) == 1
    @test JudiLing.eval_SC(Chat_val, C_val, C_train) == 1

    Chat_train = [
    1 0 0
    0 1 0
    0 0 1
    ]

    C_train = [
    1 0 0
    0 1 0
    0 0 1
    ]

    Chat_val = [
    1 0 0
    0 1 1
    1 0 1
    ]

    C_val = [
    0.9 0.1 0
    0 1 1
    1 0 1
    ]

    @test JudiLing.eval_SC(Chat_train, C_train) == 1
    @test JudiLing.eval_SC(Chat_val, C_val) == 1
    @test JudiLing.eval_SC(Chat_train, C_train, C_val) == 1
    @test JudiLing.eval_SC(Chat_val, C_val, C_train) == round(2/3, digits=4)

    latin_train =
        DataFrame(CSV.File(joinpath("data", "latin_mini.csv")))
    cue_obj_train = JudiLing.make_cue_matrix(
        latin_train,
        grams = 3,
        target_col = :Word,
        tokenized = false,
        keep_sep = false,
    )

    latin_val = latin_train[101:150, :]
    cue_obj_val = JudiLing.make_cue_matrix(
        latin_val,
        cue_obj_train,
        grams = 3,
        target_col = :Word,
        tokenized = false,
        keep_sep = false,
    )

    n_features = size(cue_obj_train.C, 2)

    S_train, S_val = JudiLing.make_S_matrix(
        latin_train,
        latin_val,
        ["Lexeme"],
        ["Person", "Number", "Tense", "Voice", "Mood"],
        ncol = n_features,
    )

    G_train = JudiLing.make_transform_matrix(S_train, cue_obj_train.C)

    Chat_train = S_train * G_train
    Chat_val = S_val * G_train

    F_train = JudiLing.make_transform_matrix(cue_obj_train.C, S_train)

    Shat_train = cue_obj_train.C * F_train
    Shat_val = cue_obj_val.C * F_train

    JudiLing.eval_SC(Chat_train, cue_obj_train.C, cue_obj_val.C)
    JudiLing.eval_SC(Chat_val, cue_obj_val.C, cue_obj_train.C)

    JudiLing.eval_SC(Chat_train, cue_obj_train.C, cue_obj_val.C, latin_train, latin_val, :Word)
    JudiLing.eval_SC(Chat_val, cue_obj_val.C, cue_obj_train.C, latin_val, latin_train, :Word)
end


text:
using GLOQ, Zygote, LinearAlgebra
using GalacticOptim,NLopt,Optim
using Plots

# System parameters for a simple two level open quantum system
N_states = 2; # number of states
freqs = [4.1] # transition frequency in GHz
omegas = 2.0*pi.*freqs # change to angular frequency
gamma1   = [1.0/(55.0*GLOQ.GLOQ_MICRO_SEC)] # Reciprocal of relaxation time - T1 (in units of ns)
gamma2   = [1.0/(15.0*GLOQ.GLOQ_MICRO_SEC)] # Reciprocal of dephasing time - T2 (in units of ns)
omr = 2.0*pi*(4.1 - 5.0e-4) # drive frequency
TC = 2.5*17.0 # total control time in ns needed for a Rabi pulse

# Initial state
initial_state = 0
state_u0 = [0.0;0.0]
state_v0 = [0.0;0.0]
state_u0[initial_state+1] = 1.0

# Duration of the Ramsey experiment, largest dark time given in Microseconds
T_Ramsey = 40.0*GLOQ.GLOQ_MICRO_SEC # convert micro-sec to nano-sec
# total number of dark time samples
N_dark_times = 401
t_dark_times = collect(range(0.0, T_Ramsey, length=N_dark_times))

# Forward solve to generate synthetic data
rho_synthetic_ramsey_u,rho_synthetic_ramsey_v = GLOQ.RamseyForwardSolve(
				 state_u0,state_v0, # initial values, u for the real part, v for the imaginary part
			     omegas,omr, # transition frequencies, drive frequency
				 gamma1,gamma2, # decay and dephasing parameters
				 initial_state, # initial state
				 TC,t_dark_times,N_states) # control time, dark time, total number of states
population_synthetic = GLOQ.get_population(rho_synthetic_ramsey_u)

# Plot the synthetic data
fig=plot(t_dark_times./GLOQ.GLOQ_MICRO_SEC,population_synthetic)
display(fig)

# Define the loss function for the GalacticOptim
# p: phyiscal parameters:
#	 p[1] = transition frequency in GHz
#	 p[2] = gamma2
# dummy_parameter: needed by GalacticOptim, one can just put [] here
function loss(p,dummy_parameter)
	_rho_ramsey_u,_rho_ramsey_v = GLOQ.RamseyForwardSolve(state_u0,state_v0,
				     (2*pi).*[p[1]],omr,
					 gamma1,[p[2]],#gamma1,gamma2,
					 initial_state, # initial state
					 TC,t_dark_times,N_states)
	_population_ramsey = GLOQ.get_population(_rho_ramsey_u)

	_loss = sum(abs2,_population_ramsey-population_synthetic)/N_dark_times
	return _loss
end

plot_callback = function(p,other_args)
	rho_ramsey_u,rho_ramsey_v = GLOQ.RamseyForwardSolve(state_u0,state_v0,
					 (2*pi).*[p[1]],omr,
					 gamma1,[p[2]],#gamma1,gamma2,
					 initial_state, # initial state
					 TC,t_dark_times,N_states)
	population_ramsey = GLOQ.get_population(rho_ramsey_u)
	fig=plot(t_dark_times./GLOQ.GLOQ_MICRO_SEC,population_synthetic,label=["Syn-0" "Syn-1"],
			 line = (:dash,0.0), marker = ([:hex :hex], 5, 0.5)  )
	plot!(fig,t_dark_times./GLOQ.GLOQ_MICRO_SEC,population_ramsey,label=["Opt-0" "Opt-1"]
		 )			#,line = (:dot, 4), marker = ([:hex :hex], 5, 0.1))
	display(fig)
    return false
end

p_true = [freqs;gamma2] # values to generate synthetic data
# initial guess for the optimization
p_initial = [freqs .- 0.5e-4;0.9.*gamma2]

# bounds for the optimization
lower_bound = (0.5).*p_true
upper_bound = (1.5).*p_true

# construct optimization object, use Zygote auto-differentiation to compute the gradient
loss_gradient = GalacticOptim.OptimizationFunction(loss, GalacticOptim.AutoZygote())
opt_prob = GalacticOptim.OptimizationProblem(loss_gradient, p_initial,
										 lb = lower_bound, ub = upper_bound)

println("Optim Fminbox(LBFGS) Optimization starts")
@time sol = GalacticOptim.solve(opt_prob ,Fminbox(LBFGS()),
								cb = plot_callback,
								outer_iterations = 20,
								iterations = 10,
								show_trace=true,
								f_tol = 1e-10,
								outer_f_tol = 1e-10)
println("Optim Fminbox(LBFGS) Optimization done")
# present the solutions
println("\nOptimized results: ",sol.u,
        "\nLoss: ",sol.minimum,
		"\nError: ",sol.u-p_true)

println("NLopt LBFGS Optimization starts")
@time sol_nlopt_LBFGS = GalacticOptim.solve(opt_prob, Opt(:LD_LBFGS,length(p_initial)),
							  maxiters=200,
							  cb = plot_callback,
							  ftol_rel=1e-7)
println("NLopt LBFGS Optimization done")

# present the solution
println("\n\n\n---------------------------\nOptimized results: \n[Transition freq. (GHz), gamma2]\n",sol.u,
        "\nLoss: ",sol.minimum,
		"\nError: ",sol.u-p_true,"\n---------------------------\n")

ext: jl


text:
# %% imports

include("../../../Modules/save_param.jl")
using Statistics
using PyCall
gym = pyimport("gym")
include("../../../Modules/Networks/net_arch.jl")
include("../../../Modules/Poisson/PoissonStateSpike.jl")
include("../../../Modules/simviz.jl")

# %% Read the results

results = load_param("Julia/CartPole/LIF_base/results/results_EI888.jld")

# %% sim functions

function count_spikes(arch::Array{Int64}, spiked::Array{Int64}, layer::Int64)
    count_layer = zeros(arch[layer])
    interval = [sum(arch[1:layer-1])+1, sum(arch[1:layer])]
    spikes_in = filter(x-> interval[1] <= x <= interval[2], spiked)
    spikes_in .-= sum(arch[1:layer-1])
    count_layer[spikes_in] .+= 1
    return count_layer
end

function left_right(count_spikes::Array{Float64}, n_left::Int64, n_right::Int64)
    sl = sum(count_spikes[1:n_left])
    sr = sum(count_spikes[n_left+1:end])
    (sl >= sr) ? action = 0 : action = 1
    return action, sl, sr
end

function borne(x,a,b)
    y = a + 0.5*(b-a)*(1+cos(pi*x/10))
    return y
end

function set_weight!(net::Network, genes_w::Array{Float64})
    w = net.s.e
    for mat_idx in 1:length(w)
        for idx in eachindex(w[mat_idx])
            offset = sum([length(mat) for mat in w[1:mat_idx-1]])
            net.s.e[mat_idx][idx] = (genes_w[idx+offset]+1)/2
        end
    end
end

function set_weight_bound!(net::Network, genes_w::Array{Float64})
    w = net.s.e
    for mat_idx in 1:length(w)
        for idx in eachindex(w[mat_idx])
            offset = sum([length(mat) for mat in w[1:mat_idx-1]])
            net.s.e[mat_idx][idx] = borne(genes_w[idx+offset], -2.0, 4.0)
        end
    end
end

function set_weight!(net::Network, genes_w::Array{Float64}, exc::Array{Int64}, inh::Array{Int64})
    w = net.s.e
    inh_count = 0
    for mat_idx in 1:length(w)
        for idx in eachindex(w[mat_idx])
            offset = sum([length(mat) for mat in w[1:mat_idx-1]]) - inh_count
            offset_neuron = sum(net.arch[1:mat_idx-1])
            (idx%size(w[mat_idx])[1]==0) ? line = size(w[mat_idx])[1] : line = idx%size(w[mat_idx])[1]
            current_neuron = line+offset_neuron
            if current_neuron in exc
                # net.s.e[mat_idx][idx] = abs(genes_w[idx+offset])
                net.s.e[mat_idx][idx] = borne(genes_w[idx+offset], 0.0, 1.1)
            elseif current_neuron in inh
                # net.s.e[mat_idx][idx] = - abs(genes_w[end])
                net.s.e[mat_idx][idx] = borne(genes_w[end], -2.0, 0.0)
                inh_count += 1
            end
        end
    end
end

function ei_idx(arch::Array{Int64})
    inh = Int64[]
    exc = Int64[]
    n_inh = [div(arch[i],4) for i in 1:length(arch)-1]
    push!(n_inh, 0)
    for l in 1:length(arch)
        for neuron in 1:arch[l]
            offset = sum(arch[1:l-1])
            if neuron <= n_inh[l]
                push!(inh, neuron+offset)
            else
                push!(exc, neuron+offset)
            end
        end
    end
    return exc, inh
end

function sig(z::Float64, alpha::Float64)
    return 1.0/(1.0+exp(-z/alpha))
end

function norm_sig(obs::Array{Float64})
    pos = obs[1] # in [-2.4,2.4]
    theta = mod(obs[3]+pi,2pi) - pi # in [-pi, pi]
    vel = [obs[2], obs[4]] # in [-15,15]
    alpha_pos = 2.4/5
    alpha_theta = pi/5
    alpha_vel = 15/5
    pos_scaled = sig(pos,alpha_pos)
    theta_scaled = sig(theta,alpha_theta)
    vel_scaled = sig.(vel,alpha_vel)
    obs_scaled = [pos_scaled, vel_scaled[1], theta_scaled, vel_scaled[2]]
    return obs_scaled
end

# %% main functions

function play_episode(net::Network, matalpha::Array{Float64,2}, render::Bool = false, n_steps::Int64 = 200)
    env = gym.make("CartPole-v1")
    obs = env.reset()
    #rng_action = Random.MersenneTwister(0)
    win_size = 40 # ms
    spike_hist = Array{Int64}[]
    tot_reward = 0
    N_rand = 0
    N_tot = 0
    @inbounds for j in 1:n_steps
        count_tot = zeros(net.arch[end])
        #obs_net = matalpha*norm_sig(obs)
        obs_net = matalpha*norm_obs(obs)
        spikes = int_spikes(obs_net) # state2spike function
        net.v .= net.param["c"]
        @inbounds for msec in 1:win_size
            spiked = step_LIF!(net, spikes, false)
            influence = 2
            if msec >= length(net.arch)-1
                count_tot .+= influence*count_spikes(net.arch, spiked, length(net.arch))
                influence *= 0.95
                push!(spike_hist,spiked)
            end
        end
        n_left = div(net.arch[end],2)
        n_right = net.arch[end] - n_left
        action, sl, sr = left_right(count_tot, n_left, n_right)
        if sl == sr
            action = rand(0:1)
            tot_reward -= 0.9
            N_rand += 1
        end
        N_tot += 1
        #push!(score_hist[1], sl)
        #push!(score_hist[2], sr)
        obs_new, reward, done, _ = env.step(action)
        if render==true
            env.render()
        end
        tot_reward += 1
        if done
            break
        end
        obs = obs_new
    end
    println("Rand = ", N_rand, "/", N_tot)
    env.close()
    env = nothing
    Base.GC.gc()
    return tot_reward, spike_hist
end

function sim(genes, render::Bool = false, n_steps::Int64 = 200)
    genes_alpha = map(x->borne(x,-1.0, 1.0), genes[1:32])
    genes_w = genes[33:end]
    matalpha = reshape(genes_alpha, (8,4))
    #params = generate(genes_p,tunekeys)
    params = YAML.load(open("Julia/CartPole/LIF_base/EI/cfglifEI.yml"))
    arch = [8,8,8]
    n = Network(arch, params)
    exc, inh = ei_idx(arch)
    set_weight!(n,genes_w, exc, inh)
    reward, spike_hist = play_episode(n, matalpha, render, n_steps)
    return reward, spike_hist, n.arch
end

# %% test robustness

function test_random(results::Dict)
    mean_reward = Float64[]
    for i in 1:length(results["fit"])
        genes = results["weights"][i]
        tmp_reward = []
        for ep in 1:30
            reward, _, _ = sim(genes)
            push!(tmp_reward, reward)
        end
        push!(mean_reward, mean(tmp_reward))
    end
    return mean_reward
end

function test_random(results::Dict, indiv::Int64)
    genes = results["weights"][indiv]
    tmp_reward = []
    for ep in 1:30
        reward, _, _ = sim(genes)
        push!(tmp_reward, reward)
    end
    mean_r = mean(tmp_reward)
    return mean_r
end

means = test_random(results)

# %% visualize the activity in the network

function activiz(results::Dict, indiv::Int64)
    genes = results["weights"][indiv]
    r, spike_hist, arch = sim(genes)
    spike_histogram(arch, spike_hist)
end

for idx in [19,20]
    activiz(results, idx)
end

# %% run one episode

genes = results["weights"][20]
r, spike_hist, arch = sim(genes, false, 500)

# %%
using Plots

for indiv in 19:20
    display(histogram(map(x->borne(x,0.0, 1.1), results["weights"][indiv][33:end-1]), xlabel = "indiv $(indiv)",  bins = collect(0.0:0.1:1.1)))
    println("Inh weight $(indiv)= $(borne(results["weights"][indiv][end], -2.0, 0.0))")
end

# %% visualize matalpha
gr()
histogram(map(x->borne(x,-1.0,1.0),results["weights"][19][1:32]), bins = collect(-1.0:0.2:1.0))

# %% definitve matalpha try
alpha_ar = map(x->borne(x,-1.0,1.0),results["weights"][19][1:32])

@show alpha_ar
# %% study matalpha

matalpha = reshape(alpha_ar, (8,4))
@show matalpha

max_stars_repo_name: SuReLI/RL-STDP

max_issues_repo_name: SuReLI/RL-STDP

max_forks_repo_name: SuReLI/RL-STDP


text:
todisplay = VERSION <= VersionNumber("1.2") ? true : false

using LinearAlgebra
using SparseArrays
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using CAGD
using Random
if todisplay
    using ViewerGL
    GL = ViewerGL
end

# Data 0
n = 3;
rngseed = 1234;
maxVisual = 100;
atol = 1e-6;

# Data 1
n = 7;
rngseed = 1111;
maxVisual = 10;
atol = 1e-4;

# Data 2
n = 8;
rngseed = 4321;
maxVisual = 10;
atol = 1e-4;

function generateRotCubes(n = 3; seed = 1234)
    rng = MersenneTwister(seed);
    V,(_,EV,FV,_) = Lar.cuboid([1,1,1],true,[-1,-1,-1])
    cube = V,FV,EV
    cubos = []
    for k = 1:n
        angles = rand(rng, 3) .* 2π
        R = Lar.r(angles[1],0,0) * Lar.r(0,angles[2],0) * Lar.r(0,0,angles[3])
        push!(cubos, Lar.apply(R, cube))
    end
    V,FV,EV = Lar.struct2lar(Lar.Struct(cubos))
    if todisplay  GL.VIEW([ GL.GLGrid(V,FV,GL.COLORS[1],0.5) ])  end

    model = CAGD.Model(V);
    CAGD.addModelCells!(model, 1, EV, signed = true)
    cFE = convert(Lar.ChainOp, Lar.coboundary_1(model.G, FV, EV))
    CAGD.addModelCells!(model, 2, cFE)

    return model
end

#== UNCOMMENT for step-by-step pipeline ==#
#==

model = generateRotCubes(n, seed = rngseed)

atol = 1e-6;

if todisplay
    GL.VIEW([
        GL.GLPol(model.G, Lar.cop2lar(model.T[1]), GL.COLORS[1])
        GL.GLFrame
    ]);
end

split_model = CAGD.pairwise_decomposition(model, atol = atol)

if todisplay
    GL.VIEW([
        GL.GLPol(model.G, Lar.cop2lar(model.T[1]), GL.COLORS[1])
        GL.GLFrame
    ]);
end

congr_model = CAGD.mergeModelVertices(split_model, err=atol, signed_merge=true)

if todisplay
    GL.VIEW([
        GL.GLPol(congr_model.G, Lar.cop2lar(congr_model.T[1]), GL.COLORS[1])
        GL.GLFrame
    ]);
end

A_model = deepcopy(congr_model);
FC, bicon_comps = CAGD.tgw(congr_model, 3)
CAGD.addModelCells!(A_model, 3, convert(Lar.ChainOp, FC'))

if todisplay
    GL.VIEW([
        GL.GLPol(A_model.G, Lar.cop2lar(A_model.T[1]), GL.COLORS[1])
        GL.GLFrame
    ]);
end
==#

model = generateRotCubes(n, seed = rngseed)
A_model = CAGD.spatial_arrangement(model, atol = atol, outerCell = true)

if todisplay
    V,CVs,FVs,EVs = Lar.pols2tria(A_model.G, A_model.T[1], A_model.T[2], A_model.T[3]);
    GL.VIEW(GL.GLExplode(V, EVs,         1.5,1.5,1.5,  99,1));
    GL.VIEW(GL.GLExplode(V, FVs,         1.5,1.5,1.5,  99,1));
    GL.VIEW(GL.GLExplode(V, CVs[2:end],  1.0,1.0,1.0,  99,1.0));
    GL.VIEW(GL.GLExplode(V, CVs[2:end],  1.0,1.0,1.0,  99,0.5));
    GL.VIEW(GL.GLExplode(V, CVs[2:end],  5.0,5.0,5.0,  99,0.5));
    # Takes the 3-cell(s) with more 2-cells
    idxs = SparseArrays.sparse(length.(CVs) .== max(length.(CVs[2:end])...)).nzind;
    GL.VIEW([
        GL.GLExplode(V,CVs[idxs], 1.0,1.0,1.0, 2, 0.9)
        GL.GLExplode(V,CVs[1:1], 1.0,1.0,1.0, 1, 0.1)
    ])
    GL.VIEW([
        GL.GLFrame
        GL.GLExplode(V, CVs[idxs], 1.0,1.0,1.0, 2, 0.9)
    ])
end

# Visualise first maxVisual 3-cells within the complete structure
if todisplay
    for i = 2 : min(length(CVs), maxVisual)
        GL.VIEW([
            GL.GLFrame
            GL.GLExplode(V,CVs[i:i], 1.0,1.0,1.0, 2, 0.9)
            GL.GLExplode(V,CVs[1:1], 1.0,1.0,1.0, 1, 0.1)
        ])
    end
end


text:
module PlmDCA
using SharedArrays,Distributed,Printf, LinearAlgebra, Statistics, DelimitedFiles
using NLopt
import DCAUtils: read_fasta_alignment, remove_duplicate_sequences, compute_weights, add_pseudocount, compute_weighted_frequencies
using LoopVectorization

export PlmOut, plmdca, plmdca_asym, plmdca_sym, plmdca_asym

include("types.jl")
include("utils.jl")
include("plmdca_asym.jl")
include("plmdca_sym.jl")
include("decimation_sym.jl")
include("mi.jl")
end

max_stars_repo_name: francescocaredda/PlmDCA

max_issues_repo_name: francescocaredda/PlmDCA

max_forks_repo_name: francescocaredda/PlmDCA


text:
using InfiniteOpt, Ipopt, Plots
using KNITRO

# Parmeters
# Time span
t0 = 0
tf = 1
# Initial values
x1 = 1

# Model
model = InfiniteModel(KNITRO.Optimizer)

# infinite_parameter
@infinite_parameter(model, t ∈ [t0, tf], num_supports = 101)

# variable
@variable(model, x ≥ 0, Infinite(t))

@variable(model, u, Infinite(t), start = 0)

# objective
@objective(model, Min, ∫(u^2, t))

# constraint
@constraint(model, x(0) == x1)
@constraint(model, x_constr, ∂(x, t) == x + u)

# Optimization

print(model)

optimize!(model)

# Visulization
ts = value(t)
u = value(u)
x = value(x)

display(plot(ts, u))
display(plot(ts, x))

using DiffEqFlux, DifferentialEquations
const solveeq = DifferentialEquations.solve

tspan = (0.0f0, 20.0f0)
ann = FastChain(FastDense(1, 32, tanh), FastDense(32, 1), (x, p) -> 0.9 * sigmoid.(x))
θ = initial_params(ann)
function sir_nn(du, u, p, t)
    S, E, I, N = u
    du[1] = b * N - d * S - c * S * I - ann([t], p)[1] * S
    du[2] = c * S * I - (e + d) * E
    du[3] = e * E - (g + a + d) * I
    du[4] = (b - d) * N - a * I
end
u0 = [s0, e0, i0, n0]
ts = Float32.(collect(0.0:0.2:tspan[2]))
prob = ODEProblem(sir_nn, u0, tspan, θ)
sol_init = solveeq(prob, Vern9(), abstol = 1e-10, reltol = 1e-10)
plot(sol_init)
function predict_adjoint(θ)
    Array(solveeq(prob, Vern9(), p = θ, saveat = ts, sensealg = InterpolatingAdjoint(autojacvec = ReverseDiffVJP(true))))
end
function loss_adjoint(θ)
    x = predict_adjoint(θ)
    mean(abs, A .* x[3, :]) + mean(abs2, [first(ann([t], θ)) for t in ts])
end
l = loss_adjoint(θ)
cb = function (θ, l)
    println(l)
    p = plot(solveeq(remake(prob, p = θ), Tsit5(), saveat = 0.2), lw = 3)
    #plot!(p, ts, [first(ann([t], θ)) for t in ts], label = "u(t)", lw = 3)
    display(p)
    return false
end
# Display the ODE with the current parameter values.
cb(θ, l)
loss1 = loss_adjoint(θ)
res1 = DiffEqFlux.sciml_train(loss_adjoint, θ, ADAM(0.05), cb = cb, maxiters = 100)
res2 = DiffEqFlux.sciml_train(loss_adjoint, res1.u,
    BFGS(initial_stepnorm = 0.01), cb = cb, maxiters = 100,
    allow_f_increases = false)

# Loss saving

opt_method_loss = 13.90369069230782
nn_method_loss = 13.711496454176164


text:
############################################################################################
#Import External Packages

############################################################################################
# Load sub-container
include("diagnostics.jl")
include("statistics.jl")
include("shuffle.jl")
include("resample.jl")
include("weights.jl")
include("logging.jl")

############################################################################################
#export


text:
@testset "Blue Style" begin
    @testset "ops" begin
        # `//` and `^` are binary ops without whitespace around them
        @test bluefmt("1 // 2 + 3 ^ 4") == "1//2 + 3^4"
        @test bluefmt("a.//10") == "a .// 10"
    end

    @testset "nest to one line" begin
        str_ = """
        var = [arg1, #com
        arg2]
        """
        str = """
        var = [
            arg1, #com
            arg2,
        ]
        """
        @test bluefmt(str_) == str

        str_ = """
        var = (arg1,
        arg2)
        """
        str = """
        var = (arg1, arg2)
        """
        @test bluefmt(str_, 4, 18) == str

        str_ = """
        var = {arg1,
        arg2}
        """
        str = """
        var = {
            arg1, arg2
        }
        """
        @test bluefmt(str_, 4, 17) == str
        @test bluefmt(str_, 4, 14) == str

        str_ = """
        var = call(arg1,
        arg2)
        """
        str = """
        var = call(
            arg1, arg2
        )
        """
        @test bluefmt(str_, 4, 14) == str

        str = """
        var = call(
            arg1,
            arg2,
        )
        """
        @test bluefmt(str_, 4, 13) == str

        str_ = """
        var = ref[arg1,
        arg2]
        """
        str = """
        var = ref[
            arg1, arg2
        ]
        """
        @test bluefmt(str_, 4, 14) == str

        str = """
        var = ref[
            arg1,
            arg2,
        ]
        """
        @test bluefmt(str_, 4, 13) == str

        str_ = """
        var = ABC{arg1,
        arg2}
        """
        str = """
        var = ABC{
            arg1,arg2
        }
        """
        @test bluefmt(str_, 4, 13) == str

        str = """
        var = ABC{
            arg1,
            arg2,
        }
        """
        @test bluefmt(str_, 4, 12) == str

        str_ = """
        var = @call(arg1,
        arg2)
        """
        str = """
        var = @call(
            arg1, arg2
        )
        """
        @test bluefmt(str_, 4, 14) == str

        str = """
        var = @call(
            arg1,
            arg2
        )
        """
        @test bluefmt(str_, 4, 13) == str

        str_ = """
        function long_name_of_function_because_i_am_writing_an_example(
            arg1, arg2, arg3, arg4, arg5, arg6
        )
            # code
        end
        """
        @test bluefmt(str_, 4, 38) == str_

        str = """
        function long_name_of_function_because_i_am_writing_an_example(
            arg1,
            arg2,
            arg3,
            arg4,
            arg5,
            arg6,
        )
            # code
        end
        """
        @test bluefmt(str_, 4, 37) == str

    end

    @testset "do not prepend return in `do` blocks" begin
        str = """
        map(arg1, arg2) do x, y
            expr1
            expr2
        end"""
        @test bluefmt(str, always_use_return = true) == str
    end

    @testset "separate kw args with semicolon" begin
        str_ = "xy = f(x, y=3)"
        str = "xy = f(x; y = 3)"
        @test bluefmt(str_) == str

        str_ = "xy = f(x=1, y=2)"
        str = "xy = f(; x = 1, y = 2)"
        @test bluefmt1(str_) == str
        @test bluefmt(str_) == str
        @test bluefmt(str) == str

        str_ = "xy = f(x=1; y=2)"
        @test bluefmt1(str_) == str
        @test bluefmt(str_) == str

        str = """
        function g(x, y = 1)
            return x + y
        end
        macro h(x, y = 1)
            nothing
        end
        shortdef1(MatrixT, VectorT = nothing) = nothing
        shortdef2(MatrixT, VectorT = nothing) where {T} = nothing
        """
        @test bluefmt1(str) == str
        @test bluefmt(str) == str
    end
end


text:
@testset "neighbourhood" begin
    @testset "LocalNeighbourhood" begin
        ln = pso.LocalNeighbourhood(5)
        @test_throws ErrorException ln[1] = [1, 2, 3]
        @test IndexStyle(ln) == IndexLinear()
        @test ln[2] == [1, 2, 3]
        @test ln[5] == [4, 5, 1]
        @test length(ln) == 5
        @test_throws AssertionError pso.LocalNeighbourhood(5, 4)
    end

    @testset "GlobalNeighbourhood" begin
        gn = pso.GlobalNeighbourhood(10)
        @test gn[1] == collect(1:10)
        @test gn[4] == collect(1:10)
        @test gn[10] == collect(1:10)
        @test size(gn) == (10,)
    end

    @testset "HirachicalNeighbourhood" begin
        hn = pso.HierachicalNeighbourhood(21, 4)
        @testset "bestchild" begin
            @test pso.bestchild(hn.graph_to_particle, hn.childs[1], vcat(1.0, 0.0, fill(1.0, 19)), <) == 2
            @test pso.bestchild(hn.graph_to_particle, hn.childs[1], vcat(0.0, 1.0, fill(0.0, 19)), <) == 3
        end
        @testset "swapifnecessary" begin
            pso.swapifnecessary!(hn, fill(0.0, 21), 2, 1, 2, 1, <)
            @test hn.particle_to_graph == collect(1:21)
            @test hn.graph_to_particle == collect(1:21)
            pso.swapifnecessary!(hn, vcat(1.0, fill(0.0, 20)), 2, 1, 2, 1, <)
            @test hn.particle_to_graph == vcat(2, 1, collect(3:21))
            @test hn.graph_to_particle == vcat(2, 1, collect(3:21))
        end
        childs = pso.getchilds(4, 2)
        @test childs == [[2, 3], [4], Int[], Int[]]
        @test pso.getparents(childs) == [1, 1, 1, 2]
        childs = pso.getchilds(8, 3)
        @test childs == [[2, 3, 4], [5, 8], [6], [7], Int[], Int[], Int[], Int[]]
        @test pso.getparents(childs) == [1, 1, 1, 1, 2, 3, 4, 2]
        hn1 = pso.HierachicalNeighbourhood(21, 4)
        @test hn1[1] == [1]
        @test hn1[3] == [1]
        @test hn1[10] == [2]
        @test hn1[15] == [3]
        pso.rearrange!(hn1, vcat(1.0, fill(0.0, 19), -1.0), <)
        @test hn1[2] == [2]
        @test hn1[3] == [2]
        @test hn1[21] == [2]
        @test hn1.graph_to_particle[hn1.particle_to_graph] == collect(1:21)
        @test hn1.particle_to_graph[hn1.graph_to_particle] == collect(1:21)
        pso.rearrange!(hn1, vcat(1.0, fill(0.0, 19), -1.0), <)
        @test hn1[21] == [21]
        @test hn1[3] == [21]
        @test hn1[2] == [21]
        @test hn1.graph_to_particle[hn1.particle_to_graph] == collect(1:21)
        @test hn1.particle_to_graph[hn1.graph_to_particle] == collect(1:21)
        @test length(hn1) == 21
        hn2 = pso.HierachicalNeighbourhoodByTreeHeight(3, 4)
        @test hn2[1] == [1]
        @test hn2[3] == [1]
        @test hn2[21] == [5]
        @test size(hn2) == (21,)
        hn3 = pso.HierachicalNeighbourhood(40, 3)
        @test hn3[9] == [3]
        @test hn3[17] == [6]
        @test hn3[39] == [10]
        @test hn3[40] == [13]
        hn4 = pso.HierachicalNeighbourhoodByTreeHeight(5, 2)
        @test hn4[9] == [6]
        @test hn4[17] == [12]
        @test hn4[23] == [15]
    end

end


text:
module JuliaToFortran

using   StaticArrays,
        MacroTools,
        Test

export  translate_expr,
        create_subroutine,
        create_module

include("julia_to_fortran.jl")

end # module


text:
function f(x :: Any, y :: Any) :: Bool return true end
function f(x :: Bool, y :: Any) :: Int64 return 0 end
function bar(x :: Any, y :: Any) return g(f(x, y)) end
function g(x :: Bool) 42 end
println(bar(0, true))
println(f(true, 0))

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
import SpecialFunctions: logbeta

@node Beta Stochastic [ out, (a, aliases = [ α ]), (b, aliases = [ β ]) ]

@average_energy Beta (q_out::Any, q_a::Any, q_b::Any) = logbeta(mean(q_a), mean(q_b)) - (mean(q_a) - 1.0) * mean(log, q_out) - (mean(q_b) - 1.0) * mean(mirrorlog, q_out)


text:
# ---
# title: The BM3D(sparse 3D transform-domain collaborative filtering) denoising algorithm.
# id: bm3d-demo
# cover: assets/bm3d_cover.png
# date: 2021-08-04
# author: "[Longhao Chen](https://github.com/Longhao-Chen)"
# description: Use the BM3D denoising algorithm to reduce gaussian noise
# ---

using ImageNoise
using TestImages, ImageShow, ImageCore, ImageQualityIndexes, ImageTransformations
using FileIO, Random #src

# First, load an image and and add gaussian noise to it

gray_img = float.(imresize(testimage("cameraman"), ratio=0.5))
n = AdditiveWhiteGaussianNoise(0.1)
noisy_img = apply_noise(gray_img, n)

# Then calling the standard `reduce_noise` API

f_bm3d = BM3D(0.1)
denoised_img = reduce_noise(noisy_img, f_bm3d)

mosaicview(gray_img, noisy_img, denoised_img; nrow=1)

# Get the PSNR using ImageQualityIndexes package:

assess_psnr(gray_img, denoised_img)

mkpath("assets") #src
save("assets/bm3d_cover.png", denoised_img) #src


text:
"""
Annealed Importance Sampling. Estimate the partition function of lnpdf

Args:
    - lnpdf     : target distribution unnomralized log prob.  Potentially
                  vectorized --- takes in z = D-dimensional tensor, outputs
                  joint log prob
    - lnpdf0    : initial distribution log likelihood --- same input and
                  output as lnpdf
    - z0        : D-dimensional initial sample
    - step      : function that takes in a current state and an
                  unnormalized log probability function, performs
                  an MCMC transition (in place)
    - init_logZ : (optional) initial log partition (unused for now)
    - schedule  : temperature schedule, monotonically increasing from 0 to 1
"""
function ais(lnpdf, lnpdf0, step, z0;
             init_logZ=nothing, schedule=nothing, verbose=false)
    # set up initial log partition function for p(z) --- should be ln(1) = 0
    D = length(z0)

    # set up and check schedule
    if schedule == nothing
        schedule = collect(linspace(0, 1, 10))   # default to 10 steps ...
    end
    @assert(isapprox(schedule[1], 0.), "schedule must start with temp = 0.")
    @assert(isapprox(schedule[end], 1.), "schedule must end with temp = 1.")

    # set up function that returns intermediate distributions
    function lnpdf_t(z, t)
        # do this to catch the 0 * -Inf = NaN cases
        if t==0.
            return lnpdf0(z)
        elseif t==1.
            return lnpdf(z)
        end
        return t*lnpdf(z) + (1.-t)*lnpdf0(z)
    end

    # run AIS: track the ratio at each step, and current z
    z = deepcopy(z0)
    llratios = zeros(length(schedule) - 1)
    accepted = 0.
    for ti in 2:length(schedule)
        # unpack previous and current temperatures
        tprev, tcurr = schedule[ti-1], schedule[ti]
        if ti % 25 == 0
            Log.info(@sprintf "  temp %2.3f (%d/%d): logpost = %2.4f  logprior = %2.4f" tcurr ti length(schedule) lnpdf(z) lnpdf0(z))
        end

        # generate zt | zt-1 using tcurr --- this leaves the distribution
        #   lnp_t = tcurr*lnpdf(z) + (1-tcurr)*lnpdf0(0) invariant
        z, _ = step(z, z -> lnpdf_t(z, tcurr))  # returns z, and aux info

        # record log ratio of distributions
        llprev = lnpdf_t(z, tprev)
        llcurr = lnpdf_t(z, tcurr)
        llratios[ti-1] = llcurr - llprev
    end

    # return z (final posterior sample) and weight
    return z, sum(llratios), llratios, float(accepted)/(length(schedule))
end

"""broadcasting logsumexp"""
function logsumexp(mat; dim=1)
    max_score = maximum(mat, dim)
    lse = log.(sum(exp.(broadcast(-, mat, max_score)), dim))
    return broadcast(+, lse, max_score)
end

""" compute a bootstrap sample of the normalizing constant """
function bootstrap_lnZ(lnZ_samps; num_bootstrap=100)
    lnZ    = zeros(num_bootstrap)
    nsamps = length(lnZ_samps)
    for n in 1:num_bootstrap
        boot_sample = [lnZ_samps[rand(1:nsamps)] for n in 1:nsamps]
        lnZ[n]      = logsumexp(boot_sample)[1] - log(nsamps)
    end
    return lnZ
end

""" copy of sigmoid schedule from https://github.com/tonywu95/eval_gen/blob/master/algorithms/ais.py
This is defined as:
      gamma_t = sigma(rad * (2t/T - 1))
      beta_t = (gamma_t - gamma_1) / (gamma_T - gamma_1),
where sigma is the logistic sigmoid. This schedule allocates more distributions near
the inverse temperature 0 and 1, since these are often the places where the distributon
changes the fastest.
"""
function sigmoid_schedule(num_steps; rad=4)
    if num_steps == 1
        return collect(linspace(0, 1, 2))
    end
    t    = collect(linspace(-rad, rad, num_steps))
    sigm = 1. ./ (1. .+ exp.(-t))
    return (sigm - minimum(sigm)) / (maximum(sigm) - minimum(sigm))
end

""" Run AIS multiple times to get multiple posterior samples, and a tighter
marginal likelihood estimate """
function ais_slicesample(logposterior::Function,
                         logprior::Function,
                         prior_sample::Function;
                         schedule::Array{Float64,1}=nothing,
                         num_temps::Int=50,
                         num_samps::Int=10,
                         num_bootstrap::Int=5000,
                         num_samples_per_step::Int=1)
    if schedule == nothing
        schedule = sigmoid_schedule(num_temps; rad=1)
    end

    function step(z, lnpdf)
        for i in 1:num_samples_per_step
            z, llh = MCMC.slicesample(z, lnpdf)
        end
        return z, 0.
    end

    D = length(prior_sample())
    zsamps = zeros(D, num_samps)
    wsamps = zeros(num_samps)
    for n in 1:num_samps
        Log.info(@sprintf "ais samp %d / %d" n num_samps)
        z0   = prior_sample()
        #step = (z, lnpdf) -> MCMC.slicesample(z, lnpdf)
        z, w, llrats, _ = ais(logposterior, logprior, step, z0; schedule=schedule)
        zsamps[:,n] = z
        wsamps[n]   = w
    end

    # estimate the partition function and get a bootstrap confidence interval
    lnZ  = logsumexp(wsamps, dim=1)[1] - log(num_samps)
    lnZs = bootstrap_lnZ(wsamps; num_bootstrap=num_bootstrap)
    res = Dict(:lnZ => lnZ, :lnZ_bootstrap => lnZs,
               :zsamps => zsamps, :lnZsamps => wsamps)
    return res
end


text:
# This file contains code taken from https://github.com/JuliaIntervals/IntervalSpecialFunctions.jl
# which is available under the following MIT license:

# > Copyright (c) 2018: David Sanders.
# >
# > Permission is hereby granted, free of charge, to any person obtaining a copy
# > of this software and associated documentation files (the "Software"), to deal
# > in the Software without restriction, including without limitation the rights
# > to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# > copies of the Software, and to permit persons to whom the Software is
# > furnished to do so, subject to the following conditions:
# >
# > The above copyright notice and this permission notice shall be included in all
# > copies or substantial portions of the Software.
# >
# > THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# > IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# > FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# > AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# > LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# > OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# > SOFTWARE.

using IntervalArithmetic
using IntervalArithmetic: @round, big53

import SpecialFunctions: erf, erfc, erfinv, erfcinv

for f in (:erf, :erfc)
    @eval function($f)(x::BigFloat, r::RoundingMode)
        setrounding(BigFloat, r) do
            ($f)(x)
        end
    end
    @eval ($f)(a::Interval{Float64}) = convert(Interval{Float64}, ($f)(big53(a)))
end

function erf(a::Interval{T}) where T
    isempty(a) && return a
    @round( erf(a.lo), erf(a.hi) )
end

function erfc(a::Interval{T}) where T
    isempty(a) && return a
    @round( erfc(a.hi), erfc(a.lo) )
end


text:
include("template.jl")

type taResults
  recovered_key::String
  templates::Array
  res_dir::String
end

#this will basically use the template to figure out the keys
function ta(templates, trace_data, cipher_type)
  unit_size = trace_data.unit_size
  num_unit = trace_data.num_unit
  plaintexts = trace_data.plaintexts
  traces = trace_data.traces
  res_dir = trace_data.res_dir
  sbox = S_boxes[cipher_type]

  pt_unit = div(num_unit,unit_size)

  #res_key[i...byte_n] contains the best guess of the ith byte using templates
  res_key_num = []
  for byte_n in 1:num_unit
    #res[k] contains the score of the subkey k
    res = zeros(2^unit_size)
    classMeans = templates[byte_n].classMeans
    classMeanCovs = templates[byte_n].classMeanCovs
    POI = templates[byte_n].POI
    for sk in 0:2^unit_size-1
      for p in 1:length(plaintexts)
        p_byte = parse(Int,plaintexts[p][(byte_n-1)*pt_unit+1:(byte_n-1)*pt_unit+2],16)
        hw = HW[sbox[(p_byte$sk)+1]]
        μ = classMeanCovs[hw][1]
        Σ = classMeanCovs[hw][2]
        pdf_val = mvgaussian_pdf(traces[p,:][POI], μ, Σ)
        res[sk+1] += log(pdf_val)
      end
    end
    max_prob_byte = indmax(res)-1
    push!(res_key_num, max_prob_byte)
  end

  recovered_key = ""
  for v in res_key_num
    recovered_key = join([recovered_key,hex(v,2)])
  end
  return taResults(recovered_key, templates, res_dir)
end

max_stars_repo_name: udion/zypher

max_issues_repo_name: udion/zypher

max_forks_repo_name: udion/zypher


text:
###### This file is part of the MendelImpute.jl package.
###### This file contains helper functions that facilitates phasing by
###### window-by-window intersections. They are called in the `phasing.jl` file.

"""
    phase_sample!(happair1, happair2, compressed_Hunique, [seen])

Phases `happair1` and `happair2` window-by-window using a heuristic strategy.
`happair1` and `happair1` index off the complete haplotype panel and can be
many other haplotypes in window `w`. We want to pick just 1 in this pool such
that overall windows the number of haplotype switches is minimzed.

# Arguments
- `happair1`: `happair1[w]` is the haplotype index for strand1 in window `w`
- `happair2`: `happair2[w]` is the haplotype index for strand2 in window `w`
- `compressed_Hunique`: A `CompressedHaplotypes` object

# Optional storage argument
- `seen`: Preallocated storage container
"""
function phase_sample!(
    happair1::AbstractVector{<:Integer},
    happair2::AbstractVector{<:Integer},
    compressed_Hunique::CompressedHaplotypes,
    # preallocated items
    survivors1::AbstractVector{<:Integer}=Int32[],
    survivors2::AbstractVector{<:Integer}=Int32[],
    )

    windows = length(happair1)
    windows == length(happair2) || error("happair1 and happair2" *
        " have different lengths.")
    haplotypes = nhaplotypes(compressed_Hunique)
    lifespan1 = lifespan2 = 1 # counter to track survival time

    # get first window's optimal haplotypes
    h1 = happair1[1]
    h2 = happair2[1]
    store!(survivors1, get(compressed_Hunique.CW_typed[1].hapmap, h1, h1))
    store!(survivors2, get(compressed_Hunique.CW_typed[1].hapmap, h2, h2))

    @inbounds for w in 2:windows
        # get current window's best haplotypes
        h1 = happair1[w]
        h2 = happair2[w]
        h1set = get(compressed_Hunique.CW_typed[w].hapmap, h1, h1)
        h2set = get(compressed_Hunique.CW_typed[w].hapmap, h2, h2)

        # heuristic to decide whether cross-over is better
        # A   B      A   B
        # |   |  or    X
        # C   D      C   D
        AC = intersect_size_sorted(survivors1, h1set)
        BD = intersect_size_sorted(survivors2, h2set)
        AD = intersect_size_sorted(survivors1, h2set)
        BC = intersect_size_sorted(survivors2, h1set)
        crossed = false
        if AC + BD < AD + BC
            crossed = true
        end

        # Prune survivors. If there are none, record last survivor into
        # happair for all previous windows and reset survivors to current
        # window's haplotypes
        if crossed
            if AD == 0 # no survivors
                happair1[(w - lifespan1):(w - 1)] .= survivors1[1]
                store!(survivors1, h2set)
                lifespan1 = 1
            else
                intersect_lange!(survivors1, h2set)
                lifespan1 += 1
            end
            if BC == 0 # no survivors
                happair2[(w - lifespan2):(w - 1)] .= survivors2[1]
                store!(survivors2, h1set)
                lifespan2 = 1
            else
                intersect_lange!(survivors2, h1set)
                lifespan2 += 1
            end
        else
            if AC == 0 # no survivors
                happair1[(w - lifespan1):(w - 1)] .= survivors1[1]
                store!(survivors1, h1set)
                lifespan1 = 1
            else
                intersect_lange!(survivors1, h1set)
                lifespan1 += 1
            end
            if BD == 0 # no survivors
                happair2[(w - lifespan2):(w - 1)] .= survivors2[1]
                store!(survivors2, h2set)
                lifespan2 = 1
            else
                intersect_lange!(survivors2, h2set)
                lifespan2 += 1
            end
        end
    end

    # treat last few windows separately since intersection may not become empty
    happair1[(windows - lifespan1 + 1):windows] .= survivors1[1]
    happair2[(windows - lifespan2 + 1):windows] .= survivors2[1]

    return nothing
end

"""
    store!(v, u)

Deletes everything in `v` and saves each element of `u` to `v`.
"""
function store!(v::AbstractVector, u)
    empty!(v)
    for i in u
        push!(v, i)
    end
    return nothing
end

"""
    intersect_size_sorted(v::AbstractVector, u::AbstractVector)

Computes the size of `v ∩ u` in place. `v` and `u` is assumed sorted.
Repeated elements is allowed.

# Arguments
- `v`: An integer vector
- `u`: An integer vector
"""
function intersect_size_sorted(
    v::AbstractVector{<:Integer},
    u::AbstractVector{<:Integer}
    )
    lv = length(v)
    lu = length(u)
    s = 0
    i = 1
    j = 1
    @inbounds while i ≤ lv && j ≤ lu
        if v[i] == u[j]
            s += 1
            i += 1
            j += 1
        elseif v[i] > u[j]
            j += 1
        else
            i += 1
        end
    end
    return s
end
intersect_size_sorted(v::AbstractVector, u::Integer) = u in v

"""
    intersect_lange!(v::AbstractVector, u::AbstractVector)

Computes `v ∩ u` in place and stores result in `v`. `v` and `u` is assumed
sorted. Repeated elements is allowed.
"""
function intersect_lange!(v::Vector{T}, u::Vector{T}) where T <: Integer
    lv = length(v)
    lu = length(u)
    s = 0
    i = 1
    j = 1
    @inbounds while i ≤ lv && j ≤ lu
        if v[i] == u[j]
            s += 1
            v[s] = v[i]
            i += 1
            j += 1
        elseif v[i] > u[j]
            j += 1
        else
            i += 1
        end
    end
    deleteat!(v, s+1:lv)
    return nothing
end

function intersect_lange!(
    v::AbstractVector{<:Integer},
    u::Integer
    )
    keep = u in v
    empty!(v)
    keep && push!(v, u)
    nothing
end

max_stars_repo_name: Hua-Zhou/MendelImpute


text:
using VDPTag2
using Test
using MCTS
using POMDPs
using ParticleFilters
using ProgressMeter
using LinearAlgebra
using Random
using POMDPModelTools
using POMDPPolicies
using POMDPModels
using POMDPSimulators
using POMCPOW

Random.seed!(1)
pomdp = VDPTagPOMDP()
gen = NextMLFirst(MersenneTwister(31))
global s = TagState([1.0, 1.0], [-1.0, -1.0])

cpomdp = VDPTagPOMDP(mdp=VDPTagMDP(barriers=CardinalBarriers(0.2, 1.8)))
adpomdp = ADiscreteVDPTagPOMDP(cpomdp=cpomdp)
dpomdp = AODiscreteVDPTagPOMDP(pomdp, 10, 5.0)

@testset "POMCPOW Test" begin
    hr = HistoryRecorder(max_steps=100)
    solver = POMCPOWSolver(
        estimate_value=FORollout(RandomSolver()),
        tree_queries=100000,
        max_time=1.0,
        criterion=MaxUCB(110.0),
        final_criterion=MaxQ(),
        max_depth=90,
        k_action=30.0,
        alpha_action=1/30,
        k_observation=3.0,
        alpha_observation=1/200,
        next_action=NextMLFirst(Random.GLOBAL_RNG),
        rng=Random.GLOBAL_RNG,
        check_repeat_obs=false,
        check_repeat_act=false,
    )
    # Discrete
    planner = solve(solver, dpomdp)
    beliefupdater =BootstrapFilter(dpomdp, 10000)
    hist = simulate(hr, dpomdp, planner, beliefupdater)
    @show r = discounted_reward(hist)
    @test r != NaN
    # Continuous
    planner = solve(solver, cpomdp)
    beliefupdater =BootstrapFilter(cpomdp, 10000)
    hist = simulate(hr, cpomdp, planner, beliefupdater)
    @show r = discounted_reward(hist)
    @test r != NaN
end

@testset "Heuristics action type test." begin
    @test typeof(action(ToNextML(mdp(adpomdp)), s)) == Float64
    @test typeof(action(ToNextML(adpomdp), s)) == TagAction
end

struct MyNode end
MCTS.n_children(::MyNode) = rand(1:10)

@inferred next_action(gen, mdp(pomdp), s, MyNode())
@inferred next_action(gen, pomdp, initialstate(pomdp), MyNode())

for a in range(0.0, stop=2*pi, length=100)
    local s = TagState(Vec2(0,0), Vec2(1,1))
    barriers = CardinalBarriers(0.2, 1.8)
    agent_speed = 1.0
    step_size = 0.5
    delta = agent_speed*step_size*Vec2(cos(a), sin(a))
    agent = barrier_stop(barriers, s.agent, delta)
    @test agent == s.agent+delta
end

pomdp = VDPTagPOMDP()
for sao in stepthrough(pomdp, RandomPolicy(pomdp), "s,a,o", max_steps=10)
    @show sao
end

dpomdp = AODiscreteVDPTagPOMDP(pomdp, 30, 0.5)
for sao in stepthrough(dpomdp, RandomPolicy(dpomdp), "s,a,o", max_steps=10)
    @show sao
    # to address #7
    rand(Random.GLOBAL_RNG, POMDPs.observation(dpomdp, sao[1], sao[2], sao[1]))
end

pomdp = VDPTagPOMDP(mdp=VDPTagMDP(barriers=CardinalBarriers(0.2, 1.8)))
filter = BootstrapFilter(pomdp, 1000)
for sao in stepthrough(pomdp, ToNextML(pomdp), filter, "s,a,o", max_steps=10)
    @show sao
end

# test to make sure it can't pass through any walls
pomdp = VDPTagPOMDP(mdp=VDPTagMDP(barriers=CardinalBarriers(0.0, 100.0)))
filter = BootstrapFilter(pomdp, 1000)
for quadrant in [Vec2(1,1), Vec2(-1,1), Vec2(1,-1), Vec2(-1,-1)]
    @showprogress for i in 1:100
        is = rand(Random.GLOBAL_RNG, initialstate(pomdp))
        is = TagState(quadrant, is.target)
        for (s, sp) in stepthrough(pomdp, ToNextML(pomdp), filter, initialstate(pomdp), is, "s,sp", max_steps=100)
            @test all(s.agent.*quadrant .>= 0.0)
            if s == sp
                println("did not move (this should not happen a bunch of times)")
            end
        end
        policy = RandomPolicy(pomdp)
        for (s, sp) in stepthrough(pomdp, policy, updater(policy), nothing, is, "s,sp", max_steps=100)
            @test all(s.agent.*quadrant .>= 0.0)
            if s == sp
                println("did not move (this should not happen a bunch of times)")
            end
        end
    end
end

# make sure it does pass through walls sometimes without any barriers
pomdp = VDPTagPOMDP()
N = 100
for quadrant in [Vec2(1,1), Vec2(-1,1), Vec2(1,-1), Vec2(-1,-1)]
    in_other = falses(N)
    @showprogress for i in 1:N
        is = rand(Random.GLOBAL_RNG, initialstate(pomdp))
        is = TagState(quadrant, is.target)
        hr = HistoryRecorder(max_steps=100)
        hist = simulate(hr, pomdp, ToNextML(pomdp), filter)
        in_other[i] = any(any(s.agent.*quadrant .< 0.0) for s in state_hist(hist))
    end
    @show sum(in_other)/length(in_other)
    println("Should be near one?")
    @test any(in_other)
end
for quadrant in [Vec2(1,1), Vec2(-1,1), Vec2(1,-1), Vec2(-1,-1)]
    in_other = falses(N)
    @showprogress for i in 1:N
        is = rand(Random.GLOBAL_RNG, initialstate(pomdp))
        is = TagState(quadrant, is.target)
        hr = HistoryRecorder(max_steps=2)
        hist = simulate(hr, pomdp, RandomPolicy(pomdp))
        in_other[i] = any(any(s.agent.*quadrant .< 0.0) for s in state_hist(hist))
    end
    @show sum(in_other)/length(in_other)
    println("Should be near 3/4?")
    @test any(in_other)
end


text:
import CLIMAParameters as CP
import RRTMGP.Parameters as RP

#=
include(joinpath(pkgdir(RRTMGP), "parameters", "create_parameters.jl"))
param_set = create_insolation_parameters(FT)
=#
function create_insolation_parameters(FT, overrides::NamedTuple = NamedTuple())
    toml_dict = CP.create_toml_dict(FT; dict_type = "alias")
    aliases = string.(fieldnames(RP.RRTMGPParameters))
    pairs = CP.get_parameter_values!(toml_dict, aliases, "RRTMGP")
    params = merge((; pairs...), overrides) # overrides
    param_set = RP.RRTMGPParameters{FT}(; params...)
    return param_set
end

max_stars_repo_name: climate-machine/JRRTMGP

max_issues_repo_name: climate-machine/JRRTMGP

max_forks_repo_name: climate-machine/JRRTMGP


text:
using DataFrames

export one_hot_encoder

"""
_given a dataframe and a column, return the dataframe with the unique values of the column added as new columns to the dataframe_
#### parameters:
    data : DataFrame
        dataframe to apply one hot encoder
    column : Symbol or String
        column to one hot encode
    prefix : String
        the name to prefix to the new columns, default to the name of the original column
    remove_original: Boolean
        whether or not to remove the original column from the dataframe
#### returns: DataFrame
    a new dataframe with one hot encoded columns added
"""
function one_hot_encoder(data::DataFrame, column; prefix=missing, remove_original=true)
    tmp = copy(data)
    column = Symbol(column)
    prefix = ismissing(prefix) ? string(column) : prefix
    possibilities = unique(tmp[!, column])
    col_names = prefix .* "_" .* string.(possibilities)
    for i in 1:length(possibilities)
        tmp[!, Symbol(col_names[i])] .= 0
        tmp[tmp[!, column] .== possibilities[i], Symbol(col_names[i])] .= 1
    end
    if remove_original select!(tmp, Not(column)) end
    return tmp
end


text:
#addprocs(4)
using JaynesCummings, PlotlyJS
JaynesCummings.load_plotfuns()

N = 16  # N should be larger than the highest resonator initial state
ω_q = 2π*6.57E9
ω_r = ω_q
g = 0.5*2π*19E6
t_f = 200E-9 # final evolution time in sec
t_samples = 200 # number of data points in time evolution
w_samples = 50 # number of points in each dimension of Wigner matrix
#ket = ["|g,1>","1im|g,3>"]
ket = ["|g,0>","3.27|g,3>","2.39|g,6>","0.85|g,9>"] # zombie cat state
#ket = ["|g,0>","0.70710678+0.90710678im|g,3>","1.60254-0.5im|g,4>","-1.6|g,6>"]

# Generate the initial state and the Hamiltonian
ρ = gen_initialstate(N,ket)
H = gen_hamiltonian(N,ω_q,ω_r,g)

# Calculate the time evolution operators for each time
times, U_array = gen_timeevoarray(H,t_f,t_samples)

# Calculate the excited probability of the qubit as a function of time
e_prob = calc_qubittimeevo(ρ,U_array)
plot_qubittimeevo(times,e_prob)

# Calculate the photon numbers in the resonator
photons = calc_photonnumbers(N,g,times,e_prob)
plot_photonnumbers(photons)

# Calculate the density matrix by repeating the above while displacing the resonator
ρ_resonator =  calc_densitymatrix_resonator(ρ,g,times,U_array)
plot_densitymatrix(ρ_resonator)

# Calculate the Wigner function (this will take some time)
W = calc_wignerfunction_resonator(ρ,g,times,U_array,w_samples)
plot_wignerfunction(W)


text:
using PhaseSlopeIndex
using Test
using Statistics: mean

@time @testset "PhaseSlopeIndex.jl" begin
    # tests of data2psi ###############################################
    # two random signals
    signal = [[randn(1000000);] [randn(1000000);]]
    psi, _ = data2psi(signal, 100, subave=true, segave=true)
    @test all(isapprox(psi, zeros(2, 2); atol=0.01))

    # induced causality should always be inferred independent of optional arguments
    for eplen_ in [0, 200], detrend_ in [true, false]
        for segave_ in [true, false], subave_ in [true, false]
            for method_ in ["jackknife", "bootstrap"]
                # channel 2 being exactly channel 1
                ch1_ = randn(1000)
                signal = [[ch1_;] [ch1_;]]
                psi, _ = data2psi(signal, 100, method=method_,
                    eplen=eplen_, detrend=detrend_, segave=segave_, subave=subave_)
                @test psi == zeros(2, 2)

                # channel 1 leading channel 2
                ch1_ = randn(100000)
                signal = [[ch1_[2:end];] [ch1_[1:(end - 1)];]]
                psi, _ = data2psi(signal, 100, method=method_,
                    eplen=eplen_, detrend=detrend_, segave=segave_, subave=subave_)
                @test psi[1, 1] == 0.0 && psi[2, 2] == 0.0
                @test psi[1, 2] > 1.0 && psi[2, 1] + psi[1, 2] == 0.0

                # channel 2 leading channel 1
                ch2_ = randn(100000)
                signal = [[ch2_[1:(end - 1)];] [ch2_[2:end];]]
                psi, _ = data2psi(signal, 100, method=method_,
                    eplen=eplen_, detrend=detrend_, segave=segave_, subave=subave_)
                @test psi[1, 1] == 0.0 && psi[2, 2] == 0.0
                @test psi[1, 2] < -1.0 && psi[2, 1] + psi[1, 2] == 0.0
            end
        end
    end

    # test freqlist ###################################################
    ch1_ = randn(100000)
    signal = [[ch1_[2:end];] [ch1_[1:(end - 1)];]]
    psi_range, _ = data2psi(signal, 100, freqlist=1:49)
    psi_default, _ = data2psi(signal, 100)  # default is based on seglen
    @test all(psi_range == psi_default)

    # auxiliary function for testing
    """
        sin_sum = v_sin(X, A, Ω, Φ)

    returns an array summed over series of Sin waves with A-amplitudes, Ω-frequencies and Φ-phase delay
    """
    v_sin(X, A, Ω, Φ) = sum([a .* sin.(ω .* X .+ ϕ) for (a, ω, ϕ) in zip(A, Ω, Φ)])

    nsamples_ = 200000
    fs_ = 100  # sampling frequency (Hz)
    f_range_ = 12:30  # beta frequency band [12-30] Hz
    Ω_ = f_range_ * 2pi
    X_ = range(0., nsamples_ / fs_, length = nsamples_)
    A_ = rand(.1:.01:1., size(f_range_, 1))
    Φ_ = rand(size(f_range_, 1))
    Y_ = v_sin(X_, A_, Ω_, Φ_)
    Y_ch1, Y_ch2 = Y_[5:end], Y_[1:end-4]
    noise_ = randn(size(Y_ch1))
    signal = [Y_ch1 Y_ch2] .+ noise_
    seglen_ = fs_
    psi_full_range, _ = data2psi(signal, seglen_)
    psi_low, _ = data2psi(signal, seglen_, freqlist = 1:11)
    psi_freq, _ = data2psi(signal, seglen_, freqlist = f_range_)
    psi_high, _ = data2psi(signal, seglen_, freqlist = 31:49)
    @test psi_full_range[1, 2] > 1.  # ch1 leading ch2
    @test psi_freq[1, 2] > psi_full_range[1, 2]  # PSI higher in target frequency range
    @test psi_freq[1, 2] > psi_low[1, 2]  # PSI higher in target frequency range
    @test psi_freq[1, 2] > psi_high[1, 2]  # PSI higher in target frequency range

    # test data2para ##################################################
    # ndims(data) should be 2
    signal = rand(100, 3, 2)
    @test_throws DimensionMismatch data2psi(signal, 100)

    # squeeze
    signal = rand(100, 1, 4)
    @test data2psi(signal, 100)[1] == data2psi(signal[:, 1, :], 100)[1]

    # size(data, 1) >!> seglen
    signal = rand(50, 3)
    @test_throws DimensionMismatch data2psi(signal, 100)

    # if eplen == 0 then no std estimation
    signal = [[randn(1000000);] [randn(1000000);]]
    _, psi_std = data2psi(signal, 100)
    @test all(isnan.(psi_std))

    # if continuous data, subave shall not be true otherwise NaN is returned
    signal = [[randn(1000000);] [randn(1000000);]]
    psi, _ = data2psi(signal, 100, subave=true)
    @test !all(isnan.(psi))

    # tests of int ####################################################
    @test PhaseSlopeIndex.int(3.14) == 3
    @test PhaseSlopeIndex.int(-2.72) == -2
    @test PhaseSlopeIndex.int(0.0) == 0

    # tests of dropmean ###############################################
    test_a = rand(13)
    @test PhaseSlopeIndex.dropmean(test_a, 1) == mean(test_a; dims=1)
    test_a = rand(3, 5)
    @test PhaseSlopeIndex.dropmean(test_a, 1) == mean(test_a; dims=1)[1, :]
    test_a = rand(7, 11)
    @test PhaseSlopeIndex.dropmean(test_a, 2) == mean(test_a; dims=2)[:, 1]

    # tests of squeeze ################################################
    test_a = rand(17)
    @test PhaseSlopeIndex.squeeze(test_a) == test_a
    test_a = rand(3, 5, 7)
    @test PhaseSlopeIndex.squeeze(test_a) == test_a
    test_a = rand(3, 1, 7)
    @test PhaseSlopeIndex.squeeze(test_a) == test_a[:, 1, :]

    # tests of detrend ################################################
    # testing the inplace
    test_x = [0:0.01:2.0;]
    zeros_x = zeros(size(test_x))
    @test all(isapprox.(PhaseSlopeIndex.detrend!(test_x, 1), zeros_x, atol=0.001))
    @test all(isapprox.(test_x, zeros_x, atol=0.001))

    test_x = [0:0.01:2.0;]
    @test all(isapprox.(PhaseSlopeIndex.detrend!(test_x, 0), [-1.0:0.01:1.0;], atol=0.001))
    test_x = [0:0.01:2.0;]
    @test all(isapprox.(PhaseSlopeIndex.detrend!(test_x, 1), zeros_x, atol=0.001))

    test_x = [0:0.01:(2pi);]
    test_y = sin.(test_x)
    @test all(isapprox.(PhaseSlopeIndex.detrend!(test_y, 0), test_y, atol=0.001))
    test_x = [0:0.01:(2pi);]
    test_y = sin.(test_x)
    @test all(isapprox.(PhaseSlopeIndex.detrend!(test_y, 1), test_y, atol=0.001))

    test_x = [0:0.01:(2pi);]
    test_y = cos.(test_x) .+ test_x
    aux_x = [0:0.01:(2pi);] .- pi
    aux_y = cos.(test_x) .+ aux_x
    @test all(isapprox.(PhaseSlopeIndex.detrend!(test_y, 0), aux_y, atol=0.001))
    test_x = [0:0.01:(2pi);]
    test_y = cos.(test_x) .+ test_x
    @test all(isapprox.(PhaseSlopeIndex.detrend!(test_y, 1), cos.(test_x), atol=0.01))

    # tests of hanning_fun ############################################
    # the arrays are from MATLAB `hanning` function.
    hann_12 = [
        0.057271987173395
        0.215967626634422
        0.439731659872338
        0.677302443521268
        0.874255374085551
        0.985470908713026
        0.985470908713026
        0.874255374085551
        0.677302443521268
        0.439731659872338
        0.215967626634422
        0.057271987173395
    ]
    hann_13 = [
        0.049515566048790
        0.188255099070633
        0.388739533021843
        0.611260466978157
        0.811744900929367
        0.950484433951210
        1.000000000000000
        0.950484433951210
        0.811744900929367
        0.611260466978157
        0.388739533021843
        0.188255099070633
        0.049515566048790
    ]
    hanning_window_12 = PhaseSlopeIndex.hanning_fun(12)
    hanning_window_13 = PhaseSlopeIndex.hanning_fun(13)
    @test all(hanning_window_12 .≈ hann_12)
    @test all(hanning_window_13 .≈ hann_13)
    @test hanning_window_12[1] == hanning_window_12[end]
    @test hanning_window_13[1] == hanning_window_13[end]
    @test hanning_window_12[1] > 0.0
    @test hanning_window_13[1] > 0.0
end


text:
function handler(session, request)
    global test_observable
    test_observable = Observable(Dict{String, Any}())
    checkbox1 = JSServe.Checkbox(true)
    checkbox2 = JSServe.Checkbox(false)

    on(checkbox1) do value
        test_observable[] = Dict{String, Any}("checkbox1" => value)
    end
    on(checkbox2) do value
        test_observable[] = Dict{String, Any}("checkbox2" => value)
    end

    return DOM.div(checkbox1, checkbox2)
end

testsession(handler, port=8558) do app
    checkbox1_jl = children(app.dom)[1]
    checkbox2_jl = children(app.dom)[2]
    @test evaljs(app, js"document.querySelectorAll('input[type=\"checkbox\"]').length") == 2
    checkbox1 = js"document.querySelectorAll('input[type=\"checkbox\"]')[0]"
    checkbox2 = js"document.querySelectorAll('input[type=\"checkbox\"]')[1]"
    @test evaljs(app, js"$(checkbox1).checked") == true
    @test evaljs(app, js"$(checkbox2).checked") == false
    val = test_value(app, js"$(checkbox1).click()")
    @test val["checkbox1"] == false
    @test checkbox1_jl[] == false
    val = test_value(app, js"$(checkbox2).click()")
    @test val["checkbox2"] == true
    @test checkbox2_jl[] == true
    # button press from Julia
    checkbox1_jl[] = true # this won't trigger onchange!!!!
    @test evaljs(app, js"$(checkbox1).checked") == true
    checkbox2_jl[] = false
    @test evaljs(app, js"$(checkbox2).checked") == false
end


text:
module FourierTransforms

using AbstractFFTs

import AbstractFFTs: Plan, ScaledPlan,
                     fft, ifft, bfft, fft!, ifft!, bfft!,
                     plan_fft, plan_ifft, plan_bfft, plan_fft!, plan_ifft!, plan_bfft!,
                     rfft, irfft, brfft, plan_rfft, plan_irfft, plan_brfft,
                     fftshift, ifftshift,
                     rfft_output_size, brfft_output_size,
                     plan_inv, normalization

import Base: show, summary, size, ndims, length, eltype,
             *, inv, \

import LinearAlgebra: mul!

##############################################################################
# Native Julia FFTs:

include("ctfft.jl")

end


text:
using JuCSBi

atoms = JuCSBi.parse_PDB("C:/Users/Serdar/Documents/GitHub/JuCSBi.jl/test/PDB/1ggg.pdb.gz")
pdb = JuCSBi.create_structure(atoms)
print(JuCSBi.get_all_atoms(JuCSBi.get_chain(pdb, "A"))[1])


text:
__precompile__()

module GRU

include("GLHelper.jl")
include("Shapes.jl")
include("Math2D.jl")
include("Math3D.jl")
include("FTFont.jl")

using LinearAlgebra

using ModernGL
using .GLHelper

import .Shapes

abstract type AbstractRenderer end
abstract type Renderable end
abstract type Resource end

function init_resource(r::Resource, renderer::AbstractRenderer, id::Symbol)
	r.renderer = renderer
	r.id = id
	add_renderer_resource(r)
	finalizer(done, r)
end

getid(r::Resource) = r.id
getrenderer(r::Resource) = r.renderer

abstract type AbstractMesh <: Resource end
abstract type AbstractTexture <: Resource end

include("GRTypes.jl")

include("GRCamera.jl")
include("GRState.jl")
include("GRRenderer.jl")
include("GRUniforms.jl")
include("GRShader.jl")
include("GRMesh.jl")
include("GRTexture.jl")
include("GRMaterial.jl")
include("GRModel.jl")
include("GRFont.jl")

export Vec2, Vec3, Vec4, Matrix4
export VertexLayout, Mesh, Shader, Texture, Material, Model
export isvalid, init, done, set_shader_preamble, apply, setuniform, getuniform, render, position_func

end


text:
# ------------ Region alias and APIs ------------ #

const Region = MLIR.API.MlirRegion

create_region() = MLIR.API.mlirRegionCreate()
destroy!(r::Region) = MLIR.API.mlirRegionDestroy(r)
is_null(r::Region) = MLIR.API.mlirRegionIsNull(r)
get_first_block(r::Region) = MLIR.API.mlirRegionGetFirstBlock(r)
insert!(r::Region, pos::Int, b::Block) = MLIR.API.mlirRegionAppendOwnedBlock(r, pos, b)
insertafter!(r::Region, ref::Block, b::Block) = MLIR.API.mlirRegionInsertOwnedBlockAfter(r, ref, b)
insertbefore!(r::Region, ref::Block, b::Block) = MLIR.API.mlirRegionInsertOwnedBlockBefore(r, ref, b)

# Constructor.
Region() = create_region()

@doc(
"""
const Region = MLIR.API.MlirRegion
""", Region)


text:
using Nemo

if VERSION < v"0.7.0-DEV.2004"
	using Base.Test
else
	using Test
end

include("../test/ModularForms-test.jl")

test_all()


text:
include("dynamic_smagorinsky.jl")

function realspace_LES_calculation!(s::AbstractSimulation)
    hasles(s) || return nothing
    if is_dynamic_les(s)
        realspace_dynamic_les_calculation!(s)
        fourierspace_dynamic_les_calculation!(s)
        coefficient_les_calculation!(s)
    end
    calc_les!(s)
end

@par function calc_les!(s::A) where {A<:@par(AbstractSimulation)}
    @mthreads for j in TRANGE
        calc_les!(s,j)
    end
    return nothing
end

@par function calc_les!(s::A,j::Integer) where {A<:@par(AbstractSimulation)}

    rhs = s.rhs.rr
    τ = s.lesmodel.tau.rr
    Δ² = s.lesmodel.Δ²

    if haspassivescalarles(A)
        φ = s.passivescalar.φ.field.data
        fφ = s.passivescalar.flux.rr
        has_les_scalar_vorticity_model(A) && (cφ = s.passivescalar.lesmodel.c*Δ²)
    end

    if hasdensityles(A)
        ρ = s.densitystratification.ρ.field.data
        f = s.densitystratification.flux.rr
        has_les_density_vorticity_model(A) && (cρ = s.densitystratification.lesmodel.c*Δ²)
    end

    if is_dynamic_les(A)
        ca = s.lesmodel.c.rr
        û = s.lesmodel.û.rr
        La = s.lesmodel.L.rr
        Ma = s.lesmodel.M.rr
        Sa = s.lesmodel.S.rr
        Δ̂² = s.lesmodel.Δ̂²
        cmin = s.lesmodel.cmin
    elseif is_Smagorinsky(A) || is_Vreman(A) || is_production_model(A)
        c = s.lesmodel.c
        α = c*c*Δ²
    end

    if is_dynP_les(A)
        Pa = s.lesmodel.P.rr
        ŵ = s.lesmodel.ŵ.rr
        cpa = s.lesmodel.cp.rr
    end

    if is_SandP(A)
        β = s.lesmodel.cb*Δ²
    end

    if is_Silvis(A)
        c = s.lesmodel.c
        cp = s.lesmodel.cp
    end

    if is_stable_nl_model(A)
        c = s.lesmodel.c
    end

    @inbounds @msimd for i in REAL_RANGES[j]
        w = rhs[i]
        S = τ[i]

        if is_dynamic_les(A)
            νt = max(cmin,min(0.4,ca[i]))*Δ²*norm(S)
            t = 2*νt*S
        elseif is_Smagorinsky(A)
            νt = α*norm(S)
            t = 2*νt*S
        elseif is_Vreman(A)
            νt = Vreman_eddy_viscosity(S,w,c,Δ²)
            t = 2*νt*S
        elseif is_production_model(A)
            νt = production_eddy_viscosity(S,AntiSymTen(-0.5w),c,Δ²)
            t = 2*νt*S
        elseif is_Silvis(A)
            fv = fvs(S,w)
            νt = Silvis_eddy_viscosity(S,fv,c,Δ²)
            νp = Silvis_P_coeff(cp,fv,Δ²)
            t = 2*νt*S + νp*Lie(S,AntiSymTen(0.5*w))
        elseif is_stable_nl_model(A)
            W = AntiSymTen(-0.5w)
            T = traceless(square(W) - Lie(S,W) - square(S))
            νt, mnut = stable_nl_eddy_viscosity(S,W,c,Δ²)
            t = c*Δ²*T + 2*mnut*S
        end

        if is_SandP(A)
            P = Lie(S,AntiSymTen(0.5*w))
            if is_FakeSmagorinsky(A)
                t = β*P
            else
                t += β*P
            end
        end

        if is_dynP_les(A)
           t += max(-0.8,min(0.8,cpa[i]))*Δ²*Lie(S,AntiSymTen(0.5*w))
        end

        if !(is_FakeSmagorinsky(A) && !is_SandP(A))
            τ[i] = t
        end

        if hasdensityles(A)
            is_FakeSmagorinsky(A) && (νt = 0.0)
            ∇ρ = f[i]
            rhsden = νt*∇ρ
            if has_les_density_vorticity_model(A)
                rhsden += (cρ/2) * (w × ∇ρ)
            end
            f[i] = rhsden
        end

        if haspassivescalarles(A)
            is_FakeSmagorinsky(A) && (νt = 0.0)
            ∇φ = fφ[i]
            rhsp = νt*fφ[i]
            if has_les_scalar_vorticity_model(A)
                rhsp += (cφ/2) * (w × ∇φ)
            end
            fφ[i] = rhsp
        end

    end

    return nothing
end


text:
using TestFunctionRunner
TestFunctionRunner.@run(packages = ["SampleMultiplePackage1", "SampleMultiplePackage2"])


text:
"""
    BangBang.Extras

`BangBang` APIs that have no counterparts in `Base`.
"""
module Extras

export mergewith!!, modify!!

using ..BangBang: delete!!, setindex!!, mergewith!!

"""
    modify!!(f, dictlike, key) -> (dictlike′, ret)

Lookup and then update, insert or delete in one go.  If supported
(e.g., when `dictlike isa Dict`), it is done without re-computing the
hash.  Immutable containers like `NamedTuple` is also supported.

The callable `f` must accept a single argument of type
`Union{Some{valtype(dictlike)}, Nothing}`.  The value
`Some(dictlike[key])` is passed to `f` if `haskey(dictlike, key)`;
otherwise `nothing` is passed.

If `f` returns `nothing`, corresponding entry in the dictionary
`dictlike` is removed.  If `f` returns non-`nothing` value `x`, `key
=> something(x)` is inserted to `dictlike` (equivalent to
`dictlike[key] = something(x)` but more efficient).

`modify!!` returns a 2-tuple `(dictlike′, ret)` where `dictlike′` is
an updated version of `dictlike` (which may be identical to
`dictlike`) and `ret` is the returned value of `f`.

This API is inspired by
[`Control.Lens.At`](http://hackage.haskell.org/package/lens-4.17.1/docs/Control-Lens-At.html)
of Haskell's lens library.

# Examples
```jldoctest
julia> using BangBang.Extras: modify!!

julia> dict = Dict("a" => 1);

julia> dict′, ret = modify!!(dict, "a") do val
           Some(something(val, 0) + 1)
       end;

julia> ret
Some(2)

julia> dict === dict′
true

julia> dict
Dict{String,Int64} with 1 entry:
  "a" => 2

julia> dict = Dict();

julia> modify!!(dict, "a") do val
           Some(something(val, 0) + 1)
       end;

julia> dict
Dict{Any,Any} with 1 entry:
  "a" => 1

julia> modify!!(_ -> nothing, dict, "a");

julia> dict
Dict{Any,Any} with 0 entries
```

# Discussion

* [Add modify! function for lookup/update/insert/delete in one go by tkf · Pull Request #33758 · JuliaLang/julia](https://github.com/JuliaLang/julia/pull/33758)
* [Support adding to named tuples? · Issue #65 · jw3126/Setfield.jl](https://github.com/jw3126/Setfield.jl/issues/65)
"""
modify!!

modify!!(f, dict, key) = generic_modify!!(f, dict, key)

function generic_modify!!(f, dict, key)
    if haskey(dict, key)
        val = f(Some(dict[key]))
    else
        val = f(nothing)
    end
    if val === nothing
        dict′ = delete!!(dict, key)
    else
        dict′ = setindex!!(dict, something(val), key)
    end
    return (dict′, val)
end

function modify!!(f, h::Dict{K1}, key0::K2) where {K1, K2}
    promote_type(K1, K2) <: K1 || return generic_modify!!(f, h, key0)
    key = convert(K1, key0)
    isequal(key, key0) || return generic_modify!!(f, h, key0)

    # Ideally, to improve performance for the case that requires
    # resizing, we should use something like `ht_keyindex` while
    # keeping computed hash value and then do something like
    # `ht_keyindex2!` if `f` returns non-`nothing`.
    idx = Base.ht_keyindex2!(h, key)

    age0 = h.age
    found = false
    if idx > 0
        found = true
        @inbounds vold = h.vals[idx]
        vnew = f(Some(vold))
    else
        vnew = f(nothing)
    end
    if h.age != age0
        idx = Base.ht_keyindex2!(h, key)
    end

    if vnew === nothing
        if idx > 0
            Base._delete!(h, idx)
        end
    else
        val = something(vnew)
        if !(K1 isa String)
            if found && val === vold
                return (h, vnew)  # nothing to do (optimization)
            end
        end
        if !(typeof(val) <: eltype(h.vals))
            return setindex!!(h, val, key0), vnew
        end
        if idx > 0
            h.age += 1
            @inbounds h.keys[idx] = key
            @inbounds h.vals[idx] = val
        else
            @inbounds Base._setindex!(h, val, key, -idx)
        end
    end
    return (h, vnew)
end

end  # module

# `Extras` was called `Experimental` before
const Experimental = Extras


text:
# Create two cubes as VTK_POLYHEDRON cells.
using WriteVTK

vpoints = Float32[
    -1, -1, -1,
     1, -1, -1,
     1,  1, -1,
    -1,  1, -1,
    -1, -1,  1,
     1, -1,  1,
     1,  1,  1,
    -1,  1,  1,
]

append!(vpoints, vpoints .+ 2)  # append points for second cube
points = reshape(vpoints, 3, :)

cells = [
    # Cube 1
    VTKPolyhedron(
        1:8,           # connectivity vector
        (1, 4, 3, 2),  # face 1
        (1, 5, 8, 4),  # face 2
        (5, 6, 7, 8),  # etc...
        (6, 2, 3, 7),
        (1, 2, 6, 5),
        (3, 4, 8, 7),
    ),
    # Cube 2
    VTKPolyhedron(
        8 .+ (1:8),         # connectivity vector
        8 .+ (1, 4, 3, 2),  # face 1
        8 .+ (1, 5, 8, 4),  # face 2
        8 .+ (5, 6, 7, 8),  # etc...
        8 .+ (6, 2, 3, 7),
        8 .+ (1, 2, 6, 5),
        8 .+ (3, 4, 8, 7),
    ),
]

@time vtk_grid("polyhedron_cube", points, cells; compress = false) do vtk
    vtk["temperature"] = vec(sum(x -> abs2(x + 1), points; dims = 1))
    vtk["pressure"] = [42.0, 84.0]  # one per cube
end


text:
export embedding,
       NgLaplacian,
       ShiMalikLaplacian,
       YuShiPopout,
       PartialGroupingConstraints

using LightGraphs
using LightGraphs.LinAlg
using ArnoldiMethod
using LinearAlgebra
using Arpack
abstract type AbstractEmbedding <: EigenvectorEmbedder
end

"""
```julia
type NgLaplacian <: AbstractEmbedding
```
# Members
- `nev::Integer`. The number of eigenvectors to obtain
- `normalize::Bool`. Wether normalize the obtained eigenvectors

Given a affinity matrix `` W \\in \\mathbb{R}^{n \\times n} ``.  Ng et al defines the laplacian as `` L =  D^{-\\frac{1}{2}} W D^{-\\frac{1}{2}} `` where `` D `` is a diagonal matrix whose (i,i)-element is the sum of W's i-th row.

The embedding function solves a relaxed version of the following optimization problem:
``\\begin{array}{crclcl}
    \\displaystyle \\max_{ U \\in \\mathbb{R}^{n\\times k} \\hspace{10pt} } & \\mathrm{Tr}(U^T L  U)  &\\\\
   \\textrm{s.a.}  {U^T U}  =   I &&
\\end{array}``

U is a matrix that contains the `nev`  largest eigevectors of `` L ``.

# References
- [On Spectral Clustering: Analysis and an algorithm. Andrew Y. Ng, Michael I. Jordan, Yair Weiss](http://ai.stanford.edu/~ang/papers/nips01-spectral.pdf)
"""
struct NgLaplacian <: AbstractEmbedding
    nev::Integer
    normalize::Bool
end
NgLaplacian(nev::Integer) = NgLaplacian(nev, true)

"""
```julia
embedding(cfg::NgLaplacian, W::CombinatorialAdjacency)
```
Performs the eigendecomposition of the laplacian matrix of the weight matrix `` W `` defined according to [`NgLaplacian`](@ref)
"""
function embedding(cfg::NgLaplacian, W::CombinatorialAdjacency)
    return embedding(cfg, NormalizedAdjacency(W))
end

"""
```julia
embedding(cfg::NgLaplacian, L::NormalizedAdjacency)
```
Performs the eigendecomposition of the laplacian matrix of the weight matrix `` W `` defined according to [`NgLaplacian`](@ref)
"""
embedding(cfg::NgLaplacian, L::NormalizedAdjacency) = embedding(cfg, sparse(L))

"""
```julia
embedding(cfg::NgLaplacian, L::SparseMatrixCSC)
```
Performs the eigendecomposition of the laplacian matrix of the weight matrix `` W `` defined according to [`NgLaplacian`](@ref)
"""
function embedding(cfg::NgLaplacian, L::AbstractMatrix)
    (vals, vec) = LightGraphs.eigs(L, nev = cfg.nev + 5, which = LM(), restarts=5000)
    vec = real(vec)
    a = (.!(isapprox.(vals, 1)))
    vec  = vec[:, a]
    vec = vec[:, 1:cfg.nev]
    if cfg.normalize
        return normalize_rows(vec)
    else
        return vec
    end
end

"""
The normalized laplacian as defined in  `` D^{-\\frac{1}{2}} (D-W) D^{-\\frac{1}{2}} ``.

## References:
- Spectral Graph Theory. Fan Chung
- Normalized Cuts and Image Segmentation. Jiambo Shi and Jitendra Malik

```
type ShiMalikLaplacian <: AbstractEmbedding
```

# Members
- `nev::Integer`. The number of eigenvector to obtain.
- `normalize::Bool`. Wether normalize the obtained eigenvectors

"""
struct ShiMalikLaplacian <: AbstractEmbedding
    nev::Integer
    normalize::Bool
end
ShiMalikLaplacian(nev::Integer) = ShiMalikLaplacian(nev, true)

"""
```
struct PartialGroupingConstraints <: AbstractEmbedding
```

# Members

- `nev::Integer`. The number of eigenvector to obtain.
- `smooth::Bool`. Whether to user Smooth Constraints
- `normalize::Bool`. Whether to normalize the rows of the obtained vectors

Segmentation Given Partial Grouping Constraints
Stella X. Yu and Jianbo Shi
"""
struct PartialGroupingConstraints <: AbstractEmbedding
    nev::Integer
    smooth::Bool
    normalize::Bool
end

function PartialGroupingConstraints(nev::Integer; smooth::Bool=true, normalize::Bool=false)
    return PartialGroupingConstraints(nev, smooth, normalize)
end

"""
```
struct PGCMatrix{T,I,F} <: AbstractMatrix{T}
```

Partial grouping constraint structure. This sturct is passed to eigs to
performe the L*x computation according to (41), (42) and (43) of
""Segmentation Given Partial Grouping Constraints""
"""
struct PGCMatrix{T} <: AbstractMatrix{T}
    W::NormalizedAdjacency{T}
    At
end
import Base.size
import LinearAlgebra.issymmetric
import Base.*
import LinearAlgebra.mul!
function size(a::PGCMatrix)
    return size(a.W)
end
function issymmetric(a::PGCMatrix)
    return true
end
function mul!(dest::AbstractVector, a::PGCMatrix, x::AbstractVector)
    z = x - (a.At * x)
    y = a.W * z
    dest[:] = y - (a.At* y)
end

function restriction_matrix(nv::Integer,  restrictions::Vector{Vector{Integer}})
    I = Integer[]
    J = Integer[]
    V = Float64[]
    k = 0
    for j = 1:length(restrictions)
        U_t = restrictions[j]
        for i=1:length(U_t)-1
            elem_i = U_t[i]
            elem_j = U_t[i+1]
            k += 1
            push!(I, elem_i)
            push!(J, k)
            push!(V, 1.0)

            push!(I, elem_j)
            push!(J, k)
            push!(V, -1.0)
        end
    end
    return sparse(I, J, V, nv, k)
end
"""
```
function embedding(cfg::PartialGroupingConstraints, L::NormalizedAdjacency, restrictions::Vector{Vector{Integer}})
```

# Arguments
- `cfg::PartialGroupingConstraints`
- `L::NormalizedAdjacency`
- `restrictions::Vector{Vector{Integer}}`
"""
function embedding(cfg::PartialGroupingConstraints, L::NormalizedAdjacency, restrictions::Vector{Vector{Integer}})
    U = restriction_matrix(size(L, 1), restrictions)
    if (cfg.smooth)
        U = sparse(L.A)' * U
    end
    DU = sparse(spdiagm(0=>prescalefactor(L)) * U)
    F = svds(DU, nsv=size(U, 2)-1)[1]
    AAt = sparse(F.U)*sparse(F.U)'
    (eigvals, eigvec) = LightGraphs.eigs(PGCMatrix(L, AAt), nev = cfg.nev, which = LM())
    eigvec = real(eigvec)
    V = spdiagm(0=>prescalefactor(L)) * eigvec
    if cfg.normalize
        return normalize_rows(V)
    else
        return V
    end
end
"""
```
embedding(cfg::PartialGroupingConstraints, gr::Graph, restrictions::Vector{Vector{Integer}})
```

# Arguments
- `cfg::PartialGroupingConstraints`
- `L::NormalizedAdjacency`
-  `restrictions::Vector{Vector{Integer}}`
"""
function embedding(cfg::PartialGroupingConstraints, gr::Graph, restrictions::Vector{Vector{Integer}})
    L = NormalizedAdjacency(CombinatorialAdjacency(adjacency_matrix(gr)))
    return embedding(cfg, L, restrictions)
end

"""
```
struct YuShiPopout <: AbstractEmbedding
```

# Members

- `nev::Integer`. The number of eigenvector to obtain.

Understanding Popout through Repulsion
Stella X. Yu and Jianbo Shi
"""

struct YuShiPopout <: AbstractEmbedding
    nev::Integer
    normalize::Bool
end
YuShiPopout(nev::Integer) = YuShiPopout(nev, true)

"""

```julia
function embedding(cfg::YuShiPopout,  grA::Graph, grR::Graph)
```

# References
- Grouping with Directed Relationships. Stella X. Yu and Jianbo Shi
- Understanding Popout through Repulsion. Stella X. Yu and Jianbo Shi
"""
function embedding(cfg::YuShiPopout, grA::Graph, grR::Graph, restrictions::Vector{Vector{Integer}})
    Wa = adjacency_matrix(grA)
    Wr = adjacency_matrix(grR)
    dr = vec(sum(Wr, 1))
    W = Wa - Wr + spdiagm(dr)
    return embedding(PartialGroupingConstraints(cfg.nev), NormalizedAdjacency(CombinatorialAdjacency(W)), restrictions)
end

"""
```julia
function embedding(cfg::YuShiPopout,  grA::Graph, grR::Graph)
```

# References
- Grouping with Directed Relationships. Stella X. Yu and Jianbo Shi
- Understanding Popout through Repulsion. Stella X. Yu and Jianbo Shi
"""
function embedding(cfg::YuShiPopout, grA::Graph, grR::Graph)
    Wa = adjacency_matrix(grA)
    Wr = adjacency_matrix(grR)
    dr = vec(sum(Wr, dims=1))
    da = vec(sum(Wa, dims=1))
    Weq = Wa - Wr + spdiagm(0=>dr)
    Deq = spdiagm(0=>da + dr)
    Wa = nothing
    da = nothing
    Wr = nothing
    dr = nothing
    (eigvals, eigvec) = Arpack.eigs(Weq, Deq, nev = cfg.nev, tol = 0.000001,  which = :LM)
    #indexes = sortperm(real(eigvals))
    return eigvec
    #if (cfg.normalize)
    #    return SpectralClustering.normalize_rows(eigvec)
    #else
    #    return eigvec
    #end
end

"""
```julia
embedding(cfg::NgLaplacian, W::CombinatorialAdjacency)
```
Performs the eigendecomposition of the laplacian matrix of the weight matrix `` W `` defined according to [`NgLaplacian`](@ref)
"""
function embedding(cfg::ShiMalikLaplacian, W::CombinatorialAdjacency)
    return embedding(cfg, NormalizedLaplacian(NormalizedAdjacency(W)))
end

"""
```julia
embedding(cfg::ShiMalikLaplacian, L::NormalizedLaplacian)
```
# Parameters
-   `cfg::ShiMalikLaplacian`. An instance of a [`ShiMalikLaplacian`](@ref)  that specify the number of eigenvectors to obtain
- `gr::Union{Graph,SparseMatrixCSC}`. The `Graph`(@ref Graph) or the weight matrix of wich is going to be computed the normalized laplacian matrix.

Performs the eigendecomposition of the normalized laplacian matrix of
the laplacian matriz `L` defined acoording to [`ShiMalikLaplacian`](@ref). Returns
the cfg.nev eigenvectors associated with the non-zero smallest
eigenvalues.
"""
function embedding(cfg::ShiMalikLaplacian, L::NormalizedLaplacian)
    (vals, V) = LightGraphs.eigs(sparse(L), nev=min(cfg.nev + 10, size(L, 1)), which = SR(), restarts=5000)
    idxs = findall(real(vals) .> 0.0000001)
    idxs = idxs[1:min(length(idxs), cfg.nev)]
    V = spdiagm(0 => L.A.A.D.^(-1 / 2)) * real(V[:,idxs])
    if cfg.normalize
        return normalize_rows(V)
    else
        return V
    end
end

"""
```
embedding{T<:AbstractEmbedding}(cfg::T, neighborhood::VertexNeighborhood, oracle::Function, data)
```
"""
function embedding(cfg::T, neighborhood::VertexNeighborhood,  oracle::Function, data) where T <: AbstractEmbedding
    graph = create(cfg.graph_creator, data)
    return embedding(cfg, graph)
end

"""
```julia
embedding(cfg::T, gr::Graph) where T<:AbstractEmbedding
```
Performs the eigendecomposition of the laplacian matrix of the weight matrix `` W `` derived from the graph `gr` defined according to [`NgLaplacian`](@ref)
"""
function embedding(cfg::T, gr::Graph) where T <: AbstractEmbedding
    return embedding(cfg, CombinatorialAdjacency(adjacency_matrix(gr, dir = :both)))
end


text:
# MetNet
mets(metnet::MetNet) = metnet.mets
mets(metnet::MetNet, ider::IDER_TYPE) = metnet.mets[metindex(metnet, ider)]

metscount(metnet::MetNet) = length(metnet.mets)

rxns(metnet::MetNet) = metnet.rxns
rxns(metnet::MetNet, ider::IDER_TYPE) = metnet.rxns[rxnindex(metnet, ider)]

rxnscount(model) = length(model.rxns)

ub(metnet::MetNet) = metnet.ub
ub(metnet::MetNet, ider::IDER_TYPE) = metnet.ub[rxnindex(metnet, ider)]

lb(metnet::MetNet) = metnet.lb
lb(metnet::MetNet, ider::IDER_TYPE) = metnet.lb[rxnindex(metnet, ider)]

b(metnet::MetNet) = metnet.b
b(metnet::MetNet, ider::IDER_TYPE) = metnet.b[rxnindex(metnet, ider)]

bounds(metnet::MetNet, ider::IDER_TYPE) = (idx = rxnindex(metnet, ider); (metnet.lb[idx], metnet.ub[idx]))

S(metnet::MetNet, metider::IDER_TYPE, rxnider::IDER_TYPE) =
    metnet.S[metindex(metnet, metider), rxnindex(metnet, rxnider)]

rxn_mets(metnet::MetNet, ider::IDER_TYPE) = findall(metnet.S[:,rxnindex(metnet, ider)] .!= 0.0)
rxn_reacts(metnet::MetNet, ider::IDER_TYPE) = findall(metnet.S[:,rxnindex(metnet, ider)] .< 0.0)
rxn_prods(metnet::MetNet, ider::IDER_TYPE) = findall(metnet.S[:,rxnindex(metnet, ider)] .> 0.0)

met_rxns(metnet::MetNet, ider::IDER_TYPE) = findall(metnet.S[metindex(metnet, ider), :] .!= 0.0)

# MetState
# Interface
av(s::AbstractMetState) = error("You must implement a 'av(s::$(typeof(s)))' method")
va(s::AbstractMetState) = error("You must implement a 'va(s::$(typeof(s)))' method")

av(s::Vector{<:Real}) = s
va(s::Vector{<:Real}) = s

# Commons getter interface
for fun_name in [:av, :va]

    @eval begin
        $(fun_name)(state::AbstractMetState, idxs) =
            $(fun_name)(state)[idxs]
        $(fun_name)(metnet::MetNet, state::AbstractMetState, ider) =
            $(fun_name)(state)[rxnindex(metnet, ider)]
        $(fun_name)(metnet::MetNet, state::AbstractMetState, iders::Vector) =
            [$(fun_name)(metnet, state, ider) for ider in iders]
        $(fun_name)(metnet::MetNet, states::Vector, ider) =
            [$(fun_name)(metnet, state, ider) for state in states]
        $(fun_name)(metnets::Vector, states::Vector, ider) =
            [$(fun_name)(metnet, state, ider) for (metnet, state) in zip(metnets, states)]
    end
end