text:
function mt_net_comparison(
      parent_folder::String,
      output_folder::String,
      objective::String;
      gen_areas::Array=[],
      area_transfer::Array=[],
      gen_zones::Array=[],
      zone_transfer::Array=[],
      enum_branches::Bool=false,
      plot_best_x::Int64=-1,
      print_results::Bool=false,
      results_folders::Array=[],
      scopf::Bool=false,
      no_converter_loss::Bool=false,
      output_to_files::Bool=true,
      regularize_f::Float64=0.0
   )

   # Run the opf for the base case then all dc configurations
   results_dict_allplots = []
   n_subnets = 0
   subnet_array = []
   idx_sorted = []
   series_labels = ["control P, Q, X(ω)"]
   if length(results_folders) == 0
      # Base case opf
      (results_dict,n_subnets,subnet_array,idx_sorted,series_output_folder,plot_best_x) = run_series(
         parent_folder,
         objective;
         gen_areas,
         area_transfer,
         gen_zones,
         zone_transfer,
         enum_branches,
         plot_best_x,
         print_results,
         suffix="",
         no_converter_loss=no_converter_loss,
         output_to_files=output_to_files,
         regularize_f=regularize_f
      )
      push!(results_dict_allplots, results_dict)

      results_dict_string = JSON.json(results_dict)
      println("Saving results at $(series_output_folder)/all_results_dict.json")
      open("$series_output_folder/all_results_dict.json", "w") do f
         write(f, results_dict_string)
      end
      push!(results_folders, series_output_folder)
   else
      # Use results data that has already been saved
      results_dict_allplots = []
      for (i,output_folder) in enumerate(results_folders)

         results_dict_parsed = JSON.parsefile("$(output_folder)/all_results_dict.json")
         if i == 1
            (n_subnets, subnet_array, idx_sorted) = get_sort_indices(results_dict_parsed, objective, gen_areas, area_transfer, gen_zones, zone_transfer)
         end

         push!(results_dict_allplots, results_dict_parsed)
      end
   end
   if output_to_files
      println("Saving results at $(output_folder)/results_dict_allplots.json")
      results_dict_string = JSON.json(results_dict_allplots)
      if !isdir(output_folder)
         mkpath(output_folder)
      end
      open("$(output_folder)/results_dict_allplots.json", "w") do f
         write(f, results_dict_string)
      end
   end
   plot_results_dicts_bar(results_dict_allplots,n_subnets,subnet_array,idx_sorted,output_folder,plot_best_x,series_labels)

end

function multi_folder(
      parent_folders::Array,
      output_folder::String,
      objective::String;
      gen_areas::Array=[],
      area_transfer::Array=[],
      gen_zones::Array=[],
      zone_transfer::Array=[],
      enum_branches::Bool=false,
      plot_best_x::Int64=-1,
      print_results::Bool=false,
      series_labels::Array=[],
      regularize_f::Float64=0.0
   )

   results_dict_allplots = []
   n_subnets = 0
   subnet_array = []
   idx_sorted = []
   for (i,parent_folder) in enumerate(parent_folders)
      (results_dict,n_subnets,subnet_array,idx_sorted,series_output_folder,plot_best_x) = run_series(
         parent_folder,
         objective;
         gen_areas,
         area_transfer,
         gen_zones,
         zone_transfer,
         enum_branches,
         plot_best_x,
         print_results,
         regularize_f=regularize_f
      )
      push!(results_dict_allplots, results_dict)

      results_dict_string = JSON.json(results_dict)
      println("Saving results at $(series_output_folder)/all_results_dict.json")
      open("$series_output_folder/all_results_dict.json", "w") do f
         write(f, results_dict_string)
      end
   end

   results_dict_allplots = []
   for (i,parent_folder) in enumerate(parent_folders)
      folder_split = split(rstrip(parent_folder,'/'),"/")
      folder_split = folder_split[folder_split .!= ""]
      toplevels = folder_split[1:end-2]
      output_folder = join([("$i/") for i in toplevels])*"results/$(folder_split[end])"

      results_dict_parsed = JSON.parsefile("$(output_folder)/all_results_dict.json")

      if i==1
         (n_subnets, subnet_array, idx_sorted) = get_sort_indices(results_dict_parsed, objective, gen_areas, area_transfer, gen_zones, zone_transfer)
      end

      push!(results_dict_allplots, results_dict_parsed)
   end

   plot_results_dicts_bar(results_dict_allplots,n_subnets,subnet_array,idx_sorted,output_folder,plot_best_x,series_labels)
end

function run_series(
      parent_folder::String,
      objective::String;
      gen_areas::Array=[],
      area_transfer::Array=[],
      gen_zones::Array=[],
      zone_transfer::Array=[],
      enum_branches::Bool=false,
      plot_best_x::Int64=-1,
      print_results::Bool=false,
      override_param::Dict=Dict(),
      suffix::String="",
      fix_f_override::Bool=false,
      direct_pq::Bool=true,
      master_subnet::Int64=1,
      no_converter_loss::Bool=false,
      output_to_files::Bool=false,
      start_vals=Dict{String, Dict}("sn"=>Dict()),
      regularize_f::Float64=0.0
   )
   folders = [f for f in readdir(parent_folder) if (isdir("$parent_folder/$f"))]

   results_dict = Dict{String, Dict{String, Any}}()

   summary_df = DataFrame()
   folder_split = split(rstrip(parent_folder,'/'),"/")
   folder_split = folder_split[folder_split .!= ""]
   toplevels = folder_split[1:end-2]
   output_folder = join([("$i/") for i in toplevels])*"results/$(folder_split[end])$suffix"
   # output_folder = join([(i=="data" ? "results/" : "$i/") for i in split(parent_folder,"/")])
   println("output location: $output_folder")
   if !isdir(output_folder)
      mkpath(output_folder)
   end
   start_vals = Dict{String, Dict}("sn"=>Dict())
   subnet_array = []
   n_subnets = 0

   # Choose which label to use for sorting all the results
   if objective == "mincost"
      first_label = "cost"
   elseif objective == "areagen"
      n_areas = length(gen_areas)
      if n_areas > 1
         gen_areas_sorted = sort(gen_areas)
         first_label = "generation in areas"
         for i in 1:(n_areas-2)
            first_label *= " $(gen_areas_sorted[i]),"
         end
         first_label *=  " $(gen_areas_sorted[end-1]) and $(gen_areas_sorted[end]) (p.u.)"
      else
         first_label = "generation in areas $(gen_areas[1]) (p.u.)"
      end
   elseif objective == "zonegen"
      n_zones = length(gen_zones)
      if n_zones > 1
         gen_zones_sorted = sort(gen_zones)
         first_label = "generation in zones"
         for i in 1:(n_zones-2)
            first_label *= " $(gen_zones_sorted[i]),"
         end
         first_label *=  " $(gen_zones_sorted[end-1]) and $(gen_zones_sorted[end]) (p.u.)"
      else
         first_label = "generation in zones $(gen_zones[1]) (p.u.)"
      end
   end

   for (i,folder) in enumerate(folders)
      println("running $parent_folder/$(folders[i])  |  $i/$(length(folders)) folders")
      override_param_tmp = deepcopy(override_param)
      (result_temp, res_summary_temp, solution_pm_temp, results_dict_temp) = run_opf(
         parent_folder, folders[i],
         objective, gen_areas, area_transfer, gen_zones, zone_transfer, results_dict,
         fix_f_override,
         direct_pq, master_subnet,
         suffix, print_results,
         override_param=override_param_tmp,
         no_converter_loss=no_converter_loss,
         output_to_files=output_to_files,
         start_vals=start_vals,
         regularize_f=regularize_f
      )

      # println("STATUS: $(res_summary_temp.status[1])")
      if !(res_summary_temp.status[1] in [LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED])
         fix_f_init = false

         println("Running OPF again without converter losses.")
         # Now run the opf again without converter losses by applying the noloss_override dictionary
         (result_init, res_summary_init, solution_pm_init, results_dict_init) = run_opf(
            parent_folder, folders[i],
            objective, gen_areas, area_transfer, gen_zones, zone_transfer, results_dict,
            fix_f_override,
            direct_pq, master_subnet,
            suffix, print_results,
            override_param=override_param_tmp,
            no_converter_loss=no_converter_loss,
            output_to_files=output_to_files,
            start_vals=start_vals,
            regularize_f=regularize_f
         )
         if (res_summary_init.status[1] in [LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED])
            println("Using no-loss solution as initialization, running OPF again with converter losses.")
            # Run the opf again with converter losses and variable frequency, now using the solution from the no-loss case as initialization
            (result_temp, res_summary_temp, solution_pm_temp, results_dict_temp) = run_opf(
               parent_folder, folders[i],
               objective, gen_areas, area_transfer, gen_zones, zone_transfer, results_dict,
               fix_f_override,
               direct_pq, master_subnet,
               suffix, print_results,
               override_param=override_param_tmp,
               start_vals=solution_pm_init,
               no_converter_loss=no_converter_loss,
               output_to_files=output_to_files,
               regularize_f=regularize_f
            )
         end
      end
      result = result_temp
      res_summary = res_summary_temp
      solution_pm = solution_pm_temp
      results_dict = results_dict_temp

      res_summary[:,:network] .= rstrip(folders[i],'/')*suffix
      if i == 1
         summary_df = res_summary
      else
         summary_df = vcat(summary_df, res_summary)
      end

   end
   label = first_label
   results = results_dict[label]

   br_results = Array{Float64,1}()
   br_indices = Array{String,1}()
   base_res = NaN

   for (folder,res) in results
      subnet_array = union(subnet_array, results_dict["subnet"][folder]) # Get maximal set with subnets represented in any folders
      if occursin("br", folder)
         push!(br_results, sum(res))
         push!(br_indices, folder)
      elseif occursin("base", folder)
         base_res = res
         println("set base_res to $base_res")
      end
   end

   br_results_tosort = copy(br_results)
   br_results_tosort[isnan.(br_results_tosort) .& (br_results_tosort.==-9999.0)] .= Inf
   perm = sortperm(br_results_tosort)
   idx_sorted = br_indices[perm]

   if output_to_files
      results_dict_string = JSON.json(results_dict)
      println("Saving results at $(output_folder)/all_results_dict.json")
      open("$output_folder/all_results_dict.json", "w") do f
         write(f, results_dict_string)
      end

      println("Saving summary at $(output_folder)/summary.csv")
      CSV.write("$(output_folder)/summary.csv", summary_df)
   end

   n_subnets = length(subnet_array)

   return results_dict, n_subnets, subnet_array, idx_sorted, output_folder, plot_best_x

end

function get_sort_indices(results_dict, objective, gen_areas::Array, area_transfer::Array, gen_zones::Array, zone_transfer::Array)
   # Choose which label to use for sorting all the results
   if objective == "mincost"
      first_label = "cost"
   elseif objective == "areagen"
      n_areas = length(gen_areas)
      if n_areas > 1
         gen_areas_sorted = sort(gen_areas)
         first_label = "generation in areas"
         for i in 1:(n_areas-2)
            first_label *= " $(gen_areas_sorted[i]),"
         end
         first_label *=  " $(gen_areas_sorted[end-1]) and $(gen_areas_sorted[end]) (p.u.)"
      else
         first_label = "generation in areas $(gen_areas[1]) (p.u.)"
      end
   elseif objective == "zonegen"
      n_zones = length(gen_zones)
      if n_zones > 1
         gen_zones_sorted = sort(gen_zones)
         first_label = "generation in zones"
         for i in 1:(n_zones-2)
            first_label *= " $(gen_zones_sorted[i]),"
         end
         first_label *=  " $(gen_zones_sorted[end-1]) and $(gen_zones_sorted[end]) (p.u.)"
      else
         first_label = "generation in zones $(gen_zones[1]) (p.u.)"
      end
   end

   label = first_label
   results = results_dict[label]

   br_results = Array{Float64,1}()
   br_indices = Array{String,1}()
   base_res = NaN

   subnet_array = []
   for (folder,res) in results
      # if results_dict["status"][folder] in [LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED]
         subnet_array = union(subnet_array, results_dict["subnet"][folder]) # Get maximal set with subnets represented in any folders
         if occursin("br", folder)
            push!(br_results, sum(res))
            push!(br_indices, folder)
         elseif occursin("base", folder)
            base_res = res
            println("set base_res to $base_res")
         end
      # end
   end

   br_results_tosort = copy(br_results)
   br_results_tosort[br_results_tosort.==-9999.0] .= Inf
   perm = sortperm(br_results_tosort)
   idx_sorted = br_indices[perm]

   n_subnets = length(subnet_array)

   return n_subnets, subnet_array, idx_sorted
end


text:
save_output = false

using DSGE, ModelConstructors, Test, JLD2, FileIO, OrderedCollections, CSV, Random, DataFrames
path = dirname(@__FILE__)

m = AnSchorfheide()

m <= Setting(:use_population_forecast, false)
m <= Setting(:cond_id, 0)
m <= Setting(:date_forecast_start, quartertodate("2015-Q4"))
m <= Setting(:date_conditional_end, quartertodate("2015-Q4"))
m <= Setting(:forecast_block_size, 5)
m <= Setting(:data_vintage, "REF", true, "vint", "Date of data")
m <= Setting(:dataroot, "$path/../reference/input_data")

estroot = normpath(joinpath(dirname(@__FILE__), "..", "reference"))
overrides = forecast_input_file_overrides(m)
overrides[:mode] = joinpath(estroot, "optimize.h5")
overrides[:full] = joinpath(estroot, "metropolis_hastings.h5")

output_vars = add_requisite_output_vars([:histpseudo, :histobs, :histstdshocks,
                                         :histutpseudo, :histutobs,
                                         :hist4qpseudo, :hist4qobs,
                                         :forecaststates, :forecastpseudo, :forecastobs, :forecaststdshocks,
                                         :forecastutpseudo, :forecastutobs,
                                         :forecast4qpseudo, :forecast4qobs,
                                         :bddforecaststates, :bddforecastshocks, :bddforecastpseudo, :bddforecastobs,
                                         :shockdecpseudo, :shockdecobs,
                                         :trendstates, :trendobs, :trendpseudo,
                                         :dettrendstates, :dettrendobs, :dettrendpseudo,
                                         :irfstates, :irfpseudo, :irfobs])

