text:
using DataFrames
using CSV
using Random

##  FUNCTION build_dict IS USED TO QUICKLY BUILD DICTIONARIES FROM DATAFRAMES
function build_dict(df_in::DataFrame)

    sz = size(df_in)
    n = sz[2]

    if n == 2
        #in two columns, simply map a to b
        dict_out = Dict([x for x in zip(df_in[:, 1], df_in[:, 2])])
    else
        #if there are more than two columns, leading columns are tuples
        fields_key = names(df_in)[1:(n - 1)]
        tups = [tuple(x...) for x in eachrow(Matrix(df_in[:, fields_key]))]
        dict_out = Dict([x for x in zip(tups, df_in[:, names(df_in)[n]])])
    end

    return dict_out
end

##  function to clean table names

function clean_field_names(nms::Array{String, 1})
    nms = replace.(lowercase.(nms), "  " => " ")
    nms = replace.(nms, " " => "_")
    nms = replace.(nms, "\$" => "")
    nms = replace.(nms, "`" => "")
    nms = replace.(nms, "-" => "_")

    for i in 1:length(nms)
        nm = nms[i]
        while nm[1] in ["_", "-", "."]
            nm = string(nm[2:end])
        end

        nms[i] = nm
    end

    return nms
end
function clean_field_names(nms::Array{Symbol, 1})
    nms = replace.(lowercase.(String.(nms)), "  " => " ")
    nms = replace.(nms, " " => "_")
    nms = replace.(nms, "\$" => "")
    nms = replace.(nms, "`" => "")
    nms = replace.(nms, "-" => "_")

    for i in 1:length(nms)
        nm = nms[i]
        while nm[1] in ["_", "-", "."]
            nm = string(nm[2:end])
        end

        nms[i] = nm
    end

    return Symbol.(nms)
end

function clean_field_names!(df_in::DataFrame)
    nms0 = names(df_in)
    nms = clean_field_names(nms0)
    rename!(df_in, Dict(zip(Symbol.(nms0), Symbol.(nms))))
end

##  functions for checking fields/keys
function check_fields!(data_frame::DataFrame, fields_check::Array{Symbol, 1})
    # verify fields
    fields_req = Set(String.(fields_check))

    if !issubset(fields_req, Set(names(data_frame)))
        str_missing_fields = setdiff(Set(names(data_frame)), fields_req)
        str_missing_fields = join("'".*str_missing_fields.*"'", ", ")
        error("Fields $(str_missing_fields) not found. Check the data frame.")
    end
end

##  check to ensure keys are present in the dictionary
function check_keys!(dict::Dict, keys_check::Array{String, 1})
    # verify fields
    keys_req = Set(String.(keys_check))
    keys_avail = Set(String.(collect(keys(dict))))

    if !issubset(keys_req, keys_avail)
        str_missing_keys = setdiff(keys_avail, keys_req)
        str_missing_keys = join("'".*str_missing_keys.*"'", ", ")
        error("Keys $(str_missing_keys) not found. Check the dictionary.")
    end
end
function check_keys!(dict::Dict, keys_check::Array{Symbol, 1})
    # verify fields
    keys_req = Set(String.(keys_check))
    keys_avail = Set(String.(collect(keys(dict))))

    if !issubset(keys_req, keys_avail)
        str_missing_keys = setdiff(keys_avail, keys_req)
        str_missing_keys = join("'".*str_missing_keys.*"'", ", ")
        error("Keys $(str_missing_keys) not found. Check the dictionary.")
    end
end

##  use check_path to return a path if it exists and throw an error otherwise; optional 'create_directory_q' can be used to create a directory
function check_path(path::String, create_directory_q::Bool)
    if ispath(path)
        return path
    else
        if create_directory_q & !occursin(".", basename(path))
            mkdir(path)
            print("Created directory '$(path)'\n")
            return path
        else
            error("Path '$(path)' not found.")
        end
    end
end

##  FUNCTION check_trailing_num_field checks if the end of a field (split by "_") is numeric; allows for the return of that number
function check_trailing_num_field(field::String, return_val::String = "query", field_delim::String = "_")
    nv = split(field, field_delim)

        tp = nothing
    if length(nv) > 1
        tp = tryparse(Int64, nv[length(nv)])
        rv_bool = (tp !== nothing)
    else
        rv_bool = false
    end

    if return_val == "query"
        return rv_bool
    elseif return_val == "num"
        return tp
    end
end

##  FUNCTION FOR EXPORTING A SUBSET OF A DATAFRAME AND CONVERTING IT TO A TAR GZ
function export_sub_file_csv(df_in::DataFrame, fields_save::Array{String, 1}, fp_csv_exp::String, conv_gz_q::Bool = true, rm_csv_q::Bool = false)

    print("\n\n###    Exporting drop fields to '$(fp_csv_exp)'... \n")

    #export unused data to CSV
    CSV.write(fp_csv_exp, df_in[:, fields_save])
    print("\n###    Export complete. ")

    if conv_gz_q
        print("Creating compressed file...\n")
        fp_targz_exp = replace(fp_csv_exp, ".csv" => ".tar.gz")
        fn_targz_exp = basename(fp_targz_exp)

        #remove existing tar.gz
        if ispath(fp_targz_exp)
            rm(fp_targz_exp)
        end

        #build tar.gz file
        dir_cur_tmp = pwd()
        comm_tgz = `tar -cvzf "$(fn_targz_exp)" "$(basename(fp_csv_exp))"`
        cd(dirname(fp_csv_exp))
        run(comm_tgz)
        cd(dir_cur_tmp)

        print("\n###    Compressed file complete.")

        if rm_csv_q
            print(" Removing CSV...\n")
            #remove the csv
            rm(fp_csv_exp)
        end

        print("\n###    Done.\n\n")
    end
end

##  define some functions for retrieving tables
function get_table(mode::String, fp_csv::String, query::String)

    # switch between csv/sql quickly with the "mode"
    if mode == "csv"
        # assume that the header is available
        df_in = read_csv(fp_csv, true)
    elseif mode == "sql"
        # placeholder for integration with sql
        df_in = DataFrame()
    else
        error("Mode $(mode) not found for get_table. Check to ensure it is specified as 'csv' or 'sql'.")
    end

    return df_in
end

## use parse_config to read a configuration file and return a dictionary
function parse_config(fp_config::String)

    if ispath(fp_config)
        # read lines from the file
        file_read = open(fp_config, "r")
        lines = readlines(file_read)
        close(file_read)

        # initialize the output dictionary
        dict_out = Dict{String, Any}()
        for i in 1:length(lines)
            line = lines[i]
            if length(line) > 0
                # check conditions
                if string(line[1]) != "#"
                    # drop comments
                    line = split(line, "#")[1]
                    # next, split on colon
                    line = split(line, ":")
                    key = strip(string(line[1]))
                    value = strip(string(line[2]))

                    # split on commas and convert to vector
                    value = split(value, ",")
                    # check first value
                    if tryparse(Float64, value[1]) != nothing
                        # convert to numeric
                        value = tryparse.(Float64, value)
                        value = filter(x -> (x != nothing), value)
                        # check for integer conversion
                        if round.(value) == value
                            value = Int64.(value)
                        end
                    else
                        # convert to character and return
                        value = string.(strip.(value))
                    end

                    # remove arrays and convert to boolean if applicable
                    if length(value) == 1
                        value = value[1]
                        if lowercase(string(value)) in ["false", "true"]
                            value = (lowercase(string(value)) == "true")
                        end
                    end

                    dict_out[key] = value
                end

            end
        end

        return dict_out
    else
        print("\nFile $(fp_config) not found. Please check for the configuration file.")
        return(-1)
    end
end

##  FUNCTION FOR PRINTING NICE HEADERS
function print_header(str_in::String)
    nchar = length(str_in)
    char_div = "#"

    n_space = 4
    n_buffer = 3

    buffer = char_div^n_buffer
    spacer = " "^n_space

    bridge = char_div^(nchar + 2*(n_space + n_buffer))
    river = buffer*(" "^(nchar + 2*n_space))*buffer
    txt = buffer*spacer*str_in*spacer*buffer

    div = "\n"*join([bridge, river, txt, river, bridge], "\n")*"\n"

    print(div)
end

##  FUNCTION printnames USED TO PRINT LONG LISTS QUICKLY
function printnames(list::Array{String})
    for nm in list
        print("$(nm)\n")
    end
end

##  print_valid_values can be used to print a list of values into a string
function print_valid_values(array::Array{String, 1})
    str_valid = join("'".*(string.(array)).*"'", ", ", ", and ")
    return str_valid
end
function print_valid_values(array::Array{Int64, 1})
    str_valid = join("'".*(string.(array)).*"'", ", ", ", and ")
    return str_valid
end
function print_valid_values(array::Array{Symbol, 1})
    str_valid = join("'".*(string.(array)).*"'", ", ", ", and ")
    return str_valid
end

##  FUNCTION read_csv IS USED TO READ IN CSVs TO DATAFRAMES
function read_csv(path_in::String, header_q::Bool, lim::Int64=-1, names_vec::Array{Symbol, 1} = Array{Symbol, 1}(), delimiter::String = ",")

    if lim == -1
        #get some data
        df_in = CSV.File(
            path_in,
            delim = delimiter,
            ignorerepeated = false,
            header = header_q
        ) |> DataFrame
    else
        #get some data
        df_in = CSV.File(
            path_in,
            delim = delimiter,
            ignorerepeated = false,
            header = header_q;
            limit = lim
        ) |> DataFrame
    end

    if (!header_q & (length(names_vec) == size(df_in)[2]))
        dict_rnm = Dict([x for x in zip(Symbol.(names(df_in)), names_vec)])
        rename!(df_in, dict_rnm)
    end

    return df_in
end

##  FUNCTION FOR TRACKING TIME TO BUILD DB
function track_timer_db(checkpoint::Int64, time_0::Float64)
    print("\n Checkpoint (" * string(checkpoint)* ") complete at about " * string(round(time() - time_0)) * " seconds.\n\n")
end

# use try_parse_float to attempt to parse objects into a float; use multiple dispatch to avoid conditionals
function try_parse_float(x::String)
    tp = tryparse(Float64, x)
    if tp == nothing
        return missing
    else
        return tp
    end
end
function try_parse_float(x::Real)
    return Float64(x)
end
function try_parse_float(x::Missing)
    return missing
end

# use zip_cols to zip columns of a data frame or matrix together for iteration
function zip_cols(mat::DataFrame)
    return zip([mat[:, i] for i in 1:size(mat)[2]]...)
end
function zip_cols(mat::Matrix)
    return zip([mat[:, i] for i in 1:size(mat)[2]]...)
end

#############################################################
#    FUNCTIONS FOR SAMPLING AND OPTIMIZATION EXPERIMENTS    #
#############################################################

## FUNCTION get_rand_blocks SELECTS RANDOM BLOCKS (STRATIFIED BY POPULATION) TO USE TO ENFORCE CONSTRAINTS
function get_rand_blocks(vec::Array{Float64, 1}, q_vec::Array{Float64, 1} = [0, 0.25, 0.5, 0.75, 1], prob::Float64 = 0.1, min_num::Int64 = 5)

    # get the quantiles for "stratification"
    q = reverse(quantile(vec, sort(q_vec)))
    #add upper limit to ensure sampling is clean
    q[1] = q[1] + 1
    n_q = length(q) - 1

    # set the number of samples
    n = Int64(minimum([maximum([min_num, prob*length(vec)]), length(vec)]))

    n_per = Int64(floor(n/n_q))
    n_extra = n%n_q

    samp_out = Array{Int64, 1}(zeros(n))

    ind = 1

    for i in range(1, n_q, step = 1)
        if i <= n_extra
            n_samp = n_per + 1
        else
            n_samp = n_per
        end

        samp_space = findall(x -> (x < q[i]) & (x >= q[i + 1]), vec)
        samp_sub = rand_sample(samp_space, n_samp)

        samp_out[ind:(ind + n_samp - 1)] = samp_sub

        ind = ind + n_samp
    end

    return samp_out
end

##  FUNCTION rand_sample GENERATES A RANDOM SUBSAMPLE OF vec OF LENGTH n
function rand_sample(vec::AbstractArray, n::Int64)
    i = 1

    m = minimum([n, length(vec)])
    vec_out = Array{Any, 1}(zeros(m))

    while i <= m
        vec_out[i] = Random.rand(vec, 1)[1]
        vec = [x for x in vec if x != vec_out[i]]
        i += 1
    end

    return vec_out
end

##  sort row ids
function sort_row_ids(df_to_sort::DataFrame, fields_id::Array{Symbol, 1})

    # check id fields
    check_fields!(df_to_sort, fields_id)
    fields_in = Symbol.(names(df_to_sort))

    # cut data frame in two and sort fields, then re-concatenate
    df_sort_ids = copy(df_to_sort[:, fields_id])
    df_sort_ids = map(x -> sort(x), eachrow(Matrix(df_sort_ids)))
    df_sort_ids = DataFrame(permutedims(hcat(df_sort_ids...)), fields_id)

    df_to_sort = hcat(select(df_to_sort, Not(fields_id)), df_sort_ids)

    return df_to_sort[:, fields_in]

end


text:
@system BulbAppearance(Stage, FloralInitiation) begin
    bulb_appearable(floral_initiateable) ~ flag

    #HACK bulbing used to begin one phyllochron after floral initiation in bolting cultivars of garlic, see Meredith 2008
    bulb_appeared(floral_initiated) ~ flag

    bulb_appearing(bulb_appearable & !bulb_appeared) ~ flag

    # #FIXME postprocess similar to @produce?
    # def finish(self):
    #     GDD_sum = self.pheno.gdd_recorder.rate
    #     print(f"* Bulbing: time = {self.time}, GDDsum = {GDD_sum}")
end

ext: jl


