text:
"""
    $(SIGNATURES)

Plot absolute values of variables and measurement surrounding fsym factor.
"""
function plotKDEofnc(fgl::AbstractDFG,
                     fsym::Symbol;
                     solveKey::Symbol=:default,
                     marg=nothing,
                     N::Int=100 )
  #
  fnc = nothing
  if haskey(fgl.fIDs, fsym)
    fnc = getfnctype( fgl, fgl.fIDs[fsym] )
  else
    error("fIDs doesn't have $(fsym)")
  end
  p = kde!( getSample(fnc, N)[1]  )
  # mmarg = length(marg) > 0 ? marg : collect(1:Ndim(p))
  plotKDE(fgl, lsf(fgl, fsym), solveKey=solveKey, addt=[p], marg=marg)
end

"""
    $(SIGNATURES)

Try plot relative values of variables and measurement surrounding fsym factor.
"""
function plotKDEresiduals(fgl::AbstractDFG,
                          fsym::Symbol;
                          N::Int=100,
                          levels::Int=3,
                          dims=nothing,
                          fill=false  )
  #
  # COLORS = ["black";"red";"green";"blue";"cyan";"deepskyblue"]
  fnc = getfnctype( fgl, fgl.fIDs[fsym] )
  # @show sxi = lsf(fgl, fsym)[1]
  # @show sxj = lsf(fgl, fsym)[2]
  fct = getFactor(fgl, fsym)
  @show sxi = getData(fct).fncargvID[1]
  @show sxj = getData(fct).fncargvID[2]
  xi = getVal(fgl, sxi)
  xj = getVal(fgl, sxj)
  measM = getSample(fnc, N)
  meas = length(measM) == 1 ? (0*measM[1], ) : (0*measM[1], measM[2])
  @show size(meas[1])
  d = size(measM[1],1)
  RES = zeros(d,N)
  for i in 1:N
    res = zeros(d)
    fnc(res, i, meas, xi, xj)
    RES[:,i] = res
    if length(measM) > 1
      if measM[2][i] == 0
        RES[:,i] = 0.5*randn(d)
      end
    end
  end
  pR = kde!(RES)
  pM = kde!(measM[1])
  fncvar = getfnctype(fct)

  COLORS = getColorsByLength()
  plotKDE([pR;pM], c=COLORS[1:2], dims=dims,levels=3, legend=["residual";"model"], fill=fill, title=string(fsym, ", ", fncvar))
end

# """
#     $(SIGNATURES)

# Plot the product of priors and incoming upward messages for variable in clique.

# plotPriorsAtCliq(tree, :x2, :x1[, marg=[1;2]])
# """
# function plotPriorsAtCliq(tree::BayesTree,
#                           lb::Symbol,
#                           cllb::Symbol;
#                           dims::Vector{Int}=Int[],
#                           levels::Int=1,
#                           fill::Bool=false  )
#   #
#   # COLORS = ["black";"red";"green";"blue";"cyan";"deepskyblue"]

#   cliq = whichCliq(tree, cllb)
#   cliqprs = cliq.attributes["debug"].priorprods[1]

#   vidx = 1
#   for lbel in cliqprs.lbls
#     if lbel == lb
#       break;
#     else
#       vidx += 1
#     end
#   end
#   dims = length(dims)>0 ? dims : collect(1:size(cliqprs.prods[vidx].prev,1))

#   arr = BallTreeDensity[]
#   push!(arr, marginal(kde!(cliqprs.prods[vidx].prev), dims)  )
#   push!(arr, marginal(kde!(cliqprs.prods[vidx].product), dims)  )
#   len = length(cliqprs.prods[vidx].potentials)
#   lg = String["p";"n"]
#   i=0
#   for pot in cliqprs.prods[vidx].potentials
#     push!(arr, marginal(pot, dims)  )
#     i+=1
#     push!(lg, cliqprs.prods[vidx].potentialfac[i])
#   end

#   COLORS = getColorsByLength(len+2)
#   # cc = COLORS[1:(len+2)]
#   plotKDE(arr, c=COLORS, legend=lg, levels=levels, fill=fill );
# end

# function plotMCMC(treel::BayesTree,
#                   lbll::Symbol;
#                   delay::Int=200,
#                   show::Bool=true,
#                   w=20cm, h=15cm,
#                   levels::Int=1,
#                   dims=nothing  )
#   #
#   cliq = whichCliq(treel, string(lbll))
#   cliqdbg = cliq.attributes["debug"]

#   vidx = 1
#   for lb in cliqdbg.mcmc[1].lbls
#     if lb == lbll
#       break;
#     else
#       vidx += 1
#     end
#   end

#   tmpfilepath = joinpath(dirname(@__FILE__),"tmpimgs")
#   ARR = BallTreeDensity[]
#   COLORS = getColorsByLength()
#   # COLORS = ["black";"red";"green";"blue";"cyan";"deepskyblue";"magenta"]
#   for i in 1:length(cliqdbg.mcmc)
#     ppr = kde!(cliqdbg.mcmc[i].prods[vidx].prev)
#     ppp = kde!(cliqdbg.mcmc[i].prods[vidx].product)
#     ARR = [ARR;ppr;ppr;cliqdbg.mcmc[i].prods[vidx].potentials]
#   end
#   rangeV = getKDERange(ARR)
#   ppp = nothing
#   for i in 1:length(cliqdbg.mcmc)
#     ppr = kde!(cliqdbg.mcmc[i].prods[vidx].prev)
#     ppp = kde!(cliqdbg.mcmc[i].prods[vidx].product)
#     arr = [ppr;ppp;cliqdbg.mcmc[i].prods[vidx].potentials]
#     len = length(cliqdbg.mcmc[i].prods[vidx].potentials)
#     lg = String["p";"n";cliqdbg.mcmc[i].prods[vidx].potentialfac] #map(string, 1:len)]
#     COLORS_l = getColorsByLength(len+2)
#     cc = plotKDE(arr, c=COLORS_l, legend=lg, levels=levels, fill=true, axis=rangeV, dims=dims );
#     Gadfly.draw(PNG(joinpath(tmpfilepath,"$(string(lbll))mcmc$(i).png"),w,h),cc)
#   end
#   # draw initial and final result
#   pp0 = kde!(cliqdbg.mcmc[1].prods[vidx].prev)
#   i = 0
#   cc = plotKDE([pp0], c=[COLORS[1]], legend=["0"], levels=levels, fill=true, axis=rangeV, dims=dims );
#   Gadfly.draw(PNG(joinpath(tmpfilepath,"$(string(lbll))mcmc$(i).png"),w,h),cc)

#   i = length(cliqdbg.mcmc)+1
#   cc = plotKDE([pp0;ppp], c=COLORS[1:2], legend=["0";"n"], levels=levels, fill=true, axis=rangeV, dims=dims );
#   Gadfly.draw(PNG(joinpath(tmpfilepath,"$(string(lbll))mcmc$(i).png"),w,h),cc)
#   # generate output
#   @warn "Maintenance required, run manually:  \"convert -delay $(delay) $(tmpfilepath)/$(string(lbll))mcmc*.png $(tmpfilepath)/$(string(lbll))mcmc.gif\""
#   # run(`convert -delay $(delay) $(tmpfilepath)/$(string(lbll))mcmc*.png $(tmpfilepath)/$(string(lbll))mcmc.gif`)
#   # !show ? nothing : (@async run(`eog $(tmpfilepath)/$(string(lbll))mcmc.gif`) )
#   # return "$(tmpfilepath)/$(string(lbll))mcmc.gif"
# end

# function plotMCMCDebug(cliq; offs=2.0)
#     println("$(cliq.attributes["label"])")
#     cliqDbg = cliq.attributes["data"].debug
#     MCd = 100.0/length(cliqDbg.mcmc)
#     MCo=0.0
#     minmax=[99999,-99999,-10,-99999.0]
#     for mc in cliqDbg.mcmc
#         drawOneMC!(mc, minmax, MCo, offs=offs)
#         MCo+=MCd
#     end

#     n = 3
#     y = linspace(minmax[1], minmax[2], n)
#     x = linspace(minmax[3],minmax[4],n)
#     xgrid = repmat(x',n,1)
#     ygrid = repmat(y,1,n)
#     z = zeros(n,n)

#     for i in 1:length(cliqDbg.mcmc[1].prods)
#       surf(xgrid,ygrid,z-(i-1)*offs,alpha=0.04, linewidth=0.0)
#     end
# end


text:
@testset "RepulsiveCore" begin

##
@info("--------------- Testing RepulsiveCore Implementation ---------------")

at = bulk(:W, cubic=true) * 3
rattle!(at, 0.03)
r0 = rnn(:W)
z = atomic_number(:W)
trans = PolyTransform(2, r0)
fcut = PolyCutoff2s(2, 0.5*r0, 1.95*r0)
B = PolyPairBasis(:W, 10, trans, fcut)
coeffs = randcoeffs(B)
pot = combine(B, coeffs)

## try out the repulsive potential
Vfit = pot

ri = 2.1
@show @D Vfit(ri)
if (@D Vfit(ri)) > 0
   Vfit = PolyPairPot(- Vfit.coeffs, Vfit.J, Vfit.zlist, Vfit.bidx0)
end
@show @D Vfit(ri)
e0 = Vfit(ri) - 1.0
Vrep = PolyPairPots.Repulsion.RepulsiveCore(Vfit, ri)

rout = range(ri+1e-15, 4.0, length=100)
println(@test all(Vfit(r) == Vrep(r,z,z) for r in rout))
rin = range(0.5, ri, length=100)
println(@test all(Vrep.Vin[1](r) == Vrep(r,z,z) for r in rin))

@info("JuLIP FD test")
println(@test JuLIP.Testing.fdtest(Vrep, at))

@info("check scaling")
println(@test energy(Vfit, at) ≈ energy(Vrep, at))

##

@info("--------------- Multi-Species RepulsiveCore ---------------")

at = bulk(:W, cubic=true) * 3
at.Z[2:3:end] .= atomic_number(:Fe)
rattle!(at, 0.03)
r0 = rnn(:W)
trans = PolyTransform(2, r0)
fcut = PolyCutoff2s(2, 0.5*r0, 1.95*r0)
B = PolyPairBasis([:W, :Fe], 10, trans, fcut)
coeffs = randcoeffs(B)
pot = combine(B, coeffs)

## try out the repulsive potential
Vfit = pot

ri = 2.1
z1 = 26
z2 = 74
@show @D Vfit(ri, z1, z2)
e0 = min(Vfit(ri, z1, z2), Vfit(ri, z1, z1), Vfit(ri, z2, z2)) - 1.0
Vrep = PolyPairPots.Repulsion.RepulsiveCore(Vfit, ri, e0)

for (z, z0) in zip([z1, z1, z2], [z1, z2, z2])
   i, i0 = JuLIP.Potentials.z2i(Vfit, z), JuLIP.Potentials.z2i(Vfit, z0)
   rout = range(ri+1e-15, 4.0, length=100)
   println(@test all(Vfit(r, z, z0) == Vrep(r, z, z0) for r in rout))
   rin = range(0.5, ri, length=100)
   println(@test all(Vrep.Vin[i,i0](r) == Vrep(r, z, z0) for r in rin))
end

@info("JuLIP FD test")
println(@test JuLIP.Testing.fdtest(Vrep, at))

@info("check scaling")
println(@test energy(Vfit, at) ≈ energy(Vrep, at))

@info("check FIO")
println(@test (Vrep == decode_dict(Dict(Vrep))))
fname = tempname()
save_json(fname, Dict(Vrep))
D1 = load_json(fname)
rm(fname)
println(@test (Vrep == decode_dict(D1)))

##

end

ext: jl


text:
"""
    LocalGeometry

The necessary local metric information defined at each node.
"""
struct LocalGeometry{C, FT, Mxξ, Mξx}
    "Coordinates of the current point"
    coordinates::C
    "Jacobian determinant of the transformation `ξ` to `x`"
    J::FT
    "Metric terms: `J` multiplied by the quadrature weights"
    WJ::FT
    "Partial derivatives of the map from `ξ` to `x`: `∂x∂ξ[i,j]` is ∂xⁱ/∂ξʲ"
    ∂x∂ξ::Mxξ
    "Partial derivatives of the map from `x` to `ξ`: `∂ξ∂x[i,j]` is ∂ξⁱ/∂xʲ"
    ∂ξ∂x::Mξx
end

const LocalGeometry1D =
    LocalGeometry{C, FT, M} where {C <: Abstract1DPoint{FT}} where {FT, M}

const LocalGeometry2D =
    LocalGeometry{C, FT, M} where {C <: Abstract2DPoint{FT}} where {FT, M}

const LocalGeometry3D =
    LocalGeometry{C, FT, M} where {C <: Abstract3DPoint{FT}} where {FT, M}

"""
    SurfaceGeometry

The necessary local metric information defined at each node on each surface.
"""
struct SurfaceGeometry{FT, N}
    "surface Jacobian determinant, multiplied by the surface quadrature weight"
    sWJ::FT
    "surface outward pointing normal vector"
    normal::N
end

undertype(::Type{LocalGeometry{C, FT, Mxξ, Mξx}}) where {C, FT, Mxξ, Mξx} = FT
undertype(::Type{SurfaceGeometry{FT, N}}) where {FT, N} = FT


text:
# breusch_godfrey.jl
# Breusch-Godfrey test for autocorrelation
#
# Copyright (C) 2017   Benjamin Born
#
# Permission is hereby granted, free of charge, to any person obtaining
# a copy of this software and associated documentation files (the
# "Software"), to deal in the Software without restriction, including
# without limitation the rights to use, copy, modify, merge, publish,
# distribute, sublicense, and/or sell copies of the Software, and to
# permit persons to whom the Software is furnished to do so, subject to
# the following conditions:
#
# The above copyright notice and this permission notice shall be
# included in all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
# LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
# OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
# WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

export BreuschGodfreyTest

struct BreuschGodfreyTest <: HypothesisTest
    n::Int              # number of observations
    lag::Int            # number of lags in test statistic
    BG::Float64         # test statistic
end