Random.seed!(47)
df = VERSION >= v"1.3" ? CSV.read("$path/../reference/input_data/data/data_dsid=00_vint=REF.csv", DataFrame) :
    CSV.read("$path/../reference/input_data/data/data_dsid=00_vint=REF.csv")
forecast_one(m, :full, :none, output_vars, df = df, verbose = :none)
compute_meansbands(m, :full, :none, [:histobs, :forecastobs, :histpseudo, :forecastpseudo,
                                     :hist4qobs, :forecast4qobs, :hist4qpseudo, :forecast4qpseudo, :shockdecobs, :irfobs], df = df, verbose = :none)

hist, fore = construct_fcast_and_hist_dfs(m, :none, [:obs_gdp, :obs_cpi, :obs_nominalrate])
hist_noT, fore_noT = construct_fcast_and_hist_dfs(m, :none, [:obs_gdp, :obs_cpi, :obs_nominalrate], include_T_in_df_forecast = false)
hist_4q, fore_4q = construct_fcast_and_hist_dfs(m, :none, [:obs_gdp, :obs_cpi, :obs_nominalrate], use_4q = true)

if save_output
    fn = VERSION >= v"1.5" ? "$path/../reference/df_to_table_out_v1p5.jld2" : "$path/../reference/df_to_table_out.jld2"
    jldopen(fn, true, true, true, IOStream) do file
        file["hist"] = hist
        file["hist_noT"] = hist_noT
        file["hist_4q"] = hist_4q
        file["fore"] = fore
        file["fore_noT"] = fore_noT
        file["fore_4q"] = fore_4q
    end
end

fn = VERSION >= v"1.5" ? "$path/../reference/df_to_table_out_v1p5.jld2" : "$path/../reference/df_to_table_out.jld2"
saved_hist = load(fn, "hist")
saved_hist_noT = load(fn, "hist_noT")
saved_hist_4q = load(fn, "hist_4q")
saved_fore = load(fn, "fore")
saved_fore_noT = load(fn, "fore_noT")
saved_fore_4q = load(fn, "fore_4q")

@testset "Test df_to_table" begin
    @test Matrix(hist) == Matrix(saved_hist)
    @test Matrix(hist_noT) == Matrix(saved_hist_noT)
    @test Matrix(hist_4q) == Matrix(saved_hist_4q)
    @test Matrix(fore) == Matrix(saved_fore)
    @test Matrix(fore_noT) == Matrix(saved_fore_noT)
    @test Matrix(fore_4q) == Matrix(saved_fore_4q)
    @test_throws AssertionError construct_fcast_and_hist_dfs(m, :none, [:obs_gdp], save_to_table = true)
    # Test that saving to tex table functionality runs
    construct_fcast_and_hist_dfs(m, :none, [:obs_gdp], save_to_table = true, table_caption = "Test Caption",
                                 filename = "test.tex", savedir = "$path")
end

@testset "Test auxiliary methods" begin
    header_mappings = DSGE.create_table_header_mappings(m, collect(keys(m.observables)))
    @test header_mappings == Dict(:obs_gdp => Symbol("Real GDP Growth"), :obs_cpi => Symbol("CPI Inflation"), :obs_nominalrate => Symbol("Nominal FFR"))
    #df = load_data(m)
    @test DSGE.unit_mappings(m, df, header_mappings) == OrderedDict(Symbol("Real GDP Growth") =>
                                                             "(Q/Q) \\% Annualized",
                                                             Symbol("CPI Inflation") =>
                                                             "(Q/Q) \\% Annualized")
end

#These really belong in analysis/io.jl MeansBands tests but put them here to avoid duplicate meansbands computations
@testset "Test read_mb" begin
    # just make sure they run
    read_mb(m, :full, :none, :forecastobs)
    read_mb(m, :full, :none, :forecastobs, bdd_and_unbdd = true)
    DSGE.get_mb_metadata(m, :full, :none, :forecastobs)
    DSGE.get_mb_metadata(m, :full, :none, :irfobs)
end

@testset "Test creating q4q4 mb" begin
    mb_4q = read_mb(m, :full, :none, :forecast4qobs)
    @test typeof(create_q4q4_mb(mb_4q)) == MeansBands
end

@testset "Test write_meansbands_tables functions" begin
    @test_throws AssertionError write_meansbands_tables_timeseries(m, :mode, :none, :irf)
    write_meansbands_tables_timeseries(m, :full, :none, :histobs)
    write_meansbands_tables_timeseries(m, :full, :none, :histforecastobs)

    write_means_tables_shockdec(m, :full, :none, :obs)

    write_meansbands_tables_all(m, :full, :none, [:histobs, :shockdecobs])
    @test_throws ErrorException write_meansbands_tables_all(m, :full, :none, [:irfobs])
end

@testset "Test load_posterior_moments works" begin
    m <= Setting(:sampling_method, :MH)
    load_posterior_moments(m)
    moment_tables(m)
    moment_tables(m, use_mode = true)

    m <= Setting(:sampling_method, :SMC)
    m <= Setting(:n_parts, 20)
    m <= Setting(:n_Φ, 10)
    m <= Setting(:adaptive_tempering_target_smc, false)
    data = df_to_matrix(m, load_data(m))
    DSGE.smc2(m, data, verbose = :none, run_csminwel = false)
    load_posterior_moments(m)

end

@testset "Test meansbands_to_matrix works" begin
    meansbands_to_matrix(m, :full, :none, [:histobs])
end

fp = dirname(@__FILE__)
if isfile(joinpath(fp, "test.tex_forecast.tex"))
    rm(joinpath(fp, "test.tex_forecast.tex"))
end
if isfile(joinpath(fp, "test.tex_history.tex"))
    rm(joinpath(fp, "test.tex_history.tex"))
end

ext: jl


text:
using Aqua
using Test
using RecordedArrays
using RecipesBase: apply_recipe
using ArrayInterface

function test_recipe(recipe, args)
    # this hack into the RecipesBase, may break future
    rds = apply_recipe(Dict{Symbol, Any}(), recipe)
    @test length(rds) == length(args)
    for rd in rds
        @test rd.args in args
    end
    return nothing
end

function test_show(x, str)
    @test repr("text/plain", x) == str
end

@testset "RecordedArrays" begin
    @testset "QA" begin
        _params(T::UnionAll) = _params(T.body)
        _params(T::DataType) = T.parameters
        @testset "Ambiguity" begin
            ambiguities = Test.detect_ambiguities(RecordedArrays)
            filter!(ambiguities) do (m1, m2)
                p1 = _params(m1.sig)
                p2 = _params(m2.sig)
                for (t1, t2) in zip(p1, p2)
                    typeintersect(t1, t2) === Union{} && return false
                end
                return true
            end
            if !isempty(ambiguities)
                for amb in ambiguities
                    println(amb[1])
                    println(amb[2])
                end
            end
            @test isempty(ambiguities)
        end
        Aqua.test_all(RecordedArrays; ambiguities=false, deps_compat=false)
    end

    @testset "Utilities" begin
        include("utils.jl")
    end

    @testset "Clock" begin
        include("clock.jl")
    end

    @testset "Entry" begin
        include("entry.jl")
    end

    @testset "RArray" begin
        include("rarray.jl")
    end

    @testset "RNumber" begin
        include("rnumber.jl")
    end
end

# vim:tw=92:ts=4:sw=4:et


text:
function irepeat(v::AbstractVector, n::Int)
    nV = length(v)
    res = similar(v, nV*n)
    @inbounds for j = 1:nV
        vj = v[j]
        base = (j-1)*n
        @inbounds @simd for i = 1:n
            res[base+i] = vj
        end
    end
    res
end

function orepeat(v::AbstractVector, n::Int)
    nV = length(v)
    res = similar(v, nV*n)
    @inbounds for i = 1:n
        base = (i-1)*nV
        @inbounds @simd for j = 1:nV
            res[base+j] = v[j]
        end
    end
    res
end

# TODO: since 0.7 transpose is different, we don't take transpose serious here.
####### kronecker product ###########
# TODO: if IMatrix{1}, do nothing
kron(A::IMatrix{Na, Ta}, B::IMatrix{Nb, Tb}) where {Na, Nb, Ta<:Number, Tb<:Number}= IMatrix{Na*Nb, promote_type(Ta, Tb)}()
kron(A::IMatrix{Na, Ta}, B::IMatrix{1, Tb}) where {Na, Ta<:Number, Tb<:Number}= IMatrix{Na, promote_type(Ta, Tb)}()
kron(A::IMatrix{1, Ta}, B::IMatrix{Nb, Tb}) where {Nb, Ta<:Number, Tb<:Number}= IMatrix{Nb, promote_type(Ta, Tb)}()
kron(A::IMatrix{Na, <:Number}, B::Diagonal{<:Number}) where Na = Diagonal(orepeat(B.diag, Na))
kron(B::Diagonal{<:Number}, A::IMatrix{Na}) where Na = Diagonal(irepeat(B.diag, Na))
for MT in [:AbstractMatrix, :PermMatrix, :SparseMatrixCSC, :Diagonal]
    @eval kron(A::IMatrix{1,<:Number}, B::$MT{<:Number}) = B
    @eval kron(B::$MT{<:Number}, A::IMatrix{1,<:Number}) = B
end

####### diagonal kron ########
kron(A::StridedMatrix{<:Number}, B::Diagonal{<:Number}) = kron(A, PermMatrix(B))
kron(A::Diagonal{<:Number}, B::StridedMatrix{<:Number}) = kron(PermMatrix(A), B)
kron(A::Diagonal{<:Number}, B::SparseMatrixCSC{<:Number}) = kron(PermMatrix(A), B)
kron(A::SparseMatrixCSC{<:Number}, B::Diagonal{<:Number}) = kron(A, PermMatrix(B))

function kron(A::AbstractMatrix{Tv}, B::IMatrix{Nb}) where {Nb, Tv<:Number}
    mA, nA = size(A)
    nzval = Vector{Tv}(undef, Nb*mA*nA)
    rowval = Vector{Int}(undef, Nb*mA*nA)
    colptr = collect(1:mA:Nb*mA*nA+1)
    @inbounds for j = 1:nA
        source = view(A,:,j)
        startbase = (j-1)*Nb*mA - mA
        for j2 = 1:Nb
            start = startbase + j2*mA
            row = j2-Nb
            @inbounds @simd for i = 1:mA
                nzval[start+i] = source[i]
                rowval[start+i] = row+Nb*i
            end
        end
    end
    SparseMatrixCSC(mA*Nb, nA*Nb, colptr, rowval, nzval)
end

function kron(A::IMatrix{Na}, B::AbstractMatrix{Tv}) where {Na, Tv<:Number}
    mB, nB = size(B)
    rowval = Vector{Int}(undef, nB*mB*Na)
    nzval = Vector{Tv}(undef, nB*mB*Na)
    @inbounds for j in 1:Na
        r0 = (j-1)*mB
        for j2 in 1:nB
            start = ((j-1)*nB+j2-1)*mB
            @inbounds @simd for i in 1:mB
                rowval[start+i] = r0+i
                nzval[start+i] = B[i,j2]
            end
        end
    end
    colptr = collect(1:mB:nB*mB*Na+1)
    SparseMatrixCSC(mB*Na, Na*nB, colptr, rowval, nzval)
end

function kron(A::IMatrix{Na}, B::SparseMatrixCSC{T}) where {Na, T<:Number}
    mB, nB = size(B)
    nV = nnz(B)
    nzval = Vector{T}(undef, Na*nV)
    rowval = Vector{Int}(undef, Na*nV)
    colptr = Vector{Int}(undef, nB*Na+1)
    nzval = Vector{T}(undef, Na*nV)
    colptr[1] = 1
    for i = 1:Na
        r0 = (i-1)*mB
        start = nV*(i-1)
        @inbounds @simd for k = 1:nV
            rowval[start+k] = B.rowval[k]+r0
            nzval[start+k] = B.nzval[k]
        end
        colbase = (i-1)*nB
        @inbounds @simd for j=2:nB+1
            colptr[colbase+j] = B.colptr[j]+start
        end
    end
    SparseMatrixCSC(mB*Na, nB*Na, colptr, rowval, nzval)
end

function kron(A::SparseMatrixCSC{T}, B::IMatrix{Nb}) where {T<:Number, Nb}
    mA, nA = size(A)
    nV = nnz(A)
    rowval = Vector{Int}(undef, Nb*nV)
    colptr = Vector{Int}(undef, nA*Nb+1)
    nzval = Vector{T}(undef, Nb*nV)
    z=1
    colptr[1] = 1
    @inbounds for i in 1:nA
        rstart = A.colptr[i]
        rend = A.colptr[i+1]-1
        colbase = (i-1)*Nb+1
        @inbounds for k in 1:Nb
            irow_Nb = k - Nb
            @inbounds @simd for r in rstart:rend
                rowval[z] = A.rowval[r]*Nb+irow_Nb
                nzval[z] = A.nzval[r]
                z+=1
            end
            colptr[colbase+k] = z
        end
    end
    SparseMatrixCSC(mA*Nb, nA*Nb, colptr, rowval, nzval)
end

function kron(A::PermMatrix{T}, B::IMatrix) where T<:Number
    nA = size(A, 1)
    nB = size(B, 1)
    vals = Vector{T}(undef, nB*nA)
    perm = Vector{Int}(undef, nB*nA)
    @inbounds for i = 1:nA
        start = (i-1)*nB
        permAi = (A.perm[i]-1)*nB
        val = A.vals[i]
        @inbounds @simd for j = 1:nB
            perm[start+j] = permAi + j
            vals[start+j] = val
        end
    end
    PermMatrix(perm, vals)
end

function kron(A::IMatrix, B::PermMatrix{Tv, Ti}) where {Tv<:Number, Ti <: Integer}
    nA = size(A, 1)
    nB = size(B, 1)
    perm = Vector{Int}(undef, nB*nA)
    vals = Vector{Tv}(undef, nB*nA)
    @inbounds for i = 1:nA
        start = (i-1)*nB
        @inbounds @simd for j = 1:nB
            perm[start+j] = start +B.perm[j]
            vals[start+j] = B.vals[j]
        end
    end
    PermMatrix(perm, vals)