text:
module TestInterface

using Phylo
using DataFrames
using Test

@testset "Build and tear down trees" begin
    @testset "For $TreeType" for TreeType in
        [NamedTree, NamedBinaryTree,
         BinaryTree{ManyRoots, DataFrame, Vector{Float64}},
         RootedTree, ManyRootTree]

        species = ["Dog", "Cat", "Human", "Potato", "Apple"]
        tree = TreeType(species)
        othernodes = ["Some 1", "Some 2"]
        nodes = createnodes!(tree, othernodes)
        @test [getnodename(tree, node) for node in nodes] == othernodes
        extra = [getnodename(tree, node) for node in createnodes!(tree, 1)]
        @test isa(extra[1], String)
        append!(othernodes, extra)
        also = getnodename(tree, createnode!(tree))
        @test also ∉ othernodes
        @test_throws Exception createnode!(tree, also)
        @test_throws Exception createnodes!(tree, [also])
        @test_throws Exception deletenode!(tree, "really not")
        @test deletenode!(tree, also)
        also = getnodename(tree, createnode!(tree))
        push!(othernodes, also)
        @test Set(othernodes) ⊆ Set(getnodenames(tree))
        innodes = append!(copy(species), othernodes)
        allnodes = reverse(innodes)
        pop!(innodes)
        @test_throws Exception getroot(tree)
        branches = map(innodes) do node
            itr = Iterators.filter(name -> hasoutboundspace(tree, name) &&
                                   name != node, allnodes)
            getbranch(tree, createbranch!(tree, first(itr), node))
        end
        @test_nowarn getroot(tree)
        branchnames = [getbranchname(tree, branch) for branch in branches]
        @test Set(branchnames) == Set(getbranchnames(tree))
        @test validate!(tree)
        @test_throws Exception createbranch!(tree, allnodes[1], allnodes[2])
        who = getparent(tree, species[1])
        b = getbranch(tree, getinbound(tree, species[1]))
        bn = getbranchname(tree, b)
        node1 = getnode(tree, species[1])
        @test !isunattached(tree, node1) &&
            hasoutboundspace(tree, node1) && !hasinboundspace(tree, node1) &&
            outdegree(tree, node1) == 0 && indegree(tree, node1) == 1
        @test hasbranch(tree, b)
        @test hasbranch(tree, bn)
        @test getnodename(tree, dst(tree, getbranch(tree, b))) == species[1]
        @test src(tree, getbranch(tree, b)) == getnode(tree, who)
        @test deletebranch!(tree, b)
        @test !hasbranch(tree, b)
        @test isunattached(tree, species[1]) &&
            hasoutboundspace(tree, species[1]) &&
            hasinboundspace(tree, species[1]) &&
            outdegree(tree, species[1]) == 0 &&
            indegree(tree, species[1]) == 0
        @test_throws Exception getbranch(tree, b)
        branches = [branch for branch in branches if branch != b]
        @test Set(branches) == Set(getbranches(tree))
        b3 = getinbound(tree, species[2])
        source = src(tree, b3)
        destination = dst(tree, b3)
        @test deletebranch!(tree, b3)
        branches = [branch for branch in branches if branch != b3]
        b3 = createbranch!(tree, who, species[1])
        @test who == getnodename(tree, src(tree, b3))
        spparent = getparent(tree, getnode(tree, species[1]))
        @test getnode(tree, spparent) == getnode(tree, src(tree, b3))
        @test species[1] == getnodename(tree, dst(tree, b3))
        b2 = getbranchname(tree, createbranch!(tree, source, destination))
        branchnames = [getbranchname(tree, branch) for branch in branches]
        push!(branchnames, b2, getbranchname(tree, b3))
        @test Set(branchnames) == Set(getbranchnames(tree))
        @test deletenode!(tree, species[1])
        @test !hasnode(tree, species[1])
        @test species[1] == getnodename(tree, createnode!(tree, species[1]))
        @test_nowarn createbranch!(tree, who, species[1])
        @test hasnode(tree, species[1])
        @test validate!(tree)
        @test all(isleaf(tree, node) for node in species)
        @test all((!isroot(tree, node) & !isunattached(tree, node) &
                   !isinternal(tree, node)) for node in species)
        @test Set(getnodename(tree, src(tree, branch))
                  for branch in getbranches(tree)) ∪
            Set(getnodename(tree, dst(tree, branch))
                for branch in getbranches(tree)) == Set(getnodenames(tree))
        createnode!(tree)
        @test (roottype(TreeType) == OneRoot) ⊻ validate!(tree)
    end
end

end


text:
module BrainForth

using MacroTools

include("brainfuck.jl")
include("bootstrap.jl")
include("forth.jl")
include("array.jl")
include("text.jl")
include("utils.jl")

end # module

max_stars_repo_name: JuliaTagBot/BrainForth

max_issues_repo_name: JuliaTagBot/BrainForth

max_forks_repo_name: JuliaTagBot/BrainForth


text:
# Autogenerated wrapper script for Spot_julia_jll for x86_64-linux-gnu-cxx03
export autcross, autfilt, dstar2tgba, genaut, genltl, libbddx, libspot, libspot_julia, libspotgen, libspotltsmin, ltl2tgba, ltl2tgta, ltlcross, ltldo, ltlfilt, ltlgrind, ltlsynt, randaut, randltl

using libcxxwrap_julia_jll
JLLWrappers.@generate_wrapper_header("Spot_julia")
JLLWrappers.@declare_executable_product(autcross)
JLLWrappers.@declare_executable_product(autfilt)
JLLWrappers.@declare_executable_product(dstar2tgba)
JLLWrappers.@declare_executable_product(genaut)
JLLWrappers.@declare_executable_product(genltl)
JLLWrappers.@declare_library_product(libbddx, "libbddx.so.0")
JLLWrappers.@declare_library_product(libspot, "libspot.so.0")
JLLWrappers.@declare_library_product(libspot_julia, "libspot_julia.so")
JLLWrappers.@declare_library_product(libspotgen, "libspotgen.so.0")
JLLWrappers.@declare_library_product(libspotltsmin, "libspotltsmin.so.0")
JLLWrappers.@declare_executable_product(ltl2tgba)
JLLWrappers.@declare_executable_product(ltl2tgta)
JLLWrappers.@declare_executable_product(ltlcross)
JLLWrappers.@declare_executable_product(ltldo)
JLLWrappers.@declare_executable_product(ltlfilt)
JLLWrappers.@declare_executable_product(ltlgrind)
JLLWrappers.@declare_executable_product(ltlsynt)
JLLWrappers.@declare_executable_product(randaut)
JLLWrappers.@declare_executable_product(randltl)
function __init__()
    JLLWrappers.@generate_init_header(libcxxwrap_julia_jll)
    JLLWrappers.@init_executable_product(
        autcross,
        "bin/autcross",
    )

    JLLWrappers.@init_executable_product(
        autfilt,
        "bin/autfilt",
    )

    JLLWrappers.@init_executable_product(
        dstar2tgba,
        "bin/dstar2tgba",
    )

    JLLWrappers.@init_executable_product(
        genaut,
        "bin/genaut",
    )

    JLLWrappers.@init_executable_product(
        genltl,
        "bin/genltl",
    )

    JLLWrappers.@init_library_product(
        libbddx,
        "lib/libbddx.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libspot,
        "lib/libspot.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libspot_julia,
        "lib/libspot_julia.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libspotgen,
        "lib/libspotgen.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libspotltsmin,
        "lib/libspotltsmin.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        ltl2tgba,
        "bin/ltl2tgba",
    )

    JLLWrappers.@init_executable_product(
        ltl2tgta,
        "bin/ltl2tgta",
    )

    JLLWrappers.@init_executable_product(
        ltlcross,
        "bin/ltlcross",
    )

    JLLWrappers.@init_executable_product(
        ltldo,
        "bin/ltldo",
    )

    JLLWrappers.@init_executable_product(
        ltlfilt,
        "bin/ltlfilt",
    )

    JLLWrappers.@init_executable_product(
        ltlgrind,
        "bin/ltlgrind",
    )

    JLLWrappers.@init_executable_product(
        ltlsynt,
        "bin/ltlsynt",
    )

    JLLWrappers.@init_executable_product(
        randaut,
        "bin/randaut",
    )

    JLLWrappers.@init_executable_product(
        randltl,
        "bin/randltl",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using DiscreteChoice
using JLD2
using ForwardDiff
using Optim
using StatsBase

@load "example/data/homo.jld2" Y X β ξ

P = zeros(10)
[P[i] = mean(Y.==i) for i in 1:10]
sum(P)

function y_prob(X, β, ξ)
    ξ = vcat(zero(eltype(ξ)), ξ)
    index = exp.(X*β .+ ξ)
    y_prob = index ./ sum(index)
    return y_prob
end

model_1 = DCModel(Y, X; hetero_preference=false)
model_2 = DCModel(Y, X; hetero_preference=false)

# Start from true value
estDCModel!(model_1; init_guess = vcat(β, ξ))
abs.(Optim.minimizer(model_1.optResults) .- vcat(β, ξ)) |> maximum
updateCoef!(model_1)

P .- y_prob(model_1.X, model_1.β, model_1.ξ)

# Start from random value
estDCModel!(model_2; init_guess = randn(length(β)+length(ξ)))
updateCoef!(model_2)

y_emp = y_prob(model_2.X, model_2.β, model_2.ξ)
P .- y_emp
# Notice that this can perfectly match the data. Why?
# simple reason is that the number of products is too small !!
# Because there are only 10 products, essentially only 10 data points, too many dimensions
# This will get changed if we include individual demographics into the data

abs.(Optim.minimizer(model_2.optResults) .- vcat(β, ξ)) |> maximum
Optim.minimizer(model_2.optResults) .- vcat(β, ξ)

β
ξ
# Why does the random starting value performs poorly? Is this some non-identification issue?

res_homo = opt_homo(Y, X, vcat(β, ξ))


text:
import LinearAlgebra: checksquare
import LinearAlgebra.BLAS: BlasFloat, BlasInt, liblapack, @blasfunc

struct LUSolverLAPACK{T<:BlasFloat} <: LinearSolver{T}
    n::BlasInt
    A::Matrix{T}
    pivots::Vector{BlasInt}
    info::BlasInt
end

function LUSolverLAPACK(A::Matrix{T}) where {T}
    n = checksquare(A)
    lu = LUSolverLAPACK{T}(n, zero(A), zeros(BlasInt, n), 0)
    factorize!(lu, A)
end

LUSolverLAPACK{T}(n::BlasInt) where {T} = LUSolverLAPACK(zeros(T, n, n))

## LAPACK LU factorization and solver for general matrices (GE)
for (getrf, getrs, elty) in
    ((:dgetrf_,:dgetrs_,:Float64),
     (:sgetrf_,:sgetrs_,:Float32),
     (:zgetrf_,:zgetrs_,:ComplexF64),
     (:cgetrf_,:cgetrs_,:ComplexF32))
    @eval begin
        function factorize!(lu::LUSolverLAPACK{$elty}, A::AbstractMatrix{$elty})
            copy!(lu.A, A)
            ccall((@blasfunc($getrf), liblapack), Nothing,
                  (Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty},
                   Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}),
                   Ref(lu.n), Ref(lu.n), lu.A, Ref(lu.n), lu.pivots, Ref(lu.info))

            if lu.info > 0
                throw(SingularException(lu.info))
            elseif lu.info < 0
                throw(ArgumentError(lu.info))
            end
            return lu
        end

        function LinearAlgebra.ldiv!(x::AbstractVector{$elty}, lu::LUSolverLAPACK{$elty}, b::AbstractVector{$elty})
            copy!(x, b)
            trans = UInt8('N')
            nrhs = BlasInt(1)
            ccall((@blasfunc($getrs), liblapack), Nothing,
                  (Ptr{UInt8}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt},
                   Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}),
                   Ref(trans), Ref(lu.n), Ref(nrhs), lu.A, Ref(lu.n), lu.pivots, x, Ref(lu.n), Ref(lu.info))

            if lu.info < 0
                throw(ArgumentError(lu.info))
            end
            return x
        end
    end
end


text:
"""
GenX: An Configurable Capacity Expansion Model
Copyright (C) 2021,  Massachusetts Institute of Technology
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
A complete copy of the GNU General Public License v2 (GPLv2) is available
in LICENSE.txt.  Users uncompressing this from an archive may not have
received this license file.  If not, see <http://www.gnu.org/licenses/>.
"""

@doc raw"""
	load_minimum_capacity_requirement(path::AbstractString, inputs::Dict, setup::Dict)

Function for reading input parameters related to mimimum capacity requirement constraints (e.g. technology specific deployment mandates)
"""
function load_minimum_capacity_requirement(path::AbstractString, inputs::Dict, setup::Dict)
	MinCapReq = DataFrame(CSV.File(joinpath(path, "Minimum_capacity_requirement.csv"), header=true), copycols=true)
	NumberOfMinCapReqs = size(collect(skipmissing(MinCapReq[!,:MinCapReqConstraint])),1)
	inputs["NumberOfMinCapReqs"] = NumberOfMinCapReqs
	inputs["MinCapReq"] = MinCapReq[!,:Min_MW]

	if setup["ParameterScale"] == 1
		inputs["MinCapReq"] = inputs["MinCapReq"]/ModelScalingFactor # Convert to GW
	end

	println("Minimum_capacity_requirement.csv Successfully Read!")

	return inputs
end

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
"""
    PODeigen(X)

Uses the eigenvalue method of snapshots to calculate the POD basis of X. Method of
snapshots is efficient when number of data points `n` > number of snapshots `m`.
"""
function PODeigen(X; subtractmean::Bool = false)

    Xcop = deepcopy(X)
    PODeigen!(Xcop, subtractmean = subtractmean)

end