"""
    BreuschGodfreyTest(X, e, lag, start0 = true)

Compute the Breusch-Godfrey test for serial correlation in the residuals
of a regression model.

`X` is the matrix of regressors from the original model and `e` the vector of residuals.
`lag` determines the number of lagged residuals included in the auxiliary regression.
Set `start0` to specify how the starting values for the lagged residuals are handled.
`start0 = true` (default) sets them to zero (as in Godfrey, 1978); `start0 = false`
uses the first `lag` residuals as starting values, i.e. shortening the sample by `lag`.

# External links

  * [Breusch-Godfrey test on Wikipedia](https://en.wikipedia.org/wiki/Breusch–Godfrey_test)
"""
function BreuschGodfreyTest(xmat::AbstractArray{T}, e::AbstractVector{T},
                            lag::Int, start0::Bool=true) where T<:Real
    n = size(e,1)
    elag = zeros(Float64,n,lag)
    for ii = 1:lag  # construct lagged residuals
        elag[ii+1:end,ii] = e[1:end-ii]
    end

    offset = start0 ? 0 : lag

    regmat = [xmat[offset+1:end,:] elag[offset+1:end,:]]
    regcoeff = regmat\e[offset+1:end]
    resid = e[offset+1:end] - regmat*regcoeff

    rsq = 1 - dot(resid,resid)/dot(e[offset+1:end],e[offset+1:end]) # uncentered R^2
    BG = (n-offset)*rsq
    BreuschGodfreyTest(n-offset,lag,BG)
end

testname(::BreuschGodfreyTest) = "Breusch-Godfrey autocorrelation test"
population_param_of_interest(x::BreuschGodfreyTest) =
    ("coefficients on lagged residuals up to lag p", "all zero", NaN)
default_tail(test::BreuschGodfreyTest) = :right

function show_params(io::IO, x::BreuschGodfreyTest, ident)
    println(io, ident, "number of observations:         ", x.n)
    println(io, ident, "number of lags:                 ", x.lag)
    println(io, ident, "T*R^2 statistic:                ", x.BG)
end

pvalue(x::BreuschGodfreyTest) = pvalue(Chisq(x.lag), x.BG; tail=:right)


text:
######### StanOptimize example  ###########

using StanOptimize

bernoulli_model = "
data {
  int<lower=1> N;
  int<lower=0,upper=1> y[N];
}
parameters {
  real<lower=0,upper=1> theta;
}
model {
  theta ~ beta(1,1);
  y ~ bernoulli(theta);
}
";

bernoulli_data = Dict("N" => 10, "y" => [0, 1, 0, 1, 0, 0, 0, 0, 0, 1])

stanmodel = OptimizeModel("bernoulli",  bernoulli_model;
  tmpdir = joinpath(@__DIR__, "tmp"));

(sample_path, log_path) = stan_sample(stanmodel, data=bernoulli_data);

if sample_path !== Nothing
  optim1, cnames = read_optimize(stanmodel)
  println()
  display(optim1)
  println()
end

# Same with saved iterations
stanmodel = OptimizeModel("bernoulli", bernoulli_model;
  method = StanOptimize.Optimize(save_iterations = true),
  tmpdir = joinpath(@__DIR__, "tmp"));

(sample_path, log_path)  = stan_optimize(stanmodel, data=bernoulli_data);

if sample_path !== Nothing
  optim2, cnames = read_optimize(stanmodel)
  println()
  display(optim2)
  println()
end


text:
using Distributions
using PyPlot

include("genw_sparse.jl")
include("run_target.jl")
include("run_train.jl")
include("run_test.jl")

# connectivity matrix
N = 200;
p = 0.3;
g = 1.0;

# neuron dynamics
b  = 1/2;

# training variables
train_duration = 100.;
nloop          = 10; #30
lambda         = 0.5; # 1, 10
learn_every    = 1.0;

# training param
stim_on         = 100.;
stim_off        = 110.;
train_time      = stim_off + train_duration;
test_time       = train_time + 100.;

dt = 0.1;

#---------- Choose target type ----------#
# target_type = "periodic";
target_type = "ratemodel";
# target_type = "ou";

#---------- Generate target trajectories ----------#
# generate M
M = g*genw_sparse(N,0,1,p)/sqrt(N*p);

# generate stim
stim = 2*(2*rand(N)-1);

# generate target patterns
net_param = [N, p, g, b, dt, test_time];
utarg = run_target(M,net_param,target_type);

#---------- Train recurrent connectivity ----------#
# Initialize variables
x0 = 2*pi*rand(N);
r0 = zeros(N);

# Run training
net_param   = [x0, r0, b];
train_param = [stim, lambda, learn_every, nloop, target_type];
time_param  = [stim_on, stim_off, train_time]
M = run_train(M,net_param,train_param,time_param,dt,utarg);

#---------- Test ----------#
# Initialize variables
x0 = 2*pi*rand(N);
r0 = zeros(N);
z0 = zeros(N);

net_param = [x0, r0, b, stim];
time_param = [stim_on, stim_off, test_time];
utest = run_test(M,net_param,time_param,dt);

#---------- Plot results ----------#
tvec = collect(dt: dt : test_time);

figure(figsize=(10,5))
nrow = 2;
ncol = 2;

for i = 1:nrow*ncol
    subplot(2,2,i)
    axvspan(10*stim_on,10*stim_off,alpha=0.5,color="dodgerblue")
    axvspan(10*stim_off,10*train_time,alpha=0.5,color="gray")

    ton = Int(stim_off/dt);
    toff = Int(train_time/dt);
    plot(tvec[ton:toff]*10,utarg[ton:toff,i],linewidth=3,color="black",label="target")
    plot(tvec*10,utest[:,i],linewidth=2,color="red",alpha=1,label="actual")

    xlim([500, 2500])
    xlabel("time (ms)",fontsize=15)
    ylabel("synaptic current",fontsize=15)

    if i == 4
      legend(loc=1,frameon=false,fontsize=15)
    end
end
tight_layout()


text:
"""
Contains all constraints and data to solve the global optimization optimization problem. Add linear constraints to a `GlobalModel` through `JuMP.@constraint(gm.model, ...)`, and nonlinear constraints using `add_nonlinear_constraint`.
"""
@with_kw mutable struct GlobalModel
    model::JuMP.Model                                            # Associated JuMP.Model
    name::String = "Model"                                       # Name
    bbls::Array{BlackBoxLearner} = BlackBoxLearner[]             # Constraints to be learned
    vars::Array{JuMP.VariableRef} = JuMP.all_variables(model)    # JuMP variables
    objective = JuMP.objective_function(model)                # Original objective function
    solution_history::DataFrame = DataFrame(string.(vars) .=> [Float64[] for i=1:length(vars)]) # Solution history
    cost::Array = []                                             # List of costs.
    soldict::Dict = Dict()                                       # For solution extraction
    params::Dict = gm_defaults()                                 # GM settings
end

function Base.show(io::IO, gm::GlobalModel)
    println(io, "GlobalModel " * gm.name * " with $(length(gm.vars)) variables: ")
    println(io, "Has $(length(gm.bbls)) BlackBoxLearners.")
    if get_param(gm, :ignore_feasibility)
        if get_param(gm, :ignore_accuracy)
            println(io, "Ignores training accuracy and data_feasibility thresholds.")
        else
            println(io, "Ignores data feasibility thresholds.")
        end
    else
        if get_param(gm, :ignore_accuracy)
            println("Ignores training accuracy thresholds.")
        end
    end
end

function set_param(gm::GlobalModel, key::Symbol, val)
    set_param(gm.params, key, val)
    for bbl in gm.bbls
        set_param(bbl.params, key, val, false)
    end
end

get_param(gm::GlobalModel, key::Symbol) = get_param(gm.params, key)

"""
    (gm::GlobalModel)(name::String)

Finds BlackBoxLearner in GlobalModel by name.
"""
function (gm::GlobalModel)(name::String)
    fn_names = getfield.(gm.bbls, :name)
    fns = gm.bbls[findall(x -> x == name, fn_names)]
    if length(fns) == 1
        return fns[1]
    elseif length(fns) == 0
        @warn("No constraint with name ", name)
        return
    else
        @warn("Multiple constraints with name ", name)
        return fns
    end
end

""" Finds active leaves of all constraints in GlobalModel."""
function active_leaves(gm::GlobalModel)
    for bbl in gm.bbls
        active_leaves(bbl)
        for ll in bbl.lls
            active_leaves(ll)
        end
    end
end

"""
    $(TYPEDSIGNATURES)

Extends JuMP.all_variables to GlobalModels and BlackBoxLearners, but only returns non-auxiliary variables
"""
JuMP.all_variables(bbo::Union{GlobalModel, BlackBoxLearner, LinkedLearner}) = bbo.vars

function JuMP.all_variables(bbls::Union{Array{BlackBoxLearner}, Array{LinkedLearner}})
    return unique(Iterators.flatten(([JuMP.all_variables(bbl) for bbl in bbls])))
end

""" Extends JuMP.set_optimizer to GlobalModels. """
function JuMP.set_optimizer(gm::GlobalModel, optimizer_factory)
    JuMP.set_optimizer(gm.model, optimizer_factory)
end

count_types(gm::GlobalModel) = count_types(gm.model)
count_variables(gm::GlobalModel) = count_variables(gm.model)
count_constraints(gm::GlobalModel) = count_constraints(gm.model)

"""
    $(TYPEDSIGNATURES)

Returns bounds of all variables.
"""
function get_bounds(model::Union{GlobalModel, JuMP.Model, BlackBoxLearner,
                                 Array{BlackBoxLearner}, LinkedLearner,
                                 Array{LinkedLearner}})
    return get_bounds(JuMP.all_variables(model))
end

""" Returns bounds of all variables, flattened in order of .vars attribute. """
function flattened_bounds(bbc::Union{GlobalModel, BlackBoxLearner, LinkedLearner})
    var_bounds = get_bounds(bbc)
    return [var_bounds[var] for var in bbc.vars]
end

"""
    $(TYPEDSIGNATURES)

Returns only unbounded variables.
"""
function get_unbounds(model::Union{GlobalModel, JuMP.Model, BlackBoxLearner,
                      Array{BlackBoxLearner}})
    return get_unbounds(JuMP.all_variables(model))
end

"""
    $(TYPEDSIGNATURES)

Takes on parsing and allocation of variables depending on user input.
"""
function determine_vars(gm::GlobalModel,
                        constraint::Union{JuMP.ConstraintRef, Expr};
                        vars::Union{Nothing, Array{JuMP.VariableRef, 1}} = nothing,
                        expr_vars::Union{Nothing, Array} = nothing)
    if constraint isa JuMP.ConstraintRef
        if isnothing(vars) && isnothing(expr_vars)
            vars = vars_from_constraint(constraint)
            return vars, vars
        elseif !isnothing(vars)
            return vars, vars
        end
    else
        if isnothing(vars) && isnothing(expr_vars)
            expr_vars = vars_from_expr(constraint, gm.model)
            vars = flat(expr_vars)
            return vars, expr_vars
        elseif isnothing(expr_vars)
            expr_vars = vars_from_expr(constraint, gm.model)
            @assert length(flat(expr_vars)) >= length(vars)
            return vars, expr_vars
        elseif isnothing(vars)
            vars = flat(expr_vars)
            return vars, expr_vars
        else
            @assert length(flat(expr_vars)) >= length(vars)
            return vars, expr_vars
        end
    end
end

"""
    add_nonlinear_constraint(gm::GlobalModel,
                     constraint::Union{JuMP.ConstraintRef, Expr};
                     vars::Union{Nothing, Array{JuMP.VariableRef, 1}} = nothing,
                     expr_vars::Union{Nothing, Array} = nothing,
                     dependent_var::Union{Nothing, JuMP.VariableRef} = nothing,
                     name::String = gm.name * " " * string(length(gm.bbls) + 1),
                     equality::Bool = false)

 Adds a new nonlinear constraint to Global Model. Standard method for adding BlackBoxClassifiers and BlackBoxRegressors.

 Note: If adding constraints via an `Expr`, please define the appropriate vector indices. Eg.
    `expr = :((x) -> -x[1:\$(N)]'*\$(Q)*x[1:\$(N)] - \$(c)'*x[1:\$(N)])`,
instead of
    `expr = :((x) -> -x'*\$(Q)*x - \$(c)'*x)`,
to avoid vectorization errors.
"""
function add_nonlinear_constraint(gm::GlobalModel,
                     constraint::Union{JuMP.ConstraintRef, Expr};
                     vars::Union{Nothing, Array{JuMP.VariableRef, 1}} = nothing,
                     expr_vars::Union{Nothing, Array} = nothing,
                     dependent_var::Union{Nothing, JuMP.VariableRef} = nothing,
                     name::String = "bbl" * string(length(gm.bbls) + 1),
                     equality::Bool = false)
    vars, expr_vars = determine_vars(gm, constraint, vars = vars, expr_vars = expr_vars)
    if constraint isa Expr
        if isnothing(dependent_var)
            new_bbl = BlackBoxClassifier(constraint = constraint, vars = vars, expr_vars = expr_vars,
                                         equality = equality, name = name)
            set_param(new_bbl, :n_samples, Int(ceil(get_param(gm, :sample_coeff)*sqrt(length(vars)))))
            push!(gm.bbls, new_bbl)
            return
        else
            new_bbl = BlackBoxRegressor(constraint = constraint, vars = vars, expr_vars = expr_vars,
                                        dependent_var = dependent_var, equality = equality, name = name)
            set_param(new_bbl, :n_samples, Int(ceil(get_param(gm, :sample_coeff)*sqrt(length(vars)))))
            push!(gm.bbls, new_bbl)
            return
        end
    elseif constraint isa JuMP.ConstraintRef
        !isnothing(dependent_var) && throw(OCTHaGOnException("Constraint " * name * " is of type $(string(typeof(constraint))) " *
                                                        "and cannot have a dependent variable " * string(dependent_var) * "."))
        new_bbl = BlackBoxClassifier(constraint = constraint, vars = vars, expr_vars = expr_vars,
                                        equality = equality, name = name)
        set_param(new_bbl, :n_samples, Int(ceil(get_param(gm, :sample_coeff)*sqrt(length(vars)))))
        push!(gm.bbls, new_bbl)
        JuMP.delete(gm.model, constraint)
    end
end