end

function kron(A::StridedMatrix{Tv}, B::PermMatrix{Tb}) where {Tv<:Number, Tb<:Number}
    mA, nA = size(A)
    nB = size(B, 1)
    perm = fast_invperm(B.perm)
    nzval = Vector{promote_type(Tv, Tb)}(undef, mA*nA*nB)
    rowval = Vector{Int}(undef, mA*nA*nB)
    colptr = collect(1:mA:nA*nB*mA+1)
    z = 1
    @inbounds for j = 1:nA
        @inbounds for j2 = 1:nB
            p2 = perm[j2]
            val2 = B.vals[p2]
            ir = p2
            @inbounds @simd for i = 1:mA
                nzval[z] = A[i, j]*val2  # merge
                rowval[z] = ir
                z += 1
                ir += nB
            end
        end
    end
    SparseMatrixCSC(mA*nB, nA*nB, colptr, rowval, nzval)
end

function kron(A::PermMatrix{Ta}, B::StridedMatrix{Tb}) where {Tb<:Number, Ta<:Number}
    mB, nB = size(B)
    nA = size(A, 1)
    perm = fast_invperm(A.perm)
    nzval = Vector{promote_type(Ta, Tb)}(undef, mB*nA*nB)
    rowval = Vector{Int}(undef, mB*nA*nB)
    colptr = collect(1:mB:nA*nB*mB+1)
    z = 0
    @inbounds for j = 1:nA
        colbase = (j-1)*nB
        p1 = perm[j]
        val2 = A.vals[p1]
        ir = (p1-1)*mB
        for j2 = 1:nB
            @inbounds @simd for i2 = 1:mB
                nzval[z+i2] = B[i2, j2]*val2
                rowval[z+i2] = ir+i2
            end
            z += mB
        end
    end
    SparseMatrixCSC(nA*mB, nA*nB, colptr, rowval, nzval)
end

function kron(A::PermMatrix{<:Number}, B::PermMatrix{<:Number})
    nA = size(A, 1)
    nB = size(B, 1)
    vals = kron(A.vals, B.vals)
    perm = Vector{Int}(undef, nB*nA)
    @inbounds for i = 1:nA
        start = (i-1)*nB
        permAi = (A.perm[i]-1)*nB
        @inbounds @simd for j = 1:nB
            perm[start+j] = permAi +B.perm[j]
        end
    end
    PermMatrix(perm, vals)
end

kron(A::PermMatrix{<:Number}, B::Diagonal{<:Number}) = kron(A, PermMatrix(B))
kron(A::Diagonal{<:Number}, B::PermMatrix{<:Number}) = kron(PermMatrix(A), B)

function kron(A::PermMatrix{Ta}, B::SparseMatrixCSC{Tb}) where {Ta<:Number, Tb<:Number}
    nA = size(A, 1)
    mB, nB = size(B)
    nV = nnz(B)
    perm = fast_invperm(A.perm)
    nzval = Vector{promote_type(Ta, Tb)}(undef, nA*nV)
    rowval = Vector{Int}(undef, nA*nV)
    colptr = Vector{Int}(undef, nA*nB+1)
    colptr[1] = 1
    @inbounds @simd for i in 1:nA
        start_row = (i-1)*nV
        start_ri = (perm[i]-1)*mB
        v0 = A.vals[perm[i]]
        @inbounds @simd for j = 1:nV
            nzval[start_row+j] = B.nzval[j]*v0
            rowval[start_row+j] = B.rowval[j] + start_ri
        end
        start_col = (i-1)*nB+1
        start_ci = (i-1)*nV
        @inbounds @simd for j = 1:nB
            colptr[start_col+j] = B.colptr[j+1] + start_ci
        end
    end
    SparseMatrixCSC(mB*nA, nB*nA, colptr, rowval, nzval)
end

function kron(A::SparseMatrixCSC{T}, B::PermMatrix{Tb}) where {T<:Number, Tb<:Number}
    nB = size(B, 1)
    mA, nA = size(A)
    nV = nnz(A)
    perm = fast_invperm(B.perm)
    rowval = Vector{Int}(undef, nB*nV)
    colptr = Vector{Int}(undef, nA*nB+1)
    nzval = Vector{promote_type(T, Tb)}(undef, nB*nV)
    z=1
    colptr[z] = 1
    @inbounds for i in 1:nA
        rstart = A.colptr[i]
        rend = A.colptr[i+1]-1
        @inbounds for k in 1:nB
            irow = perm[k]
            bval = B.vals[irow]
            irow_nB = irow - nB
            @inbounds @simd for r in rstart:rend
                rowval[z] = A.rowval[r]*nB+irow_nB
                nzval[z] = A.nzval[r]*bval
                z+=1
            end
            colptr[(i-1)*nB+k+1] = z
        end
    end
    SparseMatrixCSC(mA*nB, nA*nB, colptr, rowval, nzval)
end


text:
#   This script is not part of the PosDefManifoldML package.
#   It allows to build the PosDefManifoldML package
#   and its documentation locally from the source code,
#   without actually installing the package.
#   It is useful for developing purposes using the Julia
#   `Revise` package (that you need to have installed on your PC,
#   together with the `Documenter` package for building the documentation).
#   You won't need this script for using the package.
#
#   MIT License
#   Copyright (c) 2019, Marco Congedo, CNRS, Grenobe, France:
#   https://sites.google.com/site/marcocongedo/home
#
#   DIRECTIONS:
#   1) If you have installed PosDefManifoldML from github or Julia registry,
#      uninstall it.
#   2) Change the `juliaCodeDir` path here below to the path
#           where the PosDefManifoldML folder is located on your computer.
#   3) Under Linux, replace all '\\' with `/`
#   4) Put the cursor in this unit and hit SHIFT+CTRL+ENTER
#
#   Nota Bene: all you need for building the package is actually
#   the 'push!' line and the 'using' line.
#   You can safely delete the rest once
#   you have identified the 'srcDir' to be used in the push command.

begin
  juliaCodeDir = homedir()*"\\Documents\\Code\\julia\\"
  scrDir       = juliaCodeDir*"PosDefManifoldML\\src\\"
  docsDir      = juliaCodeDir*"PosDefManifoldML\\docs\\"

  push!(LOAD_PATH, scrDir)
  using LinearAlgebra,
        #Base.Threads, Random, Dates, Distributions, GLMNet, LIBSVM,
        Documenter, Statistics, PosDefManifold, Revise, PosDefManifoldML

  cd(docsDir)
  clipboard("""makedocs(sitename="PosDefManifoldML", modules=[PosDefManifoldML])""")
  @info("\nhit CTRL+V+ENTER on the REPL for building the documentation.");
end


text:
using Test

using MathOptInterface
const MOI = MathOptInterface
const MOIT = MathOptInterface.Test
const MOIU = MathOptInterface.Utilities
const MOIB = MathOptInterface.Bridges

include("../utilities.jl")

mock = MOIU.MockOptimizer(MOIU.Model{Float64}())
config = MOIT.TestConfig()

@testset "RSOC" begin
    bridged_mock = MOIB.Constraint.RSOC{Float64}(mock)

    MOIT.basic_constraint_tests(
        bridged_mock, config,
        include = [(F, S)
                   for F in [MOI.VectorOfVariables, MOI.VectorAffineFunction{Float64},
                             MOI.VectorQuadraticFunction{Float64}]
                   for S in [MOI.RotatedSecondOrderCone]])

    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [0.5, 1.0, 1/√2, 1/√2],
                          (MOI.SingleVariable,                MOI.EqualTo{Float64}) => [-√2, -1/√2],
                          (MOI.VectorAffineFunction{Float64}, MOI.SecondOrderCone)  => [[3/2, 1/2, -1.0, -1.0]])
    MOIT.rotatedsoc1vtest(bridged_mock, config)
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1/√2, 1/√2],
                          (MOI.VectorAffineFunction{Float64}, MOI.SecondOrderCone)  => [[3/2, 1/2, -1.0, -1.0]])
    MOIT.rotatedsoc1ftest(bridged_mock, config)
    ci = first(MOI.get(bridged_mock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.RotatedSecondOrderCone}()))
    test_delete_bridge(bridged_mock, ci, 2, ((MOI.VectorAffineFunction{Float64}, MOI.SecondOrderCone, 0),))
end


text:
using DataVoyager
using DataValues
using Electron
using Test

@testset "DataVoyager" begin

source = [(a=1,b=1), (a=2,b=2)]
source2 = [(a=DataValue(1), b=DataValue{Int}()), (a=DataValue{Int}(), b=DataValue(2))]

v = Voyager()
@test typeof(v.w) == Electron.Window

v = Voyager(source)
@test typeof(v.w) == Electron.Window

v(source)
@test typeof(v.w) == Electron.Window

source |> v
@test typeof(v.w) == Electron.Window

v = Voyager(source2)
@test typeof(v.w) == Electron.Window

v(source2)
@test typeof(v.w) == Electron.Window

source2 |> v
@test typeof(v.w) == Electron.Window

end


text:
# Sugarscape
# Chapter II
module ChapterII

using StatsBase: sample
import PyPlot
import Colors

## Types ##

# We define an abstract type AgentInfo with two subtypes: NoAgent and Agent.
# This allows for removing an agent from the board while keeping his entry in
# the agent list. Also it allows for recursive type definition, an old
# unresolved Julia issue #269.

abstract AgentInfo
type NoAgent <: AgentInfo
end

type Place
    x::Int
    y::Int
    agent::AgentInfo
    sugar::Int
    capacity::Int
end

type Agent <: AgentInfo
    place::Place
    vision::Int
    metabolism::Int
    stash::Int
end

# We only keep in this version the current lattice and agents in memory, not the
# history.
type Scape
    lattice::Array{Place}
    agents::Vector{AgentInfo}
    steps::Int
end

# convenience functions
Base.show(io::IO, scape::Scape) = println(io, "Sugarscape with $(length(agents(scape))) agents")
agents(scape::Scape) = scape.agents
Base.isempty(place::Place) = isa(place.agent, NoAgent)
alive(agent::AgentInfo) = isa(agent, Agent)
living(scape) = scape.agents[map(alive, scape.agents)]

## Initialization ##

# We define a sugarscape capacity similar to graphs in the book.
function init_capacity()
    grid = 50
    hump1 = [15, 15]
    hump2 = [40, 40]
    dist(xy, hump) = hypot(xy[1]-hump[1], xy[2]-hump[2])
    capacity = Int[max(0, 4 - floor(min(dist([i,j], hump1), dist([i,j], hump2))/5)) for i=1:grid,j=1:grid]
end

# Keyword arguments still have a performance penalty if not type defined.
function init_scape(capacity; N_agents::Int=400, init_stash=10,
                    init_vision=()->rand(1:6), init_metabolism=()->rand(1:4))
    gridx, gridy = size(capacity)
    # create empty lattice with sugar levels at full capacity
    lattice = [Place(i,j, NoAgent(), capacity[i,j], capacity[i,j]) for i=1:gridx, j=1:gridy]

    #populate with agents
    agents = AgentInfo[]
    for place in sample(lattice, N_agents, replace=false)
        vision = init_vision()
        metabolism = init_metabolism()
        stash = init_stash
        agent = Agent(place, vision, metabolism, stash)
        push!(agents, agent)

        place.agent = agent
    end

    Scape(lattice, agents, 0)
end

## Movement ##

# (function input annotations are not necessary (not even for speed), just for type checks)

#step functions to select neighbouring place with circular boundaries
 left(p::Place, lat::Array{Place}) = lat[p.x==1 ? end : p.x-1, p.y]
right(p::Place, lat::Array{Place}) = lat[p.x==end ? 1 : p.x+1, p.y]
   up(p::Place, lat::Array{Place}) = lat[p.x, p.y==1 ? end : p.y-1]
 down(p::Place, lat::Array{Place}) = lat[p.x, p.y==end ? 1 : p.y+1]

# function to move an agent to a destination place
function move(agent::Agent, dest_place::Place)
    agent.place === dest_place && return
    isempty(dest_place) || error("destination place occupied")

    agent.place.agent = NoAgent()
    dest_place.agent = agent
    agent.place = dest_place
    return nothing
end

# function to select best next step in view of agent.
function view_place(agent::Agent, step_function::Function, lattice::Array{Place})
    p = agent.place

    p_out = pstep = p #if no unoccupied place in sight, stay put
    sugar = -1 #Agent will move to empty place rather than stay put
    for _ in 1:agent.vision
        pstep = step_function(pstep, lattice)
        # Only move further if more sugar, not on equal sugar level.
        # As in book, agents will take steps of size=1 when far from sugar
        # mountains, but this could be improved to larger steps.
        # It seems from Animation II-1 that agents do not jump over other agents
        # (as a horse), but only slide like a rook.

        if isempty(pstep)
            if pstep.sugar > sugar
                sugar = pstep.sugar
                p_out = pstep
            end
        else #disallow jumps
            break
        end
    end
    p_out
end

# to sort array of places by sugar, define `isless`
Base.isless(p1::Place, p2::Place) = p1.sugar < p2.sugar

function harvest(agent)
    agent.place.sugar == 0 && return(nothing)
    agent.stash += agent.place.sugar
    agent.place.sugar = 0
    nothing
end

# move all agents and harvest
function move(scape::Scape)
    lat = scape.lattice
    direction = [left, right, up, down]

    @inbounds for agent in shuffle(living(scape))
        view_places = Array(Place, 4)
        # Randomize search directions in case of equal sugar place, as in book.
        shuffle(direction)
        for i = 1:4
            view_places[i] = view_place(agent, direction[i], lat)
        end
        select_place = last(sort!(view_places))
        move(agent, select_place)
        harvest(agent)
    end
end

## Evolve ##

function grow(scape::Scape; α::Int=1)
    @inbounds for place in scape.lattice
        place.sugar = min(place.capacity, place.sugar + α)
    end
end

function consume(scape::Scape)
    @inbounds for agent in living(scape)
        agent.stash -= agent.metabolism
    end
end