"""
    PODeigen(X,W)

Same as `PODeigen(X)` but uses weights for each data point. The weights are equal to the
cell volume for a volume mesh.
"""
function PODeigen(X,W::AbstractVector; subtractmean::Bool = false)

    Xcop = deepcopy(X)
    PODeigen!(Xcop,W, subtractmean = subtractmean)

end

"""
    PODeigen!(X,W)

Same as `PODeigen!(X)` but uses weights for each data point. The weights are equal to the
cell volume for a volume mesh.
"""
function PODeigen!(X,W::AbstractVector; subtractmean::Bool = false)

    if subtractmean
        X .-= mean(X,dims=2)
    end

    # Number of snapshots
    m = size(X,2)

    # Correlation matrix for method of snapshots
    C = X'*Diagonal(W)*X

    # Eigen Decomposition
    E = eigen!(C)
    eigVals = E.values
    eigVects = E.vectors

    # Sort the eigen vectors
    sortInd = sortperm(abs.(eigVals)/m,rev=true)
    eigVects = eigVects[:,sortInd]
    eigVals = eigVals[sortInd]

    # Diagonal matrix containing the square roots of the eigenvalues
    S = sqrt.(abs.(eigVals))

    # Construct the modes and coefficients
    Φ = X*eigVects*Diagonal(1 ./S)
    a = Diagonal(S)*eigVects'

    POD = PODBasis(a, Φ)

    return POD, S

end

"""
    PODeigen!(X)

Same as `PODeigen(X)` but overwrites memory.
"""
function PODeigen!(X; subtractmean::Bool = false)

    if subtractmean
        X .-= mean(X,dims=2)
    end

    # Number of snapshots
    m = size(X,2)

    # Correlation matrix for method of snapshots
    C = X'*X

    # Eigen Decomposition
    E = eigen!(C)
    eigVals = E.values
    eigVects = E.vectors

    # Sort the eigen vectors
    sortInd = sortperm(abs.(eigVals)/m,rev=true)
    eigVects = eigVects[:,sortInd]
    eigVals = eigVals[sortInd]

    # Diagonal matrix containing the square roots of the eigenvalues
    S = sqrt.(abs.(eigVals))

    # Construct the modes and coefficients
    Φ = X*eigVects*Diagonal(1 ./S)
    a = Diagonal(S)*eigVects'

    POD = PODBasis(a, Φ)

    return POD, S

end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct VirtualMachineScaleSetDataDisk <: SwaggerModel
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    lun::Any # spec type: Union{ Nothing, Int32 } # spec name: lun
    caching::Any # spec type: Union{ Nothing, Caching } # spec name: caching
    writeAcceleratorEnabled::Any # spec type: Union{ Nothing, Bool } # spec name: writeAcceleratorEnabled
    createOption::Any # spec type: Union{ Nothing, CreateOption } # spec name: createOption
    diskSizeGB::Any # spec type: Union{ Nothing, Int32 } # spec name: diskSizeGB
    managedDisk::Any # spec type: Union{ Nothing, VirtualMachineScaleSetManagedDiskParameters } # spec name: managedDisk
    diskIOPSReadWrite::Any # spec type: Union{ Nothing, Int64 } # spec name: diskIOPSReadWrite
    diskMBpsReadWrite::Any # spec type: Union{ Nothing, Int64 } # spec name: diskMBpsReadWrite

    function VirtualMachineScaleSetDataDisk(;name=nothing, lun=nothing, caching=nothing, writeAcceleratorEnabled=nothing, createOption=nothing, diskSizeGB=nothing, managedDisk=nothing, diskIOPSReadWrite=nothing, diskMBpsReadWrite=nothing)
        o = new()
        validate_property(VirtualMachineScaleSetDataDisk, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(VirtualMachineScaleSetDataDisk, Symbol("lun"), lun)
        setfield!(o, Symbol("lun"), lun)
        validate_property(VirtualMachineScaleSetDataDisk, Symbol("caching"), caching)
        setfield!(o, Symbol("caching"), caching)
        validate_property(VirtualMachineScaleSetDataDisk, Symbol("writeAcceleratorEnabled"), writeAcceleratorEnabled)
        setfield!(o, Symbol("writeAcceleratorEnabled"), writeAcceleratorEnabled)
        validate_property(VirtualMachineScaleSetDataDisk, Symbol("createOption"), createOption)
        setfield!(o, Symbol("createOption"), createOption)
        validate_property(VirtualMachineScaleSetDataDisk, Symbol("diskSizeGB"), diskSizeGB)
        setfield!(o, Symbol("diskSizeGB"), diskSizeGB)
        validate_property(VirtualMachineScaleSetDataDisk, Symbol("managedDisk"), managedDisk)
        setfield!(o, Symbol("managedDisk"), managedDisk)
        validate_property(VirtualMachineScaleSetDataDisk, Symbol("diskIOPSReadWrite"), diskIOPSReadWrite)
        setfield!(o, Symbol("diskIOPSReadWrite"), diskIOPSReadWrite)
        validate_property(VirtualMachineScaleSetDataDisk, Symbol("diskMBpsReadWrite"), diskMBpsReadWrite)
        setfield!(o, Symbol("diskMBpsReadWrite"), diskMBpsReadWrite)
        o
    end
end # type VirtualMachineScaleSetDataDisk

const _property_map_VirtualMachineScaleSetDataDisk = Dict{Symbol,Symbol}(Symbol("name")=>Symbol("name"), Symbol("lun")=>Symbol("lun"), Symbol("caching")=>Symbol("caching"), Symbol("writeAcceleratorEnabled")=>Symbol("writeAcceleratorEnabled"), Symbol("createOption")=>Symbol("createOption"), Symbol("diskSizeGB")=>Symbol("diskSizeGB"), Symbol("managedDisk")=>Symbol("managedDisk"), Symbol("diskIOPSReadWrite")=>Symbol("diskIOPSReadWrite"), Symbol("diskMBpsReadWrite")=>Symbol("diskMBpsReadWrite"))
const _property_types_VirtualMachineScaleSetDataDisk = Dict{Symbol,String}(Symbol("name")=>"String", Symbol("lun")=>"Int32", Symbol("caching")=>"Caching", Symbol("writeAcceleratorEnabled")=>"Bool", Symbol("createOption")=>"CreateOption", Symbol("diskSizeGB")=>"Int32", Symbol("managedDisk")=>"VirtualMachineScaleSetManagedDiskParameters", Symbol("diskIOPSReadWrite")=>"Int64", Symbol("diskMBpsReadWrite")=>"Int64")
Base.propertynames(::Type{ VirtualMachineScaleSetDataDisk }) = collect(keys(_property_map_VirtualMachineScaleSetDataDisk))
Swagger.property_type(::Type{ VirtualMachineScaleSetDataDisk }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_VirtualMachineScaleSetDataDisk[name]))}
Swagger.field_name(::Type{ VirtualMachineScaleSetDataDisk }, property_name::Symbol) =  _property_map_VirtualMachineScaleSetDataDisk[property_name]

function check_required(o::VirtualMachineScaleSetDataDisk)
    (getproperty(o, Symbol("lun")) === nothing) && (return false)
    (getproperty(o, Symbol("createOption")) === nothing) && (return false)
    true
end

function validate_property(::Type{ VirtualMachineScaleSetDataDisk }, name::Symbol, val)
end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.8.0-DEV
#     language: julia
#     name: julia-1.8
# ---

# %%
module My

struct Hello{T}
    a::T
end

end

My.Hello("Julia")

# %%
module My

struct Hello{T}
    a::T
end

function Base.show(io::IO, x::Hello)
    print(io, "Hello, ")
    show(io, x.a)
    print(io, '!')
end

end

My.Hello("Julia")

# %%
dump(My.Hello("Julia"))

# %%
str = repr(My.Hello("Julia"))

# %%
π

# %%
My.Hello(π)

# %%
module My

struct Hello{T}
    a::T
end

function Base.show(io::IO, x::Hello)
    print(io, "Hello, ")
    show(io, x.a)
    print(io, '!')
end

function Base.show(io::IO, ::MIME"text/plain", x::Hello)
    print(io, "Hello, ")
    show(io, MIME"text/plain"(), x.a)
    print(io, "!")
end

end

My.Hello("Julia")

# %%
My.Hello(π)

# %%
methods(Base.show)

# %%
methods(repr)

# %%
methods(string)

# %%
methods(print)

# %%

max_stars_repo_path: 0015/Base.show example.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0015/Base.show example.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0015/Base.show example.jl

max_forks_repo_name: genkuroki/public


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

# tests for accurate updating of method tables

using Base: get_world_counter
tls_world_age() = ccall(:jl_get_tls_world_age, UInt, ())
@test typemax(UInt) > get_world_counter() == tls_world_age() > 0

# test simple method replacement
begin
    g265a() = f265a(0)
    f265a(x::Any) = 1
    @test g265a() == 1
    @test Base.return_types(g265a, ()) == Any[Int]
    @test Core.Compiler.return_type(g265a, ()) == Int

    f265a(x::Any) = 2.0
    @test g265a() == 2.0

    @test Base.return_types(g265a, ()) == Any[Float64]
    @test Core.Compiler.return_type(g265a, ()) == Float64
end

# test signature widening
begin
    f265b(x::Int) = 1
    let ty = Any[1, 2.0e0]
        global g265b(i::Int) = f265b(ty[i])
    end
    @test g265b(1) == 1
    @test Base.return_types(g265b, (Int,)) == Any[Int]
    @test Core.Compiler.return_type(g265b, (Int,)) == Int

    f265b(x::Any) = 2.0
    @test g265b(1) == 1
    @test g265b(2) == 2.0
    @test Base.return_types(g265b, (Int,)) == Any[Union{Int, Float64}]
    @test Core.Compiler.return_type(g265b, (Int,)) == Union{Int, Float64}
end

# test signature narrowing
begin
    g265c() = f265c(0)
    f265c(x::Any) = 1
    @test g265c() == 1
    @test Base.return_types(g265c, ()) == Any[Int]
    @test Core.Compiler.return_type(g265c, ()) == Int

    f265c(x::Int) = 2.0
    @test g265c() == 2.0

    @test Base.return_types(g265c, ()) == Any[Float64]
    @test Core.Compiler.return_type(g265c, ()) == Float64
end

# test constructor narrowing
mutable struct A265{T}
    field1::T
end
A265_() = A265(1)
@test (A265_()::A265{Int}).field1 === 1
A265(fld::Int) = A265(Float64(fld))
@test (A265_()::A265{Float64}).field1 === 1.0e0

# test constructor widening
mutable struct B265{T}
    field1::T
    # dummy arg is present to prevent (::Type{T}){T}(arg) from matching the test calls
    B265{T}(field1::Any, dummy::Nothing) where {T} = new(field1) # prevent generation of outer ctor
end
  # define some constructors
B265(x::Int, dummy::Nothing) = B265{Int}(x, dummy)
let ty = Any[1, 2.0e0, 3.0f0]
    global B265_(i::Int) = B265(ty[i], nothing)
end
  # test for correct answers
@test (B265_(1)::B265{Int}).field1 === 1
@test_throws MethodError B265_(2)
@test_throws MethodError B265_(3)
@test Base.return_types(B265_, (Int,)) == Any[B265{Int}]
@test Core.Compiler.return_type(B265_, (Int,)) == B265{Int}

  # add new constructors
B265(x::Float64, dummy::Nothing) = B265{Float64}(x, dummy)
B265(x::Any, dummy::Nothing) = B265{UInt8}(x, dummy)

  # make sure answers are updated
@test (B265_(1)::B265{Int}).field1 === 1
@test (B265_(2)::B265{Float64}).field1 === 2.0e0
@test (B265_(3)::B265{UInt8}).field1 === 0x03

@test Base.return_types(B265_, (Int,)) == Any[B265]
@test Core.Compiler.return_type(B265_, (Int,)) == B265

# test oldworld call / inference
function wfunc(c1,c2)
    while true
        (f, args) = take!(c1)
        put!(c2, f(args...))
    end
end
function put_n_take!(v...)
    put!(chnls[1], v)
    take!(chnls[2])
end

g265() = [f265(x) for x in 1:3.]
wc265 = get_world_counter()
f265(::Any) = 1.0
@test wc265 + 1 == get_world_counter()
chnls, tasks = Base.channeled_tasks(2, wfunc)
t265 = tasks[1]

wc265 = get_world_counter()
@test put_n_take!(get_world_counter, ()) == wc265
@test put_n_take!(tls_world_age, ()) == wc265
f265(::Int) = 1
@test put_n_take!(get_world_counter, ()) == wc265 + 1 == get_world_counter() == tls_world_age()
@test put_n_take!(tls_world_age, ()) == wc265

@test g265() == Int[1, 1, 1]
@test Core.Compiler.return_type(f265, (Any,)) == Union{Float64, Int}
@test Core.Compiler.return_type(f265, (Int,)) == Int
@test Core.Compiler.return_type(f265, (Float64,)) == Float64

@test put_n_take!(g265, ()) == Float64[1.0, 1.0, 1.0]
@test put_n_take!(Core.Compiler.return_type, (f265, (Any,))) == Float64
@test put_n_take!(Core.Compiler.return_type, (f265, (Int,))) == Float64
@test put_n_take!(Core.Compiler.return_type, (f265, (Float64,))) == Float64
@test put_n_take!(Core.Compiler.return_type, (f265, (Float64,))) == Float64

# test that reflection ignores worlds
@test Base.return_types(f265, (Any,)) == Any[Int, Float64]
@test put_n_take!(Base.return_types, (f265, (Any,))) == Any[Int, Float64]

