text:
export run!, collect_agent_data!, collect_model_data!,
       init_agent_dataframe, init_model_dataframe, aggname,
       should_we_collect

###################################################
# Definition of the data collection API
###################################################
get_data(a, s::Symbol, obtainer::Function) = obtainer(getproperty(a, s))
get_data(a, f::Function, obtainer::Function) = obtainer(f(a))

get_data_missing(a, s::Symbol, obtainer::Function) =
    hasproperty(a, s) ? obtainer(getproperty(a, s)) : missing
function get_data_missing(a, f::Function, obtainer::Function)
    try
        obtainer(f(a))
    catch
        missing
    end
end

should_we_collect(s, model, when::AbstractVector) = s ∈ when
should_we_collect(s, model, when::Bool) = when
should_we_collect(s, model, when) = when(model, s)

"""
    run!(model, agent_step! [, model_step!], n::Integer; kwargs...) → agent_df, model_df
    run!(model, agent_step!, model_step!, n::Function; kwargs...) → agent_df, model_df

Run the model (step it with the input arguments propagated into [`step!`](@ref)) and collect
data specified by the keywords, explained one by one below. Return the data as
two `DataFrame`s, one for agent-level data and one for model-level data.

## Data-deciding keywords
* `adata::Vector` means "agent data to collect". If an entry is a `Symbol`, e.g. `:weight`,
  then the data for this entry is agent's field `weight`. If an entry is a `Function`, e.g.
  `f`, then the data for this entry is just `f(a)` for each agent `a`.
  The resulting dataframe columns are named with the input symbol (here `:weight, :f`).

* `adata::Vector{<:Tuple}`: if `adata` is a vector of tuples instead,
  data aggregation is done over the agent properties.

  For each 2-tuple, the first entry is the "key" (any entry like the ones mentioned above,
  e.g. `:weight, f`). The second entry is an aggregating function that aggregates the key,
  e.g. `mean, maximum`. So, continuing from the above example, we would have
  `adata = [(:weight, mean), (f, maximum)]`.

  It's also possible to provide a 3-tuple, with the third entry being a conditional
  function (returning a `Bool`), which assesses if each agent should be included in the
  aggregate. For example: `x_pos(a) = a.pos[1]>5` with `(:weight, mean, x_pos)` will result
  in the average weight of agents conditional on their x-position being greater than 5.

  The resulting data name columns use the function [`aggname`](@ref), and create something
  like `:mean_weight` or `:maximum_f_x_pos`.
  This name doesn't play well with anonymous functions, but you can simply use
  `DataFrames.rename!` to change the returned dataframe's column names.

  **Notice:** Aggregating only works if there are agents to be aggregated over.
  If you remove agents during model run, you should modify the aggregating functions.
  *E.g.* instead of passing `mean`, pass `mymean(a) = isempty(a) ? 0.0 : mean(a)`.

* `mdata::Vector` means "model data to collect" and works exactly like `adata`.
  For the model, no aggregation is possible (nothing to aggregate over).

By default both keywords are `nothing`, i.e. nothing is collected/aggregated.

## Mixed-Models

  For mixed-models, the `adata` keyword has some additional options & properties.
  An additional column `agent_type` will be placed in the output
  dataframe.

  In the case that data is needed for one agent type that does not exist
  in a second agent type, `missing` values will be added to the dataframe.

  **Warning:** Since this option is inherently type unstable, try to avoid this
  in a performance critical situation.

  Aggregate functions will fail if `missing` values are not handled explicitly.
  If `a1.weight` but `a2` (type: Agent2) has no `weight`, use
  `a2(a) = a isa Agent2; adata = [(:weight, sum, a2)]` to filter out the missing results.

### Other keywords
* `when=true` : at which steps `s` to perform the data collection and processing.
  A lot of flexibility is offered based on the type of `when`. If `when::Vector`,
  then data are collect if `s ∈ when`. Otherwise data are collected if `when(model, s)`
  returns `true`. By default data are collected in every step.
* `when_model = when` : same as `when` but for model data.
* `obtainer = identity` : method to transfer collected data to the `DataFrame`.
  Typically only change this to [`copy`](https://docs.julialang.org/en/v1/base/base/#Base.copy)
  if some data are mutable containers (e.g. `Vector`) which change during evolution,
  or [`deepcopy`](https://docs.julialang.org/en/v1/base/base/#Base.deepcopy) if some data are
  nested mutable containers.
  Both of these options have performance penalties.
* `replicates=0` : Run `replicates` replicates of the simulation.
* `parallel=false` : Only when `replicates>0`. Run replicate simulations in parallel.
* `agents_first=true` : Whether to update agents first and then the model, or vice versa.
"""
function run! end

run!(model::ABM, agent_step!, n::Int = 1; kwargs...) =
run!(model::ABM, agent_step!, dummystep, n; kwargs...)

function run!(model::ABM, agent_step!, model_step!, n;
  replicates::Int=0, parallel::Bool=false, kwargs...)

  r = replicates
  if r > 0
    if parallel
      return parallel_replicates(model, agent_step!, model_step!, n, r; kwargs...)
    else
      return series_replicates(model, agent_step!, model_step!, n, r; kwargs...)
    end
  else
    return _run!(model, agent_step!, model_step!, n; kwargs...)
  end
end

###################################################
# Core data collection loop
###################################################
"""
  _run!(model, agent_step!, model_step!, n; kwargs...)
Core function that loops over stepping a model and collecting data at each step.
"""
function _run!(model, agent_step!, model_step!, n;
               when = true, when_model = when,
               mdata=nothing, adata=nothing, obtainer = identity,
               agents_first=true)

    df_agent = init_agent_dataframe(model, adata)
    df_model = init_model_dataframe(model, mdata)
    if n isa Integer
        if when == true; for c in eachcol(df_agent); sizehint!(c, n); end; end
        if when_model == true; for c in eachcol(df_model); sizehint!(c, n); end; end
    end

    s = 0
    while until(s, n, model)
        if should_we_collect(s, model, when)
            collect_agent_data!(df_agent, model, adata, s; obtainer)
        end
        if should_we_collect(s, model, when_model)
            collect_model_data!(df_model, model, mdata, s; obtainer)
        end
        step!(model, agent_step!, model_step!, 1, agents_first)
        s += 1
    end
    if should_we_collect(s, model, when)
        collect_agent_data!(df_agent, model, adata, s; obtainer)
    end
    if should_we_collect(s, model, when_model)
        collect_model_data!(df_model, model, mdata, s; obtainer)
    end
    return df_agent, df_model
end

###################################################
# core data collection functions per step
###################################################
"""
    init_agent_dataframe(model, adata) → agent_df
Initialize a dataframe to add data later with [`collect_agent_data!`](@ref).
"""
init_agent_dataframe(model, properties::Nothing) = DataFrame()

"""
    collect_agent_data!(df, model, properties, step = 0; obtainer = identity)
Collect and add agent data into `df` (see [`run!`](@ref) for the dispatch rules
of `properties` and `obtainer`). `step` is given because the step number information
is not known.
"""
collect_agent_data!(df, model, properties::Nothing, step::Int=0; kwargs...) = df

function init_agent_dataframe(model::ABM{S,A}, properties::AbstractArray) where {S,A<:AbstractAgent}
    nagents(model) < 1 && throw(ArgumentError(
        "Model must have at least one agent to initialize data collection",
    ))

    utypes = union_types(A)
    std_headers = length(utypes) > 1 ? 3 : 2

    headers = Vector{Symbol}(undef, std_headers + length(properties))
    headers[1] = :step
    headers[2] = :id

    for i in 1:length(properties)
        headers[i+std_headers] = Symbol(properties[i])
    end

    types = Vector{Vector}(undef, std_headers + length(properties))
    types[1] = Int[]
    types[2] = Int[]

    if std_headers == 3
        headers[3] = :agent_type
        multi_agent_types!(types, utypes, model, properties)
    else
        single_agent_types!(types, model, properties)
    end

    DataFrame(types, headers)
end

function single_agent_types!(types::Vector{Vector{T} where T}, model::ABM, properties::AbstractArray)
    a = random_agent(model)
    for (i, k) in enumerate(properties)
        current_type = typeof(get_data(a, k, identity))
        isconcretetype(current_type) || warn(
            "Type is not concrete when using $(k) " *
            "on agents. Consider narrowning the type signature of $(k).",
        )
        types[i+2] = current_type[]
    end
end

function multi_agent_types!(types::Vector{Vector{T} where T}, utypes::Tuple, model::ABM, properties::AbstractArray)
    types[3] = Symbol[]

    for (i, k) in enumerate(properties)
        current_types = DataType[]
        for atype in utypes
            a = first(Iterators.filter(a -> a isa atype, allagents(model)))
            if k isa Symbol
                current_type =
                    hasproperty(a, k) ? typeof(get_data(a, k, identity)) : Missing
            else
                current_type = try
                    typeof(get_data(a, k, identity))
                catch
                    Missing
                end
            end
            isconcretetype(current_type) || warn(
                "Type is not concrete when using $(k) " *
                "on $(atype) agents. Consider narrowning the type signature of $(k).",
            )
            push!(current_types, current_type)
        end
        unique!(current_types)
        if length(current_types) == 1
            current_types[1] <: Missing && error("$(k) does not yield a valid agent property.")
            types[i+3] = current_types[1][]
        else
            types[i+3] = Union{current_types...}[]
        end
    end
end

function collect_agent_data!(df, model, properties::Vector, step::Int=0; kwargs...)
    alla = sort!(collect(values(model.agents)), by = a -> a.id)
    dd = DataFrame()
    dd[!, :step] = fill(step, length(alla))
    dd[!, :id] = map(a -> a.id, alla)
    if :agent_type ∈ propertynames(df)
        dd[!, :agent_type] = map(a -> Symbol(typeof(a)), alla)
    end

    for fn in properties
        _add_col_data!(dd, eltype(df[!, Symbol(fn)]), fn, alla; kwargs...)
    end
    append!(df, dd)
    return df
end

function _add_col_data!(dd::DataFrame, col::Type{T}, property, agent_iter; obtainer = identity) where {T}
    dd[!, Symbol(property)] = collect(get_data(a, property, obtainer) for a in agent_iter)
end

function _add_col_data!(dd::DataFrame, col::Type{T}, property, agent_iter; obtainer = identity) where {T>:Missing}
    dd[!, Symbol(property)] = collect(get_data_missing(a, property, obtainer) for a in agent_iter)
end

# Aggregating version
function init_agent_dataframe(model::ABM{S,A}, properties::Vector{<:Tuple}) where {S,A<:AbstractAgent}
    nagents(model) < 1 && throw(ArgumentError(
        "Model must have at least one agent to " * "initialize data collection",
    ))
    headers = Vector{String}(undef, 1 + length(properties))
    types = Vector{Vector}(undef, 1 + length(properties))

    utypes = union_types(A)

    headers[1] = "step"
    types[1] = Int[]

    if length(utypes) > 1
        multi_agent_agg_types!(types, utypes, headers, model, properties)
    else
        single_agent_agg_types!(types, headers, model, properties)
    end
    DataFrame(types, headers)
end

function single_agent_agg_types!(types::Vector{Vector{T} where T}, headers::Vector{String}, model::ABM, properties::AbstractArray)
    for (i, property) in enumerate(properties)
        k, agg = property
        headers[i+1] = aggname(property)
        # This line assumes that `agg` can work with iterators directly
        current_type =
            typeof(agg(get_data(a, k, identity) for a in Iterators.take(allagents(model), 1)))
        isconcretetype(current_type) || warn(
            "Type is not concrete when using function $(agg) " *
            "on key $(k). Consider using type annotation, e.g. $(agg)(a)::Float64 = ...",
        )
        types[i+1] = current_type[]
    end
end

function multi_agent_agg_types!(types::Vector{Vector{T} where T}, utypes::Tuple, headers::Vector{String}, model::ABM, properties::AbstractArray)
    for (i, property) in enumerate(properties)
        k, agg = property
        headers[i+1] = aggname(property)
        current_types = DataType[]
        for atype in utypes
            a = first(Iterators.filter(a -> a isa atype, allagents(model)))
            if k isa Symbol
                current_type =
                    hasproperty(a, k) ? typeof(agg(get_data(a, k, identity))) : Missing
            else
                current_type = try
                    typeof(agg(get_data(a, k, identity)))
                catch
                    Missing
                end
            end
            isconcretetype(current_type) || warn(
                "Type is not concrete when using function $(agg) " *
                "on key $(k) for $(atype) agents. Consider using type annotation, e.g. $(agg)(a)::Float64 = ...",
            )
            push!(current_types, current_type)
        end
        unique!(current_types)
        filter!(t -> !(t <: Missing), current_types) # Ignore missings here
        if length(current_types) == 1
            types[i+1] = current_types[1][]
        elseif length(current_types) > 1
            error("Multiple types found for aggregate function $(agg) on key $(k).")
        else
            error("No possible aggregation for $(k) using $(agg)")
        end
    end
end

"""
    aggname(k) → name
    aggname(k, agg) → name
    aggname(k, agg, condition) → name

Return the name of the column of the `i`-th collected data where `k = adata[i]`
(or `mdata[i]`).
`aggname` also accepts tuples with aggregate and conditional values.
"""
aggname(k, agg) = string(agg)*"_"*string(k)
aggname(k, agg, condition) = string(agg)*"_"*string(k)*"_"*string(condition)
aggname(x::Tuple{K,A}) where {K,A} = aggname(x[1], x[2])
aggname(x::Tuple{K,A,C}) where {K,A,C} = aggname(x[1], x[2], x[3])
aggname(x::Union{Function, Symbol, String}) = string(x)

function collect_agent_data!(df, model::ABM, properties::Vector{<:Tuple}, step::Int=0; kwargs...)
    alla = allagents(model)
    push!(df[!, 1], step)
    for (i, prop) in enumerate(properties)
        _add_col_data!(df[!, i+1], prop, alla; kwargs...)
    end
    return df
end

# Normal aggregates
function _add_col_data!(col::AbstractVector{T}, property::Tuple{K,A}, agent_iter; obtainer = identity) where {T,K,A}
    k, agg = property
    res::T = agg(get_data(a, k, obtainer) for a in agent_iter)
    push!(col, res)
end

# Conditional aggregates
function _add_col_data!(col::AbstractVector{T}, property::Tuple{K,A,C}, agent_iter; obtainer = identity) where {T,K,A,C}
    k, agg, condition = property
    res::T = agg(get_data(a, k, obtainer) for a in Iterators.filter(condition, agent_iter))
    push!(col, res)
end