function Base.kill(scape::Scape)
    ag = scape.agents #living(scape) will make copy

    @inbounds for i = 1:length(ag)
        if alive(ag[i]) && (ag[i].stash < 0)
            ag[i].place.agent = NoAgent()
            ag[i] = NoAgent()
        end
    end
end

function timestep(scape::Scape; α::Int=1, args...)
    grow(scape, α=α)
    move(scape) #includes harvestig
    consume(scape)
    kill(scape)
    scape.steps += 1
    nothing
end

function run(N; kwargs...)
    scape = init_scape(init_capacity(); kwargs...)
    for cnt = 1:N
        timestep(scape; kwargs...)
    end
    scape
end

## Plotting ##

function plot(scape::Scape)
    lattice = scape.lattice
    gridx, gridy = size(lattice)
    sugar = reshape([p.sugar for p in lattice], gridx, gridy)
    s_i, s_j = collect(ind2sub(size(lattice), 1:length(lattice)))
    colmap = PyPlot.ColorMap([Colors.RGB(1,1,1), Colors.RGB(1,1,0)],5,1.5)
    PyPlot.scatter(s_i, s_j, s=20*sugar, c=sugar, cmap=colmap, marker="o", lw=0)

    step_agents = living(scape)
    a_i = [a.place.x for a in step_agents]
    a_j = [a.place.y for a in step_agents]
    PyPlot.scatter(a_i, a_j, c="red", lw=0)

    PyPlot.xlim(0, gridx+1); PyPlot.ylim(0, gridy+1)
    PyPlot.title("Sugarscape step $(scape.steps)")
end

end #module


text:
## --- Load (and install if neccesary) the StatGeochem package which has the resampling functions we'll want

    using StatGeochem
    using Plots; gr();

    if VERSION>=v"0.7"
        using Statistics, DelimitedFiles, SpecialFunctions
    end

## --- Download and unzip Keller and Schoene (2012) dataset

    if ~isfile("ign.h5") # Unless it already exists
        download("https://storage.googleapis.com/statgeochem/ign.h5.gz","./ign.h5.gz")
        download("https://storage.googleapis.com/statgeochem/err2srel.csv","./err2srel.csv")
        run(`gunzip -f ign.h5.gz`) # Unzip file
    end

    # Read HDF5 file
    using HDF5
    ign = h5read("ign.h5","vars")

## --- Compute proximity coefficients (inverse weights)

    # # Compute inverse weights
    # k = invweight(ign["Latitude"] .|> Float32, ign["Longitude"] .|> Float32, ign["Age"] .|> Float32)

    # Since this is pretty computatually intensive, let's load a precomputed version instead
    k = ign["k"]

    # Probability of keeping a given data point when sampling
    p = 1.0 ./ ((k .* median(5.0 ./ k)) .+ 1.0) # Keep rougly one-fith of the data in each resampling

    # Set absolute uncertainties for each element where possible, using errors defined inerr2srel.csv
    err2srel = importdataset("err2srel.csv", ',')
    for e in ign["elements"]
        # If there's an err2srel for this variable, create a "_sigma" if possible
        if haskey(err2srel, e) && !haskey(ign, e*"_sigma")
            ign[e*"_sigma"] = ign[e] .* (err2srel[e] / 2);
        end
    end

    # Special cases: age uncertainty
    ign["Age_sigma"] = (ign["Age_Max"]-ign["Age_Min"])/2;
    t = (ign["Age_sigma"] .< 50) .| isnan.(ign["Age_sigma"]) # Find points with < 50 Ma absolute uncertainty
    ign["Age_sigma"][t] .= 50 # Set 50 Ma minimum age uncertainty (1-sigma)

    # Special cases: location uncertainty
    ign["Latitude_sigma"] = ign["Loc_Prec"]
    ign["Longitude_sigma"] = ign["Loc_Prec"]

## --- Resample a single variable

    xmin = 0 # Minimum Age
    xmax = 3900 # Maximum Age
    nbins = 39
    elem = "K2O" # Element to plot

    # Look only at samples from a specific silica range
    t = 43 .< ign["SiO2"] .< 51 # Mafic
    # t = 51 .< ign["SiO2"] .< 62  # Intermediate
    # t = 62 .< ign["SiO2"] .< 74 # Felsic
    # t = 40 .< ign["SiO2"] .< 80 # All normal igneous
    # t = trues(size(ign[elem])) # Everything

    # Resample, returning binned means and uncertainties
    # (c = bincenters, m = mean, el = lower 95% CI, eu = upper 95% CI)
    (c,m,el,eu) = bin_bsr_means(ign["Age"][t],ign[elem][t],xmin,xmax,nbins, p=p[t], x_sigma=ign["Age_sigma"][t])

    # Plot results
    plot(c,m,yerror=(el,eu),seriestype=:scatter,color=:darkblue,markerstrokecolor=:darkblue,label="")
    plot!(xlabel="Age (Ma)", ylabel="$elem (wt. %)",xlims=(xmin,xmax),framestyle=:box,grid=:off,xflip=true) # Format plot

## --- Multiple silica ranges together

    xmin = 0 # Minimum Age
    xmax = 800 # Maximum Age
    nbins = 40
    elem = "Al2O3" # Element to plot
    rsi = [43,51,62,74,80] # Ranges of silica to plot together

    t = trues(size(ign[elem]))
    h = plot(xlabel="Age (Ma)", ylabel="$elem (wt. %)",xlims=(xmin,xmax),framestyle=:box,grid=:off,xflip=true) # Format plot
    for i=1:length(rsi)-1
        t .= rsi[i] .< ign["SiO2"] .< rsi[i+1]

        # Resample, returning binned means and uncertainties
        # (c = bincenters, m = mean, el = lower 95% CI, eu = upper 95% CI)
        (c,m,el,eu) = bin_bsr_means(ign["Age"][t],ign[elem][t],xmin,xmax,nbins, p=p[t], x_sigma=ign["Age_sigma"][t])

        # Plot results
        plot!(h, c,m,yerror=(el,eu),seriestype=:scatter,color=lines[i],markerstrokecolor=lines[i],label="$(rsi[i])-$(rsi[i+1]) % SiO2")
    end
    savefig("$(elem)_$(xmin)-$(xmax) Ma.pdf")
    display(h)

## ---  Resample a ratio

    tmin = 0 # Minimum age
    tmax = 3900 # Maximum age
    nbins = 39
    num = "Sm" # Numerator
    denom = "Nd" # Denominator

    # Look only at samples from a specific silica range
    t = 43 .< ign["SiO2"] .< 51 # Mafic
    # t = 51 .< ign["SiO2"] .< 62 # Intermediate
    # t = 62 .< ign["SiO2"] .< 74 # Felsic

    # Exclude outliers
    t = t .& inpctile(ign[num], 99) .& inpctile(ign[denom], 99)

    # Resample, returning binned means and uncertainties
    # (c = bincenters, m = mean, el = lower 95% CI, eu = upper 95% CI)
    (c,m,el,eu) = bin_bsr_ratios(ign["Age"][t],ign[num][t],ign[denom][t],tmin,tmax,nbins, p=p[t],
                    x_sigma=ign["Age_sigma"][t])

    # Plot results
    h = plot(c,m,yerror=(el,eu),seriestype=:scatter,color=:darkred,markerstrokecolor=:darkred,label="")
    plot!(h, xlabel="Age (Ma)", ylabel="$(num) / $(denom)",xlims=(tmin,tmax),framestyle=:box,grid=:off,xflip=true) # Format plot
    display(h)

    # savefig(h,"$(num)$(denom)_$(tmax)-$(tmin)Ma.pdf")

## --- Single element differentiation example

    xelem = "SiO2"
    xmin = 45
    xmax = 75
    nbins = 8
    elem = "K2O"

    h = plot(xlabel=xelem, ylabel="$(elem)",xlims=(xmin,xmax),framestyle=:box,grid=:off,fg_color_legend=:white) # Format plot

    rt = [0,1,2,3,4] # Time range (Ga)
    colors = reverse(resize_colormap(viridis[1:end-20],length(rt)-1))
    for i=1:length(rt)-1
        t = rt[i]*1000 .< ign["Age"] .< rt[i+1]*1000

        # Resample, returning binned means and uncertainties
        # (c = bincenters, m = mean, el = lower 95% CI, eu = upper 95% CI)
        (c,m,el,eu) = bin_bsr_means(ign[xelem][t],ign[elem][t],xmin,xmax,nbins, p=p[t],
                        x_sigma=ign[xelem*"_sigma"][t], y_sigma=ign[elem*"_sigma"][t])

        # Plot results
        plot!(h, c,m,yerror=(el,eu),color=colors[i],mscolor=colors[i],seriestype=:scatter,label="$(rt[i])-$(rt[i+1]) Ga")
        plot!(h, c,m,style=:dot,color=colors[i],mscolor=colors[i],label="")
    end
    # savefig(h,"$(xelem)_$(num)$(denom).pdf")
    display(h)

## --- Ratio differentiation example

    xelem = "SiO2"
    xmin = 45
    xmax = 75
    nbins = 8
    num = "Sm" # Numerator
    denom = "Nd" # Denominator

    h = plot(xlabel=xelem, ylabel="$(num) / $(denom)",xlims=(xmin,xmax),framestyle=:box,grid=:off,legend=:topleft,fg_color_legend=:white) # Format plot

    rt = [0,1,2,3,4]
    colors = reverse(resize_colormap(viridis[1:end-20],length(rt)-1))
    for i=1:length(rt)-1
        t = rt[i]*1000 .< ign["Age"] .< rt[i+1]*1000

        # Resample, returning binned means and uncertainties
        # (c = bincenters, m = mean, el = lower 95% CI, eu = upper 95% CI)
        (c,m,el,eu) = bin_bsr_ratios(ign[xelem][t],ign[num][t],ign[denom][t],xmin,xmax,nbins, p=p[t],
                        x_sigma=ign[xelem*"_sigma"][t], num_sigma=ign[num*"_sigma"][t], denom_sigma=ign[denom*"_sigma"][t])
        # Plot results
        plot!(h, c,m,yerror=(el,eu),color=colors[i], mscolor=colors[i], seriestype=:scatter,label="$(rt[i])-$(rt[i+1]) Ga")
        plot!(h, c,m,style=:dot,color=colors[i],label="")
    end
    display(h)

## --- Ratio differentiation

    xelem = "SiO2"
    xmin = 40 # Minimum age
    xmax = 80 # Maximum age
    nbins = 20
    num = "Sc" # Numerator
    denom = "Yb" # Denominator

    # Exclude outliers
    t = inpctile(ign[num], 99) .& inpctile(ign[denom], 99)

    # Resample, returning binned means and uncertainties
    # (c = bincenters, m = mean, el = lower 95% CI, eu = upper 95% CI)
    (c,m,el,eu) = bin_bsr_ratios(ign[xelem][t],ign[num][t],ign[denom][t],xmin,xmax,nbins, p=p[t],
                    x_sigma=ign[xelem][t]*0.01, num_sigma=ign[num][t]*0.05, denom_sigma=ign[denom][t]*0.05)

    # Plot results
    h = plot(c,m,yerror=(el,eu),seriestype=:scatter,color=:darkblue,markerstrokecolor=:darkblue,label="")
    plot!(h, xlabel=xelem, ylabel="$(num) / $(denom)",xlims=(xmin,xmax),framestyle=:box,grid=:off) # Format plot
    display(h)

## --- Export differentiation trends

    xelem = "SiO2"
    xmin = 45
    xmax = 75
    nbins = 10

    rt = [0,1,2,3,4] # Time range (Ga)
    data = Dict()
    for elem in ["Al2O3", "MgO", "Na2O", "Fe2O3T", "K2O", "CaO"]
        for i=1:length(rt)-1
            t = rt[i]*1000 .< ign["Age"] .< rt[i+1]*1000

            # Resample, returning binned means and uncertainties
            (c,m,e) = bin_bsr(ign[xelem][t],ign[elem][t],xmin,xmax,nbins, p=p[t],
                            x_sigma=ign[xelem*"_sigma"][t], y_sigma=ign[elem*"_sigma"][t])
            data[xelem] = c
            data[elem*"_$(rt[i])-$(rt[i+1])Ga"] = m
            data[elem*"_$(rt[i])-$(rt[i+1])Ga_sigma"] = e
        end
    end
    exportdataset(data,"MajorDifferentiation.csv",',')

## --- End of File
""


text:
function NACA(number::AbstractString, x::AbstractVector; sharpTE::Bool=false)
    ndig = length(strip(number))
    ndig == 4 || error("invalid number=\"$number\": must be a 4-digit integer")
    xc = convert(Vector{Float64}, x)
    sort!(xc)
    if xc[1] != 0
        insert!(xc, 1, 0)
    end
    return NACA(Val(ndig), parse(Int, number), xc, sharpTE)
end

function NACA(number::AbstractString, np::Int=33; sharpTE::Bool=false)
    xc = collect(LinRange(0, π, np))
    @. xc = (1-cos(xc))/2
    return NACA(number, xc; sharpTE=sharpTE)
end

function NACA_T(xc::AbstractVector{Float64}, t::Float64, sharpTE::Bool)
    a0 = +0.2969
    a1 = -0.1260
    a2 = -0.3516
    a3 = +0.2843
    a4 = sharpTE ? -0.1036 : -0.1015
    yt = similar(xc)
    @. yt = 5t*(a0*sqrt(xc)+a1*xc+a2*xc^2+a3*xc^3+a4*xc^4)
    return yt
end

function NACA_Airfoil(
    name::AbstractString,
    desc::AbstractString,
    xc::AbstractVector{Float64},
    yt::AbstractVector{Float64},
    yc::AbstractVector{Float64},
    dyc::AbstractVector{Float64})

    np = length(xc)
    x = zeros(2np-1)
    y = zeros(2np-1)
    x[np] = y[np] = 0 # leading edge
    @inbounds for i in 2:np
        iu = np-i+1
        il = np+i-1
        θ = atan(dyc[i])
        x[iu] = xc[i] - yt[i]*sin(θ)
        x[il] = xc[i] + yt[i]*sin(θ)
        y[iu] = yc[i] + yt[i]*cos(θ)
        y[il] = yc[i] - yt[i]*cos(θ)
    end
    return AirfoilData(name, desc, x, y)