# test for method errors
h265() = true
loc_h265 = "$(@__FILE__):$(@__LINE__() - 1)"
@test h265()
@test_throws MethodError put_n_take!(h265, ())
@test_throws MethodError fetch(t265)
@test istaskdone(t265)
let ex = t265.exception
    @test ex.f == h265
    @test ex.args == ()
    @test ex.world == wc265
    str = sprint(showerror, ex)
    wc = get_world_counter()
    cmps = """
        MethodError: no method matching h265()
        The applicable method may be too new: running in world age $wc265, while current world is $wc."""
    @test startswith(str, cmps)
    cmps = "\n  h265() at $loc_h265 (method too new to be called from this world context.)"
    @test occursin(cmps, str)
end

# test for generated function correctness
# and min/max world computation validity of cache_method
f_gen265(x) = 1
@generated g_gen265(x) = f_gen265(x)
@generated h_gen265(x) = :(f_gen265(x))
f_gen265(x::Int) = 2
f_gen265(x::Type{Int}) = 3
@generated g_gen265b(x) = f_gen265(x)
@test h_gen265(0) == 2
@test g_gen265(0) == 1
@test f_gen265(Int) == 3
@test g_gen265b(0) == 3

# Test that old, invalidated specializations don't get revived for
# intermediate worlds by later additions to the method table that
# would have capped those specializations if they were still valid
f26506(@nospecialize(x)) = 1
g26506(x) = f26506(x[1])
z = Any["ABC"]
f26506(x::Int) = 2
g26506(z) # Places an entry for f26506(::String) in mt.name.cache
f26506(x::String) = 3
let cache = typeof(f26506).name.mt.cache
    # The entry we created above should have been truncated
    @test cache.min_world == cache.max_world
end
c26506_1, c26506_2 = Condition(), Condition()
# Captures the world age
result26506 = Any[]
t = Task(()->begin
    wait(c26506_1)
    push!(result26506, g26506(z))
    notify(c26506_2)
end)
yield(t)
f26506(x::Float64) = 4
let cache = typeof(f26506).name.mt.cache
    # The entry we created above should have been truncated
    @test cache.min_world == cache.max_world
end
notify(c26506_1)
wait(c26506_2)
@test result26506[1] == 3

max_stars_repo_name: musm/julia

max_issues_repo_name: AsForMe/julia

max_forks_repo_name: AsForMe/julia


text:
using Test: @testset, @test, @test_throws, @test_broken
using CBindingGen
using CBindingGen.Clang
using CBinding

@testset "CBindingGen" begin
	@testset "parsing C" begin
		ctx = ConverterContext(decl -> false, nothing)
		parse_headers!(ctx, ["test.h"], args = ["-I."])
		@test isempty(ctx.converted)

		ctx = ConverterContext(decl -> true, nothing)
		parse_headers!(ctx, ["test.h"], args = ["-I."])
		@test !isempty(ctx.converted)
	end

	ctx = ConverterContext(decl -> true, nothing)
	parse_headers!(ctx, ["test.h"], args = ["-I."])
	generate(ctx, pwd(), ".")

	mod = @eval module $(gensym())
		import ..CBinding

		CBinding.@cstruct _S24
		CBinding.@cstruct _S31

		include(joinpath(pwd(), ".-atdevelop.jl"))
		include(joinpath(pwd(), ".-atcompile.jl"))
		include(joinpath(pwd(), ".-atcompile_typedefs.jl"))
		include(joinpath(pwd(), ".-atcompile_bindings.jl"))
	end

	@testset "exports" begin
		@test mod isa Module

		expected = (
			(Symbol("T", i) for i in 1:26)...,       # typedefs
			(Symbol("S", i) for i in 1:32)...,       # structs
			(Symbol("E", i) for i in 1:5)...,        # enums
			(Symbol("E", i, "_1") for i in 1:5)...,  # enum values
			(Symbol("g", i) for i in 1:7)...,        # globals
			(Symbol("f", i) for i in 1:10)...,       # functions
			:E1_2, :E1_3, :E2_2, :E2_3, :E2_4, :E2_5, :E2_6, :E4_,
			:S27_, :S29_, :S29_ptr,
			:F1,
		)
		present = names(mod)
		for p in present
			@test p in expected
		end
		for e in expected
			@test e in present
		end
	end

	@testset "typedefs" begin
		@test mod.T1 === Cvoid
		@test mod.T2 === mod.T1
		@test mod.T3 === Cchar
		@test mod.T4 === Cshort
		@test mod.T5 === Cint
		@test mod.T6 === Clong
		@test mod.T7 === Clonglong
		@test mod.T8 === Cchar
		@test mod.T9 === Cshort
		@test mod.T10 === Cint
		@test mod.T11 === Clong
		@test mod.T12 === Clonglong
		@test mod.T13 === Cuchar
		@test mod.T14 === Cushort
		@test mod.T15 === Cuint
		@test mod.T16 === Culong
		@test mod.T17 === Culonglong
		@test mod.T18 === Cfloat
		@test mod.T19 === Cdouble
		@test mod.T20 === Clongdouble
		@test mod.T21 === Complex{Cfloat}
		@test mod.T22 === Complex{Cdouble}
		@test mod.T23 === Complex{Clongdouble}
		@test mod.T24 === @carray Cint[2]
		@test mod.T25 === @carray (@carray Cint[2])[5]
	end

	@testset "structs" begin
		@test mod.S1 <: Cstruct
		@test () == propertynames(mod.S1)

		@test mod.S2 <: Cstruct
		@test (:i,) == propertynames(mod.S2)
		@test (Cint,) == propertytypes(mod.S2)

		@test mod.S3 <: Cstruct
		@test (:s,) == propertynames(mod.S3)
		@test (mod.S1,) == propertytypes(mod.S3)

		@test mod.S4 <: Cstruct
		@test (:i, :s) == propertynames(mod.S4)
		@test (Cint, mod.S3) == propertytypes(mod.S4)

		@test mod.S5 <: Cstruct
		@test (:s,) == propertynames(mod.S5)

		@test mod.S6 <: Cstruct
		@test (:s, :s1, :i, :s2) == propertynames(mod.S6)

		@test mod.S7 <: Cstruct
		@test (:i, :s) == propertynames(mod.S7)

		@test mod.S8 <: Cstruct
		@test (:s,) == propertynames(mod.S8)

		@test mod.S9 <: Cstruct
		@test (:s1, :s2) == propertynames(mod.S9)

		@test mod.S10 <: Cstruct
		@test () == propertynames(mod.S10)

		@test mod.S11 <: Cstruct
		@test (:i,) == propertynames(mod.S11)
		@test (Cint,) == propertytypes(mod.S2)

		@test mod.S12 <: Cstruct
		@test (:s1, :s2, :s3, :s4) == propertynames(mod.S12)
		typ = propertytypes(mod.S12)[1]
		@test (typ, Ptr{typ}, Ptr{Ptr{typ}}, Ptr{Ptr{Ptr{typ}}}) == propertytypes(mod.S12)

		@test mod.S13 <: Cstruct
		@test (:s,) == propertynames(mod.S13)
		@test propertytypes(mod.S13)[1] <: Ptr

		@test mod.S14 <: Cstruct
		@test (:s1, :s2, :s3, :s4, :s5, :s6) == propertynames(mod.S14)
		typ = propertytypes(mod.S14)[1]
		@test (typ, Ptr{typ}, (@carray typ[2]), (@carray (@carray typ[2])[5]), (@carray (@carray Ptr{typ}[2])[5]), Ptr{(@carray Ptr{typ}[2])}) == propertytypes(mod.S14)

		@test mod.S15 <: Cstruct
		@test (:s,) == propertynames(mod.S15)
		@test eltype(propertytypes(mod.S15)[1]) <: Cstruct

		@test mod.S16 <: Cstruct
		@test (:i, :u) == propertynames(mod.S16)
		@test (Cint, Cuint) == propertytypes(mod.S16)

		@test mod.S17 <: Cstruct
		@test (:c, :i) == propertynames(mod.S17)
		@test (Cchar, Cint) == propertytypes(mod.S17)
		@test CBinding._strategy(mod.S17) === CBinding.ALIGN_NATIVE

		@test mod.S18 <: Cstruct
		@test (:c, :i) == propertynames(mod.S18)
		@test (Cchar, Cint) == propertytypes(mod.S18)
		@test CBinding._strategy(mod.S18) === CBinding.ALIGN_PACKED

		@test mod.S19 <: Cstruct
		@test (:c, :i) == propertynames(mod.S19)
		@test (Cchar, Cint) == propertytypes(mod.S19)
		@test CBinding._strategy(mod.S19) === CBinding.ALIGN_NATIVE

		@test mod.S20 <: Cstruct
		@test (:c, :i) == propertynames(mod.S20)
		@test (Cchar, Cint) == propertytypes(mod.S20)
		@test CBinding._strategy(mod.S20) === CBinding.ALIGN_PACKED

		@test mod.S21 <: Cstruct
		@test (:f,) == propertynames(mod.S21)
		@test (Ptr{Cfunction{Cvoid, Tuple{}}},) == propertytypes(mod.S21)

		@test mod.S22 <: Cstruct
		@test (:f,) == propertynames(mod.S22)
		@test (Ptr{Cfunction{Ptr{Cfunction{Cint, Tuple{Cint, Cfloat}}}, Tuple{Ptr{mod.S21}, Cint}}},) == propertytypes(mod.S22)

		@test mod.S23 <: Cstruct
		@test (:f,) == propertynames(mod.S23)
		@test (Ptr{Cfunction{Cint, Tuple{Ptr{Cchar}, Vararg}}},) == propertytypes(mod.S23)

		@test mod.S25 <: Cstruct
		@test (:s1, :s2, :s3) == propertynames(mod.S25)
		@test (Ptr{mod.S24}, Ptr{mod._S24}, Ptr{mod.S25}) == propertytypes(mod.S25)

		@test mod.S24 <: Cstruct
		@test (:s1, :s2, :s3) == propertynames(mod.S24)
		@test (Ptr{mod.S24}, Ptr{mod._S24}, mod.S25) == propertytypes(mod.S24)

		@test mod.S26 <: Cstruct
		@test (:e1, :e2, :e3, :e4, :e5) == propertynames(mod.S26)
		@test (mod.E3, Ptr{mod.E3}, mod.E3, mod.E1, Ptr{mod.E1}) == propertytypes(mod.S26)

		@test mod.S27 <: Cstruct
		@test (:i,) == propertynames(mod.S27)
		@test (Cint,) == propertytypes(mod.S27)

		@test mod.S28 <: Cstruct
		@test (:i,) == propertynames(mod.S28)
		@test (Cint,) == propertytypes(mod.S28)

		@test mod.S29 <: Cstruct
		@test (:s1, :s2) == propertynames(mod.S29)
		@test (Ptr{mod.S29}, mod.S29_ptr) == propertytypes(mod.S29)

		@test mod.S30 <: Cstruct
		@test (:f,) == propertynames(mod.S30)
		@test (Ptr{Cfunction{Cvoid, Tuple{mod.S30}}},) == propertytypes(mod.S30)

		@test mod.S32 <: Cstruct
		@test (:s,) == propertynames(mod.S32)
		@test (Ptr{mod.S31},) == propertytypes(mod.S32)
	end

	@testset "enums" begin
		@test mod.E1 <: Cenum
		@test mod.E1_1 == 0
		@test mod.E1_2 == 1
		@test mod.E1_3 == 2

		@test mod.E2 <: Cenum
		@test mod.E2_1 == 1
		@test mod.E2_2 == mod.E2_1
		@test mod.E2_3 == mod.E2_1+10
		@test mod.E2_4 == 0
		@test mod.E2_5 == 0xffffffff
		@test mod.E2_6 == 0xffffffffffff

		@test mod.E3 <: Cenum
		@test mod.E3_1 == 0

		@test mod.E4 <: Cenum
		@test mod.E4_1 == 0

		@test mod.E5 <: Cenum
		@test mod.E5_1 == 0
	end

	@testset "globals" begin
		@test mod.g1 isa Cglobalconst{Cint}
		@test mod.g2 isa Cglobalconst{mod.S24}
		@test mod.g3 isa Cglobalconst{Ptr{mod.S24}}

		@test mod.g4 isa Cglobal{Cint}
		@test mod.g5 isa Cglobal{mod.S24}
		@test mod.g6 isa Cglobal{Ptr{mod.S24}}

		@test mod.g7 isa Cglobalconst{Cint}
	end

	@testset "functions" begin
	end
end


text:
# introducing the matchin function

match(r"\.com$", "google.com") |> typeof
match(r"\.com$", "w3.org") |> typeof
#=
julia> match(r"\.com$", "google.com") |> typeof
RegexMatch

julia> match(r"\.com$", "w3.org") |> typeof
Nothing
=#

# multiple matches scenario
url = "http://google.com"
if match(r"\.com$", url) !== nothing
    # do something about .com sites
else if match(r"\.org$", url) !== nothing
    # do something about .org sites
else
    # do something different
end

# hypothetically, what would we do if `match` uses exception instead?
url = "http://google.com"
try
    match(r"\.com$", url)
    # do something about .com sites
catch ex1
    try
        match(r"\.org$", url)
        # do something about .org sites
    catch ex2
        # do something different
    end
end


