text:
adj = [0 1 0 1;
       1 0 1 0;
       0 1 0 1;
       1 0 1 0]

@testset "defaultstaterepresentation.jl" begin

    @testset "DefaultStateRepresentation structure" begin
        g = SeaPearl.CPLayerGraph()
        nodeFeatures = [1.0f0 1.0f0; 2.0f0 2.0f0]
        variableIdx = 1
        dsr = SeaPearl.DefaultStateRepresentation{SeaPearl.DefaultFeaturization,SeaPearl.DefaultTrajectoryState}(g, nodeFeatures, nothing, variableIdx, nothing)

        @test dsr.cplayergraph == g
        @test dsr.nodeFeatures == nodeFeatures
        @test dsr.variableIdx == 1

    end

    @testset "DefaultSR from CPModel" begin
        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)

        x = SeaPearl.IntVar(3, 3, "x", trailer)
        y = SeaPearl.IntVar(2, 3, "y", trailer)
        SeaPearl.addVariable!(model, x)
        SeaPearl.addVariable!(model, y)
        SeaPearl.addConstraint!(model, SeaPearl.Equal(x, y, trailer))
        SeaPearl.addConstraint!(model, SeaPearl.NotEqual(x, y, trailer))

        dsr = SeaPearl.DefaultStateRepresentation(model)
        SeaPearl.update_representation!(dsr, model, x)

        @test Matrix(LightGraphs.LinAlg.adjacency_matrix(dsr.cplayergraph)) == [0 0 1 1 0 0
                                                                                0 0 1 1 0 0
                                                                                1 1 0 0 0 1
                                                                                1 1 0 0 1 1
                                                                                0 0 0 1 0 0
                                                                                0 0 1 1 0 0]

        @test dsr.nodeFeatures == Float32[  1 1 0 0 0 0
                                        0 0 1 1 0 0
                                        0 0 0 0 1 1]
        @test dsr.variableIdx == 3
        @test SeaPearl.cpVertexFromIndex(SeaPearl.CPLayerGraph(model), dsr.variableIdx).variable == model.variables["x"]
    end

    @testset "update_representation!()" begin
        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)

        x = SeaPearl.IntVar(2, 3, "x", trailer)
        y = SeaPearl.IntVar(2, 3, "y", trailer)
        SeaPearl.addVariable!(model, x)
        SeaPearl.addVariable!(model, y)
        SeaPearl.addConstraint!(model, SeaPearl.Equal(x, y, trailer))
        SeaPearl.addConstraint!(model, SeaPearl.NotEqual(x, y, trailer))

        dsr = SeaPearl.DefaultStateRepresentation(model)
        SeaPearl.update_representation!(dsr, model, x)

        @test Matrix(LightGraphs.LinAlg.adjacency_matrix(dsr.cplayergraph)) == [0 0 1 1 0 0
                                                    0 0 1 1 0 0
                                                    1 1 0 0 1 1
                                                    1 1 0 0 1 1
                                                    0 0 1 1 0 0
                                                    0 0 1 1 0 0]

        @test dsr.nodeFeatures == Float32[   1 1 0 0 0 0
                                                        0 0 1 1 0 0
                                                        0 0 0 0 1 1]

        @test dsr.variableIdx == 3
        @test SeaPearl.cpVertexFromIndex(SeaPearl.CPLayerGraph(model), dsr.variableIdx).variable == model.variables["x"]

        SeaPearl.assign!(y, 2)
        g = SeaPearl.CPLayerGraph(model)

        SeaPearl.update_representation!(dsr, model, y)

        @test Matrix(LightGraphs.LinAlg.adjacency_matrix(dsr.cplayergraph)) == [0 0 1 1 0 0
                                                    0 0 1 1 0 0
                                                    1 1 0 0 1 1
                                                    1 1 0 0 1 0
                                                    0 0 1 1 0 0
                                                    0 0 1 0 0 0]

        @test dsr.nodeFeatures == Float32[   1 1 0 0 0 0
                                                        0 0 1 1 0 0
                                                        0 0 0 0 1 1]
        @test dsr.variableIdx == 4
        @test SeaPearl.cpVertexFromIndex(g, dsr.variableIdx).variable == model.variables["y"]

    end

    @testset "trajectoryState constructor" begin
        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)

        x = SeaPearl.IntVar(2, 3, "x", trailer)
        y = SeaPearl.IntVar(2, 3, "y", trailer)
        SeaPearl.addVariable!(model, x)
        SeaPearl.addVariable!(model, y)
        SeaPearl.addConstraint!(model, SeaPearl.Equal(x, y, trailer))
        SeaPearl.addConstraint!(model, SeaPearl.NotEqual(x, y, trailer))

        dsr = SeaPearl.DefaultStateRepresentation(model)
        SeaPearl.update_representation!(dsr, model, x) #add x as the branching variable
        dts = SeaPearl.DefaultTrajectoryState(dsr)  #creates the DefaultTrajectoryState object

        @test dts.variableIdx == SeaPearl.indexFromCpVertex(dsr.cplayergraph, SeaPearl.VariableVertex(x))

    end

    @testset "BatchedDefaultTrajectoryState constructor" begin
        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)

        x = SeaPearl.IntVar(2, 3, "x", trailer)
        y = SeaPearl.IntVar(2, 3, "y", trailer)
        SeaPearl.addVariable!(model, x)
        SeaPearl.addVariable!(model, y)
        SeaPearl.addConstraint!(model, SeaPearl.Equal(x, y, trailer))
        SeaPearl.addConstraint!(model, SeaPearl.NotEqual(x, y, trailer))

        dsr = SeaPearl.DefaultStateRepresentation(model)
        SeaPearl.update_representation!(dsr, model, x) #add x as the branching variable
        dts = SeaPearl.DefaultTrajectoryState(dsr)  #creates the DefaultTrajectoryState object

        batchedDtsSingle = dts |> cpu
        @test batchedDtsSingle.fg.graph[:,:,1]== adjacency_matrix(dts.fg)
        @test batchedDtsSingle.fg.nf[:,:,1] == dts.fg.nf
        @test batchedDtsSingle.fg.ef[:,:,1] == dts.fg.ef
        @test batchedDtsSingle.fg.gf[:,1] == dts.fg.gf
        @test batchedDtsSingle.variableIdx[1] == dts.variableIdx

        SeaPearl.assign!(x, 2)
        SeaPearl.fixPoint!(model, SeaPearl.getOnDomainChange(x))
        SeaPearl.update_representation!(dsr, model, x) #add x as the branching variable
        dts2 = SeaPearl.DefaultTrajectoryState(dsr)  #creates the DefaultTrajectoryState object

        batchedDts = [dts,dts2] |> cpu
        @test size(batchedDts.fg.graph,3) == 2
        @test size(batchedDts.fg.nf,3) == 2
        @test size(batchedDts.fg.ef,3) == 2
        @test size(batchedDts.fg.gf,2) == 2
        @test size(batchedDts.variableIdx,1) == 2
    end
end


text:
################################################################################
#
#   FUNCTIONS FOR PLOTTING OF RECIPROCAL POINTS
#   (overwriting plotSite function for AbstractSite)
#
################################################################################

# PYPLOT PLOT FUNCTION OVERWRITTEN
function PyPlot.plot(
            path :: P
            ;
            color :: Vector{<:Integer} = [250,30,0],
            site_radius :: Real = 0.5,
            kwargs...
        ) where {S<:AbstractReciprocalPoint{2}, P<:AbstractReciprocalPath{S}}

    # plot all points in the path
    for p in points(path)
        # plot the point
        plotSite(p,site_radius,color;kwargs...)
    end
    # connect the points
    PyPlot.plot([point(p)[1] for p in points(path)], [point(p)[2] for p in points(path)], color=color./255, linewidth=1, linestyle="-")
end

# PYPLOT PLOT FUNCTION OVERWRITTEN
function PyPlot.plot(
            path :: P
            ;
            color :: Vector{<:Integer} = [250,30,0],
            site_radius :: Real = 0.5,
            kwargs...
        ) where {S<:AbstractReciprocalPoint{3}, P<:AbstractReciprocalPath{S}}

    # plot all points in the path
    for p in points(path)
        # plot the point
        plotSite(p,site_radius,color;kwargs...)
    end
    # connect the points
    PyPlot.plot3D([point(p)[1] for p in points(path)], [point(p)[2] for p in points(path)], [point(p)[3] for p in points(path)], color=color./255, linewidth=1, linestyle="-")
end

# export the plot function
export plot

ext: jl


text:
using HorizonSideRobots
include("roblib.jl")

function perimetr!(r::Robot)
    vertically,horizontal=angle(r)
    for side in (Nord,Ost,Sud,West)
        while isborder(r,side)==false
            putmarkers!(r,side)
        end
    end
    moves!(r,Ost,horizontal)
    moves1!(r,Nord,vertically)
end

function moves!(r::Robot,side::HorizonSide)
    num_steps=0
    while isborder(r,side)==false
        move!(r,side)
        num_steps+=1
    end
    return num_steps
end

function moves!(r::Robot,side::HorizonSide,num_steps::Int)
    for _ in 1:num_steps
        move!(r,side)
    end
end

perimetr!(r)


text:
f(u,p,t) = 0.98u

using DifferentialEquations
f(u,p,t) = 0.98u
u0 = 1.0
tspan = (0.0,1.0)
prob = ODEProblem(f,u0,tspan)

sol = solve(prob)

using Plots; gr()
plot(sol)

plot(sol,linewidth=5,title="Solution to the linear ODE with a thick line",
     xaxis="Time (t)",yaxis="u(t) (in μm)",label="My Thick Line!") # legend=false

plot!(sol.t, t->1.0*exp(0.98t),lw=3,ls=:dash,label="True Solution!")

sol.t

sol.u

[t+u for (u,t) in tuples(sol)]

sol

sol(0.45)

sol = solve(prob,abstol=1e-8,reltol=1e-8)

plot(sol)
plot!(sol.t, t->1.0*exp(0.98t),lw=3,ls=:dash,label="True Solution!")

sol = solve(prob,saveat=0.1)

sol = solve(prob,saveat=[0.2,0.7,0.9])

sol = solve(prob,dense=false)

sol = solve(prob,save_everystep=false)

sol = solve(prob,save_everystep=false,save_start = false)

sol = solve(prob,alg_hints=[:stiff])

sol = solve(prob,Tsit5(),reltol=1e-6)

function lorenz!(du,u,p,t)
    σ,ρ,β = p
    du[1] = σ*(u[2]-u[1])
    du[2] = u[1]*(ρ-u[3]) - u[2]
    du[3] = u[1]*u[2] - β*u[3]
end

u0 = [1.0,0.0,0.0]

p = (10,28,8/3) # we could also make this an array, or any other type!

tspan = (0.0,100.0)
prob = ODEProblem(lorenz!,u0,tspan,p)

sol = solve(prob)

sol.t[10],sol[10]

sol[2,10]

A = Array(sol)

plot(sol)

plot(sol,vars=(1,2,3))

plot(sol,vars=(1,2,3),denseplot=false)

plot(sol,vars=(0,2))

function lotka_volterra!(du,u,p,t)
  du[1] = p[1]*u[1] - p[2]*u[1]*u[2]
  du[2] = -p[3]*u[2] + p[4]*u[1]*u[2]
end

using ParameterizedFunctions
lv! = @ode_def LotkaVolterra begin
  dx = a*x - b*x*y
  dy = -c*y + d*x*y
end a b c d

u0 = [1.0,1.0]
p = (1.5,1.0,3.0,1.0)
tspan = (0.0,10.0)
prob = ODEProblem(lv!,u0,tspan,p)
sol = solve(prob)
plot(sol)

lv!.Jex

A  = [1. 0  0 -5
      4 -2  4 -3
     -4  0  0  1
      5 -2  2  3]
u0 = rand(4,2)
tspan = (0.0,1.0)
f(u,p,t) = A*u
prob = ODEProblem(f,u0,tspan)
sol = solve(prob)

sol[3]

big_u0 = big.(u0)

prob = ODEProblem(f,big_u0,tspan)
sol = solve(prob)

sol[1,3]

prob = ODEProblem(f,big_u0,big.(tspan))
sol = solve(prob)

using StaticArrays
A  = @SMatrix [ 1.0  0.0 0.0 -5.0
                4.0 -2.0 4.0 -3.0
               -4.0  0.0 0.0  1.0
                5.0 -2.0 2.0  3.0]
u0 = @SMatrix rand(4,2)
tspan = (0.0,1.0)
f(u,p,t) = A*u
prob = ODEProblem(f,u0,tspan)
sol = solve(prob)

sol[3]

using DiffEqTutorials
DiffEqTutorials.tutorial_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])


text:
using DifferentialEquations, NLsolve