end

function NACA(::Val{4}, num::Int, xc::AbstractVector{Float64}, sharpTE::Bool)
    @assert num > 0 "number must be positive"
    m = (num ÷ 1000) * 0.01
    p = ((num ÷ 100) % 10) * 0.1
    t = (num % 100) * 0.01
    yt = NACA_T(xc, t, sharpTE)
    np = length(xc)
    yc = zeros(np)
    dyc = zeros(np)
    name = @sprintf "NACA %04d" num
    desc = @sprintf "Thickness=%d%%, Camber=%d%% AT %d%%" t*100 m*100 p*100
    if p > 0
        @inbounds for i in eachindex(xc)
            xx = xc[i]
            if xx <= p
                yc[i] = m * p^-2 * (2p*xx-xx^2)
                dyc[i] = 2m * p^-2 * (p-xx)
            else
                yc[i] = m * (1-p)^-2 * ((1-2p) + 2p*xx - xx^2)
                dyc[i] = 2m * (1-p)^-2 * (p-xx)
            end
        end
    end
    return NACA_Airfoil(name, desc, xc, yt, yc, dyc)
end


text:
using Makie, CatmullClark

const inputpoints = [
    [-1.0, -1.0, -1.0],
    [-1.0, -1.0, 1.0],
    [-1.0, 1.0, -1.0],
    [-1.0, 1.0, 1.0],
    [1.0, -1.0, -1.0],
    [1.0, -1.0, 1.0],
    [1.0, 1.0, -1.0],
    [1.0, 1.0, 1.0]]

const inputfaces = [
    [0, 4, 5, 1],
    [4, 6, 7, 5],
    [6, 2, 3, 7],
    [2, 0, 1, 3],
    [1, 5, 7, 3],
    [0, 2, 6, 4]]

const donutpoints = [
    [-2.0, -0.5, -2.0], [-2.0, -0.5, 2.0], [2.0, -0.5, -2.0], [2.0, -0.5, 2.0],
    [-1.0, -0.5, -1.0], [-1.0, -0.5, 1.0], [1.0, -0.5, -1.0], [1.0, -0.5, 1.0],
    [-2.0, 0.5, -2.0], [-2.0, 0.5, 2.0], [2.0,  0.5, -2.0], [2.0, 0.5, 2.0],
    [-1.0, 0.5, -1.0], [-1.0, 0.5, 1.0], [1.0, 0.5, -1.0], [1.0, 0.5, 1.0]]

const donutfaceindices = [
    [4, 5, 1, 0], [3, 1, 5, 7], [0, 2, 6, 4], [2, 3, 7, 6],
    [8,  9, 13, 12], [15, 13, 9, 11], [12, 14, 10, 8], [14, 15, 11, 10],
    [0, 1, 9, 8], [1, 3, 11, 9], [2, 0, 8, 10], [3, 2, 10, 11],
    [12, 13, 5, 4], [13, 15, 7, 5], [14, 12, 4, 6], [15, 14, 6, 7]]

const tetrapoints = [[1.0, 1.0, 1.0], [1.0, -1.0, -1.0], [-1, 1, -1], [-1, -1, 1]]

const tetrafaceindices = [[1, 2, 3], [1, 3, 4], [2, 3, 4], [1, 2, 4]]

const faces = [map(x -> Point3f0(inputpoints[x]), p .+ 1) for p in inputfaces]
const donutfaces = [map(x -> Point3f0(donutpoints[x]), p .+ 1) for p in donutfaceindices]
const tetrafaces = [map(x -> Point3f0(tetrapoints[x]), p) for p in tetrafaceindices]

# cube, rounds toward a sphere
scene = drawfaces(faces, :black)
display(scene)
setscene(scene)
sleep(1)
catmullclark(faces, 4, CatmullClark.displaycallback)

# tetrahedron
scene2 = drawfaces(tetrafaces, :black)
display(scene2)
setscene(scene2)
sleep(2)
catmullclark(tetrafaces, 3, CatmullClark.displaycallback)

# torus
scene3 = drawfaces(donutfaces, :black)
display(scene3)
setscene(scene3)
sleep(1)
catmullclark(donutfaces, 3, CatmullClark.displaycallback)

# if a face missing in cube, makes cuplike shape
sleep(3)
scene4 = drawfaces(faces[2:end], :black)
display(scene4)
setscene(scene4)
sleep(2)
catmullclark(faces[2:end], 3, CatmullClark.displaycallback)

println("Press Enter to continue", readline())


text:
import FileIO: File, filename, file_extension, @format_str, query

function iscompiled(shader::GLuint)
    success = GLint[0]
    glGetShaderiv(shader, GL_COMPILE_STATUS, success)
    return first(success) == GL_TRUE
end

abstract type AbstractShader end

struct Shader <: AbstractShader
    name    ::Symbol
    source  ::Vector{UInt8} #UInt representation of the source program string,
    typ     ::GLenum
    id      ::GLuint
    function Shader(name, source, typ, id)
        new(name, source, typ, id)
    end
end

function Shader(name, shadertype, source::Vector{UInt8})
    shaderid = glCreateShader(shadertype)::GLuint
    @assert shaderid > 0 "opengl context is not active or shader type not accepted. Shadertype: $(GLENUM(shadertype).name)"
    glShaderSource(shaderid, source)
    glCompileShader(shaderid)
    if !iscompiled(shaderid)
        print_with_lines(String(source))
        @error "shader $(name) didn't compile. \n$(getinfolog(shaderid))"
    end
    Shader(name, source, shadertype, shaderid)
end
function Shader(path::String, source_str::AbstractString)
    typ = shadertype(query(path))
    source = Vector{UInt8}(source_str)
    name = Symbol(path)
    Shader(name, typ, source)
end
Shader(path::File{format"GLSLShader"}) = load(path)

import Base: ==
(==)(a::Shader, b::Shader) = a.source == b.source && a.typ == b.typ && a.id == b.id && a.context == b.context
Base.hash(s::Shader, h::UInt64) = hash((s.source, s.typ, s.id, s.context), h)

function Base.show(io::IO, shader::Shader)
    println(io, GLENUM(shader.typ).name, " shader: $(shader.name))")
    println(io, "source:")
    print_with_lines(io, String(shader.source))
end

shadertype(s::Shader) = s.typ
function shadertype(f::File{format"GLSLShader"})
    shadertype(file_extension(f))
end
function shadertype(ext::AbstractString)
    ext == ".comp" && return GL_COMPUTE_SHADER
    ext == ".vert" && return GL_VERTEX_SHADER
    ext == ".frag" && return GL_FRAGMENT_SHADER
    ext == ".geom" && return GL_GEOMETRY_SHADER
    @error "$ext not a valid shader extension."
end
function shadertype(typ::Symbol)
    (typ == :compute  || typ == :comp) && return GL_COMPUTE_SHADER
    (typ == :vertex   || typ == :vert) && return GL_VERTEX_SHADER
    (typ == :fragment || typ == :frag) && return GL_FRAGMENT_SHADER
    (typ == :geometry || typ == :geom) && return GL_GEOMETRY_SHADER
    @error "$typ not a valid shader symbol."
end

#Implement File IO interface
function load(f::File{format"GLSLShader"})
    fname = filename(f)
    source = read(open(fname), String)
    Shader(fname, source)
end

function save(f::File{format"GLSLShader"}, data::Shader)
    s = open(f, "w")
    write(s, data.source)
    close(s)
end

# Different shader string literals- usage: e.g. frag" my shader code"
macro frag_str(source::AbstractString)
    quote
        ($source, GL_FRAGMENT_SHADER)
    end
end
macro vert_str(source::AbstractString)
    quote
        ($source, GL_VERTEX_SHADER)
    end
end
macro geom_str(source::AbstractString)
    quote
        ($source, GL_GEOMETRY_SHADER)
    end
end
macro comp_str(source::AbstractString)
    quote
        ($source, GL_COMPUTE_SHADER)
    end
end

function getinfolog(id::GLuint)
    # Get the maximum possible length for the descriptive error message
    maxlength = GLint[0]
    glGetShaderiv(id, GL_INFO_LOG_LENGTH, maxlength)
    maxlength = first(maxlength)
    # Return the text of the message if there is any
    if maxlength > 0
        buffer = zeros(GLchar, maxlength)
        sizei = GLsizei[0]
        glGetShaderInfoLog(id, maxlength, sizei, buffer)
        length = first(sizei)
        return unsafe_string(pointer(buffer), length)
    else
        return "success"
    end
end


text:
using MakieLayout
using Makie

begin
    scene = Scene(resolution = (1000, 1000));
    screen = display(scene)
    campixel!(scene);

    nrows = 1
    ncols = 1

    maingl = GridLayout(
        scene, nrows, ncols,
        colsizes = Auto(false, 1),
        rowsizes = Auto(false, 1),
        alignmode = Outside(30))

    las = [maingl[i, j] = LAxis(scene) for i in 1:nrows, j in 1:ncols]

    linkxaxes!(las...)
    linkyaxes!(las...)

    for i in 1:nrows, j in 1:ncols

        scatter!(las[i, j], rand(200, 2) .+ [i j], markersize=0.1, color=(:black, 0.3))

        i > 1 && (las[i, j].titlevisible = false)
        j > 1 && (las[i, j].ylabelvisible = false)
        j > 1 && (las[i, j].yticklabelsvisible = false)
        j > 1 && (las[i, j].yticksvisible = false)
        i < nrows && (las[i, j].xticklabelsvisible = false)
        i < nrows && (las[i, j].xticksvisible = false)
        i < nrows && (las[i, j].xlabelvisible = false)
    end

    tl = maingl[0, :] = LText(scene, text="Super Title", padding=(10, 10, 20, 20), textsize=50)
    maingl[1, :] = LRect(scene, strokevisible=false)
    stl = maingl[2:end, end+1] = LText(scene, text="Side Title", textsize=50, rotation=-pi/2)

    slgl = maingl[end+1, 1:end-1] = GridLayout(1, 2)

    slgl[1, 1] = LText(scene, text="Supertitle Size", halign=:left)
    sl1 = slgl[1, 2] = LSlider(scene, height = 30, range = 1:200)
    on(sl1.value) do val
        @async tl.attributes.textsize = val
    end

    slgl[2, 1] = LText(scene, text="Sidetitle Size", halign=:left)
    sl2 = slgl[2, 2] = LSlider(scene, height = 30, range = 1:200)
    on(sl2.value) do val
        @async stl.attributes.textsize = val
    end

    nothing
end

# save("layout.png", scene)


text:
"""
This struct contains all data for running a CEM for price forecasts or for endoogeneous Economic Dispatch.
    zones: Name of zones included in the market clearing problem.
    lines: Vector of inter-zonal transmission lines included in the market clearing problem.
    capital_cost_multiplier: Average capital cost multiplier for new builds.
    inv_periods: Vector of all markets modeled for the entire CEM planning horizon.
    projects: Vector of all projects included in the market clearing problem.
    rep_hour_weight: Weight of representative hours in the market clearing problem.
"""

struct MarketClearingProblem{Z, T}

    zones::Vector{String}
    lines::Vector{ZonalLine}
    capital_cost_multiplier::Float64
    inv_periods::Vector{MarketCollection{Z, T}}
    projects::Vector{MarketProject}
    rep_hour_weight::Vector{Float64}

    function MarketClearingProblem(
        zones::Vector{String},
        lines::Vector{ZonalLine},
        capital_cost_multiplier::Float64,
        inv_periods::Vector{MarketCollection{Z, T}},
        projects::Vector{MarketProject},
        rep_hour_weight::Vector{Float64}
    ) where {Z, T}
        @assert T == length(rep_hour_weight)
        @assert length(zones) ==  Z
        @assert allunique(getproperty.(projects, :name))
        new{Z, T}(zones, lines, capital_cost_multiplier, inv_periods, projects, rep_hour_weight)
    end

end


text:
#!/usr/bin/env julia

# Wrapper for GDIP library libGDIP.so

using Libdl

# Open the library explicitly.
lib = Libdl.dlopen("$(@__DIR__)/../gdip/lib/libGDIP")

@inline function gdip_get_length()
    return ccall( (:julia_get_length, :libGDIP), Cdouble, () )
end

@inline function gdip_sample_many(step::Float64)
    length = gdip_get_length()
    samples = Vector{Vector{Float64}}()
    sample_l = collect(0:step:length)
    push!(sample_l, length)
    for l in sample_l
        ccall( (:julia_sample_state_to_tmp, :libGDIP), Cvoid, (Cdouble,), l)
        x = ccall( (:julia_get_tmp_x, :libGDIP), Cdouble, (), )
        y = ccall( (:julia_get_tmp_y, :libGDIP), Cdouble, (), )
        theta = ccall( (:julia_get_tmp_theta, :libGDIP), Cdouble, (), )
        push!(samples, [x,y,theta])
    end
    return samples
end

@inline function gdip_init_dubins_maneuver(
    from::Array{Float64,1}, to::Array{Float64,1}, turningRadius::Float64
)
    length = ccall(
        (:julia_set_configurations, :libGDIP),                            # name of C function and library
        Cdouble,                                                          # output type
        (Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble),  # tuple of input types
        from[1], from[2], from[3], to[1], to[2], to[3], turningRadius     # names of Julia variables to pass in
    )
    return length
end

@inline function gdip_init_dip(
    p1::Array{Float64,1}, interval1::Array{Float64,1},
    p2::Array{Float64,1}, interval2::Array{Float64,1},
    turningRadius::Float64
)
    length = ccall(
        (:julia_set_configurations_dip, :libGDIP),                                                         # name of C function and library
        Cdouble,                                                                                # output type
        (Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble),                 # tuple of input types
        p1[1], p1[2], interval1[1], interval1[2], p2[1], p2[2], interval2[1], interval2[2], turningRadius  # names of Julia variables to pass in
    )
    return length
end

