text:
# D:\AppliRTE\repo\naza-mpc\powsybl\bin\itools convert-network --input-file 24Nodes.xiidm --output-file ampl --output-format AMPL
# D:\AppliRTE\repo\naza-mpc\powsybl\bin\itools convert-network --input-file RSO_5nodes_PJM_slackA.xiidm --output-file ampl --output-format AMPL
# D:\AppliRTE\repo\naza-mpc\powsybl\bin\itools convert-network --input-file RC_PPE2035_S0_deboucle.xiidm --output-file ampl --output-format AMPL
# D:\AppliRTE\repo\naza-mpc\powsybl\bin\itools convert-network --input-file RSO_5nodes_PJM_slackA_wind.xiidm --output-file ampl --output-format AMPL
module AmplTxt
    const GENERIC_HEADER = "ampl_network_";
    const GENERIC_EXTENSION = ".txt";

    const FILENAME_COLNAME = Dict(
        "rtc" => ["variant", "num", "tap", "table", "onLoad", "fault", "curative", "id"],
        "hvdc" => ["variant", "num", "type", "converterStation1", "converterStation2", "r (ohm)", "nomV (KV)", "convertersMode", "targetP (MW)", "maxP (MW)", "fault", "curative", "id", "description"],
        "limits" => ["variant", "num", "branch", "side", "limit (A)", "accept. duration (s)", "fault", "curative"],
        "generators" => ["variant", "num", "bus", "con. bus", "substation", "minP (MW)", "maxP (MW)", "minQmaxP (MVar)", "minQ0 (MVar)", "minQminP (MVar)", "maxQmaxP
        (MVar)", "maxQ0 (MVar)", "maxQminP (MVar)", "v regul.", "targetV (pu)", "targetP (MW)", "targetQ (MVar)", "fault", "curative", "id", "description", "P (MW)",
        "Q (MVar)"],
        "batteries" => ["variant", "num", "bus", "con. bus", "substation", "p0 (MW)", "q0 (MW)", "minP (MW)", "maxP (MW)", "minQmaxP (MVar)", "minQ0 (MVar)", "minQminP (MVar)", "maxQmaxP (MVar)", "maxQ0 (MVar)", "maxQminP (MVar)", "fault", "curative", "id", "description", "P (MW)", "Q (MVar)"],
        "buses" => ["variant", "num", "substation", "cc", "v (pu)", "theta (rad)", "p (MW)", "q (MVar)", "fault", "curative", "id"],
        "tct" => ["variant", "num", "tap", "var ratio", "x (pu)", "angle (rad)", "fault", "curative"],
        "loads" => ["variant", "num", "bus", "substation", "p0 (MW)", "q0 (MVar)", "fault", "curative", "id", "description", "p (MW)", "q (MVar)"],
        "branches" => ["variant", "num", "bus1", "bus2", "3wt num", "sub.1", "sub.2", "r (pu)", "x (pu)", "g1 (pu)", "g2 (pu)", "b1 (pu)", "b2 (pu)", "cst ratio (pu)", "ratio tc", "phase tc", "p1 (MW)", "p2 (MW)", "q1 (MVar)", "q2 (MVar)", "patl1 (A)", "patl2 (A)", "merged", "fault", "curative", "id", "description"],
        "shunts" => ["variant", "num", "bus", "con. bus", "substation", "minB (pu)", "maxB (pu)", "inter. points", "b (pu)", "fault", "curative", "id", "description", "P (MW)", "Q (MVar)", "sections count"],
        "substations" => ["variant", "num", "unused1", "unused2", "nomV (KV)", "minV (pu)", "maxV (pu)", "fault", "curative", "country", "id", "description"],
        "static_var_compensators" => ["variant", "num", "bus", "con. bus", "substation", "minB (pu)", "maxB (pu)", "v regul.", "targetV (pu)", "targetQ (MVar)", "fault", "curative", "id", "description", "P (MW)", "Q (MVar)"],
        "ptc" => ["variant", "num", "tap", "table", "fault", "curative", "id"],
        "vsc_converter_stations" => ["variant", "num", "bus", "con. bus", "substation", "minP (MW)", "maxP (MW)", "minQmaxP (MVar)", "minQ0 (MVar)", "minQminP (MVar)", "maxQmaxP (MVar)", "maxQ0 (MVar)", "maxQminP (MVar)", "v regul.", "targetV (pu)", "targetQ (MVar)", "lossFactor (%PDC)", "fault", "curative", "id", "description", "P (MW)", "Q (MVar)"],
        "lcc_converter_stations" => ["variant", "num", "bus", "con. bus", "substation", "lossFactor (%PDC)", "powerFactor", "fault", "curative", "id", "description",
        "P (MW)", "Q (MVar)"]
    );
    export AmplTxtDataRow;
    export split_with_space;

    mutable struct AmplTxtDataRow
        colNameIdx::Dict{String,Int64};

        data::Vector{Vector{String}};
    end

    function readNetworkData(name::String, root::String)
        file_path = abspath(root, GENERIC_HEADER * name * GENERIC_EXTENSION);
        # println("reading ", name, " in ", file_path);
        colNameIdx =  Dict{String,Int64}();
        idx = 1;
        colNames = FILENAME_COLNAME[name];
        for col in colNames
            push!(colNameIdx, col => idx);
            idx += 1;
        end

        data = Vector{String}[];
        open(file_path) do file
            for ln in eachline(file)
                # don't read commentted line
                if ln[1] != '#'
                    push!(data, split_with_space(ln));
                end
            end
            # println(data)
        end
        return AmplTxtDataRow(colNameIdx, data)
    end

    function readNetworkData(name::String)
        return readNetworkData(name, ".");
    end

    function split_with_space(str::String)
        result = String[];
        if length(str) > 0
            start_with_quote = startswith(str, "\"");
            buffer_quote = split(str, keepempty=false, "\"");
            i = 1;
            while i <= length(buffer_quote)
                if i > 1 || !start_with_quote
                    str2 = buffer_quote[i];
                    buffer_space = split(str2, keepempty=false);
                    for str3 in buffer_space
                        push!(result, str3);
                    end
                    i += 1;
                end
                if i <= length(buffer_quote)
                    push!(result, buffer_quote[i]);
                    i += 1;
                end
            end
        end
        return result;
    end
    function read()
        return read(".");
    end
    function read(root::String)
        amplTxt = Dict{String,AmplTxtDataRow}();
        for kvp in collect(FILENAME_COLNAME)
            name = kvp[1];
            data = readNetworkData(name, root);
            push!(amplTxt, name => data);
        end
        return amplTxt;
    end
end

max_stars_repo_name: azized/pscopf

max_issues_repo_name: azized/pscopf

max_forks_repo_name: azized/pscopf


text:
using DrawSimpleGraphs, SimpleGraphs
using Test
G =  Cycle(12)
DrawSimpleGraphs.draw(G)
@test 1==1

ext: jl


text:
abstract type AbstractDifferenceEquationAlgorithm <: DiffEqBase.DEAlgorithm end
struct DirectIteration <: AbstractDifferenceEquationAlgorithm end
struct KalmanFilter <: AbstractDifferenceEquationAlgorithm end

# The typical algorithm in discrete-time is DirectIteration()
# Unlike continuous time, there aren't many simple variations
default_alg(prob::AbstractStateSpaceProblem) = DirectIteration()

# If a normal prior, normal observational noise, no noise given, and observables provided then can use a kalman filter
default_alg(prob::LinearStateSpaceProblem{uType,uPriorType,tType,P,NP,F,AType,BType,CType,RType,ObsType,K}) where {uType,uPriorType<:MvNormal,tType,P,NP<:Nothing,F,AType,BType,CType,RType<:MvNormal,ObsType,K} = KalmanFilter()

# Select default algorithm if not provided
DiffEqBase.solve(prob::AbstractStateSpaceProblem; kwargs...) = DiffEqBase.solve(prob,
                                                                                default_alg(prob);
                                                                                kwargs...)
DiffEqBase.solve(prob::AbstractStateSpaceProblem, alg::Nothing, args...; kwargs...) = DiffEqBase.solve(prob,
                                                                                default_alg(prob),
                                                                                args...;
                                                                                kwargs...)


text:
using Test, JuLIP, LinearAlgebra
using JuLIP.Potentials
using JuLIP.Testing

@info("Testing One-Body")
E0 = rand()
V = OneBody(E0)
at = bulk(:Si, cubic=true) * (2,1,2)
println(@test energy(V, at) ≈ length(at) * E0)
println(@test norm(forces(V, at)) == 0)
D = Dict(V)
V1 = JuLIP.decode_dict(D)
println(@test V1 == V)

@info("Testing M-One-Body")
Vm = MOneBody(:Si => E0)
println(@test energy(Vm, at) ≈ energy(V, at))

E1 = rand()
V2 = MOneBody(:Si => E0, :C => E1)
at.Z[2:2:end] .= atomic_number(:C)
println(@test energy(V2, at) ≈ length(at) * (E0+E1) / 2)
println(@test norm(forces(V, at)) == 0)
D = Dict(V2)
V21 = JuLIP.decode_dict(D)
println(@test V2 == V21)


text:
using Winston

# plot some data
pl = plot(cumsum(rand(500) .- 0.5), "r", cumsum(rand(500) .- 0.5), "b")
# display the plot (not done automatically!)
display(pl)

# save the current figure
savefig("winston.svg")
# .eps, .pdf, & .png are also supported
# we used svg here because it respects the width and height specified above

max_stars_repo_name: inkydragon/JuliaByExample

max_issues_repo_name: inkydragon/JuliaByExample

max_forks_repo_name: inkydragon/JuliaByExample


text:
# --------------------------------------------------------------------------
# ACE1.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# Licensed under ASL - see ASL.md for terms and conditions.
# --------------------------------------------------------------------------

@assert length(ARGS) in [0, 1]

using ACE1, JuLIP, InteractiveUtils, ProgressMeter
using JuLIP.MLIPs: combine
using BenchmarkTools: @belapsed

