text:
include("Self_recognition-Compartments.jl")

function simulate_self_recognition(l,n,runin,steps,β,βeq,equil_steps,γ_steps,compartments,x_max,xsteps,reps)

    println("test1")
    filename = "../Data/Workingmemory_compartment/L_$(l)/reps_$(reps)/Scaled_temp_compartments_$(compartments)-SelfRec_Xsteps_$(xsteps)-temp_$(β)-elqutemp_$(βeq)_Xmax_$(x_max)_runin_$(runin)_eqsteps_$(equil_steps)_Patterns_$(n)_lowx.jld2"

    println("I want to do  $(filename)")
    if isfile(filename)
        println("Die daten $(filename) exestieren schon ich hoehre auf")
    else
        if !isdir("../Data")
            mkdir("../Data")
        end

        if !isdir("../Data/Workingmemory_compartment")
            mkdir("../Data/Workingmemory_compartment")
        end

        if !isdir("../Data/Workingmemory_compartment/L_$(l)")
            mkdir("../Data/Workingmemory_compartment/L_$(l)")
        end
        if !isdir("../Data/Workingmemory_compartment/L_$(l)/reps_$(reps)")
            mkdir("../Data/Workingmemory_compartment/L_$(l)/reps_$(reps)")
        end
        println("I'm starting the simulation")
        β_scaled=10^(β/10)/(1 /(l/2/n*compartments -sqrt(l/π/2))*(compartments*l)^2 /n/n)
        β_eq_scaled=10^(βeq/10)*n/compartments
        if (compartments==1)
            println("only one Compartment")
            e_rac,e_rand,overlapsload,γ,x = build_rec_for_gamma_x_One(l,n,runin,steps,β_eq_scaled,equil_steps,γ_steps,x_max,xsteps,reps)

            overlaps=zeros(γ_steps,xsteps)
            notrecognized=zeros(γ_steps,xsteps)
            optG=zeros(xsteps)
            min_gamma=zeros(xsteps)
            max_gamma=zeros(xsteps)
            for xx in 1:xsteps
                for gg in 1:γ_steps
                    overlaps[gg,xx] =mean(overlapsload[:,gg,xx])
                    for attempt in 1:length(overlapsload[:,gg,xx])
                            if (overlapsload[attempt,gg,xx]/l <0.8)
                                notrecognized[gg,xx]+=1
                            end
                    end
                end

                errorregime=sqrt(max(1,2*notrecognized[argmax(overlaps[:,xx]),xx]))
                threshEqual= 1/(2*n*reps)*errorregime
                optG[xx],min_gamma[xx],max_gamma[xx],dummi= opt_lambda_analysis(overlaps[:,xx],γ[:,xx],threshEqual)
            end

            parameters=(l,n,runin,steps,β,βeq,equil_steps,γ_steps,compartments,x_max,xsteps,reps)
            save(filename, "overlaps",overlaps, "gamma",γ,"x",x,"parameters",parameters,"gammaopt",[optG,min_gamma,max_gamma],"notrecognized",notrecognized) ## Here we don't save the energyies to minimize the data size
            #save(filename,"e_rac",e_rac, "e_rand",e_rand, "overlaps",overlaps, "gamma",γ,"x",x,"parameters",parameters)

        else
            if (compartments > n)
                println("to many compartments" )
                return -1
            end
            println("you used $(compartments) compartments for $(n) patterns" )
            e_rac,e_rand,overlapsload,γ,x = build_rec_for_gamma_x_compartment(l,n,runin,steps,β_scaled,β_eq_scaled,equil_steps,γ_steps,compartments,x_max,xsteps,reps)

            overlaps=zeros(γ_steps,xsteps)
            notrecognized=zeros(γ_steps,xsteps)
            optG=zeros(xsteps)
            min_gamma=zeros(xsteps)
            max_gamma=zeros(xsteps)
            for xx in 1:xsteps
                for gg in 1:γ_steps
                    overlaps[gg,xx] =mean(overlapsload[:,gg,xx])
                    for attempt in 1:length(overlapsload[:,gg,xx])
                            if (overlapsload[attempt,gg,xx]/l <0.8)
                                notrecognized[gg,xx]+=1
                            end
                    end
                end

                errorregime=sqrt(max(1,2*notrecognized[argmax(overlaps[:,xx]),xx]))
                threshEqual= 1/(2*n*reps)*errorregime
                optG[xx],min_gamma[xx],max_gamma[xx],dummi= opt_lambda_analysis(overlaps[:,xx],γ[:,xx],threshEqual)
            end

            parameters=(l,n,runin,steps,β,βeq,equil_steps,γ_steps,compartments,x_max,xsteps,reps)
            save(filename, "overlaps",overlaps, "gamma",γ,"x",x,"parameters",parameters,"gammaopt",[optG,min_gamma,max_gamma],"notrecognized",notrecognized) ## Here we don't save the energyies to minimize the data size
            #save(filename,"e_rac",e_rac, "e_rand",e_rand, "overlaps",overlaps, "gamma",γ,"x",x,"parameters",parameters)

        end

        println("I'm saving the file $(filename)")
    end
end

function simulate_self_recognition_alongopt(l,n,runin,steps,β,βeq,equil_steps,γ_steps,compartments,x_max,xsteps,reps)

    filename = "../Data/Workingmemory_compartment_opt/L_$(l)/reps_$(reps)/Optimal-compartments_$(compartments)-SelfRec_Xsteps_$(xsteps)-temp_$(β)-elqutemp_$(βeq)_Xmax_$(x_max)_runin_$(runin)_eqsteps_$(equil_steps)_Patterns_$(n)_lowx.jld2"

    println("I want to do  $(filename)")
    if isfile(filename)
        println("Die daten $(filename) exestieren schon ich hoehre auf")
    else

        if !isdir("../Data/Workingmemory_compartment_opt")
            mkdir("../Data/Workingmemory_compartment_opt")
        end

        if !isdir("../Data/Workingmemory_compartment_opt/L_$(l)")
            mkdir("../Data/Workingmemory_compartment_opt/L_$(l)")
        end
        if !isdir("../Data/Workingmemory_compartment_opt/L_$(l)/reps_$(reps)")
            mkdir("../Data/Workingmemory_compartment_opt/L_$(l)/reps_$(reps)")
        end
        println("Ich starte mit der siumulation")

        filenameload = "../Data/Workingmemory_compartment/L_$(l)/reps_$(reps)/Scaled_temp_compartments_$(compartments)-SelfRec_Xsteps_$(xsteps)-temp_$(β)-elqutemp_$(βeq)_Xmax_$(x_max)_runin_$(runin)_eqsteps_$(equil_steps)_Patterns_$(n)_lowx.jld2"

        β_scaled=10^(β/10)/(1 /(l/2/n*compartments -sqrt(l/π/2))*(compartments*l)^2 /n/n)
        β_eq_scaled=10^(βeq/10)*n/compartments
        if (compartments==1)
            println("only one compartment")
            optg1,dumi1,dumi2=load(filenameload,"gammaopt")

            x=load(filenameload,"x")

            e_rac,e_rand,meane,std_low,std_high,overlaps = build_rec_for_opt_gamma_One(l,n,x,optg1,runin,steps,β_eq_scaled,equil_steps*20,reps)
            parameters=(l,n,runin,steps,β,βeq,equil_steps,γ_steps,compartments,x_max,xsteps,reps)
            save(filename,"e_rac",e_rac, "e_rand",e_rand, "overlaps",overlaps, "gamma",optg1,"x",x,"meanE",meane,"std_low",std_low,"std_high",std_high,"parameters",parameters)

        else
            if (compartments > n)
                println("to many compartments" )
                return -1
            end
            overlaps=load(filenameload,"overlaps")

            x=load(filenameload,"x")

            optg1,dumi1,dumi2=load(filenameload,"gammaopt")
            pattern_used,compartment_used,e_rac,e_rand,overlaps = build_rec_for_opt_gamma_compartment(l,n,x,optg1,runin,steps,β_scaled,β_eq_scaled,equil_steps*20,reps,compartments)
            parameters=(l,n,runin,steps,β,βeq,equil_steps,γ_steps,compartments,x_max,xsteps,reps)
            save(filename,"e_rac",e_rac, "e_rand",e_rand, "overlaps",overlaps, "gamma",optg1,"x",x,"pattern_used",pattern_used,"compartment_used",compartment_used,"parameters",parameters)

        end

        println("I'm saving the file $(filename)")
    end
end

function simulate_self_recognition_alongopt_no_evolution(l,n,runin,steps,β,βeq,equil_steps,γ_steps,γ_min,compartments,x_max,xsteps,reps)

    println("test1")
    #filename = "../Data/Workingmemory_compartment/L_$(l)/reps_$(reps)/compartments_$(compartments)-SelfRec_Xsteps_$(xsteps)-temp_$(β)-eqtemp-$(βeq)_Xmax_$(x_max)_runin_$(runin)_eqsteps_$(equil_steps)_Patterns_$(n).jld2"
    filename = "../Data/Workingmemory_compartment_opt/L_$(l)/reps_$(reps)/Scaled_temp_compartments_$(compartments)-SelfRec_Xsteps_$(xsteps)-temp_$(β)-elqutemp_$(βeq)_Xmax_$(x_max)_runin_$(runin)_eqsteps_$(equil_steps)_Patterns_$(n)_small_no_evolution.jld2"

    println("I want to do  $(filename)")
    if isfile(filename)
        println("Die daten $(filename) exestieren schon ich hoehre auf")
    else

        if !isdir("../Data/Workingmemory_compartment_opt")
            mkdir("../Data/Workingmemory_compartment_opt")
        end

        if !isdir("../Data/Workingmemory_compartment_opt/L_$(l)")
            mkdir("../Data/Workingmemory_compartment_singleX_opt/L_$(l)")
        end
        if !isdir("../Data/Workingmemory_compartment_pt/L_$(l)/reps_$(reps)")
            mkdir("../Data/Workingmemory_compartment_/L_$(l)/reps_$(reps)")
        end
        println("Ich starte mit der siumulation")

        filenameload = "../Data/Workingmemory_compartment/L_$(l)/reps_$(reps)/Scaled_temp_compartments_$(compartments)-SelfRec_Xsteps_$(xsteps)-temp_$(β)-elqutemp_$(βeq)_Xmax_$(x_max)_runin_$(runin)_eqsteps_$(equil_steps)_Patterns_$(n)_lowx.jld2"

        β_scaled=10^(β/10)/(1 /(l/2/n*compartments -sqrt(l/π/2))*(compartments*l)^2 /n/n)
        β_eq_scaled=10^(βeq/10)*n/compartments
        if (compartments==1)
            println("only one compartment")
            optg1,dumi1,dumi2=load(filenameload,"gammaopt")

            x=load(filenameload,"x")

            e_rac,e_rand,meane,std_low,std_high,overlaps = build_rec_for_opt_gamma_One(l,n,0,optg1,runin,steps,β_eq_scaled,equil_steps*20,reps)
            parameters=(l,n,runin,steps,β,βeq,equil_steps,γ_steps,compartments,x_max,xsteps,reps)
            save(filename,"e_rac",e_rac, "e_rand",e_rand, "overlaps",overlaps, "gamma",optg1,"x",x,"meanE",meane,"std_low",std_low,"std_high",std_high,"parameters",parameters)

        else
            if (compartments > n)
                println("to many compartments" )
                return -1
            end
            overlaps=load(filenameload,"overlaps")

            x=load(filenameload,"x")

            optg1,dumi1,dumi2=load(filenameload,"gammaopt")
            pattern_used,compartment_used,e_rac,e_rand,overlaps = build_rec_for_opt_gamma_compartment(l,n,0,optg1,runin,steps,β_scaled,β_eq_scaled,equil_steps*20,reps,compartments)
            parameters=(l,n,runin,steps,β,βeq,equil_steps,γ_steps,compartments,x_max,xsteps,reps)
            save(filename,"e_rac",e_rac, "e_rand",e_rand, "overlaps",overlaps, "gamma",optg1,"x",x,"pattern_used",pattern_used,"compartment_used",compartment_used,"parameters",parameters)

        end

        println("I'm saving the file $(filename)")
    end
end