# Entire analysis for one pulse simulation (1 system)
function pulse_unit(par, p_length, p_strength)
    # Parameters
    u0 = [1.0, 1.5, 1.5]
    pulse_start = 200
    pulse_end = pulse_start + p_length
    pulse_event_times = union(pulse_start, pulse_end)
    t_end = 350
    t_span = (0.0, t_end)
    t_start = 175.0
    len = 100000
    t_grid = range(t_start, t_end, length = len)
    t_after = range(pulse_end, t_end, length = len)
    Ω = par.Ω

    # functions
    pulse_event(u, t, integrator) = t ∈ pulse_event_times
    # NB: for pulse p_strength is a multiplicator
    function forcing_pulse!(integrator)
        if integrator.t == pulse_start
            # multiplication for pulse
            integrator.p.K = p_strength * integrator.p.K_base
        elseif integrator.t == pulse_end
            integrator.p.K = integrator.p.K_base
        end
        return
    end
    cb_pulse = DiscreteCallback(pulse_event, forcing_pulse!)

    # run
    if par.pref == adapt_pref
        # compute ω
        # NB: This may often increases computation time as we may compute
        # results for passive omnivory twice but the code is much easier to
        # read. Here the time is not a major concern as it takes few minutes to
        # run.
        par_fix = deepcopy(par)
        par_fix.pref = fixed_pref
        prob_fix = ODEProblem(model!, u0, t_span, deepcopy(par_fix),
            tstops = pulse_event_times)
        sol_fix = solve(prob_fix, reltol = 1e-8, abstol = 1e-8,
            callback = cb_pulse)
        # sol_fix(t_grid)[1] starts at t_start
        eq = nlsolve(
            (du, u) -> model!(du, u, par_fix, 0.0),
            sol_fix(t_grid)[1]
        ).zero
        ### Solving for ω we have `ω = Ω * C^* / (Ω * C^* + (1 - Ω) * R^*)`
        ω = Ω * eq[2] / (Ω * eq[2] + (1 - Ω) * eq[1])
        par.ω = ω
    end

    par_pulse = deepcopy(par)
    par_pulse.K = p_strength * par.K_base
    prob = ODEProblem(model!, u0, t_span, deepcopy(par), tstops = pulse_event_times)
    sol = solve(prob, reltol = 1e-8, abstol = 1e-8, callback = cb_pulse)
    sol_grid = sol(t_grid)
    t_hit_eq = find_times_hit_equil_press(sol(t_after))
    # println(t_hit_eq)

    # Same equilibrium before and after
    eq_before = find_eq(sol_grid[1], par)

    # evaluate on
    λ1 = λ1_stability(cmat(eq_before, par))
    λ1_imag = λ1_stability_imag(cmat(eq_before, par))
    react = ν_stability(cmat(eq_before, par))

    predator_OS = overshoot(sol, eq_before, 3, t_hit_eq[3], t_end)
    consumer_OS = overshoot(sol, eq_before, 2, t_hit_eq[2], t_end)
    resource_OS = overshoot(sol, eq_before, 1, t_hit_eq[1], t_end)

    # Calculate max-min metric
    predator_mm = min_max(sol, 3, t_hit_eq[3], t_end)
    consumer_mm = min_max(sol, 2, t_hit_eq[2], t_end)
    resource_mm = min_max(sol, 1, t_hit_eq[1], t_end)

    push!(
        [],
        λ1,
        -1 ./ λ1,
        λ1_imag,
        react,
        predator_OS,
        consumer_OS,
        resource_OS,
        predator_mm,
        consumer_mm,
        resource_mm,
        eq_before,  # for consistency in output size with `press()`
        eq_before,
        par,
        sol
    )

end

# Perform the pulse for the 3 systems:
# - Food chain: "chain";C OP (fixed)
# - Passive omnivory: "passive";
# - Active omnivory: "active";
function pulse(base_par, Ω, p_length, p_strength)

    par_chain = deepcopy(base_par)
    par_chain.Ω = 0.0
    par_passive = deepcopy(base_par)
    par_passive.Ω = Ω
    par_passive.pref = fixed_pref
    par_active = deepcopy(par_passive)
    par_active.pref = adapt_pref

    push!(
        [],
        pulse_unit(par_chain, p_length, p_strength),
        pulse_unit(par_passive, p_length, p_strength),
        pulse_unit(par_active, p_length, p_strength)
    )

end

max_stars_repo_name: McCannLab/Labmnivory

max_issues_repo_name: McCannLab/Labmnivory

max_forks_repo_name: McCannLab/Labmnivory


text:
@testset "fit normal SnpLinAlg" begin
	#simulat data with k true predictors, from distribution d and with link l.
	n = 1000
	p = 10000
	k = 10
	d = Normal
	l = canonicallink(d())

	#construct SnpArraym, snpmatrix, and non genetic covariate (intercept)
	x = simulate_random_snparray(undef, n, p)
	xla = SnpLinAlg{Float64}(x, model=ADDITIVE_MODEL, center=true, scale=true)
	z = ones(n)

	# simulate response, true model b, and the correct non-0 positions of b
	y, true_b, correct_position = simulate_random_response(xla, k, d, l)

	# run IHT
	result = fit_iht(y, xla, z, J=1, k=k, d=d(), l=l)
	show(result)

	@test length(result.beta) == 10000
	@test count(!iszero, result.beta) == k
	@test result.c[1] != 0
	@test result.k == 10
end

@testset "fit Bernoulli SnpLinAlg" begin
	#simulat data with k true predictors, from distribution d and with link l.
	n = 1000
	p = 10000
	k = 10
	d = Bernoulli
	l = canonicallink(d())

	#set random seed
	Random.seed!(1111)

	#construct SnpArraym, snpmatrix, and non genetic covariate (intercept)
	x = simulate_random_snparray(undef, n, p)
	xla = SnpLinAlg{Float64}(x, model=ADDITIVE_MODEL, center=true, scale=true)
	z = ones(n)

	# simulate response, true model b, and the correct non-0 positions of b
	y, true_b, correct_position = simulate_random_response(xla, k, d, l)

	# run IHT
	result = fit_iht(y, xla, z, J=1, k=k, d=d(), l=l)

	@test length(result.beta) == 10000
	@test count(!iszero, result.beta) == k
	@test result.c[1] != 0
	@test result.k == 10
end

@testset "fit Poisson SnpLinAlg" begin
	#simulat data with k true predictors, from distribution d and with link l.
	n = 1000
	p = 10000
	k = 10
	d = Poisson
	l = canonicallink(d())

	#construct SnpArraym, snpmatrix, and non genetic covariate (intercept)
	x = simulate_random_snparray(undef, n, p)
	xla = SnpLinAlg{Float64}(x, model=ADDITIVE_MODEL, center=true, scale=true)
	z = ones(n)

	# simulate response, true model b, and the correct non-0 positions of b
	y, true_b, correct_position = simulate_random_response(xla, k, d, l)

	# run IHT
	result = fit_iht(y, xla, z, J=1, k=k, d=d(), l=l)

	@test length(result.beta) == 10000
	@test count(!iszero, result.beta) == k
	@test result.c[1] != 0
	@test result.k == 10
end

@testset "fit NegativeBinomial SnpLinAlg" begin
	#simulat data with k true predictors, from distribution d and with link l.
	n = 1000
	p = 10000
	k = 10
	d = NegativeBinomial
	l = LogLink()

	#construct SnpArraym, snpmatrix, and non genetic covariate (intercept)
	x = simulate_random_snparray(undef, n, p)
	xla = SnpLinAlg{Float64}(x, model=ADDITIVE_MODEL, center=true, scale=true)
	z = ones(n)

	# simulate response, true model b, and the correct non-0 positions of b
	y, true_b, correct_position = simulate_random_response(xla, k, d, l)

	#run result
	result = fit_iht(y, xla, z, J=1, k=k, d=d(), l=l)

	@test length(result.beta) == 10000
	@test count(!iszero, result.beta) == k
	@test result.c[1] != 0
	@test result.k == 10
end

@testset "fit with >1 non-genetic covariates SnpLinAlg" begin
	#simulat data with k true predictors, from distribution d and with link l.
	n = 1000
	p = 10000
	k = 10
	d = Normal
	l = canonicallink(d())

	#construct SnpArraym, snpmatrix, and non genetic covariate (intercept)
	x = simulate_random_snparray(undef, n, p)
	xla = SnpLinAlg{Float64}(x, model=ADDITIVE_MODEL, center=true, scale=true)
	z = ones(n, 2) # 1st column intercept
	z[:, 2] .= randn(n)

	#define true_b and true_c
	true_b = zeros(p)
	true_b[1:k-2] = randn(k-2)
	shuffle!(true_b)
	correct_position = findall(!iszero, true_b)
	true_c = [3.0; 3.5]

	#simulate phenotype
	prob = GLM.linkinv.(l, xla * true_b .+ z * true_c)
	y = [rand(d(i)) for i in prob]

	# run IHT
	result = fit_iht(y, xla, z, J=1, k=k, d=d(), l=l)

	@test length(result.beta) == 10000
	@test length(result.c) == 2
	@test length(findall(!iszero, result.beta)) == 10
	@test findall(!iszero, result.c) == [1;2]
	@test result.k == 10
end

@testset "model selection on non-genetic covariates" begin
	#simulat data with k true predictors, from distribution d and with link l.
	n = 1000
	p = 10000
	k = 10
	d = Normal
	l = canonicallink(d())

	#construct SnpArraym, snpmatrix, and non genetic covariate (intercept)
	x = simulate_random_snparray(undef, n, p)
	xla = SnpLinAlg{Float64}(x, model=ADDITIVE_MODEL, center=true, scale=true)
	z = ones(n, 2) # 1st column intercept
	z[:, 2] .= randn(n)

	# define true_b and true_c
	true_b = zeros(p)
	true_b[1:k-2] = randn(k-2)
	shuffle!(true_b)
	correct_position = findall(!iszero, true_b)
	true_c = [3.0; 0]

	# simulate phenotype
	prob = GLM.linkinv.(l, xla * true_b .+ z * true_c)
	y = [rand(d(i)) for i in prob]

	#run result, keeping only intercept
	zkeep = convert(BitVector, [true, false])
	result = fit_iht(y, xla, z, J=1, k=k, d=d(), l=l, zkeep = zkeep)

	@test length(result.beta) == 10000
	@test length(result.c) == 2
	@test length(findall(!iszero, result.beta)) == 10
	@test count(!iszero, result.c) == 1 # the 2nd covariate should be excluded
	@test result.k == 10
end

@testset "Correlated predictors and double sparsity Float64 matrix" begin
	#simulat data with k true predictors, from distribution d and with link l.
    n = 1000
    p = 10000
    d = Normal
	l = canonicallink(d())
    block_size = 20
    num_blocks = Int(p / block_size)

    #set random seed
    Random.seed!(1111)

    # assign group membership
    membership = collect(1:num_blocks)
    g = zeros(Int64, p)
    for i in 1:length(membership)
        for j in 1:block_size
            cur_row = block_size * (i - 1) + j
            g[block_size*(i - 1) + j] = membership[i]
        end
    end

    #simulate correlated snparray
    x = simulate_correlated_snparray(undef, n, p)
    z = ones(n) # the intercept
    x_float = convert(Matrix{Float64}, x, model=ADDITIVE_MODEL, center=true, scale=true)

    #simulate true model, where 5 groups each with 3 snps contribute
    true_b = zeros(p)
    true_groups = randperm(num_blocks)[1:5]
    within_group = [randperm(block_size)[1:3] randperm(block_size)[1:3] randperm(block_size)[1:3] randperm(block_size)[1:3] randperm(block_size)[1:3]]
    correct_position = zeros(Int64, 15)
    for i in 1:5
        cur_group = block_size * (true_groups[i] - 1)
        cur_group_snps = cur_group .+ within_group[:, i]
        correct_position[3*(i-1)+1:3i] .= cur_group_snps
    end
    for i in 1:15
        true_b[correct_position[i]] = rand(-1:2:1) * 0.2
    end
    sort!(correct_position)

    # simulate phenotype
    if d == Normal || d == Bernoulli || d == Poisson
        prob = GLM.linkinv.(l, x_float * true_b)
        clamp!(prob, -20, 20)
        y = [rand(d(i)) for i in prob]
    elseif d == NegativeBinomial
        nn = 10
        μ = GLM.linkinv.(l, x_float * true_b)
        clamp!(μ, -20, 20)
        prob = 1 ./ (1 .+ μ ./ nn)
        y = [rand(d(nn, i)) for i in prob] #number of failtures before nn success occurs
    end
    y = Float64.(y)

    #run IHT without groups
    k = 15
    ungrouped = fit_iht(y, x_float, z, J=1, k=k, d=d(), l=l)

    #run IHT with 5 groups each with 3 predictor
    J = 5
    k = 3
    grouped = fit_iht(y, x_float, z, J=J, k=k, d=d(), l=l, group=g)

    @test length(findall(!iszero, ungrouped.beta)) == 15
    @test length(findall(!iszero, grouped.beta)) == 15
end