text:
"""
* `H : number of grid points inside wavelength
	* choose 5 for 4th order scheme
* `epsilon` : Courant number
"""
function check_stability(pa::PFdtd, verbose=true; H=div(20,_fd.order), epsilon=inv(sqrt(ndims(pa.c.medium))))
	mgrid=pa.c.medium.mgrid
	freqmin=pa.c.fc[:freqmin]
	freqmax=pa.c.fc[:freqmax]
	attrib_mod=pa.c.attrib_mod
	bounds=pa.c.medium.bounds

	δ=step.(mgrid)
	δt=pa.c.fc[:dt]

	if(isa(attrib_mod, FdtdElastic))
		vmin=minimum(filter(x->x ≠ 0, pa.c.medium[:vs])) # vs condition overrides (other than zero?)
		vmax=sqrt(bounds[:vp][2]^2 + bounds[:vs][2]^2) # see Virieux (1986)

	else
		vmin=bounds[:vp][1]
		vmax=bounds[:vp][2]
	end

	# check spatial sampling, i.e., number of grid points in minimum wavelength
	if(isa(attrib_mod, FdtdElastic))
		δs_temp=round(vmin/H/freqmax,digits=2);
	else
		δs_temp=round(vmin/H/freqmax,digits=2);
	end
	δs_max = maximum(δ)
	str1=@sprintf("%0.2e",δs_max)
	str2=@sprintf("%0.2e",δs_temp)
	if(str1 ≠ str2)
		if(all(δs_max .> δs_temp))
			@warn "decrease maximum spatial sampling ($str1) below $str2"
		else
			if(verbose)
				@info "spatial sampling ($str1) can be as high as $str2"
			end
		end
	end

	# check time sampling, i.e., Courant criteria
	δs_min = minimum(δ)
	δt_temp=epsilon*δs_min/vmax
	str1=@sprintf("%0.2e", δt)
	str2=@sprintf("%0.2e", δt_temp)
	if(str1 ≠ str2)
		if(all(δt .> δt_temp))
			@warn "decrease time sampling ($str1) below $str2"
		else
			if(verbose)
				@info "time sampling ($str1) can be as high as $str2"
			end
		end
	end
end


text:
using CSV
using DataFrames
using Dates
include("../src/TrueSkill.jl")
using .TrueSkill
global const ttt = TrueSkill

data = CSV.read("input/history.csv", DataFrame)
times = Dates.value.(data[:,"time_start"] .- Date("1900-1-1"))
composition = [ r.double == "t" ? [[r.w1_id,r.w2_id],[r.l1_id,r.l2_id]] : [[r.w1_id],[r.l1_id]] for r in eachrow(data) ]

sigmas = [i for i in 1.0:0.05:2.0]
gammas = [i for i in 0.033:0.0005:0.042]

df = DataFrame()
df[!,:rownames] = gammas
for s in sigmas
    df[!,string(s)] .= 0.0
end

for s in 1:length(sigmas)
    for g in 1:length(gammas)
        h = ttt.History(composition=composition, times = times, sigma = sigmas[s], gamma = gammas[g])
        ttt.convergence(h,epsilon=0.01, iterations=5)
        df[g,string(sigmas[s])] = ttt.log_evidence(h)
    end
end

CSV.write("output/atp_optimization.csv", df; header=true)

composition = [ r.double == "t" ? [[r.w1_id,r.w1_id*r.ground,r.w2_id,r.w2_id*r.ground],[r.l1_id,r.l1_id*r.ground,r.l2_id,r.l2_id*r.ground]] : [[r.w1_id,r.w1_id*r.ground],[r.l1_id,r.l1_id*r.ground]] for r in eachrow(data) ]

players = Set(vcat(([ r.double == "t" ? [[r.w1_id,r.w2_id],[r.l1_id,r.l2_id]] : [[r.w1_id],[r.l1_id]] for r in eachrow(data) ]...)...))

sigmas = [i for i in 0.25:0.05:0.45]
gammas = [i for i in 0.000:0.001:0.002]

df = DataFrame()
df[!,:rownames] = gammas
for s in sigmas
    df[!,string(s)] .= 0.0
end

#s=0.35;g=0.0; sp=1.6; gp=0.036
function scenario(s,g, sp, gp)
    priors = Dict([(p, ttt.Player(ttt.Gaussian(0., sp), 1.0, gp) ) for p in players])
    h = ttt.History(composition=composition, times = times, sigma = s, gamma = g, beta = 0.0, priors = priors)
    ttt.convergence(h,epsilon=0.01, iterations=2)
    return ttt.log_evidence(h,online=false,agents=["n409"])
end

#scenario(0.35, 0.000, 1.6, 0.036)
for s in 1:length(sigmas)
    for g in 1:length(gammas)
        GC.gc()
        df[g,string(sigmas[s])] = scenario(sigmas[s],gammas[g], 1.6, 0.036)
        println("s:",sigmas[s], ", g:", gammas[g], ". ", df[g,string(sigmas[s])])
    end
end

CSV.write("output/atp_ground_optimization.csv", df; header=true)


text:
println("Calculating discharges")
p1,Q1 = computedischarge(d1,csd)
p2,Q2 = computedischarge(d2,csd)


text:
mutable struct DeldirArguments
    # The variables for the master Fortran subroutine
    # The variable names are copied from the R package
    x::Vector{Float64}
    y::Vector{Float64}
    rw::Vector{Float64}
    npd::Vector{Int32}
    ntot::Vector{Int32}
    nadj::Vector{Int32}
    madj::Vector{Int32}
    tx::Vector{Float64}
    ty::Vector{Float64}
    epsilon::Float64
    delsgs::Vector{Float64}
    ndel::Vector{Int32}
    delsum::Vector{Float64}
    dirsgs::Vector{Float64}
    ndir::Vector{Int32}
    dirsum::Vector{Float64}
    nerror::Vector{Int32}

    # The variables for sorting the points
    indices::Vector{Int32}
    reverse_indices::Vector{Int32}

    function DeldirArguments(x, y, rw, npd, ntot, nadj, madj, tx, ty, epsilon, delsgs, ndel,
                             delsum, dirsgs, ndir, dirsum, nerror, indices, reverse_indices)
        if length(x) != length(y)
            throw(DimensionMismatch("Coordinate vectors must be of equal length"))
        end

        if epsilon < eps(Float64)
            throw(DomainError(epsilon, "Must be at least `eps(Float64)`"))
        end

        min_x, max_x = extrema(x)
        min_y, max_y = extrema(y)
        if min_x < rw[1] || max_x > rw[2] || min_y < rw[3] || max_y > rw[4]
            throw(DomainError(rw, "Boundary window is too small"))
        end

        new(x, y, rw, npd, ntot, nadj, madj, tx, ty, epsilon, delsgs,
            ndel, delsum, dirsgs, ndir, dirsum, nerror, indices, reverse_indices)
    end
end

function DeldirArguments(x, y, rw, epsilon)
    unique_x, unique_y = remove_duplicates(x, y)
    num_points = length(unique_x)

    # Dummy points: Ignored!
    ndm = 0
    npd = num_points + ndm

    # The total number of points
    ntot = npd + 4

    indices, reverse_indices = sortperm_points!(unique_x, unique_y, rw)
    X = [zeros(4); unique_x; zeros(4)]
    Y = [zeros(4); unique_y; zeros(4)]

    # Set up fixed dimensioning constants
    ntdel = 4*npd
    ntdir = 3*npd

    # Set up dimensioning constants which might need to be increased
    madj_val = max(20, ceil(Int32, 3*sqrt(ntot)))
    madj = Int32[madj_val]
    tadj = (madj_val + 1)*(ntot + 4)
    ndel = Int32[madj_val*(madj_val + 1)/2]
    tdel = 6*ndel[]
    ndir = copy(ndel)
    tdir = 10*ndir[]

    nadj   = zeros(Int32, tadj)
    tx     = zeros(Float64, npd)
    ty     = zeros(Float64, npd)
    delsgs = zeros(Float64, tdel)
    delsum = zeros(Float64, ntdel)
    dirsgs = zeros(Float64, tdir)
    dirsum = zeros(Float64, ntdir)
    nerror = Int32[1]

    DeldirArguments(X, Y, rw, [Int32(npd)], [Int32(ntot)], nadj, madj, tx, ty, epsilon,
                    delsgs, ndel, delsum, dirsgs, ndir, dirsum, nerror, indices, reverse_indices)
end

function sortperm_points!(x, y, rw)
    n = length(x)
    if n != length(y)
        DimensionMismatch("x and y must have the same length")
    end

    tx = similar(x)
    ty = similar(y)

    indices = Vector{Int32}(undef, n)
    reverse_indices = similar(indices)
    ilst = similar(indices)
    nerror = Int32[1]

    ccall((:binsrt_, Deldir_jll.libdeldir), Cvoid,
          (Ref{Float64}, Ref{Float64}, Ref{Float64}, Ref{Int32}, Ref{Int32},
           Ref{Int32}, Ref{Float64}, Ref{Float64}, Ref{Int32}, Ref{Int32}),
          x, y, rw, n, indices,
          reverse_indices, tx, ty, ilst, nerror
    )

    if nerror[] > 0
        error("Mismatch between number of points and number of sorted points")
    end

    return indices, reverse_indices
end

function error_handling!(da::DeldirArguments)
    error_number = da.nerror[]

    if error_number == 4
        madj_val = ceil(Int32, 1.2*da.madj[])
        da.madj = Int32[madj_val]

        tadj = (madj_val + 1)*(da.ntot[] + 4)
        resize!(da.nadj, tadj)

        ndel_val = max(da.ndel[], div(madj_val*(madj_val + 1), 2))
        da.ndel = Int32[ndel_val]
        da.ndir = copy(da.ndel)

        tdel = 6*ndel_val
        resize!(da.delsgs, tdel)

        tdir = 8*ndel_val
        resize!(da.dirsgs, tdir)

        @info "Fortran error $(error_number). Increasing madj to $(da.madj[])"
    elseif error_number == 12
        error("Vertices of triangle are collinear")
    elseif error_number == 14 || error_number == 15
        ndel_val = ceil(Int32, 1.2*ndel[])
        da.ndel = Int32[ndel_val]
        da.ndir = copy(ndel)

        tdel = 6*ndel_val
        resize!(da.delsgs, tdel)

        tdir = 8*ndel_val
        resize!(da.dirsgs, tdir)

        @info "Fortran error $(error_number). Increasing ndel & ndir to $(da.ndel[])"
    elseif error_number > 1
        error("From `deldir` Fortran, nerror = ", error_number)
    end

    return(da)
end

function get_delaunay(da::DeldirArguments)
    num_del = Int64(da.ndel[])
    delsgs  = reshape(da.delsgs[1:6*num_del], 6, num_del) |> transpose

    del_df = DataFrames.DataFrame(
        [Float64, Float64, Float64, Float64, Int64, Int64],
        ["x1", "y1", "x2", "y2", "ind1", "ind2"],
        num_del
    )

    del_df[!, "x1"] = delsgs[:, 1]
    del_df[!, "y1"] = delsgs[:, 2]
    del_df[!, "x2"] = delsgs[:, 3]
    del_df[!, "y2"] = delsgs[:, 4]

    ind1 = round.(Int, delsgs[:, 5])
    del_df[!, "ind1"] = da.reverse_indices[ind1]

    ind2 = round.(Int, delsgs[:, 6])
    del_df[!, "ind2"] = da.reverse_indices[ind2]

    return del_df
end

function get_voronoi(da::DeldirArguments)
    num_dir = Int64(da.ndir[])
    vor  = reshape(da.dirsgs[1:10*num_dir], 10, num_dir) |> transpose

    vor_df = DataFrames.DataFrame(
        [Float64, Float64, Float64, Float64, Int64, Int64, Bool, Bool, Int64, Int64],
        ["x1", "y1", "x2", "y2", "ind1", "ind2", "bp1", "bp2", "thirdv1", "thirdv2"],
        num_dir
    )

    vor_df[!, "x1"] = vor[:, 1]
    vor_df[!, "y1"] = vor[:, 2]
    vor_df[!, "x2"] = vor[:, 3]
    vor_df[!, "y2"] = vor[:, 4]
    vor_df[!, "bp1"] = vor[:, 7] .== 1
    vor_df[!, "bp2"] = vor[:, 8] .== 1

    ind1 = round.(Int, vor[:, 5])
    vor_df[!, "ind1"] = da.reverse_indices[ind1]

    ind2 = round.(Int, vor[:, 6])
    vor_df[!, "ind2"] = da.reverse_indices[ind2]

    thirdv1 = Int.(vor[:, 9])
    idx1 = thirdv1 .>= 0
    thirdv1[idx1] = da.reverse_indices[thirdv1[idx1]]
    vor_df[!, "thirdv1"] = thirdv1

    thirdv2 = Int.(vor[:, 10])
    idx2 = thirdv2 .>= 0
    thirdv2[idx2] = da.reverse_indices[thirdv2[idx2]]
    vor_df[!, "thirdv2"] = thirdv2

    return vor_df
end

function get_summary(da::DeldirArguments)
    npd = Int64(da.npd[1])

    summary_df = DataFrames.DataFrame(
        [Float64, Float64, Int64, Float64, Int64, Int64, Float64],
        ["x", "y", "ntri", "del_area", "n_tside", "nbpt", "vor_area"],
        npd
    )

    delsum = reshape(da.delsum, npd, 4)
    delsum = delsum[da.indices, :]
    summary_df[!, "x"] = delsum[:, 1]
    summary_df[!, "y"] = delsum[:, 2]
    summary_df[!, "ntri"] = round.(Int, delsum[:, 3])
    summary_df[!, "del_area"] = delsum[:, 4]

    dirsum = reshape(da.dirsum, npd, 3)
    dirsum = dirsum[da.indices, :]
    summary_df[!, "n_tside"] = round.(Int, dirsum[:, 1])
    summary_df[!, "nbpt"] = round.(Int, dirsum[:, 2])
    summary_df[!, "vor_area"] = dirsum[:, 3]

    return summary_df
end

"""
deldirwrapper!(da::DeldirArguments)

Wrapper for the Fortran code that returns the output undigested.
"""
function deldirwrapper!(da::DeldirArguments)
    # Call Fortran routine
    while da.nerror[] >= 1
        ccall((:master_, Deldir_jll.libdeldir), Cvoid,
            (Ref{Float64}, Ref{Float64}, Ref{Float64}, Ref{Int32}, Ref{Int32},
             Ref{Int32}, Ref{Int32}, Ref{Float64}, Ref{Float64},
             Ref{Float64}, Ref{Float64}, Ref{Int32},
             Ref{Float64}, Ref{Float64}, Ref{Int32}, Ref{Float64}, Ref{Int32}),
            da.x, da.y, da.rw, da.npd, da.ntot,
            da.nadj, da.madj, da.tx, da.ty,
            da.epsilon, da.delsgs, da.ndel,
            da.delsum, da.dirsgs, da.ndir, da.dirsum, da.nerror
        )

        error_handling!(da)
    end

    return da