function simulate_self_recognition_singlex_alongopts_stats(l,n,runin,steps,β,βeq,equil_steps,γ_steps,γ_min,compartments,x_max,xsteps,reps,this_x)

    println("test1")
        filename = "../Data/Workingmemory_compartment_singleX_opt/L_$(l)/reps_$(reps)/Scaled_temp_compartments_$(compartments)-SelfRec_Xsteps_$(xsteps)-temp_$(β)-elqutemp_$(βeq)_Xmax_$(x_max)_runin_$(runin)_eqsteps_$(equil_steps)_Patterns_$(n)_thisX-$(this_x)_small_stats.jld2"

    println("I want to do  $(filename)")
    if isfile(filename)
        println("Die daten $(filename) exestieren schon ich hoehre auf")
    else

        if !isdir("../Data/Workingmemory_compartment_singleX_opt")
            mkdir("../Data/Workingmemory_compartment_singleX_opt")
        end

        if !isdir("../Data/Workingmemory_compartment_singleX_opt/L_$(l)")
            mkdir("../Data/Workingmemory_compartment_singleX_opt/L_$(l)")
        end
        if !isdir("../Data/Workingmemory_compartment_singleX_opt/L_$(l)/reps_$(reps)")
            mkdir("../Data/Workingmemory_compartment_singleX_opt/L_$(l)/reps_$(reps)")
        end
        println("Ich starte mit der siumulation")

        filenameload = "../Data/Workingmemory_compartment_singleX/L_$(l)/reps_50/Scaled_temp_compartments_$(compartments)-SelfRec_Xsteps_$(xsteps)-temp_$(β)-elqutemp_$(βeq)_Xmax_$(x_max)_runin_$(runin)_eqsteps_$(equil_steps)_Patterns_$(n)_thisX-$(this_x)_small.jld2"
        println("Ich lade datei:")
        println(filenameload)
        β_scaled=10^(β/10)/(50)
        β_eq_scaled=10^(βeq/10)*n/compartments

        if (compartments==1)
            println("only one compartment")
            overlaps=load(filenameload,"overlaps")
            γ=load(filenameload,"gamma")
            x=load(filenameload,"x")

            #m,s,opt,optg1,opt2,optg2=give_stats_Overlaps(overlaps./l,γ);
            optg1=γ[argmax(overlaps)]
            runin_e=max.(10*n,min.(floor.(Int64,log(0.00001)./ log.(optg1) .+1),runin))
            steps_e=min.(runin_e,steps)

            dist_start,dist_next,open_paths,     dist_start_rand,dist_next_rand,open_paths_rand,   partision_mean_vecs,partision_std_vecs,partision_first_vecs,partision_last_vecs,     partision_mean_patts,partision_std_patts,partision_first_patts,partision_last_patts,best_overlaps,mean_current_direction,best_next_direciton = get_statistics_One_allstats(l,n,x[this_x],optg1[this_x],runin_e,steps_e,reps,equil_steps*20,β_eq_scaled)
            parameters=(l,n,runin,steps,β,βeq,equil_steps,γ_steps,γ_min,compartments,x_max,xsteps,reps,this_x)
            save(filename,"dist_start",dist_start, "dist_next",dist_next,"open_paths" ,open_paths,   "dist_start_rand",dist_start_rand,"dist_next_rand",dist_next_rand,"open_paths_rand",open_paths_rand,   "partision_mean_vecs",partision_mean_vecs,"partision_std_vecs",partision_std_vecs,"partision_first_vecs",partision_first_vecs,"partision_last_vecs",partision_last_vecs,     "partision_mean_patts",partision_mean_patts,"partision_std_patts",partision_std_patts,"partision_first_patts",partision_first_patts,"partision_last_patts",partision_last_patts,"best_overlaps",best_overlaps,"mean_current_direction",mean_current_direction,"best_next_direciton",best_next_direciton      ,"gamma",optg1,"x",x,"parameters",parameters)

        else
            println("sorry this only works for one compartment at this point try again another time")
        end

        println("I'm saving the file $(filename)")
    end
end

function Make_eigenvalue_eigenvectors_stats_along_x(l,n,runin,steps,equil_steps,β,βeq,compartments,x_max,xsteps,reps,this_x)

    println("test1")
    filename = "../Data/Workingmemory_compartment_singleX_opt/EV_stats/ev_stats_L-$(l)_N-$(n)_xsteps$(xsteps)_reps=$(reps)_thisX_$(this_x).jld2"

    println("I want to do  $(filename)")
    if isfile(filename)
        println("Die daten $(filename) exestieren schon ich hoehre auf")
    else

        if !isdir("../Data/Workingmemory_compartment_singleX_opt")
            mkdir("../Data/Workingmemory_compartment_singleX_opt")
        end

        if !isdir("../Data/Workingmemory_compartment_singleX_opt/EV_stats")
            mkdir("../Data/Workingmemory_compartment_singleX_opt/EV_stats")
        end

            println("Ich starte mit der siumulation")

            filenameload = "../Data/Workingmemory_compartment_singleX/L_$(l)/reps_50/Scaled_temp_compartments_$(compartments)-SelfRec_Xsteps_$(xsteps)-temp_$(β)-elqutemp_$(βeq)_Xmax_$(x_max)_runin_$(runin)_eqsteps_$(equil_steps)_Patterns_$(n)_thisX-$(this_x)_small.jld2"
            println("Ich lade datei:")
            println(filenameload)
            β_scaled=10^(β/10)/(1 /(l/2/n*compartments -sqrt(l/π/2))*(compartments*l)^2 /n/n)
            β_eq_scaled=10^(βeq/10)*n/compartments

            if (compartments==1)
                println("only one compartment")
                overlaps=load(filenameload,"overlaps")
                γ=load(filenameload,"gamma")
                x=load(filenameload,"x")

                #m,s,opt,optg1,opt2,optg2=give_stats_Overlaps(overlaps./l,γ);
                optg1=γ[argmax(overlaps)]
                runin_e=max.(10*n,min.(floor.(Int64,log(0.00001)./ log.(optg1) .+1),runin))

                data =    distribution_only_allEV(l,n,x[this_x],optg1,runin_e,reps)
                info=" data[1][rep][ii] is the iits eigenvalue of the data with ii in 1:n;  data[2][ii] is the mean of the evs"
                save(filename,"../Data",data,"gamma",optg1, "x",x,"this_x",this_x, "info",info)

            else
                println("sorry this only works for one compartment at this point try again a nother time")
            end

            println("I'm saving the file $(filename)")

    end

    println("Im now doing the non evolving one")
    #filename = "../Data/Workingmemory_compartment/L_$(l)/reps_$(reps)/compartments_$(compartments)-SelfRec_Xsteps_$(xsteps)-temp_$(β)-eqtemp-$(βeq)_Xmax_$(x_max)_runin_$(runin)_eqsteps_$(equil_steps)_Patterns_$(n).jld2"
    filename = "../Data/Workingmemory_compartment_singleX_opt/EV_stats/ev_stats_L-$(l)_N-$(n)_xsteps$(xsteps)_reps=$(reps)_thisX_$(this_x)_noevo.jld2"

    println("I want to do  $(filename)")
    if isfile(filename)
        println("Die daten $(filename) exestieren schon ich hoehre auf")
    else

        if !isdir("../Data/Workingmemory_compartment_singleX_opt")
            mkdir("../Data/Workingmemory_compartment_singleX_opt")
        end

        if !isdir("../Data/Workingmemory_compartment_singleX_opt/EV_stats")
            mkdir("../Data/Workingmemory_compartment_singleX_opt/EV_stats")
        end

            println("Ich starte mit der siumulation")

            filenameload = "../Data/Workingmemory_compartment_singleX/L_$(l)/reps_50/Scaled_temp_compartments_$(compartments)-SelfRec_Xsteps_$(xsteps)-temp_$(β)-elqutemp_$(βeq)_Xmax_$(x_max)_runin_$(runin)_eqsteps_$(equil_steps)_Patterns_$(n)_thisX-$(this_x)_small.jld2"
            println("Ich lade datei:")
            println(filenameload)
            β_scaled=10^(β/10)/(1 /(l/2/n*compartments -sqrt(l/π/2))*(compartments*l)^2 /n/n)
            β_eq_scaled=10^(βeq/10)*n/compartments

            if (compartments==1)
                println("only one compartment")
                overlaps=load(filenameload,"overlaps")
                γ=load(filenameload,"gamma")
                x=load(filenameload,"x")

                #m,s,opt,optg1,opt2,optg2=give_stats_Overlaps(overlaps./l,γ);
                optg1=γ[argmax(overlaps)]
                runin_e=max.(10*n,min.(floor.(Int64,log(0.00001)./ log.(optg1) .+1),runin))

                data =    distribution_only_allEV(l,n,0,optg1,runin_e,reps)
                info=" data[1][rep][ii] is the iits eigenvalue of the data with ii in 1:n;  data[2][ii] is the mean of the evs"
                save(filename,"../Data",data,"gamma",optg1, "x",x,"this_x",this_x, "info",info)

            else
                println("sorry this only works for one compartment at this point try again a nother time")
            end

            println("I'm saving the file $(filename)")

    end

end

function Make_eigenvalue_eigenvectors_stats(l,n,xsteps,xmin,xmax,runin_max,reps)

    println("test1")
    #filename = "../Data/Workingmemory_compartment/L_$(l)/reps_$(reps)/compartments_$(compartments)-SelfRec_Xsteps_$(xsteps)-temp_$(β)-eqtemp-$(βeq)_Xmax_$(x_max)_runin_$(runin)_eqsteps_$(equil_steps)_Patterns_$(n).jld2"
    filename = "../Data/Workingmemory_compartment_singleX_opt/EV_stats/ev_stats_L-$(l)_N-$(n)_xsteps$(xsteps)_reps=$(reps).jld2"

    println("I want to do  $(filename)")
    if isfile(filename)
        println("Die daten $(filename) exestieren schon ich hoehre auf")
    else

        if !isdir("../Data/Workingmemory_compartment_singleX_opt")
            mkdir("../Data/Workingmemory_compartment_singleX_opt")
        end

        if !isdir("../Data/Workingmemory_compartment_singleX_opt/EV_stats")
            mkdir("../Data/Workingmemory_compartment_singleX_opt/EV_stats")
        end

        println("Ich starte mit der siumulation")

        data =distribution_EV_for_x(l,n,xsteps,xmin,xmax,runin_max,reps)
        parameters=(l,n,xsteps,xmin,xmax,runin_max,reps)
        x=exp.(range(log(xmin),stop=log(0.01),length=4))
        if (xsteps>5)
            append!(x,range(0.015,stop=xmax,length=xsteps-4))
        end
        γ=1 .- sqrt.(8 .*x./(n.-1))
        info="../Data[g] is the data for gamma[g]; data[g][1][rep][ii] is the iits eigenvalue of the data with ii in 1:n;data[g][2][rep][ee] is the entropy stats of the data with ee= 1,2,3 entorpyevec(mean,1 ,last) , 456= same but normed 789 fror eneergies, where the mean,come with mean,std,and lower-uper std "
        save(filename,"../Data",data,"gamma",γ, "parameters",parameters,"info",info)

        println("I'm saving the file $(filename)")
    end
end


text:
using ReinforcementLearningEnvironmentAtari, Flux, GR, JLD2,
ReinforcementLearning
const withgpu = false
const game = "alien"
const gpudevice = 1

"""
AtariPreprocessor(; gpu = false, croptosquare = false,
                    cropfromrow = 34, color = false,
                    outdim = (80, croptosquare ? 80 : 105))
"""
function AtariPreprocessor(; gpu = false, croptosquare = false,
                             cropfromrow = 34, color = false,
                             outdim = (80, croptosquare ? 80 : 105))
    chain = []
    if croptosquare
        push!(chain, ImageCrop(1:160, cropfromrow:cropfromrow + 159))
    end
    if !((croptosquare && outdim == (160, 160)) || outdim == (160, 210))
        push!(chain, ImageResizeNearestNeighbour(outdim))
    end
    if !color
        push!(chain, x -> reshape(x, (outdim[1], outdim[2], 1)))
    end
    if gpu
        push!(chain, togpu)
    end
    ImagePreprocessor(color ? (3, 160, 210) : (160, 210), chain)
end

if withgpu
    using CUDAnative
    device!(gpudevice)
    using CuArrays
    const inputdtype = Float32
else
    const inputdtype = Float64
end
env = AtariEnv(game)
model = Chain(x -> x./inputdtype(255), Conv((8, 8), 4 => 16, relu, stride = (4, 4)),
                         Conv((4, 4), 16 => 32, relu, stride = (2, 2)),
                         x -> reshape(x, :, size(x, 4)),
                         Dense(2592, 256, relu),
                         Dense(256, length(env.actions)));
learner = DQN(model, opttype = x -> Flux.ADAM(x, .0001),
              loss = huberloss, doubledqn = true,
              updatetargetevery = 2500, nsteps = 10,
              updateevery = 4, replaysize = 10^6, nmarkov = 4,
              startlearningat = 200000);
x = RLSetup(learner,
            env,
            ConstantNumberSteps(2000),
            preprocessor = AtariPreprocessor(gpu=withgpu, outdim = (84, 84)),
            callbacks = [Visualize(wait = 0)])
GR.inline("mov")
beginprint("before$game.mov")
run!(x)
endprint()
learner.t = 0
x.callbacks = [Progress(5*10^2), EvaluationPerEpisode(TotalReward()),
               LinearDecreaseEpsilon(5 * 10^4, 10^6, 1, .01)];
x.stoppingcriterion = ConstantNumberSteps(4 * 10^6)
@info "start learning."
@time learn!(x)
@save "model$game.jld2" model
x.callbacks = [Visualize(wait = 0)]
x.stoppingcriterion = ConstantNumberSteps(2000)
beginprint("after$game.mov")
run!(x)
endprint()

ext: jl


text:
#Julia PyPlot Plotting

#Pkg.add("PyPlot")
using PyPlot
precipitation = [0,0,0,0,0,0,0,0,0,0,0.12,0.01,0,0,0,0.37,0,0,0,0,0.01,0,0,0,0.01,0.01,0,0.17,0.01,0.11,0.31]
date = collect(1:31)
fig = figure(1, figsize=(4, 4))
plot(date, precipitation, ".")
title("Boston Precipitation")
xlabel("May 2013")
ylabel("Precipitation")

max_stars_repo_name: PacktPublishing/Learning-Jupyter

max_issues_repo_name: PacktPublishing/Learning-Jupyter

max_forks_repo_name: PacktPublishing/Learning-Jupyter


text:
using Test, Downloads
using SDMX

const ecb_url = "https://sdw-wsrest.ecb.europa.eu/service/"
const flowref = "CBD2/"
const BS = "Q.IE.W0.67._Z._Z.A.F.._X..._Z.LE._T.EUR"
const PL = "Q.IE.W0.67._Z._Z.A.F.._X.ALL.._Z.T._T.EUR"

io = IOBuffer()
flat = Downloads.download(ecb_url*"data/"*flowref*BS*"?startPeriod=2020&dimensionAtObservation=AllDimensions&format=jsondata", io) |> take!
series = Downloads.download(ecb_url*"data/"*flowref*BS*"?startPeriod=2020&format=jsondata", io) |> take!

dt_flat = SDMX.read(flat)
dt_series = SDMX.read(series)

@test SDMX.name(dt_flat) == SDMX.name(dt_series)
@test SDMX.headers(dt_flat) == SDMX.headers(dt_series)
@test SDMX.data(dt_flat) == SDMX.data(dt_series)

max_stars_repo_name: lsaenzt/SDMX