@testset "Negative binomial nuisance parameter Float32 matrix" begin
	n = 1000
	p = 10000
	k = 10
	d = NegativeBinomial
	l = LogLink()

	# simulate SnpArrays data
	x = simulate_correlated_snparray(undef, n, p)
	xla = SnpLinAlg{Float64}(x, model=ADDITIVE_MODEL, center=true, scale=true)
	z = ones(n, 1)
	y, true_b, correct_position = simulate_random_response(xla, k, d, l, r=10);

	@time newton = fit_iht(y, xla, z, J=1, k=k, d=d(), l=l, est_r=:Newton)
	@test typeof(newton.d) == NegativeBinomial{Float64}
	@test newton.d.p == 0.5 # p parameter not used
	@test newton.d.r ≥ 1 # r converges to 10 faster

	@time mm = fit_iht(y, xla, z, J=1, k=k, d=d(), l=l, est_r=:MM)
	@test typeof(mm.d) == NegativeBinomial{Float64}
	@test mm.d.p == 0.5
	@test mm.d.r ≥ 1 # r converges to 10 slower

	# simulate floating point data
	T = Float32
	x = randn(T, n, p)
	z = ones(T, n)

	# simulate true model b, and the correct non-0 positions of b
	true_b = zeros(T, p)
	true_b[1:k] .= collect(0.1:0.1:1.0)
	true_c = [T.(4.0)]
	shuffle!(true_b)

	# simulate response
	r = 20
    μ = GLM.linkinv.(l, xla * true_b)
    clamp!(μ, -20, 20)
    prob = 1 ./ (1 .+ μ ./ r)
    y = [rand(d(r, Float64(i))) for i in prob]
    y = T.(y)

	d = d(r, Float32(0.5)) # need Float32 for eltype of d
	@time newton = fit_iht(y, x, z, J=1, k=k, d=d, l=l, est_r=:Newton)
	@test typeof(newton.d) == NegativeBinomial{Float32}
	@test newton.d.p == 0.5
	@test newton.d.r ≥ 0

	@time mm = fit_iht(y, x, z, J=1, k=k, d=d, l=l, est_r=:MM)
	@test typeof(mm.d) == NegativeBinomial{Float32}
	@test mm.d.p == 0.5
	@test mm.d.r ≥ 0
end

@testset "initialze beta" begin
	#simulat data with k true predictors, from distribution d and with link l.
	n = 1000
	p = 10000
	k = 10
	d = Normal
	l = canonicallink(d())

	#set random seed
	Random.seed!(1111)

	#construct SnpArraym, snpmatrix, and non genetic covariate (intercept)
	x = simulate_random_snparray(undef, n, p)
	xla = SnpLinAlg{Float64}(x, model=ADDITIVE_MODEL, center=true, scale=true)
	z = ones(n)

	# simulate response, true model b, and the correct non-0 positions of b
	y, true_b, correct_position = simulate_random_response(xla, k, d, l)

	#run with and without initializing beta
	@time result = fit_iht(y, xla, z, J=1, k=k, d=d(), l=l, init_beta=true)
	@test count(!iszero, result.beta) == 10
end


text:
# Arquivo contendo comparação cinemática, neste arquivo farei a geração de trajetória no domínio cartesiano e no domínio das juntas para comparação.
push!(LOAD_PATH, expanduser("/home/phelipe/Documentos"))
using Trajectory
using Plots
pyplot()

m=[23.902, 1.285]
l= [.45, .45]
r= [0.091, 0.048]
I= [1.266, 0.093]
g = 9.81
t0 = 0.
tend = 4.
Ts = 0.01

function cinematica_direta(θ::T) where T<:AbstractArray
    x = l[1] * sin(θ[1]) + l[2] * sin(θ[1] + θ[2])
    y = -(l[1] * cos(θ[1]) + l[2] * cos(θ[1] + θ[2]))
    [x,y]
end

function cinematica_inversa(cord::Vector{T}) where T<:AbstractFloat
    x = cord[1]
    y = cord[2]
    t = (x^2 + y^2 - l[1]^2 - l[2]^2) / (2 * l[1] * l[2])
    θ2 = atan2( (1 - t^2), t)
    k1 = l[1] + l[2] * cos(θ2)
    k2 = l[2] * sin(θ2)
    γ = atan2(k2,k1)
    θ1 = atan2(x,-y) - γ
    [θ1, θ2]
end

θ_start = [0.,0.]
#θ_end = [deg2rad(10),deg2rad(90)]
θ_end = [1.8, 2.0]
position_start = cinematica_direta(θ_start)
position_end = cinematica_direta(θ_end)

#teta1 = deg2rad(45)
#teta2 = deg2rad(0)
#println("Valores na cinemática direta $(cord = cinematica_direta([teta1,teta2])) metros")
#println("Valores na cinemática inversa $(rad2deg.(cinematica_inversa(cord))) deg")

#mínimo jerk cartesiano
x1,v1,a1,j1 = minimumjerkf(position_start[1],0.,0.,t0,position_end[1],0.,0.,tend)
x2,v2,a2,j2 = minimumjerkf(position_start[2],0.,0.,t0,position_end[2],0.,0.,tend)

time = collect(t0:Ts:tend)
x1data = map(a-> x1(a),time)
x2data = map(a-> x2(a),time)
function plotcartesiano()
    p1 = plot(time,x1data,xlabel="tempo",ylabel="posição x")
    p2 = plot(time,x2data,xlabel="tempo",ylabel="posição y")
    plot(p1,p2, tittle="Cartesiano")
end

θ = map((x,y)-> cinematica_inversa([x,y]),x1data,x2data)
θ1 = map(x->x[1],θ)
θ2 = map(x->x[2],θ)
vθ1 = diff(θ1)./Ts
vθ2 = diff(θ2)./Ts
aθ1 = diff(vθ1)./Ts
aθ2 = diff(vθ2)./Ts
jθ1 = diff(aθ1)./Ts
jθ2 = diff(aθ2)./Ts

function plotjuntas()
    p1 = plot(time,θ1,xlabel="tempo",ylabel="posição θ1")
    p2 = plot(time,θ2,xlabel="tempo",ylabel="posição θ2")
    plot(p1,p2, tittle="Jerk cartesiano nas juntas")
end

function plotvjuntas()
    p1 = plot(time[1:length(vθ1)],vθ1,xlabel="tempo",ylabel="velocidade θ1")
    p2 = plot(time[1:length(vθ1)],vθ2,xlabel="tempo",ylabel="velocidade θ2")
    plot(p1,p2, tittle="Jerk cartesiano nas juntas")
end

function plotajuntas()
    p1 = plot(time[1:length(aθ1)],aθ1,xlabel="tempo",ylabel="aceleração θ1")
    p2 = plot(time[1:length(aθ1)],aθ2,xlabel="tempo",ylabel="aceleração θ2")
    plot(p1,p2, tittle="Jerk cartesiano nas juntas")
end

function plotjjuntas()
    p1 = plot(time[1:length(jθ1)],jθ1,xlabel="tempo",ylabel="tranco θ1")
    p2 = plot(time[1:length(jθ2)],jθ2,xlabel="tempo",ylabel="tranco θ2")
    plot(p1,p2, tittle="Jerk cartesiano nas juntas")
end

# domínio das juntas
xj1,vj1,aj1,jj1 = minimumjerkf(θ_start[1],0.,0.,t0,θ_end[1],0.,0.,tend)
xj2,vj2,aj2,jj2 = minimumjerkf(θ_start[2],0.,0.,t0,θ_end[2],0.,0.,tend)

θ_junta1 = map(a-> xj1(a),time)
θ_junta2 = map(a-> xj2(a),time)
v_junta1 = map(a-> vj1(a),time)
v_junta2 = map(a-> vj2(a),time)
a_junta1 = map(a-> aj1(a),time)
a_junta2 = map(a-> aj2(a),time)
j_junta1 = map(a-> jj1(a),time)
j_junta2 = map(a-> jj2(a),time)

cart = map((x,y)-> cinematica_direta([x,y]),θ_junta1,θ_junta2)
cart_x = map(x->x[1],cart)
cart_y = map(x->x[2],cart)

function plotjerks()
    p1 = plot(time,[jθ1,j_junta1],xlabel="tempo",ylabel="tranco θ1")
    p2 = plot(time,[jθ2,j_junta2],xlabel="tempo",ylabel="tranco θ2")
    plot(p1,p2)
end

function plotpositions()
    p1 = plot(time,[θ1,θ_junta1],xlabel="tempo",ylabel="posição θ1")
    p2 = plot(time,[θ2,θ_junta2],xlabel="tempo",ylabel="posição θ2")
    plot(p1,p2)
end

function plotvelocitys()
    p1 = plot(time,[vθ1,v_junta1],xlabel="tempo",ylabel="velocidade θ1")
    p2 = plot(time,[vθ2,v_junta2],xlabel="tempo",ylabel="velocidade θ2")
    plot(p1,p2)
end

function plotacelerations()
    p1 = plot(time,[aθ1,a_junta1],xlabel="tempo",ylabel="aceleração θ1")
    p2 = plot(time,[aθ2,a_junta2],xlabel="tempo",ylabel="aceleração θ2")
    plot(p1,p2)
end

function plotxy()
    p1 = plot(x1data,x2data)
    p2 = plot(cart_x,cart_y)
    p = plot(p1,p2)
end

#Essas funções mostram o que acontece nas juntas quando o jerk é feito no domínio cartesiano
#plotjuntas()
#plotvjuntas()
#plotajuntas()
#plotjjuntas()

# Essas funções comparam as ações nas juntas
plotxy()
plotjerks()


text:
# function for converting back and forth between xyz and barycentric coordinates

"""
  This function calculates xyz coordinates given the barycentric corodinates
  and the vertices of the triangle xyz coordinates

  Inputs:
    bary: a dim x numpoints array of barycentric coordinates
    vtx: a nvertex x dim array of xyz coordinates of the vertices
         of the triangle (or tetrahedron

  Outputs:
    coords: a dim x numpoints array of xyz coordinates corresponding to
            the barycentric coordinates
"""
function baryToXY(bary::AbstractMatrix{T2}, vtx::AbstractMatrix) where T2

  @assert (size(bary, 1) >= 1 && size(bary, 1) <= 3)
  @assert (size(vtx, 2) >= 1 && size(vtx, 2) <= 3)
  @assert (size(vtx, 1) >= 2 && size(vtx, 1) <= 4)

  npoints = size(bary, 2)
  dim = size(bary, 1)
  coords = zeros(bary)

  T = zeros(dim, dim)
  r1 = vtx[1, :]

  for i=1:dim
    T[i, :] = vtx[i+1, :] - r1
  end

  # compute xyz coordinates
  for i=1:npoints
    coords[:, i] = T*bary[:, i] + r1
  end

  return coords
end

#=
function calc3dnodes(vtx::Array{Float64, 2})
# the columns of vtx are the x, y, and z coordinates of the nodes

  v1 = vtx[:, 2] - vtx[:, 1]
  v2 = vtx[:, 3] - vtx[:, 1]
  v3 = vtx[:, 4] - vtx[:, 1]
  xi = [1/8 5/8 1/8 1/8;
        1/8 1/8 5/8 1/8;
        1/8 1/8 1/8 5/8]
  T = zeros(3,3)
  T[:, 1] = v1
  T[:, 2] = v2
  T[:, 3] = v3

  coords = zeros(3, 4)
  r1 = vtx[:, 1]
  for i=1:4
    coords[:, i] = T*xi[:, i] + r1
  end

  return coords
end
#=
vtx = [1. 2 1 1;
       1  1 2 1;
       1  1 1 2]
coords = calcnodes(vtx)
println("vtx = \n", vtx)
println("coords = \n", coords)
=#
=#


text:
using Test, Random
import Flux: activations