"""
    function add_nonlinear_or_compatible(gm::GlobalModel,
        constraint::Union{JuMP.ConstraintRef, Expr};
        vars::Union{Nothing, Array{JuMP.VariableRef, 1}} = nothing,
        expr_vars::Union{Nothing, Array} = nothing,
        dependent_var::Union{Nothing, JuMP.VariableRef} = nothing,
        name::String = gm.name * "_" * string(length(gm.bbls) + 1),
        equality::Bool = false)

Extends `add_nonlinear_constraint` to recognize JuMP compatible linear constraints and add them as linear JuMP constraints instead.
"""
function add_nonlinear_or_compatible(gm::GlobalModel,
                     constraint::Union{JuMP.ConstraintRef, Expr};
                     vars::Union{Nothing, Array{JuMP.VariableRef, 1}} = nothing,
                     expr_vars::Union{Nothing, Array} = nothing,
                     dependent_var::Union{Nothing, JuMP.VariableRef} = nothing,
                     name::String = gm.name * "_" * string(length(gm.bbls) + 1),
                     equality::Bool = false)
    vars, expr_vars = determine_vars(gm, constraint, vars = vars, expr_vars = expr_vars)
    fn = functionify(constraint)
    if fn isa Function
        try
            constr_expr = Base.invokelatest(fn, expr_vars...)
            if constr_expr isa JuMP.GenericAffExpr || get_param(gm, :convex_constrs)
                if equality && !isnothing(dependent_var)
                    @constraint(gm.model, dependent_var == constr_expr)
                elseif equality
                    @constraint(gm.model, constr_expr == 0)
                elseif !isnothing(dependent_var)
                    @constraint(gm.model, dependent_var >= constr_expr)
                else
                    @constraint(gm.model, constr_expr >= 0)
                end
            else
                throw(ErrorException())
            end
        catch
            add_nonlinear_constraint(gm, constraint, vars = vars, expr_vars = expr_vars,
                                     dependent_var = dependent_var,
                                     name = name, equality = equality)
        end
    else
        add_nonlinear_constraint(gm, constraint, vars = vars, expr_vars = expr_vars,
                                    dependent_var = dependent_var,
                                    name = name, equality = equality)
    end
end

"""
    add_linked_constraint(gm::GlobalModel, bbc::BlackBoxClassifier, linked_vars::Array{JuMP.Variable})
    add_linked_constraint(gm::GlobalModel, bbr::BlackBoxRegressor, linked_vars::Array{JuMP.Variable}, linked_dependent::JuMP.Variable)

Adds a linked constraint of the same structure as the BlackBoxLearner.
When a nonlinear constraint is repeated more than once, this function allows the underlying
approximator to be replicated without retraining trees for each constraint.
Note that the bounds used for sampling are for the original variables of the BlackBoxLearner, so be careful!
"""
function add_linked_constraint(gm::GlobalModel, bbc::BlackBoxClassifier, vars::Array{JuMP.VariableRef})
    length(vars) == length(bbc.vars) || throw(OCTHaGOnException("BBC $(bbc.name) does not" *
    " have the same number of variables as linked variables $(vars)."))
    if !isempty(bbc.mi_constraints)
        clear_tree_constraints!(gm, bbc)
        @info "Cleared constraints from BBC $(bbc.name) since it was relinked."
    end
    lc = LinkedClassifier(vars = vars, equality = bbc.equality)
    push!(bbc.lls, lc)
    return
end

function add_linked_constraint(gm::GlobalModel, bbr::BlackBoxRegressor, vars::Array{JuMP.VariableRef},
                                dependent_var::JuMP.VariableRef)
    length(vars) == length(bbr.vars) || throw(OCTHaGOnException("BBR $(bbr.name) does not" *
    " have the same number of variables as linked variables $(vars)"))
    if !isempty(bbr.mi_constraints)
        clear_tree_constraints!(gm, bbr)
        @info "Cleared constraints from BBR $(bbr.name) since it was relinked."
    end
    lr = LinkedRegressor(vars = vars, dependent_var = dependent_var,
                                   equality = bbr.equality)
    push!(bbr.lls, lr)
    return
end

"""
    nonlinearize!(gm::GlobalModel, bbls::Array{BlackBoxLearner})
    nonlinearize!(gm::GlobalModel)

Turns gm.model into the nonlinear representation of global optimization problem.
NOTE: to get back to MI-compatible forms, must regenerate GlobalModel from scratch.
"""
function nonlinearize!(gm::GlobalModel, bbls::Array{BlackBoxLearner})
    for (i, bbl) in enumerate(bbls)
        if bbl.constraint isa JuMP.ConstraintRef
            JuMP.add_constraint(gm.model, bbl.constraint)
        elseif bbl.constraint isa Expr
            symb = Symbol(bbl.name)
            vars = flat(bbl.expr_vars) # We want flattening of dense vars.
            var_ranges = get_var_ranges(bbl.expr_vars)
            expr = bbl.constraint
            flat_expr = :((x...) -> $(expr)([x[i] for i in $(var_ranges)]...))
            fn = eval(flat_expr)
            JuMP.register(gm.model, symb, length(vars), fn; autodiff = true)
            expr = Expr(:call, symb, vars...)
            if bbl.equality && bbl isa BlackBoxRegressor
                JuMP.add_NL_constraint(gm.model, :($(expr) == $(bbl.dependent_var)))
            elseif bbl.equality
                JuMP.add_NL_constraint(gm.model, :($(expr) == 0))
            elseif bbl isa BlackBoxRegressor
                JuMP.add_NL_constraint(gm.model, :($(expr) <= $(bbl.dependent_var)))
            else
                JuMP.add_NL_constraint(gm.model, :($(expr) >= 0))
            end
        end
    end
    return
end

function nonlinearize!(gm::GlobalModel)
    nonlinearize!(gm, gm.bbls)
end

function bound!(model::GlobalModel, bounds::Union{Pair,Dict})
    bound!(model.model, bounds)
end

""" Separates and returns linear and nonlinear constraints in a model. """
function classify_constraints(model::Union{GlobalModel, JuMP.Model})
    jump_model = model
    if model isa GlobalModel
        jump_model = model.model
    end
    all_types = list_of_constraint_types(jump_model)
    nl_constrs = [];
    l_constrs = [];
    l_vartypes = [JuMP.VariableRef, JuMP.GenericAffExpr{Float64, VariableRef}]
    l_constypes = [MOI.GreaterThan{Float64}, MOI.LessThan{Float64}, MOI.EqualTo{Float64}]
    for (vartype, constype) in all_types
        constrs_of_type = JuMP.all_constraints(jump_model, vartype, constype)
        if any(vartype .== l_vartypes) && any(constype .== l_constypes)
            append!(l_constrs, constrs_of_type)
        else
            append!(nl_constrs, constrs_of_type)
        end
    end
    if !isnothing(jump_model.nlp_data)
        append!(nl_constrs, jump_model.nlp_data.nlconstr)
    end
    return l_constrs, nl_constrs
end

""" Returns the feasibility of data points in a bbl or GM. """
feasibility(bbc::BlackBoxClassifier) = bbc.feas_ratio
feasibility(bbr::BlackBoxRegressor) = size(bbr.X, 1) / (size(bbr.X, 1) + size(bbr.infeas_X, 1))
feasibility(bbls::BlackBoxLearner) = feasibility.(bbls)
feasibility(gm::GlobalModel) = feasibility.(gm.bbls)

""" Returns the accuracy of learners in a bbl or GM. """
function evaluate_accuracy(bbc::BlackBoxClassifier)
    if bbc.convex
        return 1.
    elseif bbc.feas_ratio in [1., 0]
        @warn(string("Accuracy of BlackBoxClassifier ", bbc.name, " is tautological."))
        return 1.
    elseif isempty(bbc.learners)
        throw(OCTHaGOnException(string("BlackBoxClassifier ", bbc.name, " has not been trained yet.")))
    else
        return bbc.accuracies[end]
    end
end

function evaluate_accuracy(bbr::BlackBoxRegressor)
    if bbr.convex && !bbr.equality
        return 1.
    elseif isempty(bbr.learners)
        throw(OCTHaGOnException(string("BlackBoxRegressor ", bbr.name, " has not been trained yet.")))
    else
        return bbr.accuracies[active_lower_tree(bbr)]
    end
end

evaluate_accuracy(gm::GlobalModel) = evaluate_accuracy.(gm.bbls)

"""
    JuMP.optimize!(gm::GlobalModel; kwargs...)

Applies JuMP.optimize! to GlobalModels, and saves solution history.
"""
function JuMP.optimize!(gm::GlobalModel; kwargs...)
    JuMP.optimize!(gm.model, kwargs...) # Solves the optimization problem.
    append!(gm.solution_history, solution(gm), cols=:intersect) # Saves the solution.
    feas_gap(gm) # Computes the feasibility gaps of all constraints.
    push!(gm.cost, JuMP.getobjectivevalue(gm.model)) # Updates the final cost.
    active_leaves(gm) # Updates the active leaves of all approximations.
    gm.soldict = Dict(key => JuMP.getvalue.(gm.model[key]) for (key, value) in gm.model.obj_dict)
    return
end