max_issues_repo_name: lsaenzt/SDMX

max_forks_repo_name: lsaenzt/SDMX


text:
__precompile__(true)
module ADCME

    export tf,
    tfp,
    AUTO_REUSE,
    GLOBAL_VARIABLES,
    TRAINABLE_VARIABLES,
    UPDATE_OPS

    using PyCall
    using Random
    using LinearAlgebra
    using SparseArrays
    using LibGit2
    using Libdl

    tf = PyNULL()
    tfp = PyNULL()
    libadcme = missing

    #------------------------------------------------------------------------------------------
    # Global Storage
    DTYPE = Dict{Type, PyObject}()
    LIBADCME = abspath(joinpath("$(@__DIR__)", "../deps/CustomOps/build", "libadcme.$dlext"))
    if Sys.iswindows()
        LIBADCME = abspath(joinpath("$(@__DIR__)", "../deps/CustomOps/build", "adcme.dll"))
    end
    LIBPLUGIN = abspath(joinpath("$(@__DIR__)", "../deps/Plugin"))

    if isfile("$(@__DIR__)/../deps/deps.jl")
        include("$(@__DIR__)/../deps/deps.jl")
        if Sys.iswindows()
            ENV["PATH"] = LIBDIR*";"*ENV["PATH"]
        else
            if haskey(ENV, "LD_LIBRARY_PATH")
                ENV["LD_LIBRARY_PATH"] = LIBDIR*":"*ENV["LD_LIBRARY_PATH"]
            else
                ENV["LD_LIBRARY_PATH"] = LIBDIR
            end
        end
    else
        error("ADCME is not properly built; run `Pkg.build(\"ADCME\")` to fix the problem.")
    end
    run_metadata = nothing
    STORAGE = Dict{String, Any}()

    function __init__()
        global AUTO_REUSE, GLOBAL_VARIABLES, TRAINABLE_VARIABLES, UPDATE_OPS, DTYPE, libadcme
        copy!(tf, pyimport("tensorflow"))
        tf.compat.v1.logging.set_verbosity(tf.compat.v1.logging.ERROR)
        try
            copy!(tfp, pyimport("tensorflow_probability"))
        catch
        end

        DTYPE = Dict(Float64=>tf.float64,
            Float32=>tf.float32,
            Int64=>tf.int64,
            Int32=>tf.int32,
            Bool=>tf.bool,
            ComplexF64=>tf.complex128,
            ComplexF32=>tf.complex64,
            String=>tf.string,
            Char=>tf.string)
        AUTO_REUSE = tf.compat.v1.AUTO_REUSE
        GLOBAL_VARIABLES = tf.compat.v1.GraphKeys.GLOBAL_VARIABLES
        TRAINABLE_VARIABLES = tf.compat.v1.GraphKeys.TRAINABLE_VARIABLES
        UPDATE_OPS = tf.compat.v1.GraphKeys.UPDATE_OPS
        global options = Options()

        try
            libadcme = tf.load_op_library(LIBADCME)
        catch
            @warn "Cannot load $LIBADCME. Please recompile the shared library by `ADCME.precompile()` for using custom operators."
        end
    end

    include("core.jl")
    include("io.jl")
    include("optim.jl")
    include("run.jl")
    include("variable.jl")
    include("ops.jl")
    include("layers.jl")
    include("extra.jl")
    include("install.jl")
    include("sparse.jl")
    include("random.jl")
    include("gan.jl")
    include("ot.jl")
    include("ode.jl")
    include("flow.jl")
    include("options.jl")
    include("mpi.jl")
    include("toolchain.jl")
    include("kit.jl")
    include("rbf.jl")
    include("pcl.jl")
end


text:
"""
    df(sem_fit::SemFit)
    df(model::AbstractSem)

Return the degrees of freedom.
"""
function df end

df(sem_fit::SemFit) = df(sem_fit.model)

df(model::AbstractSem) = n_dp(model) - n_par(model)

function n_dp(model::AbstractSemSingle)
    nman = n_man(model)
    ndp = 0.5(nman^2 + nman)
    if !isnothing(model.imply.μ)
        ndp += n_man(model)
    end
    return ndp
end

n_dp(model::SemEnsemble) = sum(n_dp.(model.sems))


text:
using Pkg
Pkg.activate(".")
Pkg.instantiate()

using Pluto
using PDELib

Pkg.update()

notebooks=["examples/Pluto-GridsAndVisualization.jl",
           "examples/Pluto-MultECatGrids.jl"]

for notebook in notebooks
    println("Updating packages in $(notebook):")
    Pluto.activate_notebook_environment(notebook)
    Pkg.update()
    Pkg.status()
    run(`git diff $(notebook)`)
    println("Updating of  $(notebook) done\n")
    Pkg.activate(".")
end


text:
#=
core.jl
ImageGeom type, constructors, and core methods
=#

export ImageGeom
export downsample, oversample, expand_nz
export axis, grids, axisf, axesf, gridf

#using ImageGeoms: downsample, upsample
using LazyGrids: ndgrid
using FillArrays: Trues, Falses, Zeros, Ones

## Type

"""
    ImageGeom{D,S,M}

Image geometry struct with essential grid parameters.

- `dims::Dims{D}` image dimensions
- `deltas::S` where `S <: NTuple{D}` pixel sizes,
  where each Δ is usually `Real` or `Unitful.Length`
- `offsets::NTuple{D,Float32}` unitless
- `mask::M` where `M <: AbstractArray{Bool,D}` logical mask, often `FillArrays.Trues(dims)`.
"""
struct ImageGeom{D,S,M} # S <: NTuple{D,RealU}, M <: AbstractArray{Bool,D}}
    dims::Dims{D} # image dimensions
    deltas::S # pixel sizes
    offsets::NTuple{D,Float32} # unitless
    mask::M # logical mask for image support

    function ImageGeom{D,S,M}(
        dims::Dims{D},
        deltas::S,
        offsets::NTuple{D,Real},
        mask::M
    ) where {D,
            S <: NTuple{D,RealU}, # i.e., Union{Real,Unitful.Length}
            M <: AbstractArray{Bool,D},
        }
        any(<=(zero(Int)), dims) && throw("dims must be positive")
        any(iszero, deltas) && throw("deltas must be nonzero")
        size(mask) == dims ||
            throw(DimensionMismatch("mask size $(size(mask)) vs dims $dims"))
        new{D,S,M}(dims, deltas, Float32.(offsets), mask)
    end
end

## Constructors

"""
    ig = ImageGeom(dims::Dims, deltas, offsets, [, mask])
Constructor for `ImageGeom` of dimensions `dims`.
* The `deltas` elements (tuple of grid spacings)
  should each be `Real` or a `Unitful.Length`; default `(1,…,1)`.
* The `offsets` (tuple of grid offsets) must be unitless;  default `(0,…,0)`.
* The `dims`, `deltas` and `offsets` tuples must be same length.
* Default `mask` is `FillArrays.Trues(dims)` which is akin to `trues(dims)`.

# Example

```jldoctest
julia> ImageGeom((5,7), (2.,3.))
ImageGeom{2, NTuple{2,Float64}}
 dims::NTuple{2,Int64} (5, 7)
 deltas::NTuple{2,Float64} (2.0, 3.0)
 offsets::NTuple{2,Float32} (0.0f0, 0.0f0)
 mask: 5×7 Ones{Bool} {35 of 35}
```
"""
function ImageGeom(
    dims::Dims{D},
    deltas::S,
    offsets::NTuple{D,Real},
    mask::M,
) where {D, S <: NTuple{D,RealU}, M <: AbstractArray{Bool,D}}
    ImageGeom{D,S,M}(dims, deltas, offsets, mask)
end

function ImageGeom(
    dims::Dims{D},
    deltas::NTuple{D,RealU} = ntuple(i -> 1, D),
    offsets::NTuple{D,Real} = ntuple(i -> 0, D),
) where {D}
    ImageGeom(dims, deltas, offsets, Trues(dims))
end

"""
    ig = ImageGeom( ; dims=(nx,ny), deltas=(1,1), offsets=(0,0), mask=Trues )
Constructor using named keywords.
"""
function ImageGeom( ;
    dims::Dims{D} = (128,128),
    deltas::NTuple{D,RealU} = ntuple(i -> 1.0f0, length(dims)),
    offsets::NTuple{D,Real} = ntuple(i -> 0.0f0, length(dims)),
    mask::M = Trues(dims),
) where {D, M <: AbstractArray{Bool}}
    D == ndims(mask) || throw(DimensionMismatch("D=$D mask size $(size(mask))"))
    ImageGeom(dims, deltas, offsets, mask)
end

"""
    ig = ImageGeom(masktype::Mask ; kwargs...)
Constructor with specified mask, e.g., `ImageGeom(MaskCircle() ; dims=(6,8))`.
"""
function ImageGeom(masktype::Mask ; kwargs...)
    ig = ImageGeom( ; kwargs...)
    mask = makemask(ig, masktype)
    return ImageGeom( ; mask, kwargs...)
end

## Methods

# Simplify type display of Tuple to NTuple if possible
function _typeof(x::Tuple)
    D = length(x)
    return x isa NTuple{D, typeof(x[1])} ?
        "NTuple{$D,$(typeof(x[1]))}" : typeof(x)
end

"""
    show(io::IO, ::MIME"text/plain", ig::ImageGeom)
"""
function Base.show(io::IO, ::MIME"text/plain", ig::ImageGeom{D}) where D
    t = typeof(ig)
    t = replace(string(t), string(typeof(ig.deltas)) => string(_typeof(ig.deltas)))
    println(io, t)
    for f in (:dims, :deltas, :offsets)
        p = getproperty(ig, f)
        t = _typeof(p)
        println(io, " ", f, "::", t, " ", p)
    end
    f = :mask
    mask = getproperty(ig, f)
    print(io, " ", f, ":", " ", summary(mask))
    println(io, " {", count(mask), " of ", length(mask), "}")
end

Base.ndims(ig::ImageGeom{D}) where D = D
Base.size(ig::ImageGeom) = ig.dims
Base.size(ig::ImageGeom, d::Int) = ig.dims[d]

Base.zeros(T::DataType, ig::ImageGeom) = Zeros{T}(ig.dims...)
Base.zeros(ig::ImageGeom) = zeros(Float32, ig)
Base.ones(T::DataType, ig::ImageGeom) = Ones{T}(ig.dims...)
Base.ones(ig::ImageGeom) = ones(Float32, ig)

Base.trues(ig::ImageGeom) = Trues(ig.dims...)
Base.falses(ig::ImageGeom) = Falses(ig.dims...)

Base.zero(ig::ImageGeom{D,S}) where {D, S <: NTuple{D,Real}} = zero(Float32)
Base.zero(ig::ImageGeom{D,S}) where {D, S <: NTuple{D,Any}} = zero(Int32) # alert!
Base.zero(ig::ImageGeom{D,S}) where {D, S <: NTuple{D,T}} where {T <: Number} = zero(T)
Base.zero(ig::ImageGeom{D,S}) where {D, S <: NTuple{D,T}} where {T <: Real} = zero(T)

"""
    ig_down = downsample(ig, down::Tuple{Int})
Down sample an image geometry by the factor `down`.
Image size `ig.dims` must be ≥ `down`.
cf `image_geom_downsample`.
"""
function downsample(ig::ImageGeom{D,S,M}, down::NTuple{D,Int}) where {D,S,M}

    any(ig.dims .< down) && throw("dims=$(ig.dims) < down=$down")
    down_dim = ig.dims .÷ down # round down (if needed)
    deltas = ig.deltas .* down
    down_offsets = ig.offsets ./ down # adjust offsets to new "pixel" units

    # carefully down-sample the mask
    if ig.mask isa Trues # || all(ig.mask)
        down_mask = Trues(down_dim)
    else
        v = view(ig.mask, Base.OneTo.(down_dim .* down)...)
        down_mask = downsample(v, down)
    end

    return ImageGeom{D,S,M}(down_dim, deltas, down_offsets, down_mask)
end

downsample(ig::ImageGeom{D}, down::Int) where {D} =
    downsample(ig, ntuple(i -> down, D))

"""
    ig_over = oversample(ig::ImageGeom, over::Int or NTuple)
Over-sample an image geometry by the factor `over`.
"""
function oversample(ig::ImageGeom{D,S}, over::NTuple{D,Int}) where {D,S}
    if ig.mask isa Trues # || all(ig.mask)
        mask_over = Trues(ig.dims .* over)
    else
        mask_over = upsample(ig.mask, over)
    end
    return ImageGeom(ig.dims .* over, ig.deltas ./ over,
        ig.offsets .* over, mask_over,
    )
end

oversample(ig::ImageGeom{D}, over::Int) where {D} =
    oversample(ig, ntuple(i -> over, D))

"""
    ig_new = expand_nz(ig::ImageGeom{3}, nz_pad::Int)
Pad both ends along z.
"""
function expand_nz(ig::ImageGeom{3,S,M}, nz_pad::Int) where {S,M}
    out_nz = size(ig,3) + 2*nz_pad
    if ig.mask isa Trues
        out_mask = Trues(ig.dims[1:2]..., out_nz)
    else
        out_mask = cat(dims=3, repeat(ig.mask[:,:,1], 1, 1, nz_pad),
            ig.mask, repeat(ig.mask[:,:,end], 1, 1, nz_pad),
        )
    end
    return ImageGeom{3,S,M}((ig.dims[1], ig.dims[2], out_nz),
        ig.deltas, ig.offsets, out_mask,
    )
end

"""
    plot(ig, how ; kwargs...)
The `how` argument should be `MIRTjim.jim` to be useful.
"""
plot(ig::ImageGeom{2}, how::Function ; kwargs...) =
    how(axes(ig)..., ig.mask, "(nx,ny)=$(ig.dims)" ; kwargs...)
plot(ig::ImageGeom{3}, how::Function ; kwargs...) =
    how(axis(ig,1), axis(ig,2), mask_or(ig.mask),
        "(dx,dy,dz)=$(ig.deltas)" ; kwargs...)