@inline function gdip_init_gdip(
    p1::Array{Float64,1}, interval1::Array{Float64,1}, radius1::Float64,
    p2::Array{Float64,1}, interval2::Array{Float64,1}, radius2::Float64,
    turningRadius::Float64
)
    length = ccall(
        (:julia_set_configurations_gdip, :libGDIP),                                                                          # name of C function and library
        Cdouble,                                                                                # output type
        (Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble),                 # tuple of input types
        p1[1], p1[2], interval1[1], interval1[2], radius1, p2[1], p2[2], interval2[1], interval2[2], radius2,  turningRadius # names of Julia variables to pass in
    )
    return length
end

@inline function gdip_dubins_circle_intersection(
    circle_origin::Array{Float64,1}, circle_radius::Float64
)
    isintersecting = ccall(
        (:julia_dubins_circle_intersection, :libGDIP),                                                                          # name of C function and library
        Cuchar,                                                                                # output type
        (Cdouble, Cdouble, Cdouble),                        # circle_origin_x, circle_origin_y, circle_radius
        circle_origin[1], circle_origin[2], circle_radius   # names of Julia variables to pass in
    )
    return isintersecting > 0
end

@inline function gdip_dubins_closest(point::Array{Float64,1})
    ccall(
        (:julia_dubins_closest, :libGDIP),                                                                          # name of C function and library
        Cvoid,                                                                                # output type
        (Cdouble, Cdouble),  # point_x, point_y
        point[1], point[2]   # names of Julia variables to pass in
    )
    x = ccall( (:julia_get_tmp_x, :libGDIP), Cdouble, (), )
    y = ccall( (:julia_get_tmp_y, :libGDIP), Cdouble, (), )
    theta = ccall( (:julia_get_tmp_theta, :libGDIP), Cdouble, (), )
    return [x,y,theta]
end

max_stars_repo_name: petvana/gdip

max_issues_repo_name: petvana/gdip

max_forks_repo_name: petvana/gdip


text:
bike_monday_env = {
    "coverage"  => [3, 4, 5 * ones(Int, 6), 4 * ones(Int, 2), 2 * ones(Int, 2)],
    "cycle_length" => 12,
}

# Coverage requires 49 man hours

bike_monday_employees = {
    "Cody" => {
        "min" => 3,
        "max" => 8,
        "availability" => ones(Int, 12),
    },
    "Damian" => {
        "min" => 3,
        "max" => 8,
        "availability" => ones(Int, 12),
    },
    "Douglas" => {
        "min" => 3,
        "max" => 8,
        "availability" => ones(Int, 12),
    },
    "Desmond" => {
        "min" => 3,
        "max" => 8,
        "availability" => ones(Int, 12),
    },
    "Dylan" => {
        "min" => 3,
        "max" => 8,
        "availability" => [
            ones(Int, 7), zeros(Int, 5),
        ],
    },
    "Jaeger" => {
        "min" => 3,
        "max" => 8,
        "availability" => ones(Int, 12),
    },
    "Leonardo" => {
        "min" => 3,
        "max" => 8,
        "availability" => [
            ones(Int, 7), zeros(Int, 5),
        ],
    },
    "Shaun" => {
        "min" => 3, # hours worked over simulation
        "max" => 8, # hours worked over simulation
        "availability" => [
            ones(Int, 12),
        ],
    },
}

max_stars_repo_name: shdann/autoshift

max_issues_repo_name: shdann/autoshift

max_forks_repo_name: shdann/autoshift


text:
function setup_filedataset_test()
    files = [
        "root/p1/f1.csv",
        "root/p2/f2.csv",
        "root/p2/p2p1/f2.csv",
        "root/p3/p3p1/f1.csv"
    ]

    for (i, file) in enumerate(files)
        paths = splitpath(file)[1:(end - 1)]
        root = ""
        for p in paths
            fullp = joinpath(root, p)
            isdir(fullp) || mkdir(fullp)
            root = fullp
        end

        open(file, "w") do io
            write(io, "a,b,c\n")
            write(io, join(i .* [1, 2, 3], ","))
        end
    end

    return files
end
cleanup_filedataset_test() = rm("root"; recursive = true)

@testset "FileDataset" begin
    files = setup_filedataset_test()
    dataset = FileDataset(f -> CSV.read(f, DataFrame), "root", "*.csv")
    @test numobs(dataset) == length(files)
    for (i, file) in enumerate(files)
        true_obs = CSV.read(file, DataFrame)
        @test getobs(dataset, i) == true_obs
    end
    cleanup_filedataset_test()
end


text:
using DroneSurveillance
using POMDPs
# import a solver from POMDPs.jl e.g. SARSOP
using SARSOP
# for visualization
using POMDPGifs
import Cairo, Fontconfig

pomdp = DroneSurveillancePOMDP() # initialize the problem
solver = SARSOPSolver(precision=0.1, timeout=10.0) # configure the solver
policy = solve(solver, pomdp) # solve the problem

sim = GifSimulator(filename="out.gif", max_steps=30)
simulate(sim, pomdp, policy);


text:
using JuliaGHA
using Test

using Logging: global_logger, current_logger
using LoggingExtras
using GitHubActions
get(ENV, "GITHUB_ACTIONS", "false") == "true" &&
    global_logger(TeeLogger(current_logger(), GitHubActions.GitHubActionsLogger()))

using Pkg
Pkg.develop(url = "https://github.com/adamslc/LoggingTestSets.jl")
using LoggingTestSets

@testset LoggingTestSet "JuliaGHA.jl" begin
    @testset "add" begin
        @test add(1, 1) == 2
        @test add(2, 2) == 4
        @test add(3, 3) == 6
    end

    @testset "mul" begin
        @test mul(1, 1) == 1
        @test mul(2, 2) == 4
        @test mul(3, 3) == 9
    end
end


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

module TestBoundsCheck

using Base.Test

@enum BCOption bc_default bc_on bc_off
bc_opt = BCOption(Base.JLOptions().check_bounds)

# test for boundscheck block eliminated at same level
@inline function A1()
    r = 0
    @boundscheck r += 1
    return r
end

@noinline function A1_noinline()
    r = 0
    @boundscheck r += 1
    return r
end

function A1_inbounds()
    r = 0
    @inbounds begin
        @boundscheck r += 1
    end
    return r
end

if bc_opt == bc_default
    @test A1() == 1
    @test A1_inbounds() == 0
elseif bc_opt == bc_on
    @test A1() == 1
    @test A1_inbounds() == 1
else
    @test A1() == 0
    @test A1_inbounds() == 0
end

# test for boundscheck block eliminated one layer deep, if the called method is inlined
@inline function A2()
    r = A1()+1
    return r
end

function A2_inbounds()
    @inbounds r = A1()+1
    return r
end

function A2_notinlined()
    @inbounds r = A1_noinline()+1
    return r
end

Base.@propagate_inbounds function A2_propagate_inbounds()
    r = A1()+1
    return r
end

if bc_opt == bc_default
    @test A2() == 2
    @test A2_inbounds() == 1
    @test A2_notinlined() == 2
    @test A2_propagate_inbounds() == 2
elseif bc_opt == bc_on
    @test A2() == 2
    @test A2_inbounds() == 2
    @test A2_notinlined() == 2
    @test A2_propagate_inbounds() == 2
else
    @test A2() == 1
    @test A2_inbounds() == 1
    @test A2_notinlined() == 1
    @test A2_propagate_inbounds() == 1
end

# test boundscheck NOT eliminated two layers deep, unless propagated

function A3()
    r = A2()+1
    return r
end

function A3_inbounds()
    @inbounds r = A2()+1
    return r
end

function A3_inbounds2()
    @inbounds r = A2_propagate_inbounds()+1
    return r
end

if bc_opt == bc_default
    @test A3() == 3
    @test A3_inbounds() == 3
    @test A3_inbounds2() == 2
elseif bc_opt == bc_on
    @test A3() == 3
    @test A3_inbounds() == 3
    @test A3_inbounds2() == 3
else
    @test A3() == 2
    @test A3_inbounds() == 2
    @test A3_inbounds2() == 2
end

# swapped nesting order of @boundscheck and @inbounds
function A1_nested()
    r = 0
    @boundscheck @inbounds r += 1
    return r
end

if bc_opt == bc_default || bc_opt == bc_on
    @test A1_nested() == 1
else
    @test A1_nested() == 0
end

# elide a throw
cb(x) = x > 0 || throw(BoundsError())

function B1()
    y = [1,2,3]
    @inbounds begin
        @boundscheck cb(0)
    end
    return 0
end

if bc_opt == bc_default || bc_opt == bc_off
    @test B1() == 0
else
    @test_throws BoundsError B1()
end

# elide a simple branch
cond(x) = x > 0 ? x : -x

function B2()
    y = [1,2,3]
    @inbounds begin
        @boundscheck cond(0)
    end
    return 0
end

@test B2() == 0

# Make sure type inference doesn't incorrectly optimize out
# `Expr(:inbounds, false)`
# Simply `return a[1]` doesn't work due to inlining bug
@inline function f1(a)
    # This has to be an arrayget / arrayset since these currently have a
    # implicit `Expr(:boundscheck)` that's not visible to type inference
    x = a[1]
    return x
end
# second level
@inline function g1(a)
    x = f1(a)
    return x
end
function k1(a)
    # This `Expr(:inbounds, true)` shouldn't affect `f1`
    @inbounds x = g1(a)
    return x
end
if bc_opt != bc_off
    @test_throws BoundsError k1(Int[])
end

# Ensure that broadcast doesn't use @inbounds when calling the function
if bc_opt != bc_off
    let A = zeros(3,3)
        @test_throws BoundsError broadcast(getindex, A, 1:3, 1:3)
    end
end

# issue #19554
function f19554(a)
    a[][3]
end
function f19554_2(a, b)
    a[][3] = b
    return a
end
a19554 = Ref{Array{Float64}}([1 2; 3 4])
@test f19554(a19554) === 2.0
@test f19554_2(a19554, 1) === a19554
@test a19554[][3] === f19554(a19554) === 1.0

# Ensure unsafe_view doesn't check bounds
function V1()
    A = rand(10,10)
    B = view(A, 4:7, 4:7)
    C = Base.unsafe_view(B, -2:7, -2:7)
    @test C == A
    nothing
end

if bc_opt == bc_default || bc_opt == bc_off
    @test V1() === nothing
else
    @test_throws BoundsError V1()
end

end

max_stars_repo_name: liveinjs/julia

max_issues_repo_name: liveinjs/julia

max_forks_repo_name: liveinjs/julia


text:
using AbstractDifferentiation
using Test, FiniteDifferences, LinearAlgebra
using ForwardDiff
using AbstractDifferentiation: ForwardDiffBackend
using Zygote
using Random
Random.seed!(1234)

const FDM = FiniteDifferences

## FiniteDifferences
struct FDMBackend1{A} <: AD.AbstractFiniteDifference
    alg::A
end
FDMBackend1() = FDMBackend1(central_fdm(5, 1))
const fdm_backend1 = FDMBackend1()
# Minimal interface
AD.@primitive function jacobian(ab::FDMBackend1, f, xs...)
    return FDM.jacobian(ab.alg, f, xs...)
end

struct FDMBackend2{A} <: AD.AbstractFiniteDifference
    alg::A
end
FDMBackend2() = FDMBackend2(central_fdm(5, 1))
const fdm_backend2 = FDMBackend2()
AD.@primitive function pushforward_function(ab::FDMBackend2, f, xs...)
    return function (vs)
        FDM.jvp(ab.alg, f, tuple.(xs, vs)...)
    end
end

struct FDMBackend3{A} <: AD.AbstractFiniteDifference
    alg::A
end
FDMBackend3() = FDMBackend3(central_fdm(5, 1))
const fdm_backend3 = FDMBackend3()
AD.@primitive function pullback_function(ab::FDMBackend3, f, xs...)
    return function (vs)
        # Supports only single output
        if vs isa AbstractVector
            return FDM.j′vp(ab.alg, f, vs, xs...)
        else
            @assert length(vs) == 1
            return FDM.j′vp(ab.alg, f, vs[1], xs...)

        end
    end
end
##

## ForwardDiff
const forwarddiff_backend = ForwardDiffBackend()
struct ForwardDiffBackend1 <: AD.AbstractForwardMode end
const forwarddiff_backend1 = ForwardDiffBackend1()
AD.@primitive function jacobian(ab::ForwardDiffBackend1, f, xs)
    if xs isa Number
        return (ForwardDiff.derivative(f, xs),)
    elseif xs isa AbstractArray
        out = f(xs)
        if out isa Number
            return (adjoint(ForwardDiff.gradient(f, xs)),)
        else
            return (ForwardDiff.jacobian(f, xs),)
        end
    elseif xs isa Tuple
        error(typeof(xs))
    else
        error(typeof(xs))
    end
end
AD.primal_value(::ForwardDiffBackend1, ::Any, f, xs) = ForwardDiff.value.(f(xs...))

struct ForwardDiffBackend2 <: AD.AbstractForwardMode end
const forwarddiff_backend2 = ForwardDiffBackend2()
AD.@primitive function pushforward_function(ab::ForwardDiffBackend2, f, xs...)
    # jvp = f'(x)*v, i.e., differentiate f(x + h*v) wrt h at 0
    return function (vs)
        if xs isa Tuple
            @assert length(xs) <= 2
            if length(xs) == 1
                (ForwardDiff.derivative(h->f(xs[1]+h*vs[1]),0),)
            else
                ForwardDiff.derivative(h->f(xs[1]+h*vs[1], xs[2]+h*vs[2]),0)
            end
        else
            ForwardDiff.derivative(h->f(xs+h*vs),0)
        end
    end
end
AD.primal_value(::ForwardDiffBackend2, ::Any, f, xs) = ForwardDiff.value.(f(xs...))
##

## Zygote
struct ZygoteBackend1 <: AD.AbstractReverseMode end
const zygote_backend1 = ZygoteBackend1()
AD.@primitive function pullback_function(ab::ZygoteBackend1, f, xs...)
    return function (vs)
        # Supports only single output
        _, back = Zygote.pullback(f, xs...)
        if vs isa AbstractVector
            back(vs)
        else
            @assert length(vs) == 1
            back(vs[1])
        end
    end
end
##

fder(x, y) = exp(y) * x + y * log(x)
dfderdx(x, y) = exp(y) + y * 1/x
dfderdy(x, y) = exp(y) * x + log(x)