end

"""
deldir(x::Vector, y::Vector; ...)

Compute the Delaunay triangulation and Voronoi tesselation of the 2D points with x-coordinates `x` and y-coordinates `y`.

Optional arguments are

- `rw`: Boundary rectangle specified as a vector with `[xmin, xmax, ymin, ymax]`. By default, `rw` is the unit rectangle.
- `epsilon`: A value of epsilon used in testing whether a quantity is zeros, mainly in the context of whether points are collinear.
If anomalous errors arise, it is possible that these may averted by adjusting the value of `epsilon` upward or downward.
By default, `epsilon = 1e-9`.

The output are three `DataFrame`s:

###### `delsgs`

- The `x1`, `y1`, `x2` & `y2` entires are the coordinates of the points joined by an edge of a Delaunay triangle.
- The `ind1` and `ind2` entries are the indices of the two points which are joined.

###### `vorsgs`

- The `x1`, `y1`, `x2` & `y2` entires are the coordinates of the endpoints of one the edges of a Voronoi cell.
- The `ind1` and `ind2` entries are the indices of the two points, in the set being triangulated, which are separated by that edge
- The `bp1` entry indicates whether the first endpoint of the corresponding edge of a Voronoi cell is a boundary point (a point on the boundary of the rectangular window).
Likewise for the `bp2` entry and the second endpoint of the edge.
- The `thirdv1` and `thirdv2` columns are the indices of the respective third vertices of the Delaunay triangle whose circumcentres constitute the corresponding endpoints of the edge under consideration.

###### `summary`

- The `x` and `y` entries of each row are the coordinates of the points in the set being triangulated.
- The `ntri` entry are the number of Delaunay triangles emanating from the point.
- The `del_area` entry is `1/3` of the total area of all the Delaunay triangles emanating from the point.
- The `n_tside` entry is the number of sides — within the rectangular window — of the Voronoi cell surrounding the point.
- The `nbpt` entry is the number of points in which the Voronoi cell intersects the boundary of the rectangular window.
- The `vor_area` entry is the area of the Voronoi cell surrounding the point.
"""
function deldir(x::Vector{Float64}, y::Vector{Float64}, rw::Vector = [0.0; 1.0; 0.0; 1.0], epsilon = 1e-9)
    da = DeldirArguments(x, y, rw, epsilon)
    deldirwrapper!(da)

    del = get_delaunay(da)
    vor = get_voronoi(da)
    summary = get_summary(da)

    return del, vor, summary
end


text:
using Clustering, RDatasets, Random, Measures, Plots; pyplot()
Random.seed!(0)