# Model data
"""
    init_model_dataframe(model, mdata) → model_df
Initialize a dataframe to add data later with [`collect_model_data!`](@ref).
"""
function init_model_dataframe(model::ABM, properties::Vector)
    headers = Vector{Symbol}(undef, 1+length(properties))
    headers[1] = :step
    for i in 1:length(properties); headers[i+1] = Symbol(properties[i]); end

    types = Vector{Vector}(undef, 1+length(properties))
    types[1] = Int[]
    for (i,k) in enumerate(properties)
        types[i+1] =
            if typeof(k) <: Symbol
                typeof(model.properties[k])[]
            else
                current_type = typeof(k(model))
                isconcretetype(current_type) || warn("Type is not concrete when using $(k)"*
                "on the model. Considering narrowing the type signature of $(k).")
                current_type[]
            end
    end
    DataFrame(types, headers)
end

init_model_dataframe(model::ABM, properties::Nothing) = DataFrame()

"""
    collect_model_data!(df, model, properties, step = 0, obtainer = identity)
Same as [`collect_agent_data!`](@ref) but for model data instead.
"""
function collect_model_data!(df, model, properties::Vector, step::Int=0; obtainer = identity)
  push!(df[!, :step], step)
  for fn in properties
    push!(df[!, Symbol(fn)], get_data(model, fn, obtainer))
  end
  return df
end

collect_model_data!(df, model, properties::Nothing, step::Int=0; kwargs...) = df

###################################################
# Parallel / replicates
###################################################
function replicate_col!(df, rep)
  df[!, :replicate] = [rep for i in 1:size(df, 1)]
end

"Run replicates of the same simulation"
function series_replicates(model, agent_step!, model_step!, n, replicates; kwargs...)

  df_agent, df_model = _run!(deepcopy(model), agent_step!, model_step!, n; kwargs...)
  replicate_col!(df_agent, 1)
  replicate_col!(df_model, 1)

  for rep in 2:replicates
    df_agentTemp, df_modelTemp = _run!(deepcopy(model), agent_step!, model_step!, n; kwargs...)
    replicate_col!(df_agentTemp, rep)
    replicate_col!(df_modelTemp, rep)

    append!(df_agent, df_agentTemp)
    append!(df_model, df_modelTemp)
  end
  return df_agent, df_model
end

"Run replicates of the same simulation in parallel"
function parallel_replicates(model::ABM, agent_step!, model_step!, n, replicates; kwargs...)

  all_data = pmap(j -> _run!(deepcopy(model), agent_step!, model_step!, n; kwargs...),
                  1:replicates)

  df_agent = DataFrame()
  df_model = DataFrame()
  for (rep, d) in enumerate(all_data)
    replicate_col!(d[1], rep)
    replicate_col!(d[2], rep)
    append!(df_agent, d[1])
    append!(df_model, d[2])
  end

  return df_agent, df_model
end


text:
using ArgParse
using util
using ranking
using preprocessing
using Method
using data
using Metrics
using preprocessing
using JSON
using DataStructures

# Include run_method file
dir_name, script_name = splitdir(@__FILE__)
include(joinpath(dir_name,"run_method.jl"))

s = ArgParseSettings()
@add_arg_table s begin
    "bug_dataset"
        help="File that contains all bug report data(categorical, summary,...)"
    "evaluation"
        help="Evaluation dataset"
    "method_name"
    "search_space_script"
    "--max_evals"
            help = ""
            arg_type = Int
            default = 100
    "--test"
        help="Test dataset"
end

import_settings(s, get_basic_args())

parsed_args = parse_args(s, as_symbols=true)

@info parsed_args

# Load validation set. This set is used to find the best parameters.
validationDataset = readDataset(parsed_args[:evaluation])

# Load search space for the experiment
using PyCall

hyperopt = pyimport("hyperopt")
space_script_path = parsed_args[:search_space_script]

py"""
space_script = $(space_script_path)
f = open(space_script, 'r')
code_str = f.read()
f.close()

exec(compile(code_str, space_script, 'exec'))
"""

fixed_values = py"""fixed_values"""
space = py"""space"""

# Objective Function
function objective(params...)
    exp_parameter = Dict(parsed_args)

    for (k,v) in  params[1]
        exp_parameter[Symbol(k)] = v
    end

    for (k,v) in  fixed_values
        exp_parameter[Symbol(k)] = v
    end

    results = Dict()

    timeScore = @elapsed begin
        try
            results = run_method(validationDataset.bugIds, exp_parameter)
        catch e
            @error exception=(e, catch_backtrace())
            throw(e)
        end
    end

    map = 0.0
    auc = 0.0
    rr1 = 0.0
    rr5 = 0.0

    print(typeof(results))
    for r in results
        if r["label"] == "MAP_RecallRate"
            map = r["map"]::Float64
            rr1 = r["rr"][1]::Float64
            rr5 = r["rr"][5]::Float64
        elseif r["label"] == "BinaryPredictionROC"
            auc = r["auc"]::Float64
        end
    end

    loss = 2.0 - (map + auc)

    new_trial = OrderedDict(:x=> params[1],:loss => loss, :map=> map, :auc=>auc,:rr1 => rr1, :rr5=>rr5 , :args => exp_parameter)
    @info new_trial
    @info "Finish a new trial in $(timeScore)"

    return Dict("loss"=> 2.0 - (map + auc), "status"=> hyperopt.STATUS_OK, "results"=> results, "map"=> map, "auc"=> auc,
                "args"=> exp_parameter)
end

# Find the best parameters
trials = hyperopt.Trials()
best = hyperopt.fmin(objective,
            space=space,
            algo=hyperopt.tpe.suggest,
            max_evals=parsed_args[:max_evals],
            trials=trials)

@info "Paremeter search is over!"
@info "Best: $(best)"

@info "Sorted trial:"

all_trials = [ OrderedDict{Symbol, Any}(:vals => t["misc"]["vals"], :loss => t["result"]["loss"], :map => t["result"]["map"], :auc => t["result"]["auc"], :result => t["result"]["results"], :args => t["result"]["args"]) for t in trials]

sort!(all_trials, by=x-> x[:loss])

for t in all_trials
    @info "$(t)"
end

best_result = all_trials[1]
println(best_result)

if !isnothing(get(parsed_args, :test, nothing))
    @info "Start test:"
    testDataset = readDataset(parsed_args[:test])
    best_args = Dict{Symbol, Any}(( Symbol(k) => v for (k, v) in best_result[:args]))

    @info "$(best_args)"

    results = run_method(testDataset.bugIds, best_args)

    @info  JSON.json(results)
end

ext: jl

max_stars_repo_name: irving-muller/FaST

max_issues_repo_name: irving-muller/FaST

max_forks_repo_name: irving-muller/FaST


text:
"""
State machine for these tests:

@startuml
state TestDeepHistoryMachine {
    state TestDeepHistoryStartState
    state TestDeepHistoryL1State {
        state TestDeepHistoryL2State {
            state TestDeepHistoryL3State
        }
    }
}
@enduml
"""

using Logging, Test, HSM
include("state_transition_test_utils.jl")

@transition_tests_state TestDeepHistoryMachine
@transition_tests_state TestDeepHistoryStartState
@transition_tests_state TestDeepHistoryL1State
@transition_tests_state TestDeepHistoryL2State
@transition_tests_state TestDeepHistoryL3State

function build_deep_history_machine()
    machine = TestDeepHistoryMachine(nothing)
    start_state = TestDeepHistoryStartState(machine)
    l1_state = TestDeepHistoryL1State(machine)
    l2_state = TestDeepHistoryL2State(l1_state)
    l3_state = TestDeepHistoryL3State(l2_state)

    HSM.transition_to_state!(machine, start_state)

    return machine, start_state, l1_state, l2_state, l3_state
end

#
# Deep history transition to a state that has not been entered before is the
# transitioned-to state.
#
@testset "transition_to_deep_history!() -- start -deep-> l1 => l1" begin
    machine, start_state, l1_state, l2_state, l3_state = build_deep_history_machine()
    HSM.transition_to_deep_history!(machine, l1_state)
    @test HSM.active_state(machine) == l1_state
end

#
# Deep history transition to a state whose substate's child has been entered
# before is the substate.
#
@testset "transition_to_deep_history!() -- start -> l3 -> start -deep-> l1 => l3" begin
    machine, start_state, l1_state, l2_state, l3_state = build_deep_history_machine()
    HSM.transition_to_state!(machine, l3_state)
    HSM.transition_to_state!(machine, start_state)
    HSM.transition_to_deep_history!(machine, l1_state)
    @test HSM.active_state(machine) == l3_state
end


text:
# This file is a part of BitOperations.jl, licensed under the MIT License (MIT).

using BitOperations
using Test

@testset "zigzag encoding" begin
    @testset "zigzagenc" begin
        @test zigzagenc(0) == 0x0000000000000000
        @test zigzagenc(-1) == 0x0000000000000001
        @test zigzagenc(1) == 0x0000000000000002
        @test zigzagenc(- signed(0x8000000000000000)) == 0xffffffffffffffff
        @test zigzagenc(+ signed(0x7fffffffffffffff)) == 0xfffffffffffffffe
    end

    @testset "zigzagdec" begin
        @test zigzagdec(0x0000000000000000) == 0
        @test zigzagdec(0x0000000000000001) == -1
        @test zigzagdec(0x0000000000000002) == 1
        @test zigzagdec(0xffffffffffffffff) == - signed(0x8000000000000000)
        @test zigzagdec(0xfffffffffffffffe) == + signed(0x7fffffffffffffff)
    end
end


text:
module Finance
    export fortune

    function fortune(f0, p, c0, n, i)::Bool
        addpercent(current,percent) = current+(current*percent)÷100
        f = f0
        c = c0
        for year in 2:n
            f = addpercent(f,p)-c
            c = addpercent(c,i)
        end
        f ≥ 0
    end

end

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
# Main module for BayesianMixtures package
module BayesianMixtures

include("MFM.jl")
include("RandomNumbers.jl")

include("Normal.jl")
include("MVN.jl")
include("MVNaaC.jl")
include("MVNaaN.jl")
include("MVNaaRJ.jl")
include("NormalNonoptimized.jl")

using Random
using Statistics
using SpecialFunctions
lgamma_(x) = logabsgamma(x)[1]

# ===================================================================
# ===================================================================
# ================== Functions to generate results ==================
# ===================================================================
# ===================================================================

# Create an options object to specify model, data, and MCMC parameters.
function options(
        mode, # "Normal", "MVN", "MVNaaC", "MVNaaN", or "MVNaaRJ"
        model_type, # "MFM" or "DPM"
        x, # data
        n_total; # total number of MCMC sweeps to run the sampler
        n_keep=n_total, # number of MCMC sweeps to keep after thinning
        n_burn=round(Int,n_total/10), # number of MCMC sweeps (out of n_total) to discard as burn-in
        verbose=true, # display information or not
        use_hyperprior=true, # put hyperprior on the base distn parameters or not
        t_max=40, # a guess at an upper bound on # of clusters that will be encountered during MCMC

        # MFM options:
        gamma=1.0, # Dirichlet_k(gamma,...,gamma)
        log_pk="k -> log(0.1)+(k-1)*log(0.9)", # string representation of log(p(k))
            # (where p(k) is the log of the prior on # of components, K)

        # DPM options:
        alpha_random=true, # put prior on alpha (DPM concentration parameter) or not
        alpha=1.0, # value of alpha (initial value if alpha_random=true)

        # Jain-Neal split-merge options:
        use_splitmerge=true, # use split-merge or not
        n_split=5, # number of intermediate sweeps for split launch state
        n_merge=5,  #                 "         "       merge    "     "

        # RJMCMC options:
        k_max=t_max # a guess at an upper bound on # of components that will be encountered during MCMC
    )

    # Compute partition distribution values
    n = length(x)
    if model_type=="MFM"
        lpk = eval(Meta.parse(log_pk))
        log_pk_fn(k) = Base.invokelatest(lpk,k)
        log_v = MFM.coefficients(log_pk_fn,gamma,n,t_max+1)
        a = b = gamma
    elseif model_type=="DPM"
        log_v = float(1:t_max+1)*log(alpha) .- lgamma_(alpha+n) .+ lgamma_(alpha)
        a,b = 1.,0.
    else
        error("Invalid model_type: $model_type.")
    end
    if mode=="MVNaaRJ"; @assert(model_type=="MFM", "RJMCMC is not implemented for DPMs."); end

    n_keep = min(n_keep,n_total)
    module_ = getfield(BayesianMixtures,Symbol(mode))
    return module_.Options(mode, model_type, x, n_total, n_keep, n_burn, verbose,
                           use_hyperprior, t_max, gamma, log_pk, alpha_random, alpha,
                           use_splitmerge, n_split, n_merge, k_max, a, b, log_v, n)
end

# Run the MCMC sampler with the specified options.
function run_sampler(options)
    o = options
    n,n_total,n_keep = o.n,o.n_total,o.n_keep
    module_ = getfield(BayesianMixtures,Symbol(o.mode))

    # Short run to precompile (calling precompile doesn't seem to work...)
    module_.sampler(o,1,1)

    if o.verbose
        println(o.mode, " ", o.model_type)
        println("n = $n, n_total = $n_total, n_keep = $n_keep")
        print("Running... ")
    end

    # Main run
    elapsed_time = (@elapsed t_r,N_r,z_r,theta_r,keepers = module_.sampler(o,n_total,n_keep))
    time_per_step = elapsed_time/(n_total*n)

    if o.verbose
        println("complete.")
        println("Elapsed time = $elapsed_time seconds")
        println("Time per step ~ $time_per_step seconds")
    end

    # # Run profiler
    # @profile sampler(x,n_total,n_keep,o)
    # Profile.print(format = :flat)
    # Profile.clear()

    return module_.Result(o,t_r,N_r,z_r,theta_r,keepers,elapsed_time,time_per_step)
end

# ===================================================================
# ===================================================================
# ================== Functions to analyze results ===================
# ===================================================================
# ===================================================================

# Compute running average of posterior CDF on t (# of clusters)
function t_running(result)
    o = result.options
    n_total = o.n_total
    T = zeros(n_total,o.t_max)
    for i = 1:n_total; T[i,result.t[i]] = 1; end
    cdfs = cumsum(cumsum(T,dims=1),dims=2)./(1:n_total)
    return cdfs
end