fgrad(x, y) = prod(x) + sum(y ./ (1:length(y)))
dfgraddx(x, y) = prod(x)./x
dfgraddy(x, y) = one(eltype(y)) ./ (1:length(y))
dfgraddxdx(x, y) = prod(x)./(x*x') - Diagonal(diag(prod(x)./(x*x')))
dfgraddydy(x, y) = zeros(length(y),length(y))

function fjac(x, y)
    x + -3*y + [y[2:end];zero(y[end])]/2# Bidiagonal(-ones(length(y)) * 3, ones(length(y) - 1) / 2, :U) * y
end
function dfjacdx(x, y)
    if VERSION < v"1.3"
        return Matrix{Float64}(I, length(x), length(x))
    else
        return I(length(x))
    end
end
dfjacdy(x, y) = Bidiagonal(-ones(length(y)) * 3, ones(length(y) - 1) / 2, :U)

# Jvp
jxvp(x,y,v) = dfjacdx(x,y)*v
jyvp(x,y,v) = dfjacdy(x,y)*v

# vJp
vJxp(x,y,v) = dfjacdx(x,y)'*v
vJyp(x,y,v) = dfjacdy(x,y)'*v

const xscalar = rand()
const yscalar = rand()

const xvec = rand(5)
const yvec = rand(5)

# to check if vectors get mutated
xvec2 = deepcopy(xvec)
yvec2 = deepcopy(yvec)

function test_higher_order_backend(backends...)
    ADbackends = AD.HigherOrderBackend(backends)
    @test backends[end] == AD.lowest(ADbackends)
    @test backends[end-1] == AD.second_lowest(ADbackends)

    for i in length(backends):-1:1
        @test backends[i] == AD.lowest(ADbackends)
        ADbackends = AD.reduce_order(ADbackends)
    end
    backends[1] == AD.reduce_order(ADbackends)
end

function test_derivatives(backend; multiple_inputs=true)
    # test with respect to analytical solution
    der_exact = (dfderdx(xscalar,yscalar), dfderdy(xscalar,yscalar))
    if multiple_inputs
        der1 = AD.derivative(backend, fder, xscalar, yscalar)
        @test minimum(isapprox.(der_exact, der1, rtol=1e-10))
        valscalar, der2 = AD.value_and_derivative(backend, fder, xscalar, yscalar)
        @test valscalar == fder(xscalar, yscalar)
        @test der2 .- der1 == (0, 0)
    end
    # test if single input (no tuple works)
    valscalara, dera = AD.value_and_derivative(backend, x -> fder(x, yscalar), xscalar)
    valscalarb, derb = AD.value_and_derivative(backend, y -> fder(xscalar, y), yscalar)
    @test fder(xscalar, yscalar) == valscalara
    @test fder(xscalar, yscalar) == valscalarb
    @test isapprox(dera[1], der_exact[1], rtol=1e-10)
    @test isapprox(derb[1], der_exact[2], rtol=1e-10)
end

function test_gradients(backend; multiple_inputs=true)
    # test with respect to analytical solution
    grad_exact = (dfgraddx(xvec,yvec), dfgraddy(xvec,yvec))
    if multiple_inputs
        grad1 = AD.gradient(backend, fgrad, xvec, yvec)
        @test minimum(isapprox.(grad_exact, grad1, rtol=1e-10))
        valscalar, grad2 = AD.value_and_gradient(backend, fgrad, xvec, yvec)
        @test valscalar == fgrad(xvec, yvec)
        @test norm.(grad2 .- grad1) == (0, 0)
    end
    # test if single input (no tuple works)
    valscalara, grada = AD.value_and_gradient(backend, x -> fgrad(x, yvec), xvec)
    valscalarb, gradb = AD.value_and_gradient(backend, y -> fgrad(xvec, y), yvec)
    @test fgrad(xvec, yvec) == valscalara
    @test fgrad(xvec, yvec) == valscalarb
    @test isapprox(grada[1], grad_exact[1], rtol=1e-10)
    @test isapprox(gradb[1], grad_exact[2], rtol=1e-10)
    @test xvec == xvec2
    @test yvec == yvec2
end

function test_jacobians(backend; multiple_inputs=true)
    # test with respect to analytical solution
    jac_exact = (dfjacdx(xvec, yvec), dfjacdy(xvec, yvec))
    if multiple_inputs
        jac1 = AD.jacobian(backend, fjac, xvec, yvec)
        @test minimum(isapprox.(jac_exact, jac1, rtol=1e-10))
        valvec, jac2 = AD.value_and_jacobian(backend, fjac, xvec, yvec)
        @test valvec == fjac(xvec, yvec)
        @test norm.(jac2 .- jac1) == (0, 0)
    end

    # test if single input (no tuple works)
    valveca, jaca = AD.value_and_jacobian(backend, x -> fjac(x, yvec), xvec)
    valvecb, jacb = AD.value_and_jacobian(backend, y -> fjac(xvec, y), yvec)
    @test fjac(xvec, yvec) == valveca
    @test fjac(xvec, yvec) == valvecb
    @test isapprox(jaca[1], jac_exact[1], rtol=1e-10)
    @test isapprox(jacb[1], jac_exact[2], rtol=1e-10)
    @test xvec == xvec2
    @test yvec == yvec2
end

function test_hessians(backend; multiple_inputs=false)
    if multiple_inputs
        # ... but
        error("multiple_inputs=true is not supported.")
    else
        # explicit test that AbstractDifferentiation throws an error
        # don't support tuple of Hessians
        @test_throws AssertionError H1 = AD.hessian(backend, fgrad, (xvec, yvec))
        @test_throws MethodError H1 = AD.hessian(backend, fgrad, xvec, yvec)
    end

    # @test dfgraddxdx(xvec,yvec) ≈ H1[1] atol=1e-10
    # @test dfgraddydy(xvec,yvec) ≈ H1[2] atol=1e-10

    # test if single input (no tuple works)
    fhess = x -> fgrad(x, yvec)
    hess1 = AD.hessian(backend, fhess, xvec)
    # test with respect to analytical solution
    @test dfgraddxdx(xvec,yvec) ≈ hess1[1] atol=1e-10

    valscalar, hess2 = AD.value_and_hessian(backend, fhess, xvec)
    @test valscalar == fgrad(xvec, yvec)
    @test norm.(hess2 .- hess1) == (0,)
    valscalar, grad, hess3 = AD.value_gradient_and_hessian(backend, fhess, xvec)
    @test valscalar == fgrad(xvec, yvec)
    @test norm.(grad .- AD.gradient(backend, fhess, xvec)) == (0,)
    @test norm.(hess3 .- hess1) == (0,)

    @test xvec == xvec2
    @test yvec == yvec2
    fhess2 = x-> dfgraddx(x, yvec)
    hess4 = AD.jacobian(backend, fhess2, xvec)
    @test minimum(isapprox.(hess4, hess1, atol=1e-10))
end

function test_jvp(backend; multiple_inputs=true)
    v = (rand(length(xvec)), rand(length(yvec)))

    if multiple_inputs
        if backend isa Union{FDMBackend2,ForwardDiffBackend,ForwardDiffBackend2} # augmented version of v
            identity_like = AD.identity_matrix_like(v)
            vaug = map(identity_like) do identity_like_i
                identity_like_i .* v
            end

            pf1 = map(v->AD.pushforward_function(backend, fjac, xvec, yvec)(v), vaug)
            ((valvec1, pf2x), (valvec2, pf2y)) = map(v->AD.value_and_pushforward_function(backend, fjac, xvec, yvec)(v), vaug)
        else
            pf1 = AD.pushforward_function(backend, fjac, xvec, yvec)(v)
            valvec, pf2 = AD.value_and_pushforward_function(backend, fjac, xvec, yvec)(v)
            ((valvec1, pf2x), (valvec2, pf2y)) = (valvec, pf2[1]), (valvec, pf2[2])
        end

        @test valvec1 == fjac(xvec, yvec)
        @test valvec2 == fjac(xvec, yvec)
        @test norm.((pf2x,pf2y) .- pf1) == (0, 0)
        # test with respect to analytical solution
        @test minimum(isapprox.(pf1, (jxvp(xvec,yvec,v[1]), jyvp(xvec,yvec,v[2])), atol=1e-10))
        @test xvec == xvec2
        @test yvec == yvec2
    end

    valvec1, pf1 = AD.value_and_pushforward_function(backend, x -> fjac(x, yvec), xvec)(v[1])
    valvec2, pf2 = AD.value_and_pushforward_function(backend, y -> fjac(xvec, y), yvec)(v[2])

    if backend isa Union{FDMBackend2}
      pf1 = (pf1,)
      pf2 = (pf2,)
    end
    @test valvec1 == fjac(xvec, yvec)
    @test valvec2 == fjac(xvec, yvec)
    @test minimum(isapprox.((pf1[1],pf2[1]), (jxvp(xvec,yvec,v[1]), jyvp(xvec,yvec,v[2])), atol=1e-10))
end

function test_j′vp(backend; multiple_inputs=true)
    # test with respect to analytical solution
    w = rand(length(fjac(xvec, yvec)))
    if multiple_inputs
        pb1 = AD.pullback_function(backend, fjac, xvec, yvec)(w)
        valvec, pb2 = AD.value_and_pullback_function(backend, fjac, xvec, yvec)(w)
        @test valvec == fjac(xvec, yvec)
        @test norm.(pb2 .- pb1) == (0, 0)
        @test minimum(isapprox.(pb1, (vJxp(xvec,yvec,w), vJyp(xvec,yvec,w)), atol=1e-10))
        @test xvec == xvec2
        @test yvec == yvec2
    end

    valvec1, pb1 = AD.value_and_pullback_function(backend, x -> fjac(x, yvec), xvec)(w)
    valvec2, pb2 = AD.value_and_pullback_function(backend, y -> fjac(xvec, y), yvec)(w)
    @test valvec1 == fjac(xvec, yvec)
    @test valvec2 == fjac(xvec, yvec)
    @test minimum(isapprox.((pb1[1],pb2[1]), (vJxp(xvec,yvec,w), vJyp(xvec,yvec,w)), atol=1e-10))
end

function test_lazy_derivatives(backend; multiple_inputs=true)
    # single input function
    der1 = AD.derivative(backend, x->fder(x, yscalar), xscalar)
    lazyder = AD.LazyDerivative(backend, x->fder(x, yscalar), xscalar)

    # multiplication with scalar
    @test lazyder*yscalar == der1.*yscalar
    @test lazyder*yscalar isa Tuple

    @test yscalar*lazyder == yscalar.*der1
    @test yscalar*lazyder isa Tuple

    # multiplication with array
    @test lazyder*yvec == (der1.*yvec,)
    @test lazyder*yvec isa Tuple

    @test yvec*lazyder == (yvec.*der1,)
    @test yvec*lazyder isa Tuple

    # multiplication with tuple
    @test lazyder*(yscalar,) == lazyder*yscalar
    @test lazyder*(yvec,) == lazyder*yvec

    @test (yscalar,)*lazyder == yscalar*lazyder
    @test (yvec,)*lazyder == yvec*lazyder

    # two input function
    if multiple_inputs
        der1 = AD.derivative(backend, fder, xscalar, yscalar)
        lazyder = AD.LazyDerivative(backend, fder, (xscalar, yscalar))

        # multiplication with scalar
        @test lazyder*yscalar == der1.*yscalar
        @test lazyder*yscalar isa Tuple

        @test yscalar*lazyder == yscalar.*der1
        @test yscalar*lazyder isa Tuple

        # multiplication with array
        @test lazyder*yvec == (der1[1]*yvec, der1[2]*yvec)
        @test lazyder*yvec isa Tuple

        @test yvec*lazyder == (yvec*der1[1], yvec*der1[2])
        @test lazyder*yvec isa Tuple

        # multiplication with tuple
        @test_throws AssertionError lazyder*(yscalar,)
        @test_throws AssertionError lazyder*(yvec,)

        @test_throws AssertionError (yscalar,)*lazyder
        @test_throws AssertionError (yvec,)*lazyder
    end
end

function test_lazy_gradients(backend; multiple_inputs=true)
    # single input function
    grad1 = AD.gradient(backend, x->fgrad(x, yvec), xvec)
    lazygrad = AD.LazyGradient(backend, x->fgrad(x, yvec), xvec)

    # multiplication with scalar
    @test norm.(lazygrad*yscalar .- grad1.*yscalar) == (0,)
    @test lazygrad*yscalar isa Tuple

    @test norm.(yscalar*lazygrad .- yscalar.*grad1) == (0,)
    @test yscalar*lazygrad isa Tuple

    # multiplication with tuple
    @test lazygrad*(yscalar,) == lazygrad*yscalar
    @test (yscalar,)*lazygrad == yscalar*lazygrad

    # two input function
    if multiple_inputs
        grad1 = AD.gradient(backend, fgrad, xvec, yvec)
        lazygrad = AD.LazyGradient(backend, fgrad, (xvec, yvec))

        # multiplication with scalar
        @test norm.(lazygrad*yscalar .- grad1.*yscalar) == (0,0)
        @test lazygrad*yscalar isa Tuple

        @test norm.(yscalar*lazygrad .- yscalar.*grad1) == (0,0)
        @test yscalar*lazygrad isa Tuple

        # multiplication with tuple
        @test_throws AssertionError lazygrad*(yscalar,) == lazygrad*yscalar
        @test_throws AssertionError (yscalar,)*lazygrad == yscalar*lazygrad
    end
end