# spatial axes
_center(n::Int, offset::Real) = (n - 1)/2 + offset
_axis(n::Int, Δ::RealU, offset::Real) = ((0:n-1) .- _center(n, offset)) * Δ
axis(ig::ImageGeom{D}, d::Int) where D =
    _axis(ig.dims[d], ig.deltas[d], ig.offsets[d])
Base.axes(ig::ImageGeom{D}) where D = ntuple(d -> axis(ig, d), Val(D))

grids(ig::ImageGeom) = ndgrid(axes(ig)...) # spatial grids

# spatial frequency axes
_axisf(n::Int, Δ::RealU) = (-n/2:n/2-1) / (n*Δ)
axisf(ig::ImageGeom{D}, d::Int) where {D} =
    _axisf(ig.dims[d], ig.deltas[d])
axesf(ig::ImageGeom{D}) where {D} = ntuple(d -> axisf(ig, d), Val(D))

gridf(ig::ImageGeom) = ndgrid(axesf(ig)...) # spectral grids


text:
"""
There's a lot of noise if we just go through all parameters of a `RenderObject`.
This function filters out the internal values
"""
function is_editable(k, v_v)
    v = value(v_v)
    !(
        k == :objectid ||
        k == :is_fully_opaque ||
        k == :instances ||
        k == Symbol("position.multiplicator") ||
        k == Symbol("position.dims") ||
        k == Symbol("resolution") ||
        k == :visible ||
        startswith(string(k), "boundingbox") ||
        (k == Symbol("color") && isa(v, AbstractArray)) ||
        k in fieldnames(PerspectiveCamera) ||
        k == :instances ||
        isa(v, Symbol) ||
        isa(v, Void) ||
        isa(v, NativeMesh) ||
        isa(v, Int) ||
        isa(v, Int32) ||
        isa(v, UInt32) ||
        isa(v, UInt) ||
        (isa(v, FixedVector) && eltype(v) <: Integer)
    )
end

makesignal2(s::Signal) = s
makesignal2(v) = Signal(v)
makesignal2(v::GPUArray) = v

function extract_edit_menu(robj::RenderObject, edit_screen, isvisible)
    extract_edit_menu(robj.uniforms, edit_screen, isvisible)
end
function extract_edit_menu(edit_dict::Dict, edit_screen, isvisible)
    lines = Point2f0[]
    screen_w = edit_screen.area.value.w-10mm
    labels = String[]
    glyph_scale = GLVisualize.glyph_scale!('X')
    pos = 1mm
    scale = 2mm ./ glyph_scale
    widget_text = scale .* 1.2f0
    glyph_height = round(Int, glyph_scale[2]*scale[2])
    atlas = GLVisualize.get_texture_atlas()
    font = GLVisualize.DEFAULT_FONT_FACE
    textpositions = Point2f0[]
    for (k,v) in edit_dict
        is_editable(k, v) || continue
        s = makesignal2(v)
        if applicable(widget, s, edit_screen)
            sig, vis = widget(s, edit_screen,
                visible=isvisible, text_scale=widget_text,
                area=(screen_w, value(icon_size)),
                knob_scale = 1.6mm
            )
            edit_dict[k] = sig
            bb = value(boundingbox(vis))
            height = widths(bb)[2]
            mini = minimum(bb)
            to_origin = -Vec3f0(mini[1], mini[2], 0)
            GLAbstraction.transform!(vis, translationmatrix(Vec3f0(2mm,pos,0)+to_origin))
            _view(vis, edit_screen, camera=:fixed_pixel)
            pos += round(Int, height) + 1mm

            label = replace(string(k), "_", " ")*":"
            push!(labels, label)
            append!(textpositions,
                GLVisualize.calc_position(label, Point2f0(1mm, pos), scale, font, atlas)
            )
            pos += glyph_height + 4mm
            push!(lines, Point2f0(0, pos-2mm), Point2f0(screen_w, pos-2mm))

        end
    end
    _view(visualize(
            join(labels), position=textpositions,
            color=RGBA{Float32}(0.8, 0.8, 0.8, 1.0),
            relative_scale=scale
        ), edit_screen, camera=:fixed_pixel
    )
    _view(visualize(
        lines, :linesegment, thickness=0.25mm, color=RGBA{Float32}(0.9, 0.9, 0.9, 1.0)
    ), edit_screen, camera=:fixed_pixel)

    pos
end


text:
export getGhostIndices

"""
function getGhotIndiced(n)

get indices of ghost points, boundary points and interior points for
cell-centered grid with n cells

Input:

   n::Array{Int}    -  number of cells

Output:

   idg   - indices of ghost points
   idgn  - indices of ghost point neighbors
   idint - indices of interior nodes

"""
function getGhostIndices(n::Array{Int})
	dim = length(n)

	if dim==2
		idg   = Array{Array{Int,1}}(undef, 4)
		idgn  = Array{Array{Int,1}}(undef, 4)

		idn    = reshape(collect(1:prod(n.+2)),tuple(n.+2...))::Array{Int,2}
		idg[1] = vec(idn[1,2:end-1])
		idg[2] = vec(idn[2:end-1,1])
		idg[3] = vec(idn[end,2:end-1])
		idg[4] = vec(idn[2:end-1,end])

		idgn[1] = vec(idn[2,2:end-1])
		idgn[2] = vec(idn[2:end-1,2])
		idgn[3] = vec(idn[end-1,2:end-1])
		idgn[4] = vec(idn[2:end-1,end-1])
		idint   = vec(idn[2:end-1,2:end-1])
	else
		idg   = Array{Array{Int,1}}(undef, 6)
		idgn  = Array{Array{Int,1}}(undef, 6)

		idn   = reshape(collect(1:prod(n.+2)),tuple(n.+2...))::Array{Int,3}
		idg[1] = vec(idn[1,2:end-1,2:end-1])
		idg[2] = vec(idn[2:end-1,1,2:end-1])
		idg[3] = vec(idn[2:end-1,2:end-1,1])
		idg[4] = vec(idn[end,2:end-1,2:end-1])
		idg[5] = vec(idn[2:end-1,end,2:end-1])
		idg[6] = vec(idn[2:end-1,2:end-1,end])

		idgn[1] = vec(idn[2,2:end-1,2:end-1])
		idgn[2] = vec(idn[2:end-1,2,2:end-1])
		idgn[3] = vec(idn[2:end-1,2:end-1,2])
		idgn[4] = vec(idn[end-1,2:end-1,2:end-1])
		idgn[5] = vec(idn[2:end-1,end-1,2:end-1])
		idgn[6] = vec(idn[2:end-1,2:end-1,end-1])

		idint = vec(idn[2:end-1,2:end-1,2:end-1])
	end
	return idg,idgn,idint
end

max_stars_repo_name: JuliaInv/ConvDiffMIPDECO

max_issues_repo_name: JuliaInv/ConvDiffMIPDECO

max_forks_repo_name: JuliaInv/ConvDiffMIPDECO


text:
#=
  Structs for the approxChol solver
=#

"""
  LLp elements are all in the same column.
  row tells us the row, and val is the entry.
  val is set to zero for some edges that we should remove.
  next gives the next in the column.  It points to itself to terminate.
  reverse is the index into lles of the other copy of this edge,
  since every edge is stored twice as we do not know the order of elimination in advance.
"""
mutable struct LLp{Tind,Tval}
    row::Tind
    val::Tval
    next::LLp{Tind,Tval}
    reverse::LLp{Tind,Tval}

    LLp{Tind,Tval}() where {Tind,Tval} = (x = new(zero(Tind), zero(Tval)); x.next = x; x.reverse = x)
    LLp{Tind,Tval}(row, val, next, rev) where {Tind,Tval} = new(row, val, next, rev)
    LLp{Tind,Tval}(row, val) where {Tind,Tval} = (x = new(row, val); x.next = x; x.reverse = x)
    LLp{Tind,Tval}(row, val, next) where {Tind,Tval} = (x = new(row, val, next); x.reverse = x)
end

"""
  LLmatp is the data structure used to maintain the matrix during elimination.
  It stores the elements in each column in a singly linked list (only next ptrs)
  Each element is an LLp (linked list pointer).
  The head of each column is pointed to by cols.

  We probably can get rid of degs - as it is only used to store initial degrees.
"""
struct LLmatp{Tind,Tval}
    n::Int64
    degs::Array{Tind,1}
    cols::Array{LLp{Tind,Tval},1}
    lles::Array{LLp{Tind,Tval},1}
end

# The followin are for an array-based structure
struct LLa{Tind,Tval}
    row::Tind
    val::Tval
    next::Tind
    reverse::Tind
end

struct LLmata{Tind,Tval}
    n::Int64
    degs::Vector{Tind}
    cols::Vector{Tind}
    lles::Vector{LLa{Tind,Tval}}
end

struct LLcol{Tind,Tval}
    row::Tind
    ptr::Tind
    val::Tval
end

# the following code might not work
bdd_type(B::BlockCSC{D,Tval} where {D, Tval}) = (D, Tval)
bdd_type(B::SparseMatrixCSC{Tval,Tind} where {Tval,Tind}) = (1, Tval)

#=============================================================

LDLinv

=============================================================#

struct IJop{Ti,Tv}
    i::Ti
    j::Ti
    wii::Tv
    wji::Tv
end

"""
Records that col i was used to elim entry in row j,
and that it was swapped into an entry in j, k.
"""
struct IJKop{Ti,Tv}
    i::Ti
    j::Ti
    k::Ti
    prob::Tv
end

"""
  LDLinv contains the information needed to solve the Laplacian systems.
  It does it by applying Linv, then Dinv, then Linv (transpose).
  But, it is specially constructed for this particular solver.
  It does not explicitly make the matrix triangular.
  Rather, col[i] is the name of the ith col to be eliminated
"""
mutable struct LDLinv{Tind,Tval}
    col::Array{Tind,1}
    colptr::Array{Tind,1}
    rowval::Array{Tind,1}
    fval::Array{Tval,1}
    d::Array{Tval,1}
end

abstract type AbsLTCSC{Tind, Tval} end

"""
A data type for holding a permutation of a lower-triangular matrix,
in a CSC-ish format.  The fields are colptr`, `rowval`, `nzval`.
The entries of the ith columns in the order are stored in rowval
and nzval in entries `colptr[i]`` through `colptr[i+1]-1`.  The identity
of the column is in `roval[colptr[i]]`, and its diagonal entry is in nzval.
No guarantees are made about the orders of other entries in that column.
"""
struct PermLTCSC{Tind,Tval} <: AbsLTCSC{Tind, Tval}
    colptr::Vector{Tind}
    rowval::Vector{Tind}
    nzval::Vector{Tval}
end

# similar type, but lower triangular although not nec sorted
struct LTCSC{Tind,Tval} <: AbsLTCSC{Tind, Tval}
    colptr::Vector{Tind}
    rowval::Vector{Tind}
    nzval::Vector{Tval}
end

#=============================================================

ApproxCholPQ
the data strcture we use to keep track of degrees

=============================================================#

struct acPQ_elem{Tind}
    prev::Tind
    next::Tind
    key::Tind
end

"""
  An approximate priority queue.
  Items are bundled together into doubly-linked lists with all approximately the same key.
  minlist is the min list we know to be non-empty.
  It should always be a lower bound.
  keyMap maps keys to lists
"""
mutable struct acPQ{Tind}
    elems::Array{acPQ_elem{Tind},1} # indexed by node name
    lists::Array{Tind,1}
    minlist::Int
    nitems::Int
    n::Int
end

mutable struct acPQ_rev{Tind}
    elems::Array{acPQ_elem{Tind},1} # indexed by node name
    heads::Array{Tind,1}  # pop from heads: prev points to head
    tails::Array{Tind,1}  # add to tails, next points to tail
    minlist::Int
    nitems::Int
    n::Int
end

mutable struct randishPQ{Tind}
    elems::Array{acPQ_elem{Tind},1} # indexed by node name
    heads::Array{Tind,1}  # pop from heads: prev points to head
    tails::Array{Tind,1}  # add to tails, next points to tail
    minlist::Int
    nitems::Int
    n::Int
end

struct NodeDeg{Tind}
    node::Tind
    deg::Tind
end

mutable struct RandPQ{Tind}
    elems::Vector{NodeDeg{Tind}}
    starts::Vector{Tind}
    nodeptr::Vector{Tind}
    cur::Tind
end

struct Compression
    col::Int
    into::Int
    from::Int
end

struct IJKop2{Ti,Tv}
    i::Ti
    j::Ti
    k::Ti
    prob::Tv
    ptr_ji::Ti
    ptr_ki::Ti
    ptr_kj::Ti
    #ptr_jk::Ti equals ptr_ji
end

struct PCG_params
    verbose::Bool
    tol::Float64
    maxits::Float64
    maxtime::Float64
    pcgIts::Vector{Int}
    stag_test::Integer
end

PCG_params() = PCG_params(false, 1e-6, Inf, Inf, Int[], 5)


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "libfreenect"
version = v"0.6.1"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/OpenKinect/libfreenect.git", "edad69caabe662e7996db9425c8b5c568efa275d"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/libfreenect/
atomic_patch -p1 ../patches/libusb_header.patch
mkdir build && cd build/
FLAGS=()
if [[ "${target}" == *-mingw* ]]; then
    # Give a hint about where to find pthread library
    FLAGS+=(-DTHREADS_PTHREADS_WIN32_LIBRARY="/opt/${target}/${target}/sys-root/bin/libwinpthread-1.dll")
fi
cmake -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_BUILD_TYPE=Release \
    -DBUILD_CPP=OFF \
    -DBUILD_EXAMPLES=OFF \
    -DBUILD_FAKENECT=OFF \
    -DBUILD_REDIST_PACKAGE=OFF \
    "${FLAGS[@]}" \
    ..