@testset "basic" begin
  @testset "helpers" begin
    @testset "activations" begin
      dummy_model = Chain(x->x.^2, x->x .- 3, x -> tan.(x))
      x = randn(10)
      @test activations(dummy_model, x)[1] == x.^2
      @test activations(dummy_model, x)[2] == (x.^2 .- 3)
      @test activations(dummy_model, x)[3] == tan.(x.^2 .- 3)

      @test activations(Chain(), x) == ()
      @test activations(Chain(identity, x->:foo), x)[2] == :foo # results include `Any` type
    end
  end

  @testset "Chain" begin
    @test_nowarn Chain(Dense(10, 5, σ), Dense(5, 2))(randn(10))
    @test_throws DimensionMismatch Chain(Dense(10, 5, σ),Dense(2, 1))(randn(10))
    # numeric test should be put into testset of corresponding layer
  end

  @testset "Activations" begin
    c = Chain(Dense(3,5,relu), Dense(5,1,relu))
    X = Float32.([1.0; 1.0; 1.0])
    @test_nowarn gradient(()->Flux.activations(c, X)[2][1], params(c))
  end

  @testset "Dense" begin
    @testset "constructors" begin
      @test size(Dense(10, 100).weight) == (100, 10)
      @test size(Dense(10, 100).bias) == (100,)
      @test Dense(rand(100,10), rand(100)).σ == identity
      @test Dense(rand(100,10)).σ == identity

      @test Dense(rand(100,10), false).σ == identity
      @test Dense(rand(100,10), false, tanh).σ == tanh
      @test Dense(rand(100,10), rand(100)).σ == identity
      @test Dense(rand(Float16, 100,10), true).bias isa Vector{Float16}  # creates matching type
      @test_skip Dense(rand(Float16, 100,10), rand(100)).bias isa Vector{Float16}  # converts to match

      @test Dense(3,4; init=Base.randn, bias=true).bias isa Vector{Float64}
      @test_skip Dense(3,4; init=Base.randn, bias=[1,2,3,4]).bias isa Vector{Float64}

      @test_throws MethodError Dense(10, 10.5)
      @test_throws MethodError Dense(10, 10.5, tanh)
      @test_throws DimensionMismatch Dense(3,4; bias=rand(5))
      @test_throws DimensionMismatch Dense(rand(4,3), rand(5))
      @test_throws MethodError Dense(rand(5))
      @test_throws MethodError Dense(rand(5), rand(5))
      @test_throws MethodError Dense(rand(5), rand(5), tanh)
    end
    @testset "dimensions" begin
      @test  length(Dense(10, 5)(randn(10))) == 5
      @test_throws DimensionMismatch Dense(10, 5)(randn(1))
      @test_throws MethodError Dense(10, 5)(1) # avoid broadcasting
      @test_throws MethodError Dense(10, 5).(randn(10)) # avoid broadcasting
      @test size(Dense(10, 5)(randn(10))) == (5,)
      @test size(Dense(10, 5)(randn(10,2))) == (5,2)
      @test size(Dense(10, 5)(randn(10,2,3))) == (5,2,3)
      @test size(Dense(10, 5)(randn(10,2,3,4))) == (5,2,3,4)
      @test_throws DimensionMismatch Dense(10, 5)(randn(11,2,3))
    end
    @testset "zeros" begin
      @test Dense(10, 1, identity, init = ones)(ones(10,1)) == 10*ones(1, 1)
      @test Dense(10, 1, identity, init = ones)(ones(10,2)) == 10*ones(1, 2)
      @test Dense(10, 2, identity, init = ones)(ones(10,1)) == 10*ones(2, 1)
      @test Dense(10, 2, identity, init = ones)([ones(10,1) 2*ones(10,1)]) == [10 20; 10 20]
      @test Dense(10, 2, identity, init = ones, bias = false)([ones(10,1) 2*ones(10,1)]) == [10 20; 10 20]
    end
  end

  @testset "Diagonal" begin
    @test length(Flux.Diagonal(10)(randn(10))) == 10
    @test length(Flux.Diagonal(10)(1)) == 10
    @test length(Flux.Diagonal(10)(randn(1))) == 10
    @test_throws DimensionMismatch Flux.Diagonal(10)(randn(2))

    @test Flux.Diagonal(2)([1 2]) == [1 2; 1 2]
    @test Flux.Diagonal(2)([1,2]) == [1,2]
    @test Flux.Diagonal(2)([1 2; 3 4]) == [1 2; 3 4]

    @test Flux.Diagonal(2)(rand(2,3,4)) |> size == (2, 3, 4)
    @test Flux.Diagonal(2,3)(rand(2,3,4)) |> size == (2, 3, 4)
    @test Flux.Diagonal(2,3,4)(rand(2,3,4)) |> size == (2, 3, 4)
    @test Flux.Diagonal(2,3)(rand(2,1,4)) |> size == (2, 3, 4)
  end

  @testset "Maxout" begin
    # Note that the normal common usage of Maxout is as per the docstring
    # These are abnormal constructors used for testing purposes

    @testset "Constructor" begin
      mo = Maxout(() -> identity, 4)
      input = rand(40)
      @test mo(input) == input
    end

    @testset "simple alternatives" begin
      mo = Maxout((x -> x, x -> 2x, x -> 0.5x))
      input = rand(40)
      @test mo(input) == 2*input
    end

    @testset "complex alternatives" begin
      mo = Maxout((x -> [0.5; 0.1]*x, x -> [0.2; 0.7]*x))
      input = [3.0 2.0]
      target = [0.5, 0.7].*input
      @test mo(input) == target
    end

    @testset "params" begin
      mo = Maxout(()->Dense(32, 64), 4)
      ps = params(mo)
      @test length(ps) == 8  #4 alts, each with weight and bias
    end
  end

  @testset "SkipConnection" begin
    @testset "zero sum" begin
      input = randn(10, 10, 10, 10)
      @test SkipConnection(x -> zeros(size(x)), (a,b) -> a + b)(input) == input
    end

    @testset "concat size" begin
      input = randn(10, 2)
      @test size(SkipConnection(Dense(10,10), (a,b) -> cat(a, b, dims = 2))(input)) == (10,4)
    end
  end

  @testset "Bilinear" begin
    @testset "SkipConnection recombinator" begin
      d = Dense(10, 10)
      b = Flux.Bilinear(10, 10, 5)
      x = randn(Float32,10,9)
      sc = SkipConnection(d, b)
      @test size(sc(x)) == (5,9)
    end

    @testset "Two-streams zero sum" begin
      x = zeros(Float32,10,9)
      y = zeros(Float32,2,9)
      b = Flux.Bilinear(10, 2, 3)
      @test size(b(x,y)) == (3,9)
      @test sum(abs2, b(x,y)) == 0f0
    end

    @testset "Inner interactions" begin
      x = randn(Float32,11,7)
      b = Flux.Bilinear(11, 11, 3)
      @test size(b(x)) == (3,7)
      @test_nowarn gs = gradient(() -> sum(abs2.(b(x))), params(b))
    end

    @testset "constructors" begin
      b1 = Flux.Bilinear(randn(3,4,5))
      @test b1.bias isa Vector{Float64}
      @test b1.σ == identity

      b2 = Flux.Bilinear(randn(3,4,5), false)
      @test b2.bias == Flux.Zeros()

      b3 = Flux.Bilinear(randn(Float16, 3,4,5), true, tanh)
      @test b3.σ == tanh
      @test b3.bias isa Vector{Float16}
      @test size(b3(rand(4), rand(5))) == (3,)

      b4 = Flux.Bilinear(3,3,7; bias=1:7, init=Flux.zeros)
      @test_skip  b4.bias isa Vector{Float32}

      @test_throws ArgumentError Flux.Bilinear(rand(3)) # expects a 3-array
      @test_throws ArgumentError Flux.Bilinear(rand(3,4), false, tanh)
      @test_throws DimensionMismatch Flux.Bilinear(rand(3,4,5), rand(6), tanh) # wrong length bias
    end
  end

  @testset "Parallel" begin
    @testset "zero sum" begin
      input = randn(10, 10, 10, 10)
      @test Parallel(+, x -> zeros(size(x)), identity)(input) == input
    end

    @testset "concat size" begin
      input = randn(10, 2)
      @test size(Parallel((a, b) -> cat(a, b; dims=2), Dense(10, 10), identity)(input)) == (10, 4)
    end

    @testset "vararg input" begin
      inputs = randn(10), randn(5), randn(4)
      @test size(Parallel(+, Dense(10, 2), Dense(5, 2), Dense(4, 2))(inputs)) == (2,)
    end
  end
end


text:
"""
    Lasso(; alpha = 0.1, tol = 1e-4, mi = 1e+8)
Lasso Regression structure.

# Parameters
- `alpha` : leaarning rate.
- `th` : Threshold to judge that it has converged.
- `max` : Maximum number of learning.

# Example
```jldoctest regression
julia> model = Lasso()
Lasso(Float64[], 0.1, 0.0001, 100000000)

julia> fit!(model, x, t)
3-element Vector{Float64}:
   0.0
   0.5022766549841176
 154.43624186616267

julia> model(x)
20-element Vector{Float64}:
 188.64541774549468
 187.30088075704523
 181.5191726873298
 184.15748544986084
 185.72158909964372
 176.33798923891868
 185.80014722707335
 184.71565381947883
 189.20524796663838
 189.67502263476888
 189.50409373058318
 188.39535519538825
 188.481083670683
 188.88872347085172
 182.8477136378307
 188.64156231429416
 181.43996475587224
 188.9400571253936
 179.39836073711297
 185.6065850765288
```
"""
mutable struct Lasso
    w::Vector{Float64}
    α::Float64
    max::Int64
    th::Float64
    Lasso(; alpha = 0.1, max = 1e+8, th = 1e-4) = new(Array{Float64}(undef, 0), alpha, max, th)
end

soft_threshold(y, α) = sign.(y) .* max.(abs.(y) .- α, 0)

supermum_eigen(x) = maximum(sum(abs.(x), dims=1))