function test_lazy_jacobians(backend; multiple_inputs=true)
    # single input function
    jac1 = AD.jacobian(backend, x->fjac(x, yvec), xvec)
    lazyjac = AD.LazyJacobian(backend, x->fjac(x, yvec), xvec)

    # multiplication with scalar
    @test norm.(lazyjac*yscalar .- jac1.*yscalar) == (0,)
    @test lazyjac*yscalar isa Tuple

    @test norm.(yscalar*lazyjac .- yscalar.*jac1) == (0,)
    @test yscalar*lazyjac isa Tuple

    w = rand(length(fjac(xvec, yvec)))
    v = (rand(length(xvec)),rand(length(xvec)))

    # vjp
    pb1 = (vJxp(xvec,yvec,w),)
    res = w'*lazyjac
    @test minimum(isapprox.(pb1, res, atol=1e-10))
    @test res isa Tuple

    # jvp
    pf1 = (jxvp(xvec,yvec,v[1]),)
    res = lazyjac*v[1]
    @test minimum(isapprox.(pf1, res, atol=1e-10))
    @test res isa Tuple

    # two input function
    if multiple_inputs
        jac1 = AD.jacobian(backend, fjac, xvec, yvec)
        lazyjac = AD.LazyJacobian(backend, fjac, (xvec, yvec))

        # multiplication with scalar
        @test norm.(lazyjac*yscalar .- jac1.*yscalar) == (0,0)
        @test lazyjac*yscalar isa Tuple

        @test norm.(yscalar*lazyjac .- yscalar.*jac1) == (0,0)
        @test yscalar*lazyjac isa Tuple

        # vjp
        pb1 = (vJxp(xvec,yvec,w), vJyp(xvec,yvec,w))
        res = w'lazyjac
        @test minimum(isapprox.(pb1, res, atol=1e-10))
        @test res isa Tuple

        # jvp
        pf1 = (jxvp(xvec,yvec,v[1]), jyvp(xvec,yvec,v[2]))

        if backend isa Union{FDMBackend2,ForwardDiffBackend,ForwardDiffBackend2} # augmented version of v
            identity_like = AD.identity_matrix_like(v)
            vaug = map(identity_like) do identity_like_i
                identity_like_i .* v
            end

            res = map(v->(lazyjac*v)[1], vaug)
        else
            res = lazyjac*v
        end
        @test minimum(isapprox.(pf1, res, atol=1e-10))
        @test res isa Tuple
    end
end

function test_lazy_hessians(backend; multiple_inputs=true)
    # fdm_backend not used here yet..
    # single input function
    fhess = x -> fgrad(x, yvec)
    hess1 = (dfgraddxdx(xvec,yvec),)
    lazyhess = AD.LazyHessian(backend, fhess, xvec)

    # multiplication with scalar
    @test minimum(isapprox.(lazyhess*yscalar, hess1.*yscalar, atol=1e-10))
    @test lazyhess*yscalar isa Tuple

    # multiplication with scalar
    @test minimum(isapprox.(yscalar*lazyhess, yscalar.*hess1, atol=1e-10))
    @test yscalar*lazyhess isa Tuple

    w = rand(length(xvec))
    v = rand(length(xvec))

    # Hvp
    Hv = map(h->h*v, hess1)
    res = lazyhess*v
    @test minimum(isapprox.(Hv, res, atol=1e-10))
    @test res isa Tuple

    # H′vp
    wH = map(h->h'*w, hess1)
    res = w'*lazyhess
    @test minimum(isapprox.(wH, res, atol=1e-10))
    @test res isa Tuple
end

@testset "AbstractDifferentiation.jl" begin
    @testset "Utils" begin
        test_higher_order_backend(fdm_backend1, fdm_backend2, fdm_backend3, zygote_backend1, forwarddiff_backend, forwarddiff_backend2)
    end
    @testset "FiniteDifferences" begin
        @testset "Derivative" begin
            test_derivatives(fdm_backend1)
            test_derivatives(fdm_backend2)
            test_derivatives(fdm_backend3)
        end
        @testset "Gradient" begin
            test_gradients(fdm_backend1)
            test_gradients(fdm_backend2)
            test_gradients(fdm_backend3)
        end
        @testset "Jacobian" begin
            test_jacobians(fdm_backend1)
            test_jacobians(fdm_backend2)
            test_jacobians(fdm_backend3)
        end
        @testset "Hessian" begin
            test_hessians(fdm_backend1)
            test_hessians(fdm_backend2)
            test_hessians(fdm_backend3)
        end
        @testset "jvp" begin
            test_jvp(fdm_backend1)
            test_jvp(fdm_backend2)
            test_jvp(fdm_backend3)
        end
        @testset "j′vp" begin
            test_j′vp(fdm_backend1)
            test_j′vp(fdm_backend2)
            test_j′vp(fdm_backend3)
        end
        @testset "Lazy Derivative" begin
            test_lazy_derivatives(fdm_backend1)
            test_lazy_derivatives(fdm_backend2)
            test_lazy_derivatives(fdm_backend3)
        end
        @testset "Lazy Gradient" begin
            test_lazy_gradients(fdm_backend1)
            test_lazy_gradients(fdm_backend2)
            test_lazy_gradients(fdm_backend3)
        end
        @testset "Lazy Jacobian" begin
            test_lazy_jacobians(fdm_backend1)
            test_lazy_jacobians(fdm_backend2)
            test_lazy_jacobians(fdm_backend3)
        end
        @testset "Lazy Hessian" begin
            test_lazy_hessians(fdm_backend1)
            test_lazy_hessians(fdm_backend2)
            test_lazy_hessians(fdm_backend3)
        end
    end
    @testset "ForwardDiff" begin
        @testset "Derivative" begin
            test_derivatives(forwarddiff_backend)
            test_derivatives(forwarddiff_backend1; multiple_inputs=false)
            test_derivatives(forwarddiff_backend2)
        end
        @testset "Gradient" begin
            test_gradients(forwarddiff_backend)
            test_gradients(forwarddiff_backend1; multiple_inputs=false)
            test_gradients(forwarddiff_backend2)
        end
        @testset "Jacobian" begin
            test_jacobians(forwarddiff_backend)
            test_jacobians(forwarddiff_backend1; multiple_inputs=false)
            test_jacobians(forwarddiff_backend2)
        end
        @testset "Hessian" begin
            test_hessians(forwarddiff_backend)
            test_hessians(forwarddiff_backend1; multiple_inputs=false)
            test_hessians(forwarddiff_backend2)
        end
        @testset "jvp" begin
            test_jvp(forwarddiff_backend)
            test_jvp(forwarddiff_backend1; multiple_inputs=false)
            test_jvp(forwarddiff_backend2)
        end
        @testset "j′vp" begin
            test_j′vp(forwarddiff_backend)
            test_j′vp(forwarddiff_backend1; multiple_inputs=false)
            test_j′vp(forwarddiff_backend2)
        end
        @testset "Lazy Derivative" begin
            test_lazy_derivatives(forwarddiff_backend)
            test_lazy_derivatives(forwarddiff_backend1; multiple_inputs=false)
            test_lazy_derivatives(forwarddiff_backend2)
        end
        @testset "Lazy Gradient" begin
            test_lazy_gradients(forwarddiff_backend)
            test_lazy_gradients(forwarddiff_backend1; multiple_inputs=false)
            test_lazy_gradients(forwarddiff_backend2)
        end
        @testset "Lazy Jacobian" begin
            test_lazy_jacobians(forwarddiff_backend)
            test_lazy_jacobians(forwarddiff_backend1; multiple_inputs=false)
            test_lazy_jacobians(forwarddiff_backend2)
        end
        @testset "Lazy Hessian" begin
            test_lazy_hessians(forwarddiff_backend)
            test_lazy_hessians(forwarddiff_backend1; multiple_inputs=false)
            test_lazy_hessians(forwarddiff_backend2)
        end
    end
    @testset "Zygote" begin
        @testset "Derivative" begin
            test_derivatives(zygote_backend1)
        end
        @testset "Gradient" begin
            test_gradients(zygote_backend1)
        end
        @testset "Jacobian" begin
            test_jacobians(zygote_backend1)
        end
        @testset "Hessian" begin
            # Zygote over Zygote problems
            backends = AD.HigherOrderBackend((forwarddiff_backend,zygote_backend1))
            test_hessians(backends)
            backends = AD.HigherOrderBackend((forwarddiff_backend2,zygote_backend1))
            test_hessians(backends)
            if VERSION >= v"1.3"
                # backends = AD.HigherOrderBackend((zygote_backend1,forwarddiff_backend))
                # test_hessians(backends)
                backends = AD.HigherOrderBackend((zygote_backend1,forwarddiff_backend1))
                test_hessians(backends)
            end
            # fails:
            # backends = AD.HigherOrderBackend((zygote_backend1,forwarddiff_backend2))
            # test_hessians(backends)
        end
        @testset "jvp" begin
            test_jvp(zygote_backend1)
        end
        @testset "j′vp" begin
            test_j′vp(zygote_backend1)
        end
        @testset "Lazy Derivative" begin
            test_lazy_derivatives(zygote_backend1)
        end
        @testset "Lazy Gradient" begin
            test_lazy_gradients(zygote_backend1)
        end
        @testset "Lazy Jacobian" begin
            test_lazy_jacobians(zygote_backend1)
        end
        @testset "Lazy Hessian" begin
            # Zygote over Zygote problems
            backends = AD.HigherOrderBackend((forwarddiff_backend,zygote_backend1))
            test_lazy_hessians(backends)
            backends = AD.HigherOrderBackend((forwarddiff_backend2,zygote_backend1))
            test_lazy_hessians(backends)
            if VERSION >= v"1.3"
                # backends = AD.HigherOrderBackend((zygote_backend1,forwarddiff_backend))
                # test_lazy_hessians(backends)
                backends = AD.HigherOrderBackend((zygote_backend1,forwarddiff_backend1))
                test_lazy_hessians(backends)
            end
        end
    end
end


text:
# This file is a part of SimilaritySearch.jl

export neardup

"""
    neardup(idx::AbstractSearchContext, X::AbstractVector, ϵ)

Find nearest duplicates in database `X` using the empty index `idx`. The algorithm iteratively try to index elements in `X`,
and items being near than `ϵ` to some element in `idx` will be ignored.

The function returns a named tuple `(idx, map, nn, dist)` where:
- `idx`: it is the index of the non duplicated elements
- `map`: a mapping from `1-|idx|` to its positions in `X`
- `nn`: an array where each element in ``x \\in X`` points to its covering element (previously indexed element `u` such that ``d(u, x_i) \\leq ϵ``)
- `dist`: an array of distance values to each covering element (correspond to each element in `nn`)

# Arguments
- `idx`: An empty index (i.e., a `SearchGraph`)
- `X`: The input dataset
- `ϵ`: Real value to cut

# Notes
- The index `idx` must support incremental construction with `push!`
- You can access the set of elements being 'ϵ'-non duplicates (the ``ϵ-net``) using `idx.db` or where `nn[i] == i`
"""
function neardup(idx::AbstractSearchContext, X::AbstractDatabase, ϵ::Real)
    res = KnnResult(3)  # should be 1, but index's setups work better on larger `k` values
    L = zeros(Int32, length(X))
    D = zeros(Float32, length(X))
    M = Int32[1]
    L[1] = 1
    D[1] = 0.0
    push!(idx, X[1])
    @inbounds for i in 2:length(X)
        res = reuse!(res)
        x = X[i]
        search(idx, x, res)
        if length(res) == 0 || minimum(res) > ϵ
            push!(idx, x)
            L[i] = i
            D[i] = 0.0
            push!(M, i)
        else
            L[i] = M[argmin(res)]
            D[i] = minimum(res)

        end
    end

    (idx=idx, map=M, nn=L, dist=D)
end

max_stars_repo_name: sadit/SimilaritySearch

max_issues_repo_name: sadit/NNS

max_forks_repo_name: sadit/SimilaritySearch


text:
abstract type UMRRuleModel <: MixingRule end

struct UMRRule{γ} <: UMRRuleModel
    components::Array{String,1}
    activity::γ
end

@registermodel UMRRule
export UMRRule
function UMRRule(components::Vector{String}; activity = Wilson, userlocations::Vector{String}=String[],activity_userlocations::Vector{String}=String[], verbose::Bool=false)
    init_activity = activity(components;userlocations = activity_userlocations,verbose)
    model = UMRRule(components, init_activity)
    return model
end

function ab_premixing(::Type{PR},mixing::UMRRuleModel,Tc,pc,kij)
    Ωa, Ωb = ab_consts(PR)
    _Tc = Tc.values
    _pc = pc.values
    a = epsilon_LorentzBerthelot(SingleParam(pc, @. Ωa*R̄^2*_Tc^2/_pc),kij)
    bi = @. Ωb*R̄*_Tc/_pc
    bij = ((bi.^(1/2).+bi'.^(1/2))/2).^2
    b = PairParam("b",Tc.components,bij)
    return a,b
end

function mixing_rule(model::PRModel,V,T,z,mixing_model::UMRRuleModel,α,a,b,c)
    n = sum(z)
    #x = z./n
    invn = (one(n)/n)
    invn2 = invn^2
    lnγ_SG_  = lnγ_SG(mixing_model.activity,1e5,T,z)
    lnγ_res_ = lnγ_res(mixing_model.activity,1e5,T,z)
    g_E = sum(z[i]*R̄*T*(lnγ_res_[i]+lnγ_SG_[i]) for i ∈ @comps)*invn
    #b = Diagonal(b).diag
    #b = ((b.^(1/2).+b'.^(1/2))/2).^2
    b̄ = dot(z,Symmetric(b),z) * invn2
    c̄ = dot(z,c)*invn
    Σab = sum(z[i]*a[i,i]*α[i]/b[i,i]/(R̄*T) for i ∈ @comps)*invn
    ā = b̄*R̄*T*(Σab-1/0.53*g_E/(R̄*T))
    return ā,b̄,c̄
end


text:
function center_of_mass(
    pos::Vector{T},
    box::MolecularBox{T,N,P};
    weights = ones(Float64,length(pos))
) where {T,N,P}
    if sum(P) != N
        error("Calculating the center of mass is currently only supported for"*
            "fully-periodic systems")
    end
    if length(pos)==1 ; return pos[1] ; end

    L = box.lengths

    length(weights) == length(pos) || error("Weights and positions mismatch")

    invtotweight = 1.0/sum(weights)
    i2pi = 0.5/pi

    a = zero(eltype(pos))
    b = zero(eltype(pos))

    for i in eachindex(pos)
        t = pos[i]./L*2pi
        a -= cos.(t) * weights[i]
        b -= sin.(t) * weights[i]
    end
    a = a*invtotweight
    b = b*invtotweight

    tcom = atan.(b,a)+pi
    com = tcom.*i2pi.*L
end