make -j${nproc}
make install
install_license ../APACHE20 ../GPL2
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct("libfreenect_sync", :libfreenect_sync),
    LibraryProduct("libfreenect", :lifreenect)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="libusb_jll", uuid="a877fdc9-fe69-5ed6-b93d-11ecd0dc2d49"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: dmillard/Yggdrasil

max_issues_repo_name: dmillard/Yggdrasil

max_forks_repo_name: dmillard/Yggdrasil


text:
using MatrixNetworks

"""

Ear Decomposition
-----------------

An ear decomposition of an undirected graph G is a partition of its set of edges into a sequence of ears, such that the one or two endpoints of each ear
belong to earlier ears in the sequence and such that the internal vertices of each ear do not belong to any earlier ear.

An open ear decomposition or a proper ear decomposition is an ear decomposition in which the two endpoints of each ear after the first are distinct from each other.

The code is based on the algorithm given by Schmidt et.al in  "A simple test on 2-vertex- and 2-edge-connectivity" (Information Processing Letters Journal'
2013).

This approach is an important subroutine in many graph problems like Planarity Testing, Identifying the Factor critical graphs, recognizing series parallel
graphs, identifying Tri connected components etc.

Input: An undirected graph with i--->j and j--->i present (NO self loops!)

Observed:

      1) Slower on random graphs generated using GT-graph Suite ( Erdos Reyni model - with 20K vertices and a sparse quotient of 10 ) by a factor 5-8%

TODO: 1) Add a function to test for bi-connectivity and 2-edge connectivity.
      2) Compare the Performance with C ( SNAP DATA set )

"""

I = [1,2,2,3,3,4,1,4,3,5,5,6,6,7,7,8,6,8,5,9,9,10,5,10,1,3,2,5]
J = [2,1,3,2,4,3,4,1,5,3,6,5,7,6,8,7,8,6,9,5,10,9,10,5,3,1,5,2]
V = [1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1]

A = sparse(I,J,V,10,10)

function ear(A::MatrixNetwork,u::Int64,full::Int64,target::Int64)

    (rp,ci) = (A.rp,A.ci)
    n=length(rp)-1
    d=-1*ones(Int64,n)
    dt=-1*ones(Int64,n)
    rev_dt= -1*ones(Int64,n)
    ft=-1*ones(Int64,n)
    pred=zeros(Int64,n)
    rs=zeros(Int64,2*n)
    rss=0 # recursion stack holds two nums (v,ri)
    be = Any[] #back edges array
    ears = zeros(Int64,2 * length(A.ci))
    ear_number = zeros(Int64,length(A.ci))
    ear_index = 0
    ear_visited = zeros(Int64, A.n)

    for i=1:n
        push!(be, [])
    end

    #start dfs at u
    t=1
    targethit=0
    for i=1:n
        if i==1
            v=u
        else
            v=mod(u+i-1,n)+1
            if d[v]>0
                continue
            end
        end
        d[v]=1
        dt[v]=t
        rev_dt[t]=v;
	t=t+1
        ri=rp[v]
        rss=rss+1
        rs[2*rss-1]=v
        rs[2*rss]=ri # add v to the stack

        while rss>0
            v=rs[2*rss-1]
            ri=rs[2*rss]
            rss=rss-1 # pop v from the stack
            if v==target || targethit == 1
                ri=rp[v+1]
                targethit=1 # end the algorithm if v is the target
            end
            while ri<rp[v+1]
                w=ci[ri]
                ri=ri+1
                if d[w]<0
                    d[w]=d[v]+1
                    pred[w]=v
                    rss=rss+1
                    rs[2*rss-1]=v
                    rs[2*rss]=ri # add v to the stack
                    v=w
                    ri=rp[w]
                    dt[v]=t
		    rev_dt[t]=v
                    t=t+1
                #    continue #discover a new vertex!
                elseif (pred[v]!=w)
		    if(dt[w]< dt[v])
		       push!(be[w], v)
		    end
		end
            end
        end
        if full == 0
            break
        end
    end
    cn=1
    for i=1:n
	startVertex = rev_dt[i]
	ear_visited[startVertex] = 1
	while !isempty(be[startVertex])  # A disjoint Path (or edge) is identified.
	    nextVertex = pop!(be[startVertex])
            ear_index += 1
            ears[2*ear_index-1] = startVertex
            ears[2*ear_index] = nextVertex
            ear_number[ear_index] = cn
            while (ear_visited[nextVertex] != 1) #Identify the ear (Chain) and number it
                ear_visited[nextVertex] = 1
                startVertex_n = nextVertex
                nextVertex = pred[startVertex_n]
                ear_index += 1
                ears[2*ear_index-1] = startVertex_n # Store the ear (Chain).
                ears[2*ear_index] = nextVertex
                ear_number[ear_index] = cn  # Provide numbering to that Edge (#ear)
            end
            cn += 1
	end
    end
    print(ears)
    return (ears,ear_number)
end

function dfs_test()
#print(A);
B = MatrixNetwork(A);
#print(B)
ear(B, 1,1,0)
end

dfs_test();


text:
# Dimensionality operations for 4D Tensors
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: January 2020

export squeeze, unsqueeze, wavelet_squeeze, wavelet_unsqueeze, Haar_squeeze, invHaar_unsqueeze, tensor_split, tensor_cat
export cat_states, split_states

####################################################################################################
# Squeeze and unsqueeze
function patch_inds(N::NTuple{n, Int}, d::Integer) where n
    indsX = [ix*N[1]+1:(ix+1)*N[1] for ix=[(i+1)%2 for i=1:2^(d-2)]]
    indsY = [iy*N[2]+1:(iy+1)*N[2] for iy=[div(i-1,2)%2 for i=1:2^(d-2)]]
    d == 4 && (return zip(indsX, indsY))
    indsZ = [iz*N[3]+1:(iz+1)*N[3] for iz=[div(i-1,4)%2 for i=1:2^(d-2)]]
    return zip(indsX, indsY, indsZ)
end

function checkboard_inds(N::NTuple{n, Int}, d::Integer) where n
    indsX = [ix+1:2:N[1] for ix=[(i+1)%2 for i=1:2^(d-2)]]
    indsY = [iy+1:2:N[2] for iy=[div(i-1,2)%2 for i=1:2^(d-2)]]
    d == 4 && (return zip(indsX, indsY))
    indsZ = [iz:2:N[3] for iz=[div(i-1,4)%2 for i=1:2^(d-2)]]
    return zip(indsX, indsY, indsZ)
end

"""
    Y = squeeze(X; pattern="column")

 Reshape input image such that each spatial dimension is reduced by a factor
 of 2, while the number of channels is increased by a factor of 4.

 *Input*:

 - `X`: 4D input tensor of dimensions `nx` x `ny` x `n_channel` x `batchsize`

 - `pattern`: Squeezing pattern

        1 2 3 4        1 1 3 3        1 3 1 3
        1 2 3 4        1 1 3 3        2 4 2 4
        1 2 3 4        2 2 4 4        1 3 1 3
        1 2 3 4        2 2 4 4        2 4 2 4

        column          patch       checkerboard

 *Output*:

 - `Y`: Reshaped tensor of dimensions `nx/2` x `ny/2` x `n_channel*4` x `batchsize`

 See also: [`unsqueeze`](@ref), [`wavelet_squeeze`](@ref), [`wavelet_unsqueeze`](@ref)
"""
function squeeze(X::AbstractArray{T, N}; pattern="column") where {T, N}

    # Dimensions
    nc_in, batchsize = size(X)[N-1:N]
    if any([mod(nn, 2) == 1 for nn=size(X)[1:N-2]])
        throw("Input dimensions must be multiple of 2")
    end
    N_out = Tuple(nn÷2 for nn=size(X)[1:N-2])
    nc_out = size(X, N-1) * 2^(N-2)
    cinds = Tuple((:) for i=1:N-2)

    if pattern == "column"
        Y = reshape(X, N_out..., nc_out, batchsize)
    elseif pattern == "patch"
        Y = cuzeros(X, N_out..., nc_out, batchsize)
        iX = patch_inds(N_out, N)
        for (i, ix)=enumerate(iX)
            Y[cinds..., (i-1)*nc_in+1:i*nc_in, :] = X[ix..., :, :]
        end
    elseif pattern == "checkerboard"
        Y = cuzeros(X, N_out..., nc_out, batchsize)
        iX = checkboard_inds(size(X), N)
        for (i, ix)=enumerate(iX)
            Y[cinds..., (i-1)*nc_in+1:i*nc_in, :] = X[ix..., :, :]
        end
    else
        throw("Specified pattern not defined.")
    end
    return Y
end

"""
    X = unsqueeze(Y; pattern="column")

 Undo squeezing operation by reshaping input image such that each spatial dimension is
 increased by a factor of 2, while the number of channels is decreased by a factor of 4.

 *Input*:

 - `Y`: 4D input tensor of dimensions `nx` x `ny` x `n_channel` x `batchsize`

 - `pattern`: Squeezing pattern

            1 2 3 4        1 1 3 3        1 3 1 3
            1 2 3 4        1 1 3 3        2 4 2 4
            1 2 3 4        2 2 4 4        1 3 1 3
            1 2 3 4        2 2 4 4        2 4 2 4

            column          patch       checkerboard

 *Output*:

 - `X`: Reshaped tensor of dimensions `nx*2` x `ny*2` x `n_channel/4` x `batchsize`

 See also: [`squeeze`](@ref), [`wavelet_squeeze`](@ref), [`wavelet_unsqueeze`](@ref)
"""
function unsqueeze(Y::AbstractArray{T,N}; pattern="column") where {T, N}

    # Dimensions
    batchsize = size(Y, N)
    if any([mod(nn, 2) == 1 for nn=size(Y)[1:N-2]])
        throw("Input dimensions must be multiple of 2")
    end
    N_out = Tuple(nn*2 for nn=size(Y)[1:N-2])
    nc_out = size(Y, N-1) ÷ 2^(N-2)
    cinds = Tuple((:) for i=1:N-2)

    if pattern == "column"
        X = reshape(Y, N_out..., nc_out, batchsize)
    elseif pattern == "patch"
        X = cuzeros(Y, N_out..., nc_out, batchsize)
        iX = patch_inds(size(Y), N)
        for (i, ix)=enumerate(iX)
            X[ix..., :, :] = Y[cinds..., (i-1)*nc_out+1:i*nc_out, :]
        end
    elseif pattern == "checkerboard"
        X = cuzeros(Y, N_out..., nc_out, batchsize)
        iX = checkboard_inds(N_out, N)
        for (i, ix)=enumerate(iX)
            X[ix..., :, :] = Y[cinds..., (i-1)*nc_out+1:i*nc_out, :]
        end
    else
        throw("Specified pattern not defined.")
    end
    return X
end

function unsqueeze(X::AbstractArray{T,N}, Y::AbstractArray{T,4}; pattern="column") where {T,N}
    return unsqueeze(X; pattern=pattern), unsqueeze(Y; pattern=pattern)
end

####################################################################################################
# Squeeze and unsqueeze using the wavelet transform

"""
    Y = wavelet_squeeze(X; type=WT.db1)

 Perform a 1-level channelwise 2D wavelet transform of X and squeeze output of each
 transform into 4 channels (per 1 input channel).

 *Input*:

 - `X`: 4D input tensor of dimensions `nx` x `ny` x `n_channel` x `batchsize`

 - `type`: Wavelet filter type. Possible values are `WT.haar` for Haar wavelets,
    `WT.coif2`, `WT.coif4`, etc. for Coiflet wavelets, or `WT.db1`, `WT.db2`, etc.
    for Daubechies wavetlets. See *https://github.com/JuliaDSP/Wavelets.jl* for a
    full list.

 *Output*:

 - `Y`: Reshaped tensor of dimensions `nx/2` x `ny/2` x `n_channel*4` x `batchsize`

 See also: [`wavelet_unsqueeze`](@ref), [`squeeze`](@ref), [`unsqueeze`](@ref)
"""
function wavelet_squeeze(X::AbstractArray{T, N}; type=WT.db1) where {T, N}
    batchsize = size(X, N)
    N_in = size(X)[1:N-2]
    N_out = Tuple(nn÷2 for nn=size(X)[1:N-2])
    nd = 2^(N-2)
    nc_out = size(X, N-1) * nd
    cinds = Tuple((:) for i=1:N-2)

    Y = cuzeros(X, N_out..., nc_out, batchsize)
    for i=1:batchsize
        for j=1:size(X, N-1)
            Ycurr = dwt(X[cinds..., j, i], wavelet(type), 1)
            Y[cinds..., (j-1)*nd + 1: j*nd, i] = squeeze(reshape(Ycurr, N_in..., 1, 1); pattern="patch")
        end
    end
    return Y
end

"""
    X = wavelet_unsqueeze(Y; type=WT.db1)

 Perform a 1-level inverse 2D wavelet transform of Y and unsqueeze output.
 This reduces the number of channels by 4 and increases each spatial
 dimension by a factor of 2. Inverse operation of `wavelet_squeeze`.

 *Input*:

 - `Y`: 4D input tensor of dimensions `nx` x `ny` x `n_channel` x `batchsize`

 - `type`: Wavelet filter type. Possible values are `haar` for Haar wavelets,
  `coif2`, `coif4`, etc. for Coiflet wavelets, or `db1`, `db2`, etc. for Daubechies
  wavetlets. See *https://github.com/JuliaDSP/Wavelets.jl* for a full list.

 *Output*:

 - `X`: Reshaped tensor of dimenions `nx*2` x `ny*2` x `n_channel/4` x `batchsize`

 See also: [`wavelet_squeeze`](@ref), [`squeeze`](@ref), [`unsqueeze`](@ref)
"""
function wavelet_unsqueeze(Y::AbstractArray{T,N}; type=WT.db1) where {T, N}
    N_out = Tuple(nn*2 for nn=size(Y)[1:N-2])
    nd = 2^(N-2)
    nc_out = size(Y, N-1) ÷ nd
    cinds = Tuple((:) for i=1:N-2)

    X = cuzeros(Y, N_out..., nc_out, size(Y, N))
    for i=1:size(Y, N)
        for j=1:nc_out
            Ycurr = unsqueeze(Y[cinds..., (j-1)*nd + 1: j*nd, i:i]; pattern="patch")[cinds..., 1, 1]
            Xcurr = idwt(Ycurr, wavelet(type), 1)
            X[cinds..., j, i] = reshape(Xcurr, N_out..., 1, 1)
        end
    end
    return X