# Compute histogram with the specified bin edges,
# where x[i] is in bin j if edges[j] < x[i] <= edges[j+1].
function histogram(x, edges=[]; n_bins=50, weights=ones(length(x)))
    if isempty(edges)
        mn,mx = minimum(x),maximum(x)
        r = mx-mn
        edges = range(mn-r/n_bins, stop=mx+r/n_bins, length=n_bins)
    else
        n_bins = length(edges)-1
    end
    counts = zeros(Float64,n_bins)
    for i=1:length(x)
        for j=1:n_bins
            if (edges[j] < x[i] <= edges[j+1])
                counts[j] += weights[i]
                break
            end
        end
    end
    return counts,edges
end

# Compute posterior on t (# of clusters)
function t_posterior(result)
    o = result.options
    use = o.n_burn+1:o.n_total
    counts,edges = histogram(result.t[use], 0:o.t_max) # bins are: (0,1],(1,2],...,(t_max-1,t_max]
    return counts/length(use)
end

# Compute posterior on k (# of components)
function k_posterior(result; upto=result.options.t_max)
    o = result.options
    @assert(o.model_type=="MFM", "The posterior on k is not defined for the DPM.")
    lpk = eval(Meta.parse(o.log_pk))
    log_pk_fn(k) = Base.invokelatest(lpk,k)
    p_kt = MFM.p_kt(log_pk_fn,o.gamma,o.n,upto,o.t_max)
    return p_kt*t_posterior(result)
end

# Compute density estimate at a point x (using Green's pseudo-predictive density).
function density_estimate(x,result)
    o = result.options
    n,b = o.n,o.b
    t,N,theta,keepers = result.t, result.N, result.theta, result.keepers
    n_keep = length(keepers)
    loglik = getfield(BayesianMixtures,Symbol(o.mode)).log_likelihood
    r = 0.0
    n_used = 0
    for ik = 1:n_keep
        i = keepers[ik]
        if i>o.n_burn
            for c in findall(N[:,i].>0)
                r += exp(loglik(x,theta[c,ik])) * (N[c,i]+b)/(n + b*t[i])
            end
            n_used += 1
        end
    end
    return r/n_used
end

rand_categorical(p) = (s=rand(); i=0; while s>0; i+=1; s-=p[i]; end; i)
rand_dirichlet(a) = (w=RandomNumbers.gamma.(a,1); w/sum(w))
counts(a,k) = (N=zeros(Int,k); for ai in a; N[ai]+=1; end; N)
logsumexp(x) = (m = maximum(x); m == -Inf ? -Inf : log.(sum(exp.(x.-m))) + m)

# Convert MFM result into posterior samples of k,w,a,theta,
# where k = number of components
#       w = mixture weights
#       a = assignments of data points to mixture components
#       theta = mixture component parameters
function sample_mixture_parameters(result,kmax)
    r,o = result,result.options
    @assert(o.model_type=="MFM", "This function is not defined for the DPM.")
    lpk = eval(Meta.parse(o.log_pk))
    log_pk_fn(k) = Base.invokelatest(lpk,k)
    p_kt = MFM.p_kt(log_pk_fn,o.gamma,o.n,kmax,o.t_max)

    module_ = getfield(BayesianMixtures,Symbol(o.mode))
    Theta = module_.Theta
    H = module_.construct_hyperparameters(o)

    k = zeros(Int16,o.n_keep)
    a = zeros(Int16,o.n,o.n_keep)
    w = Array{Float64,1}[]
    theta = Array{Theta,1}[]
    for i = 1:o.n_keep
        # sample k from p(k|t)
        ikeep = r.keepers[i]
        k[i] = rand_categorical(p_kt[:,r.t[ikeep]])

        # sample a from p(a|k,C) where C = clustering
        uz = sort(unique(r.z[:,i]))
        assign = zeros(Int,o.t_max+3)
        assign[uz] = randperm(k[i])[1:length(uz)]
        a[:,i] = assign[r.z[:,i]]

        # sample w from p(w|k,a)
        push!(w, rand_dirichlet(o.gamma .+ counts(a[:,i],k[i])))

        # sample from p(theta|k,a,phi)
        theta_i = [module_.new_theta(H)::Theta for j=1:k[i]]
        for j = 1:k[i]
            c_j = findfirst(isequal(j),assign)
            if c_j==nothing
                module_.prior_sample!(theta_i[j],H)
            else
                theta_i[j] = deepcopy(r.theta[c_j,i])
            end
        end
        push!(theta,theta_i)
    end
    return k,a,w,theta
end

# Compute log-likelihood of the data for a sequence of mixture weights w and component parameters theta.
function log_likelihoods(x,w,theta,mode)
    module_ = getfield(BayesianMixtures,Symbol(mode))
    log_f = module_.log_likelihood
    m = length(w)
    n = length(x)
    loglik = zeros(m)
    for iter = 1:m
        k = length(w[iter])
        ll = 0.0
        for i = 1:n
            ll += logsumexp(Float64[log(w[iter][j]) + log_f(x[i],theta[iter][j]) for j = 1:k])
        end
        loglik[iter] = ll
    end
    return loglik
end

# Compute the posterior similarity matrix (probability that i and j are in same cluster).
function similarity_matrix(result)
    o = result.options
    n = o.n
    z = result.z
    C = zeros(Int,n,n)
    n_used = 0
    for index = 1:o.n_keep
        if result.keepers[index]>o.n_burn
            for i = 1:n, j = 1:n
                if z[i,index]==z[j,index]; C[i,j] += 1; end
            end
            n_used += 1
        end
    end
    return C/n_used
end

# Compute cross-correlation between two vectors
function xcorr_(x,y)
    m,n = length(x),length(y)
    z = zeros(m+n-1)
    for k = 1:(m+n-1)
        s = 0.0
        for i = max(1,k-n+1):min(k,m)
            s += x[i]*y[n-k+i]
        end
        z[k] = s
    end
    return z
end

# Compute autocorrelation curve for the sequence y=[y_1,...,y_T]
function autocorrelation(y)
    T = length(y)
    u = ones(typeof(y[1]),T)
    d = [1:T; T-1:-1:1]
    autocov = xcorr_(y,y)./d - (xcorr_(y,u)./d).*(xcorr_(u,y)./d)
    a = autocov[T:-1:1]
    return a/a[1]
end
#autocorrelation_time(auto,c) = (tau=1+2*cumsum(auto[2:end]); M=findfirst([1:length(tau)] .>= c*tau); tau[M])

end # module BayesianMixtures


text:
module Lorenz63

export L63

include("rhs.jl")
include("time_integration.jl")
include("L63.jl")

end


text:
using GrupoBeta
using Test

@testset "Algoritmo_QR" begin
  E = [1 1 1; 1 2 1; 1 1 3]
  @test isapprox(autovalores(E), [4.214319 0.0 0.0; 0.0 1.460811 0.0; 0.0 0.0 0.324869], rtol = 1e-4)
  @test isapprox(autovetores(E), [0.39711 0.23319 -0.88765; 0.52065 0.73923 0.42713; 0.75578 -0.63178 0.17214], rtol = 1e-4)
  @test isapprox(iteracoes_QR(E), 24, rtol = 1e-1)
end

@testset "Metodo_potencias" begin
  A = [1 2; 3 2]
  @test isapprox(autovalor_dominante(A), 4, rtol = 1e-4)
  @test isapprox(autovetor_dominante(A), [0.55470; 0.83204], rtol = 1e-4)
  @test isapprox(iteracoes_potencias(A), 10, rtol = 1e-1)
end

max_stars_repo_name: EleonoraAvello/GrupoBetaPacote

max_issues_repo_name: EleonoraAvello/GrupoBetaPacote

max_forks_repo_name: EleonoraAvello/GrupoBetaPacote


text:
using Statistics
using DataStructures

function p1(lines)
  match = Dict('('=>')','['=>']','{'=>'}','<'=>'>')
  scores = Dict(')'=>3, ']'=>57, '}'=>1197, '>'=>25137)
  t = 0
  for line in lines
    Stack = []
    for c in line
      if c ∈ keys(match)
        push!(Stack, c)
      else
        last = pop!(Stack)
        if match[last] != c
          t += scores[c]
          break
        end
      end
    end
  end
  t
end

function p2(lines)
  match = Dict('('=>')','['=>']','{'=>'}','<'=>'>')
  score = Dict('('=>1, '['=>2, '{'=>3, '<'=>4)
  scores = []
  for line in lines
    corrupt = false
    stack = []
    for c in line
      if c ∈ keys(match)
        push!(stack, c)
      elseif match[pop!(stack)] != c
        corrupt = true
        break
      end
    end
    push!(scores, foldl((x,y) -> 5x+score[y], stack, init=0))
  end
  trunc(Int,median(scores))
end

const match = Dict('('=>')','['=>']','{'=>'}','<'=>'>')
const score = Dict('('=>1, '['=>2, '{'=>3, '<'=>4)

function p2helper(line::String)
  corrupt = false
  stack = Stack{Char}()
  for c in line
    if c ∈ keys(match)
      push!(stack, c)
    elseif match[pop!(stack)] != c
      corrupt = true
      break
    end
  end
  s = -1
  if !corrupt
    s = 0
    for x in stack
      s = 5s + score[x]
    end
  end
  s
end

# best speed singlethread
function p2fast(lines::Vector{String})
  scores = Int64[]
  for line in lines
    s = p2helper(line)
    s != -1 && push!(scores, s)
  end
  sort(scores)[trunc(Int,length(lines)/2)]
end

# best speed multithread; works best at nthreads = 3 on my machine
function p2fastMulti(lines::Vector{String})
  scores = Int64[]
  lk = ReentrantLock()
  Threads.@threads for line in lines
    s = p2helper(line)
    s != -1 && lock(lk) do
      push!(scores, s)
    end
  end
  sort(scores)[trunc(Int,length(lines)/2)]
end

using BenchmarkTools

lines = readlines("inputs/6-10/10.txt")

@btime p2fastMulti(lines)


text:
using Fermi
using Test

tol = 1E-8
#include("test_Options.jl")
@testset "Fermi" begin
    include("test_RHF.jl")
    include("test_CC.jl")
    include("test_MP.jl")
    include("test_Options.jl")
    include("test_CI.jl")
end


text:
"""
    PostImputingMatrix{T <: Number, U <: AbstractMatrix{T}, V <: AbstractVector{T}} <: AbstractMatrix{T}

A parametrized matrix that fills in a default vector of parameters whenever a \"missing\" column
is encountered during multiplication.

Supports multiplication with [`NGramMatrix`](@ref), [`MaybeHotMatrix`](@ref) and [`MaybeHotVector`](@ref).
For any other `AbstractMatrix` it falls back to standard multiplication.

# Examples
```jldoctest
julia> A = PostImputingMatrix(ones(2, 2), -ones(2))
2×2 PostImputingMatrix{Float64, Matrix{Float64}, Vector{Float64}}:
W:
 1.0  1.0
 1.0  1.0

ψ:
 -1.0
 -1.0

julia> A * maybehotbatch([1, missing], 1:2)
2×2 Matrix{Float64}:
 1.0  -1.0
 1.0  -1.0
```

See also: [`PreImputingMatrix`](@ref).
"""
struct PostImputingMatrix{T <: Number, U <: AbstractMatrix{T}, V <: AbstractVector{T}} <: AbstractMatrix{T}
    W::U
    ψ::V
end

Flux.@functor PostImputingMatrix

"""
    PostImputingMatrix(W::AbstractMatrix{T}, ψ=zeros(T, size(W, 1))) where T

Construct a [`PostImputingMatrix`](@ref) with multiplication parameters `W` and default parameters `ψ`.

# Examples
```jldoctest
julia> PostImputingMatrix([1 2; 3 4])
2×2 PostImputingMatrix{Int64, Matrix{Int64}, Vector{Int64}}:
W:
 1  2
 3  4

ψ:
 0
 0
```

See also: [`PreImputingMatrix`](@ref).
"""
PostImputingMatrix(W::AbstractMatrix{T}) where T = PostImputingMatrix(W, zeros(T, size(W, 1)))

Flux.@forward PostImputingMatrix.W Base.size, Base.getindex, Base.setindex!, Base.firstindex, Base.lastindex

Base.vcat(As::PostImputingMatrix...) = PostImputingMatrix(vcat((A.W for A in As)...), vcat((A.ψ for A in As)...))
function Base.hcat(As::PostImputingMatrix...)
    ArgumentError("It doesn't make sense to hcat PostImputingMatrices") |> throw
end

A::PostImputingMatrix * b::AbstractVector = (_check_mul(A, b); _mul(A, b))
Zygote.@adjoint A::PostImputingMatrix * b::AbstractVector = (_check_mul(A, b); Zygote.pullback(_mul, A, b))
A::PostImputingMatrix * b::MaybeHotVector = (_check_mul(A, b); _mul(A, b))
Zygote.@adjoint A::PostImputingMatrix * b::MaybeHotVector = (_check_mul(A, b); Zygote.pullback(_mul, A, b))
A::PostImputingMatrix * B::AbstractMatrix = (_check_mul(A, B); _mul(A, B))
Zygote.@adjoint A::PostImputingMatrix * B::AbstractMatrix = (_check_mul(A, B); Zygote.pullback(_mul, A, B))
A::PostImputingMatrix * B::MaybeHotMatrix = (_check_mul(A, B); _mul(A, B))
Zygote.@adjoint A::PostImputingMatrix * B::MaybeHotMatrix = (_check_mul(A, B); Zygote.pullback(_mul, A, B))
A::PostImputingMatrix * B::NGramMatrix = (_check_mul(A, B); _mul(A, B))
Zygote.@adjoint A::PostImputingMatrix * B::NGramMatrix = (_check_mul(A, B); Zygote.pullback(_mul, A, B))

_mul(A::PostImputingMatrix, B::AbstractVecOrMat) = A.W * B
_mul(A::PostImputingMatrix, b::MaybeHotVector{Missing}) = A.ψ
_mul(A::PostImputingMatrix, b::MaybeHotVector{<:Integer}) = A.W * b
_mul(A::PostImputingMatrix, B::MaybeHotMatrix{Missing}) = repeat(A.ψ, 1, size(B, 2))
_mul(A::PostImputingMatrix, B::MaybeHotMatrix{<:Integer}) = A.W * B
_mul(A::PostImputingMatrix, B::MaybeHotMatrix{Maybe{T}}) where {T <: Integer} =
    _mul_maybe_hot(A.W, A.ψ, B.I)
_mul(A::PostImputingMatrix, B::NGramMatrix{<:Sequence}) = A.W * B
_mul(A::PostImputingMatrix, B::NGramMatrix{Missing}) = repeat(A.ψ, 1, size(B, 2))
_mul(A::PostImputingMatrix, B::NGramMatrix{Maybe{T}}) where {T <: Sequence} =
    _mul_ngram(A.W, A.ψ, B.S, B.n, B.b, B.m)