function fit!(model::Lasso, x, t)
    coverge(x, th) = @. abs(x) < th
    function update(w, x, t, α, rho)
        res = t - x * w
        return soft_threshold(w + (x' * res) / rho, α / rho)
    end
    x = expand(x)
    α = model.α * size(x, 1)
    w = zeros(size(x, 2))
    rho = supermum_eigen(x' * x)
    for _ in 1 : model.max
        w_new = update(w, x, t, α, rho)
        if coverge.(w_new - w, model.th) == trues(size(w)...)
            model.w = w_new
            return
        end
        w = w_new
    end
    @warn "Not Converged!"
    model.w = w
end

(model::Lasso)(x) = expand(x) * model.w


text:
#To show plot in browser
#using Gadfly

#To save plot as png
using Compose

#To plot graph
using Cairo, Fontconfig
using GraphPlot

#To create graphs
using Graphs
using Colors

function saveGraph(g, params, savePath, edgeMembership)

	#G = LightGraphs.SimpleGraph(g)

	colorsList = [RGBA(0, 0.2*(t-1), 0.1667*t, 1) for t in 1:params["l"]]

	colors = [colorsList[i] for i in edgeMembership]

	#Png
	draw(PNG(savePath, 8cm, 8cm), gplot(g, nodelabel=0:params["n"], edgestrokec=colors))

	#Web
	#gplot(g)
end

function saveCircuits(params, circuits, savePath)
	g = Graph(params["n"] + 1)

	edgeMembership = Int64[]

	currentIndex = 1

	for circuit in circuits

		for (i, j) in zip(circuit[1:end-1], circuit[2:end])
			add_edge!(g, i+1, j+1)
			push!(edgeMembership, currentIndex)
		end

		add_edge!(g, circuit[end]+1, circuit[1]+1)
		push!(edgeMembership, currentIndex)

	end

	saveGraph(g, params, savePath, edgeMembership)

end

function saveMultiCircuits(params, multiCircuits, savePath)

	g = Graph(params["n"] + 1)

	edgeMembership = Int64[]

	edgeToId = Dict{Tuple{Int, Int}, Int}()

	for (currentIndex, circuits) in enumerate(multiCircuits)

		for circuit in circuits

			for (i, j) in zip(circuit[1:end-1], circuit[2:end])
				if add_edge!(g, i+1, j+1)
					edgeToId[(i+1, j+1)] = currentIndex
					edgeToId[(j+1, i+1)] = currentIndex
				# else
				# 	edgeToId[(i+1, j+1)] += 1
				# 	edgeToId[(j+1, i+1)] += 1
				end
			end

			if add_edge!(g, circuit[end]+1, circuit[1]+1)
				edgeToId[(circuit[end]+1, circuit[1]+1)] = currentIndex
				edgeToId[(circuit[1]+1, circuit[end]+1)] = currentIndex
			# else
			# 	edgeToId[(circuit[end]+1, circuit[1]+1)] += 1
			# 	edgeToId[(circuit[1]+1, circuit[end]+1)] += 1
			end

		end

	end

	for e in edges(g)
		push!(edgeMembership, edgeToId[src(e), dst(e)])
		#edgeMembership
	end

	saveGraph(g, params, savePath, edgeMembership)

end


text:
using JLD

export savetraces, loadtraces

function savetraces(; traces=traceitems(), path="traces.jld", varname="traces")
    JLD.save(path, varname, traces)
end

function loadtraces(; path="traces.jld", varname="traces", asdefault=true)
    traces = JLD.load(path, varname)
    asdefault && (global happysnaps = traces)
    traces
end

max_stars_repo_name: JobJob/RickTracy


text:
function servercanlisten(server::WebsocketServer, port::Int = 8080)
    ended = Condition()

    listen(server, :connectError) do  err
        notify(ended, err.msg)
    end
    listen(server, :listening) do detail
        close(server)
    end
    listen(server, :client) do client
    end
    listen(server, :closed) do detail
        notify(ended, detail)
    end

    @async serve(server, port)

    wait(ended)
end

function echoserver(server::WebsocketServer, port::Int = 8080)
    started = Condition()

    listen(server, :client) do ws
        listen(ws, :message) do message
            send(ws, message)
        end
    end
    listen(server, :listening) do detail
        notify(started, server)
    end

    @async serve(server, port)

    wait(started)
end

max_stars_repo_name: citkane/SimpleWebsockets

max_issues_repo_name: citkane/SimpleWebsockets

max_forks_repo_name: citkane/SimpleWebsockets


text:
###
  Name     : c6_14_remove_uudate_package.jl
  Book     : Hands-on Data Science with Anaconda )
  Publisher: Packt Publishing Ltd.
  Author   : Yuxing Yan and James Yan
  Date     : 1/25/2018
  email    : yany@canisius.edu
             paulyxy@hotmail.com
###

Pkg.rm("AbstractTable")

Pkg.update()


text:
module QJuliaFieldUtils

using QJuliaGrid
using QJuliaFields
#using QJuliaUtils
using QJuliaGaugeUtils
using QJuliaEnums
using QJuliaInterface

@inline function gen_random_spinor!(field::QJuliaFields.QJuliaGenericField_qj)

  if(field.field_desc.geom != QJuliaEnums.QJULIA_SCALAR_GEOMETRY); error("Cannot apply on fields with non-scalar geometry.");end
  if(field.field_desc.grid.register_type != ComplexF64 && field.field_desc.grid.register_type != ComplexF32 )
    error("Register type ", field.field_desc.register_type, " is currently not supported.")
  end

  for i in LinearIndices(field.v); field.v[i] = field.field_desc.register_type(rand(), rand()); end

end #gen_random_spinor!

function construct_gauge_field!(field::QJuliaFields.QJuliaGenericField_qj, gtype, param::QJuliaInterface.QJuliaGaugeParam_qj)

  if(field.field_desc.geom != QJuliaEnums.QJULIA_VECTOR_GEOMETRY); error("Cannot apply on fields with non-vector geometry.");end
  if(field.field_desc.grid.register_type != ComplexF64 && field.field_desc.grid.register_type != ComplexF32 )
    error("Register type ", field.field_desc.register_type, " is currently not supported.")
  end
  QJuliaGaugeUtils.construct_gauge_field!(field.v, gtype, param)

end

end #QJuliaFieldUtils

max_stars_repo_name: alexstrel/QJulia

max_issues_repo_name: alexstrel/QJulia

max_forks_repo_name: alexstrel/QJulia


text:
push!(LOAD_PATH,"/Users/jannitta/.julia/dev/MyPkg/src/")
using Documenter, MyPkg

makedocs(sitename="My Documentation")

max_stars_repo_name: jannittayao/MyPkg

max_issues_repo_name: jannittayao/MyPkg

max_forks_repo_name: jannittayao/MyPkg


text:
function InterpolatingTransformProvider(arg0::TransformProvider, arg1::CartesianDerivativesFilter, arg2::AngularDerivativesFilter, arg3::jint, arg4::jdouble, arg5::jint, arg6::jdouble, arg7::jdouble)
    return InterpolatingTransformProvider((TransformProvider, CartesianDerivativesFilter, AngularDerivativesFilter, jint, jdouble, jint, jdouble, jdouble), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7)
end

function get_grid_points(obj::InterpolatingTransformProvider)
    return jcall(obj, "getGridPoints", jint, ())
end

function get_raw_provider(obj::InterpolatingTransformProvider)
    return jcall(obj, "getRawProvider", TransformProvider, ())
end

function get_step(obj::InterpolatingTransformProvider)
    return jcall(obj, "getStep", jdouble, ())
end

function get_transform(obj::InterpolatingTransformProvider, arg0::AbsoluteDate)
    return jcall(obj, "getTransform", Transform, (AbsoluteDate,), arg0)
end

function get_transform(obj::InterpolatingTransformProvider, arg0::FieldAbsoluteDate)
    return jcall(obj, "getTransform", FieldTransform, (FieldAbsoluteDate,), arg0)
end


text:
# The code to generate Fig 1A
using PyPlot
using CSV
using SurNoR_2020
using Statistics
using HypothesisTests

PyPlot.svg(true)
rcParams = PyPlot.PyDict(PyPlot.matplotlib."rcParams")
rcParams["svg.fonttype"] = "none"
rcParams["pdf.fonttype"] = 42

# ------------------------------------------------------------------------------
# Behavioral data
# ------------------------------------------------------------------------------
BehavData = CSV.read("./data/BehavData_SwitchState.csv")

Sub_Num = 12
Sub_Set = [1,2,3,4,5,6,7,8,9,11,12,13]
Inputs = Array{SurNoR_2020.Str_Input,3}(undef,Sub_Num,2,5)
for n = 1:Sub_Num
  for Block = 1:2
    for Epi = 1:5
      Inputs[n,Block,Epi] = SurNoR_2020.Func_data_to_input(BehavData; Sub=Sub_Set[n],
                                                           Epi=Epi, Block=Block)
    end
  end
end

# ------------------------------------------------------------------------------
# Lenght of episodes
# ------------------------------------------------------------------------------
Len_Mean, Len_Std, Len_Matrix = SurNoR_2020.Func_episode_lenght(Inputs)

println("The t-test result for comparing average number of actiopns in B1E1 and B2E2:")
@show OneSampleTTest(Len_Matrix[:,1,1], Len_Matrix[:,1,2])

fig = figure(); ax = gca()

x = ((1:1:5))
y=Len_Mean[:,:,1]
dy = Len_Std[:,:,1]/sqrt(length(Len_Matrix[:,1,1]))
bp1 = ax.bar(x, y[:] , yerr=dy[:], color="#167CB7")

x = ((6:1:10))
y=Len_Mean[:,:,2]
dy = Len_Std[:,:,2]/sqrt(length(Len_Matrix[:,1,1]))
bp2 = ax.bar(x, y[:] , yerr=dy[:], color="#B71616")
ax.set_xticks(1:1:10)
ax.set_xticklabels(["E1","E2","E3","E4","E5","E1","E2","E3","E4","E5"])
legend(["Block 1","Block 2"])

x = Array(1:1:10)
σ = 0.1
for j = 1:Sub_Num
  x_1 = x .+ 2*σ*(rand() - 0.5)
  y = Len_Matrix[j,:,:]
  ax.plot(x_1, y[:],".",color="k",alpha=0.5)
end
ax.set_ylim([0,250])

ax.set_title("Number of Steps of Subjects")
ax.set_xlabel("Episode")
ax.set_ylabel("Number of Steps")


text:
@testset "D" begin
    @testset "dafac/dafdc/dafec" begin
        try
            dafpath = tempname()
            handle = ckopn(dafpath, "TEST_ex_dafac", 140)
            cmnts = ["a", "bc", "def", "ghij"]
            dafac(handle, cmnts)
            dafcls(handle)
            handle = dafopr(dafpath)
            cmnts_out = dafec(handle)
            @test cmnts_out == cmnts
            cmnts_out = dafec(handle, bufsiz=1)
            @test cmnts_out == cmnts
            dafcls(handle)
            handle = dafopw(dafpath)
            dafac(handle, cmnts)
            cmnts_out = dafec(handle)
            @test cmnts_out == repeat(cmnts, 2)
            dafdc(handle)
            cmnts_out = dafec(handle)
            @test cmnts_out == String[]
        finally
            kclear()
        end
    end
    @testset "dafbbs/daffpa" begin
        try
            handle = dafopr(path(CORE, :spk))
            dafbbs(handle)
            found = daffpa()
            @test found
        finally
            kclear()
        end
    end
    @testset "dafbfs/daffna" begin
        try
            handle = dafopr(path(CORE, :spk))
            dafbfs(handle)
            found = daffna()
            @test found
        finally
            kclear()
        end
    end
    @testset "dafcs" begin
        try
            handle = dafopr(path(CORE, :spk))
            dafbbs(handle)
            dafcs(handle)
            found = daffpa()
            @test found
        finally
            kclear()
        end
    end
    @testset "dafgda" begin
        try
            handle = dafopr(path(CORE, :spk))
            elements = dafgda(handle, 1026, 1027)
            @test elements == [691200.0 / 2]
        finally
            kclear()
        end
    end
    @testset "dafgh" begin
        try
            handle = dafopr(path(CORE, :spk))
            dafbbs(handle)
            dafcs(handle)
            search_handle = dafgh()
            @test search_handle == handle
        finally
            kclear()
        end
    end
    @testset "dafgn/dafgs" begin
        try
            handle = dafopr(path(CORE, :spk))
            dafbfs(handle)
            found = daffna()
            @test found
            out = dafgs(2)
            @test out ≈ [-9.46511378160646408796e+07, 3.15662463183953464031e+08]
            outname = dafgn(100)
            @test outname == "DE-405"
        finally
            kclear()
        end
    end
    @testset "dafgsr/dafrfr/dafus" begin
        try
            handle = dafopr(path(CORE, :spk))
            nd, ni, ifname, fward, bward, free = dafrfr(handle)
            @test nd == 2
            @test ni == 6
            ss = nd + ((ni+1) >> 1)

            while fward > 0
                recno = fward
                fward, bward, n_ss = drec = Int.(dafgsr(handle, recno, 1, 3))
                @test recno == 7
                @test fward == 0
                @test bward == 0
                @test n_ss == 15
                first_word = 4
                last_end_word = 1024
                for _ in 1:n_ss
                    drec = dafgsr(handle, recno, first_word, first_word + ss - 1)
                    dc, ic = dafus(drec, nd, ni)
                    body, center, frame, spk_type, start_word, end_word = ic
                    @test dc ≈ [-9.46511378160646408796e+07, 3.15662463183953464031e+08]
                    @test body ÷ 100 == center
                    @test frame == 1
                    @test spk_type == 2
                    @test last_end_word + 1 == start_word
                    first_word += ss
                    last_end_word = end_word
                end
                @test fward == 0
            end
        finally
            kclear()
        end
    end
    @testset "dafps/dafrs" begin
        try
            dafpath = tempname()
            ifname = "Test CK type 1 segment created by cspice_ckw01"
            handle = ckopn(dafpath, ifname, 10)
            ckw01(handle, 1.1, 4.1, -77701, "J2000", "Test type 1 CK segment",
                  [1.1, 4.1], [[1.0, 1.0, 1.0, 1.0], [2.0, 2.0, 2.0, 2.0]],
                  [[0.0, 0.0, 1.0], [0.0, 0.0, 2.0]])
            ckcls(handle)
            kclear()

            handle = dafopw(dafpath)
            dafbfs(handle)
            found = daffna()
            @test found
            out = dafgs(124)
            dc, ic = dafus(out, 2, 6)
            ic[1] = -1999
            ic[2] = -2999
            summ = dafps(dc, ic)
            dafrs(summ)
            dafcls(handle)
            kclear()

            handle = dafopr(dafpath)
            dafbfs(handle)
            found = daffna()
            @test found
            out = dafgs(124)
            dc, ic = dafus(out, 2, 6)
            @test ic[1] == -1999
            @test ic[2] == -2999
        finally
            kclear()
        end
    end
    @testset "dasac/dascl/dasopr/dasopw/dasec/dasdc" begin
        try
            daspath = tempname()
            handle = dskopn(daspath, "ex_dasac", 140)
            dasac(handle, ["spice", "naif", "julia"])
            dascls(handle)
            kclear()

            handle = dasopr(daspath)
            @test dashfn(handle) == daspath
            comments = dasec(handle, bufsiz=3)
            @test length(comments) == 3
            @test comments  == ["spice", "naif", "julia"]
            dascls(handle)

            handle = dasopr(daspath)
            idword, ifname, nresvr, nresvc, ncomr, ncomc = dasrfr(handle)
            @test idword == "DAS/DSK"
            @test ifname == "ex_dasac"
            @test nresvr == 0
            @test nresvc == 0
            @test ncomr  == 1
            @test ncomc  == 17
            dascls(handle)

            handle = dasopw(daspath)
            dasdc(handle)
            dascls(handle)
            handle = dasopr(daspath)
            comments = dasec(handle)
            @test comments == []
            dascls(handle)
        finally
            kclear()
        end
    end
    @testset "dcyldr" begin
        output = dcyldr(1.0, 0.0, 0.0)
        expected = [1.0 0.0 0.0;
                    0.0 1.0 0.0;
                    0.0 0.0 1.0]
        @test output ≈ expected
        @test_throws SpiceError dcyldr(0.0, 0.0, 1.0)
    end
    @testset "deltet" begin
        try
            furnsh(path(CORE, :lsk))
            UTC_1997 = "Jan 1 1997"
            UTC_2004 = "Jan 1 2004"
            et_1997 = str2et(UTC_1997)
            et_2004 = str2et(UTC_2004)
            delt_1997 = deltet(et_1997, "ET")
            delt_2004 = deltet(et_2004, "ET")
            @test delt_1997 ≈ 62.1839353
            @test delt_2004 ≈ 64.1839116
        finally
            kclear()
        end
    end
    @testset "dgeodr" begin
        try
            furnsh(path(CORE, :pck))
            radii = bodvrd("EARTH", "RADII", 3)
            flat = (radii[1] - radii[3]) / radii[1]
            lon = deg2rad(118.0)
            lat = deg2rad(32.0)
            alt = 0.0
            rec = latrec(lon, lat, alt)
            output = dgeodr(rec[1], rec[2], rec[3], radii[1], flat)
            expected = [-0.25730624850202866 0.41177607401581356 0.0;
                        -0.019818463887750683 -0.012383950685377182 0.0011247386599188864;
                        0.040768073853231314 0.02547471988726025 0.9988438330394612]
            @test output ≈ expected
        finally
            kclear()
        end
    end
    @testset "diags2" begin
        mat = [1.0 4.0; 4.0 -5.0]
        diag, rot = diags2(mat)
        expected_diag = [3.0 0.0; 0.0 -7.0]
        expected_rot = [0.89442719 -0.44721360; 0.44721360 0.89442719]
        @test diag ≈ expected_diag
        @test rot ≈ expected_rot
    end
    @testset "diff" begin
        cell1 = SpiceIntCell(8)
        cell2 = SpiceIntCell(8)
        insrti!(cell1, 1)
        insrti!(cell1, 2)
        insrti!(cell1, 3)
        insrti!(cell2, 2)
        insrti!(cell2, 3)
        insrti!(cell2, 4)
        out = diff(cell1, cell2)
        @test out == [1]
        out = diff(cell2, cell1)
        @test out == [4]

        cell1 = SpiceCharCell(8, 8)
        cell2 = SpiceCharCell(8, 8)
        insrtc!(cell1, "1")
        insrtc!(cell1, "2")
        insrtc!(cell1, "3")
        insrtc!(cell2, "2")
        insrtc!(cell2, "3")
        insrtc!(cell2, "4")
        out = diff(cell1, cell2)
        @test out == ["1"]
        out = diff(cell2, cell1)
        @test out == ["4"]

        cell1 = SpiceDoubleCell(8)
        cell2 = SpiceDoubleCell(8)
        insrtd!(cell1, 1.0)
        insrtd!(cell1, 2.0)
        insrtd!(cell1, 3.0)
        insrtd!(cell2, 2.0)
        insrtd!(cell2, 3.0)
        insrtd!(cell2, 4.0)
        out = diff(cell1, cell2)
        @test out == [1.0]
        out = diff(cell2, cell1)
        @test out == [4.0]
    end
    @testset "dlabfs" begin
        try
            handle = dasopr(path(EXTRA, :phobos_dsk))
            current = dlabfs(handle)
            @test current.dsize == 1300
            @test dlafns(handle, current) === nothing
            dascls(handle)
        finally
            kclear()
        end
    end
    @testset "dlabbs" begin
        try
            handle = dasopr(path(EXTRA, :phobos_dsk))
            current = dlabbs(handle)
            @test current.dsize == 1300
            @test dlafps(handle, current) === nothing
            dascls(handle)
        finally
            kclear()
        end
    end
    @testset "dlatdr" begin
        output = dlatdr(1.0, 0.0, 0.0)
        expected = [1.0 0.0 0.0;
                    0.0 1.0 0.0;
                    0.0 0.0 1.0]
        @test output ≈ expected
    end
    @testset "dp2hx" begin
        @test dp2hx(2.0e-9) == "89705F4136B4A8^-7"
        @test dp2hx(1.0) == "1^1"
        @test dp2hx(-1.0) == "-1^1"
        @test dp2hx(1024.0) == "4^3"
        @test dp2hx(-1024.0) == "-4^3"
        @test dp2hx(521707.0) == "7F5EB^5"
        @test dp2hx(27.0) == "1B^2"
        @test dp2hx(0.0) == "0^0"
    end
    @testset "dpgrdr" begin
        try
            furnsh(path(CORE, :pck))
            radii = bodvrd("MARS", "RADII", 3)
            re = radii[1]
            rp = radii[3]
            f = (re - rp) / re
            output = dpgrdr("Mars", deg2rad(90.0), deg2rad(45.0), 300, re, f)
            expected = [0.25464790894703276 -0.5092958178940655 -0.0;
                        -0.002629849831988239 -0.0013149249159941194 1.5182979166821334e-05;
                        0.004618598844358383 0.0023092994221791917 0.9999866677515724]
            @test output ≈ expected
        finally
            kclear()
        end
    end
    @testset "dpmax" begin
        @test SPICE._dpmax() == prevfloat(typemax(Float64))
    end
    @testset "dpmin" begin
        @test SPICE._dpmin() == nextfloat(typemin(Float64))
    end
    @testset "dpr" begin
        @test SPICE._dpr() == rad2deg(1.0)
    end
    @testset "drdcyl" begin
        output = drdcyl(1.0, deg2rad(180.0), 1.0)
        expected = [-1.0 0.0 0.0;
                    0.0 -1.0 0.0;
                    0.0 0.0 1.0]
        @test output ≈ expected
    end
    @testset "drdgeo" begin
        try
            furnsh(path(CORE, :pck))
            radii = bodvrd("EARTH", "RADII", 3)
            flat = (radii[1] - radii[3]) / radii[1]
            lon = deg2rad(118.0)
            lat = deg2rad(32.0)
            alt = 0.0
            output = drdgeo(lon, lat, alt, radii[1], flat)
            expected = [-4780.329375996193 1580.5982261675397 -0.3981344650201568;
                        -2541.7462156656084 -2972.6729150327574 0.7487820251299121;
                        0.0 5387.9427815962445 0.5299192642332049]
            @test output ≈ expected
        finally
            kclear()
        end
    end
    @testset "drdlat" begin
        output = drdlat(1.0, deg2rad(90.0), 0.0)
        expected = [0.0 -1.0 -0.0;
                    1.0 0.0 -0.0;
                    0.0 0.0 1.0]
        @test output ≈ expected
    end
    @testset "drdpgr" begin
        try
            furnsh(path(CORE, :pck))
            radii = bodvrd("MARS", "RADII", 3)
            re = radii[1]
            rp = radii[3]
            f = (re - rp) / re
            output = drdpgr("MARS", deg2rad(90.0), deg2rad(45), 300, re, f)
            expected = [-2620.6789148181783 0.0 0.0;
                        0.0 2606.460468253308 -0.7071067811865476;
                        -0.0 2606.460468253308 0.7071067811865475]
            @test output ≈ expected
        finally
            kclear()
        end
    end
    @testset "drdsph" begin
        output = drdsph(1.0, π/2, π)
        expected = [-1.0 0.0 0.0;
                    0.0 0.0 -1.0;
                    0.0 -1.0 0.0]
        @test output ≈ expected
    end
    @testset "dskgtl/dskstl" begin
        SPICE_DSK_KEYXFR = 1
        @test dskgtl(SPICE_DSK_KEYXFR) ≈ 1.0e-10
        dskstl(SPICE_DSK_KEYXFR, 1.0e-8)
        @test dskgtl(SPICE_DSK_KEYXFR) ≈ 1.0e-8
        dskstl(SPICE_DSK_KEYXFR, 1.0e-10)
        @test dskgtl(SPICE_DSK_KEYXFR) ≈ 1.0e-10
    end
    @testset "dskobj/dsksrf" begin
        try
            bodyids = SpiceIntCell(100)
            dskobj!(bodyids, path(EXTRA, :phobos_dsk))
            @test 401 in bodyids
            srfids = SpiceIntCell(100)
            dsksrf!(srfids, path(EXTRA, :phobos_dsk), 401)
            @test 401 in srfids
        finally
            kclear()
        end
    end
    @testset "dskb02" begin
        try
            handle = dasopr(path(EXTRA, :phobos_dsk))
            dladsc = dlabfs(handle)
            res = dskb02(handle, dladsc)
            nv, nump, nvxtot, vtxbds, voxsiz, voxori, vgrext, cgscal, vtxnpl, voxnpt, voxnpl = res
            @test nv == 422
            @test nump == 840
            @test nvxtot == 8232
            @test cgscal == 7
            @test vtxnpl == 0
            @test voxnpt == 2744
            @test voxnpl == 3257
            @test voxsiz ≈ 3.320691339664286
        finally
            kclear()
        end
    end
    @testset "dskd02" begin
        try
            handle = dasopr(path(EXTRA, :phobos_dsk))
            dladsc = dlabfs(handle)
            values = dskd02(handle, dladsc, 19, 1, 3)
            expected = [5.12656957900699912362e-16,  -0.0, -8.37260000000000026432e+00]
            @test values ≈ expected
        finally
            kclear()
        end
    end
    @testset "dskgd" begin
        try
            handle = dasopr(path(EXTRA, :phobos_dsk))
            dladsc = dlabfs(handle)
            dskdsc = dskgd(handle, dladsc)
            @test dskdsc.surfce == 401
            @test dskdsc.center == 401
            @test dskdsc.dclass == 1
            @test dskdsc.dtype  == 2
            @test dskdsc.frmcde == 10021
            @test dskdsc.corsys == 1
            @test collect(dskdsc.corpar) ≈ zeros(10)
            @test dskdsc.co1min ≈ -3.141593 rtol=1e-6
            @test dskdsc.co1max ≈ 3.141593 rtol=1e-6
            @test dskdsc.co2min ≈ -1.570796 rtol=1e-6
            @test dskdsc.co2max ≈ 1.570796 rtol=1e-6
            @test dskdsc.co3min ≈ 8.181895873588292 rtol=1e-6
            @test dskdsc.co3max ≈ 13.89340000000111 rtol=1e-6
            @test dskdsc.start  ≈ -1577879958.816059 rtol=1e-6
            @test dskdsc.stop   ≈ 1577880066.183913 rtol=1e-6
        finally
            kclear()
        end
    end
    @testset "dski02" begin
        try
            handle = dasopr(path(EXTRA, :phobos_dsk))
            dladsc = dlabfs(handle)
            plates = dski02(handle, dladsc, 2, 1, 3)
            @test length(plates) > 0
        finally
            kclear()
        end
    end
    @testset "dskn02" begin
        try
            handle = dasopr(path(EXTRA, :phobos_dsk))
            dladsc = dlabfs(handle)
            normal = dskn02(handle, dladsc, 1)
            expected = [0.20813166897151150203,  0.07187012861854354118, -0.97545676120650637309]
            @test normal ≈ expected
        finally
            kclear()
        end
    end
    @testset "dskp02" begin
        try
            handle = dasopr(path(EXTRA, :phobos_dsk))
            dladsc = dlabfs(handle)
            plates = dskp02(handle, dladsc, 2, 2)
            @test plates[1] == [1, 9, 2]
            @test plates[2] == [1, 2, 3]
        finally
            kclear()
        end
    end
    @testset "dskv02" begin
        try
            handle = dasopr(path(EXTRA, :phobos_dsk))
            dladsc = dlabfs(handle)
            vrtces = dskv02(handle, dladsc, 2, 1)
            expected = [5.12656957900699912362e-16, -0.0, -8.37260000000000026432]
            @test vrtces[1] ≈ expected
        finally
            kclear()
        end
    end
    @testset "dskw02/dskrb2/dskmi2" begin
        try
            dskpath = tempname()
            handle = dasopr(path(EXTRA, :phobos_dsk))
            dladsc = dlabfs(handle)
            finscl = 5.0
            corscl = 4
            center = 401
            surfid = 1
            dclass = 2
            frame = "IAU_PHOBOS"
            first = -50 * jyear()
            last  =  50 * jyear()
            SPICE_DSK02_MAXVRT = 16000002 ÷ 128
            SPICE_DSK02_MAXPLT = 2 * (SPICE_DSK02_MAXVRT - 2)
            SPICE_DSK02_MAXVXP = SPICE_DSK02_MAXPLT ÷ 2
            SPICE_DSK02_MAXCEL = 60000000 ÷ 128
            SPICE_DSK02_MXNVLS = SPICE_DSK02_MAXCEL + (SPICE_DSK02_MAXVXP ÷ 2)
            SPICE_DSK02_MAXCGR = 100000 ÷ 128
            SPICE_DSK02_IXIFIX = SPICE_DSK02_MAXCGR + 7
            SPICE_DSK02_MAXNPV = 3 * (SPICE_DSK02_MAXPLT ÷ 2) + 1
            SPICE_DSK02_SPAISZ = (SPICE_DSK02_IXIFIX + SPICE_DSK02_MAXVXP + SPICE_DSK02_MXNVLS
                                  + SPICE_DSK02_MAXVRT + SPICE_DSK02_MAXNPV)
            worksz = SPICE_DSK02_MAXCEL
            voxpsz = SPICE_DSK02_MAXVXP
            voxlsz = SPICE_DSK02_MXNVLS
            spaisz = SPICE_DSK02_SPAISZ
            vrtces = dskv02(handle, dladsc, 2, 422)
            plates = dskp02(handle, dladsc, 2, 840)
            dskcls(handle)
            kclear()

            handle = dskopn(dskpath, "TESTdskw02.dsk/AA/29-SEP-2017", 0)
            spaixd, spaixi = dskmi2(vrtces, plates, finscl, corscl, worksz, voxpsz, voxlsz, false, spaisz)
            corsys = 1
            mncor1 = -π
            mxcor1 = π
            mncor2 = -π/2
            mxcor2 = π/2
            corpar = zeros(10)
            mncor3, mxcor3 = dskrb2(vrtces, plates, corsys, corpar)
            dskw02(handle, center, surfid, dclass, frame, corsys, corpar,
                   mncor1, mxcor1, mncor2, mxcor2, mncor3, mxcor3, first,
                   last, vrtces, plates, spaixd, spaixi)
            dskcls(handle, true)
        finally
            kclear()
        end
    end
    @testset "dskx02" begin
        try
            handle = dasopr(path(EXTRA, :phobos_dsk))
            dladsc = dlabfs(handle)
            dskdsc = dskgd(handle, dladsc)
            r = 2.0 * dskdsc.co3max
            vertex = latrec(r, 0.0, 0.0)
            raydir = -vertex
            plid, xpt = dskx02(handle, dladsc, vertex, raydir)
            # TODO: The results of this function are inconsistent for different versions of SPICE
            @test_skip plid == 421
            @test xpt ≈ [12.36679999999999957083, 0.0, 0.0]
        finally
            kclear()
        end
    end
    @testset "dskxsi" begin
        try
            furnsh(path(EXTRA, :phobos_dsk))
            dsk1, filtyp, source, handle = kdata(1, "DSK")
            dladsc = dlabfs(handle)
            dskdsc = dskgd(handle, dladsc)
            target = bodc2n(dskdsc.center)
            fixref = frmnam(dskdsc.frmcde)
            r = 1.0e10
            vertex = latrec(r, 0.0, 0.0)
            raydir = -vertex
            srflst = [dskdsc.surfce]
            xpt, handle, dladsc2, dskdsc2, dc, ic  = dskxsi(false, target, srflst, 0.0,
                                                            fixref, vertex, raydir)
            # TODO: The results of this function are inconsistent for different versions of SPICE
            @test_skip ic[1] == 420
            @test dc[1] == 0.0
            @test xpt ≈ [12.36679999999999957083, 0.0, 0.0]
        finally
            kclear()
        end
    end
    @testset "dskxv" begin
        try
            furnsh(path(EXTRA, :phobos_dsk))
            dsk1, filtyp, source, handle = kdata(1, "DSK")
            dladsc = dlabfs(handle)
            dskdsc = dskgd(handle, dladsc)
            target = bodc2n(dskdsc.center)
            fixref = frmnam(dskdsc.frmcde)
            r = 1.0e10
            vertex = latrec(r, 0.0, 0.0)
            raydir = -vertex
            srflst = [dskdsc.surfce]
            xpt, foundarray = dskxv(false, target, srflst, 0.0, fixref, [vertex], [raydir])
            @test length(xpt) == 1
            @test length(foundarray) == 1
            @test foundarray[1]
            @test xpt[1] ≈ [12.36679999999999957083, 0.0, 0.0]
        finally
            kclear()
        end
    end
    @testset "dskz02" begin
        try
            handle = dasopr(path(EXTRA, :phobos_dsk))
            dladsc = dlabfs(handle)
            nv, nplates = dskz02(handle, dladsc)
            @test nv > 0
            @test nplates > 0
        finally
            kclear()
        end
    end
    @testset "dsphdr" begin
        output = dsphdr(-1.0, 0.0, 0.0)
        expected = [-1.0 0.0 0.0;
                    0.0 0.0 -1.0;
                    0.0 -1.0 0.0]
        @test output ≈ expected
    end
    @testset "dtpool" begin
        try
            lmpoolNames = ["DELTET/DELTA_T_A", "DELTET/K", "DELTET/EB",
                           "DELTET/M", "DELTET/DELTA_AT"]
            lmpoolLens = [1, 1, 1, 2, 46]
            textbuf = ["DELTET/DELTA_T_A = 32.184",
                       "DELTET/K         = 1.657D-3",
                       "DELTET/EB        = 1.671D-2",
                       "DELTET/M         = ( 6.239996 1.99096871D-7 )",
                       "DELTET/DELTA_AT  = ( 10, @1972-JAN-1",
                       "                     11, @1972-JUL-1",
                       "                     12, @1973-JAN-1",
                       "                     13, @1974-JAN-1",
                       "                     14, @1975-JAN-1",
                       "                     15, @1976-JAN-1",
                       "                     16, @1977-JAN-1",
                       "                     17, @1978-JAN-1",
                       "                     18, @1979-JAN-1",
                       "                     19, @1980-JAN-1",
                       "                     20, @1981-JUL-1",
                       "                     21, @1982-JUL-1",
                       "                     22, @1983-JUL-1",
                       "                     23, @1985-JUL-1",
                       "                     24, @1988-JAN-1",
                       "                     25, @1990-JAN-1",
                       "                     26, @1991-JAN-1",
                       "                     27, @1992-JUL-1",
                       "                     28, @1993-JUL-1",
                       "                     29, @1994-JUL-1",
                       "                     30, @1996-JAN-1",
                       "                     31, @1997-JUL-1",
                       "                     32, @1999-JAN-1 )"]
            lmpool(textbuf)
            for (var, expectLen) in zip(lmpoolNames, lmpoolLens)
                n, vartype = dtpool(var)
                @test expectLen == n
                @test vartype == :N
            end
        finally
            kclear()
        end
    end
    @testset "ducrss" begin
        try
            furnsh(path(CORE, :lsk), path(CORE, :pck), path(CORE, :spk))
            z_earth = [0.0, 0.0, 1.0, 0.0, 0.0, 0.0]
            et = str2et("Jan 1, 2009")
            trans = sxform("IAU_EARTH", "J2000", et)
            z_j2000 = trans * z_earth
            state, ltime = spkezr("SUN", et, "J2000", "LT+S", "EARTH")
            z_new = ducrss(state, z_j2000)
            z_expected = [-0.9798625180326394, -0.1996715076226282, 0.0008572038510904833,
                          4.453114222872359e-08, -2.1853106962531453e-07, -3.6140021238340607e-11]
            @test z_new ≈ z_expected
        finally
            kclear()
        end
    end
    @testset "dvcrss" begin
        try
            furnsh(path(CORE, :lsk), path(CORE, :pck), path(CORE, :spk))
            z_earth = [0.0, 0.0, 1.0, 0.0, 0.0, 0.0]
            et = str2et("Jan 1, 2009")
            trans = sxform("IAU_EARTH", "J2000", et)
            z_j2000 = trans * z_earth
            state, ltime = spkezr("SUN", et, "J2000", "LT+S", "EARTH")
            z = dvcrss(state, z_j2000)
            expected = [-1.32672690582546606660e+08,  -2.70353812480484284461e+07,
                        1.16064793997540167766e+05,   5.12510726479525757782e+00,
                        -2.97732415336074147660e+01,  -4.10216496370272454969e-03]
            @test z ≈ expected
        finally
            kclear()
        end
    end
    @testset "dvdot" begin
        @test dvdot([1.0, 0.0, 1.0, 0.0, 1.0, 0.0], [0.0, 1.0, 0.0, 1.0, 0.0, 1.0]) == 3.0
    end
    @testset "dvhat" begin
        try
            furnsh(path(CORE, :lsk), path(CORE, :pck), path(CORE, :spk))
            et = str2et("Jan 1, 2009")
            state, ltime = spkezr("Sun", et, "J2000", "LT+S", "Earth")
            x_new = dvhat(state)
            expected = [0.1834466376334262, -0.9019196633282948, -0.39100927360200305,
                        2.0244976750658316e-07, 3.4660106111045445e-08, 1.5033141925267006e-08]
            @test expected ≈ x_new
        finally
            kclear()
        end
    end
    @testset "dvnorm" begin
        mag = [-4.0, 4, 12]
        x = [1.0, sqrt(2.0), sqrt(3.0)]
        s1 = [x .* 10.0^mag[1]; x]
        s2 = [x .* 10.0^mag[2]; -x]
        s3 = [zeros(3); x .* 10.0^mag[3]]
        @test dvnorm(s1) ≈ 2.4494897 rtol=1e-6
        @test dvnorm(s2) ≈ -2.4494897 rtol=1e-6
        @test dvnorm(s3) ≈ 0.0
    end
    @testset "dvpool" begin
        try
            pdpool("DTEST_VAL", [3.1415, 186.0, 282.397])
            @test expool("DTEST_VAL")
            dvpool("DTEST_VAL")
            @test !expool("DTEST_VAL")
        finally
            clpool()
            kclear()
        end
    end
    @testset "dvsep" begin
        try
            furnsh(path(CORE, :lsk), path(CORE, :pck), path(CORE, :spk))
            et = str2et("JAN 1 2009")
            state_e, eltime = spkezr("EARTH", et, "J2000", "NONE", "SUN")
            state_m, mltime = spkezr("MOON", et, "J2000", "NONE", "SUN")
            dsept = dvsep(state_e, state_m)
            @test dsept ≈ 3.8121194e-09 rtol=1e-6
        finally
            kclear()
        end
    end
end


text:
module NavierStokes2D

using SparseArrays, LinearAlgebra
using Future: copy!

include("grids.jl")
include("parameters.jl")
export Grid, Parameters

include("models.jl")
export DiffusionModel

include("caches.jl")
export ExplicitDiffusionCache, CNDiffusionCache, ADIDiffusionCache

include("methods_diffusion.jl")
export ExplicitDiffusion, CrankNickolsonDiffusion, ADIDiffusion
export perform_step!

include("cache_initialisers.jl")
export init_cache

include("matrix_solvers.jl")
export init_CN, init_ADI

include("boundary_conditions.jl")
export apply_BCs!

include("problems.jl")
export DiffusionProblem

include("solutions.jl")
export DiffusionSolution
export get_analytic, solve

end


text:
using SubstitutionTilings
using Test

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


text:
struct Augment{B}
    noise::B
end

struct Augmented{A, B}
    P::A
    noise::B
end

Augmented(P, B::Augment) = Augmented(P, B.noise)

Base.size(A::Augmented) = (size(A.P, 1), size(A.P, 2) + size(A.noise, 2))
Base.size(A::Augmented, d::Integer) = size(A)[d]

struct AugmentedSigmaPoints{T, V <: AbstractVector{T}, L <: LowerTriangular{T}, W <: AbstractWeightingParameters} <: AbstractSigmaPoints{T}
    x0::V
    P_chol::L
    noise_chol::L
    weight_params::W
    AugmentedSigmaPoints{T, V, L, W}(x0, P_chol, noise_chol, weight_params) where {T<:Real, V<:AbstractVector{T}, L<:LowerTriangular{T}, W<:AbstractWeightingParameters} =
        size(x0, 1) == size(P_chol, 1) == size(P_chol, 2) ?
        new{T, V, L, W}(x0, P_chol, noise_chol, weight_params) :
        error("The length of the first dimension must be equal to the size of P_chol")
end

Base.size(S::AugmentedSigmaPoints) = (length(S.x0) + size(S.noise_chol, 1), 2 * size(S.P_chol, 2) + 2 * size(S.noise_chol, 2) + 1)

Base.getindex(S::AugmentedSigmaPoints{T}, inds::Vararg{Int,2}) where {T} =
    @inbounds if inds[2] == 1
        inds[1] <= length(S.x0) ? S.x0[inds[1]] : zero(T)
    elseif 1 < inds[2] <= size(S.P_chol, 2) + 1
        inds[1] <= length(S.x0) ? S.x0[inds[1]] + S.P_chol[inds[1], inds[2] - 1] : zero(T)
    elseif size(S.P_chol, 2) + 1 < inds[2] <= size(S.P_chol, 2) + size(S.noise_chol, 2) + 1
        inds[1] <= length(S.x0) ? S.x0[inds[1]] : S.noise_chol[inds[1] - length(S.x0), inds[2] - size(S.P_chol, 2) - 1]
    elseif size(S.P_chol, 2) + size(S.noise_chol, 2) + 1 < inds[2] <= 2 * size(S.P_chol, 2) + size(S.noise_chol, 2) + 1
        inds[1] <= length(S.x0) ? S.x0[inds[1]] - S.P_chol[inds[1], inds[2] - size(S.P_chol, 2) - size(S.noise_chol, 2) - 1] : zero(T)
    else
        inds[1] <= length(S.x0) ? S.x0[inds[1]] : -S.noise_chol[inds[1] - length(S.x0), inds[2] - 2 * size(S.P_chol, 2) - size(S.noise_chol, 2) - 1]
    end

AugmentedSigmaPoints(x0::V, P_chol::L, noise_chol::L, weight_params::W) where {T<:Real, V<:AbstractVector{T}, L<:LowerTriangular{T}, W<:AbstractWeightingParameters} =
    AugmentedSigmaPoints{T, V, L, W}(x0, P_chol, noise_chol, weight_params)
AugmentedSigmaPoints(x0::V, P_chol::L, noise_chol::Cholesky{T}, weight_params::W) where {T<:Real, V<:AbstractVector{T}, L<:LowerTriangular{T}, W<:AbstractWeightingParameters} =
    AugmentedSigmaPoints{T, V, L, W}(x0, P_chol, noise_chol.L, weight_params)
AugmentedSigmaPoints(x0::V, P_chol::Cholesky{T}, noise_chol::L, weight_params::W) where {T<:Real, V<:AbstractVector{T}, L<:LowerTriangular{T}, W<:AbstractWeightingParameters} =
    AugmentedSigmaPoints{T, V, L, W}(x0, P_chol.L, noise_chol, weight_params)
AugmentedSigmaPoints(x0::V, P_chol::Cholesky{T}, noise_chol::Cholesky{T}, weight_params::W) where {T<:Real, V<:AbstractVector{T}, W<:AbstractWeightingParameters} =
    AugmentedSigmaPoints(x0, P_chol.L, noise_chol.L, weight_params)

function calc_sigma_points(
    x::AbstractVector{T},
    P::Augmented{<:AbstractMatrix{T}, <:AbstractMatrix{T}},
    weight_params::W
) where {T, W<:AbstractWeightingParameters}
    weight = calc_cholesky_weight(weight_params, P)
    P_chol = cholesky(Hermitian(P.P .* weight, :L))
    noise_chol = cholesky(Hermitian(P.noise .* weight, :L))
    AugmentedSigmaPoints(x, P_chol.L, noise_chol.L, weight_params)
end

function calc_sigma_points(
    x::AbstractVector{T},
    P::Augmented{<:Cholesky{T}, <:Cholesky{T}},
    weight_params::W
) where {T, W<:AbstractWeightingParameters}
    weight = calc_cholesky_weight(weight_params, P)
    AugmentedSigmaPoints(x, P.P.L * sqrt(weight), P.noise.L * sqrt(weight), weight_params)
end

function calc_sigma_points!(
    P_chol_temp::Augmented{<:AbstractMatrix{T}, <:AbstractMatrix{T}},
    x::AbstractVector{T},
    P::Augmented{<:AbstractMatrix{T}, <:AbstractMatrix{T}},
    weight_params::W
) where {T, W<:AbstractWeightingParameters}
    weight = calc_cholesky_weight(weight_params, P)
    P_chol_temp.P .= P.P .* weight
    P_chol = cholesky!(Hermitian(P_chol_temp.P, :L))
    P_chol_temp.noise .= P.noise .* weight
    P_chol_noise = cholesky!(Hermitian(P_chol_temp.noise, :L))
    AugmentedSigmaPoints(x, P_chol.L, P_chol_noise.L, weight_params)
end

function calc_sigma_points!(
    P_chol_temp::Augmented{<:AbstractMatrix{T}, <:AbstractMatrix{T}},
    x::AbstractVector{T},
    P::Augmented{<:Cholesky{T}, <:Cholesky{T}},
    weight_params::W
) where {T, W<:AbstractWeightingParameters}
    weight = calc_cholesky_weight(weight_params, P)
    P_chol_temp.P .= (P.P.uplo === 'U' ? transpose(P.P.U) : P.P.L) .* sqrt(weight)
    P_chol_temp.noise .= (P.noise.uplo === 'U' ? transpose(P.noise.U) : P.noise.L) .* sqrt(weight)
    AugmentedSigmaPoints(x, LowerTriangular(P_chol_temp.P), LowerTriangular(P_chol_temp.noise), weight_params)
end

function transform(F, χ::AugmentedSigmaPoints{T}) where T
    𝓨_x0 = F(χ.x0)
    num_x = length(χ.x0)
    𝓨_xi = Matrix{T}(undef, length(𝓨_x0), 2 * size(χ.P_chol, 2) + 2 * size(χ.noise_chol, 2))
    xi_temp = copy(χ.x0)
    @inbounds for i = size(χ.P_chol, 2):-1:1
        xi_temp[i:num_x] .= @view(χ.x0[i:num_x]) .+ @view(χ.P_chol.data[i:num_x, i])
        𝓨_xi[:, i] = F(xi_temp)
        xi_temp[i:num_x] .= @view(χ.x0[i:num_x]) .- @view(χ.P_chol.data[i:num_x, i])
        𝓨_xi[:, i + size(χ.P_chol, 2) + size(χ.noise_chol, 2)] = F(xi_temp)
    end
    @inbounds for i = 1:size(χ.noise_chol, 2)
        𝓨_xi[:, i + size(χ.P_chol, 2)] = F(χ.x0, @view(χ.noise_chol[:, i]))
        𝓨_xi[:, i + 2 * size(χ.P_chol, 2) + size(χ.noise_chol, 2)] = F(χ.x0, -@view(χ.noise_chol[:, i]))
    end
    TransformedSigmaPoints(𝓨_x0, 𝓨_xi, χ.weight_params)
end

function transform!(𝓨::TransformedSigmaPoints{T}, xi_temp, F!, χ::AugmentedSigmaPoints{T}) where T
    F!(𝓨.x0, χ.x0)
    num_x = length(χ.x0)
    xi_temp.P .= χ.x0
    @inbounds for i = size(χ.P_chol, 2):-1:1
        xi_temp.P[i:num_x] .= @view(χ.x0[i:num_x]) .+ @view(χ.P_chol.data[i:num_x, i])
        F!(@view(𝓨.xi[:, i]), xi_temp.P)
        xi_temp.P[i:num_x] .= @view(χ.x0[i:num_x]) .- @view(χ.P_chol.data[i:num_x, i])
        F!(@view(𝓨.xi[:, i + size(χ.P_chol, 2) + size(χ.noise_chol, 2)]), xi_temp.P)
    end
    @inbounds for i = 1:size(χ.noise_chol, 2)
        F!(@view(𝓨.xi[:, i + size(χ.P_chol, 2)]), χ.x0, @view(χ.noise_chol[:, i]))
        xi_temp.noise[:] .= -1 .* @view(χ.noise_chol[:, i])
        F!(@view(𝓨.xi[:, i + 2 * size(χ.P_chol, 2) + size(χ.noise_chol, 2)]), χ.x0, xi_temp.noise)
    end
    TransformedSigmaPoints(𝓨.x0, 𝓨.xi, χ.weight_params)
end

function cov(χ::AugmentedSigmaPoints, unbiased_𝓨::TransformedSigmaPoints)
    weight_0, weight_i = calc_cov_weights(χ)
    num_states = length(χ.x0)
    num_noise_states = size(χ.noise_chol, 2)
    (χ.P_chol * (@view(unbiased_𝓨.xi[:, 1:num_states]))' .-
        χ.P_chol * (@view(unbiased_𝓨.xi[:, num_states + num_noise_states + 1:2 * num_states + num_noise_states]))') .* weight_i
end

function cov!(P, χ::AugmentedSigmaPoints, unbiased_𝓨::TransformedSigmaPoints)
    weight_0, weight_i = calc_cov_weights(χ)
    num_states = length(χ.x0)
    num_noise_states = size(χ.noise_chol, 2)
    P .= @~ χ.P_chol * (@view(unbiased_𝓨.xi[:, 1:num_states]))'
    P .-= @~ χ.P_chol * (@view(unbiased_𝓨.xi[:, num_states + num_noise_states + 1:2 * num_states + num_noise_states]))'
    P .*= weight_i
end

function cov(unbiased_𝓨::TransformedSigmaPoints, Q::Augment)
    weight_0, weight_i = calc_cov_weights(unbiased_𝓨)
    unbiased_𝓨.x0 * unbiased_𝓨.x0' .* weight_0 .+ unbiased_𝓨.xi * unbiased_𝓨.xi' .* weight_i
end

function cov!(P, unbiased_𝓨::TransformedSigmaPoints, Q::Augment)
    weight_0, weight_i = calc_cov_weights(unbiased_𝓨)
    P .= @~ unbiased_𝓨.x0 * unbiased_𝓨.x0' .* weight_0 .+ unbiased_𝓨.xi * unbiased_𝓨.xi' .* weight_i
end


text:
import SequentialMonteCarlo.SMCModel
import SimpleSMC.SMCIO

struct CCSMCIO{Particle, ParticleScratch}
  N::Int64
  n::Int64
  smcio1::SMCIO{Particle, ParticleScratch}
  smcio2::SMCIO{Particle, ParticleScratch}
  ref1::Vector{Particle}
  ref2::Vector{Particle}
  # normalized weights
  nws1::Vector{Float64}
  nws2::Vector{Float64}
  # common weights & probability
  cws::Vector{Float64}
  # residual weights
  rws1::Vector{Float64}
  rws2::Vector{Float64}
  # backward sampling weights
  bws1::Vector{Float64}
  bws2::Vector{Float64}

  extraRNG::RNG # for rng coupled sampling
  idxs::Vector{Int64} # for JLS
  σ::Vector{Int64} # for JLS
end

function CCSMCIO{Particle, ParticleScratch}(N::Int64, n::Int64) where
  {Particle, ParticleScratch}
  smcio1 = SMCIO{Particle, ParticleScratch}(N, n)
  smcio2 = SMCIO{Particle, ParticleScratch}(N, n)
  ref1 = Vector{Particle}(undef, n)
  ref2 = Vector{Particle}(undef, n)
  for i in 1:n
    ref1[i] = Particle()
    ref2[i] = Particle()
  end
  nws1 = Vector{Float64}(undef, N)
  nws2 = Vector{Float64}(undef, N)
  cws = Vector{Float64}(undef, N)
  rws1 = Vector{Float64}(undef, N)
  rws2 = Vector{Float64}(undef, N)
  bws1 = Vector{Float64}(undef, N)
  bws2 = Vector{Float64}(undef, N)

  extraRNG = deepcopy(getRNG())
  idxs = Vector{Int64}(undef, N)
  σ = Vector{Int64}(undef, N-1)

  return CCSMCIO{Particle, ParticleScratch}(N, n, smcio1, smcio2, ref1, ref2,
    nws1, nws2, cws, rws1, rws2, bws1, bws2, extraRNG, idxs, σ)
end


text:
export pileup

using Base.Collections

type Pileuper{T}
    reads::T
    window::PriorityQueue{Read, Int32, Base.Order.ForwardOrdering} # Read -> ref pos of last matching base
    muts::PriorityQueue{Mut, Int32, Base.Order.ForwardOrdering}    # Mut  -> pos
    chr::Int32
    Pileuper(x::T) = new(x, PriorityQueue(Read, Int32), PriorityQueue(Mut, Int32), -2)
end

pileup{T}(x::T) = pileup(Pileuper{T}(x))

function pileup{T}(p::Pileuper{T})
    for r in p.reads
        if r.refID != p.chr
            flush_muts!(p)
            p.chr = r.refID
        end
        add_muts!(p, r)
    end
    flush_muts!(p)
end

function produce_muts!{T}(p::Pileuper{T})
    mut = dequeue!(p.muts)

    while !isempty(p.window) && peek(p.window).second < mut.pos
        dequeue!(p.window)
    end

    produce((keys(p.window), p.chr, mut))
end

function flush_muts!{T}(p::Pileuper{T})
    while !isempty(p.muts)
        produce_muts!(p)
    end

    # clear p.window by hack; not sure if this is faster than just allocating a new one
    empty!(p.window.xs)
    empty!(p.window.index)
end

function add_muts!{T}(p::Pileuper{T}, r::Read)
    while !isempty(p.muts) && peek(p.muts).second < r.pos
        produce_muts!(p)
    end

    enqueue!(p.window, r, r.pos + calc_distance(r) - 1)

    for mut in r.muts
        mut = map_to_ref(mut, r)
        haskey(p.muts, mut) || enqueue!(p.muts, mut, mut.pos)
    end
end


text:
"""
    Triplets(D::EDM{T}) where T

Compute the triplets from D.
"""
function TripletEmbeddings.Triplets(D::EDM{T}) where T <: Real

    n = size(D,1)
    triplets = Vector{Tuple{Int,Int,Int}}(undef, n*binomial(n-1, 2))
    counter = 0

    for k = 1:n, j = 1:k-1, i = 1:n
        if i != j && i != k
            if D[i,j] < D[i,k]
                counter +=1
                @inbounds triplets[counter] = (i, j, k)
            elseif D[i,j] > D[i,k]
                counter += 1
                @inbounds triplets[counter] = (i, k, j)
            end
        end
    end

    return Triplets(triplets[1:counter])
end

"""
    Triplets(D::MaskedEDM{T}) where T

Compute the triplets from D, using only its unmasked (or non-zero) values.
"""
function TripletEmbeddings.Triplets(D::MaskedEDM{T}) where T <: Real

    n = size(D,1)
    Dm = masked(D)
    triplets = Vector{Tuple{Int,Int,Int}}(undef, n*binomial(n-1, 2))
    counter = 0

    for k = 1:n, j = 1:k-1, i = 1:n
        if i != j && i != k
            if Dm[i,j] != 0 && Dm[i,k] != 0
                if Dm[i,j] < Dm[i,k]
                    counter +=1
                    @inbounds triplets[counter] = (i, j, k)
                elseif Dm[i,j] > Dm[i,k]
                    counter += 1
                    @inbounds triplets[counter] = (i, k, j)
                end
            end
        end
    end

    return Triplets(triplets[1:counter])
end


text:
module TakingBroadcastSeriously
import Base: ^
using Base.Broadcast: broadcast_c, containertype

broadcast_(f, A, Bs...) = broadcast_c(f, containertype(A, Bs...), A, Bs...)

struct Broadcasted{T}
  x::T
end

unwrap(w::Broadcasted) = w.x

# We must hack each function we want to use with un-fused broadcasting.
for f in :[sin, cos, +, -, *, /, %, ^, <, <=, >, >=, !=, ==].args
  @eval Base.$f(a::Broadcasted...) = Broadcasted(broadcast_($f, unwrap.(a)...))

  #sometimes literals "resist" wrapping, so catch them too
  @eval Base.$f(a::Broadcasted, b) = Broadcasted(broadcast_($f, unwrap(a), b))
  @eval Base.$f(b, a::Broadcasted) = Broadcasted(broadcast_($f, b, unwrap(a)))
end

#Avoid ambiguitity with Base
^(a::Broadcasted, b::Integer) = Broadcasted(broadcast_(^, unwrap(a), b))

macro unfuse(T)
  T = esc(T)
  quote
    Base.broadcast(f, A::$T, Bs...) = f(Broadcasted(A), Broadcasted.(Bs)...) |> unwrap
    Base.broadcast(f, A, B::$T, Cs...) = f(Broadcasted(A), Broadcasted(B), Broadcasted.(Cs)...) |> unwrap
    Base.broadcast(f, A::$T, B::$T, Cs...) = f(Broadcasted(A), Broadcasted(B), Broadcasted.(Cs)...) |> unwrap
  end
end

end # module


text:
cum_sum = 0

for x1 = 0:1
    for x2 = 0:1
        for x3 = 0:1
            global cum_sum
            # cum_sum += .25*((i == 1 && j == 1)*(.5*(k == 1) + .1*(k != 1)) +
            # 1/6*(i != 1))
            cum_sum += (1/4)^(2)*(2.0)^(x2*x1)*(3.0)^(x3*x1*(1-x2-x3)+x3)
        end
    end
end

println(cum_sum)

max_stars_repo_name: seanrattana/courses

max_issues_repo_name: seanrattana/courses

max_forks_repo_name: seanrattana/courses


text:
abstract type AbstractTurbineDefinition end

struct TurbineDefinition{TI,AF,CTM,PM} <: AbstractTurbineDefinition
    id::TI
    rotor_diameter::AF
    hub_height::AF
    cut_in_speed::AF
    rated_speed::AF
    cut_out_speed::AF
    rated_power::AF
    generator_efficiency::AF
    ct_model::CTM
    power_model::PM
end