end

########## Haaar wavelet, GPU supported #####################

function HaarLift(x::AbstractArray{T,N}, dim) where {T, N}
    dim > N - 2 && return (x,)
    #Haar lifting
    inds = [i==dim ? (1:2:size(x, dim)) : (:) for i=1:N];
    # Splitting
    H = x[inds...]
    inds[dim] = 2:2:size(x, dim)
    L = x[inds...]

    # predict
    H .-= L

    #update
    L .+= H ./ T(2.0)

    #normalize
    H ./= sqrt(T(2.0))
    L .*= sqrt(T(2.0))

    return L, H
end

function invHaarLift(x::AbstractArray{T,N}, dim) where {T, N}
    dim > N - 2 && return x
    return invHaarLift(tensor_split(x)..., dim)
end

function invHaarLift(L::AbstractArray{T,N}, H::AbstractArray{T,N},dim) where {T, N}
    #inverse Haar lifting

    #inv normalize
    H .*= sqrt(T(2.0))
    L ./= sqrt(T(2.0))

    #inv update & predict
    L .-= H ./ T(2.0)
    H .+= L

    #allocate output:
    x = cat(cuzeros(L,size(L)...), cuzeros(H,size(H)...), dims=dim)
    inds = [i==dim ? (1:2:size(x, dim)) : (:) for i=1:N];
    x[inds...] .= H
    inds[dim] = 2:2:size(x, dim)
    x[inds...] .= L

    return x
end

"""
    Y = Haar_squeeze(X)

 Perform a 1-level channelwise 2D/3D (lifting) Haar transform of X and squeeze output of each
 transform into 8 channels (per 1 input channel).

 *Input*:

 - `X`: 4D/5D input tensor of dimensions `nx` x `ny` (x `nz`) x `n_channel` x `batchsize`

 *Output*:

 - `Y`: Reshaped tensor of dimensions `nx/2` x `ny/2` (x `nz/2`) x `n_channel*8` x `batchsize`

 See also: [`wavelet_unsqueeze`](@ref), [`Haar_unsqueeze`](@ref), [`HaarLift`](@ref), [`squeeze`](@ref), [`unsqueeze`](@ref)
"""
function Haar_squeeze(x::AbstractArray{T, N}) where {T, N}

        L, H = HaarLift(x, 2)

        a, h = HaarLift(L ,1)
        v, d = HaarLift(H, 1)

        a = HaarLift(a, 3)
        v = HaarLift(v, 3)
        h = HaarLift(h, 3)
        d = HaarLift(d, 3)

        return cat(a..., v..., h..., d...,dims=N-1)
end

"""
    X = invHaar_unsqueeze(Y)

 Perform a 1-level inverse 2D/3D Haar transform of Y and unsqueeze output.
 This reduces the number of channels by 8 and increases each spatial
 dimension by a factor of 2. Inverse operation of `Haar_squeeze`.

 *Input*:

 - `Y`: 4D/5D input tensor of dimensions `nx` x `ny` (x `nz`) x `n_channel` x `batchsize`

 *Output*:

 - `X`: Reshaped tensor of dimenions `nx*2` x `ny*2` (x `nz*2`) x `n_channel/8` x `batchsize`

 See also: [`wavelet_unsqueeze`](@ref), [`Haar_unsqueeze`](@ref), [`HaarLift`](@ref), [`squeeze`](@ref), [`unsqueeze`](@ref)
"""
function invHaar_unsqueeze(x::AbstractArray{T, N}) where {T, N}

        s = size(x, N-1)
        a, h = tensor_split(x)
        a, v, h, d = tensor_split(a)..., tensor_split(h)...

        a = invHaarLift(a, 3)
        v = invHaarLift(v, 3)
        h = invHaarLift(h, 3)
        d = invHaarLift(d, 3)

        L = invHaarLift(a, h, 1)
        H = invHaarLift(v, d, 1)

        x = invHaarLift(L, H, 2)

        return x
end

####################################################################################################
# Split and concatenate

"""
    Y, Z = tensor_split(X)

 Split ND input tensor in half along the channel dimension. Inverse operation
 of `tensor_cat`.

 *Input*:

 - `X`: ND input tensor of dimensions `nx` [x `ny` [x `nz`]] x `n_channel` x `batchsize`

 *Output*:

 - `Y`, `Z`: ND output tensors, each of dimensions `nx` [x `ny` [x `nz`]] x `n_channel/2` x `batchsize`

 See also: [`tensor_cat`](@ref)
"""
function tensor_split(X::AbstractArray{T, N}; split_index=nothing) where {T, N}
    d = max(1, N-1)
    if isnothing(split_index)
        k = Int(round(size(X, d)/2))
    else
        k = split_index
    end

    indsl = [i==d ? (1:k) : (:) for i=1:N]
    indsr = [i==d ? (k+1:size(X, d)) : (:) for i=1:N]

    return X[indsl...], X[indsr...]
end

"""
    X = tensor_cat(Y, Z)

 Concatenate ND input tensors along the channel dimension. Inverse operation
 of `tensor_split`.

 *Input*:

 - `Y`, `Z`: ND input tensors, each of dimensions `nx` [x `ny` [x `nz`]] x `n_channel` x `batchsize`

 *Output*:

 - `X`: ND output tensor of dimensions `nx` [x `ny` [x `nz`]] x `n_channel*2` x `batchsize`

 See also: [`tensor_split`](@ref)
"""
function tensor_cat(X::AbstractArray{T, N}, Y::AbstractArray{T, N}) where {T, N}
    d = max(1, N-1)
    if size(X, d) == 0
        return Y
    elseif size(Y, d) == 0
        return X
    else
        return cat(X, Y; dims=d)
    end
end

tensor_cat(X::Tuple{AbstractArray{T,N}, AbstractArray{T,N}}) where {T, N} = tensor_cat(X[1], X[2])

@inline xy_dims(dims::Tuple, ::Val{false}) = dims
@inline xy_dims(dims::Tuple, ::Val{true}) = Int.(dims .* (.5, .5, 4, 1))

# Concatenate states Zi and final output
function cat_states(XY_save::AbstractMatrix{<:AbstractArray}, X::AbstractArray{T, 4}, Y::AbstractArray{T, 4}) where T
    return cat_states(XY_save[:, 1], X), cat_states(XY_save[:, 2], Y)
end
# Concatenate states Zi and final output
function cat_states(Z_save::Vector{<:AbstractArray}, X::AbstractArray{T, N}) where {T, N}
    return cat([[vec(Z_save[j]) for j=1:length(Z_save)]..., vec(X)]..., dims=1)
end

# Split and reshape 1D vector Y in latent space back to states Zi
# where Zi is the split tensor at each multiscale level.
function split_states(Y::AbstractVector{T}, Z_dims) where {T, N}
    L = length(Z_dims) + 1
    inds = cumsum([1, [prod(Z_dims[j]) for j=1:L-1]...])
    Z_save = [reshape(Y[inds[j]:inds[j+1]-1], xy_dims(Z_dims[j], Val(j==L))) for j=1:L-1]
    X = reshape(Y[inds[L]:end], xy_dims(Z_dims[end], Val(true)))
    return Z_save, X
end

# Split and reshape 1D vector X_full and Y_vull in latent space back to states Zi
# where Zi is the split tensor at each multiscale level.
function split_states(XY_dims::AbstractArray{Tuple, 1}, X_full::AbstractArray{T, 1}, Y_full::AbstractArray{T, 1}) where T
    c1, X = split_states(X_full, XY_dims)
    c2, Y = split_states(Y_full, XY_dims)
    return hcat(c1, c2), X, Y
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    HierarchicalPartition(first, second)