_mul_maybe_hot(W, ψ, I) = _impute_maybe_hot(W, ψ, I)[1]
function ChainRulesCore.rrule(::typeof(_mul_maybe_hot), W, ψ, I)
    C, m = _impute_maybe_hot(W, ψ, I)
    function ∇W(Δ)
        dW = zero(W)
        @inbounds for (k, j) in enumerate(I)
            if !ismissing(j)
                for i in 1:size(dW, 1)
                    dW[i, j] += Δ[i, k]
                end
            end
        end
        dW
    end
    C, Δ -> (NoTangent(), @thunk(∇W(Δ)), @thunk(vec(sum(view(Δ, :, m), dims=2))), NoTangent())
end
function _impute_maybe_hot(W, ψ, I)
    m = trues(length(I))
    C = similar(ψ, size(W, 1), length(I))
    C .= ψ
    @inbounds for (j, k) in enumerate(I)
        if !ismissing(k)
            m[j] = false
            for i in 1:size(C, 1)
                C[i, j] = W[i, k]
            end
        end
    end
    C, m
end

# TODO rewrite this to less parameters once Zygote allows for composite grads
_mul_ngram(W, ψ, S, n, b, m) = _impute_ngram(W, ψ, S, n, b, m)
function ChainRulesCore.rrule(::typeof(_mul_ngram), W, ψ, S, n, b, m)
    function ∇W(Δ)
        dW = zero(W)
        z = _init_z(n, b)
        for (k, s) in enumerate(S)
            _∇A_mul_vec!(Δ, k, dW, z, s, n, b, m)
        end
        dW
    end
    _impute_ngram(W, ψ, S, n, b, m), Δ -> (NoTangent(), @thunk(∇W(Δ)),
                                           @thunk(vec(sum(view(Δ, :, ismissing.(S)), dims=2))),
                                  fill(NoTangent(), 4)...)
end

function _impute_ngram(W, ψ, S, n, b, m)
    C = zeros(eltype(ψ), size(W, 1), length(S))
    z = _init_z(n, b)
    for (k, s) in enumerate(S)
        _mul_vec!(C, k, W, z, s, n, b, m, ψ)
    end
    C
end


text:
######### StanSample Gaussian example  ###########

using StanSample, Distributions

gaussian_model = "
data {
  int<lower=0> N;
  vector[N] y;
}
parameters {
  real mu;
  real<lower=0> sigma;
}
model {
  mu ~ normal(0,1);
  sigma ~ cauchy(0,5);
  y ~ normal(mu,sigma);
}
"

simulateGaussian(; μ=0, σ=1, Nd, kwargs...) = (y=rand(Normal(μ, σ), Nd), N=Nd)

count = 0;
iter = 0