function run_regression_tests!(D::Dict, label::AbstractString)
   # store meta-data: save the label
   if label in D["labels"]
      error("label $(label) already exists in the regression tests")
   end
   push!(D["labels"], label)
   # julia version and machine
   buf = IOBuffer()
   InteractiveUtils.versioninfo(buf)
   D["versioninfo"][label] = String(take!(buf))
   # read the benchmark configuration
   at = read_dict(D["at"])
   # run the actual benchmarks
   @showprogress for test in D["results"]
      basis = read_dict(test["basis"])
      tBE, tBF = run_regtest(basis, at)
      V = combine(basis, Float64.(test["coeffs"]))
      tVE, tVF = run_regtest(V, at)
      test[label] = Dict("evalbasis" => tBE, "gradbasis" => tBF,
                         "evalpot" => tVE, "gradpot" => tVF)
   end
end

function run_regtest(calc, at)
   return (@belapsed energy($calc, $at)),
          (@belapsed forces($calc, $at))
end

#---

testfile = @__DIR__() * "/regtest.json"

if !(isfile(testfile))
   error("No test file exists. Generating...")
end

# load the parameter file - keep it open for the rest of this
D = JuLIP.FIO.load_dict(testfile)

# create a label for the tests
if length(ARGS) == 1
   label = ARGS[1]
else
   # get the current git hash and use it as label
   buf = IOBuffer()
   run(pipeline(`git rev-parse --short HEAD`, stdout = buf))
   label = String(take!(buf))[1:end-1]
end

# run and save the actual tests
run_regression_tests!(D, label)
JuLIP.save_dict(testfile, D)


text:
function splitinvariants!(W::Vector{Float64},Q::Vector{Float64},A::Vector{Float64},
    Qo::Vector{Float64},Ao::Vector{Float64},beta::Vector{Float64},c0::Vector{Float64},
    rho::Float64,k::Vector{Float64},h::Float64,diff::Float64,mu::Float64,children::Vector{Int64})
    lf = Float64;
    lb = Float64;
    W1 = Vector{Float64}[];
    dW1 = Float64;
    W2 = Vector{Float64}[];
    dW2 = Float64;

    # forward wave speed for parent's distal end
    lf = (Q[1]./A[1] .+ sqrt.(0.5.*beta[1]./rho).*A[1].^0.25);
    # parent's right running invariant at current time
    W1 = ([Qo[1]./Ao[1] .+ 4.*(sqrt.(0.5.*beta[1]./rho).*Ao[1]^0.25 .- c0[1])]);
    push!(W1,Q[1]./A[1] .+ 4.*(sqrt.(0.5.*beta[1]./rho).*A[1].^0.25 .- c0[1]));
    dW1 = (W1[2] - W1[1])/k[1];
    dQ = (Q[1] .- Qo[1])./k[1];
    dA = (A[1] .- Ao[1])./k[1];
    # advance parent's distal boundary W1 by linear interp.
    Qint = Q[1] .- dQ.*lf.*h;
    Aint = A[1] .- dA.*lf.*h;
    W[1] = W1[2] .- dW1.*lf.*h .- diff.*pi.*mu./rho.*Qint./Aint.^2.*h;
    for j = 2:(length(children)+1)
        lb = (Q[j]./A[j] .- sqrt.(0.5.*beta[j]./rho).*A[j].^0.25);
        W2 = ([Q[j]./A[j] .- 4.*(sqrt.(0.5.*beta[j]./rho).*A[j].^0.25 .- c0[j])]);
        push!(W2,Qo[j]./Ao[j] .- 4.*(sqrt.(0.5.*beta[j]./rho).*Ao[j].^0.25 .- c0[j]))
        dW2 = (W2[2] .- W2[1])./k[j];
        dQ = (Qo[j] .- Q[j])./k[j];
        dA = (Ao[j] .- A[j])./k[j];
        # advance child's proximal boundary invariant
        Qint = Q[j] .- dQ.*lb.*h;
        Aint = A[j] .- dA.*lb.*h;
        W[j] = W2[1] .- dW2.*lb.*h .- diff.*pi.*mu./rho.*Qint./Aint.^2.*h;
    end
end

max_stars_repo_name: dcanuto/CardioModeling

max_issues_repo_name: dcanuto/CardioModeling

max_forks_repo_name: dcanuto/CardioModeling


text:
include("rbm.jl")

include("data.jl")

##

using Knet: norm

using Plots: plot, plot!

train(;rbm            = nothing,
       hidden_size    = 10,
       learning_rate  = 1,
       hm_epochs      = 1,
       batch_size     = length(data_train),
       do_print       = true
      ) =
begin

    rbm == nothing ?
        rbm = RBM(in_size,hidden_size) :
            hidden_size = length(rbm.hiddens)

    do_print ? (@info "Training started. \nhidden size   $(hidden_size) \nlearning rate $(learning_rate) \nbatch size    $(batch_size)") : ()

    grad_norms      = []
    grad_sums       = []

    dev_grad_norms = []
    dev_grad_sums  = []

    do_print ? (begin

        dev_grads = batch_grads(rbm, data_dev)
        @info "initial dev stats \nnorm $(norm(dev_grads)./sqrt(length(rbm.weights))) \nsum $(sum(abs.(dev_grads))./length(rbm.weights))"

    end) : ()

    for ep in 1:hm_epochs

        total_grads = nothing

        for batch in batchify(shuffle(data_train),batch_size)

            grads = batch_grads(rbm, batch)

            update_weights!(rbm, grads, learning_rate)

            total_grads == nothing ? total_grads = grads : total_grads += grads

        end

        dev_grads = batch_grads(rbm, data_dev)

        push!(grad_norms, norm(total_grads)./sqrt(length(rbm.weights)))
        push!(grad_sums, sum(abs.(total_grads))./length(rbm.weights))
        push!(dev_grad_norms, norm(dev_grads)./sqrt(length(rbm.weights)))
        push!(dev_grad_sums, sum(abs.(dev_grads))./length(rbm.weights))

        do_print ? println("Epoch $ep, train_sum $(round(grad_sums[end],digits=3)), dev_sum $(round(dev_grad_sums[end],digits=3))") : ()

    end

    min_train_norm = argmin(grad_norms)
    min_train_sum = argmin(grad_sums)
    min_dev_norm = argmin(dev_grad_norms)
    min_dev_sum = argmin(dev_grad_sums)

    p1 = plot(1:hm_epochs, grad_norms,     title="train_norm_$(hidden_size)_$(learning_rate)_$(batch_size)",xlabel="$(grad_norms[min_train_norm]) / $(min_train_norm)")
    p2 = plot(1:hm_epochs, grad_sums,      title="train_sum_$(hidden_size)_$(learning_rate)_$(batch_size)",xlabel="$(grad_sums[min_train_sum]) / $(min_train_sum)")
    p3 = plot(1:hm_epochs, dev_grad_norms, title="dev_norm_$(hidden_size)_$(learning_rate)_$(batch_size)",xlabel="$(dev_grad_norms[min_dev_norm]) / $(min_dev_norm)")
    p4 = plot(1:hm_epochs, dev_grad_sums,  title="dev_sum_$(hidden_size)_$(learning_rate)_$(batch_size)",xlabel="$(dev_grad_sums[min_dev_sum]) / $(min_dev_sum)")

    display(plot(p1,p2,p3,p4,layout=(2,2)))

rbm, [grad_norms,grad_sums,dev_grad_norms,dev_grad_sums]
end

##

using Images: Gray