A partitioning method in which a `first` partition is applied
and then a `second` partition is applied to each subset of the
`first`.
"""
struct HierarchicalPartition <: PartitionMethod
  first::PartitionMethod
  second::PartitionMethod
end

function partition(object, method::HierarchicalPartition)
  result = Vector{Int}[]

  # use first partition method
  p = partition(object, method.first)

  # use second method to partition the first
  s = indices(p)
  for (i, d) in Iterators.enumerate(p)
    q = partition(d, method.second)

    for js in indices(q)
      push!(result, s[i][js])
    end
  end

  Partition(object, result)
end


text:
file = joinpath(@__DIR__, "documents", "jupyter_test.jmd")
using IJulia, Conda

Conda.add("nbconvert") # should be the same as IJulia.JUPYTER, i.e. the miniconda Python

Weave.notebook(file, jupyter_path = IJulia.JUPYTER)
@test "jupyter_test.ipynb" ∈ readdir(@__DIR__) # test if the result was weaved


text:
struct Extrapolator

    Proj :: SparseMatrixCSC
    R :: SparseMatrixCSC
    flip_factors :: Matrix{Int64}

end

function Extrapolator(M::Mirror, grd::Grid)

    chnl = RemoteChannel(()->Channel{Bool}(), 1)
    p = Progress(grd.Nk, desc="Finding mirror images... ")
    @async while take!(chnl)
        next!(p)
    end

    map_out = pmap(1:grd.Nk) do k

        _, proj_js, _ = findnz(grd.Proj)
        k_cart = proj_js[k]
        h, sindex = distance_to_mirror(grd.points[:,k]..., M)
        flip_row = ones(Int64, size(M.verts)[2])

        if h < 0

            flip_row[sindex] = -1

            j = Int(ceil(k_cart / grd._Nx))
            i = k_cart - (j - 1) * grd._Nx
            x, y, di, dj = mirror_image(grd.points[:,k]..., grd.dr, grd.dr, M)
            i += di
            j += dj
            # @assert distance_to_mirror(x, y, M)[1] > 0
            @assert !((di == 0) && (dj == 0))
            k_cart = i + (j-1)*grd._Nx
            k = 0
        end

        # @info k
        put!(chnl, true)
        Int64[k, k_cart, flip_row...]
    end

    C = hcat(map_out...)

    proj_j = C[1,:]
    proj_j = proj_j[proj_j .!= 0]
    r_j = C[2,:]
    flip_factors = C[3:end,:]

    Mirror = sparse([1:grd.Nk...], r_j, ones(grd.Nk), grd.Nk, grd._Nx*grd._Ny) * transpose(grd.Proj)
    Proj = sparse([1:length(proj_j)...], proj_j, ones(length(proj_j)), length(proj_j), grd.Nk)

    return Extrapolator(Proj, Mirror, flip_factors)
end

function make_dirichlet(a::Extrapolator, inds)

    factors = a.flip_factors[inds,:]
    N = size(factors)[2]
    v = ones(N)
    for i=1:size(factors)[1]
        v = v .* factors[i,:]
    end

    R = Diagonal(v)
    return Extrapolator(a.Proj, R, a.flip_factors)
end

function (a::Extrapolator)(x::Vector{Float64}, xb::Vector{Float64})

    if length(x) < length(xb)
        v = transpose(a.Proj)*x
    else
        v = x[:]
    end
    return a.R*v + (I - a.R)*xb
end

function (a::Extrapolator)(A::SparseMatrixCSC, b::Vector{Float64}, xb::Vector{Float64})

    P = a.Proj
    Pt = transpose(P)

    Anew = P * A * a.R * Pt
    bnew = P*b + P*A*(a.R - I)*xb

    return Anew, bnew
end


text:
if VERSION > v"0.7-"
    push!(LOAD_PATH, Base.NamedEnv("SegmentTrees.jl"))
    using InteractiveUtils
    using Random
    using Profile
end
versioninfo()

using SegmentTrees
using CodecZlib

function readgff3(input::IO, chrom::String)
    segments = Segment{Int}[]
    for line in eachline(input)
        if startswith(line, "#")
            continue
        end
        values = split(line, '\t')
        if values[1] == chrom
            chromstart = parse(Int, values[4])
            chromend = parse(Int, values[5])
            push!(segments, Segment(chromstart, chromend))
        end
    end
    sort!(segments)
    return SegmentTree{Int,Segment{Int}}(segments)
end

function readgff3(filepath::String, chrom::String)
    stream = open(filepath)
    if endswith(filepath, ".gz")
        stream = GzipDecompressorStream(stream)
    end
    try
        return readgff3(stream, chrom)
    finally
        close(stream)
    end
end

isdir("data") || mkdir("data")
cd("data")
filename="Danio_rerio.GRCz11.92.gff3.gz"
isfile(filename)||
    download("ftp://ftp.ensembl.org/pub/release-92/gff3/danio_rerio/Danio_rerio.GRCz11.92.gff3.gz", filename)
cd("..")

tree = readgff3("data/Danio_rerio.GRCz11.92.gff3.gz", "1")

println("random query")
function random_query(tree, queries)
    n = 0
    for q in queries
        for _ in intersect(tree, q)
            n += 1
        end
    end
    return n
end
srand(1234)
segments = shuffle!(collect(tree))
for _ in 1:5
    @time random_query(tree, segments)
end
if "-p" in ARGS
    @profile random_query(tree, segments)
    Profile.print()
end


text:
module parsim
function simulate()
end
end


text:
export ±

"""
```julia
±(a,b)
```
- returns (a-b,a+b)
```julia
# Examples

1 ± 0.5  # returns (0.5,1.5)
```
"""
±(a,b) = [a-b,a+b]


text:
""
function constraint_voltage_magnitude_link(pm::_PM.AbstractACRModel, n_1::Int, n_2::Int, i::Int)
    vr_1 = var(pm, n_1, :vr, i)
    vi_1 = var(pm, n_1, :vi, i)

    vr_2 = var(pm, n_2, :vr, i)
    vi_2 = var(pm, n_2, :vi, i)

    JuMP.@constraint(pm.model, vr_1^2 + vi_1^2 == vr_2^2 + vi_2^2)
end

""
function constraint_power_balance_shunt_dispatch(pm::_PM.AbstractACRModel, n::Int, i::Int, bus_arcs, bus_arcs_dc, bus_arcs_sw, bus_gens, bus_storage, bus_shunts_var, bus_pd, bus_qd, bus_gs_const, bus_bs_const)
    vi = var(pm, n, :vi, i)
    vr = var(pm, n, :vr, i)
    p    = get(var(pm, n),    :p, Dict()); _PM._check_var_keys(p, bus_arcs, "active power", "branch")
    q    = get(var(pm, n),    :q, Dict()); _PM._check_var_keys(q, bus_arcs, "reactive power", "branch")
    pg   = get(var(pm, n),   :pg, Dict()); _PM._check_var_keys(pg, bus_gens, "active power", "generator")
    qg   = get(var(pm, n),   :qg, Dict()); _PM._check_var_keys(qg, bus_gens, "reactive power", "generator")
    ps   = get(var(pm, n),   :ps, Dict()); _PM._check_var_keys(ps, bus_storage, "active power", "storage")
    qs   = get(var(pm, n),   :qs, Dict()); _PM._check_var_keys(qs, bus_storage, "reactive power", "storage")
    psw  = get(var(pm, n),  :psw, Dict()); _PM._check_var_keys(psw, bus_arcs_sw, "active power", "switch")
    qsw  = get(var(pm, n),  :qsw, Dict()); _PM._check_var_keys(qsw, bus_arcs_sw, "reactive power", "switch")
    p_dc = get(var(pm, n), :p_dc, Dict()); _PM._check_var_keys(p_dc, bus_arcs_dc, "active power", "dcline")
    q_dc = get(var(pm, n), :q_dc, Dict()); _PM._check_var_keys(q_dc, bus_arcs_dc, "reactive power", "dcline")

    bs = get(var(pm, n), :bs, Dict()); _PM._check_var_keys(bs, bus_shunts_var, "reactive power", "shunt")

    # possibly can save 2x in function eval, but no the dominant runtime at this moment
    #vm_sqr = @variable(pm.model, start=1.0, base_name="$(0)_vm_sqr_$(i)")

    #JuMP.@constraint(pm.model, vm_sqr == vi^2 + vr^2)
    #cstr_p = JuMP.@constraint(pm.model, 0 == - sum(p[a] for a in bus_arcs) + sum(pg[g] for g in bus_gens) - sum(pd for pd in values(bus_pd)) - sum(gs for gs in values(bus_gs_const))*vm_sqr)
    #cstr_q = JuMP.@constraint(pm.model, 0 == - sum(q[a] for a in bus_arcs) + sum(qg[g] for g in bus_gens) - sum(qd for qd in values(bus_qd)) + sum(bs for bs in values(bus_bs_const))*vm_sqr + sum(bs[s]*vm_sqr for s in bus_shunts_var))

    cstr_p = JuMP.@constraint(pm.model, 0 == - sum(p[a] for a in bus_arcs) + sum(pg[g] for g in bus_gens) - sum(pd for pd in values(bus_pd)) - sum(gs for gs in values(bus_gs_const))*(vi^2 + vr^2))
    cstr_q = JuMP.@NLconstraint(pm.model, 0 == - sum(q[a] for a in bus_arcs) + sum(qg[g] for g in bus_gens) - sum(qd for qd in values(bus_qd)) + sum(bs for bs in values(bus_bs_const))*(vi^2 + vr^2) + sum(bs[s]*(vi^2 + vr^2) for s in bus_shunts_var))

    if _IM.report_duals(pm)
        sol(pm, n, :bus, i)[:lam_kcl_r] = cstr_p
        sol(pm, n, :bus, i)[:lam_kcl_i] = cstr_q
    end
end

function constraint_ohms_yt_from_goc(pm::_PM.AbstractACRModel, n::Int, f_bus, t_bus, f_idx, t_idx, g, b, g_fr, b_fr, tr, ti, tm)
    p_fr = var(pm, n, :p, f_idx)
    q_fr = var(pm, n, :q, f_idx)
    vr_fr = var(pm, n, :vr, f_bus)
    vr_to = var(pm, n, :vr, t_bus)
    vi_fr = var(pm, n, :vi, f_bus)
    vi_to = var(pm, n, :vi, t_bus)

    JuMP.@constraint(pm.model, p_fr ==  (g/tm^2+g_fr)*(vr_fr^2 + vi_fr^2) + (-g*tr+b*ti)/tm^2*(vr_fr*vr_to + vi_fr*vi_to) + (-b*tr-g*ti)/tm^2*(vi_fr*vr_to - vr_fr*vi_to) )
    JuMP.@constraint(pm.model, q_fr == -(b/tm^2+b_fr)*(vr_fr^2 + vi_fr^2) - (-b*tr-g*ti)/tm^2*(vr_fr*vr_to + vi_fr*vi_to) + (-g*tr+b*ti)/tm^2*(vi_fr*vr_to - vr_fr*vi_to) )
end

""
function expression_branch_power_ohms_yt_from_goc(pm::_PM.AbstractACRModel, n::Int, f_bus, t_bus, f_idx, t_idx, g, b, g_fr, b_fr, tr, ti, tm)
    vr_fr = var(pm, n, :vr, f_bus)
    vr_to = var(pm, n, :vr, t_bus)
    vi_fr = var(pm, n, :vi, f_bus)
    vi_to = var(pm, n, :vi, t_bus)

    var(pm, n, :p)[f_idx] =  (g/tm^2+g_fr)*(vr_fr^2 + vi_fr^2) + (-g*tr+b*ti)/tm^2*(vr_fr*vr_to + vi_fr*vi_to) + (-b*tr-g*ti)/tm^2*(vi_fr*vr_to - vr_fr*vi_to)
    var(pm, n, :q)[f_idx] = -(b/tm^2+b_fr)*(vr_fr^2 + vi_fr^2) - (-b*tr-g*ti)/tm^2*(vr_fr*vr_to + vi_fr*vi_to) + (-g*tr+b*ti)/tm^2*(vi_fr*vr_to - vr_fr*vi_to)
end

""
function expression_bus_withdrawal(pm::_PM.AbstractACRModel, n::Int, i::Int, bus_storage, bus_pd, bus_qd, bus_gs, bus_bs)
    vr = var(pm, n, :vr, i)
    vi = var(pm, n, :vi, i)
    ps = get(var(pm, n), :ps, Dict()); _PM._check_var_keys(ps, bus_storage, "active power", "storage")
    qs = get(var(pm, n), :qs, Dict()); _PM._check_var_keys(ps, bus_storage, "reactive power", "storage")

    ps_total = 0.0
    qs_total = 0.0
    if length(bus_storage) > 0
        ps_total = sum(ps[s] for s in bus_storage)
        qs_total = sum(qs[s] for s in bus_storage)
    end

    pd_total = 0.0
    qd_total = 0.0
    if length(bus_pd) > 0
        pd_total = sum(pd for pd in values(bus_pd))
    end
        if length(bus_qd) > 0
        qd_total = sum(qd for qd in values(bus_qd))
    end

    gs_total = 0.0
    bs_total = 0.0
    if length(bus_gs) > 0
        gs_total = sum(gs for gs in values(bus_gs))*(vr^2 + vi^2)
    end
    if length(bus_bs) > 0
        bs_total = sum(bs for bs in values(bus_bs))*(vr^2 + vi^2)
    end

    var(pm, n, :bus_wdp)[i] = ps_total + pd_total + gs_total
    var(pm, n, :bus_wdq)[i] = qs_total + qd_total + bs_total
end


text:
using LabelledArrays

export rossler_rossler

"""
    eom_rossler_rossler(u, p, t) → Function

Equations of motions for a 6D dynamical system consisting
of two 3D Rössler attractors.
"""
function eom_rossler_rossler(u, p, t)
    a, b, c, ϵ₁, ϵ₂, ω₁, ω₂ = ([p[i] for i = 1:7]...,)
    x₁, y₁, z₁, x₂, y₂, z₂ = ([u[i] for i = 1:6]...,)
    # First Rössler system
    dx₁ = ω₁*(-y₁) - z₁ + ϵ₂*(x₁ - x₂)
    dy₁ = ω₁*x₁ + a*y₁
    dz₁ = b + z₁*(x₁ - c)

    # Second Rössler system
    dx₂ = ω₂*(-y₂) - z₂ + ϵ₁*(x₂ - x₁)
    dy₂ = ω₂*x₂ + a*y₂
    dz₂ = b + z₂*(x₂ - c)
    return SVector{6}(dx₁, dy₁, dz₁, dx₂, dy₂, dz₂)
end

"""
	rossler_rossler(u₀, a, b, c, ϵ₁, ϵ₂, ω₁, ω₂) → ContinuousDynamicalSystem

A coupled 6D Rossler system from Krakovská et al. (2018).
The system consists of two separate subsystems, each being a 3D Rossler
attractor. The subsystems are bidirectionally coupled, influencing each other
through variables `x₁` and `x₂`.

The parameters and default values (if calling the keyword version of `rrbi`)
are as follows:

```
a = 0.1      # parameter in the Rössler attractor
b = 0.1      # parameter in the Rössler attractor
c = 14.0     # parameter in the Rössler attractor
ϵ₁ = 0.0     # influence of subsystem 1 on subsystem 2
ϵ₂ = 0.0     # influence of subsystem 2 on subsystem 1
ω₁ = 1+0.015 # the frequency of the first system
ω₂ = 1-0.015 # the frequency of the second system
```

The dynamics is generated by the vector field:

```math
\\begin{aligned}
\\dot x₁ = ω₁ (-y₁) - z₁ + ϵ₂(x₁ - x₂)
\\dot y₁ = ω₁ x₁ + ay₁
\\dot z₁ = b + z₁(x₁ - c)

\\dot x₂ = ω₂(-y₂) - z₂ + ϵ₁(x₂ - x₁)
\\dot y₂ = ω₂x₂ + ay₂
\\dot z₂ = b + z₂(x₂ - c)
\\end{aligned}
```

with the coupling constant ``c \\geq 0``.

# References
Krakovská, A., Jakubík, J., Chvosteková, M., Coufal, D., Jajcay, N., & Paluš, M. (2018). Comparison of six methods for the detection of causality in a bivariate time series. Physical Review E, 97(4), 042207.
"""
function rossler_rossler(u₀, a, b, c, ϵ₁, ϵ₂, ω₁, ω₂)
    p = @LArray [a, b, c, ϵ₁, ϵ₂, ω₁, ω₂] (:a, :b, :c, :ϵ₁, :ϵ₂, :ω₁, :ω₂)
    return ContinuousDynamicalSystem(eom_rossler_rossler, u₀, p)
end

"""
	rossler_rossler(;u₀ = rand(6), a = 0.1, b = 0.1, c = 14.0, ϵ₁ = 0.0,
		ϵ₂ = 0.0, ω₁ = 1 + 0.015, ω₂ = 1 - 0.015) → ContinuousDynamicalSystem

A coupled 6D Rossler system from Krakovská et al. (2018).
The system consists of two separate subsystems, each being a 3D Rossler
attractor. The subsystems are bidirectionally coupled, influencing each other
through variables `x₁` and `x₂`.

The parameters and default values (if calling the keyword version of `rrbi`)
are as follows:

```
a = 0.1      # parameter in the Rössler attractor
b = 0.1      # parameter in the Rössler attractor
c = 14.0     # parameter in the Rössler attractor
ϵ₁ = 0.0     # influence of subsystem 1 on subsystem 2
ϵ₂ = 0.0     # influence of subsystem 2 on subsystem 1
ω₁ = 1+0.015 # the frequency of the first system
ω₂ = 1-0.015 # the frequency of the second system
```

The dynamics is generated by the vector field:

```math
\\begin{aligned}
dx₁ = ω₁*(-y₁) - z₁ + ϵ₂*(x₁ - x₂)
dy₁ = ω₁*x₁ + a*y₁
dz₁ = b + z₁*(x₁ - c)

dx₂ = ω₂*(-y₂) - z₂ + ϵ₁*(x₂ - x₁)
dy₂ = ω₂*x₂ + a*y₂
dz₂ = b + z₂*(x₂ - c)
\\end{aligned}
```

with the coupling constant ``c \\geq 0``.