"""
    solution(gm::GlobalModel)
    solution(m::JuMP.Model)

Returns the optimal solution of the GlobalModel/JuMP.Model in a DataFrame.
"""
function solution(gm::GlobalModel)
    vals = getvalue.(gm.vars)
    return DataFrame(vals', string.(gm.vars))
end

function solution(m::JuMP.Model)
    variables = JuMP.all_variables(m)
    vals = getvalue.(variables)
    return DataFrame(vals', string.(variables))
end

"""
    feas_gap(gm::GlobalModel, soln = solution(gm))

Evaluates relative feasibility gap of each nonlinear constraint at the given solution.
Negative values -> constraint violation for BBCs,
                    regression underestimation for BBRs.
Positive values -> constraint violation for BBC equalities,
                    regression overestimation for BBRs.
"""
function feas_gap(gm::GlobalModel, soln = solution(gm))
    for bbl in gm.bbls
        if bbl isa BlackBoxClassifier && !isnothing(bbl.constraint)
            for ll in bbl.lls # LL feas_gaps evaluated first, for descent function
                eval!(bbl, DataFrame(string.(bbl.vars) .=> values(soln[1, string.(ll.vars)])))
                push!(ll.feas_gap, bbl.Y[end] ./ (bbl.max_Y - bbl.min_Y))
            end
            eval!(bbl, soln)
            push!(bbl.feas_gap, bbl.Y[end] ./ (bbl.max_Y - bbl.min_Y))
        elseif bbl isa BlackBoxRegressor && !isnothing(bbl.constraint)
            for ll in bbl.lls # LL feas_gaps evaluated first, for descent function
                eval!(bbl, DataFrame(string.(bbl.vars) .=> values(soln[1, string.(ll.vars)])))
                optimum = soln[:, string(ll.dependent_var)][1]
                actual = bbl.Y[end]
                push!(ll.optima, optimum)
                push!(ll.actuals, actual)
                push!(ll.feas_gap, (optimum-actual) / ((bbl.max_Y - bbl.min_Y)))
            end
            eval!(bbl, soln)
            optimum = soln[:, string(bbl.dependent_var)][1]
            actual = bbl.Y[end]
            push!(bbl.optima, optimum)
            push!(bbl.actuals, actual)
            push!(bbl.feas_gap, (optimum-actual) / ((bbl.max_Y - bbl.min_Y)))
        elseif bbl isa BlackBoxClassifier && isnothing(bbl.constraint)
            for ll in bbl.lls # LL feas_gaps evaluated first, for descent function
                push!(ll.feas_gap, 0)
            end
            push!(bbl.feas_gap, 0) # data constraints are always feasible
        elseif bbl isa BlackBoxRegressor && isnothing(bbl.constraint)
            for ll in bbl.lls # LL feas_gaps evaluated first, for descent function
                optimum = soln[:, string(ll.dependent_var)][1]
                push!(ll.optima, optimum)
                push!(ll.feas_gap, 0) # data constraints are always feasible
            end
            optimum = soln[:, string(bbl.dependent_var)][1]
            push!(bbl.optima, optimum)
            push!(bbl.feas_gap, 0) # data constraints are always feasible
        end
    end
    return
end

""" Prints the last feasibility gap of each constraint. """
function print_feas_gaps(gm::GlobalModel)
    @info "Feasibility gaps:"
    for bbl in gm.bbls
        println("$(bbl.name): $(bbl.feas_gap[end])")
        for ll in bbl.lls
            println("Linked: $(ll.feas_gap[end])")
        end
    end
    return
end

""" Shows feasibility of last solution w.r.t. each approximated constraint. """
function is_feasible(bbl::Union{BlackBoxLearner, LinkedLearner}, tighttol = 1e-8)
    if bbl.equality
        return abs(bbl.feas_gap[end]) <= tighttol
    else
        return bbl.feas_gap[end] >= -tighttol
    end
end

""" Returns the feasibility of the GlobalModel. """
function is_feasible(gm::GlobalModel)
    tighttol = get_param(gm, :tighttol)
    for bbl in gm.bbls
        is_feasible(bbl, tighttol) || return false
        for ll in bbl.lls
            is_feasible(ll, tighttol) || return false
        end
    end
    return true
end

""" Checks whether a BBL or GM is sampled. """
is_sampled(bbl::BlackBoxLearner) = size(bbl.X,1) != 0

is_sampled(gm::GlobalModel) = all(is_sampled(bbl) for bbl in gm.bbls)

""" Clears all sampling, training and optimization data from GlobalModel. """
function clear_data!(gm::GlobalModel)
    clear_tree_constraints!(gm, gm.bbls)
    clear_data!.(gm.bbls)
    gm.solution_history = DataFrame(string.(gm.vars) .=> [Float64[] for i=1:length(gm.vars)])
end

update_vexity(gm::GlobalModel) = update_vexity.(gm.bbls)

""" Prints relevant details of a GlobalModel. """
function print_details(gm::GlobalModel)
    @info "GlobalModel $(gm.name) has:"
    n_vars = length(gm.vars)
    @info "$(n_vars) variables,"
    all_types = list_of_constraint_types(gm.model)
    l_vartypes = [JuMP.VariableRef, JuMP.GenericAffExpr{Float64, VariableRef}]
    l_constrs = []
    nl_constrs = []
    l_constypes = [MOI.GreaterThan{Float64}, MOI.LessThan{Float64}, MOI.EqualTo{Float64}]
    for (vartype, constype) in all_types
        constrs_of_type = JuMP.all_constraints(gm.model, vartype, constype)
        if any(vartype .== l_vartypes) && any(constype .== l_constypes)
            append!(l_constrs, constrs_of_type)
        else
            append!(nl_constrs, constrs_of_type)
        end
    end
    @info "$(length(l_constrs)) linear constraints."
    @info "$(length(nl_constrs)) nonlinear constraints."

    n_eqs = length(findall(x -> x.equality, gm.bbls))
    n_ineqs = length(gm.bbls) - n_eqs
    obj_bbl = filter(x -> x.dependent_var == gm.objective,
        [bbl for bbl in gm.bbls if bbl isa BlackBoxRegressor])
    @info "$(n_ineqs - length(obj_bbl)) black box inequalities,"
    @info "$(n_eqs) black box equalities,"
    n_lls = sum(length(bbl.lls) for bbl in gm.bbls)
    if n_lls > 0
        @info "$(n_lls) linked constraints,"
    end
    if isempty(obj_bbl)
        @info "And a linear objective."
    else
        @info "And a nonlinear objective."
    end
    return
end

""" Sets whether a BlackBoxLearner should be represented with big-M or big-M free constraints. """
function set_bigM(bbl::BlackBoxLearner, b::Bool)
    bbl.bigM = b
    return
end

set_bigM(gm::GlobalModel, b::Bool) = set_bigM.(gm.bbls, b)

"""
    $(TYPEDSIGNATURES)

Detects BlackBoxLearner with similar structure, for speeding up constraint learning.
"""
function detect_linked_constraints!(gm::GlobalModel)
    bounds = get_bounds(gm.bbls)
    vars = all_variables(gm.bbls)
    isnothing(check_bounds(bounds)) || throw(OCTHaGOnException("Not all variables in nonlinear constraints are bounded in $(gm.name)."))
    bins = Dict()
    queue = copy(gm.bbls)
    while !isempty(queue)
        bbl = popfirst!(queue)
        n_vars = length(bbl.vars)
        reg_bool = bbl isa BlackBoxRegressor
        similars = [q for q in queue if (n_vars == length(q.vars)
                        && reg_bool == (q isa BlackBoxRegressor) && bbl.equality == q.equality)]
        if !isempty(similars)
            samples = lh_sample(bbl, lh_iterations = 5, n_samples = 10)
            vals = evaluate(bbl, samples)
            act_sims = []
            for i = 1:length(similars)
                sim = similars[i]
                rename!(samples, string.(sim.vars))
                qs = evaluate(sim, samples)
                if all(isapprox.(vals, qs, atol = 1e-8))
                    push!(act_sims, sim)
                end
            end
            if !isempty(act_sims)
                queue = filter!(e -> !(e in act_sims), queue)
                bbl_bounds = flattened_bounds(bbl)
                ll_bounds = [flattened_bounds(sim) for sim in act_sims]
                outer_bounds = []
                for i = 1:length(bbl_bounds)
                    var_bound = [
                        minimum([minimum(bbl_bounds[i]),
                                [minimum(ll_bound[i]) for ll_bound in ll_bounds]...])
                        maximum([maximum(bbl_bounds[i]),
                                [maximum(ll_bound[i]) for ll_bound in ll_bounds]...])]
                    push!(outer_bounds, var_bound)
                end
                bins[bbl] = act_sims
            end
        else
            continue
        end
    end
    return bins
end

"""
    $(TYPEDSIGNATURES)

Merges detected linked constraints via LinkedLearners.
"""
function merge_linked_constraints!(gm::GlobalModel, bins::Dict)
    @info "Merging $(length(flat(values(bins)))) LinkedLearners into $(length(keys(bins))) BlackBoxLearners. "
    for (bbl, linked_bbls) in bins
        for linked_bbl in linked_bbls
            if linked_bbl isa BlackBoxClassifier
                add_linked_constraint(gm, bbl, linked_bbl.vars)
            else
                add_linked_constraint(gm, bbl, linked_bbl.vars, linked_bbl.dependent_var)
            end
        end
    end
    gm.bbls = filter!(x -> !(x in flat(values(bins))), gm.bbls)
    return
end


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

"""
    localanisotropies(data, angles, ranges, convention=:GSLIB)

Import a `LocalAnisotropy` object from an outer source. The `data` is a
georeferenced object. It must have the rotation angles and ranges as
properties. These property names are informed as vectors in `angles` and
`ranges`for the conversion. Check out the available rotation conventions at
[`RotationRule`](@ref) docstring.

## Example

```julia
localanisotropies(data, [:rot1], [:range1, :range2], convention=:EulerIntr) # 2D
localanisotropies(data, [:rot1, :rot2, :rot3], [:range1, :range2, :range3]) # 3D
```
"""
function localanisotropies(data::SpatialData, angles::AbstractVector,
                           ranges::AbstractVector, convention=:GSLIB)
    ranges = Symbol.(ranges)
    angles = Symbol.(angles)
    tab  = Tables.columns(values(data))
    cols = string.(Tables.columnnames(tab))
    len  = nelements(data)
    @assert string.(angles) ⊆ cols "angle column name do not exist"
    @assert string.(ranges) ⊆ cols "range column name do not exist"
    dim = length(ranges)
    rules[convention].main == :y && (ranges = ranges[reverse(1:dim,1,2)])

    q = Array{Quaternion}(undef,len)
    m = Array{Vector}(undef,len)

    for i in 1:len
        xranges = [Tables.getcolumn(tab,x)[i] for x in ranges] ./ Tables.getcolumn(tab,ranges[1])[i]
        xangles = [Tables.getcolumn(tab,x)[i] for x in angles]

        P, Λ = rotmat(xranges, xangles, convention; rev=false)
        size(P,1) == 2 && (P=DCM([P[1,1] P[1,2] 0; P[2,1] P[2,2] 0; 0 0 1]))
        q[i] = dcm_to_quat(P)
        m[i] = xranges
    end

    LocalAnisotropy(q, reduce(hcat,m))
end

#
"""
    convertangles(angles, convention1, convention2)

Convert angles between different convention. Check out the available rotation
conventions at [`RotationRule`](@ref) docstring.

## Example

Converting 3D rotation [30, 30 30] from GSLIB convention to Datamine convention

```julia
new_angs = convertangles([30,30,30], :GSLIB, :Datamine)
```
"""
function convertangles(angles::AbstractVector, convention1, convention2)
    P, _  = rotmat([1,1,1], angles, convention1)
    rotmat2angles(P, convention2)
end

"""
    convertangles(quaternion, convention)

Convert a `ReferenceFrameRotations.Quaternion` to angles in a given convention.
Check out the available rotation conventions at [`RotationRule`](@ref) docstring.
"""
function convertangles(quat::Quaternion, convention)
    dcm    = quat_to_dcm(quat)
    rotmat2angles(dcm, convention)
end

"""
    exportpars(filename, localaniso, convention)

Export local anisotropies to a CSV table. The angles are exported in the form of
the given convention. Check out the available rotation conventions at
[`RotationRule`](@ref) docstring.

## Example

```julia
exportpars("path/localaniso.csv", localaniso, :GSLIB)
```
"""
function exportpars(out, lpars::LocalAnisotropy, convention=:GSLIB)
    table = convertpars(lpars, convention)
    CSV.write(out, table)
end

# convert LocalAnisotropy into angles + ranges
function convertpars(lpars::LocalAnisotropy, convention=:GSLIB)
    pars = []
    len  = nvals(lpars)
    for i in 1:len
        dcm    = rotmat(lpars, i)
        angles = rotmat2angles(dcm, convention)
        ranges = magnitude(lpars, i)
        push!(pars, append!(angles,ranges))
    end
    pars = reduce(hcat, pars)
    r3   = size(pars,1) == 6 ? pars[6,:] : [0 for i in 1:len]
    ((ang1=pars[1,x],ang2=pars[2,x],ang3=pars[3,x],r1=pars[4,x],r2=pars[5,x],r3=r3[x]) for x in 1:len)
end

# reverse transformation of rotation matrix to angles
function rotmat2angles(dcm::AbstractMatrix, convention)
  N = size(dcm, 1)
  P = N == 2 ? DCM([P[1,1] P[1,2] 0; P[2,1] P[2,2] 0; 0 0 1]) : dcm

  rule = convention isa Symbol ? rules[convention] : convention
  !rule.extrinsic && (P = P')
  preangs = dcm_to_angle(P, rule.order)
  angles  = [preangs.a1, preangs.a2, preangs.a3]

  intr = @. (rule.motion == :CW)  & !rule.extrinsic
  extr = @. (rule.motion == :CCW) & rule.extrinsic
  angles[intr .| extr] *= -1

  !rule.radian && (angles = rad2deg.(angles))
  angles
end


text:
module TestArguments

end # module


text:
const ∅ = BK_EMPTY

function PATH(s, t, i, source, sink, PARENT, INDEX)
    path = Int[]
    idxs = Int[]
    sᵢ = tᵢ = i
    while s != source
        pushfirst!(path, s)
        pushfirst!(idxs, sᵢ)
        s, sᵢ = PARENT[s], INDEX[s]
    end
    pushfirst!(path, s)
    pushfirst!(idxs, sᵢ)
    pop!(idxs)  # pop i out since we gonna push it again
    while t != sink
        push!(path, t)
        push!(idxs, tᵢ)
        t, tᵢ = PARENT[t], INDEX[t]
    end
    push!(path, t)
    push!(idxs, tᵢ)
    return path, idxs
end

function boykov_kolmogorov(source::Int, sink::Int, neighbors::Vector{Vector{Tuple{Int,Int}}}, weights::AbstractMatrix)
    # initialize: S = {s}, T = {t},  A = {s,t}, O = ∅
    vertexNum = length(neighbors)
    PARENT = zeros(Int, vertexNum)
    INDEX = zeros(Int, vertexNum)
    STATUS = fill(BK_FREE, vertexNum)
    STATUS[source] = BK_S_ACTIVE
    STATUS[sink] = BK_T_ACTIVE
    ORPHAN = fill(false, vertexNum)
    A = Int[source, sink]  # this queue could also contains inactive nodes which will be automatically skipped in the growth stage
    O = Int[]
    residual = copy(weights)
    flow = zero(Float64)
    # Boykov, Yuri, and Vladimir Kolmogorov. "An experimental comparison of min-cut/max-flow
    # algorithms for energy minimization in vision." Pattern Analysis and Machine Intelligence,
    # IEEE Transactions on 26.9 (2004): 1124-1137.
    while true
        # “growth” stage: search trees S and T grow until they touch giving an s → t path
        # grow S or T to find an augmenting path P from s to t
        P, IDX = growth_stage!(source, sink, neighbors, residual, A, STATUS, PARENT, INDEX)
        isempty(P) && break
        # “augmentation” stage: the found path is augmented, search tree(s) break into forest(s)
        flow += augmentation_stage!(neighbors, residual, P, IDX, O, STATUS, PARENT, INDEX)
        # “adoption” stage: trees S and T are restored
        adoption_stage!(source, sink, neighbors, residual, O, A, ORPHAN, STATUS, PARENT, INDEX)
    end
    return flow, STATUS
end

boykov_kolmogorov(source, sink, neighbors, weights::AbstractVector) = boykov_kolmogorov(source, sink, neighbors, reshape(weights, 2, :))

function growth_stage!(source, sink, neighbors, residual, A, STATUS, PARENT, INDEX)
    TREE(x) = STATUS[x] & (BK_S | BK_T)
    while !isempty(A)
        p = last(A)  # pick an active node p ∈ A ("First-In-First-Out"): enqueue -> queue -> dequeue
        STATUS[p] & BK_ACTIVE == BK_ACTIVE || (pop!(A); continue)  # automatically skip inactive node
        @inbounds for (q,qᵢ) in neighbors[p]
            tree_cap = TREE(p)==BK_S ? (p < q ? residual[1,qᵢ] : residual[2,qᵢ]) : (p < q ? residual[2,qᵢ] : residual[1,qᵢ])
            tree_cap > 0 || continue
            if TREE(q) == ∅
                # then add q to search tree as an active node
                STATUS[q] = STATUS[p]
                PARENT[q] = p
                pushfirst!(A, q)  # "First-In-First-Out": enqueue -> queue -> dequeue
                INDEX[q] = qᵢ   # cache index for referencing residual in augmentation stage
            end
            TREE(q)≠∅ && TREE(q)≠TREE(p) && return TREE(p)==BK_S ? PATH(p,q,qᵢ,source,sink,PARENT,INDEX) : PATH(q,p,qᵢ,source,sink,PARENT,INDEX)
        end
        # remove p from A ("First-In-First-Out"): enqueue -> queue -> dequeue
        STATUS[p] &= ~BK_ACTIVE
        pop!(A)
    end
    return Int[], Int[]
end

function augmentation_stage!(neighbors, residual, P, IDX, O, STATUS, PARENT, INDEX)
    TREE(x) = STATUS[x] & (BK_S | BK_T)
    # find the bottleneck capacity Δ on P
    Δ = Inf
    @inbounds for i = 1:length(P)-1
        p, q = P[i], P[i+1]
        if p < q
            Δ > residual[1,IDX[i]] && (Δ = residual[1,IDX[i]];)
        else
            Δ > residual[2,IDX[i]] && (Δ = residual[2,IDX[i]];)
        end
    end
    # update the residual graph by pushing flow Δ through P
    @inbounds for i = 1:length(P)-1
        p, q = P[i], P[i+1]
        if p < q
            residual[1,IDX[i]] -= Δ
            residual[2,IDX[i]] += Δ
            residual[1,IDX[i]] == 0 || continue
        else
            residual[2,IDX[i]] -= Δ
            residual[1,IDX[i]] += Δ
            residual[2,IDX[i]] == 0 || continue
        end
        if TREE(p) == TREE(q) == BK_S
            PARENT[q] = 0
            pushfirst!(O, q)
            INDEX[q] = 0  # clean up index cache
        end
        if TREE(p) == TREE(q) == BK_T
            PARENT[p] = 0
            pushfirst!(O, p)
            INDEX[p] = 0  # clean up index cache
        end
    end
    return Δ
end

function adoption_stage!(source, sink, neighbors, residual, O, A, ORPHAN, STATUS, PARENT, INDEX)
    TREE(x) = STATUS[x] & (BK_S | BK_T)
    fill!(ORPHAN, false)
    while !isempty(O)
        # pick an orphan node p ∈ O and remove it from O
        p = pop!(O)
        # find a new valid parent for p among its neighbors
        has_valid_parent = false
        @inbounds for (q,qᵢ) in neighbors[p]
            ORPHAN[q] && continue
            TREE(q) == TREE(p) || continue
            tree_cap = TREE(p)==BK_S ? (q < p ? residual[1,qᵢ] : residual[2,qᵢ]) : (q < p ? residual[2,qᵢ] : residual[1,qᵢ])
            tree_cap > 0 || continue
            # the “origin” of q should be either source or sink, it should not originates from orphan
            x = q
            while PARENT[x] ≠ 0
                x = PARENT[x]
                ORPHAN[x] && break
                ORPHAN[x] = true
            end
            if x == source || x == sink
                PARENT[p] = q
                INDEX[p] = qᵢ
                has_valid_parent = true
                fill!(ORPHAN, false)
                break
            end
        end
        has_valid_parent && continue
        @inbounds for (q,qᵢ) in neighbors[p]
            TREE(q) == TREE(p) || continue
            tree_cap = TREE(p)==BK_S ? (q < p ? residual[1,qᵢ] : residual[2,qᵢ]) : (q < p ? residual[2,qᵢ] : residual[1,qᵢ])
            tree_cap > 0   && (STATUS[q] |= BK_ACTIVE; pushfirst!(A, q);)
            PARENT[q] == p && (PARENT[q] = 0; INDEX[q] = 0; pushfirst!(O, q);)
        end
        # TREE(p) := ∅, A := A - {p}
        STATUS[p] = BK_FREE   # note that this also marks p as inactive node
    end
end

max_stars_repo_name: Gnimuc/BKmaxflow

max_issues_repo_name: Gnimuc/BKmaxflow

max_forks_repo_name: Gnimuc/BKmaxflow


text:
module DocumentFormat
using CSTParser
import CSTParser.Tokenize.Tokens

struct Edit{T}
    loc::T
    text::String
end

mutable struct State{T}
    offset::Int
    edits::T
end

function format(text)
    state = State(0, Edit[])
    x = CSTParser.parse(text, true)
    pass(x, state, operator_pass)
    state.offset = 0
    pass(x, state, tuple_pass)
    state.offset = 0
    pass(x, state, curly_pass)
    state.offset = 0
    pass(x, state, call_pass)
    sort!(state.edits, lt = (a,b) -> first(a.loc) < first(b.loc), rev = true)
    for i = 1:length(state.edits)
        text = apply(text, state.edits[i])
    end
    text = indents(text)
    return text

end

function pass(x::CSTParser.LeafNode, state, f = (x,edits)->nothing)
    state.offset += x.fullspan
end

function pass(x, state, f = (x,state)->nothing)
    f(x, state)
    for a in x
        pass(a, state, f)
    end
    state
end

function ensure_single_space_after(x, state, offset)
    if x.fullspan == last(x.span)
        if x isa CSTParser.OPERATOR
            if length(x.span) > 1 && length(String(Expr(x))) == 1
                push!(state.edits, Edit(offset + 1, " "))
            else
                push!(state.edits, Edit(offset + x.fullspan, " "))
            end
        else
            push!(state.edits, Edit(offset + x.fullspan, " "))
        end
    end
end

function ensure_no_space_after(x, state, offset)
    if x.fullspan != last(x.span)
        push!(state.edits, Edit(offset .+(last(x.span)+1:x.fullspan), ""))
    end
end

function apply(text, edit::Edit{Int})
    string(text[1:edit.loc], edit.text, text[nextind(text, edit.loc):end])
end
function apply(text, edit::Edit{UnitRange{Int}})
    string(text[1:prevind(text, first(edit.loc))], edit.text, text[nextind(text, last(edit.loc)):end])
end
include("passes.jl")
include("indents.jl")
end


text:
using SparseArrays
using LinearAlgebra

include("utils.jl")
include("GloVe.jl")
include("Bias.jl")

target = get(ARGS, 1, "C0-V15-W8-D75-R0.05-E300")
pert_dir = get(ARGS, 2, "results/perturbations")
embedding_dir = get(ARGS, 3, "embeddings")

function pre_compute(M, X, word_indices)
    num_words = length(word_indices)
    H = Dict{Int64,Array{Float64,2}}()
    G = Dict{Int64,Array{Float64,1}}()
    for i in word_indices
        H[i] = inv(GloVe.∇²Li(M.U, X, i))
        G[i] = GloVe.∇Li(M.W, M.b_w, M.U, M.b_u, X, i)
    end
    return H, G
end

function predict_bias(pert_path, M, X, inv_hessians, gradients, weat_idx_set)

    Y = GloVe.load_cooc(pert_path, M.V)
    # Make the IF approximation
    target_indices = unique([i for inds in weat_idx_set for i in inds])
    deltas = GloVe.compute_IF_deltas(Y, M, X, target_indices, inv_hessians, gradients)
    # Compute the bias change
    return Bias.effect_size(M.W, weat_idx_set, deltas)
end

function main()
    vocab_path = abspath(joinpath(embedding_dir, "vocab-$(split(target, "-W")[1]).txt"))
    println("Vocab: $vocab_path")
    vocab, ivocab = GloVe.load_vocab(vocab_path)

    weat_idx_sets = [Bias.get_weat_idx_set(set, vocab) for set in Bias.WEAT_WORD_SETS]
    all_weat_indices = unique([i for set in weat_idx_sets for inds in set for i in inds])

    cooc_path = abspath(joinpath(embedding_dir, "cooc-$(split(target, "-D")[1]).bin"))
    println("Cooc: $cooc_path")
    X = GloVe.load_cooc(cooc_path, length(vocab), all_weat_indices)

    n = 0
    for model_filename in readdir(embedding_dir)
        if occursin(Regex("$target-S[0-9]+.bin"), model_filename)
            n += 1
            println("Model: $model_filename")
            M = GloVe.load_model(joinpath(embedding_dir, model_filename))
            @assert M.vocab_path == vocab_path
            H, G = pre_compute(M, X, all_weat_indices)
            effect_sizes = [Bias.effect_size(M.W, set) for set in weat_idx_sets]

            for i in 1:length(Bias.WEAT_WORD_SETS)
                target_dir = joinpath(pert_dir, target * "-B$i")
                if (n == 1)
                    open(joinpath(target_dir, "predicted_change.csv"), "w") do out_io
                        headers = ["model_file","pert_file", "pert_type", "pert_size", "pert_run", "seed", "B", "IFB̃"]
                        println(out_io, join(headers,","))
                    end
                end
                open(joinpath(target_dir, "predicted_change.csv"), "a") do out_io
                    for pert_filename in readdir(target_dir)
                        if startswith(pert_filename, "pert-")
                            println("Pert: $pert_filename")
                            fields = split(pert_filename[6:end-4], "_")
                            pert_type = fields[1]
                            pert_size = parse(Int, fields[2])
                            pert_run = parse(Int, fields[3])
                            seed = parse(Int, match(r"-S[0-9]+", model_filename).match[3:end])
                            pert_path = joinpath(target_dir, pert_filename)
                            IFB̃ = predict_bias(pert_path, M, X, H, G, weat_idx_sets[i])
                            data = [model_filename, pert_filename, pert_type, pert_size, pert_run, seed, effect_sizes[i], IFB̃]
                            println(out_io, join(data, ","))
                        end
                    end
                end
            end
        end
    end

end

main()


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

"""
    BisectPointPartitioner(normal, point)

A method for partitioning spatial data into two half spaces
defined by a `normal` direction and a reference `point`.
"""
struct BisectPointPartitioner{T,N} <: AbstractPartitioner
  normal::SVector{N,T}
  point::SVector{N,T}
end

BisectPointPartitioner(normal::NTuple{N,T},
                  point::NTuple{N,T}=ntuple(i->zero(T), N)) where {T,N} =
  BisectPointPartitioner{T,N}(normalize(SVector(normal)), SVector(point))

function partition(object::AbstractSpatialObject{T,N},
                   partitioner::BisectPointPartitioner{T,N}) where {T,N}
  n = partitioner.normal
  p = partitioner.point

  x = MVector{N,T}(undef)

  left  = Vector{Int}()
  right = Vector{Int}()
  for location in 1:npoints(object)
    coordinates!(x, object, location)
    if (x - p) ⋅ n < zero(T)
      push!(left, location)
    else
      push!(right, location)
    end
  end

  SpatialPartition(object, [left,right])
end


text:
######
# This file is part of the BivariateCopulas.jl package.
#
#   Plotting joint cdfs, densities and samples
#
#           University of Liverpool, Institute for Risk and Uncertainty
#
#                                           Author: Ander Gray
#                                           Email:  ander.gray@liverpool.ac.uk
######

wait_for_key() = (print(stdout, "press enter to continue"); read(stdin, 1); nothing);

function samplePlot(C :: AbstractCopula, N = 1)     # Generates a plot showing conditional sampling

    x = rand(N);    y = rand(N);
    ux = x;         uy = zeros(N);
    m = n;
    e = 1/m;        ygrid = range(0,stop=1,length = m);

    if !ismissing(C.func) e = sqrt(eps()); end

    for i=1:N
        conditional = (C(x[i] + e/2,ygrid) - C(x[i] - e/2,ygrid))./e
        uy[i] = findlast(conditional .<= y[i]) / m;
        println("-------------------------------------------")
        println("N: $i.")
        println("Independent ~U(0,1)'s:")
        println(hcat(x[i],y[i]))
        println();

        plot2(C.cdf,conditional,ux[i], uy[i], y[i])

        println("Samples from copula:");
        println(hcat(ux[i],uy[i]))
        println("-------------------------------------------")

        wait_for_key()
    end
    return hcat(ux,uy);
end

plotDensity(x :: AbstractCopula) = plot1(x.density);
plotCdf(x::AbstractCopula) = plot1(x.cdf);

function plotCdf(j::AbstractJoint; pn = 50)

    A = j.cdf
    m = size(A)[1];
    if m < pn; ppn = m; else ppn = pn; end

    nm = round(m/ppn);
    A = A[1:Int(nm):end,1:Int(nm):end]
    plot3(A, j.xRange[1:Int(nm):end],j.yRange[1:Int(nm):end])

end

function plotDensity(j::AbstractJoint)

    A = j.density
    m = size(A)[1];
    if m < pn; ppn = m; else ppn = pn; end

    nm = round(m/ppn);
    A = A[1:Int(nm):end,1:Int(nm):end]
    plot3(A, j.xRange[1:Int(nm):end],j.yRange[1:Int(nm):end])

end

plot(j::AbstractJoint) = plot4(j)
plotDen(j::AbstractJoint) = plot4(j,false)

plotContourCdf(j::AbstractJoint) = plot5(j)
plotContourDen(j::AbstractJoint) = plot5(j,false)

function plot(x :: AbstractCopula;title = "SurfacePlots", pn = 50, fontsize=18)
    A = x.cdf
    m = size(A)[1];
    if m < pn; ppn = m; else ppn = pn; end

    x = y = range(0,stop = 1,length=ppn)
    xgrid = repeat(x',ppn,1)
    ygrid = repeat(y,1,ppn)

    nm = round(m/ppn);

    z = A[1:Int(nm):end,1:Int(nm):end]

    fig = figure(title,figsize=(10,10))
    ax = fig.add_subplot(1,1,1,projection="3d")
    #ax = fig.add_subplot(2,1,1)
    plot_surface(xgrid, ygrid, z, rstride=2,edgecolors="k", cstride=2, alpha=1, linewidth=2, cmap=ColorMap("RdGy"))
    ax.view_init(45-27, 180+ 26)
    xlabel("y",fontsize=fontsize)
    ylabel("x",fontsize=fontsize)
    ax.zaxis.set_rotate_label(false);
    zlabel("C(x,y)", rotation = 0,fontsize=fontsize)
    #xticks(fontsize = fontsize); yticks(fontsize = fontsize);
    #PyPlot.title(title,fontsize=fontsize)
    tight_layout()

end

function plotDen(x :: AbstractCopula; title = "SurfacePlots", pn = 50, fontsize=18)
    A = x.density
    m = size(A)[1];
    if m < pn; ppn = m; else ppn = pn; end

    x = y = range(0,stop = 1,length=ppn)
    xgrid = repeat(x',ppn,1)
    ygrid = repeat(y,1,ppn)

    nm = round(m/ppn);

    z = A[1:Int(nm):end,1:Int(nm):end]

    fig = figure(title,figsize=(10,10))
    ax = fig.add_subplot(1,1,1,projection="3d")
    #ax = fig.add_subplot(2,1,1)
    plot_surface(xgrid, ygrid, z, rstride=2,edgecolors="k", cstride=2, alpha=1, linewidth=2, cmap=ColorMap("RdGy"))
    ax.view_init(45-27, 180+ 26)
    xlabel("y",fontsize=fontsize)
    ylabel("x",fontsize=fontsize)
    ax.zaxis.set_rotate_label(false);
    zlabel("C(x,y)", rotation = 0,fontsize=fontsize)
    #xticks(fontsize = fontsize); yticks(fontsize = fontsize);
    #PyPlot.title(title,fontsize=fontsize)
    tight_layout()

end

function plotContourDen(x; title = "SurfacePlots", pn = pn, fontsize=18)
    A = x.density
    m = size(A)[1];
    if m < pn; ppn = m; else ppn = pn; end

    x = y = range(0,stop = 1,length=ppn)
    xgrid = repeat(x',ppn,1)
    ygrid = repeat(y,1,ppn)

    nm = round(m/ppn);

    z = A[1:Int(nm):end,1:Int(nm):end]
    fig = figure(title,figsize=(10,10))

    cp = contour(xgrid, ygrid, z,cmap=ColorMap("coolwarm"),levels = 20)
    xlabel("y",fontsize=fontsize)
    ylabel("x",fontsize=fontsize)

end

function plotContourCdf(x :: AbstractCopula; title = "SurfacePlots",fontsize=18)
    A = x.cdf
    m = size(A)[1];
    if m < pn; ppn = m; else ppn = pn; end

    x = y = range(0,stop = 1,length=ppn)
    xgrid = repeat(x',ppn,1)
    ygrid = repeat(y,1,ppn)

    nm = round(m/ppn);

    z = A[1:Int(nm):end,1:Int(nm):end]
    fig = figure(title,figsize=(10,10))

    cp = contour(xgrid, ygrid, z,cmap=ColorMap("coolwarm"),levels = 20)
    xlabel("y",fontsize=fontsize)
    ylabel("x",fontsize=fontsize)

end

function plot1(A :: Array{<:Real,2})

#   https://github.com/gizmaa/Julia_Examples/blob/master/pyplot_surfaceplot.jl
    m = size(A)[1];
    if m < pn; ppn = m; else ppn = pn; end

    x = y = range(0,stop=1,length=ppn)
    xgrid = repeat(x',ppn,1)
    ygrid = repeat(y,1,ppn)

    nm = round(m/ppn);

    z = A[1:Int(nm):end,1:Int(nm):end]

    fig = figure("SurfacPlots",figsize=(10,10))
    ax = fig.add_subplot(2,1,1,projection="3d")
    plot_surface(xgrid, ygrid, z, rstride=2,edgecolors="k", cstride=2, alpha=0.8, linewidth=0.25,cmap=ColorMap("coolwarm"))
    xlabel("X")
    ylabel("Y")
    ax.set_zlim3d(0, max(maximum(A),1))
    PyPlot.title("Surface Plot")

    subplot(212)
    ax = fig.add_subplot(2,1,2)
    cp = contour(xgrid, ygrid, z,cmap=ColorMap("coolwarm"),levels = 15)
    ax.clabel(cp, inline=1, fontsize=10)
    xlabel("X")
    ylabel("Y")
    #ax.set_zlim3d(0, max(maximum(A),1))
    PyPlot.title("Contour Plot")
    tight_layout()

end

function plot2(A :: Array{<:Real,2}, cond:: Array{<:Real,1}, ux = missing, uy = missing, uuy = missing)

    if !isempty(get_fignums()) PyPlot.clf(); end

    m = size(A)[1];
    if m < pn; ppn = m; else ppn = pn; end

    x = y = range(0,stop=1,length=ppn)
    xgrid = repeat(x',ppn,1)
    ygrid = repeat(y,1,ppn)

    nm = m/ppn;

    z = A[1:Int(nm):end, 1:Int(nm):end]
    conditional = cond[1:Int(nm):end];

    fig = figure("ConditionalPlots",figsize=(10,10))
    ax = fig.add_subplot(2,1,1)
    cp = contour(xgrid, ygrid, z)
    ax.clabel(cp, inline=1, fontsize=10)
    xlabel("X")
    ylabel("Y")

    if !ismissing(ux)
        xx = ones(pn) * ux
        PyPlot.plot(xx,x, color = "red",linewidth=2,label = "Rand1: $ux");
    end

    subplot(212)
    ax = fig.add_subplot(2,1,2)
    PyPlot.plot(y, conditional, color = "black")
    if !ismissing(uy) && !ismissing(uuy)
        nums1 = Int(round(pn*uy))
        nums2 = Int(round(pn*uuy))
        yy = ones(nums1) * uuy
        uyy = ones(nums2) * uy
        PyPlot.plot(y[1:nums1],yy,color = "red",linewidth=2, label= "Rand2: $uy");
        PyPlot.plot(uyy,y[1:nums2],color = "blue",linewidth=2);
    end
    xlabel("Y")
    ylabel("conditional CDF")
    tight_layout()
end

function plot3(A :: Array{<:Real,2}, x = range(0,1,length = pn), y = range(0,1,length = pn))    # Won't work if n!=200

#   https://github.com/gizmaa/Julia_Examples/blob/master/pyplot_surfaceplot.jl

    #m1 = size(A)[1];
    #m2 = size(A)[2];
    #if m1 < pn; ppn1 = m1; else ppn1 = pn; end      # ppn1 is the number of points to plot in X
    #if m2 < pn; ppn2 = m2; else ppn2 = pn; end      # ppn2 is the number of points to plot in Y

    #nm = round(m/ppn);                              # The descritization of the domain

    #x = range(rangeX[1],rangeX[2],length=pn)
    #y = range(rangeY[1],rangeY[1],length=pn)
    xgrid = repeat(x',pn,1)
    ygrid = repeat(y,1,pn)
    z = A'

    fig = figure("SurfacePlots",figsize=(10,10))
    ax = fig.add_subplot(2,1,1,projection="3d")
    plot_surface(xgrid, ygrid, z, rstride=2,edgecolors="k", cstride=2, cmap=ColorMap("coolwarm"), alpha=0.8, linewidth=0.25)
    #plot_surface(xgrid, ygrid, z, rstride=2,edgecolors="k", cstride=2, cmap=ColorMap("Gray"), alpha=0.8, linewidth=0.25)
    xlabel("X")
    ylabel("Y")
    PyPlot.title("Surface Plot")

    #subplot(212)
    ax = fig.add_subplot(2,1,2)
    cp = contour(xgrid, ygrid, z, cmap=ColorMap("coolwarm"))
    ax.clabel(cp, inline=1, fontsize=10)
    xlabel("X")
    ylabel("Y")
    PyPlot.title("Contour Plot")
    tight_layout()

end

function plot4(J :: AbstractJoint, CDF = true)    # Won't work if n!=200

#   https://github.com/gizmaa/Julia_Examples/blob/master/pyplot_surfaceplot.jl

    x = J.xRange;       y = J.yRange;
    M1 = J.marginal1;   M2 = J.marginal2;

    if CDF; z = J.cdf; title = "cdf" else z = J.density; title = "pdf";end

    m = size(z)[1];
    if m < pn; ppn = m; else ppn = pn; end

    nm = round(m/ppn);
    z = z[1:Int(nm):end,1:Int(nm):end]

    x = x[1:Int(nm):end]
    y = y[1:Int(nm):end]

    xgrid = repeat(x',pn,1)
    ygrid = repeat(y,1,pn)

    fig1 = figure("SurfacePlots",figsize=(10,10))
    #ax = fig.add_subplot(2,1,1,projection="3d")
    plot_surface(xgrid, ygrid, z, rstride=2,edgecolors="k", cstride=2, cmap=ColorMap("coolwarm"), alpha=0.8, linewidth=0.25)
    xlabel("X")
    ylabel("Y")

    PyPlot.title("$title surface plot")

    tight_layout()

end

function plot5(J :: AbstractJoint, CDF = true)    # Won't work if n!=200

    #   https://github.com/gizmaa/Julia_Examples/blob/master/pyplot_surfaceplot.jl

        x = J.xRange;       y = J.yRange;
        M1 = J.marginal1;   M2 = J.marginal2;

        if CDF; z = J.cdf; title = "cdf" else z = J.density; title = "pdf";end

        m = size(z)[1];
        if m < pn; ppn = m; else ppn = pn; end

        nm = round(m/ppn);
        z = z[1:Int(nm):end,1:Int(nm):end]

        x = x[1:Int(nm):end]
        y = y[1:Int(nm):end]

        xgrid = repeat(x',pn,1)
        ygrid = repeat(y,1,pn)

        fig = plt.figure("ContourPlots", figsize=(10, 10))
        grid = plt.GridSpec(4, 4, hspace=0.2, wspace=0.2)
        main_ax = fig.add_subplot(get(grid, (slice(1,4),slice(0,3))))
        yMarg  = fig.add_subplot(get(grid, (slice(1,4),3)), xticklabels=[], sharey=main_ax)
        xMarg  = fig.add_subplot(get(grid, (0,slice(0,3))), yticklabels=[], sharex=main_ax)

        cp = main_ax.contour(xgrid, ygrid, z, cmap=ColorMap("coolwarm"), levels =15)
        main_ax.clabel(cp, inline=1, fontsize=10)
        #xlabel("X")
        #ylabel("Y")
        if (CDF) xMarg.plot(x,cdf.(M1,x));else xMarg.plot(x,pdf.(M1,x));end
        if (CDF) yMarg.plot(cdf.(M2,y),y);else yMarg.plot(pdf.(M2,y),y);end

        PyPlot.title("joint $title and marginals")
        tight_layout()

    end

function plotDist(D :: AbstractMarginal; name = "marginal", col = "red", fontsize = 12)

    fig = figure(name,figsize=(10,10))

    PyPlot.plot(D.range, D.cdf, color = col)

    xticks(fontsize = fontsize); yticks(fontsize = fontsize)
    xlabel("Distribution range",fontsize = fontsize); ylabel("CDF",fontsize=fontsize);

end

slice(x,y) = pycall(pybuiltin("slice"), PyObject, x, y)         # For performing python array slicing for

function scatter(a :: Array{Float64,2}; title = "samples",fontsize = 12, fontsizeT = 12, xlab = "U", ylab = "V", nbins = 50)

    x = a[:,1];
    y = a[:,2];

    fig = plt.figure(title, figsize=(10, 10))
    grid = plt.GridSpec(4, 4, hspace=0.2, wspace=0.2)
    main_ax = fig.add_subplot(get(grid, (slice(1,4),slice(0,3))))

    # scatter points on the main axes
    main_ax.plot(x, y, "o", markersize=3, alpha=0.2)

    xticks(fontsize=fontsizeT); yticks(fontsize=fontsizeT)
    xlabel(xlab, fontsize =  fontsize); ylabel(ylab, fontsize = fontsize);

    x_hist  = fig.add_subplot(get(grid, (0,slice(0,3))), xticklabels=[], yticklabels=[])#, sharex=main_ax)

    # histogram on the attached axes
    x_hist.hist(x, nbins, histtype="stepfilled", orientation="vertical", color="gray")
    #x_hist.invert_yaxis()
    xticks(fontsize=fontsizeT); yticks(fontsize=fontsizeT)

    y_hist  = fig.add_subplot(get(grid, (slice(1,4),3)), xticklabels=[], yticklabels=[])#, sharey=main_ax)

    y_hist.hist(y, nbins, histtype="stepfilled", orientation="horizontal", color="gray")
    #y_hist.invert_xaxis()
    xticks(fontsize=fontsizeT); yticks(fontsize=fontsizeT)

end


text:
import Pkg; Pkg.add(url = "https://github.com/tobega/RelationalData")
using RelationalData

rules = Relation((container = String, contained = String, amount = Int64),
  [(container = match(r"(\w+ \w+) bags contain", line).captures[1],
    contained = m.captures[2],
    amount =  parse(Int64, m.captures[1]))
    for line in eachline("a7.txt")
      for m in eachmatch(r"(\d+) (\w+ \w+)", line)
        ])

#part 1
function containers(b)
  isempty(b) && return b
  # c = rename(project(rules ⨝ b, :container), :container => :contained)
  c = rules ⨝ b |> r -> project(r, :container) |> r -> rename(r, :container => :contained)
  return c ∪ containers(c)
end

length(containers(Relation((contained = "shiny gold",)))) |> println

#alt part 1
function containers2(b)
  isempty(b) && return b
  c = rules ⨝ b |> r -> map(t -> (contained = t.container,), r)
  return c ∪ containers2(c)
end

length(containers2(Relation((contained = "shiny gold",)))) |> println

#part 2
function contents(b)
  isempty(b) && return b
  c = (project(b, :amount, :contained) |> r -> rename(r, :amount => :multiplier, :contained => :container)
    ⨝ rules) |>
    r -> extend(r, :tally, t -> t.amount * t.multiplier) |>
    r -> project(r, :container, :contained, :tally) |>
    r -> rename(r, :tally => :amount)
  return c ∪ contents(c)
end

sum(t -> t.amount, contents(Relation((container="x", contained="shiny gold", amount=1)))) |> println

#alternatively
function contents2(b)
  isempty(b) && return b
  c = (map(t -> (multiplier = t.amount, container = t.contained), b) ⨝ rules) |>
    r -> map(t -> (container = t.container, contained = t.contained, amount = t.amount * t.multiplier), r)
  return c ∪ contents2(c)
end

sum(t -> t.amount, contents2(Relation((container="x", contained="shiny gold", amount=1)))) |> println


text:
# generateur d'Instances pour la RODD

n = 10
P = 6

# Ouvrir le fichier "output.txt" dans lequel on pourra écrire
fout = open("outputn10P6.txt", "w")
println(fout, "title")

println(fout, "n = " * string(n))
println(fout, "P = " * string(P))

for i in 1:n
  for j in 1:n
    println(fout, string(i) * " " * string(j) * " " * string(rand(1:10)))
  end
end

for k in 1:P
  for i in 1:n
    for j in 1:n
      if rand(1:10) <= 9
        println(fout, string(k) * " " * string(i) * " " * string(j) * " " * string((rand(Float64,1)/2)[1]))
      end
    end
  end
end

close(fout)

max_stars_repo_name: boumpteryx/RODD

max_issues_repo_name: boumpteryx/RODD

max_forks_repo_name: boumpteryx/RODD


text:
export SeisAddNoise,
SeisLinearEvents,
SeisParabEvents,
SeisHypEvents,
AcousticSetup,
SeisAcousticWave

include("SeisAddNoise.jl")
include("SeisLinearEvents.jl")
include("SeisParabEvents.jl")
include("SeisHypEvents.jl")
include("SeisAcousticWave.jl")

max_stars_repo_name: fercarozzi/myseismicjulia

max_issues_repo_name: fercarozzi/myseismicjulia

max_forks_repo_name: fercarozzi/myseismicjulia


text:
n = parse(Int64, readline())

dict = Dict{String, Int64}()

for i in 1:n
    A = split(readline(), " ")
    dict[A[1]] = parse(Int64, A[2])
end

for i in 1:n
    key = readline()
    val = get(dict,key,"-1")
    if val == "-1"
        println("Not found")
    else
        println(key,"=",val)
    end
end

max_stars_repo_path: Hackerrank/30 Days of Code/Julia/day 08.jl

max_stars_repo_name: Next-Gen-UI/Code-Dynamics

max_issues_repo_name: Next-Gen-UI/Code-Dynamics

max_forks_repo_name: Next-Gen-UI/Code-Dynamics


text:
#!/bin/bash
# -*- mode: julia -*-
#=
exec julia --project=@. --compile=min -O0 --startup-file=no -e "include(popfirst!(ARGS))" "${BASH_SOURCE[0]}" "$@"
=#

using Raw2CSV

if isempty(ARGS)
	println("usage: ./$PROGRAM_FILE <list of all files>")
	exit()
end

csv = gather_csv_from_files(ARGS, [
	("instancePath", key_equals_extractor("instance_path", NoDefault{String}())),
	("instanceName", wrap(basename,
		key_equals_extractor("instance_path", NoDefault{String}()))),
	("rotationOfItems", wrap(Int,
		p_args_key_extractor("PPG2KP-allow-rotation", NoDefault{Bool}()))),
	("exportMPS", wrap(isempty,
		p_args_key_extractor("save-model", NoDefault{String}()))),
	("numberOfVariables", key_equals_extractor("num_vars", 0)),
	("numberOfConstraints", key_equals_extractor("num_constrs", 0)),
	("preprocessing_time", key_equals_extractor("enumeration_time", NaN)),
	("model_building_time", key_equals_extractor("build_time", NaN)),
	("is_faithful",
		p_args_key_extractor("PPG2KP-faithful2furini2016", NoDefault{Bool}())),
	("use_mirror_plates",
		p_args_key_extractor("PPG2KP-mirror-plates", false)),
	key_equals("qt_pevars_after_preprocess", -1),
	key_equals("qt_cmvars_after_preprocess", -1),
	key_equals("qt_plates_after_preprocess", -1),
	key_equals("qt_pevars_after_purge", -1),
	key_equals("qt_cmvars_after_purge", -1),
	key_equals("qt_plates_after_purge", -1),
	key_equals("build_stop_reason", "REASON_NOT_FOUND"),
	key_equals("save_model_time", NaN),
	key_equals("run_total_time", NaN),
	("exception_detected",
		key_equals_extractor("run_ended_by_exception", false)),
])

print(csv)


text:
#this file contains: read_bin, convert2array, convert2gcmfaces

## convert2array method:

"""
    convert2array(fld::MeshArray)

Convert MeshArray to Array (or vice versa otherwise)
"""
function convert2array(fld::MeshArray)

if fld.grid.class=="llc";
    tmp1=cat(fld.f[1],fld.f[2],rotr90(fld.f[4]),rotr90(fld.f[5]);dims=1);
    tmp2=cat(rotl90(fld.f[3]),NaN*fld.f[3],NaN*fld.f[3],NaN*fld.f[3];dims=1);
    arr=cat(tmp1,tmp2;dims=2);
elseif fld.grid.class=="cs";
    tmp1=cat(fld.f[1],fld.f[2],rotr90(fld.f[4]),rotr90(fld.f[5]);dims=1);
    tmp2=cat(rotl90(fld.f[3]),NaN*fld.f[3],NaN*fld.f[3],NaN*fld.f[3];dims=1);
    tmp0=cat(NaN*fld.f[3],NaN*fld.f[3],NaN*fld.f[3],rotr90(fld.f[6]);dims=1);
    arr=cat(tmp0,tmp1,tmp2;dims=2);
elseif fld.grid.class=="ll";
  arr=fld.f[1];
else;
  error("unknown grTopo case");
end;

return arr;

end

## convert2array ##

function convert2array(fld::Array{T,N},grid::gcmgrid) where {T,N}

grTopo=grid.class
nFaces=grid.nFaces
facesSize=grid.fSize

v1=Array{Array{T,N},1}(undef,nFaces);
N>2 ? error("N>2 case not implemented yet") : nothing

if grTopo=="llc";
    (n1,n2)=facesSize[1];
    v1[1]=fld[1:n1,1:n2];
    v1[2]=fld[n1+1:n1*2,1:n2];
    v1[3]=rotr90(fld[1:n1,n2+1:n2+n1]);
    v1[4]=rotl90(fld[n1*2+1:n1*3,1:n2]);
    v1[5]=rotl90(fld[n1*3+1:n1*4,1:n2]);
elseif grTopo=="cs";
    (n1,n2)=facesSize[1];
    v1[1]=fld[1:n1,n1+1:n1+n2];
    v1[2]=fld[n1+1:2*n1,n1+1:n1+n2];
    v1[3]=rotr90(fld[1:n1,n1+n2+1:n2+n1*2]);
    v1[4]=rotl90(fld[n1*2+1:n1*3,n1+1:n1+n2]);
    v1[5]=rotl90(fld[n1*3+1:n1*4,n1+1:n1+n2]);
    v1[6]=rotl90(fld[n1*3+1:n1*4,1:n1]);
elseif grTopo=="ll";
    v1[1]=fld;
else;
  error("unknown grTopo case");
end;

MeshArray(grid,v1);

end

## convert2gcmfaces method:

"""
    convert2gcmfaces(fld::MeshArray)

Convert mitgcm output to MeshArray (or vice versa otherwise)
"""
function convert2gcmfaces(fld::MeshArray)

    grTopo=fld.grid.class
    nFaces=fld.grid.nFaces
    (n1,n2)=fld.grid.ioSize
    facesSize=fld.grid.fSize

    aa=0;bb=0;
    for iFace=1:nFaces;
        aa=aa+prod(facesSize[iFace]);
        bb=bb+prod(size(fld.f[iFace]));
    end;
    n3=Int64(bb/aa);
    v11=NaN*zeros(aa,n3);

    i0=0; i1=0;
    for iFace=1:nFaces;
        i0=i1+1;
        nn=facesSize[iFace][1];
        mm=facesSize[iFace][2];
        i1=i1+nn*mm;
        v11[i0:i1,:]=reshape(fld.f[iFace],(nn*mm,n3));
    end;
    v1=reshape(v11,(n1,n2,n3));
    n3==1 ? v1=dropdims(v1,dims=3) : nothing

    return v1;

end

## convert2gcmfaces ##

function convert2gcmfaces(fld::Array,grid::gcmgrid)

grTopo=grid.class
nFaces=grid.nFaces
(n1,n2)=grid.ioSize
facesSize=grid.fSize

n3=Int64(prod(size(fld))/n1/n2);

v0=reshape(fld,(n1*n2,n3));
v1=Array{Array{eltype(fld),ndims(fld)}}(undef,nFaces);
i0=0; i1=0;
for iFace=1:nFaces
  i0=i1+1;
  nn=facesSize[iFace][1]; mm=facesSize[iFace][2];
  i1=i1+nn*mm;
  if n3>1;
    v1[iFace]=reshape(v0[i0:i1,:],(nn,mm,n3));
  else;
    v1[iFace]=reshape(v0[i0:i1,:],(nn,mm));
  end;
end

MeshArray(grid,v1);

end

## read_bin function with full list of argument

"""
    read_bin(fil::String,kt::Union{Int,Missing},kk::Union{Int,Missing},prec::DataType,mygrid::gcmgrid)

Read model output from binary file and convert to MeshArray. Other methods:

```
read_bin(fil::String,prec::DataType,mygrid::gcmgrid)
read_bin(fil::String,mygrid::gcmgrid)
```
"""
function read_bin(fil::String,kt::Union{Int,Missing},kk::Union{Int,Missing},prec::DataType,mygrid::gcmgrid)

  if ~ismissing(kt);
    error("non-empty kt option not implemented yet");
  end;

  if ~ismissing(kk);
    error("non-empty kk option not implemented yet");
  end;

  (n1,n2)=mygrid.ioSize

  if prec==Float64;
    reclen=8;
  else;
    reclen=4;
  end;
  tmp1=stat(fil);
  n3=Int64(tmp1.size/n1/n2/reclen);

  fid = open(fil);
  fld = Array{prec,1}(undef,(n1*n2*n3));
  read!(fid,fld);
  fld = hton.(fld);
  close(fid)

  n3>1 ? s=(n1,n2,n3) : s=(n1,n2)
  v0=reshape(fld,s);

  convert2gcmfaces(v0,mygrid)

end

## read_bin with reduced list of argument

# read_bin(fil::String,prec::DataType,mygrid::gcmgrid)
function read_bin(fil::String,prec::DataType,mygrid::gcmgrid)
  read_bin(fil,missing,missing,prec,mygrid::gcmgrid)
end

# read_bin(fil::String,mygrid::gcmgrid)
function read_bin(fil::String,mygrid::gcmgrid)
  read_bin(fil,missing,missing,mygrid.ioPrec,mygrid::gcmgrid)
end

## read_bin with alternative arguments

# read_bin(fil::String,x::MeshArray)
function read_bin(fil::String,x::MeshArray)
  read_bin(fil,missing,missing,eltype(x),x.grid::gcmgrid)
end

# read_bin(tmp::Array,mygrid::gcmgrid)
function read_bin(tmp::Array,mygrid::gcmgrid)
  convert2gcmfaces(tmp,mygrid)
end

# read_bin(tmp::Array,x::MeshArray)
function read_bin(tmp::Array,x::MeshArray)
  convert2gcmfaces(tmp,x.grid)
end


text:
using AbstractAlgebra
using Logging

global_logger(Logging.ConsoleLogger(Logging.Info))

struct RewriteSystem
    tails  # vector of tails
    rules  # vector of rewriters
end

struct Rewriter
    idx::Integer
    monom
    terminal
end

function addrule!(RS, i, m)
    push!(RS.rules, Rewriter(i, m, isone(m)))
end

function addtail!(RS, t)
    push!(RS.tails, (leading_monomial(t), tail(t)))
    addrule!(RS, length(RS.tails), one(t))
end

function findrewriter(monom, RS)
    for rule in reverse(RS.rules)
        t, d = divides(monom, rule.monom*first(RS.tails[rule.idx]))
        if t
            @info "found rule" rule
            addrule!(RS, rule.idx, rule.monom*d)
            return d, rule.monom*last(RS.tails[rule.idx])
        end
    end
    return monom, zero(monom)
end

function rewrite(mult, f, RS)
    newmonoms = []
    @warn "rewriting" mult f
    for m in monomials(f)
        mm = mult*m
        mult1, rewriter = findrewriter(mm, RS)
        @info m mm mult1 rewriter
        if !iszero(rewriter)
            append!(newmonoms, .- rewrite(mult1, rewriter, RS))
        else
            push!(newmonoms, mm)
        end
    end
    newmonoms
end

function rewrite(f, RS)
    sum(rewrite(one(f), f, RS))
end

"""

"""
function interreduce(gens::Vector{Poly}) where {Poly}
    sort!(gens, by=total_degree ∘ leading_monomial)

    RS = RewriteSystem([], [])

    # rewriters[1] = x1 + x2 + x3

    @warn "interreducing" gens

    newgens = similar(gens)
    newgens[1] = map_coefficients(c -> c // leading_coefficient(gens[1]), gens[1])

    addtail!(RS, newgens[1])

    i = 2

    while i <= length(gens)
        newgens[i] = rewrite(gens[i], RS)
        newgens[i] = map_coefficients(c -> c // leading_coefficient(newgens[i]), newgens[i])
        addtail!(RS, newgens[i])
        @error "End of $i" RS
        i += 1
    end

    newgens
end


text:
# Magic smoke starts leaking from the submarine's arithmetic logic unit (ALU).
# Without the ability to perform basic arithmetic and logic functions, the
# submarine can't produce cool patterns with its Christmas lights!
#
# It also can't navigate. Or run the oxygen system.
#
# Don't worry, though - you probably have enough oxygen left to give you enough
# time to build a new ALU.
#
# The ALU is a four-dimensional processing unit: it has integer variables w, x,
# y, and z. These variables all start with the value 0. The ALU also supports
# six instructions:
#
# - inp a - Read an input value and write it to variable a.
# - add a b - Add the value of a to the value of b, then store the result in
#   variable a.
# - mul a b - Multiply the value of a by the value of b, then store the result
#   in variable a.
# - div a b - Divide the value of a by the value of b, truncate the result to
#   an integer, then store the result in variable a. (Here, "truncate" means to
#   round the value toward zero.)
# - mod a b - Divide the value of a by the value of b, then store the remainder
#   in variable a. (This is also called the modulo operation.)
# - eql a b - If the value of a and b are equal, then store the value 1 in
#   variable a. Otherwise, store the value 0 in variable a.
#
# In all of these instructions, a and b are placeholders; a will always be the
# variable where the result of the operation is stored (one of w, x, y, or z),
# while b can be either a variable or a number. Numbers can be positive or
# negative, but will always be integers.
#
# The ALU has no jump instructions; in an ALU program, every instruction is run
# exactly once in order from top to bottom. The program halts after the last
# instruction has finished executing.
#
# (Program authors should be especially cautious; attempting to execute div
# with b=0 or attempting to execute mod with a<0 or b<=0 will cause the program
# to crash and might even damage the ALU. These operations are never intended
# in any serious ALU program.)
#
# For example, here is an ALU program which takes an input number, negates it,
# and stores it in x:
#
# inp x
# mul x -1
#
# Here is an ALU program which takes two input numbers, then sets z to 1 if the
# second input number is three times larger than the first input number, or
# sets z to 0 otherwise:
#
# inp z
# inp x
# mul z 3
# eql z x
#
# Here is an ALU program which takes a non-negative integer as input, converts
# it into binary, and stores the lowest (1's) bit in z, the second-lowest (2's)
# bit in y, the third-lowest (4's) bit in x, and the fourth-lowest (8's) bit in
# w:
#
# inp w
# add z w
# mod z 2
# div w 2
# add y w
# mod y 2
# div w 2
# add x w
# mod x 2
# div w 2
# mod w 2
#
# Once you have built a replacement ALU, you can install it in the submarine,
# which will immediately resume what it was doing when the ALU failed:
# validating the submarine's model number. To do this, the ALU will run the
# MOdel Number Automatic Detector program (MONAD, your puzzle input).
#
# Submarine model numbers are always fourteen-digit numbers consisting only of
# digits 1 through 9. The digit 0 cannot appear in a model number.
#
# When MONAD checks a hypothetical fourteen-digit model number, it uses
# fourteen separate inp instructions, each expecting a single digit of the
# model number in order of most to least significant. (So, to check the model
# number 13579246899999, you would give 1 to the first inp instruction, 3 to
# the second inp instruction, 5 to the third inp instruction, and so on.) This
# means that when operating MONAD, each input instruction should only ever be
# given an integer value of at least 1 and at most 9.
#
# Then, after MONAD has finished running all of its instructions, it will
# indicate that the model number was valid by leaving a 0 in variable z.
# However, if the model number was invalid, it will leave some other non-zero
# value in z.
#
# MONAD imposes additional, mysterious restrictions on model numbers, and
# legend says the last copy of the MONAD documentation was eaten by a tanuki.
# You'll need to figure out what MONAD does some other way.
#
# To enable as many submarine features as possible, find the largest valid
# fourteen-digit model number that contains no 0 digits. What is the largest
# model number accepted by MONAD?

# read input
program = map(eachline("input.txt")) do l
  (opcode, operands...) = split(l, ' ')
  op1 = operands[1][1]
  op2 = nothing
  if length(operands) > 1
    op2 = operands[2]
    if '0' <= op2[1] <= '9' || op2[1] == '-'
      op2 = parse(Int, op2)
    else
      op2 = op2[1]
    end
  end
  (opcode, op1, op2)
end

# location of all inp instructions
inpip = accumulate((i,_) -> findnext(((o,_,_),) -> o == "inp", program, i+1), 1:14)

# You could brute force today's puzzle by writing code to run the program with
# all possible inputs, but it will take forever to run. Instead, today's puzzle
# is mostly about figuring out what the program does. If you inspect it, you'll
# see it has 14 blocks of code that are basically identical:
#
# inp w
# mul x 0
# add x z
# mod x 26
# div z AAA
# add x BBB
# eql x w
# eql x 0
# mul y 0
# add y 25
# mul y x
# add y 1
# mul z y
# mul y 0
# add y w
# add y CCC
# mul y x
# add z y
#
# If you work that out, it becomes:
#
# w = input
# x = int((z % 26 + BBB) != w)               # always 0 or 1
# z = (z / AAA) * (25x + 1) + (w + CCC)x
#
# When BBB is positive, it's always >= 10 and AAA is always 1, so:
# * int(z % 26 + BBB != w) will always be 1 since z % 26 is always
#   non-negative, BBB >= 10, and w < 10;
# * z reduces to z = 26z + w + CCC
#
# You can think of this operation as "pushing" a digit (namely w + CCC) onto
# the end of z in base 26.
#
# When BBB is negative, AAA is always 26. AAA is 1 seven times, and 26 seven
# times. So if we're "pushing" a digit when AAA is 1, we must be popping it
# when AAA is 26. For that to be true, we want the z equation to equal exactly
# z / 26 which happens when x is zero. x is zero when (z % 26 + BBB) == w.
#
# So, basically, matching pushes with pops gives us constraints such that:
# * wpush + CCCpush == wpop - BBBpop
result = Vector{Char}(undef, 14)
stack = Vector{NTuple{2,Int}}()
sizehint!(stack, 7)
for i in 1:14
  ip = inpip[i]
  (bbb,ccc) = (program[ip+5][3], program[ip+15][3])
  if bbb > 0
    push!(stack, (i, ccc))
  else
    (j,delta) = pop!(stack)
    delta += bbb
    result[i] = min(9, 9 + delta) + '0'
    result[j] = min(9, 9 - delta) + '0'
  end
end

println("Model: ", String(result))


text:
using ObjectDetectors
using GeometryBasics

@testset "StaticObjectDetector" begin
    vertices = [Point2(2.0, 2.0), Point2(4.0, 2.0), Point2(3.0, 3.0)]
    core_args = ("class", 0.7, 256, 256, MACHINE, "detector")
    detection = PolygonDetection(vertices, core_args...)
    detector = StaticObjectDetector{PolygonDetection{Float64}}([detection])
    img = Array{UInt8}(undef, 10, 10, 3)
    detections = detect(detector, img)

    @test length(detections) == 1
    @test class_name(detections[1]) == "class"
    @test confidence(detections[1]) == 0.7f0
    @test centroid(detections[1]) == [3, 7 / 3]
end


text:
module PihVfi890

# This makes writing docstrings a bit easier.
# Docstrings are used to automatically generate documentation or help in the REPL.
using DocStringExtensions, Interpolations, Roots
using UtilityFunctions890

export Model
export init_test_model, solve, solve_one_point, euler_dev_one_point, euler_dev_one_point2
export kprime_min, kprime_max

struct Model
    wageV :: Vector{Float64}
    R :: Float64
    beta :: Float64
    u :: AbstractUtility
    # Number of capital grid points
    nk :: Int64
    # Max k grid point
    kMax :: Float64
end

## ----------  Helpers

c_min(m :: Model) = 1e-4;
betar(m :: Model) = m.beta * m.R;
n_periods(m :: Model) = length(m.wageV);
wage(m :: Model, t) = m.wageV[t];

# g + g^2 + ... + g^(T-1)
pvfactor(g, T) = (g ^ T - 1.0) / (g - 1.0);

present_value(xV, R) = sum(xV ./ (R .^ (0 : (length(xV)-1))));

euler_dev(m :: Model, ctV) = UtilityFunctions890.euler_dev(m.u, ctV, betar(m));

budget_kprime(w, R, k, c) = w .+ R .* k .- c;
budget_c(w, R, k, kPrime) = w .+ R .* k .- kPrime;

make_k_grid(m :: Model, t) = LinRange(0.0, m.kMax, m.nk);
kprime_max(m :: Model, t) = m.kMax;
# Max kPrime consistent with c > c_min
kprime_max(m :: Model, t, k) = budget_kprime(wage(m, t), m.R, k, c_min(m));
kprime_min(m :: Model, t) = 0.0;

function init_test_model()
    T = 5;
    wageV = collect(LinRange(1.0, 2.0, T));
    kMax = 20.0;
    u = UtilityCRRA(2.0);
    return Model(wageV, 1.04, 0.98, u, 20, kMax)
end

## ---------- Solve by Policy Fct Iteration

"""
	$(SIGNATURES)

Solve a model with policy function iteration.
Returns `Vector`s with interpolated policy functions k'(k) and c(k).
"""
function solve(m :: Model)
    T = n_periods(m);

    # Results are stored as interpolation objects
    kPrime_tV = Vector{Any}(undef, T);
    c_tV = Vector{Any}(undef, T);

    # Last period: easy
    kPrime_tV[T], c_tV[T] = solve_last_period(m);

    # Backward induction
    for t = (T-1) : -1 : 1
        kGridV, kPrime_tV[t], c_tV[t] = solve_one_period(m, t, kPrime_tV[t+1]);
    end

    isValid = check_solution(m, kPrime_tV, c_tV);
    if !isValid
        @warn "Solution fails checks"
    end

    return kPrime_tV, c_tV
end

function solve_last_period(m :: Model)
    T = n_periods(m);
    kGridV = make_k_grid(m, T);
    kPrimeFct = interpolate_kprime(kGridV, zeros(length(kGridV)));
    cV = budget_c(wage(m, T), m.R, kGridV, 0.0);
    cFct = interpolate_kprime(kGridV, cV);
    return kPrimeFct, cFct
end

function solve_one_period(m :: Model, t, kPrimeTomorrow)
    kGridV = make_k_grid(m, t);
    kPrimeV = zeros(length(kGridV));
    for (ik, k) in enumerate(kGridV)
        kPrimeV[ik] = solve_one_point(m, t, k, kPrimeTomorrow);
    end
    cV = budget_c(wage(m, t), m.R, kGridV, kPrimeV);
    kPrimeFct = interpolate_kprime(kGridV, kPrimeV);
    cFct = interpolate_kprime(kGridV, cV);
    return kGridV, kPrimeFct, cFct
end

function solve_one_point(m :: Model, t, k, kPrimeTomorrow)
    e_dev(kPrime) = euler_dev_one_point(m, t, k, kPrime, kPrimeTomorrow);
    kPrimeMin = 0.0;
    kPrimeMax = min(kprime_max(m, t, k), kprime_max(m, t));
    @assert kPrimeMax > kPrimeMin

    # Try corners
    eDev0 = e_dev(kPrimeMin);
    if eDev0 <= 0.0
        kPrime = kPrimeMin;
    else
        eDevMax = e_dev(kPrimeMax);
        if eDevMax >= 0.0
            kPrime = kPrimeMax
        else
            kPrime = find_zero(e_dev, (kPrimeMin, kPrimeMax, Bisection()));
        end
    end
    @assert kPrimeMin <= kPrime <= kprime_max(m, t)
    return kPrime
end

# Euler equation deviation.
# `dev > 0` means that `kPrime` implies a `c` today that is too high.
function euler_dev_one_point(m :: Model, t, k, kPrime, kPrimeTomorrow)
    cPrime = budget_c(wage(m, t+1), m.R, kPrime, kPrimeTomorrow(kPrime));
    betaRUprime = betar(m) * marg_utility(m.u, cPrime);
    dev = budget_c(wage(m, t), m.R, k, kPrime) - inv_marg_utility(m.u, betaRUprime);
    return dev
end

# Alternative Euler equation deviation.
# Easier to understand, but more non-linear.
function euler_dev_one_point2(m :: Model, t, k, kPrime, kPrimeTomorrow)
    c = budget_c(wage(m, t), m.R, k, kPrime);
    kPrimePrime = kPrimeTomorrow(kPrime);
    cPrime = budget_c(wage(m, t+1), m.R, kPrime, kPrimePrime);
    eDev = euler_dev(m, [c, cPrime]);
    return only(eDev)
end

function interpolate_kprime(kGridV :: AbstractVector{F}, kPrimeV :: AbstractVector{F}) where F
    interp = LinearInterpolation(kGridV, kPrimeV);
    return interp
end

function check_solution(m :: Model, kPrime_tV, c_tV)
    bcValid = true;
    eulerValid = true;

    T = n_periods(m);
    for t = 1 : (T-1)
        kGridV = make_k_grid(m, t);
        kp_fct = kPrime_tV[t];
        c_fct = c_tV[t];
        wage_t = wage(m, t);
        for (ik, k) in enumerate(kGridV)
            kPrime = kp_fct(k);
            c = budget_c(wage_t, m.R, k, kPrime);
            bcValid = bcValid && isapprox(c, c_fct(k), atol = 1e-4);

            if kprime_min(m, t) < kPrime < kprime_max(m, t)
                kPrimePrime = kPrime_tV[t+1](kPrime);
                cPrime = budget_c(wage(m, t+1), m.R, kPrime, kPrimePrime);
                eDev = only(euler_dev(m, [c, cPrime]));
                eulerValid = eulerValid && abs(eDev < 1e-4);
            end
        end
    end

    if !bcValid
        @warn "Budget constraint violated"
    end
    if !eulerValid
        @warn "Euler equation violated"
    end

    return bcValid  &&  eulerValid
end

end # module


text:
# This type contains measurement data (time, values and a counter)
type Measurements{T}
    i::Int64                # measurement counter
    t::Array{Float64}       # time data (when it was received by this recorder)
    t_msg::Array{Float64}   # time that the message was sent
    z::Array{T}             # measurement values
end

# This function cleans the zeros from the type above once the simulation is finished
function clean_up(m::Measurements)
    m.t     = m.t[1:m.i-1]
    m.t_msg = m.t_msg[1:m.i-1]
    m.z     = m.z[1:m.i-1,:]
end

function Quat2Euler(q::Array{Float64})
    sol = zeros(Float64,3)
    sol[1]   = atan2(2*(q[1]*q[2]+q[3]*q[4]),1-2*(q[2]^2+q[3]^2))
    sol[2]   = asin(2*(q[1]*q[3]-q[4]*q[2]))
    sol[3]   = atan2(2*(q[1]*q[4]+q[2]*q[3]),1-2*(q[3]^2+q[4]^2))
    return sol
end

function ECU_callback(msg::ECU,cmd_log::Measurements)
    cmd_log.t[cmd_log.i] = to_sec(get_rostime())
    cmd_log.t_msg[cmd_log.i] = to_sec(msg.header.stamp)
    cmd_log.z[cmd_log.i,:] = convert(Array{Float64,1},[msg.motor;msg.servo])
    cmd_log.i += 1
    nothing
end

function ECU_PWM_callback(msg::ECU,cmd_pwm_log::Measurements)
    cmd_pwm_log.t[cmd_pwm_log.i] = to_sec(get_rostime())
    cmd_pwm_log.t_msg[cmd_pwm_log.i] = to_sec(msg.header.stamp)
    cmd_pwm_log.z[cmd_pwm_log.i,:] = convert(Array{Float64,1},[msg.motor;msg.servo])
    cmd_pwm_log.i += 1
    nothing
end

function IMU_callback(msg::Imu,imu_meas::Measurements)
    imu_meas.t[imu_meas.i]      = to_sec(get_rostime())
    imu_meas.t_msg[imu_meas.i]  = to_sec(msg.header.stamp)
    imu_meas.z[imu_meas.i,:]    = [msg.angular_velocity.x;msg.angular_velocity.y;msg.angular_velocity.z;
                                    Quat2Euler([msg.orientation.w;msg.orientation.x;msg.orientation.y;msg.orientation.z]);
                                    msg.linear_acceleration.x;msg.linear_acceleration.y;msg.linear_acceleration.z]::Array{Float64}
    imu_meas.i += 1
    nothing
end

function GPS_callback(msg::hedge_pos,gps_meas::Measurements)
    gps_meas.t[gps_meas.i]      = to_sec(get_rostime())
    gps_meas.t_msg[gps_meas.i]  = msg.timestamp_ms
    gps_meas.z[gps_meas.i,:]    = [msg.x_m;msg.y_m]
    gps_meas.i += 1
    nothing
end

function pos_info_callback(msg::pos_info,pos_info_log::Measurements)
    pos_info_log.t[pos_info_log.i]      = to_sec(get_rostime())
    pos_info_log.t_msg[pos_info_log.i]  = to_sec(msg.header.stamp)
    pos_info_log.z[pos_info_log.i,:]    = [msg.s;msg.ey;msg.epsi;msg.v;msg.s_start;msg.x;msg.y;msg.v_x;msg.v_y;
                                            msg.psi;msg.psiDot;msg.x_raw;msg.y_raw;msg.psi_raw;msg.v_raw;msg.psi_drift;
                                            msg.a_x;msg.a_y;msg.a_x_raw;msg.a_y_raw]
    pos_info_log.i += 1
    nothing
end

function vel_est_callback(msg::Vel_est,vel_est_log::Measurements)
    vel_est_log.t[vel_est_log.i]      = to_sec(get_rostime())
    vel_est_log.t_msg[vel_est_log.i]  = to_sec(msg.header.stamp)
    vel_est_log.z[vel_est_log.i,:]    = [msg.vel_est,msg.vel_fl,msg.vel_fr,msg.vel_bl,msg.vel_br]
    vel_est_log.i += 1
    nothing
end

max_stars_repo_name: RoyaFiroozi/barc

max_issues_repo_name: RoyaFiroozi/barc

max_forks_repo_name: RoyaFiroozi/barc