while iter < 3
  (y, N) = simulateGaussian(; Nd=5000)
  gaussian_data = Dict("N" => N, "y" => y)

  # Keep tmpdir across multiple runs to prevent re-compilation
  stanmodel = SampleModel("gaussian", gaussian_model)

  rc = stan_sample(stanmodel; data=gaussian_data)

  if success(rc)
    global iter += 1
    samples = read_samples(stanmodel, :dimarray)

    # Show the cmdstan summary output in DataFrame format
    sdf = StanSample.read_summary(stanmodel)
    if sdf[sdf.parameters .== :mu, :ess][1] >= 3999.0 ||
        sdf[sdf.parameters .== :sigma, :ess][1] >= 3999.0
      global count += 1
      println("$(iter), $(count) : mu_ess=$(sdf[sdf.parameters .== :mu, :ess][1]),
          sigma_ess=$(sdf[sdf.parameters .== :sigma, :ess][1])")
    end
  end
end


text:
__precompile__()

module Sundials

using Reexport, DataStructures, Roots
@reexport using DiffEqBase
using DiffEqBase: check_error!
import DiffEqBase: solve

const warnkeywords =
    (:save_idxs, :d_discontinuities, :isoutofdomain, :unstable_check,
     :calck, :progress, :timeseries_steps,
     :internalnorm, :gamma, :beta1, :beta2, :qmax, :qmin, :qoldinit)

function __init__()
    const global warnlist = Set(warnkeywords)
    const global warnida  = union(warnlist, Set((:dtmin,)))
end

const depsfile = joinpath(dirname(dirname(@__FILE__)),"deps","deps.jl")
if isfile(depsfile)
    include(depsfile)
else
    error("Sundials is not properly installed. Please run Pkg.build(\"Sundials\")")
end

export solve, SundialsODEAlgorithm, SundialsDAEAlgorithm, ARKODE, CVODE_BDF, CVODE_Adams, IDA

# some definitions from the system C headers wrapped into the types_and_consts.jl
const DBL_MAX = prevfloat(Inf)
const DBL_MIN = nextfloat(-Inf)
const DBL_EPSILON = eps(Cdouble)

const FILE = Void
const __builtin_va_list = Ptr{Void}

include("wrapped_api/types_and_consts.jl")
include("types_and_consts_additions.jl")

include("handle.jl")
include("nvector_wrapper.jl")

include("wrapped_api/nvector.jl")
include("wrapped_api/libsundials.jl")
include("wrapped_api/sunmatrix.jl")
include("wrapped_api/sunlinsol.jl")
if isdefined(:libsundials_cvodes)
    include("wrapped_api/cvodes.jl")
else
    include("wrapped_api/cvode.jl")
end
include("wrapped_api/arkode.jl")
if isdefined(:libsundials_idas)
    include("wrapped_api/idas.jl")
else
    include("wrapped_api/ida.jl")
end
include("wrapped_api/kinsol.jl")

include("simple.jl")
include("common_interface/verbosity.jl")
include("common_interface/integrator_types.jl")
include("common_interface/integrator_utils.jl")
include("common_interface/callbacks.jl")
include("common_interface/algorithms.jl")
include("common_interface/solve.jl")

##################################################################
# Deprecations
##################################################################

end # module


text:
"""
    Response(Y)

Type for storing response matrix

"""
mutable struct Response

    Y::AbstractArray{Float64,2} # Response matrix Y
end

"""
    Predictors(X, Z, hasXIntercept, hasZIntercept)

Type for storing predictor (covariate) matrices. Also stores boolean
variables hasXIntercept and hasZIntercept (if they are not supplied, they
default to false).

"""
mutable struct Predictors

    # Assumes all categorical variables are coded in appropriate contrasts
    # Matrix of row covariates
    X::AbstractArray{Float64,2}
    # Matrix of column covariates
    Z::AbstractArray{Float64,2}

    # Boolean flag indicating whether X has an intercept
    hasXIntercept::Bool
    # Boolean flag indicating whether Z has an intercept
    hasZIntercept::Bool

    # Usual constructor
    Predictors(X, Z, hasXIntercept, hasZIntercept) =
        new(X, Z, hasXIntercept, hasZIntercept)

    # Modified constructor that sets hasXIntercept and hasZIntercept to false
    # by default
    Predictors(X, Z) = new(X, Z, false, false)
end

"""
    RawData(response, predictors)

Type for storing response and predictor matrices

Also stores dimensions of matrices as n, m, p, and q.
- n = number of rows of X = number of rows of Y
- m = number of rows of Z = number of columns of Y
- p = number of columns of X
- q = number of columns of Z

The constructor will compute n, m, p, and q based on the response and
predictor matrices and assert that they are consistent.

"""
mutable struct RawData

    # Response
    response::Response
    # Predictors
    predictors::Predictors

    # Data dimensions.
    # These will be generated and should not be supplied as input to the
    # constructor.
    n::Int64 # Rows of X = rows of Y
    m::Int64 # Rows of Z = columns of Y
    p::Int64 # Columns of X
    q::Int64 # Columns of Z

    RawData(response, predictors) =
        # Assert that the dimensions in response and predictors match
        (size(predictors.X,1) != size(response.Y,1)) ||
            (size(predictors.Z,1) != size(response.Y,2)) ?
        error("Dimension mismatch.") :
        # Construct RawData object
        new(response, predictors, size(predictors.X,1), size(predictors.Z,1),
            size(predictors.X,2), size(predictors.Z,2))
end

"""
    get_X(data)

Extract X matrix from RawData object

# Arguments

- data = RawData object

# Value

2d array

"""
function get_X(data::RawData)

    return data.predictors.X
end

"""
    get_Z(data)

Extract Z matrix from RawData object

# Arguments

- data = RawData object

# Value

2d array

"""
function get_Z(data::RawData)

    return data.predictors.Z
end

"""
    get_Y(data)

Extract Y matrix from RawData object

# Arguments

- data = RawData object

# Value

2d array

"""
function get_Y(data::RawData)

    return data.response.Y
end


text:
# ---
# title: 988. Smallest String Starting From Leaf
# id: problem988
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Tree, Depth-first Search
# link: <https://leetcode.com/problems/smallest-string-starting-from-leaf/description/>
# hidden: true
# ---
#
# Given the `root` of a binary tree, each node has a value from `0` to `25`
# representing the letters `'a'` to `'z'`: a value of `0` represents `'a'`, a
# value of `1` represents `'b'`, and so on.
#
# Find the lexicographically smallest string that starts at a leaf of this tree
# and ends at the root.
#
# _(As a reminder, any shorter prefix of a string is lexicographically smaller:
# for example,`"ab"` is lexicographically smaller than `"aba"`.  A leaf of a
# node is a node that has no children.)_
#
#
#
# **Example 1:**
#
# **![](https://assets.leetcode.com/uploads/2019/01/30/tree1.png)**
#
#
#
#     Input: [0,1,2,3,4,3,4]
#     Output: "dba"
#
#
# **Example 2:**
#
# **![](https://assets.leetcode.com/uploads/2019/01/30/tree2.png)**
#
#
#
#     Input: [25,1,3,1,3,0,2]
#     Output: "adz"
#
#
# **Example 3:**
#
# **![](https://assets.leetcode.com/uploads/2019/02/01/tree3.png)**
#
#
#
#     Input: [2,2,1,null,1,0,null,0]
#     Output: "abc"
#
#
#
#
# **Note:**
#
#   1. The number of nodes in the given tree will be between `1` and `8500`.
#   2. Each node in the tree will have a value between `0` and `25`.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using Documenter, Gadfly

load_dir(x) = map(file -> joinpath("lib", x, file), readdir(joinpath(Base.source_dir(), "src", "lib", x)))

makedocs(
    modules = [Gadfly],
    clean = false,
    format = :html,
    sitename = "Gadfly.jl",
    pages = Any[
        "Home" => "index.md",
        "Tutorial" => "tutorial.md",
        "Manual" => Any[
            "Plotting" => "man/plotting.md",
            "Layers and Stacks" => "man/layers.md",
            "Backends" => "man/backends.md",
            "Themes" => "man/themes.md",
        ],
        "Library" => Any[
            hide("Geometries" => "lib/geometries.md", load_dir("geoms")),
            hide("Guides" => "lib/guides.md", load_dir("guides")),
            hide("Statistics" => "lib/stats.md", load_dir("stats")),
            hide("Coords" => "lib/coords.md", load_dir("coords")),
            hide("Scales" => "lib/scales.md", load_dir("scales")),
            "Shapes" => "lib/shapes.md",
        ],
        "Development" => Any[
            "Rendering Pipeline" => "dev/pipeline.md",
            "Regression Testing" => "dev/regression.md",
        ]
    ]
)

deploydocs(
    repo   = "github.com/GiovineItalia/Gadfly.jl.git",
    julia  = "0.6",
    osname = "linux",
    deps = nothing,
    make = nothing,
    target = "build",
)


text:
sucesses = @parallel (+) for i in 1:10^8
    Int(any(bitrand(3)))
end

max_stars_repo_name: prcastro/IntroJulia

max_issues_repo_name: prcastro/IntroJulia

max_forks_repo_name: prcastro/IntroJulia


text:
"""
Comparison of means: plot posterior densities and compute odds against difference of 0.

Returns:
        diff_mean, draws_m1, draws_m2, qs, tst       ## out-of-date: plt,diff_mean,qs,tst

        plt = plot of densities
        diff_mean = MC sample from posterior of difference in means
        qs = HPD quantile intervals from diff_mean
        tst = posterior odds against zero difference and area in tail to one side of zero.

Enter either:
        compare_mean(x,y)
        x, y are two samples to compare_means

        compare_means(m1, m2, s1, s2, n1, n2)
        m1, m2, are sample means, s1, s2, are sample SDs, n1, n2 are sample number of obs.

Optional arguments:
        M = MC sample size (larger for improvement numerical accuracy)
        lbl = ["mu 1 label" "mu 1 label"],  use lbl = ["" ""] to eliminate labels
        lgd = :topright (choose legend position, :topleft to move to top left of figure, etc.)
"""
function compare_means(m1, m2, s1, s2, n1, n2; M = 1000000,  lbl=["mu 1" "mu 2"],lgd = :topright)
    draws_m1 = m1 .+ (s1/sqrt(n1)).*rand(TDist(n1-1),M)
    draws_m2 = m2 .+ (s2/sqrt(n2)).*rand(TDist(n2-1),M)
    diff_mean = draws_m1 - draws_m2
    plot(draws_m1,st=:histogram,normalize=true,bins=100)
    l = @layout([a b])
    plt1 = plot(draws_m1,st=:density,fill=(0,0.4,:red),alpha=0.4,label=lbl[1],legend=lgd,title="Posteriors for means")
    plot!(draws_m2,st=:density,fill=(0,0.4,:blue),alpha=0.4,label=lbl[2])
    plt2 = plot(diff_mean,st=:density,fill=(0,0.4,:green),alpha=0.4,label="",title="Posterior difference")
    vline!([0.0],color=:black,label="")
    plt3 = plot(plt1, plt2, layout=l)
    qs = quantile(diff_mean,[0.005 0.025 0.05 0.5 0.95 0.975 0.995])
    @show(qs)
    tst = [mcodds(diff_mean) bayespval(diff_mean)]
    @show(tst)
    return diff_mean, draws_m1, draws_m2, qs, tst, plt3
end

# function overload for using two samples directly
function compare_means(x,y; M = 1000000, lbl=["mu 1" "mu 2"],lgd = :topright)
    m1 = mean(x); m2 = mean(y)
    s1 = std(x); s2 = std(y)
    n1 = length(x); n2 = length(y)
    draws_m1 = m1 .+ (s1/sqrt(n1)).*rand(TDist(n1-1),M)
    draws_m2 = m2 .+ (s2/sqrt(n2)).*rand(TDist(n2-1),M)
    diff_mean = draws_m1 - draws_m2
    plot(draws_m1,st=:histogram,normalize=true,bins=100)
    l = @layout([a b])
    plt1 = plot(draws_m1,st=:density,fill=(0,0.4,:red),alpha=0.4,label=lbl[1],legend=lgd,title="Posteriors for means")
    plot!(draws_m2,st=:density,fill=(0,0.4,:blue),alpha=0.4,label=lbl[2])
    plt2 = plot(diff_mean,st=:density,fill=(0,0.4,:green),alpha=0.4,label="",title="Posterior difference")
    vline!([0.0],color=:black,label="")
    plt3 = plot(plt1, plt2, layout=l)
    qs = quantile(diff_mean,[0.005 0.025 0.05 0.5 0.95 0.975 0.995])
    @show(qs)
    tst = [mcodds(diff_mean) bayespval(diff_mean)]
    @show(tst)
    return diff_mean, draws_m1, draws_m2, qs, tst, plt3
end

"""
 2x2 contigency table (categorical data) comparison of proportions
 """
function compare_proportions(s1, s2, n1, n2; a =1, b = 1,M = 1000000, lbl=["mu 1" "mu 2"],lgd = :topright)
    draws_m1 = rand(Beta((s1+a),(n1-s1+b)),M)
    draws_m2 = rand(Beta((s2+a),(n2-s2+b)),M)
    diff_mean = draws_m1 - draws_m2
    plot(draws_m1,st=:histogram,normalize=true,bins=100)
    l = @layout([a b])
    plt1 = plot(draws_m1,st=:density,fill=(0,0.4,:red),alpha=0.4,label=lbl[1],legend=lgd,title="Posteriors for proportions",ylab="Density",xlab="Proportion")
    plot!(draws_m2,st=:density,fill=(0,0.4,:blue),alpha=0.4,label=lbl[2])
    plt2 = plot(diff_mean,st=:density,fill=(0,0.4,:green),alpha=0.4,label="",title="Posterior difference",xlab="Difference in proportion")
    vline!([0.0],color=:black,label="",linewidth=2)
    plt3 = plot(plt1, plt2, layout=l)
    qs = quantile(diff_mean,[0.005 0.025 0.05 0.5 0.95 0.975 0.995])
    @show(qs)
    tst = [mcodds(diff_mean) bayespval(diff_mean)]
    @show(tst)
    return diff_mean, draws_m1, draws_m2, qs, tst, plt3
end

function compare_proportions(x, y; a =1, b = 1,M = 1000000, lbl=["mu 1" "mu 2"],lgd = :topright)
    n1 = length(x); s1 = sum(x)
    n2 = length(y); s2 = sum(y)
    draws_m1 = rand(Beta((s1+a),(n1-s1+b)),M)
    draws_m2 = rand(Beta((s2+a),(n2-s2+b)),M)
    diff_mean = draws_m1 - draws_m2
    plot(draws_m1,st=:histogram,normalize=true,bins=100)
    l = @layout([a b])
    plt1 = plot(draws_m1,st=:density,fill=(0,0.4,:blue),alpha=0.4,label=lbl[1],legend=lgd,title="Posteriors for proportions")
    plot!(draws_m2,st=:density,fill=(0,0.4,:red),alpha=0.4,label=lbl[2])
    plt2 = plot(diff_mean,st=:density,fill=(0,0.4,:green),alpha=0.4,label="",title="Posterior difference")
    vline!([0.0],color=:black,label="")
    plt3 = plot(plt1, plt2, layout=l)
    qs = quantile(diff_mean,[0.005 0.025 0.05 0.5 0.95 0.975 0.995])
    @show(qs)
    tst = [mcodds(diff_mean) bayespval(diff_mean)]
    @show(tst)
    return diff_mean, draws_m1, draws_m2, qs, tst, plt3
end

"""
Function to draw from Beta(s+a, n-s+b) posterior

     Default prior is uniform: a = b = 1
"""
function beta_posterior(s,n; a=1, b=1, M=100000)
    theta_draws = rand(Beta((s+a),(n-s+b)),M)
    return theta_draws
end

"""
Function to update Beta(s+a, n-s+b) posterior given two samples

     Default prior is uniform: a = b = 1 for each sample distribution
"""
function beta_update(s1,n1,s2,n2; a1=1,b1=1,a2=1,b2=1,M=100000)
    theta_draws = rand(Beta((s1+s2+a1+a2),(n1+n2-s1-s2+b1+b2)),M)
    return theta_draws
end

"""
equiv_test = Equivalence test from MC sample of difference in ATE

    mc_sample = posterior simulation sample of difference in treatments

    interval = 1/2 tolerance interval  = distance from 0.0 to one side of interval

    see TADS_equiv.jl for example of use.
"""
function equiv_test(mc_sample,interval;h0=0.0)
  p1 = p2 = p3 = p4 = p5 = 0.0
  for th1 in mc_sample
      if th1 < (h0-3*interval)
          p1 += 1
      elseif th1 < (h0-interval)
          p2 += 1
      elseif th1 < (h0+interval)
          p3 += 1
      elseif th1 < (h0+3*interval)
          p4 += 1
      else
          p5 += 1
      end
  end
  ps = [p1; p2; p3; p4; p5]
  normps = ps./sum(p1 + p2 + p3 + p4 + p5)
  x = [-4*interval; -2*interval; 0.0; 2*interval; 4*interval]
#  plt = plot(mc_sample,st=:density,label="",ylab="",alpha=0.4,fill=(0,0.4,:green),ylim=(0.0,0.4))
#  plot!(twinx(),x,normps, st=:bar, color=:grey,alpha=0.4,bar_width=2*interval,label="",ylim=(0.0,1.0),ylab="Probability",title = "Equivalence interval probabilities") #,axis=:right
  plt = plot(mc_sample,st=:density,label="",ylab="",alpha=0.2,fill=(0,0.2,:blue)) #,ylim=(0.0,0.4))
  plot!(twinx(),x,normps, st=:bar, color=:green,alpha=0.4,bar_width=2*interval,label="",ylim=(0.0,1.0),ylab="",title = "Equivalence interval probabilities") #,axis=:right
  oddsf35 = normps[3]/normps[5]
  oddsf34 = normps[3]/normps[4]
  oddsf32 = normps[3]/normps[2]
  oddsf31 = normps[3]/normps[1]
  eqodds = [oddsf31 oddsf32 oddsf34  oddsf35]
#  @show(eqodds)
#  @show(normps)
  return eqodds, normps, plt, x
end

"""
Equivalance test with 7 intervals (equiv_test has 5)
"""
function equiv_test7(mc_sample,interval;h0=0.0)
    p1 = p2 = p3 = p4 = p5 = p6 = p7 = 0.0
    for th1 in mc_sample
        if th1 < (h0-5*interval)
            p1 += 1
        elseif th1 < (h0-3*interval)
            p2 += 1
        elseif th1 < (h0-interval)
            p3 += 1
        elseif th1 < (h0+interval)
            p4 += 1
        elseif th1 < (h0+3*interval)
            p5 += 1
        elseif th1 < (h0+5*interval)
            p6 += 1
        else
            p7 += 1
        end
    end
    ps = [p1; p2; p3; p4; p5; p6; p7]
    normps = ps./sum(p1 + p2 + p3 + p4 + p5 + p6 + p7)
    x = [-6*interval; -4*interval; -2*interval; 0.0; 2*interval; 4*interval; 6*interval;]
    #  plt = plot(mc_sample,st=:density,label="",ylab="",alpha=0.4,fill=(0,0.4,:green),ylim=(0.0,0.4))
    #  plot!(twinx(),x,normps, st=:bar, color=:grey,alpha=0.4,bar_width=2*interval,label="",ylim=(0.0,1.0),ylab="Probability",title = "Equivalence interval probabilities") #,axis=:right
    plt = plot(mc_sample,st=:density,label="",ylab="",alpha=0.2,fill=(0,0.2,:blue)) #,ylim=(0.0,0.4))
    plot!(twinx(),x,normps, st=:bar, color=:green,alpha=0.4,bar_width=2*interval,label="",ylim=(0.0,1.0),ylab="",title = "Equivalence interval probabilities") #,axis=:right
    oddsf47 = normps[4]/normps[7]
    oddsf46 = normps[4]/normps[6]
    oddsf45 = normps[4]/normps[5]
    oddsf43 = normps[4]/normps[3]
    oddsf42 = normps[4]/normps[2]
    oddsf41 = normps[4]/normps[1]
    eqodds = [oddsf41 oddsf42 oddsf43 oddsf45 oddsf46 oddsf47]
    #  @show(eqodds)
    #  @show(normps)
    return eqodds, normps, plt, x
end


text:
function countSymbols(intNumber)
	return length(matchall(r"1", bin(intNumber,2)))
end

max_stars_repo_name: driesverhees/dobble

max_issues_repo_name: driesverhees/dobble

max_forks_repo_name: driesverhees/dobble


text:
import Base: eltype

# AtomGaussian

struct AtomGaussian{N,T<:Real} <: AbstractIsotropicGaussian{N,T}
    μ::SVector{N,T}
    σ::T
    ϕ::T
    dirs::Vector{SVector{N,T}}
    node::SDFileAtom
    nodeidx::Int
end

"""
"""
function AtomGaussian(atom::SDFileAtom, nodeidx, σ, ϕ)
    N = length(atom.coords)
    T = eltype(atom.coords)
    return AtomGaussian{N,T}(SVector{N,T}(atom.coords), T(σ), T(ϕ), SVector{N,T}[], atom, nodeidx)
end

AtomGaussian(mol::Union{UndirectedGraph, SubgraphView}, nodeidx, σdict, ϕdict
    ) = AtomGaussian(nodeattr(mol,nodeidx), nodeidx, σdict[nodeidx], ϕdict[nodeidx])

function AtomGaussian(μ::AbstractArray, σ::Real, ϕ::Real, dirs::AbstractArray=SVector{length(μ),eltype(μ)}[], node=SDFileAtom(), nodeidx=0)
    t = promote_type(eltype(μ), typeof(σ), typeof(ϕ), eltype(eltype(dirs)))
    return AtomGaussian{length(μ),t}(SVector{length(μ),t}(μ), t(σ), t(ϕ), SVector{length(μ),t}[SVector{length(μ),t}(dir/norm(dir)) for dir in dirs], node, nodeidx)
end

convert(::Type{AtomGaussian{N,T}}, g::AtomGaussian) where {N,T} = AtomGaussian(SVector{N,T}(g.μ), T(g.σ), T(g.ϕ), Vector{SVector{N,T}}(g.dirs), g.node, g.nodeidx)

# MolGMM

struct MolGMM{N,T<:Real} <: AbstractIsotropicGMM{N,T}
    gaussians::Vector{AtomGaussian{N,T}}
    graph::Union{UndirectedGraph,SubgraphView}
    nodes::Set{Int}
    σfun::Function
    ϕfun::Function
end

eltype(::Type{MolGMM{N,T}}) where {N,T} = AtomGaussian{N,T}

"""
    model = MolGMM(mol, σfun=ones, ϕfun=ones, nodes=nodeset(mol))

Creates a Gaussian mixture model from a molecule or subgraph `mol`, which consists of spherical Gaussian distributions
with means `μ` equal to atom coordinates.

Optionally, functions `σfun` and `ϕfun` can be provided, which take `mol` as input and return dictionaries mapping
mapping node indicies to variances `σ` and scaling coefficients `ϕ`, respectively.

If `nodes` is provided, the Gaussian mixture model will be constructed only from atoms corresponding to the node
indexes of the molecule's graph.
"""
function MolGMM(mol::UndirectedGraph,
                nodes=nodeset(mol);
                σfun = vdwvolume_sigma,
                ϕfun = ones)
    N = length(nodeattrs(mol)[1].coords)
    T = eltype(nodeattrs(mol)[1].coords)
    σdict, ϕdict = Dict{Int, T}(σfun(mol)), Dict{Int, T}(ϕfun(mol))
    atoms = [AtomGaussian(nodeattr(mol, i), i, σdict[i], ϕdict[i]) for i in nodes]
    # remove atoms with ϕ = 0
    filter!(atom->atom.ϕ≠0,atoms)
    return MolGMM{N,T}(atoms, mol, nodes, σfun, ϕfun)
end

MolGMM(submol::SubgraphView; kwargs...) = MolGMM(submol.graph, nodeset(submol); kwargs...)

struct PharmacophoreGMM{N,T<:Real,K} <: AbstractIsotropicMultiGMM{N,T,K}
    gmms::Dict{K, IsotropicGMM{N,T}}
    graph::Union{UndirectedGraph,SubgraphView}
    nodes::Set{Int}
    σfun::Function
    ϕfun::Function
    features::Vector{K}
    directional::Bool
end

eltype(::Type{PharmacophoreGMM{N,T,K}}) where {N,T,K} = Pair{K, IsotropicGMM{N,T}}

"""
    model = PharmacophoreGMM(mol, σfun=vdwvolume_sigma, ϕfun=ones, nodes=nodeset(mol), features=pubchem_features)

Creates a set Gaussian mixture models from a molecule or subgraph `mol`, with each model corresponding to a
particular type of molecular feature (e.g. ring structures)

Optionally, functions `σfun` and `ϕfun` can be provided, which take `mol` as input and return dictionaries
mapping node indicies to variances `σ` and scaling coefficients `ϕ`, respectively.

If `nodes` is provided, the Gaussian mixture models will be constructed only from atoms corresponding to the node
indexes of the molecule's graph.

`features` specifies the keys of allowed molecular features. Features with other keys will be ignored.

`directional` specifies whether or not geometric constraints for ring structures and hydrogen bond donors will be included.
"""
function PharmacophoreGMM(mol::UndirectedGraph,
                          nodes = nodeset(mol);
                          σfun = vdwvolume_sigma,
                          ϕfun = ones,
                          features = pubchem_features,
                          directional = true)
    dim = length(nodeattrs(mol)[1].coords)
    valtype = eltype(nodeattrs(mol)[1].coords)
    # add a GMM for each type of feature
    gmms = Dict{Symbol,IsotropicGMM{dim,valtype}}()
    for (key, nodesets) in pharmfeatures!(mol)
        if key in features
            feats = IsotropicGaussian{dim,valtype}[]
            for set in nodesets
                # check if all atoms in the feature are represented by allowed nodes
                if set ∩ nodes == set
                    # if it is a ring feature, add normal vectors
                    if (key == :rings || key == :aromaticrings) && directional
                        center = sum([nodeattr(mol,idx).coords for idx in set])/length(set)
                        coordmat = fill(zero(valtype), 3, length(set))
                        for (i,idx) in enumerate(set)
                            coordmat[:,i] = nodeattr(mol,idx).coords .- center
                        end
                        n = GaussianMixtureAlignment.planefit(coordmat)[1]
                        dirs = [n, -n]
                    # if it is a hydrogen bond donor, point outward from hydrogens
                    elseif key == :donor && directional
                        if length(set)==1
                            dirs = Vector{valtype}[]
                            idx = collect(set)[1]
                            for (edge,neigh) in neighbors(mol,idx)
                                if nodeattr(mol,neigh).symbol == :H
                                    push!(dirs, nodeattr(mol,neigh).coords - nodeattr(mol,idx).coords)
                                end
                            end
                        else
                            dirs = nothing
                        end
                    # if it is a hydrogen bond acceptor, point along the lone pair orbitals
                    elseif key == :acceptor && directional
                        # the acceptor should be a single atom
                        if length(set) != 1
                            dirs = nothing
                        else
                            idx = collect(set)[1]
                            neighs = collect(neighbors(mol,idx))
                            if hybridization(mol)[idx] == :sp3      # tetrahedral
                                if length(neighs) == 3
                                    # if there are 3 bonds
                                    bonddir1 = nodeattr(mol,neighs[1].second).coords .- nodeattr(mol,idx).coords
                                    bonddir2 = nodeattr(mol,neighs[2].second).coords .- nodeattr(mol,idx).coords
                                    bonddir3 = nodeattr(mol,neighs[3].second).coords .- nodeattr(mol,idx).coords
                                    dirs = Vector{valtype}[-normalize(bonddir1.+bonddir2.+bonddir3)]
                                elseif length(neighs) == 2
                                    # if there are 2 bonds
                                    bonddir1 = nodeattr(mol,neighs[1].second).coords .- nodeattr(mol,idx).coords
                                    bonddir2 = nodeattr(mol,neighs[2].second).coords .- nodeattr(mol,idx).coords
                                    axis1 = normalize(bonddir1.+bonddir2)
                                    axis2 = normalize(cross(bonddir1,bonddir2))
                                    R = AngleAxis(π,axis2...) * AngleAxis(π/2,axis1...)
                                    dirs = Vector{valtype}[R*bonddir1, R*bonddir2]
                                else
                                    dirs = nothing
                                end
                            elseif hybridization(mol)[idx] == :sp2  # trigonal planar
                                # if there are two bonds, use them to find the lone pair's direction
                                if length(neighs) != 2
                                    dirs = nothing
                                else
                                    dirs = Vector{valtype}[]
                                    bonddir1 = nodeattr(mol,neighs[1].second).coords .- nodeattr(mol,idx).coords
                                    bonddir2 = nodeattr(mol,neighs[2].second).coords .- nodeattr(mol,idx).coords
                                    axis = normalize(cross(bonddir1, bonddir2))
                                    push!(dirs, AngleAxis(4π/3, axis...)*bonddir1)
                                end
                            elseif hybridization(mol)[idx] == :sp   # linear
                                # direction points along the axis of the atom's single bond
                                if length(neighs) != 1
                                    dirs = nothing
                                else
                                    dirs = Vector{valtype}[nodeattr(nodeattr(mol,idx) .- mol,neighs[1].second.coords)]
                                end
                            end
                        end
                    # otherwise there is no geometric constraint
                    else
                        dirs = nothing
                    end
                    push!(feats, atoms_to_feature(mol, set, dirs; σfun=σfun, ϕfun=ϕfun))
                end
            end
            push!(gmms, Pair(key, IsotropicGMM(feats)))
        end
    end
    gmms = gmms
    return PharmacophoreGMM(gmms, mol, nodes, σfun, ϕfun, features, directional)
end

PharmacophoreGMM(submol::SubgraphView; kwargs...) = PharmacophoreGMM(submol.graph, nodeset(submol); kwargs...)

# descriptive display

Base.show(io::IO, molgmm::MolGMM) = println(io,
    "MolGMM from molecule with formula $(typeof(molgmm.graph)<:GraphMol ? molecularformula(molgmm.graph) : molecularformula(molgmm.graph.graph))",
    " with $(length(molgmm)) Gaussians."
)

Base.show(io::IO, pgmm::PharmacophoreGMM) = println(io,
    summary(pgmm),
    " from molecule with formula $(typeof(pgmm.graph)<:GraphMol ? molecularformula(pgmm.graph) : molecularformula(pgmm.graph.graph))",
    " with $(sum(length(gmm.second) for gmm in pgmm)) Gaussians in $(length(pgmm)) GMMs with labels:\n",
    "$([label for (label, gmm) in pgmm.gmms])"
)


text:
const SpatialDim = Union{XDim,YDim,ZDim}

"""
    Band <: Dimension

    Band(val=:)

Band [`Dimension`]($DDdimdocs) for multi-band rasters.

## Example:
```julia
banddim = Band(10:10:100)
# Or
val = A[Band(1)]
# Or
mean(A; dims=Band)
```
"""
@dim Band

"""
    mappedbounds(x)

Get the bounds converted to the [`mappedcrs`](@ref) value.

Whithout ArchGDAL loaded, this is just the regular bounds.
"""
function mappedbounds end

mappedbounds(dims::Tuple) = map(mappedbounds, dims)
mappedbounds(dim::Dimension) = mappedbounds(mode(dim), dim)
mappedbounds(::IndexMode, dim) = bounds(dim)
mappedbounds(mode::Projected, dim) = mappedbounds(mappedcrs(mode), mode, dim)
mappedbounds(mappedcrs::Nothing, mode::Projected, dim) =
    error("No mappedcrs attached to $(name(dim)) dimension")
mappedbounds(mappedcrs::GeoFormat, mode::Projected, dim) =
    _sort(reproject(crs(mode), mappedcrs, dim, bounds(dim)))

projectedbounds(dims::Tuple) = map(projectedbounds, dims)
projectedbounds(dim::Dimension) = projectedbounds(mode(dim), dim)
projectedbounds(::IndexMode, dim) = bounds(dim)
projectedbounds(mode::Mapped, dim) = projectedbounds(crs(mode), mode, dim)
projectedbounds(crs::Nothing, mode::Mapped, dim) =
    error("No projection crs attached to $(name(dim)) dimension")
projectedbounds(crs::GeoFormat, mode::Mapped, dim) =
    _sort(reproject(mappedcrs(mode), crs, dim, bounds(dim)))

_sort((a, b)) = a <= b ? (a, b) : (b, a)

"""
    mappedindex(x)

Get the index value of a dimension converted to the `mappedcrs` value.

Whithout ArchGDAL loaded, this is just the regular dim value.
"""
function mappedindex end

mappedindex(dims::Tuple) = map(mappedindex, dims)
mappedindex(dim::Dimension) = _mappedindex(mode(dim), dim)

_mappedindex(::IndexMode, dim::Dimension) = index(dim)
_mappedindex(mode::Projected, dim::Dimension) = _mappedindex(mappedcrs(mode), mode, dim)
_mappedindex(mappedcrs::Nothing, mode::Projected, dim) =
    error("No mappedcrs attached to $(name(dim)) dimension")
_mappedindex(mappedcrs::GeoFormat, mode::Projected, dim) =
    reproject(crs(dim), mappedcrs, dim, index(dim))

projectedindex(dims::Tuple) = map(projectedindex, dims)
projectedindex(dim::Dimension) = _projectedindex(mode(dim), dim)

_projectedindex(::IndexMode, dim::Dimension) = index(dim)
_projectedindex(mode::Mapped, dim::Dimension) = _projectedindex(crs(mode), mode, dim)
_projectedindex(crs::Nothing, mode::Mapped, dim::Dimension) =
    error("No projection crs attached to $(name(dim)) dimension")
_projectedindex(crs::GeoFormat, mode::Mapped, dim::Dimension) =
    reproject(mappedcrs(dim), crs, dim, index(dim))


text:
type edb
  events::DataFrame
  event_times::Vector
  cd::ASCIIString
end

"""
Generate an event database, where all individuals are suscpetible at time 0.
Specify whether for an SI model or an SIR model.
"""
function create_event_db(pop_db, cd="discrete", gamma=Inf)
  if gamma == Inf
    return edb(DataFrame(ind_id = pop_db[:,1], s = fill(0, size(pop_db)[1]), i = fill(NaN, size(pop_db)[1])), fill(Inf, (size(pop_db)[1])), cd)
  end
  if 0 < gamma < Inf
    return edb(DataFrame(ind_id = pop_db[:,1], s = fill(0, size(pop_db)[1]), i = fill(NaN, size(pop_db)[1]), r = fill(NaN, size(pop_db)[1])), fill(Inf, 2*(size(pop_db)[1])), cd)
  end
end

"""
Update the ordered sequence of event times based on a new event at
`time`
"""
function event_time_update!(eventtime, event_db::edb)
  maxevents=length(event_db.event_times)
  if length(size(eventtime)) > 1
    error("Error: Incorrect dimensions of event times (should be 1 dimensional)")
  end
  if length(eventtime) > size(event_db.events)[1]
    error("Error: Impossible number of event times (exceeds population size)")
  end
  if sum(event_db.event_times .== Inf) < length(eventtime)
    warn("Warning: More event times inputted than there are events left to occur")
  end
  if 1 < length(eventtime) < size(event_db.events)[1]
    for i = 1:length(eventtime)
      event_db.event_times=[event_db.event_times[(1:maxevents)[eventtime[i] .>= event_db.event_times]]; eventtime[i]; event_db.event_times[(1:(maxevents-1))[eventtime[i] .< (event_db.event_times[1:(maxevents-1)])]]]
    end
  end
  if 1 == length(eventtime)
    event_db.event_times=[event_db.event_times[(1:maxevents)[eventtime .>= event_db.event_times]]; eventtime; event_db.event_times[(1:(maxevents-1))[eventtime .< (event_db.event_times[1:(maxevents-1)])]]]
  end
end

"""
Randomly infect 1 individual at time = 1.0, A random infection is required
at the beginning of all simulations. Specify whether ILM is continuous
("c") or discrete ("d"), which selects between a exponential or geometric
recovery time functions. If `gamma` (the mean recovery time) is > 0, then
a recovery for this infection will also be probabilistically generated.
"""
function initial_infect!(event_db, gamma=Inf, initial=0)
  if initial==0
    chosen_one=sample(event_db.events[:,1], 1)
  else
    chosen_one=event_db.events[initial, 1]
  end
  event_db.events[chosen_one, 3] = 1.0
  event_time_update!(1.0, event_db)
  if 0 < gamma < Inf && size(event_db.events)[2] == 4
    if event_db.cd == "continuous"
      recovery_dist = Exponential(gamma)
    end
    if event_db.cd == "discrete"
      recovery_dist = Geometric(1/gamma)
    end
    recovery_time = 1.0 + rand(recovery_dist)
    event_db.events[chosen_one, 4] = recovery_time
    event_time_update!(recovery_time, event_db)
  end
end

"""
Find the individuals falling into `state` (S, I, or R), at `time`
"""
function find_state(event_db::edb, time, state)
  state_index = fill(false, size(event_db.events)[1])
  if state == "S"
    for i = 1:length(state_index)
      state_index[i] = isnan(event_db.events[i,3]) || event_db.events[i,3] > time
    end
  end
  if state == "I"
    if size(event_db.events)[2] == 4
      for i = 1:length(state_index)
        state_index[i] = event_db.events[i,3] <= time && ((isnan(event_db.events[i,4])) || (time < event_db.events[i,4]))
      end
    end
    if size(event_db.events)[2] == 3
      for i = 1:length(state_index)
        state_index[i] = event_db.events[i,3] <= time
      end
    end
  end
  if state == "R"
    for i = 1:length(state_index)
      state_index[i] = event_db.events[i,4] <= time
    end
  end
  return state_index
end

"""
Find the alpha*(distance matrix^-beta)
"""
function dist_ab_mtx(distance_mat, alpha, beta)
  return alpha .* (distance_mat .^ -beta)
end

"""
Determine infection probabilities for all susceptible individuals
"""
function infection_probabilities(distance_mat_alphabeta, infectious, susceptible)
  return 1 - exp(-sum(distance_mat_alphabeta[susceptible, infectious], 2))
end

"""
Generate infection times (exponentially distributed) based on current
infectious and susceptible. The minimum time will become the next
infected individual, remaining times will need to be recalculated
"""
function infection_times(distance_mat_alphabeta, infectious, susceptible)
  exponential_rate = fill(Inf, length(susceptible))
  exponential_rate[susceptible] = sum(distance_mat_alphabeta[susceptible, infectious], 2).^-1.0
  infect_times = fill(Inf, length(susceptible))
  for i = 1:length(susceptible)
    if 0 < exponential_rate[i] < Inf
      infect_times[i]= rand(Exponential(exponential_rate[i]))
    end
  end
  return infect_times
end

"""
Propagate infection through a population, default is for a discrete model.
Use of an SIR or SI model inferred from dimensions of `event_db` and
specifcation of a `gamma` > 0. Default time of previous time step is based
on the maximum infection time - which may not work well in discrete SIR
models especially
"""
function infect_recover!(distance_mat_alphabeta,
                         event_db::edb,
                         time=1.0::Float64,
                         gamma=Inf::Float64)
  susceptible = find_state(event_db, time, "S")
  infectious = find_state(event_db, time, "I")
  if event_db.cd == "discrete"
    infect_probs = fill(0.0, length(susceptible))
    infect_probs[susceptible]=infection_probabilities(distance_mat_alphabeta, infectious, susceptible)
    infected = fill(false, length(susceptible))
    for i in 1:length(infect_probs)
      if 0.0 < infect_probs[i] <=1.0
        infected[i] = rand() < infect_probs[i]
      end
    end
    if sum(infected) > 0
      event_db.events[infected, 3] = time + 1.0
      event_time_update!(fill(time + 1.0, sum(infected)), event_db)
      if 0 < gamma < Inf && size(event_db.events)[2] == 4
        recovery_times = rand(Geometric(1/gamma), sum(infected))
        event_db.events[infected, 4] = time + 2.0 + recovery_times
        event_time_update!(time + 2.0 + recovery_times, event_db)
      end
    end
  end
  if event_db.cd == "continuous"
    if gamma == Inf && size(event_db.events)[2] == 3
      infect_times = infection_times(distance_mat_alphabeta, infectious, susceptible)
      which_min = infect_times .== minimum(infect_times)
      event_db.events[which_min, 3] = time + infect_times[which_min]
      event_time_update!(time + infect_times[which_min], event_db)
    end
    if 0 < gamma < Inf && size(event_db.events)[2] == 4
      infect_times = infection_times(distance_mat_alphabeta, infectious, susceptible)
      which_min = infect_times .== minimum(infect_times)
      while (time .< event_db.event_times) != ((time + infect_times[which_min][1]) .<= event_db.event_times)
        time = event_db.event_times[time .< event_db.event_times][1]
        susceptible = find_state(event_db, time, "S")
        infectious = find_state(event_db, time, "I")
        infect_times = infection_times(distance_mat_alphabeta, infectious, susceptible)
        if minimum(infect_times) == Inf
          break
        end
        which_min = infect_times .== minimum(infect_times)
      end
      if minimum(infect_times) != Inf
        event_db.events[which_min,3] = (time + infect_times[which_min])
        recovery_time = rand(Exponential(gamma),sum(which_min))
        event_db.events[which_min,4] = event_db.events[which_min,3] + recovery_time
        event_time_update!(time + infect_times[which_min], event_db)
        event_time_update!(time + infect_times[which_min] + recovery_time, event_db)
      end
    end
  end
end

"""
Function to generate initial infection then loop infect_recover function as appropriate
for continuous and discrete SI and SIR models
"""
function infect_recover_loop(pop_db::DataFrame,
                             cd="discrete",
                             ilm="SI",
                             alpha=1.::Float64,
                             beta=1.::Float64,
                             gamma=Inf::Float64,
                             limit=100::Float64,
                             initial=0)
  distance_mat_alphabeta = dist_ab_mtx(create_dist_mtx(pop_db), alpha, beta)
  event_db = create_event_db(pop_db, cd, gamma)
  initial_infect!(event_db, gamma, initial)
  time=1.0
  if event_db.cd == "discrete"
    if gamma == Inf && size(event_db.events)[2] == 3
      while sum(find_state(event_db, time, "S")) > 0 && time < limit-1
        infect_recover!(distance_mat_alphabeta, event_db, time, Inf)
        time += 1.0
      end
    end
    if 0 < gamma < Inf && size(event_db.events)[2] == 4
      while sum(find_state(event_db, time, "I")) > 0 && time < limit-1
        infect_recover!(distance_mat_alphabeta, event_db, time, gamma)
        time += 1.0
      end
    end
    if time == limit-1
        warn("Warning: Simulation was halted after limit reached")
    end
  end
  if event_db.cd == "continuous"
    if gamma == Inf && size(event_db.events)[2] == 3
      while sum(find_state(event_db, time, "S")) > 0 && time < limit
        infect_recover!(distance_mat_alphabeta, event_db, time, Inf)
        time = maximum(event_db.events[:,3])
      end
    end
    if 0 < gamma < Inf && size(event_db.events)[2] == 4
      while sum(find_state(event_db, time, "I")) > 0 && time < limit
        infect_recover!(distance_mat_alphabeta, event_db, time, gamma)
        if time == maximum(event_db.events[:,3])
          break
        end
        time = maximum(event_db.events[:,3])
      end
    end
    if time >= limit
        warn("Warning: Simulation was halted after limit reached")
    end
  end
  return event_db
end


text:
# Author: Ali Siahkoohi, alisk@gatech.edu
# Date: September 2020
# Copyright: Georgia Institute of Technology, 2020

export load_experiment

    """
    load_experiment(sim_name::String)

    A function that uses the input variables to the experiment, passed by dictionaty
    `parsed_args`, generated via `ArgParse.jl` package, to determine which experiment
    data to load from `datadir(sim_name)`, where `sim_name` is the simulation name.

    In case there are multiple experiments with the same variable valus, the last
    experiment will be returned

 *Input*:

 - `sim_name`: Simulation name to upload
 - `collect_results`: Simulation name to upload

    """

function load_experiment(parsed_args::Dict; return_path=false)

    # Upload data/sim_name
    sim_name = parsed_args["sim_name"]
    experiments = collect_results(datadir(sim_name), black_list=["Params", "fval"])

    # Drop experiments that miss one or some of the input arguments
    dropmissing!(experiments)

    # Keep experiment with same variables
    for (key, value) in parsed_args
        experiments = experiments[experiments[!, string(key)] .== value, :]
    end

    @assert size(experiments, 1) > 0 "No saved experiments found with such input values"

    println(size(experiments, 1), " experiment(s) found — loading the latest")
    experiment_path = experiments[1, :path]
    Params = wload(experiment_path)["Params"]
    fval = wload(experiment_path)["fval"]

    return_path ? (return Params, fval, experiment_path) : (return Params, fval)

end


text:
using OrdinaryDiffEq, DiffEqSensitivity

function growth(du, u, p, t)
    @. du = p * u * (1 - u)
end
u0 = [0.1]
tspan = (0.0, 2.0)
prob = ODEProblem(growth, u0, tspan, [1.0])
sol = solve(prob, Tsit5(), reltol = 1e-8, abstol = 1e-8)

savetimes = [0.0, 1.0, 1.9]

function f(a)
    _prob = remake(prob,p=[a[1]],saveat=savetimes)
    predicted = solve(_prob, Tsit5(), sensealg=InterpolatingAdjoint(), abstol=1e-12, reltol=1e-12)
    sum(predicted[end])
end

function f2(a)
    _prob = remake(prob,p=[a[1]],saveat=savetimes)
    predicted = solve(_prob, Tsit5(), sensealg=ForwardDiffSensitivity(), abstol=1e-12, reltol=1e-12)
    sum(predicted[end])
end

using Zygote
a = ones(3)
@test Zygote.gradient(f,a)[1][1] ≈ Zygote.gradient(f2,a)[1][1]
@test Zygote.gradient(f,a)[1][2] == Zygote.gradient(f2,a)[1][2] == 0
@test Zygote.gradient(f,a)[1][3] == Zygote.gradient(f2,a)[1][3] == 0


text:
module CuNFFT

using Reexport
@reexport using AbstractNFFTs
using NFFT
using LinearAlgebra
using AbstractFFTs
using CUDA

if CUDA.functional()
  using CUDA.CUSPARSE
  include("implementation.jl")
end

end


text:
begin
	import Pkg
    Pkg.activate(mktempdir())
	Pkg.add(url="https://github.com/brendanjohnharris/Catch22.jl#main")
	using Catch22
	Pkg.add(url="https://github.com/tk3369/YeoJohnsonTrans.jl")
	Pkg.develop(path="C:\\Users\\Brendan\\OneDrive - The University of Sydney (Students)\\Honours\\Code\\ParameterInference.jl")
	Pkg.add(path="C:\\Users\\Brendan\\OneDrive - The University of Sydney (Students)\\Honours\\Code\\ParameterInference.jl")
	using ParameterInference, Plots, StatsBase, Statistics, LinearAlgebra, ParameterInference.NonstationaryProcesses, MultivariateStats, StatsPlots
	gr(); fourseas!(); # Plots initialisation
	nothing
end

nₛ = 1000;
𝜂 = 0.25;
𝛼, 𝛽 = randn(nₛ, 2) |> eachcol;

# ? Case 1
𝑓ₕ₁ = hcat(𝛼, 𝛼, 𝛼, 𝛼, 𝛽)';
𝑓₁ = hcat(𝜂.*[𝛼, 𝛼, 𝛼, 𝛼]..., 𝛽)';
𝛴² = cov(𝑓₁');

𝑓̂₁ = 𝑓₁./std(𝑓₁, dims=2);
𝛴̂² = cov(𝑓̂₁');
𝑃̂ =  fit(MultivariateStats.PCA, 𝑓̂₁).proj; #eigvecs(𝛴̂²);
𝑓̂′₁ = 𝑃̂'*𝑓₁;
cov(𝑓̂′₁, dims=2)
sgn1 = sign(cor(𝑓̂′₁[1, :], 𝛼))
sgn2 = sign(cor(𝑓̂′₁[2, :], 𝛽))

𝛴²ₕ = cov(𝑓ₕ₁');
m1 = fit(MultivariateStats.PCA, 𝑓ₕ₁; maxoutdim=5, pratio=1.0)
𝑃ₕ = m1.proj; # eigvecs(𝛴²ₕ)
𝛬ₕ = sqrt(diagm(m1.prinvars)); # Diagonal(sqrt.(abs.(eigvals(𝛴²ₕ))));
𝑓′′₁ = (𝑃ₕ*inv(𝛬ₕ))'*𝑓₁;
sgn3 = sign(cor(𝑓′′₁[1, :], 𝛼))
sgn4 = sign(cor(𝑓′′₁[2, :], 𝛽))
cov(𝑓′′₁, dims=2)

p = scatter(sgn2.*𝑓̂′₁[2, :], sgn1.*𝑓̂′₁[1, :], ylabel="𝛼", xlabel="𝛽", color=:crimson, markersize=2, label=nothing)
# k = kde((𝑓̂′₁[2, :], 𝑓̂′₁[1, :]); bandwidth=(0.5, 0.5))
# plot!(k, levels=3, color=:crimson)
scatter!(sgn4.*𝑓′′₁[2, :], sgn3.*𝑓′′₁[1, :],  color=ParameterInference.cucumber, markersize=2, label=nothing)
p = scatter!(p, 𝛽, 𝜂.*𝛼,  markersize=2, markercolor=:cornflowerblue, label=nothing)
# k = kde((𝛽, 𝜂.*𝛼); bandwidth=(0.5, 0.5))
# plot!(k, levels=3, color=:cornflowerblue)
covellipse!(p, zeros(2), cov(hcat(𝑓̂′₁[2, :], 𝑓̂′₁[1, :])), n_std=2.25, fillcolor=XYZA(1, 1, 1, 0.0), linecolor=:crimson, linealpha=0.8, label="\$\\textrm{Standardised\\ PCA}: \\rho_{\\alpha} = $(abs(round(cor(𝑓̂′₁[1, :], 𝛼), digits=2))), \\rho_\\beta = $(abs(round(cor(𝑓̂′₁[2, :], 𝛽), digits=2)))\$")
covellipse!(p, zeros(2), cov(hcat(sgn4.*𝑓′′₁[2, :], sgn3.*𝑓′′₁[1, :])), n_std=2.25, fillcolor=XYZA(1, 1, 1, 0.0), linecolor=ParameterInference.cucumber, linealpha=0.8, label="\$\\textrm{Baseline\\ PCA}: \\rho_{\\alpha} = $(abs(round(cor(𝑓′′₁[1, :], 𝛼), digits=2))), \\rho_\\beta = $(abs(round(cor(𝑓′′₁[2, :], 𝛽), digits=2)))\$")
covellipse!(p, zeros(2), cov(hcat(𝛽, 𝜂.*𝛼)), n_std=2.25, fillcolor=XYZA(1, 1, 1, 0.0), linecolor=:cornflowerblue, linealpha=0.8, size=(500, 400), xlims=(-3.6, 3.6), ylims=(-2.4, 2.4), label="True parameters", legend=true, foreground_color_legend=nothing, background_color_legend=nothing)
savefig("../figures/$(rand(UInt16)).svg")

# ? Case 2

# ╔═╡ 46e7d2a0-8b62-4112-92df-74538f5f16d8
𝑓ₕ₂ = hcat(𝛼+𝛽/4, 𝛼+𝛽/4, 𝛼+𝛽/4, 𝛼+𝛽/4, 𝛽)';
𝑓₂ = hcat((𝜂.*[𝛼, 𝛼, 𝛼, 𝛼] .+ [𝛽, 𝛽, 𝛽, 𝛽]./4)..., 𝛽)';

# ╔═╡ feeaa835-28e0-4a24-905b-53347560ec9a
𝑓̂₂ = 𝑓₂./std(𝑓₂, dims=2);

# ╔═╡ 2358f79f-f44d-4e93-b4da-75e7be715882
𝛴̂²₂ = cov(𝑓̂₂')

# ╔═╡ d9fca85a-63b6-4ab5-ab0e-6cd5e68c4e06
𝑃̂₂ = fit(MultivariateStats.PCA, 𝑓̂₂).proj; # eigvecs(𝛴̂²₂);

# ╔═╡ d700b4ae-8bd5-4437-b515-c21d3da70543
𝑓̂′₂ = 𝑃̂₂'*𝑓̂₂;
cov(𝑓̂′₂, dims=2)

sgn2 = sign(cor(𝑓̂′₂[2, :], 𝛼))
sgn1 = sign(cor(𝑓̂′₂[1, :], 𝛽))
std(𝑓̂′₂[1, :])./std(𝑓̂′₂[2, :])
std(𝛽)./std(𝜂.*𝛼 .+ 𝛽./4)
# ╔═╡ 06e226a6-0f0f-4ef6-b916-04d51a3789c8
scatter(sgn1.*𝑓̂′₂[1, :], sgn2.*𝑓̂′₂[2, :], aspect_ratio=:equal, ylabel="𝛼′", xlabel="𝛽′", title="Standardised PCA estimate")

# ╔═╡ ea513127-e9bb-40a7-9ec0-2beee8fa6179
scatter!(𝛽, 𝜂.*𝛼 .+ 𝛽./4, markersize=2, markercolor=:gray)

# ╔═╡ f946e3d7-15b9-45a9-951b-b3d77fc589f8
𝛴²ₕ₂ = cov(𝑓ₕ₂');

# ╔═╡ 652b998d-a6b0-4bb7-908e-c3dc7a82764c
m2 = fit(MultivariateStats.PCA, 𝑓ₕ₂);

# ╔═╡ d192149c-1ef4-4761-9e78-879a97abee57
𝑃ₕ₂ = m2.proj; #eigvecs(𝛴²ₕ₂);

𝛬ₕ₂ = sqrt(diagm(m2.prinvars)); # Diagonal(sqrt.(abs.(eigvals(𝛴²ₕ₂))))

𝑓′′₂ = (𝑃ₕ₂*inv(𝛬ₕ₂))'*𝑓₂;

# Now pca againt o get underlying params:
𝑓′′₂ = fit(MultivariateStats.PCA, 𝑓′′₂).proj'*𝑓′′₂;
sgn3 = sign(cor(𝑓′′₂[1, :], 𝛽))
sgn4 = sign(cor(𝑓′′₂[2, :], 𝛼))
scatter(sgn4.*𝑓′′₂[2, :], sgn3.*𝑓′′₂[1, :], aspect_ratio=:equal, ylabel="𝛼′", xlabel="𝛽′");

scatter!(𝛽, 𝜂.*𝛼 .+ 𝛽./4, markersize=2, markercolor=:gray)

p = scatter(sgn1.*𝑓̂′₂[1, :], sgn2.*𝑓̂′₂[2, :], ylabel="𝛼", xlabel="𝛽", color=:crimson, markersize=2, label=nothing)
# k = kde((𝑓̂′₁[2, :], 𝑓̂′₁[1, :]); bandwidth=(0.5, 0.5))
# plot!(k, levels=3, color=:crimson)
scatter!(sgn3.*𝑓′′₂[1, :], sgn4.*𝑓′′₂[2, :], color=ParameterInference.cucumber, markersize=2, label=nothing)
p = scatter!(p, 𝛽, 𝜂.*𝛼, markersize=2, markercolor=:cornflowerblue, label=nothing)
# k = kde((𝛽, 𝜂.*𝛼); bandwidth=(0.5, 0.5))
# plot!(k, levels=3, color=:cornflowerblue)
covellipse!(p, zeros(2), cov(hcat(sgn1.*𝑓̂′₂[1, :], sgn2.*𝑓̂′₂[2, :])), n_std=2.25, fillcolor=XYZA(1, 1, 1, 0.0), linecolor=:crimson, linealpha=0.8, label="\$\\textrm{Standardised\\ PCA}: \\rho_{\\alpha} = $(abs(round(cor(𝑓̂′₂[2, :], 𝛼), digits=2))), \\rho_\\beta = $(abs(round(cor(𝑓̂′₂[1, :], 𝛽), digits=2)))\$")
covellipse!(p, zeros(2), cov(hcat(sgn3.*𝑓′′₂[1, :], sgn4.*𝑓′′₂[2, :])), n_std=2.25, fillcolor=XYZA(1, 1, 1, 0.0), linecolor=ParameterInference.cucumber, linealpha=0.8, label="\$\\textrm{Baseline\\ PCA}: \\rho_{\\alpha} = $(round(cor(𝑓′′₂[2, :], 𝛼), digits=2)), \\rho_\\beta = $(round(cor(𝑓′′₂[1, :], 𝛽), digits=2))\$")
covellipse!(p, zeros(2), cov(hcat(𝛽, 𝜂.*𝛼)), n_std=2.25, fillcolor=XYZA(1, 1, 1, 0.0), linecolor=:cornflowerblue, linealpha=0.8, size=(500, 400), xlims=(-6, 6), ylims=(-2.4, 2.4), label="True parameters", legend=true, foreground_color_legend=nothing, background_color_legend=nothing)
savefig("../figures/$(rand(UInt16)).svg")

# ? Case 3
𝜂₃ = 0.1;
𝑓ₕ₃ = hcat(repeat([𝛼.+𝜂₃.*randn(nₛ)],3)..., 𝛽.+𝜂₃.*randn(nₛ), (𝜂₃).*randn(nₛ))';
𝑓₀ = hcat(repeat([𝜂₃.*randn(nₛ)], 3)..., 𝜂₃.*randn(nₛ), 𝜂₃.*randn(nₛ))';
𝑓₃ = hcat(repeat([𝜂.*𝛼 .+ 𝜂₃.*randn(nₛ)], 3)..., 𝛽.+𝜂₃.*randn(nₛ), 𝜂₃.*randn(nₛ))';

# * Standardised PCA

# ╔═╡ a4c8677f-6bdc-4bfe-adcc-3130f45b7aaa
𝑓̂₃ = 𝑓₃./std(𝑓₃, dims=2);

# ╔═╡ 5ad01546-7daa-441f-bf2c-60e2a6e6f0c7
𝛴̂²₃ = cov(𝑓̂₃');

# ╔═╡ 36039419-6df3-45f2-ae2c-77cf6e3b405d
𝑃̂₃ = fit(MultivariateStats.PCA, 𝑓̂₃).proj;# eigvecs(𝛴̂²₃);

# ╔═╡ d046eb57-d4cb-4eda-9a81-00a234af1d57
𝑓̂′₃ = 𝑃̂₃'*𝑓̂₃;
sgn2 = sign(cor(𝑓̂′₃[2, :], 𝛽))
sgn1 = sign(cor(𝑓̂′₃[1, :], 𝛼))

cov(𝑓̂′₃, dims=2) # So actually 3 dim

# ╔═╡ 4e9b0004-da7a-4cb3-afe9-c3365d37471a
scatter(𝑓̂′₃[1, :], 𝑓̂′₃[2, :], aspect_ratio=:equal, ylabel="𝛼′", xlabel="𝛽′", title="Standardised PCA estimate")
cor(𝑓̂′₃[1, :], 𝛼)
std(𝑓̂′₃[1, :])./std(𝑓̂′₃[2, :])
std(𝛽)./std(𝜂.*𝛼)
# ╔═╡ 5e2281e9-360b-41dd-9f4f-38e14d8aea86
scatter!(𝛽, 𝜂.*𝛼, markersize=2, markercolor=:gray)

# ╔═╡ bcd85e11-8bf0-4968-a528-8b56ce14db5b
# * Baseline PCA

# ╔═╡ f771887a-77f0-4a4d-9e4c-7e988201197f
𝛴²ₕ₃ = cov(𝑓ₕ₃');

# ╔═╡ b636df98-f935-46b4-bb08-be96cee18884
𝛴²₀ = cov(𝑓₀');

# ╔═╡ c43d6053-4eeb-450a-bca4-bd5d58b179c1
𝛴²₃ = cov(𝑓₃');

# ╔═╡ beafaf38-3753-4e9f-bf6d-69b98439fe18
m3 = fit(MultivariateStats.PCA, 𝑓ₕ₃)

# ╔═╡ fd00b9e4-f2fc-4571-8596-b1e53f8d7286
𝑃ₕ₃ = m3.proj;#eigvecs(𝛴²ₕ₃);

# ╔═╡ 9d13988d-e625-42b5-8e89-86d519027912
𝛬ₕ₃ = sqrt(diagm(m3.prinvars)); #Diagonal(sqrt.(abs.(eigvals(𝛴²ₕ₃))))

# ╔═╡ 4f5c167d-dbde-4fed-bd2d-63b0100e2d1a
𝛴²ₕ′ = diagm(m3.prinvars)

# ╔═╡ 82841686-b078-4ba2-b5ce-d748070437d2
𝛴²₀′ = cov(𝑃ₕ₃'*𝑓₀, dims=2)

# ╔═╡ 8dded5c7-f4ea-4ee5-b0b9-59b21c7f3097
𝛴²′ = cov(𝑃ₕ₃'*𝑓₃, dims=2)

# ╔═╡ cb77a71b-7cba-489d-a51b-690cdd281a4c
𝑉 = I - inv(sqrt(Diagonal(𝛴²′)))*sqrt(Diagonal(𝛴²₀′));

# ╔═╡ b00f965a-18a1-4ccf-85a9-ae79d588cab6
𝑆 = 𝑉/(𝛬ₕ₃ - sqrt(Diagonal(𝛴²₀′)))

# ╔═╡ 0ae23619-0ce3-4e32-9a11-f9a444235c70
𝑓′′₃ = 𝑆'*𝑃ₕ₃'*𝑓₃;
# Now pca againt o get underlying params:
𝑓′′₃ = fit(MultivariateStats.PCA, 𝑓′′₃).proj'*𝑓′′₃;
cov(𝑓′′₃, dims=2)
sgn3 = sign(cor(𝑓′′₃[1, :], 𝛽))
sgn4 = sign(cor(𝑓′′₃[2, :], 𝛼))

# ╔═╡ ddb9528f-b57c-4333-9347-8a56aa6890c4
#f′′₃ = (𝑓′′₃./std(𝑓′′₃, dims=2))*((sqrt(Diagonal(𝛴²′)) - sqrt(Diagonal(𝛴²₀′)))/())

# ╔═╡ c3067e6d-dfaf-4938-8bbf-4b51514fd286
𝑓′′₀ = (𝑃ₕ₃*inv(𝛬ₕ₃)*𝑉)'*𝛴²₀*(𝑃ₕ₃*inv(𝛬ₕ₃)*𝑉); # Not 0 because the transform is a strong approximation

# ╔═╡ a2f4c245-031b-4f1c-9f44-b3947cbc9795
cov(𝑓′′₃, dims=2)

# ╔═╡ ad7aac3b-303a-4ca6-a44f-f11cba8d115f
𝑓′′₀

# ╔═╡ dbdcd6c4-102d-41cc-b218-100f74355b85
scatter(𝑓′′₃[2, :], .-𝑓′′₃[1, :], aspect_ratio=:equal, ylabel="𝛼′", xlabel="𝛽′", title="Baseline PCA estimate", left_margin=5Plots.mm);

# ╔═╡ 9d6334c1-7107-4fbb-b383-9a37fc56de74
scatter!(.-𝛽, .-𝜂.*𝛼, markersize=1, markercolor=:gray)

# ╔═╡ 07e9a417-5adb-4182-a98d-fd618feaf675
scatter(𝛽, 𝑓′′₃[2, :]) # Compairwise to true parameters

# ╔═╡ 008512c1-3a44-42b4-8695-2a9f4d88ef01
scatter(𝛼, 𝑓′′₃[1, :]) # Compairwise to true parameters

p = scatter(sgn2.*𝑓̂′₃[2, :], sgn1.*𝑓̂′₃[1, :], ylabel="𝛼", xlabel="𝛽", color=:crimson, markersize=2, label=nothing)
# k = kde((𝑓̂′₁[2, :], 𝑓̂′₁[1, :]); bandwidth=(0.5, 0.5))
# plot!(k, levels=3, color=:crimson)
scatter!(sgn3.*𝑓′′₃[1, :], sgn4.*𝑓′′₃[2, :], color=ParameterInference.cucumber, markersize=2, label=nothing)
p = scatter!(p, 𝛽, 𝜂.*𝛼, markersize=2, markercolor=:cornflowerblue, label=nothing)
# k = kde((𝛽, 𝜂.*𝛼); bandwidth=(0.5, 0.5))
# plot!(k, levels=3, color=:cornflowerblue)
covellipse!(p, zeros(2), cov(hcat(sgn2.*𝑓̂′₃[2, :], sgn1.*𝑓̂′₃[1, :])), n_std=2.25, fillcolor=XYZA(1, 1, 1, 0.0), linecolor=:crimson, linealpha=0.8, label="\$\\textrm{Standardised\\ PCA}: \\rho_{\\alpha} = $(abs(round(cor(𝑓̂′₃[1, :], 𝛼), digits=2))), \\rho_\\beta = $(abs(round(cor(𝑓̂′₃[2, :], 𝛽), digits=2)))\$")
covellipse!(p, zeros(2), cov(hcat(sgn3.*𝑓′′₃[1, :], sgn4.*𝑓′′₃[2, :])), n_std=2.25, fillcolor=XYZA(1, 1, 1, 0.0), linecolor=ParameterInference.cucumber, linealpha=0.8, label="\$\\textrm{Baseline\\ PCA}: \\rho_{\\alpha} = $(round(cor(𝑓′′₃[2, :], 𝛼), digits=2)), \\rho_\\beta = $(round(cor(𝑓′′₃[1, :], 𝛽), digits=2))\$")
covellipse!(p, zeros(2), cov(hcat(𝛽, 𝜂.*𝛼)), n_std=2.25, fillcolor=XYZA(1, 1, 1, 0.0), linecolor=:cornflowerblue, linealpha=0.8, size=(500, 400), xlims=(-3.6, 3.6), ylims=(-4.8, 4.8), label="True parameters", legend=true, foreground_color_legend=nothing, background_color_legend=nothing)
savefig("../figures/$(rand(UInt16)).svg")


text:
"""
    AnisotropicBiharmonicDiffusivity{FT, KH, KV}

Parameters for anisotropic biharmonic diffusivity models.
"""
struct AnisotropicBiharmonicDiffusivity{FT, KH, KV} <: TensorDiffusivity{FT}
    νh :: FT
    νv :: FT
    κh :: KH
    κv :: KV
    function AnisotropicBiharmonicDiffusivity{FT}(νh, νv, κh, κv) where FT
        return new{FT, typeof(κh), typeof(κv)}(νh, νv, convert_diffusivity(FT, κh), convert_diffusivity(FT, κv))
    end
end

"""
    AnisotropicBiharmonicDiffusivity(; νh, νv, κh, κv)

Returns parameters for a fourth-order, anisotropic biharmonic diffusivity closure with
constant horizontal and vertical biharmonic viscosities `νh`, `νv` and constant horizontal
and vertical tracer biharmonic diffusivities `κh`, `κv`.
`κh` and `κv` may be `NamedTuple`s with fields corresponding
to each tracer, or a single number to be a applied to all tracers.
The tracer flux divergence associated with an anisotropic biharmonic diffusivity is, for example

```math
    ∂ᵢ κᵢⱼ ∂ⱼc = (κh ∇h⁴ + κv ∂z⁴) c
```

"""
AnisotropicBiharmonicDiffusivity(FT=Float64; νh=0, νv=0, κh=0, κv=0) =
    AnisotropicBiharmonicDiffusivity{FT}(νh, νv, κh, κv)

function with_tracers(tracers, closure::AnisotropicBiharmonicDiffusivity{FT}) where FT
    κh = tracer_diffusivities(tracers, closure.κh)
    κv = tracer_diffusivities(tracers, closure.κv)
    return AnisotropicBiharmonicDiffusivity{FT}(closure.νh, closure.νv, κh, κv)
end

calculate_diffusivities!(K, arch, grid, closure::AnisotropicBiharmonicDiffusivity, args...) = nothing

@inline ∂ⱼ_2ν_Σ₁ⱼ(i, j, k, grid, closure::AnisotropicBiharmonicDiffusivity, U, args...) = (
    - closure.νh * ∇h⁴_fca(i, j, k, grid, U.u)
    - closure.νv * ∂z⁴_aac(i, j, k, grid, U.u)
    )

@inline ∂ⱼ_2ν_Σ₂ⱼ(i, j, k, grid, closure::AnisotropicBiharmonicDiffusivity, U, args...) = (
    - closure.νh * ∇h⁴_cfa(i, j, k, grid, U.v)
    - closure.νv * ∂z⁴_aac(i, j, k, grid, U.v)
    )

@inline ∂ⱼ_2ν_Σ₃ⱼ(i, j, k, grid, closure::AnisotropicBiharmonicDiffusivity, U, args...) = (
    - closure.νh * ∇h⁴_cca(i, j, k, grid, U.w)
    - closure.νv * ∂z⁴_aaf(i, j, k, grid, U.w)
    )

@inline function ∇_κ_∇c(i, j, k, grid, closure::AnisotropicBiharmonicDiffusivity,
                        c, ::Val{tracer_index}, args...) where tracer_index

    @inbounds κh = closure.κh[tracer_index]
    @inbounds κv = closure.κv[tracer_index]

    return (- κh * ∇h⁴_cca(i, j, k, grid, c)
            - κv * ∂z⁴_aac(i, j, k, grid, c)
           )
end