# References
Krakovská, A., Jakubík, J., Chvosteková, M., Coufal, D., Jajcay, N., & Paluš, M. (2018). Comparison of six methods for the detection of causality in a bivariate time series. Physical Review E, 97(4), 042207.
"""
rossler_rossler(;u₀ = rand(6), a = 0.1, b = 0.1, c = 14.0,
            ϵ₁ = 0.0, ϵ₂ = 0.0, ω₁ = 1 + 0.015, ω₂ = 1 - 0.015) =
    rossler_rossler(u₀, a, b, c, ϵ₁, ϵ₂, ω₁, ω₂)


text:
# signal_dft.jl
##########################
fs=Int(1000);
Ts=1/fs;
N=Int(12000);
T=N/fs;
#t=(0:Int(N-1))*Ts;

x=zeros(Float64, N)
for m=1:N
    x[m]=sin(2*Ts*(m-1)*(Ts*(m-1)-3)*(Ts*(m-1)-6)*(Ts*(m-1)-9)*(Ts*(m-1)-12))
end

using FFTW, DSP

hatx=fft(x);
k = (0:length(hatx)-1)*50/length(hatx);
G=broadcast(abs, hatx);

using Plots

plot(G,
    xaxis="k = 0, 1, 2, ..., N-1",
    yaxis="|F[x](k)|",
    title="DFT of x(t)=sin(2t(t-3)(t-6)(t-9)(t-12))",
    grid=false,
    legend=false)
plot!(xticks = ([0 5000 10000;], [0, 5000, 10000]))
plot!(yticks = ([0 200 400;], [0 200 400]))
savefig("signal_dft.png")


text:
using SmoothingSplines

function includeall(dir)
    for (root,dirs,files) in walkdir(dir)
        for file in files
            if file[end-2:end] == ".jl" && file[1:4]!="Test"
                include(joinpath(root,file))
            end
        end
    end
end

@inline function evalpoly(x::N,coefs::T) where {T<:AbstractArray,N<:Number}
    out = 0.0
    for i in length(coefs):-1:1
        @inbounds out += coefs[i]
        if i != 1
            @fastmath out *= x
        end
    end
    return out
end

export evalpoly

@inline function getBoundingIndsSorted(el::Q,x::T) where {T<:AbstractArray,Q<:Any}
    if el <= x[1]
        return (1,-1)
    end
    for i in 1:length(x)
        if @inbounds x[i] >= el
            if @inbounds x[i] == el
                return (i,-1)
            else
                return (i-1,i)
            end
        end
    end
    return (length(x),-1)
end
export getBoundingIndsSorted

"""
fit a cubic spline to data and return a function evaluating that spline
"""
function getspline(xs,vals;s=1e-10)
    smspl = fit(SmoothingSpline,xs,vals,s)
    F(x::T) where {T} = _predict(smspl,x)
    return F
end

# obtained from SmoothingSpline.jl and modified the input typing for ReverseDiff TrackedReal
function _predict(spl::SmoothingSpline{T}, x::T1) where {T<:Union{Float32, Float64},T1}
    n = length(spl.Xdesign)
    idxl = searchsortedlast(spl.Xdesign, x)
    idxr = idxl + 1
    if idxl == 0 # linear extrapolation to the left
        gl = spl.g[1]
        gr = spl.g[2]
        γ  = spl.γ[1]
        xl = spl.Xdesign[1]
        xr = spl.Xdesign[2]
        gprime = (gr-gl)/(xr-xl) - 1/6*(xr-xl)*γ
        val = gl - (xl-x)*gprime
    elseif idxl == n # linear extrapolation to the right
        gl = spl.g[n-1]
        gr = spl.g[n]
        γ  = spl.γ[n-2]
        xl = spl.Xdesign[n-1]
        xr = spl.Xdesign[n]
        gprime = (gr-gl)/(xr-xl) +1/6*(xr-xl)*γ
        val = gr + (x - xr)*gprime
    else # cubic interpolation
        xl = spl.Xdesign[idxl]
        xr = spl.Xdesign[idxr]
        γl = idxl == 1 ? zero(T) : spl.γ[idxl-1]
        γr = idxl == n-1 ? zero(T) : spl.γ[idxr-1]
        gl = spl.g[idxl]
        gr = spl.g[idxr]
        h = xr-xl
        val = ((x-xl)*gr + (xr-x)*gl)/h
        val -=  1/6*(x-xl)*(xr-x)*((1 + (x-xl)/h)*γr + (1+ (xr-x)/h)*γl)
    end
    val
end


text:
# Correlation
# TODO: Implement `CorrelationTest` for two arguments

using Statistics: clampcor

export CorrelationTest

"""
    CorrelationTest(x, y)

Perform a t-test for the hypothesis that ``\\text{Cor}(x,y) = 0``, i.e. the correlation
of vectors `x` and `y` is zero.

    CorrelationTest(x, y, Z)

Perform a t-test for the hypothesis that ``\\text{Cor}(x,y|Z=z) = 0``, i.e. the partial
correlation of vectors `x` and `y` given the matrix `Z` is zero.

Implements `pvalue` for the t-test.
Implements `confint` using an approximate confidence interval based on Fisher's
``z``-transform.

See also `partialcor` from StatsBase.

# External resources

* [Partial correlation on Wikipedia](https://en.wikipedia.org/wiki/Partial_correlation#As_conditional_independence_test)
  (for the construction of the confidence interval)
* [Section testing using Student's t-distribution from Pearson correlation coefficient on Wikipedia](https://en.wikipedia.org/wiki/Pearson_correlation_coefficient#Testing_using_Student's_t-distribution)
* [K.J. Levy and S.C. Narula (1978): Testing Hypotheses concerning Partial Correlations: Some Methods and Discussion. International Statistical Review 46(2).](https://www.jstor.org/stable/1402814)
"""
struct CorrelationTest{T<:Real} <: HypothesisTest
    r::T
    n::Int
    k::Int
    t::T

    # Error checking is done in `cor`
    function CorrelationTest(x::AbstractVector, y::AbstractVector)
        r = cor(x, y)
        n = length(x)
        t = r * sqrt((n - 2) / (1 - r^2))
        return new{typeof(r)}(r, n, 0, t)
    end

    # Error checking is done in `partialcor`
    function CorrelationTest(x::AbstractVector, y::AbstractVector, Z::AbstractMatrix)
        r = partialcor(x, y, Z)
        n, k = size(Z)
        t = r * sqrt((n - 2 - k) / (1 - r^2))
        return new{typeof(r)}(r, n, k, t)
    end

    function CorrelationTest(x::AbstractVector, y::AbstractVector, z::AbstractVector)
        r = partialcor(x, y, z)
        n = length(z)
        t = r * sqrt((n - 3) / (1 - r^2))
        return new{typeof(r)}(r, n, 1, t)
    end
end

testname(p::CorrelationTest) =
    string("Test for nonzero ", p.k != 0 ? "partial " : "", "correlation")

function population_param_of_interest(p::CorrelationTest)
    param = p.k != 0 ? "Partial correlation" : "Correlation"
    (param, zero(p.r), p.r)
end

StatsBase.nobs(p::CorrelationTest) = p.n
StatsBase.dof(p::CorrelationTest) = p.n - 2 - p.k

function StatsBase.confint(test::CorrelationTest{T}, level::Float64=0.95) where T
    dof(test) > 1 || return (-one(T), one(T))  # Otherwise we can get NaNs
    q = quantile(Normal(), 1 - (1-level) / 2)
    fisher = atanh(test.r)
    bound = q / sqrt(dof(test) - 1)
    elo = clampcor(tanh(fisher - bound))
    ehi = clampcor(tanh(fisher + bound))
    return (elo, ehi)
end

default_tail(::CorrelationTest) = :both
pvalue(test::CorrelationTest; tail=:both) = pvalue(TDist(dof(test)), test.t, tail=tail)

function show_params(io::IO, test::CorrelationTest, indent="")
    println(io, indent, "number of observations:          ", nobs(test))
    println(io, indent, "number of conditional variables: ", test.k)
    println(io, indent, "t-statistic:                     ", test.t)
    println(io, indent, "degrees of freedom:              ", dof(test))
end


text:
#=
Copyright (c) 2015, Intel Corporation

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

    * Redistributions of source code must retain the above copyright notice,
      this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of Intel Corporation nor the names of its contributors
      may be used to endorse or promote products derived from this software
      without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=#

@doc """ The whole function region. Not really used so far. """
type FunctionRegion <: Region
    func_ast         :: Expr
    members          :: Set{BasicBlockIndex} # All members in the function
    entry            :: Any #BasicBlock
    exit             :: Any #BasicBlock
    symbol_property  :: Symexpr2PropertiesMap
    property_proxies :: Any

    FunctionRegion(_func_ast) = new(_func_ast, Set{BasicBlockIndex}(), nothing,
        nothing, Symexpr2PropertiesMap(), nothing)
end

@doc """
An exit edge of a loop from a block in the loop to another block outside the loop.
"""
type LoopExit
    from_bb :: BasicBlock
    to_bb   :: BasicBlock
end

@doc """
A region that is composed of a loop. Speeding up loops is the focus of
Sparse Accelerator. A loop region may have a parent region (the immediate
enclosing outer loop region, or the function region), and may have some exits,
and immediate members (the indexes of the basic blocks that belong to the
loop, but not to any inner loop of it).
"""
type LoopRegion <: Region
    parent            :: Region
    loop              :: Loop
    exits             :: Set{LoopExit}
    symbol_property   :: Symexpr2PropertiesMap
    property_proxies  :: Any
    immediate_members :: Set{BasicBlockIndex}
end

@doc """
Form a region with the loop.
"""
function loop_region_formation(
    parent  :: Region,
    L       :: Loop,
    cfg     :: CFG,
    members :: Set{BasicBlockIndex}
)
    region = LoopRegion(parent, L, Set{LoopExit}(), Symexpr2PropertiesMap(), nothing, members)
    blocks = cfg.basic_blocks
    for bb_index in L.members
        bb = blocks[bb_index]
        for succ in bb.succs
            if !in(succ.label, L.members)
                exit = LoopExit(bb, succ)
                push!(region.exits, exit)
            end
        end
    end
    region
end

@doc """
Form regions for all the outermost loops in the control flow graph.
"""
function loop_region_formation(
    parent    :: Region,
    cfg       :: CFG,
    loop_info :: DomLoops
)
    regions = Vector{LoopRegion}()
    for L in loop_info.loops
        is_outermost      = true
        members = copy(L.members)
        for L1 in loop_info.loops
            assert(L == L1 || L.members != L1.members)
            if L != L1 && L.members < L1.members
                is_outermost = false
                break
            end
            if L != L1 && L1.members < L.members
                setdiff!(members, L1.members)
            end
        end

        if is_outermost
            region = loop_region_formation(parent , L, cfg, members)
            push!(regions, region)
        end
    end
    return regions
end

@doc """
Form regions.
"""
function map_blocks_to_loop_nesting_depths(
    cfg       :: CFG,
    loop_info :: DomLoops
)
    bb2depth = Dict{BasicBlock, Int}()
    for (bb_idx, bb) in cfg.basic_blocks
        bb2depth[bb] = 0
    end
    for L in loop_info.loops
        members = copy(L.members)
        depth   = 1
        for L1 in loop_info.loops
            if L != L1 && L.members < L1.members
                depth += 1
                setdiff!(members, L1.members)
            end
        end
        for bb_idx in members
            bb = cfg.basic_blocks[bb_idx]
            bb2depth[bb] = depth
        end
    end
    bb2depth
end

@doc """
Form regions.
"""
function region_formation(
    func_region :: FunctionRegion,
    cfg         :: CFG,
    loop_info   :: DomLoops
)
    for (bb_idx, bb) in cfg.basic_blocks
        push!(func_region.members, bb_idx)
        if bb_idx == -1
            func_region.entry = bb
        end
        if bb_idx == -2
            func_region.exit = bb
        end
    end

    # So far, form only loop regions.
    # TODO: (1) Connect multiple loop regions together into a bigger region
    #       (2) Extend a region to include non-loop code
    regions = loop_region_formation(func_region, cfg, loop_info)

    # Create a map from blocks to loop nesting depth. A block outside any loop has
    # a depth of 0.
    bb2depth = map_blocks_to_loop_nesting_depths(cfg, loop_info)

    return regions, bb2depth
end

max_stars_repo_name: IntelLabs/Sparso

max_issues_repo_name: IntelLabs/Sparso

max_forks_repo_name: IntelLabs/Sparso


text:
function compute(first::Int, modulo::Int)::Int
    min_value = max_value = result = first
    candidate = 0
    while candidate ≠ 1
        candidate = (min_value + max_value) % modulo
        if candidate > max_value
            max_value = candidate
        elseif candidate < min_value
            min_value = candidate
            result += candidate
        end
    end
    return result
end

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
using .NavAbilityCaesarExtensions
using NavAbilitySDK
using DistributedFactorGraphs
using RoME
using JSON
using Unmarshal

# TODO: Reconcile these differences
function gqlVariableToDfg(dfg, gqlVariable)
  if haskey(gqlVariable,"label")
    ppes = map(ppe -> begin
      return unmarshal(MeanMaxPPE, ppe)
    end, gqlVariable["ppes"])
    solverData = map(sd -> begin
      packedSolverData = unmarshal(PackedVariableNodeData, sd)
      return unpackVariableNodeData(dfg, packedSolverData)
    end,gqlVariable["solverData"])
  end
  dfgVariable = unpackVariable(dfg, gqlVariable, unpackPPEs=false, unpackSolverData=false, unpackBigData=false)
  [(dfgVariable.solverDataDict[sd.solveKey] = sd) for sd in solverData]
  [(dfgVariable.ppeDict[ppe.solveKey] = ppe) for ppe in ppes]
  return dfgVariable
end

function addVariable!(dfg::NavAbilityDFG, variable::AbstractDFGVariable; maxWaitSeconds::Int=150)
  context = Client(dfg.userId,dfg.robotId,dfg.sessionId)
  packedVariable = packVariable(dfg, variable)
  requestId = addPackedVariable(dfg.navabilityClient, context, packedVariable)
  label = string(getLabel(variable))
  for i in 1:maxWaitSeconds
    @info i
    if exists(dfg, label)
      break
    end
    sleep(1)
  end
  return requestId
end

function getVariable(dfg::NavAbilityDFG, label::Union{Symbol, String})
  context = Client(dfg.userId,dfg.robotId,dfg.sessionId)
  gqlVariable = NavAbilitySDK.getVariable(dfg.navabilityClient, context, string(label)) |> fetch
  dfgVariable = gqlVariableToDfg(dfg, gqlVariable)
  return dfgVariable
end

function getVariables(dfg::NavAbilityDFG, regexFilter::Union{Nothing, Regex}=nothing; tags::Vector{Symbol}=Symbol[], solvable::Int=0)
  context = Client(dfg.userId,dfg.robotId,dfg.sessionId)
  gqlVariables = NavAbilitySDK.getVariables(dfg.navabilityClient, context) |> fetch
  dfgVariables = map(v -> gqlVariableToDfg(dfg, v),gqlVariables)
  # TODO: Implement regexFilter
  # TODO: Implement tags
  # TODO: Implement solver filter
  return dfgVariables
end