generate(rbm;hidden_state=nothing,converge=true) =
begin

    if hidden_state == nothing

        random_states = randn(1,length(rbm.hiddens))

        binary ? random_states = binarize_state.(random_states) : ()

        hidden_state = random_states

    end

    rbm.hiddens = hidden_state

    rbm()

    converge ? propogate_until_convergence!(rbm, rbm.visibles) : ()

    std_transform ? (begin

        rbm.visibles = StatsBase.reconstruct(transform3, rbm.visibles)

        pca_transform ? (begin

            rbm.visibles = MultivariateStats.reconstruct(transform2, rbm.visibles')'

        end) : (begin

            rbm.visibles = StatsBase.reconstruct(transform1, rbm.visibles)

        end)

        # reconstruct!(std2, rbm.visibles)
        # reconstruct!(std1, rbm.visibles)

    end) : ()

Gray.(reshape((rbm.visibles.+1)./2, (28,28))')
end

##

train2(;rbm            = nothing,
        hidden_size    = 10,
        learning_rate  = 1,
        hm_epochs      = 1,
        batch_size     = 10,
        do_print       = true
       ) =
begin

    rbm == nothing ?
        rbm = RBM(in_size,hidden_size) :
            hidden_size = length(rbm.hiddens)

    do_print ? (@info "Training started. \nhidden size   $(hidden_size) \nlearning rate $(learning_rate) \nbatch size    $(batch_size)") : ()

    grad_norms      = []
    grad_sums       = []

    dev_grad_norms = []
    dev_grad_sums  = []

    do_print ? (begin

        dev_grads = batch_grads(rbm, data_dev)
        @info "initial dev stats \nnorm $(norm(dev_grads)./sqrt(length(rbm.weights))) \nsum $(sum(abs.(dev_grads))./length(rbm.weights))"

    end) : ()

    make_equalified_batches = ()->
    begin

        batch_size = int(batch_size/10) * 10
        class_sizes = [length(class_data) for class_data in data_train2]
        hm_batches = int(class_sizes[argmin(class_sizes)]/batch_size)
        batchify(vcat([sample for batch in collect(zip([choices(data_train2[class_ctr], hm_batches*int(batch_size/10)) for class_ctr in 1:10]...)) for sample in batch]), batch_size)

    end

    for ep in 1:hm_epochs

        total_grads = nothing

        for batch in make_equalified_batches()

            grads = batch_grads(rbm, batch)

            update_weights!(rbm, grads, learning_rate)

            total_grads == nothing ? total_grads = grads : total_grads += grads

        end

        dev_grads = batch_grads(rbm, data_dev)

        push!(grad_norms, norm(total_grads)./sqrt(length(rbm.weights)))
        push!(grad_sums, sum(abs.(total_grads))./length(rbm.weights))
        push!(dev_grad_norms, norm(dev_grads)./sqrt(length(rbm.weights)))
        push!(dev_grad_sums, sum(abs.(dev_grads))./length(rbm.weights))

        do_print ? println("Epoch $ep, train_sum $(round(grad_sums[end],digits=3)), dev_sum $(round(dev_grad_sums[end],digits=3))") : ()

    end

    min_train_norm = argmin(grad_norms)
    min_train_sum = argmin(grad_sums)
    min_dev_norm = argmin(dev_grad_norms)
    min_dev_sum = argmin(dev_grad_sums)

    p1 = plot(1:hm_epochs, grad_norms,     title="train_norm_$(hidden_size)_$(learning_rate)_$(batch_size)",xlabel="$(grad_norms[min_train_norm]) / $(min_train_norm)")
    p2 = plot(1:hm_epochs, grad_sums,      title="train_sum_$(hidden_size)_$(learning_rate)_$(batch_size)",xlabel="$(grad_sums[min_train_sum]) / $(min_train_sum)")
    p3 = plot(1:hm_epochs, dev_grad_norms, title="dev_norm_$(hidden_size)_$(learning_rate)_$(batch_size)",xlabel="$(dev_grad_norms[min_dev_norm]) / $(min_dev_norm)")
    p4 = plot(1:hm_epochs, dev_grad_sums,  title="dev_sum_$(hidden_size)_$(learning_rate)_$(batch_size)",xlabel="$(dev_grad_sums[min_dev_sum]) / $(min_dev_sum)")

    display(plot(p1,p2,p3,p4,layout=(2,2)))

rbm, [grad_norms,grad_sums,dev_grad_norms,dev_grad_sums]
end


text:
module IntPrint

export intprint

"""
`intprint(n)` converts an integer `n` to a `String` that has commas
inserted appropriately. For example:
```
julia> intprint(1234567)
"1,234,567"
```
Two named optional arguments may be passed:
+ `sep` can be set to give the separator string (default is `","`).
  For example, `intprint(12345,sep="_")` yields `"12_345"`.
+ `except4` can be set to `true` (default is `false`) so that four-digit
  numbers print without a separator. For example,
  `intprint(1234,except4=true)` yields `"1234"`.
"""
function intprint(n::Integer;sep::String=",",except4::Bool=false)::String
  # hand 0 as a special case (do we need to?)
  # if n==0
  #   return "0"
  # end

  prefix = ""
  if n < 0
    prefix = "-"
    n = -n
  end

  result = ""

  if except4 && n < 10_000
    return prefix * string(n)
  end

  while n >= 1000
    digs = lead0(string(n%1000))  # next 3 digits
    result = sep * digs * result
    n = div(n,1000)
  end

  digs = string(n)
  result = prefix * digs * result

  return result
end

"""
`lead0(s::String)` adds additions 0s on the left
to make this a length-3 string.
"""
function lead0(s::String)
  n = length(s)
  if n < 3
    s = "0"^(3-n) * s
  end
  return s
end

end # end of module


text:
# At each push event, the stage will be merge and
# it will be aknowlaged by synching the tokens
function _is_stage_up_to_day(gl::GitLink)

    stage_token = _get_stage_token(gl)
    push_token = _get_push_token(gl)

    isempty(stage_token) && return false
    isempty(push_token) && return false

    return stage_token == push_token
end

function _merge_stage(gl::GitLink)
    rdir = repo_dir(gl)
    sdir = stage_dir(gl)

    for src in _readdir(sdir; join = true)
        dest = replace(src, sdir => rdir)
        _cp(src, dest)
    end
end


text:
module ClickHouse

using CSV
using HTTP

import HTTP.URIs: URI

export
    # core.jl
    ClickHouseClient,
    query

include("core.jl")

end


text:
using BinaryProvider # requires BinaryProvider 0.3.0 or later
using CMakeWrapper: cmake_executable

include("compile.jl")

# Parse some basic command-line arguments
const verbose = "--verbose" in ARGS
const prefix = Prefix(joinpath(@__DIR__, "usr"))
products = [
    LibraryProduct(prefix, "libpng16", :libpng)
]

# Download binaries from hosted location
source_prefix = "https://sourceforge.net/projects/libpng/files/libpng16/1.6.37"

# Listing of files generated by BinaryBuilder:
download_info = Dict(
    Linux(:x86_64, libc=:glibc) => ("$source_prefix/libpng-1.6.37.tar.gz", "daeb2620d829575513e35fecc83f0d3791a620b9b93d800b763542ece9390fb4")
)

# Install unsatisfied or updated dependencies:
unsatisfied = any(!satisfied(p; verbose=verbose) for p in products)
dl_info = choose_download(download_info, platform_key_abi())
if dl_info === nothing && unsatisfied
    # If we don't have a compatible .tar.gz to download, complain.
    # Alternatively, you could attempt to install from a separate provider,
    # build from source or something even more ambitious here.
    error("Your platform (\"$(Sys.MACHINE)\", parsed as \"$(triplet(platform_key_abi()))\") is not supported by this package!")
end

# If we have a download, and we are unsatisfied (or the version we're
# trying to install is not itself installed) then load it up!
if unsatisfied || !isinstalled(dl_info...; prefix=prefix)
    # Download and install binaries
    compile(basename(dl_info[1]), dl_info[1], dl_info[2], prefix=prefix, verbose=verbose)
end

# check again whether the dependency is satisfied, which
# may not be true if dlopen fails due to a libc++ incompatibility
any(!satisfied(p; verbose=verbose) for p in products) && error("Unable to install libpng!")

# Write out a deps.jl file that will contain mappings for our products
write_deps_file(joinpath(@__DIR__, "deps.jl"), products, verbose=verbose)

max_stars_repo_name: FugroRoames/LibPNG

max_issues_repo_name: FugroRoames/LibPNG

max_forks_repo_name: FugroRoames/LibPNG


text:
using Test
using Random
using MarketEquilibria

@testset "Homogeneous Fisher model" begin
    samp = 5

    @testset "Linear" begin
        for _ in 1:samp
            (n, m) = rand(10:50, 2)

            endowments =  1 .+ randexp(n)
            A = .25 .+ 0.5 .* rand(n, m)    # Yields better conditioning
            supplies = 1 .+ rand(m)

            X, prices = homogeneousfisher(endowments, A, supplies, :linear)
            @test -endowments ≈ X * prices atol=1e-5
        end
    end

    @testset "Cobb–Douglas" begin
        for _ in 1:samp
            (n, m) = rand(10:50, 2)

            endowments =  1 .+ randexp(n)
            A = rand(n, m)
            supplies = 1 .+ rand(m)

            A ./= sum(A, dims=2)         # Required by Cobb–Douglas family
            X, prices = homogeneousfisher(endowments, A, supplies, :cobb_douglas)
            @test -endowments ≈ X * prices atol=1e-4
        end
    end

    @testset "Leontief" begin
        for _ in 1:samp
            (n, m) = rand(10:50, 2)

            endowments =  1 .+ randexp(n)
            A = .25 .+ 0.5 .* rand(n, m)
            supplies = 1 .+ rand(m)

            A ./= sum(A, dims=2)
            X, prices = homogeneousfisher(endowments, A, supplies, :leontief)
            @test -endowments ≈ X * prices atol=1e-4
        end
    end

    @testset "CES partial substitutes: ρ ∈ (0, 1)" begin
        for _ in 1:samp
            (n, m) = rand(10:50, 2)

            endowments =  1 .+ randexp(n)
            A = rand(n, m)
            supplies = 1 .+ rand(m)

            A ./= sum(A, dims=2)
            ρ = .97 * rand()        # Fails consistently for ρ  > .99, so the model uses
                                    # linear whenever ρ > .97.
            X, prices = homogeneousfisher(endowments, A, supplies, :CES, ρ)
            @test -endowments ≈ X * prices atol=1e-5
        end
    end

    @testset "CES complementarity fx: ρ < 0" begin
        for _ in 1:samp
            (n, m) = rand(10:50, 2)

            endowments =  1 .+ randexp(n)
            A = rand(n, m)
            supplies = 1 .+ rand(m)

            A ./= sum(A, dims=2)
            X, prices = homogeneousfisher(endowments, A, supplies, :CES, -randexp())
            @test -endowments ≈ X * prices atol=1e-5
        end
    end
end

@testset "Tâtonnement process for exchange economies" begin
    samp = 10

    for _ in 1:samp
        (n, m) = rand(5:10, 2)

        W = rand(n, m)      # Endowments
        A = rand(n, m)      # Utility params
        ρ = 1 .- randexp(n)

        prices, demands = tatonnement(W, A, ρ, nothing)

        @test W * prices ≈ demands * prices atol=1e-4
    end
end

@testset "Convex programs for exchange economies" begin
    samp = 10

    @testset "Linear" begin
        for _ in 1:samp
            (n, m) = rand(20:50, 2)

            endowments =  1 .+ randexp(n, m)
            A = randexp(n, m)
            A .*= (rand(n, m) .< .75)
            A ./= sum(A, dims=2)

            prices, demands = exchange(endowments, A, :linear)

            @test all(demands * prices .≥ endowments * prices .- 10e-2)
        end
    end

    @testset "CES" begin
        for _ in 1:samp
            (n, m) = rand(10:50, 2)

            endowments =  1 .+ randexp(n, m)
            A = rand(n, m)
            ρ = -rand(n)

            prices, demands = exchange(endowments, A, ρ)

            @test all(demands * prices .≥ endowments * prices .- 10e-4)
        end
    end
end

@testset "CES production and exchange" begin
    samp = 10

    for i in 1:samp
        n = rand(5:20)
        m = rand(4:n)       # m > n case doesn't converge well
        l = rand(3:m)       # l > m makes no sense

        endowments =  1 .+ randexp(n, m)
        A_consumers = rand(n, m)
        A_consumers ./= sum(A_consumers, dims=2)

        A_producers = rand(l, m)
        A_producers ./= sum(A_producers, dims=2)

        ρ_consumers = .1 .+ .8 * rand(n)
        ρ_producers = .1 .+ .8 * rand(l)

        # Index of what producers produce. Use perm here so none are redundant.
        o_producers = randperm(m)[1:l]

        # Can't use your output as an input. Unclear if this affects convergence.
        for (k, j) in enumerate(o_producers)
            A_producers
            A_producers[k, j] = 0
        end

        prices, demands, inputs = production(endowments, A_consumers, A_producers,
                                             ρ_consumers, ρ_producers, o_producers)

        @test all(inputs * prices .≈ endowments * prices)
    end
end

max_stars_repo_name: maxkapur/MarketEquilibria

max_issues_repo_name: maxkapur/MarketEquilibria

max_forks_repo_name: maxkapur/MarketEquilibria


text:
# Wrapper for current backend as `Symbol`
mutable struct CurrentADBackend
    backend::Symbol
end

const _current_adbackend = CurrentADBackend(:none)
const _adbackends = Symbol[]

"""
    adbackend(backend::Symbol = :default) -> Symbol

Check if autodiff `backend` is valid (see [`adbackends`](@ref)) and return
formatted. Note that `backend=:default` returns the current backend.
"""
adbackend() = _current_adbackend.backend
adbackend(backend::Symbol) = adbackend(Val(backend))
adbackend(backend::Val{:default}) = adbackend()

function adbackend(backend::Val{T}) where {T}
    T in _adbackends && return T
    throw(ArgumentError("Invalid AD backend $(T). Valid options are $(adbackends())."))
end

"""
    adbackend!(backend::Symbol)

Set current backend for autodiff to `backend`.
"""
function adbackend!(backend)
    backend = adbackend(backend)
    _current_adbackend.backend = backend
    return backend
end

"""
    adbackends() -> Vector{Symbol}

Get vector of currently valid AD backends.
"""
adbackends() = _adbackends

"""
    _gradient(f, x::Number, backend = adbackend()) -> Number
    _gradient(f, x::Array, backend = adbackend()) -> Array

Compute gradient of function `f` at point `x` using AD backend `backend`.
"""
_gradient(f, x) = _gradient(f, x, Val(adbackend()))
_gradient(f, x, backend::Symbol) = _gradient(f, x, Val(adbackend(backend)))

"""
    _jacobian(f, x::Array, backend = adbackend()) -> Array

Compute Jacobian matrix of function `f` at point `x` using AD backend `backend`.
Inputs and outputs of `f` are vectorized.
"""
_jacobian(f, x) = _jacobian(f, x, Val(adbackend()))
_jacobian(f, x, backend::Symbol) = _jacobian(f, x, Val(adbackend(backend)))

# Finite differences

push!(_adbackends, :finitedifferences)
adbackend!(:finitedifferences)

_default_fdm() = central_fdm(5, 1)

_gradient(f, x, backend::Val{:finitedifferences}) = _gradient(f, x, _default_fdm())

function _gradient(f, x, fdm::FiniteDifferences.FiniteDifferenceMethod)
    return FiniteDifferences.grad(fdm, f, x)[1]
end

_jacobian(f, x, backend::Val{:finitedifferences}) = _jacobian(f, x, _default_fdm())

function _jacobian(f, x, fdm::FiniteDifferences.FiniteDifferenceMethod)
    return FiniteDifferences.jacobian(fdm, f, x)[1]
end


text:
function notnull(s::AbstractString)
   if s == ""
       return false
   else
       return true
   end
end

max_stars_repo_name: atantos/RegTools

max_issues_repo_name: atantos/RegTools

max_forks_repo_name: atantos/RegTools


text:
export Sphere

# ---------- #
# - Sphere - #
# ---------- #

"""
    Sphere

Sphere is a primitive object.

### Fields:

* `center`   - Center of the Sphere in 3D world space
* `radius`   - Radius of the Sphere
* `material` - Material of the Sphere
"""
struct Sphere{C, P, Q, R, S, T, U} <: Object
    center::Vec3{C}
    radius::C
    material::Material{P, Q, R, S, T, U}
end

show(io::IO, s::Sphere) =
    print(io, "Sphere Object:\n    Center - ", s.center, "\n    Radius - ", s.radius[],
          "\n    ", s.material)

@diffops Sphere

function intersect(s::Sphere, origin, direction)
    b = dot(direction, origin - s.center)  # direction is a vec3 with array
    c = l2norm(s.center) .+ l2norm(origin) .- 2 .* dot(s.center, origin) .- (s.radius .^ 2)
    disc = (b .^ 2) .- c
    function get_intersections(x, y)
        t = bigmul(x + y) # Hack to split the 0.0 gradient to both. Otherwise one gets nothing
        if y > 0
            sqrty = sqrt(y)
            z1 = -x - sqrty
            z2 = -x + sqrty
            if z1 <= 0 && z2 > 0
                t = z2
            elseif z1 > 0
                t = z1
            end
        end
        return t
    end
    result = broadcast(get_intersections, b, disc)
    return result
end

get_normal(s::Sphere, pt, dir) = normalize(pt - s.center)


text:
# This file includes the testset for modulation pulses

@testset "PulseTestset" begin
    # Fields
    for type in [RectangularPulse, RaisedCosinePulse]
        @test :amplitude in fieldnames(type)
        @test :duration in fieldnames(type)
    end

    # Defaults
    pulse = RectangularPulse()
    @test pulse.amplitude ≈ 1.
    @test pulse.duration  ≈ 1.
    pulse = RaisedCosinePulse()
    @test pulse.amplitude ≈ 1.
    @test pulse.duration  ≈ 1.

    # Energy
    pulse = RectangularPulse(2., 4.)
    @test energy(pulse) ≈ 16.
    pulse = RaisedCosinePulse(2., 4.)
    @test energy(pulse) ≈ 4. * 3 /

    # Waveforms
    pulse = RectangularPulse(4., 1.)
    @test pulse(-1) ≈ 0.
    @test pulse(2.) ≈ 0.
    @test pulse(1/2) ≈ pulse.amplitude

    pulse = RaisedCosinePulse(4., 1.)
    @test pulse(-1.) ≈ 0.
    @test pulse(1/2) ≈ pulse.amplitude / 2 * (1 - cos(2π / pulse.duration * 1/2))
    @test pulse(2.) ≈ 0.
end # Test set.


text:
using CounterTools
using Test

# Function for compiling read/write tests
function compile_test_program(kernel; ntimes = 10, array_size = 100000000)
    commands = [
        "g++",
        "-march=native",
        "-mtune=native",
        "-mcmodel=large",
        "-DSTREAM_ARRAY_SIZE=$array_size",
        "-DFUNCTION=$kernel",
        "-DNTIMES=$ntimes",
        "-O3",
        "-fopenmp",
        joinpath(@__DIR__, "programs", "stream.cpp"),
        "-o",
        joinpath(@__DIR__, "programs", "stream"),
    ]
    run(`$commands`)
end

function run_test_program(program = "stream"; cpunodebind=0, membind=0)
    path = joinpath(@__DIR__,  "programs", program)
    run(`numactl --cpunodebind=$cpunodebind --membind=$membind $path`)
    return nothing
end

include("record.jl")
include("indexzero.jl")
include("utils.jl")
include("uncore/handle.jl")
include("uncore/cha.jl")
include("core.jl")
#include("uncore.jl")


text:
module Fracciones

export Fraccion, numerador, denominador, fraccion, reciproco, @fraccion

include("fraccion.jl")

end # module


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ 538f05be-5053-11eb-19a0-959e34e1c2a1
using Pkg, DrWatson

# ╔═╡ 5d84f90c-5053-11eb-076b-5f30fc9685e3
begin
	#@quickactivate "StatisticalRethinkingStan"
	using StanSample, StanOptimize
	using StatisticalRethinking
end

# ╔═╡ 235c5298-5053-11eb-0608-435d1aa4716c
md" ## Clip-07-13-14s.jl"

# ╔═╡ 861accd4-5053-11eb-0432-81db212f4f38
begin
	sppnames = [:afarensis, :africanus, :hapilis, :boisei, :rudolfensis, :ergaster, :sapiens]
	brainvol = [438, 452, 612, 521, 752, 871, 1350]
	masskg = [37, 35.5, 34.5, 41.5, 55.5, 61, 53.5]
	df = DataFrame(species = sppnames, brain = brainvol, mass = masskg)
	scale!(df, [:mass, :brain])
	df.brain_std = df.brain/maximum(df.brain)
end;

# ╔═╡ 290d19ae-59cf-11eb-1636-772efccc7cb9
df

# ╔═╡ 24328f84-541d-11eb-0ec3-4df5b8a8cb19
begin
	scatter(df.mass, df.brain, xlab="body mass [kg]", ylab="brain vol [cc]",
		lab="Observations")
	for (ind, species) in pairs(df.species)
		annotate!([(df[ind, :mass] + 1, df[ind, :brain] + 30,
			Plots.text(df[ind, :species], 6, :red, :right))])
	end
	plot!()
end

# ╔═╡ 5710f362-5a99-11eb-1282-ad61ea6dde7e
md"
!!! note

	Below results need to be further studied. By manipulating `sigma` the overall results have a similar trend as in the book, but I'm not sure I trust this approach.
"

# ╔═╡ 9457b384-5a97-11eb-3a64-39599bd271a5
sig = 0.1

# ╔═╡ ee275e7a-5067-11eb-325b-7760b758e85e
stan7_2 = "
data {
	 int < lower = 1 > N; 			// Sample size
	 int < lower = 1 > K;			// Degree of polynomial
	 vector[N] brain; 				// Outcome
	 matrix[N, K] mass; 			// Predictor
}

parameters {
	real a;                        // Intercept
	vector[K] b;                  // K slope(s)
	real log_sigma;
}

transformed parameters {
    vector[N] mu;
    mu = a + mass * b;
}

model {
	a ~ normal(0.5, 1);
	b ~ normal(0, 10);
	brain ~ normal(mu , $(sig));
}
generated quantities {
	vector[N] log_lik;
	real sigma;
	for (i in 1:N)
		log_lik[i] = normal_lpdf(brain[i] | mu[i], $(sig));
	sigma = $(sig);
}
";

# ╔═╡ e1bd7ef2-5a98-11eb-19a8-f32fcb3d7fb6
sig6 = 0.001

# ╔═╡ f2898fe6-5a98-11eb-10f0-db1c4fe19bff
stan7_6 = "
data {
	 int < lower = 1 > N; 			// Sample size
	 int < lower = 1 > K;			// Degree of polynomial
	 vector[N] brain; 				// Outcome
	 matrix[N, K] mass; 			// Predictor
}

parameters {
	real a;                        // Intercept
	vector[K] b;                  // K slope(s)
	real log_sigma;
}

transformed parameters {
    vector[N] mu;
    mu = a + mass * b;
}

model {
	a ~ normal(0.5, 1);
	b ~ normal(0, 10);
	brain ~ normal(mu , $(sig6));
}
generated quantities {
	vector[N] log_lik;
	real sigma;
	for (i in 1:N)
		log_lik[i] = normal_lpdf(brain[i] | mu[i], $(sig6));
	sigma = $(sig6);
}
";

# ╔═╡ f59ecd72-60ab-11eb-2fe2-c3f7159e2080
begin
	loo = Vector{Float64}(undef, 6)
	loos = Vector{Vector{Float64}}(undef, 6)
	pk = Vector{Vector{Float64}}(undef, 6)
	deviance = Vector{Float64}(undef, 6)
end;

# ╔═╡ 95cbb456-520e-11eb-2e0a-1d173ddc9dec
begin
	tmpdir = joinpath(projectdir(), "tmp")
	lppd_res = Matrix{Float64}(undef, 6, 7)

	for K in 1:6
		N = size(df, 1)
		mass = create_observation_matrix(df.mass_s, K)
		data = (N = N, K = K, brain = df.brain_std, mass = mass)

		# `sigma` should really be `exp(log_sigma)`!

		if K < 6
			m7_2s = SampleModel("m7.2s", stan7_2; tmpdir)
		else
			m7_2s = SampleModel("m7.2s", stan7_6; tmpdir)
		end
		rc7_2s = stan_sample(m7_2s; data=data)

		if success(rc7_2s)
			nt7_2s = read_samples(m7_2s, :namedtuple)
			post7_2s_df = read_samples(m7_2s, :dataframe)
		end

		log_lik = nt7_2s.log_lik'
		lppd_res[K, :] = lppd(log_lik)
		loo[K], loos[K], pk[K] = psisloo(log_lik)

		lp = logprob(post7_2s_df, mass, df.brain_std, K)
		deviance[K] = -2sum(lppd(lp))
	end
	deviance
end

# ╔═╡ 58eaed70-6809-11eb-192f-2b1edb51a17e
-deviance/2

# ╔═╡ 1ab3e174-59e0-11eb-3db8-4d91bc007f0a
lppd_res'

# ╔═╡ 6441d73c-666b-11eb-2fa8-a920bcbf2380
sum(lppd_res', dims=1)

# ╔═╡ 946d5e90-5fe5-11eb-1ec8-ed64a60550e3
md"
!!! note

	Take a look at these runs using `psisloo()` from PSIS.jl. This is explained later in chapter 7.
"

# ╔═╡ ef3da47a-60ac-11eb-39b9-f3fe46d314d5
loo

# ╔═╡ 19541136-5fe5-11eb-00a2-db4ef2435713
sum.(loos)

# ╔═╡ 1ed9cb28-5fe5-11eb-0e3e-29e68eb14885
pk_qualify(pk[1])

# ╔═╡ 3ae2831e-5fe5-11eb-24ec-d36a56204b4e
pk_plot(pk[1]; title="PSIS diagnostic plot (K = 1)")

# ╔═╡ 5246f7ba-5f7f-11eb-2a0b-b77c47d6bdcc
begin
	fig = Vector{Plots.Plot{Plots.GRBackend}}(undef, 6)
	for i in 1:6
		fig[i] = pk_plot(pk[i], title="K = $i", leg=false)
	end
	plot(fig..., layout=(3,2))
end

# ╔═╡ 91e9af0a-5065-11eb-212c-f751fd114263
md" ## End of clip-07-13-14s.jl"

# ╔═╡ Cell order:
# ╟─235c5298-5053-11eb-0608-435d1aa4716c
# ╠═538f05be-5053-11eb-19a0-959e34e1c2a1
# ╠═5d84f90c-5053-11eb-076b-5f30fc9685e3
# ╠═861accd4-5053-11eb-0432-81db212f4f38
# ╠═290d19ae-59cf-11eb-1636-772efccc7cb9
# ╠═24328f84-541d-11eb-0ec3-4df5b8a8cb19
# ╟─5710f362-5a99-11eb-1282-ad61ea6dde7e
# ╠═9457b384-5a97-11eb-3a64-39599bd271a5
# ╠═ee275e7a-5067-11eb-325b-7760b758e85e
# ╠═e1bd7ef2-5a98-11eb-19a8-f32fcb3d7fb6
# ╠═f2898fe6-5a98-11eb-10f0-db1c4fe19bff
# ╠═f59ecd72-60ab-11eb-2fe2-c3f7159e2080
# ╠═95cbb456-520e-11eb-2e0a-1d173ddc9dec
# ╠═58eaed70-6809-11eb-192f-2b1edb51a17e
# ╠═1ab3e174-59e0-11eb-3db8-4d91bc007f0a
# ╠═6441d73c-666b-11eb-2fa8-a920bcbf2380
# ╟─946d5e90-5fe5-11eb-1ec8-ed64a60550e3
# ╠═ef3da47a-60ac-11eb-39b9-f3fe46d314d5
# ╠═19541136-5fe5-11eb-00a2-db4ef2435713
# ╠═1ed9cb28-5fe5-11eb-0e3e-29e68eb14885
# ╠═3ae2831e-5fe5-11eb-24ec-d36a56204b4e
# ╠═5246f7ba-5f7f-11eb-2a0b-b77c47d6bdcc
# ╟─91e9af0a-5065-11eb-212c-f751fd114263


text:
module FastIO
const FastReadInt = Union{Type{Int16}, Type{UInt16}, Type{Int32}, Type{UInt32}, Type{Int64}, Type{UInt64}}
const FastReads   = Union{FastReadInt, Type{Float16}, Type{Float32}, Type{Float64}}

#=TO DO:
fastread Char
IBM-Float
=#

# =====================================================================
#=
  This section: file position commands

  No rewrite needed in 1.3: mark, reset, seekstart
=#

# fastpos
fastpos(io::IO) = io.ptr-1
fastpos(io::IOStream) = ccall(:ios_pos, Int64, (Ptr{Cvoid},), io.ios)

# fasteof
fasteof(io::IO) = (io.ptr-1 == io.size)
fasteof(io::IOStream) = Bool(ccall(:ios_eof_blocking, Cint, (Ptr{Cvoid},), io.ios))

# fastseek
fastseek(io::IO, p::Integer) = seek(io, p)
function fastseek(io::IOStream, n::Integer)
  ccall(:ios_seek, Int64, (Ptr{Cvoid}, Int64), io.ios, n)
  return nothing
end

# fastskip
fastskip(io::IO, n::Integer) = skip(io, n)
function fastskip(io::IOStream, n::Integer)
  ccall(:ios_skip, Int64, (Ptr{Cvoid}, Int64), io.ios, n)
  return nothing
end

# fastseekend
fastseekend(io::IO) = (io.ptr = io.size+1)
fastseekend(io::IOStream) = ccall(:ios_seek_end, Int64, (Ptr{Cvoid},), io.ios)

# =====================================================================
#= This section: file read commands

fastread(io) returns a UInt8
fastread(io, T) returns a T

=#
fastread!(io::IO, buf::Array{T}) where T = read!(io, buf)
fastread!(io::IOStream, buf::Array{T}) where T = @GC.preserve buf ccall(:ios_readall, Csize_t, (Ptr{Cvoid}, Ptr{Cvoid}, Csize_t), io, pointer(buf), sizeof(buf))
# fastread!(io::IOStream, buf::Array{UInt8,1}) = ccall(:ios_readall, Csize_t, (Ptr{Cvoid}, Ptr{Cvoid}, Csize_t), io.ios, pointer(buf, 1), sizeof(buf))

# fastread
fastread(io::IO) = read(io, UInt8)
fastread(io::IO, T::FastReadInt) = read(io, T)
fastread(io::IO, n::Integer) = read(io, n)

# IOStream methods avoid locking in 1.3
fastread(io::IOStream) = ccall(:ios_getc, Cint, (Ptr{Cvoid},), io.ios) % UInt8
function fastread(io::IOStream, n::Integer)
  buf = Vector{UInt8}(undef, n)
  ccall(:ios_readall, Csize_t, (Ptr{Cvoid}, Ptr{Cvoid}, Csize_t), io.ios, pointer(buf, 1), n)
  return buf
end

# Working solution 2020-02-10T16:00:00
function fastread(io::IOStream, T::FastReadInt)
  if VERSION > v"1.2.0"
    # ccall(:ios_flush, Cint, (Ptr{Cvoid},), io.ios)
    ccall(:ios_readprep, UInt64, (Ptr{Cvoid}, Csize_t), io.ios, sizeof(T))
  end
  return (ccall(:jl_ios_get_nbyte_int, UInt64, (Ptr{Cvoid}, Csize_t), io.ios, sizeof(T)) % T)
end
#= I *don't* like this. jl_ios_get_nbyte_int seems more likely to change than
   ios_getc, and version-dependent "if" loop is undesirable =#

# function fastread1(io::IOStream, ::Type{Int16})
#   ccall(:ios_flush, Cint, (Ptr{Cvoid},), io.ios)
#   x  = ccall(:ios_getc, Cint, (Ptr{Cvoid},), io.ios) | (ccall(:ios_getc, Cint, (Ptr{Cvoid},), io.ios) << 0x08)
#   return signed(x % UInt16)
# end
#
fastread(io::IO, ::Type{Bool})    = (fastread(io) != 0x00)
fastread(io::IO, ::Type{UInt8})   = fastread(io)
fastread(io::IO, ::Type{Int8})    = signed(fastread(io))
fastread(io::IO, ::Type{Float16}) = Base.reinterpret(Float16, fastread(io, Int16))
fastread(io::IO, ::Type{Float32}) = Base.reinterpret(Float32, fastread(io, Int32))
fastread(io::IO, ::Type{Float64}) = Base.reinterpret(Float64, fastread(io, Int64))

# use if creating a new array during the call; I have yet to read anything but a 1d array in any format
function fastread(io::IO, T::Type, n::Integer)
  a = Array{T, 1}(undef, n)
  fastread!(io, a)
  return a
end

# TO DO
fastread(io::IO, ::Type{Char}) = read(io, Char)

# =====================================================================
# fast_readbytes!
fast_readbytes!(io::IO, buf::Array{UInt8,1}, nb::Integer) = readbytes!(io, buf, nb)
function fast_readbytes!(io::IOStream, buf::Array{UInt8,1}, nb::Integer)
    olb = lb = length(buf)
    nr = 0
    GC.@preserve buf while nr < nb
        if lb < nr+1
            lb = max(65536, (nr+1) * 2)
            resize!(buf, lb)
        end
        nr += Int(ccall(:ios_readall, Csize_t, (Ptr{Cvoid}, Ptr{Cvoid}, Csize_t),
                        io.ios, pointer(buf, nr+1), min(lb-nr, nb-nr)))
        fasteof(io) && break
    end
    if lb > olb && lb > nr
        resize!(buf, nr) # shrink to just contain input data if was resized
    end
    return nr
end
# was (and crashed with)
# = @GC.preserve buf ccall(:ios_readall, Csize_t, (Ptr{Cvoid}, Ptr{Cvoid}, Csize_t), io.ios, pointer(buf, 1), nb)

# fast_readline
fast_readline(io::IO) = readline(io)
fast_readline(io::IOStream) = ccall(:jl_readuntil, Ref{String}, (Ptr{Cvoid}, UInt8, UInt8, UInt8), io.ios, '\n', 1, 2)

# fast_unsafe_read
fast_unsafe_read(io::IO, p::Ptr{UInt8}, nb::Integer) = unsafe_read(io, p, nb)
fast_unsafe_read(io::IOStream, p::Ptr{UInt8}, nb::Integer) = ccall(:ios_readall, Csize_t, (Ptr{Cvoid}, Ptr{Cvoid}, Csize_t), io, p, nb)

# =====================================================================
export
  FastReadInt,
  FastReads,
  fast_readbytes!,
  fast_readline,
  fast_unsafe_read,
  fasteof,
  fastpos,
  fastread!,
  fastread,
  fastseek,
  fastseekend,
  fastskip
end


text:
"""
    This revert function reverses the internal implementation of the
`ScaleTransform` operation from `KernelFunctions.jl`. This is to preserve the
previous convention of this package.
"""
_revert(x) = sqrt(2.0 * x)

"""
    _build_kernel_matrix(x; kwargs...) -> AbstractMatrix
    _build_kernel_matrix(x, y; kwargs...) -> AbstractMatrix

It builds a kernel matrix using a user specified kernel. When just one matrix is passed,
`x`, the resulting matrix has the same size as `x`.
When both `x` and `y` are passed, the resulting matrix has the size
(`size(x, 2)`, `size(y, 2)`).

The matrix is built using the `kernelmatrix!` functions from `KernelFunctions.jl`.

# Arguments
- `x::AbstractMatrix`: The data matrix with the training instances.
- `y::AbstractMatrix`: Another data matrix, normally used for the prediction step in the classification implementation.

# Keywords
- `kernel::Symbol=:rbf`: The kernel to be used. The available options are `:rbf` for a RBF type kernel; `:linear`, for a classic linear kernel; and `:poly` for a polynomial kernel. The degree of the polynomial is handled with the `degree` keyword argument.
- `sigma::Float64`: The hyperparameter for the RBF kernel.
- `degree::Float64`: The degree of the polynomial kernel if "poly" is chosen.

# Returns
- `kern_mat`: The kernel matrix.
"""
function _build_kernel_matrix(x; kwargs...)
    # Extract the matrix for the keyword arguments
    kernel = kwargs[:kernel]

    if kernel == :rbf
        # Create the kernel with the corresponding scale
        t = ScaleTransform(_revert(kwargs[:sigma]))
        κ = transform(SqExponentialKernel(), t)
        kern_mat = kernelmatrix(κ, x)
    elseif kernel == :linear
        κ = LinearKernel()
        kern_mat = kernelmatrix(κ, x)
    elseif kernel == :poly
        # Create the kernel with the corresponding degree
        κ = PolynomialKernel(; degree=kwargs[:degree], c=0.0)
        kern_mat = kernelmatrix(κ, x)
    end

    return kern_mat
end

function _build_kernel_matrix(x, y; kwargs...)
    # Check that the first dimension is the same
    @assert size(x, 1) == size(y, 1)

    # Extract the matrix for the keyword arguments
    kernel = kwargs[:kernel]

    if kernel == :rbf
        # Create the kernel with the corresponding scale
        t = ScaleTransform(_revert(kwargs[:sigma]))
        κ = transform(SqExponentialKernel(), t)
        kern_mat = kernelmatrix(κ, x, y)
    elseif kernel == :linear
        κ = LinearKernel()
        kern_mat = kernelmatrix(κ, x, y)
    elseif kernel == :poly
        # Create the kernel with the corresponding degree
        κ = PolynomialKernel(; degree=kwargs[:degree], c=0.0)
        kern_mat = kernelmatrix(κ, x, y)
    end

    return kern_mat
end

"""
    Takes a Support Vector Machine type and converts some of its attributes to a
dictionary that makes it easier to handle as keyword arguments.
"""
_kwargs2dict(svm) = Dict(:kernel => svm.kernel, :sigma => svm.σ, :degree => svm.degree)

"""
    Function to find all the instances in an array `y` that are
equal to some value `k`. It returns a copy of the array and the indices where the
condition is met.
"""
function _find_and_copy(k, y)
    indices = findall(isequal(k), y)

    return (copy(y[indices]), indices)
end

"""
    This function takes in a matrix `x` and does the following logic:
- Obtains the unique elements from the matrix.
- Counts how many occurrences of each element happen in the array.
- Using `argmax`, the indices where this condition is met are extracted.
- Finally, we only need the first dimension index, so we extract it as such.

Essentially, this is a voting scheme for the multiclass classification problem.
In the case of a tie, the smallest index is always chosen, i.e. 1. This is not the best
strategy, but it is after the following paper:

Chih-Wei Hsu and Chih-Jen Lin (2002) ‘A comparison of methods for multiclass support vector machines’, IEEE Transactions on Neural Networks, 13(2), pp. 416. doi: 10.1109/72.991427.

Where it says the following quote:

> [...] Then we predict is in the class with the largest vote. The
> voting approach described above is also called the “Max Wins”
> strategy. In case that two classes have identical votes, though it
> may not be a good strategy, now we simply select the one with
> the smaller index. [...]
"""
function _predictions_by_votes(x)
    unique_elements = unique(x) |> sort
    predictions = zeros(size(x, 2))
    counts = map(z -> count(==(z), x, dims=1), unique_elements)
    largest_values = argmax(vcat(counts...), dims=1)

    for (i, l) in enumerate(largest_values)
        predictions[i] = l[1]
    end

    return predictions
end

max_stars_repo_name: gmagannaDevelop/Elysivm

max_issues_repo_name: gmagannaDevelop/Elysivm

max_forks_repo_name: gmagannaDevelop/Elysivm


text:
# vim: set et ts=2 sw=2;
# Non-parametric optimizers
println("# peek.jl ...")
include("lib.jl")

#### Config
@with_kw mutable struct It
  data = (file="auto93.csv", dir="data",some=128,best=.75)
  char = (skip='?',less='-',more='+',num=':',klass='!')
  str  = (skip="?")
  some = (max=256,bins=.5, cohen=0.3, trivial=1.05)
  divs = (few=126)
  seed = 1
end

#### Globals
it=It()
Random.seed!(it.seed)
no = nothing

#-------------------------------------------------------------------
#### Columns
# Count all the symbols. Keep a sample of the numbers.
@with_kw mutable struct Some pos=0;txt="";w=1;n=0;_all=[];ok=true end
@with_kw mutable struct Sym  pos=0;txt="";w=1;n=0;seen=Dict() end
@with_kw mutable struct Skip pos=0;txt="";w=1;n=0 end

inc!(i,x::Array) = begin [inc!(i,y) for y in x]; i end
inc!(i,x) = begin
  x==it.char.skip ? x : (i.n += 1; inc1!(i,x)); x end

inc1!(i::Skip, x) = i
inc1!(i::Sym,  x) = begin
  new = i.seen[x] = 1 + get(i.seen,x,0)
  if new > i.n i.n, i.mode = now, x end end
inc1!(i::Some, x) = begin
  m = length(i._all)
  if m < it.some.max
    i.ok=false; push!(i._all, x);
  elseif rand() < m/i.n
    i.ok=false; i._all[int(m*rand())+1]=x end  end

mid(i::Sym)   = i.mode
mid( i::Some) = per(all(i),.5)

var(i::Sym) = sum(-v/i.n*log(v/i.n)/log(2) for (_,v) in i.seen)
var(i::Some,a=all(i)) = (per(a,.9) - per(a,.1)) / 2.56

norm(i::Some,x, a=all(i)) =
  x==it.char.skip ? x : (x-a[1])/(a[end]-a[1]+1E-32)

all(i::Some)  = begin
  i._all = i.ok ? i._all : sort(i._all)
  i.ok = true
  i._all end

combine(i::Some,j::Some, small,lo) = begin
  mi, mj = mid(i), mid(j)
  if (abs(mi) - abs(mj) < small) || (mi < lo && mi < lo)
    inc!(Some(), [i._all;j._all]) end end

#-------------------------------------------------------------------
#### Table
# Load rows, Summarize the columns.
@with_kw mutable struct Table ys=[]; xs=[]; rows=[]; cols=[] end
@with_kw mutable struct Row   has=[]; gt=0; tag=no end

function data(file; t=Table())
  for a in csv(it.data.dir * "/" * file) inc!(t,a)
  t end

inc!( t::Table,a) = length(t.cols)==0 ? t.cols=cols(t,a) : push!(t.rows, cells(t,a))
cells(t::Table,a) = Row(has= [inc!(c, a[c.pos]) for c in t.cols])
cols(t::Table,a)  = [col(txt=txt, pos=pos) for (pos,txt) in enumerate(a)]

function  col(t::Table;txt="", pos=0, c=it.char) = begin
  what = c.less in txt||c.more in txt||c.num in txt ? Some : Sym
  what = c.skip in txt ? Skip : what
  now  = what(txt=txt, pos=pos, w= c.less in txt ? -1 : 1)
  at   = c.less in txt||c.more in txt||c.klass in txt ? t.ys : t.xs
  push!(at, now)
  now
end

#-------------------------------------------------------------------
#### Classify
# Score rows by how they are better than others.
# Tag the `best` percent rows.
function tag(t::Table)
  function gt(row1, row2)
    "Zitler's continous domination predicate (from IBEA, 2005)."
    s1, s2, n = 0, 0, length(t.ys)
    for col in t.ys
      a, b = row1.has[col.pos], row2.has[col.pos]
      a, b = norm(col, a), norm(col, b)
      s1  -= ℯ ^ (col.w * (a - b) / n)
      s2  -= ℯ ^ (col.w * (b - a) / n) end
    s1/n < s2/n
  end
  n = it.data.some
  print(":::",any(t.rows))
  for row in t.rows
    row.gt = sum(gt(row, any(t.rows)) for _ in 1:n)/n  end
  for (n, row) in enumerate(sort(t.rows, by = x -> x.gt))
    row.tag = n > length(t.rows) * it.data.best end
  t end

#-------------------------------------------------------------------
#### Discretize
# Score rows by how they are better than others.
# Tag the `best` percent rows.
@with_kw mutable struct Span
   lo=typemax(Int); hi=typemin(Int); _has=[] end

o(i::Span) =
  if     i.lo <= typemin(Int) "<= $(i.hi)"
  elseif i.hi <= typemax(Int) ">= $(i.lp)"
  else   "[$(i.lo)..$(i.hi)]" end

#function discretize(x::Some, t:Table)
  # function div(xsmall,ysmall,ymin,xy)
  #   n = length(xy) * the data.xchop
  #   #while n < 4 && n <= length(xy)/2 n+= 1.2 end
  #   n, now, tmp = int(n), int(n), []
  #   b4, span    = 1, Span(lo=xy[1][1])
  #   while now < length(xy) - n
  #     x=xy[now][1]
  #     span.hi = x
  #     now += 1
  #     if (now - b4 => n    && now < len(xy) - 2 &&
  #         x != xy[now][1]  && span.hi - span.lo > xsmall)
  #       span._has= [z[2] for z in xy[b4:now]]
  #       tmp += [span]
  #       span = Span(lo=xy[now][1])
  #       b4   = now
  #       now += n end end
  #    tmp += [Span(lo=xy[b4][1],hi=xy[end][1],
  #                 _has=[z[2] for z in xy[b4:]])]
  #    merge(tmp,ysmall,ymin) end
  #
  # function merge(b4, ysmall, ymin)
  #   j, now = 0, []
  #   while j < length(b4)
  #     a=b4[j]
  #     if j < length(b4)
  #       b = b4[j+1]
  #       c  = combine(a._has,b._has, ysmall, ymin)
  #       if c != nothing
  #         now += [Span(lo=a.lo,hi=b.hi,_has=c)] end
  #   length(now) == length(b4) ? b4 : merge(now,ymin,small) end
  #
  # for row in rows
  #   inc!(y, row.gt)
  #   push!(xy, (row.has[x.pos], row.gt)) end
  #   tmp = div(var(x)*it.some.xsmall,
  #           var(y)*it.some.ysmall,
  #           per(all(y), it.data.best),
  #           sort(xy)) end
  # tmp[1].lo = typemin(Int)
  # tmp[end].hi = typemax(Int)
  # tmp end

max_stars_repo_name: timm/peek

max_issues_repo_name: timm/peek

max_forks_repo_name: timm/peek


text:
function get_inits(::MIMIXNoFactors, params::Dict{Symbol, Any}, data::Dict{Symbol, Any})
    θ_init = iclr(proportionalize(data[:Y]))
    ω_init = params[:ω] == "ones" ? ones(data[:K], data[:p]) : zeros(data[:K], data[:p])
    π_init = params[:π] == "even" ? [0.5 for _ in 1:data[:p]] : params[:π]
    β_init = params[:β] == "zeros" ? zeros(data[:K], data[:p]) : params[:β]
    γ_init = params[:γ] == "zeros" ? zeros(data[:K], data[:q]) : params[:γ]
    β_var_init = params[:β_var] == "ones" ? ones(data[:p]) : params[:β_var]
    γ_var_init = params[:γ_var] == "ones" ? ones(data[:num_blocking_factors]) : params[:γ_var]

    inits = Dict{Symbol, Any}(
        :Y => data[:Y],
        :θ => θ_init,
        :θ_var => vec(var(θ_init, dims=1)),
        :μ => vec(mean(θ_init, dims=1)),
        :μ_var => params[:μ_var],
        :ω => ω_init,
        :π => π_init,
        :β_full => β_init,
        :γ => γ_init,
        :β_var => β_var_init,
        :γ_var => γ_var_init
    )
    return inits
end

function get_model(::MIMIXNoFactors, monitor::Dict{Symbol, Any}, hyper::Dict{Symbol, Any})

    model = Model(
        # High-dimensional counts
        Y = Stochastic(2,
            (m, ϕ, N) -> MultivariateDistribution[
                Multinomial(m[i], ϕ[i, :]) for i in 1:N
            ],
            monitor[:Y]
        ),

        ϕ = Logical(2,
            (θ) -> clr(θ),
            monitor[:ϕ]
        ),

        # Sample-specific random effects

        θ = Stochastic(2,
            (θ_mean, θ_var, N) -> MultivariateDistribution[
                MvNormal(θ_mean[i, :], sqrt.(θ_var)) for i in 1:N
            ],
            monitor[:θ]
        ),

        θ_mean = Logical(2,
            (μ, Xβ, sumZγ) -> Xβ + sumZγ .+ transpose(μ),
            monitor[:θ_mean]
        ),

        θ_var = Stochastic(1,
            () -> InverseGamma(1, 1),
            monitor[:θ_var]
        ),

        # Population mean

        μ = Stochastic(1,
            (K, μ_var) -> MvNormal(K, sqrt(μ_var)),
            monitor[:μ]
        ),

        μ_var = Stochastic(
            () -> InverseGamma(1, 1),
            monitor[:μ_var]
        ),

        # Blocking factor random effects

        sumZγ = Logical(2,
            (Zγ) -> transpose(squeezesum(Zγ, 3)),
            false
        ),

        Zγ = Logical(3,
            (Z, γ) -> γ[:, Z],
            false
        ),

        γ = Stochastic(2,
            (q, K, γ_var, blocking_factor) -> UnivariateDistribution[
                Normal(0.0, sqrt(γ_var[blocking_factor[r]])) for k in 1:K, r in 1:q
            ],
            monitor[:γ]
        ),

        γ_var = Stochastic(1,
            () -> InverseGamma(1, 1),
            monitor[:γ_var]
        ),

        # Fixed treatment effects w/ spike-and-slab variable selection

        Xβ = Logical(2,
            (X, β) -> X * transpose(β),
            false
        ),

        β = Logical(2,
            (β_full, ω) -> ω .* β_full,
            monitor[:β]
        ),

        β_full = Stochastic(2,
            (p, K, β_var) -> MultivariateDistribution[
                MvNormal(zeros(p), sqrt.(β_var)) for k in 1:K
            ],
            monitor[:β_full]
        ),

        β_var = Stochastic(1,
            () -> InverseGamma(1, 1),
            monitor[:β_var]
        ),

        ω = Stochastic(2,
            (p, K, π) -> UnivariateDistribution[
                Bernoulli(π[j]) for k in 1:K, j in 1:p
            ],
            monitor[:ω]
        ),

        π = Stochastic(1,
            (K) -> Beta(1, K),
            monitor[:π]
        )
    )

    samplers = [
        Sampler(:θ, (Y, m, θ, θ_mean, θ_var, N, K) ->
            begin
                acc = 0
                for i in 1:N
                    eps = rand(Exponential(hyper[:hmc_epsilon]))
                    yi = Y[i, :]
                    mi = m[i]
                    u1 = θ[i, :]
                    θi_mean = θ_mean[i, :]
                    logf0, grad0 = logf1, grad1 = logfgrad(u1, θi_mean, θ_var, yi, mi)
                    v0 = v1 = randn(length(u1))
                    v1 += 0.5 * eps * grad0
                    for step in 1:hyper[:hmc_steps]
                        u1 += eps * v1
                        logf1, grad1 = logfgrad(u1, θi_mean, θ_var, yi, mi)
                        v1 += eps * grad1
                    end
                    v1 -= 0.5 * eps * grad1
                    v1 *= -1.0
                    Kv0 = 0.5 * sum(abs2, v0)
                    Kv1 = 0.5 * sum(abs2, v1)
                    if rand() < exp((logf1 - Kv1) - (logf0 - Kv0))
                        acc += 1
                        θ[i, :] = u1
                    end
                end
                if hyper[:hmc_verbose]
                    println("HMC accept rate: $(round(100acc/N, digits=1))%, ")
                end
                θ
            end
        ),

        Sampler(:μ, (θ, θ_var, μ, μ_var, Xβ, sumZγ, N, K) ->
            begin
                Sigma = 1 ./ (N ./ θ_var .+ 1 / μ_var)
                mu = Sigma .* vec(sum(θ - Xβ - sumZγ, dims=1)) ./ θ_var
                for k in 1:K
                    μ[k] = rand(Normal(mu[k], sqrt(Sigma[k])))
                end
                μ
            end
        ),

        Sampler([:γ], (θ, θ_var, Xβ, μ, γ, γ_var, Z, sumZγ, q, K) ->
            begin
                A = θ - Xβ - sumZγ .- transpose(μ)
                for j in 1:size(Z, 2)
                    z = Z[:, j]
                    A += transpose(γ[:, z])
                    for r in unique(z)
                        in_block = z .== r
                        B = vec(sum(A[in_block, :], dims=1))
                        Sigma = 1 ./ (sum(in_block) ./ θ_var .+ 1 / γ_var[j])
                        mu = Sigma .* B ./ θ_var
                        for k in 1:K
                            γ[k, r] = rand(Normal(mu[k], sqrt(Sigma[k])))
                        end
                    end
                    A -= transpose(γ[:, z])
                end
                γ
            end
        ),

        Sampler(:β_full, (θ, θ_var, sumZγ, μ, β_full, β_var, X, ω, K, p) ->
            begin
                A = θ - sumZγ .- transpose(μ)
                for k in 1:K
                    Xtilde = X * Diagonal(ω[k, :])
                    Sigma = inv(cholesky(Hermitian(transpose(Xtilde) * Xtilde ./ θ_var[k] + Diagonal(1 ./ β_var))))
                    mu = Sigma * ((transpose(Xtilde) * A[:, k]) ./ θ_var[k])
                    β_full[k, :] = rand(MvNormal(mu, Sigma))
                end
                β_full
            end
        ),

        Sampler(:ω, (ω, θ, θ_var, sumZγ, X, Xβ, μ, β, β_full, π, p, K) ->
            begin
                D = θ - sumZγ - Xβ .- transpose(μ)
                for k in 1:K
                    d = D[:, k]
                    for j in 1:p
                        x = X[:, j]
                        d += x .* β[k, j]
                        lprob_in = log(π[j]) - 0.5sum(abs2, d - x .* β_full[k, j]) / θ_var[k]
                        lprob_out = log(1 - π[j]) - 0.5sum(abs2, d) / θ_var[k]
                        diff = min(lprob_in - lprob_out, 10.0)
                        prob_in = exp(diff) / (1 + exp(diff))
                        ω[k, j] = rand(Bernoulli(prob_in))
                        d -= x .* (ω[k, j] * β_full[k, j])
                    end
                end
                ω
            end
        ),

        Sampler(:π, (π, ω, K, p) ->
            begin
                S = vec(sum(ω, dims=1))
                c0, d0 = params(π.distr)
                c = c0 .+ S
                d = d0 + K .- S
                [rand(Beta(c[j], d[j])) for j in 1:p]
            end
        ),

        Sampler(:μ_var, (μ, μ_var, K) ->
            begin
                u = 0.5K + shape(μ_var.distr)
                v = 0.5sum(abs2, μ) + scale(μ_var.distr)
                rand(InverseGamma(u, v))
            end
        ),

        Sampler(:γ_var, (γ, γ_var, K, q, blocking_factor) ->
            begin
                num_blocking_factor = length(γ_var)
                u = fill(shape(γ_var.distr), num_blocking_factor)
                v = fill(scale(γ_var.distr), num_blocking_factor)
                G = sum(abs2, γ, dims=1)
                for r in 1:q
                    u[blocking_factor[r]] += 0.5K
                    v[blocking_factor[r]] += 0.5G[r]
                end
                for i in 1:num_blocking_factor
                    γ_var[i] = rand(InverseGamma(u[i], v[i]))
                end
                γ_var
            end
        ),

        Sampler(:β_var, (β, β_var, ω, p) ->
            begin
                u = 0.5vec(sum(ω, dims=1)) .+ shape(β_var.distr)
                v = 0.5vec(sum(abs2, β, dims=1)) .+ scale(β_var.distr)
                [rand(InverseGamma(u[j], v[j])) for j in 1:p]
            end
        ),

        Sampler(:θ_var, (θ, θ_mean, θ_var, N, K) ->
            begin
                u = 0.5N + shape(θ_var.distr)
                v = 0.5vec(sum(abs2, θ - θ_mean, dims=1)) .+ scale(θ_var.distr)
                [rand(InverseGamma(u, v[k])) for k in 1:K]
            end
        )
    ]

    setsamplers!(model, samplers)
    return model
end

max_stars_repo_name: reich-group/MIMIX

max_issues_repo_name: reich-group/MIMIX

max_forks_repo_name: reich-group/MIMIX


text:
function my_pmap(f, lst)
    np = nprocs()  # determine the number of processes available
    n = length(lst)
    results = Vector{Any}(n)
    i = 1
    # function to produce the next work item from the queue.
    # in this case it's just an index.
    nextidx() = (idx=i; i+=1; idx)
    @sync begin
        for p=1:np
            if p != myid() || np == 1
                @async begin
                    while true
                        idx = nextidx()
                        if idx > n
                            break
                        end
                        println("p=",p)
                        results[idx] = remotecall_fetch(f, p, lst[idx])
                    end
                end
            end
        end
    end
    results
end

@everywhere f(x) = 2x

A = rand(100)
my_pmap(f, A)


text:
################################################################################
# Utils
################################################################################

# rotation matrix rotating unit vector a onto unit vector b
function rot(a, b)
	v = cross(a, b)
	s = sqrt(transpose(v) * v)
	c = transpose(a) * b

	R = Diagonal(@SVector ones(3)) + skew(v) + 1.0 / (1.0 + c) * skew(v) * skew(v)
end

function skew(x)
	SMatrix{3,3}([0.0 -x[3] x[2];
	               x[3] 0.0 -x[1];
				   -x[2] x[1] 0.0])
end

function default_background!(vis)
    setvisible!(vis["/Background"], true)
    setprop!(vis["/Background"], "top_color", RGBA(1.0, 1.0, 1.0, 1.0))
    setprop!(vis["/Background"], "bottom_color", RGBA(1.0, 1.0, 1.0, 1.0))
    setvisible!(vis["/Axes"], false)
end

function get_material(;α=1.0)
    orange_col = [1,153/255,51/255]
    blue_col = [51/255,1,1]
    black_col = [0,0,0]
    orange_mat = MeshPhongMaterial(color=MeshCat.RGBA(orange_col...,α))
    blue_mat = MeshPhongMaterial(color=MeshCat.RGBA(blue_col...,α))
    black_mat = MeshPhongMaterial(color=MeshCat.RGBA(black_col...,α))
    return orange_mat, blue_mat, black_mat
end

function get_line_material(size::Real; α=1.0)
    orange_col = [1,153/255,51/255]
    blue_col = [51/255,1,1]
    black_col = [0,0,0]
    orange_mat = LineBasicMaterial(color=color=RGBA(orange_col...,α), linewidth=size)
    blue_mat = LineBasicMaterial(color=color=RGBA(blue_col...,α), linewidth=size)
    black_mat = LineBasicMaterial(color=color=RGBA(black_col...,α), linewidth=size)
    return orange_mat, blue_mat, black_mat
end

function clean!(vis::Visualizer)
    delete!(vis["/Grid"])
    delete!(vis["/Axes"])
    return nothing
end

function ModifiedMeshFileObject(obj_path::String, material_path::String; scale::T=1.0) where {T}
    obj = MeshFileObject(obj_path)
    rescaled_contents = rescale_contents(obj_path, scale=scale)
    material = select_material(material_path)
    mod_obj = MeshFileObject(
        rescaled_contents,
        obj.format,
        material,
        obj.resources,
        )
    return mod_obj
end

function rescale_contents(obj_path::String; scale::T=1.0) where T
    lines = readlines(obj_path)
    rescaled_lines = copy(lines)
    for (k,line) in enumerate(lines)
        if length(line) >= 2
            if line[1] == 'v'
                stringvec = split(line, " ")
                vals = []
                try
                    vals = map(x->parse(Float64,x),stringvec[end-2:end])
                catch e
                    vals = map(x->parse(Float64,x),stringvec[end-1:end])
                end
                rescaled_vals = vals .* scale
                rescaled_lines[k] = join([stringvec[1]; string.(rescaled_vals)], " ")
            end
        end
    end
    rescaled_contents = join(rescaled_lines, "\r\n")
    return rescaled_contents
end

function select_material(material_path::String)
    mtl_file = open(material_path)
    mtl = read(mtl_file, String)
    return mtl
end

function transRot(model::AbstractGameModel, x::SVx, i::Int) where {SVx}
    return transRot(standardize(model, x[model.pz[i]]))
end

function transRot(x::VehicleState)
    t = Translation(x.x, x.y, 0.)
    r = LinearMap(AngleAxis(x.θ, 0, 0, 1.))
    tr = compose(t,r)
    return compose(t,r)
end


text:
#  MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
#  Copyright (C) 2020-2021, UChicago Argonne, LLC. All rights reserved.
#  Released under the modified BSD license. See COPYING.md for more details.

using Cbc
using JuMP
using MIPLearn

@testset "LearningSolver" begin
    @testset "Model with annotations" begin
        model = build_knapsack_model()
        solver = LearningSolver(Cbc.Optimizer)
        instance = JuMPInstance(model)
        stats = solve!(solver, instance)
        @test stats["mip_lower_bound"] == 11.0
        @test length(instance.samples) == 1
        fit!(solver, [instance])
        solve!(solver, instance)
    end

    @testset "Model without annotations" begin
        model = build_knapsack_model()
        solver = LearningSolver(Cbc.Optimizer)
        instance = JuMPInstance(model)
        stats = solve!(solver, instance)
        @test stats["mip_lower_bound"] == 11.0
    end

    @testset "Save and load" begin
        solver = LearningSolver(Cbc.Optimizer)
        solver.py.components = "Placeholder"
        filename = tempname()
        save(filename, solver)
        @test isfile(filename)
        loaded = load_solver(filename)
        @test loaded.py.components == "Placeholder"
    end

    # @testset "Discard output" begin
    #     instance = build_knapsack_file_instance()
    #     solver = LearningSolver(Cbc.Optimizer)
    #     solve!(solver, instance, discard_output = true)
    #     loaded = load_instance(instance.filename)
    #     @test length(loaded.samples) == 0
    # end
end