K = 3
df = dataset("cluster", "xclara")
data = copy(convert(Array{Float64}, df)')

seeds = initseeds(:rand, data, K)
xclaraKmeans = kmeans(data, K, init = seeds)

println("Number of clusters: ", nclusters(xclaraKmeans))
println("Counts of clusters: ", counts(xclaraKmeans))

df.Group  = assignments(xclaraKmeans)

p1 = scatter(df[:, :V1], df[:, :V2], c=:blue, msw=0)
     scatter!(df[seeds, :V1], df[seeds, :V2], markersize=12, c=:red, msw=0)

p2 = scatter( df[df.Group .== 1, :V1], df[df.Group .== 1, :V2], c=:blue, msw=0)
     scatter!( df[df.Group .== 2, :V1], df[df.Group .== 2, :V2], c=:red, msw=0)
     scatter!( df[df.Group .== 3, :V1], df[df.Group .== 3, :V2], c=:green, msw=0)

plot(p1,p2,legend=:none,ratio=:equal,
    size=(800,400), xlabel="V1", ylabel="V2", margin = 5mm)

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
"""
Holds all data needed for a documentation example... header, description, and plotting expression (Expr)
"""
mutable struct PlotExample
  header::AbstractString
  desc::AbstractString
  exprs::Vector{Expr}
end

# the _examples we'll run for each
const _examples = PlotExample[

PlotExample("Lines",
    "A simple line plot of the columns.",
    [:(begin
        plot(Plots.fakedata(50,5), w=3)
    end)]
),

PlotExample("Functions, adding data, and animations",
"""
Plot multiple functions.  You can also put the function first, or use the form `plot(f,
xmin, xmax)` where f is a Function or AbstractVector{Function}.\n\nGet series data:
`x, y = plt[i]`.  Set series data: `plt[i] = (x,y)`. Add to the series with
`push!`/`append!`.\n\nEasily build animations.  (`convert` or `ffmpeg` must be available
to generate the animation.)  Use command `gif(anim, filename, fps=15)` to save the
animation.
""",
    [:(begin
        p = plot([sin,cos], zeros(0), leg=false)
        anim = Animation()
        for x in range(0, stop=10π, length=100)
            push!(p, x, Float64[sin(x), cos(x)])
            frame(anim)
        end
    end)]
),

PlotExample("Parametric plots",
    "Plot function pair (x(u), y(u)).",
    [:(begin
        plot(sin, x->sin(2x), 0, 2π, line=4, leg=false, fill=(0,:orange))
    end)]
),

PlotExample("Colors",
"""
Access predefined palettes (or build your own with the `colorscheme` method).
Line/marker colors are auto-generated from the plot's palette, unless overridden.  Set
the `z` argument to turn on series gradients.
""",
    [:(begin
y = rand(100)
plot(0:10:100,rand(11,4),lab="lines",w=3,palette=:grays,fill=0, α=0.6)
scatter!(y, zcolor=abs.(y.-0.5), m=(:heat,0.8,Plots.stroke(1,:green)), ms=10*abs.(y.-0.5).+4,
         lab="grad")
    end)]
),

PlotExample("Global",
"""
Change the guides/background/limits/ticks.  Convenience args `xaxis` and `yaxis` allow
you to pass a tuple or value which will be mapped to the relevant args automatically.
The `xaxis` below will be replaced with `xlabel` and `xlims` args automatically during
the preprocessing step. You can also use shorthand functions: `title!`, `xaxis!`,
`yaxis!`, `xlabel!`, `ylabel!`, `xlims!`, `ylims!`, `xticks!`, `yticks!`
""",
    [:(begin
using Statistics
y = rand(20,3)
plot(y, xaxis=("XLABEL",(-5,30),0:2:20,:flip), background_color = RGB(0.2,0.2,0.2),
     leg=false)
hline!(mean(y, dims = 1)+rand(1,3), line=(4,:dash,0.6,[:lightgreen :green :darkgreen]))
vline!([5,10])
title!("TITLE")
yaxis!("YLABEL", :log10)
    end)]
),

# PlotExample("Two-axis",
#             "Use the `axis` arguments.",
#             [
#               :(plot(Vector[randn(100), randn(100)*100], axis = [:l :r], ylabel="LEFT", yrightlabel="RIGHT", xlabel="X", title="TITLE"))
#             ]),

PlotExample("Images",
    "Plot an image.  y-axis is set to flipped",
    [:(begin
    import FileIO
    path = download("http://juliaplots.org/PlotReferenceImages.jl/Plots/pyplot/0.7.0/ref1.png")
    img = FileIO.load(path)
        plot(img)
    end)]
),

PlotExample("Arguments",
"""
Plot multiple series with different numbers of points.  Mix arguments that apply to all
series (marker/markersize) with arguments unique to each series (colors).  Special
arguments `line`, `marker`, and `fill` will automatically figure out what arguments to
set (for example, we are setting the `linestyle`, `linewidth`, and `color` arguments with
`line`.)  Note that we pass a matrix of colors, and this applies the colors to each
series.
""",
    [:(begin
        ys = Vector[rand(10), rand(20)]
        plot(ys, color=[:black :orange], line=(:dot,4), marker=([:hex :d],12,0.8,Plots.stroke(3,:gray)))
    end)]
),

PlotExample("Build plot in pieces",
    "Start with a base plot...",
    [:(begin
        plot(rand(100)/3, reg=true, fill=(0,:green))
    end)]
),

PlotExample("",
    "and add to it later.",
    [:(begin
        scatter!(rand(100), markersize=6, c=:orange)
    end)]
),

PlotExample("Histogram2D",
    "",
    [:(begin
        histogram2d(randn(10000), randn(10000), nbins=20)
    end)]
),

PlotExample("Line types",
    "",
    [:(begin
        linetypes = [:path :steppre :steppost :sticks :scatter]
        n = length(linetypes)
        x = Vector[sort(rand(20)) for i in 1:n]
        y = rand(20,n)
        plot(x, y, line=(linetypes,3), lab=map(string,linetypes), ms=15)
    end)]
),

PlotExample("Line styles",
    "",
    [:(begin
styles = filter(s -> s in Plots.supported_styles(),
                [:solid, :dash, :dot, :dashdot, :dashdotdot])
styles = reshape(styles, 1, length(styles)) # Julia 0.6 unfortunately gives an error when transposing symbol vectors
n = length(styles)
y = cumsum(randn(20,n), dims = 1)
plot(y, line = (5, styles), label = map(string,styles), legendtitle = "linestyle")
             end)]
),

PlotExample("Marker types",
    "",
    [:(begin
        markers = filter(m -> m in Plots.supported_markers(), Plots._shape_keys)
        markers = reshape(markers, 1, length(markers))
        n = length(markers)
        x = range(0, stop=10, length=n+2)[2:end-1]
        y = repeat(reshape(reverse(x),1,:), n, 1)
        scatter(x, y, m=(8,:auto), lab=map(string,markers), bg=:linen, xlim=(0,10), ylim=(0,10))
    end)]
),

PlotExample("Bar",
    "`x` is the midpoint of the bar. (todo: allow passing of edges instead of midpoints)",
    [:(begin
        bar(randn(99))
    end)]
),

PlotExample("Histogram",
    "",
    [:(begin
        histogram(randn(1000), bins = :scott, weights = repeat(1:5, outer = 200))
    end)]
),

PlotExample("Subplots",
"""
Use the `layout` keyword, and optionally the convenient `@layout` macro to generate
arbitrarily complex subplot layouts.
""",
    [:(begin
l = @layout([a{0.1h}; b [c;d e]])
plot(randn(100,5), layout=l, t=[:line :histogram :scatter :steppre :bar], leg=false,
     ticks=nothing, border=:none)
    end)]
),

PlotExample("Adding to subplots",
"""
Note here the automatic grid layout, as well as the order in which new series are added
to the plots.
""",
    [:(begin
plot(Plots.fakedata(100,10), layout=4, palette=[:grays :blues :heat :lightrainbow],
     bg_inside=[:orange :pink :darkblue :black])
    end)]
),

PlotExample("",
    "",
    [:(begin
        using Random
        Random.seed!(111)
        plot!(Plots.fakedata(100,10))
    end)]
),

PlotExample("Open/High/Low/Close",
"""
Create an OHLC chart.  Pass in a list of (open,high,low,close) tuples as your `y`
argument.  This uses recipes to first convert the tuples to OHLC objects, and
subsequently create a :path series with the appropriate line segments.
""",
    [:(begin
n=20
hgt=rand(n).+1
bot=randn(n)
openpct=rand(n)
closepct=rand(n)
y = OHLC[(openpct[i]*hgt[i]+bot[i], bot[i]+hgt[i], bot[i],
          closepct[i]*hgt[i]+bot[i]) for i in 1:n]
ohlc(y)
    end)]
),

PlotExample("Annotations",
"""
The `annotations` keyword is used for text annotations in data-coordinates.  Pass in a
tuple (x,y,text) or a vector of annotations.  `annotate!(ann)` is shorthand for `plot!(;
annotation=ann)`.  Series annotations are used for annotating individual data points.
They require only the annotation... x/y values are computed.  A `PlotText` object can be
build with the method `text(string, attr...)`, which wraps font and color attributes.
""",
    [:(begin
y = rand(10)
plot(y, annotations = (3,y[3], Plots.text("this is #3",:left)), leg=false)
annotate!([(5, y[5], Plots.text("this is #5",16,:red,:center)),
          (10, y[10], Plots.text("this is #10",:right,20,"courier"))])
scatter!(range(2, stop=8, length=6), rand(6), marker=(50,0.2,:orange),
         series_annotations = ["series","annotations","map","to","series",
                               Plots.text("data",:green)])
    end)]
),

PlotExample("Custom Markers",
"""A `Plots.Shape` is a light wrapper around vertices of a polygon.  For supported
backends, pass arbitrary polygons as the marker shapes.  Note: The center is (0,0) and
the size is expected to be rougly the area of the unit circle.
""",
    [:(begin
verts = [(-1.0,1.0),(-1.28,0.6),(-0.2,-1.4),(0.2,-1.4),(1.28,0.6),(1.0,1.0),
         (-1.0,1.0),(-0.2,-0.6),(0.0,-0.2),(-0.4,0.6),(1.28,0.6),(0.2,-1.4),
         (-0.2,-1.4),(0.6,0.2),(-0.2,0.2),(0.0,-0.2),(0.2,0.2),(-0.2,-0.6)]
x = 0.1:0.2:0.9
y = 0.7rand(5).+0.15
plot(x, y, line = (3,:dash,:lightblue), marker = (Shape(verts),30,RGBA(0,0,0,0.2)),
     bg=:pink, fg=:darkblue, xlim = (0,1), ylim=(0,1), leg=false)
    end)]
),

PlotExample("Contours",
"""
Any value for fill works here.  We first build a filled contour from a function, then an
unfilled contour from a matrix.
""",
    [:(begin
        x = 1:0.5:20
        y = 1:0.5:10
        f(x,y) = (3x+y^2)*abs(sin(x)+cos(y))
        X = repeat(reshape(x,1,:), length(y), 1)
        Y = repeat(y, 1, length(x))
        Z = map(f, X, Y)
        p1 = contour(x, y, f, fill=true)
        p2 = contour(x, y, Z)
        plot(p1, p2)
    end)]
),

PlotExample("Pie",
    "",
    [:(begin
        x = ["Nerds", "Hackers", "Scientists"]
        y = [0.4, 0.35, 0.25]
        pie(x, y, title="The Julia Community", l=0.5)
    end)]
),

PlotExample("3D",
    "",
    [:(begin
        n = 100
        ts = range(0, stop=8π, length=n)
        x = ts .* map(cos,ts)
        y = 0.1ts .* map(sin,ts)
        z = 1:n
        plot(x, y, z, zcolor=reverse(z), m=(10,0.8,:blues,Plots.stroke(0)), leg=false, cbar=true, w=5)
        plot!(zeros(n),zeros(n),1:n, w=10)
    end)]
),

PlotExample("DataFrames",
    "Plot using DataFrame column symbols.",
    [:(using StatsPlots), # can't be inside begin block because @df gets expanded first
     :(begin
        import RDatasets
        iris = RDatasets.dataset("datasets", "iris")
        @df iris scatter(:SepalLength, :SepalWidth, group=:Species,
            title = "My awesome plot", xlabel = "Length", ylabel = "Width",
            marker = (0.5, [:cross :hex :star7], 12), bg=RGB(.2,.2,.2))
    end)]
),

PlotExample("Groups and Subplots",
    "",
    [:(begin
        group = rand(map(i->"group $i",1:4),100)
        plot(rand(100), layout=@layout([a b;c]), group=group,
            linetype=[:bar :scatter :steppre], linecolor = :match)
    end)]
),

PlotExample("Polar Plots",
    "",
    [:(begin
        Θ = range(0, stop=1.5π, length=100)
        r = abs.(0.1randn(100)+sin.(3Θ))
        plot(Θ, r, proj=:polar, m=2)
    end)]
),

PlotExample("Heatmap, categorical axes, and aspect_ratio",
    "",
    [:(begin
        xs = [string("x",i) for i=1:10]
        ys = [string("y",i) for i=1:4]
        z = float((1:4)*reshape(1:10,1,:))
        heatmap(xs, ys, z, aspect_ratio=1)
    end)]
),

PlotExample("Layouts, margins, label rotation, title location",
    "",
    [:(begin
        using Plots.PlotMeasures # for Measures, e.g. mm and px
        plot(rand(100,6),layout=@layout([a b; c]),title=["A" "B" "C"],
                        title_location=:left, left_margin=[20mm 0mm],
                        bottom_margin=10px, xrotation=60)
    end)]
),

PlotExample("Boxplot and Violin series recipes",
    "",
    [:(using StatsPlots), # can't be inside begin block because @df gets expanded first
     :(begin
        import RDatasets
        singers = RDatasets.dataset("lattice", "singer")
        @df singers violin(:VoicePart, :Height, line = 0, fill = (0.2, :blue))
        @df singers boxplot!(:VoicePart, :Height, line = (2,:black), fill = (0.3, :orange))
    end)]
),

PlotExample("Animation with subplots",
    "The `layout` macro can be used to create an animation with subplots.",
    [:(begin
        l = @layout([[a; b] c])
        p = plot(plot([sin,cos],1,leg=false),
        scatter([atan,cos],1,leg=false),
        plot(log,1,xlims=(1,10π),ylims=(0,5),leg=false),layout=l)

        anim = Animation()
        for x = range(1, stop=10π, length=100)
          plot(push!(p,x,Float64[sin(x),cos(x),atan(x),cos(x),log(x)]))
          frame(anim)
        end
    end)]
),

PlotExample("Spy",
"""
For a matrix `mat` with unique nonzeros `spy(mat)` returns a colorless plot. If `mat` has
various different nonzero values, a colorbar is added. The colorbar can be disabled with
`legend = nothing`.
""",
    [:(begin
    using SparseArrays
    a = spdiagm(0 => ones(50), 1 => ones(49), -1 => ones(49), 10 => ones(40), -10 => ones(40))
    b = spdiagm(0 => 1:50, 1 => 1:49, -1 => 1:49, 10 => 1:40, -10 => 1:40)
    plot(spy(a), spy(b), title = ["Unique nonzeros" "Different nonzeros"])
    end)]
),

PlotExample("Magic grid argument",
"""
The grid lines can be modified individually for each axis with the magic `grid` argument.
""",
    [:(begin
    x = rand(10)
    p1 = plot(x, title = "Default looks")
    p2 = plot(x, grid = (:y, :olivedrab, :dot, 1, 0.9), title = "Modified y grid")
    p3 = plot(deepcopy(p2), title = "Add x grid")
    xgrid!(p3, :on, :cadetblue, 2, :dashdot, 0.4)
    plot(p1, p2, p3, layout = (1, 3), label = "", fillrange = 0, fillalpha = 0.3)
    end)]
),

PlotExample("Framestyle",
"""
The style of the frame/axes of a (sub)plot can be changed with the `framestyle`
attribute. The default framestyle is `:axes`.
""",
    [:(begin
    scatter(fill(randn(10), 6), fill(randn(10), 6),
        framestyle = [:box :semi :origin :zerolines :grid :none],
        title = [":box" ":semi" ":origin" ":zerolines" ":grid" ":none"],
        color = permutedims(1:6), layout = 6, label = "", markerstrokewidth = 0,
        ticks = -2:2)
    end)]
),

PlotExample("Lines and markers with varying colors",
"""
You can use the `line_z` and `marker_z` properties to associate a color with
each line segment or marker in the plot.
""",
    [:(begin
        t = range(0, stop=1, length=100)
        θ = 6π .* t
        x = t .* cos.(θ)
        y = t .* sin.(θ)
        p1 = plot(x, y, line_z=t, linewidth=3, legend=false)
        p2 = scatter(x, y, marker_z=(x,y)->x+y, color=:bluesreds, legend=false)
        plot(p1, p2)
    end)]
),

PlotExample("Portfolio Composition maps",
"""
see: http://stackoverflow.com/a/37732384/5075246
""",
    [:(begin
    using Random
    Random.seed!(111)
    tickers = ["IBM", "Google", "Apple", "Intel"]
    N = 10
    D = length(tickers)
    weights = rand(N,D)
    weights ./= sum(weights, dims = 2)
    returns = sort!((1:N) + D*randn(N))

    portfoliocomposition(weights, returns, labels = permutedims(tickers))
    end)]
),

PlotExample("Ribbons",
    """
    Ribbons can be added to lines via the `ribbon` keyword;
    you can pass a tuple of arrays (upper and lower bounds),
    a single Array (for symmetric ribbons), a Function, or a number.
    """,
    [:(begin
        plot(
            plot(0:10; ribbon = (LinRange(0, 2, 10), LinRange(0, 1, 10))),
            plot(0:10; ribbon = 0:0.5:5),
            plot(0:10; ribbon = sqrt),
            plot(0:10; ribbon = 1),
        )
    end)]
),

PlotExample("Histogram2D (complex values)",
    "",
    [:(begin
        n = 10_000
        x = exp.(0.1randn(n) .+ randn(n).*(im))
        histogram2d(x, nbins=(20,40), show_empty_bins=true,
                    normed=true, aspect_ratio=1)
    end)]
),

PlotExample("Unconnected lines using `missing` or `NaN`",
"""
Missing values and non-finite values, including `NaN`, are not plotted.
Instead, lines are separated into segments at these values.
""",
    [:(begin
        x,y = [1,2,2,1,1], [1,2,1,2,1]
        plot(
              plot([rand(5); NaN; rand(5); NaN; rand(5)]),
              plot([1,missing,2,3], marker=true),
              plot([x; NaN; x.+2], [y; NaN; y.+1], arrow=2),
              plot([1, 2+3im, Inf, 4im, 3, -Inf*im, 0, 3+3im], marker=true),
              legend=false
        )
    end)]
),

]

# Some constants for PlotDocs and PlotReferenceImages
_animation_examples = [2, 30]
_backend_skips = Dict(
    :gr => [25, 30],
    :pyplot => [25, 30],
    :plotlyjs => [2, 21, 24, 25, 30, 31],
    :pgfplots => [2, 5, 6, 10, 16, 20, 22, 23, 25, 28, 30],
)

# ---------------------------------------------------------------------------------

# make and display one plot
function test_examples(pkgname::Symbol, idx::Int; debug = false, disp = true)
  Plots._debugMode.on = debug
  @info("Testing plot: $pkgname:$idx:$(_examples[idx].header)")
  backend(pkgname)
  backend()

  # prevent leaking variables (esp. functions) directly into Plots namespace
  m = Module(:PlotExampleModule)
  Base.eval(m, :(using Plots))
  map(exprs -> Base.eval(m, exprs), _examples[idx].exprs)

  plt = current()
  if disp
    gui(plt)
  end
  plt
end

# generate all plots and create a dict mapping idx --> plt
"""
test_examples(pkgname[, idx]; debug = false, disp = true, sleep = nothing,
                                        skip = [], only = nothing

Run the `idx` test example for a given backend, or all examples if `idx`
is not specified.
"""
function test_examples(pkgname::Symbol; debug = false, disp = true, sleep = nothing,
                                        skip = [], only = nothing)
  Plots._debugMode.on = debug
  plts = Dict()
  for i in 1:length(_examples)
    only !== nothing && !(i in only) && continue
    i in skip && continue
    try
      plt = test_examples(pkgname, i, debug=debug, disp=disp)
      plts[i] = plt
    catch ex
      # TODO: put error info into markdown?
      @warn("Example $pkgname:$i:$(_examples[i].header) failed with: $ex")
    end
    if sleep !== nothing
        Base.sleep(sleep)
    end
  end
  plts
end


text:
using QuadGG
using Test

@testset "QuadGG" begin
    @testset "GaussLobatto" begin
        @testset "Continuous" begin
            @test adaptlob(sqrt,0,1;tol=1e-8) ≈ 2/3 atol=1e-8
            @test !isapprox(adaptlob(sqrt,0,1;tol=1e-3),2/3, atol=1e-8)
            @test adaptlob(x -> exp(im*x), 0,1) ≈ (exp(1im)-1)/im
        end

        @testset "Discontinuous" begin
            fun1(x) = ifelse(x<1/2,0,1)
            @test adaptlob(fun1,0,1;tol=1e-8) ≈ 1/2 atol=1e-8
        end
    end
    @testset "Simpson" begin
        @testset "Continuous" begin
            @test adaptsim(sqrt,0,1;tol=1e-8) ≈ 2/3 atol=1e-7
            @test !isapprox(adaptsim(sqrt,0,1;tol=1e-3),2/3, atol=1e-8)
            @test adaptsim(x -> exp(im*x), 0,1) ≈ (exp(1im)-1)/im
        end

        @testset "Discontinuous" begin
            fun1(x) = ifelse(x<1/2,0,1)
            @test adaptsim(fun1,0,1;tol=1e-8) ≈ 1/2 atol=1e-5

            fun2(x) = ifelse(x < 1,x+1, ifelse(x>3, 2,3-x))
            @test adaptsim(fun2,0,5;tol=1e-6,trace=true) ≈ 7.5 atol=1e-5
        end
    end
end


text:
"""
complex bulk modulus K(ω)
then Q is given by Re(K(ω))/Im(K(ω))
* `nsls` : number of standard linear solids
*
"""
function get_relaxation_times(Q, nsls, fmin, fmax)
	tau_epsilon=zeros(nsls, size(Q)...)
	tau_sigma=zeros(nsls, size(Q)...)

	nom=100
	om=range(fmin, stop=fmax, length=nom);

	# unifrmly distribute tau_sigma over the pass band (eq. 27)
	τ_σ0=inv.(exp.(range(log(fmin), stop=log(fmax), length=nsls)))

	for ix in 1:size(Q,2)
		for iz in 1:size(Q,1)
			Q1=Q[iz,ix]

			for i in 1:nsls
				tau_sigma[i, iz, ix] = τ_σ0[i]
			end

			# see eq. 24
			A=hcat([om .* inv.(1. .+ om.*om .* tau_sigma[i,iz,ix].*tau_sigma[i,iz,ix]) for i in 1:nsls]...)

			# our objective is to optimize Q such that it is close to Q_0, for all ω
			τ_ϵ0 = A\fill(inv(Q1),nom)

			for i in 1:nsls
				tau_epsilon[i,iz,ix] = tau_sigma[i,iz,ix] .+ τ_ϵ0[i]
			end
		end
    end
    return tau_sigma, tau_epsilon
end

"""
Compute that factor that converts unrelaxed bulk modulus to relaxed bulk modulus.
"""
function relaxed_conversion_factor(tau_sigma, tau_epsilon)
	@assert size(tau_epsilon) == size(tau_sigma)
	nsls, nz, nx = size(tau_sigma)
	K_relaxation_factor=zeros(nz,nx)
	for ix in 1:nx
		for iz in 1:nz
			stau=0.0
			for i in 1:nsls
				stau  += (tau_epsilon[i,iz,ix]*inv(tau_sigma[i,iz,ix]) - 1)
			end
			K_relaxation_factor[iz,ix] = 1.0 - stau
		end
	end
	return K_relaxation_factor
end

# return complex bulk modulus, for a frequency domain solver (currently implemented for homogeneous media)
function complexK(Ku, freqs, tau_sigma, tau_epsilon)
    fac=relaxed_conversion_factor(tau_sigma, tau_epsilon)[1]
    nsls=size(tau_sigma,1)

    Kc=complex.(zeros(length(freqs)), zeros(length(freqs)))
    for (ifreq,freq) in enumerate(freqs)
        s=complex(0.0,0.0)
        for isls in 1:nsls
            s+=(1.0 + im*freq*tau_epsilon[isls,1,1])*inv(1.0 + im*freq*tau_sigma[isls,1,1])
        end
		Kc[ifreq]=(1.0 + inv(nsls)*(-nsls + s))*Ku#*inv(fac)
		s=0.0
		for isls in 1:nsls
			s+=tau_epsilon[isls,1,1]*inv(tau_sigma[isls,1,1])
		end
		Kc[ifreq] *= inv(1 + inv(nsls)*(-nsls + s))
	end
	return Kc

end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export dal_dashaladinlite

"""
    dal_dashaladinlite(;kwargs...)

A DashAladinLite component.
ExampleComponent is an example component.
It takes a property, `label`, and
displays it.
It renders an input with the property `value`
which is editable by the user.
Keyword arguments:
- `id` (String; optional): The ID used to identify this component in Dash callbacks.
- `label` (String; required): A label that will be printed when this component is rendered.
- `value` (String; optional): The value displayed in the input.
"""
function dal_dashaladinlite(; kwargs...)
        available_props = Symbol[:id, :label, :value]
        wild_props = Symbol[]
        return Component("dal_dashaladinlite", "DashAladinLite", "dash_aladin_lite", available_props, wild_props; kwargs...)
end


max_stars_repo_path: test/4 Trig functions/4.1 Sine/4.1.9 trig^m (a+b sin^n+c sin^(2 n))^p.jl


text:
module TestRadTrans

using Test
using Dates: DateTime
using Canopy.RadTrans.Blackbody
using Canopy.RadTrans.Solar

@time @testset "Canopy.RadTrans.Blackbody" begin
    # ephoton
    @test isapprox(ephoton(450e-9), 4.414e-19, rtol=1e-3)  # blue
    @test isapprox(ephoton(700e-9), 2.838e-19, rtol=1e-3)  # red
    # energy2photon
    @test isapprox(energy2photon(600, 0.5e-6), 2508e-6, rtol=1e-3)
    # planck
    @test isapprox(planck(500e-9, 5772), 2.624e13, rtol=1e-3)  # sun
    @test isapprox(planck(10e-6, 288), 8.114e6, rtol=1e-3)  # earth
    # stefan_boltzmann
    @test isapprox(stefan_boltzmann(5772), 6.294e7, rtol=1e-3)  # sun
    @test isapprox(stefan_boltzmann(255), 239.8, rtol=1e-3)  # earth
    # blackbody_temp
    @test isapprox(blackbody_temp(240), 255, rtol=1e-3)  # earth
    @test isapprox(blackbody_temp(149.6), 226.6, rtol=1e-3)  # venus
end

@time @testset "Canopy.RadTrans.Solar" begin
    # eccentricity
    @test isapprox(eccentricity(-1), 0.01712, rtol=1e-3)  # 1000 AD
    @test isapprox(eccentricity(DateTime(1000)), eccentricity(-1), rtol=1e-6)
    @test isapprox(eccentricity(1), 0.01627, rtol=1e-3)  # 3000 AD
    @test isapprox(eccentricity(DateTime(3000)), eccentricity(1), rtol=1e-6)

    # atmos_refrac
    @test atmos_refrac(89.) == 0.0
    angle = 40.
    # compare with a simpler approximation by Bennett (1982) J. Navigat.
    # accuracy = 0.1 arcmin
    @test isapprox(atmos_refrac(angle),
                   cotd(angle + 7.31 / (angle + 4.4)) / 60., atol=0.1 / 60.)

    # solar_angle
    solar_angle_result = solar_angle(
        DateTime(2018, 9, 20, 12), 34.069444, -118.445278, -7)
    sunrise = solar_angle_result["sunrise"]
    sunset = solar_angle_result["sunset"]
    @test isapprox(Solar.partial_day(DateTime(2018, 9, 20, 6, 40)),
                   sunrise, atol=2. / 1440.)  # tolerance = 2 min
    @test isapprox(Solar.partial_day(DateTime(2018, 9, 20, 18, 52)),
                   sunset, atol=2. / 1440.)  # tolerance = 2 min
end

end  # module


text:
using DataFrames, NamedArrays
import DocOpt
using JLD

doc = """CSV Process Utility

Usage:
  process_csv.jl check <csvfile> [--exclude=<excludefile>]
  process_csv.jl statuscounts <csvfile> <csvfileout> [--exclude=<excludefile>] [--bestof=<list>]...
  process_csv.jl geomeans <csvfile> [--exclude=<excludefile>] [--bestof=<list>]...
  process_csv.jl perfprofile <csvfile> <output_jld> <solver>... [--exclude=<excludefile>] [--bestof=<list>]...
  process_csv.jl -h | --help

Options:
  -h --help               Show this screen.
  --exclude=<excludefile> A list of runs to exclude.
  --bestof=<list>         Create an artificial solver by taking the best of multiple solvers. Format: "A B C" creates a new solver A from best of B and C.
"""

arguments = DocOpt.docopt(doc)

# process the csv produced by process_output.jl
results = readtable(arguments["<csvfile>"])
println()

# file listing manually excluded results
# format by line: solver instance ignored comment
# example: PAJ_CPLEX_MOSEK clay0304h.cbf # CPLEX gave wrong answer
if isa(arguments["--exclude"], String)
    excluded = [split(strip(l))[1:2] for l in split(readstring(arguments["--exclude"]),'\n') if strip(l) != ""]
else
    excluded = []
end

function isexcluded(solver, inst)
    for k in 1:length(excluded)
        if excluded[k][1] == solver && excluded[k][2] == inst
            return true
        end
    end
    return false
end

# exclude any in the exclude file
for i in 1:size(results,1)
    if isexcluded(results[i, :solver], results[i, :instance])
        results[i, :newstatus] = "excl"
    end
end

for solverlist in arguments["--bestof"]
    newsolver = split(solverlist)[1]
    best_of_solvers = split(solverlist)[2:end]
    newrows = []
    for rowlist in groupby(results, :instance)
        besttime = Inf
        beststatus = "blank"
        inst = rowlist[1, :instance]
        for i in 1:size(rowlist, 1)
            if rowlist[i, :solver] in best_of_solvers
                if rowlist[i, :newstatus] == "conv"
                    besttime = min(besttime, rowlist[i, :total_time])
                    beststatus = "conv"
                elseif rowlist[i, :newstatus] == "lim"
                    if beststatus != "conv"
                        besttime = min(besttime, rowlist[i, :total_time])
                        beststatus = "lim"
                    end
                end
            end
        end
        newrow = Dict{Symbol,Any}(n => rowlist[1, n] for n in names(results))
        newrow[:instance] = inst
        newrow[:solver] = newsolver
        newrow[:newstatus] = beststatus
        newrow[:total_time] = besttime
        push!(newrows, newrow)
    end
    for r in newrows
        push!(results, r)
    end
end

conv_runs = results[results[:newstatus] .== "conv", :]

if arguments["check"]
    # check for disagreements in objective value
    println("Objective disagreements")
    for optval_by_instance in groupby(conv_runs, :instance)
        first_optval = optval_by_instance[1, :objval]
        sense = optval_by_instance[1, :sense]
        inst = optval_by_instance[1, :instance]
        for i in 2:size(optval_by_instance, 1)
            solver = optval_by_instance[i, :solver]
            optval = optval_by_instance[i, :objval]
            if abs(optval - first_optval)/(abs(first_optval) + 1e-5) > 1e-4
                println()
                println("$inst (sense = $sense)")
                show(optval_by_instance[:,[:solver, :objval, :objbound]])
                println()
                break
            end
        end
    end

    # # check for duplicated runs
    # for g1 in groupby(results, :instance)
    #     for g2 in groupby(g1, :solver)
    #         if size(g2,1) > 1
    #             s = g2[1, :solver]
    #             inst = g2[1, :instance]
    #             println("Multiple results for solver $s instance $inst")
    #         end
    #     end
    # end

elseif arguments["statuscounts"]
    # status counts by solver
    statuses = sort(collect(unique(results[:newstatus])))
    all_solvers = sort(collect(unique(results[:solver])))
    status_table = NamedArray(zeros(Int, length(all_solvers), length(statuses)+1))
    setnames!(status_table, all_solvers, 1)
    setnames!(status_table, [statuses; "Total"], 2)
    for i in 1:size(results,1)
        status_table[results[i, :solver], results[i, :newstatus]] += 1
        status_table[results[i, :solver], "Total"] += 1
    end

    df2 = DataFrame(convert(Array, status_table))
    names!(df2, convert(Vector{Symbol}, names(status_table,2)))
    df2[:solver] = names(status_table, 1)
    writetable(arguments["<csvfileout>"], df2)

elseif arguments["geomeans"]
    function shifted_geomean(table, field, shift, groupby)
        solvers = unique(table[groupby])
        sum_by = Dict(g => 0.0 for g in solvers)
        counter_by = Dict(g => 0 for g in solvers)
        # prod(x_i + s)^(1/n) - s
        # compute prod(x_i+s)^(1/n) as
        # exp(sum(log(x_i+s))/n)
        for i in 1:size(table, 1)
            if ismissing(table[i, field])
                continue
            end
            sum_by[table[i, groupby]] += log(table[i, field] + shift)
            counter_by[table[i, groupby]] += 1
        end
        for g in solvers
            r = exp(sum_by[g]/counter_by[g]) - shift
            if !isnan(r)
                @printf("%50s %f\n", g, r)
                # println("$g $r")
            end
        end
    end

    timeshift = 10.0
    itershift = 1
    nodeshift = 10

    # set any missing times to time limit
    for i in 1:size(results, 1)
        if ismissing(results[i, :total_time])
            results[i, :total_time] = results[i, :timelimit]
        else
            @assert !isnan(results[i, :total_time]) && results[i, :total_time] > 0
        end
    end

    all_solvers = unique(results[:solver])
    all_conv = []
    for g1 in groupby(conv_runs, :instance)
        hassolver = falses(length(all_solvers))
        for i in 1:size(g1, 1)
            solver = g1[i, :solver]
            idx = findnext(all_solvers, solver, 1)
            hassolver[idx] = true
        end
        if sum(hassolver) == length(all_solvers)
            push!(all_conv, g1[1, :instance])
        end
    end

    # subset of runs where all solvers converged
    runs_all_conv = conv_runs[find(t-> t in all_conv, conv_runs[:instance]), :]

    println("$(length(all_conv)) instances solved optimally by all solvers:\n")
    println(all_conv)
    println()

    println("\nShifted ($timeshift) geomean of solve time...")
    println("\n...on all (before termination)")
    shifted_geomean(results, :total_time, timeshift, :solver)
    println("\n...on conv by this solver")
    shifted_geomean(conv_runs, :total_time, timeshift, :solver)
    println("\n...on conv by all solvers")
    shifted_geomean(runs_all_conv, :total_time, timeshift, :solver)
    println()

    println("\nShifted ($itershift) geomean of conic subproblem count...")
    println("\n...on all (before termination)")
    shifted_geomean(results, :conic_count, itershift, :solver)
    println("\n...on conv by this solver")
    shifted_geomean(conv_runs, :conic_count, itershift, :solver)
    println("\n...on conv by all solvers")
    shifted_geomean(runs_all_conv, :conic_count, itershift, :solver)
    println()

    println("\nShifted ($itershift) geomean of iteration count...")
    println("\n...on all (before termination)")
    shifted_geomean(results, :iter_count, itershift, :solver)
    println("\n...on conv by this solver")
    shifted_geomean(conv_runs, :iter_count, itershift, :solver)
    println("\n...on conv by all solvers")
    shifted_geomean(runs_all_conv, :iter_count, itershift, :solver)
    println()

    println("\nShifted ($nodeshift) geomean of node count...")
    println("\n...on all (before termination)")
    shifted_geomean(results, :node_count, nodeshift, :solver)
    println("\n...on conv by this solver")
    shifted_geomean(conv_runs, :node_count, nodeshift, :solver)
    println("\n...on conv by all solvers")
    shifted_geomean(runs_all_conv, :node_count, nodeshift, :solver)
    println()

elseif arguments["perfprofile"]
    solvers = arguments["<solver>"]

    # need a table of [time] X [solver] where each row is a solver
    time_rows = []
    itndcount_rows = []
    instance_names = String[]
    for by_instance in groupby(conv_runs, :instance)
        time_row = fill(Inf, length(solvers))
        itndcount_row = fill(Inf, length(solvers))
        for i in 1:size(by_instance, 1)
            push!(instance_names, by_instance[1, :instance])
            if by_instance[i, :newstatus] == "conv"
                whichsolver = indexin([by_instance[i, :solver]], solvers)[1]
                if whichsolver != 0
                    time_row[whichsolver] = by_instance[i, :total_time]
                    if !ismissing(by_instance[i, :iter_count])
                        itndcount_row[whichsolver] = by_instance[i, :iter_count]
                    else
                        @assert !ismissing(by_instance[i, :node_count])
                        itndcount_row[whichsolver] = by_instance[i, :node_count]
                    end
                end # else not plotting this solver
            end # else not solved, Inf by default
        end
        push!(time_rows, time_row')
        push!(itndcount_rows, itndcount_row')
    end
    time_mat = vcat(time_rows...)
    itndcount_mat = vcat(itndcount_rows...)

    @assert endswith(arguments["<output_jld>"],".jld")
    save(arguments["<output_jld>"], "time_table", time_mat, "itndcount_table", itndcount_mat, "solvers", solvers, "instance_names", instance_names)
end
println()

max_stars_repo_name: chriscoey/PajaritoSupplement

max_issues_repo_name: chriscoey/PajaritoSupplement

max_forks_repo_name: chriscoey/PajaritoSupplement