text:
__precompile__()

module EZLearn

using OBOParse
using SQLite
using JSON

abstract type ClassifierView end
abstract type Belief end

immutable EZLearnTask
    views::Vector{ClassifierView}
    beliefs::Dict{String, Vector{Belief}}
    models::Dict{String, Vector{Any}}
    ontology::Ontology
    params::Dict
    all_samples::Vector{String}
    cache
end

function ezlearn_step(task::EZLearnTask)
    for v in task.views
        @show typeof(v)
        labels = construct_labels(v, task.beliefs, task)
        belief, model = train_and_predict(v, labels, task)
        push!(task.beliefs[v.id], BeliefDict(belief))
        push!(task.models[v.id], model)
    end
end

train_and_predict(v::ClassifierView, labels, task) = error("not implemented")
construct_labels(v::ClassifierView, beliefs, task) = error("not implemented")

include("beliefs.jl")

export EZLearnTask, BeliefDict, BeliefSQLite, get_beliefs, intersect_labels_core, intersect_labels_simple, ezlearn_step, store_beliefs, get_all_samples

end


text:
const QUINTIC_ZERO_ONE_INTERPOLATOR = fit_quintic(x0=0.0, xf=1.0, y0=0.0, yd0=0.0, ydd0=0.0, yf=1.0, ydf=0.0, yddf=0.0)

function interpolated_orientation_trajectory(t0::Number, tf::Number, rot0::Rotation, rotf::Rotation)
    Interpolated(t0, tf, Quat(rot0), Quat(rotf), QUINTIC_ZERO_ONE_INTERPOLATOR, min_num_derivs=Val(2))
end

function init_support!(end_effector_controller::SE3PDController; t0::Number, tf::Number)
    # TODO: frame lookup is kind of nasty
    bodyframe = end_effector_controller.trajectory[].body
    baseframe = end_effector_controller.trajectory[].base
    end_effector_controller.gains[] = SE3PDGains(FramePDGains(bodyframe, PDGains(0.0, 15.0)), FramePDGains(bodyframe, PDGains(0.0, 0.0)))
    angulartraj = interpolated_orientation_trajectory(t0, tf, one(Quat{Float64}), one(Quat{Float64}))
    lineartraj = convert(BasicFootTrajectory{Float64}, SVector(0.0, 0.0, 0.0), tf)
    end_effector_controller.trajectory[] = SE3Trajectory(bodyframe, baseframe, angulartraj, lineartraj)
end

function init_swing!(end_effector_controller::SE3PDController, pose0::Transform3D, posef::Transform3D;
        t0::Number, tf::Number, Δzmid::Number=0.1, zdf::Number=-0.1)
    # TODO: frame lookup is kind of nasty
    bodyframe = end_effector_controller.trajectory[].body
    baseframe = end_effector_controller.trajectory[].base
    end_effector_controller.gains[] = SE3PDGains(
        FramePDGains(bodyframe, critically_damped_gains(100.)),
        FramePDGains(bodyframe, critically_damped_gains(100.))
    )
    angulartraj = interpolated_orientation_trajectory(t0, tf, rotation(pose0), rotation(posef))
    lineartraj = BasicFootTrajectory(t0, tf, translation(pose0), Δzmid, translation(posef), zdf)
    end_effector_controller.trajectory[] = SE3Trajectory(bodyframe, baseframe, angulartraj, lineartraj)
end

ext: jl


text:
"""
    This struct contains all the required data, parameters and solvers
    for creating and running AgentSimulation.
        base_dir: Directory where Simulation input data for markets and investors are stored.
        sys_dir: Test system directory.
        solver: Solvers used for optimization problems. (The solver should be able to solve QP for price prediction and MILP for SIIP production cost model)
        siip_market_clearing: Whether SIIP production cost model is to be used for energy market clearing. If false, the endogenous Economic Dispatch model will be used for market clearing.
        start_year: Start year for the simulation (default is set to 2020)
        total_horizon: Number of years of data available for the simulation.
        rolling_horizon: Number of years to be used for price prediction. If end of rolling horizon exceeds the years of available data, a receding horizon approach is used.
        simulation_years: Number of years to be simulated.
        num_rep_days: Number of representative days used for price prediction.
        heterogeneity: Whether investors' heterogeneous financial characteristics and technology preferences are modeled.
        forecast_type: "Perfect" or "imperfect" forecasts used for price prediction.
        info_symmetry: Whether investors have symmetric information about forecast parameters.
        belief_update: Whether investors' beliefs are updated each year after actual market clearing.
        uncertainty: Whether multiple probability weighted scenarios are used instead of a deterministic forecast.
        risk_aversion: Whether investors are risk averse.
        parallel_investors: Whether investors' price prediction is to be parallelized.
        parallel_scenarios: Whether each investor's price prediction scenarios are to be parallelized.
"""

struct CaseDefinition
    base_dir::String
    sys_dir::String
    solver::JuMP.MOI.OptimizerWithAttributes
    siip_market_clearing::Bool
    start_year::Int64
    total_horizon::Int64
    rolling_horizon::Int64
    simulation_years::Int64
    num_rep_days::Int64
    heterogeneity::Bool
    forecast_type::String
    info_symmetry::Bool
    belief_update::Bool
    uncertainty::Bool
    risk_aversion::Bool
    parallel_investors::Bool
    parallel_scenarios::Bool

    function CaseDefinition(base_dir,
                            sys_dir,
                            solver,
                            siip_market_clearing,
                            start_year,
                            total_horizon,
                            rolling_horizon,
                            simulation_years,
                            num_rep_days,
                            heterogeneity,
                            forecast_type,
                            info_symmetry,
                            belief_update,
                            uncertainty,
                            risk_aversion,
                            parallel_investors,
                            parallel_scenarios)

        @assert total_horizon >= simulation_years

        forecast_type = lowercase(forecast_type)
        @assert forecast_type == "perfect" || lowercase(forecast_type) == "imperfect"

        if forecast_type == "perfect"
            @assert info_symmetry == true
            @assert belief_update == false
            @assert uncertainty == false
            @assert risk_aversion == false
        end

        return new(base_dir,
                   sys_dir,
                   solver,
                   siip_market_clearing,
                   start_year,
                   total_horizon,
                   rolling_horizon,
                   simulation_years,
                   num_rep_days,
                   heterogeneity,
                   forecast_type,
                   info_symmetry,
                   belief_update,
                   uncertainty,
                   risk_aversion,
                   parallel_investors,
                   parallel_scenarios)
    end
end

function CaseDefinition(base_dir::String,
                        sys_dir::String,
                        solver::JuMP.MOI.OptimizerWithAttributes;
                        siip_market_clearing::Bool = true,
                        start_year::Int64 = 2020,
                        total_horizon::Int64 = 20,
                        rolling_horizon::Int64 = 10,
                        simulation_years::Int64 = 10,
                        num_rep_days::Int64 = 12,
                        heterogeneity::Bool = false,
                        forecast_type::String = "perfect",
                        info_symmetry::Bool = true,
                        belief_update::Bool = false,
                        uncertainty::Bool = false,
                        risk_aversion::Bool = false,
                        parallel_investors::Bool = false,
                        parallel_scenarios::Bool = false)

    CaseDefinition(base_dir,
                   sys_dir,
                   solver,
                   siip_market_clearing,
                   start_year,
                   total_horizon,
                   rolling_horizon,
                   simulation_years,
                   num_rep_days,
                   heterogeneity,
                   forecast_type,
                   info_symmetry,
                   belief_update,
                   uncertainty,
                   risk_aversion,
                   parallel_investors,
                   parallel_scenarios)
end

get_base_dir(case::CaseDefinition) = case.base_dir
get_sys_dir(case::CaseDefinition) = case.sys_dir
get_solver(case::CaseDefinition) = case.solver
get_siip_market_clearing(case::CaseDefinition) = case.siip_market_clearing
get_start_year(case::CaseDefinition) = case.start_year
get_total_horizon(case::CaseDefinition) = case.total_horizon
get_rolling_horizon(case::CaseDefinition) = case.rolling_horizon
get_simulation_years(case::CaseDefinition) = case.simulation_years
get_num_rep_days(case::CaseDefinition) = case.num_rep_days
get_heterogeneity(case::CaseDefinition) = case.heterogeneity
get_info_symmetry(case::CaseDefinition) = case.info_symmetry
get_belief_update(case::CaseDefinition) = case.belief_update
get_forecast_type(case::CaseDefinition) = case.forecast_type
get_uncertainty(case::CaseDefinition) = case.uncertainty
get_risk_aversion(case::CaseDefinition) = case.risk_aversion
get_parallel_investors(case::CaseDefinition) = case.parallel_investors
get_parallel_scenarios(case::CaseDefinition) = case.parallel_scenarios

function get_name(case::CaseDefinition)
    if get_heterogeneity(case)
        investors = "Heterogeneous"
    else
        investors = "Homogeneous"
    end

    if get_info_symmetry(case)
        information = "InfoSym"
    else
        information = "InfoASym"
    end

    if get_belief_update(case)
        update = "UpdateBelief"
    else
        update = "NoUpdate"
    end

    if get_uncertainty(case)
        uncertainty = "Uncertain"
    else
        uncertainty = "Deterministic"
    end

    if get_risk_aversion(case)
        risk = "RiskAverse"
    else
        risk = "RiskNeutral"
    end

    case_name = "$(investors)_$(information)_Forecast-$(get_forecast_type(case))_$(uncertainty)_$(update)_$(risk)_$(get_simulation_years(case))years"

    return case_name
end

function get_data_dir(case::CaseDefinition)
    base_dir = get_base_dir(case)
    case_dir = joinpath(base_dir, get_name(case))

    return case_dir
end


text:
# ps4000openunit.jl

export ps4000openunit

function ps4000openunit()
  handle = Ref{UInt16}()
  ps_status = ccall((:ps4000OpenUnit, ps4000driver),
                  PicoStatus,
                  (Ref{UInt16},),
                  handle)
  return handle[]
end


text:
function checkout_branch(dir::AbstractString,
                         branch::AbstractString;
                         git_command::AbstractString = "git")
    original_working_directory = pwd()
    cd(dir)
    Base.run(`$(git_command) checkout $(branch)`)
    cd(original_working_directory)
end

clone_repo(repo::GitHub.Repo) = clone_repo(repo_url(repo))

function clone_repo(url::AbstractString)
    parent_dir = mktempdir()
    atexit(() -> rm(parent_dir; force = true, recursive = true))
    repo_dir = joinpath(parent_dir, "REPO")
    my_retry(() -> _clone_repo_into_dir(url, repo_dir))
    @info("Clone was successful")
    return repo_dir
end

function _clone_repo_into_dir(url::AbstractString, repo_dir)
    @info("Attempting to clone...")
    rm(repo_dir; force = true, recursive = true)
    mkpath(repo_dir)
    LibGit2.clone(url, repo_dir)
    return repo_dir
end

function _comment_disclaimer()
    result = string("\n\n",
                    "Note that the guidelines are only required for the pull request ",
                    "to be merged automatically. However, it is **strongly recommended** ",
                    "to follow them, since otherwise the pull request needs to be ",
                    "manually reviewed and merged by a human.")
    return result
end

function _comment_noblock()
    result = string("\n\n---\n",
                    "If you want to prevent this pull request from ",
                    "being auto-merged, simply leave a comment. ",
                    "If you want to post a comment without blocking ",
                    "auto-merging, you must include the text ",
                    "`[noblock]` in your comment.")
    return result
end

function comment_text_pass(::NewVersion,
                           suggest_onepointzero::Bool,
                           version::VersionNumber)
    result = string("Your `new version` pull request met all of the ",
                    "guidelines for auto-merging and is scheduled to ",
                    "be merged in the next round.",
                    _comment_noblock(),
                    _onepointzero_suggestion(suggest_onepointzero, version),
                    "\n<!-- [noblock] -->")
    return result
end

function comment_text_pass(::NewPackage,
                           suggest_onepointzero::Bool,
                           version::VersionNumber)
    result = string("Your `new package` pull request met all of the ",
                    "guidelines for auto-merging and is scheduled to ",
                    "be merged when the mandatory waiting period has elapsed.",
                    _comment_noblock(),
                    _onepointzero_suggestion(suggest_onepointzero, version),
                    "\n<!-- [noblock] -->")
    return result
end

function comment_text_fail(::NewPackage,
                           reasons::Vector{String},
                           suggest_onepointzero::Bool,
                           version::VersionNumber)
    reasons_formatted = join(string.("- ", reasons), "\n")
    result = string("Your `new package` pull request does not meet ",
                    "the following guidelines for auto-merging:\n\n",
                    reasons_formatted,
                    _comment_disclaimer(),
                    _comment_noblock(),
                    _onepointzero_suggestion(suggest_onepointzero, version),
                    "\n<!-- [noblock] -->")
    return result
end

function comment_text_fail(::NewVersion,
                           reasons::Vector{String},
                           suggest_onepointzero::Bool,
                           version::VersionNumber)
    reasons_formatted = join(string.("- ", reasons), "\n")
    result = string("Your `new version` pull request does not meet ",
                    "the following guidelines for auto-merging:\n\n",
                    reasons_formatted,
                    _comment_disclaimer(),
                    _comment_noblock(),
                    _onepointzero_suggestion(suggest_onepointzero, version),
                    "\n<!-- [noblock] -->")
    return result
end

function comment_text_merge_now()
    result = string("The mandatory waiting period has elapsed.\n\n",
                    "Your pull request is ready to merge.\n\n",
                    "I will now merge this pull request.",
                    "\n<!-- [noblock] -->")
    return result
end

is_julia_stdlib(name) = name in julia_stdlib_list()

function julia_stdlib_list()
    return readdir(Pkg.Types.stdlib_dir())
end

function now_utc()
    # my_local_timezone = TimeZones.localzone()
    utc_timezone = TimeZones.TimeZone("UTC")
    _now_utc = Dates.now(utc_timezone)
    return _now_utc
end

function _onepointzero_suggestion(suggest_onepointzero::Bool,
                                  version::VersionNumber)
    if suggest_onepointzero && version < v"1.0.0"
        result = string("\n\n---\n",
                        "On a separate note, I see that you are registering ",
                        "a release with a version number of the form ",
                        "`v0.X.Y`.\n\n",
                        "Does your package have a stable public API? ",
                        "If so, then it's time for you to register version ",
                        "`v1.0.0` of your package. ",
                        "(This is not a requirement. ",
                        "It's just a recommendation.)\n\n",
                        "If your package does not yet have a stable public ",
                        "API, then of course you are not yet ready to ",
                        "release version `v1.0.0`.")
        return result
    else
        return ""
    end
end

function time_is_already_in_utc(dt::Dates.DateTime)
    # my_local_timezone = TimeZones.localzone()
    utc_timezone = TimeZones.TimeZone("UTC")
    result = TimeZones.ZonedDateTime(dt, utc_timezone; from_utc = true)
    return result
end


text:
@testset "Neighbour joining implementation" begin

    D = DataFrame(CSV.File(joinpath(@__DIR__,"..","examples","caudata_dist.txt")); copycols=false)
    tree = nj(D)
    # read results from ape implementation of nj
    apetree = readTopology(joinpath(@__DIR__, "..", "examples", "caudata_dist_nj.txt"))
    @test hardwiredClusterDistance(tree, apetree, false) == 0
    # also check branch lengths (more or less)
    @test sort!([e.length for e in tree.edge]) ≈ sort!([e.length for e in apetree.edge])

    # example where Ints are converted to Floats, and there's a < 0 edge length
    df = DataFrame(s1=[0, 5,9,9,5], s2=[5,0,10,10,6], s3=[9,10,0,8,4],
                   s4=[9,10,8,0,0], s5=[5,6, 4, 0,0])
    tree = (@test_logs (:info, r"have negative lengths") nj(df))
    @test writeTopology(tree) == "(((s1:2.0,s2:3.0):3.0,s3:4.0):2.0,s4:2.0,s5:-2.0);"
    # example with no names argument, force_nonnegative_edges
    D = [0 5 9 9 6.5; 5 0 10 10 7.5; 9 10 0 8 5.5; 9 10 8 0 1.5; 6.5 7.5 5.5 1.5 0]
    tree = (@test_logs (:info, r"reset to 0") PhyloNetworks.nj!(D, force_nonnegative_edges=true))
    @test writeTopology(tree) == "(((1:2.0,2:3.0):3.0,3:4.0):2.0,4:2.0,5:0.0);"
end


text:
function readMap(fileName)
    result = Dict()
    i, j = 0, 0

    for line in eachline(fileName)
        j = 0
        for c in line
            result[(i, j)] = c
            j += 1
        end
        i += 1
    end

    result
end

function nextState(stateMap, i, j)
    current = stateMap[(i, j)]

    if (current == '.')
        return '.'
    end

    neighbours = countNeighbours(stateMap, i, j)

    if (neighbours == 0 && current == 'L')
        return '#'
    end

    if (current == '#' && neighbours >= 5)
        return 'L'
    end

    return current
end

function countNeighbours(stateMap, i, j)
    result = 0

    for (dx, dy) in [(1, 1), (1, 0), (1, -1), (0, 1), (0, -1), (-1, 1), (-1, 0), (-1, -1)]
        ii, jj = i, j
        while (true)
            val = get(stateMap, (ii + dx, jj + dy), 'x')

            ii += dx
            jj += dy

            if (val == '#')
                result += 1
                break
            end

            if (val == 'x' || val == 'L')
                break
            end
        end
    end

    result
end

function solve(fileName)
    seatMap = readMap(fileName)
    while (true)
        newState = Dict()

        stateChanges = 0

        for ((i, j), value) in pairs(seatMap)
            newState[(i, j)] = nextState(seatMap, i, j)
            if (seatMap[(i, j)] != newState[(i, j)])
                stateChanges += 1
            end
        end

        seatMap = newState

        if (stateChanges == 0)
            break;
        end
    end

    count(c->(c == '#'), values(seatMap))
end

const fileName = "input.txt"
println(solve(fileName))


text:
module MadNLPLapackGPU

import ..MadNLPGPU:
    @kwdef, Logger, @debug, @warn, @error,
    AbstractOptions, AbstractLinearSolver, set_options!,
    SymbolicException,FactorizationException,SolveException,InertiaException,
    introduce, factorize!, solve!, improve!, is_inertia, inertia, MadNLPLapackCPU, tril_to_full!,
    CUBLAS, CUSOLVER, CuVector, CuMatrix, version, R_64F
import .CUSOLVER:
    cusolverDnDsytrf_bufferSize, cusolverDnDsytrf,
    cusolverDnDpotrf_bufferSize, cusolverDnDpotrf, cusolverDnDpotrs,
    cusolverDnDgetrf_bufferSize, cusolverDnDgetrf, cusolverDnDgetrs,
    cusolverDnDgeqrf_bufferSize, cusolverDnDgeqrf, cusolverDnDgeqrf_bufferSize,
    cusolverDnDormqr_bufferSize, cusolverDnDormqr,
    libcusolver, cusolverStatus_t, CuPtr, cudaDataType, cublasFillMode_t, cusolverDnHandle_t, dense_handle
import .CUBLAS: cublasDtrsm_v2, handle, CUBLAS_DIAG_NON_UNIT,
    CUBLAS_FILL_MODE_LOWER, CUBLAS_FILL_MODE_UPPER, CUBLAS_SIDE_LEFT, CUBLAS_OP_N, CUBLAS_OP_T

const INPUT_MATRIX_TYPE = :dense

@enum(Algorithms::Int, BUNCHKAUFMAN = 1, LU = 2, QR = 3, CHOLESKY = 4)
@kwdef mutable struct Options <: AbstractOptions
    lapackgpu_algorithm::Algorithms = BUNCHKAUFMAN
end

mutable struct Solver{MT} <: AbstractLinearSolver
    dense::MT
    fact::CuMatrix{Float64}
    rhs::CuVector{Float64}
    work::CuVector{Float64}
    lwork
    work_host::Vector{Float64}
    lwork_host
    info::CuVector{Int32}
    etc::Dict{Symbol,Any} # throw some algorithm-specific things here
    opt::Options
    logger::Logger
end

function Solver(dense::MT;
                option_dict::Dict{Symbol,Any}=Dict{Symbol,Any}(),
                opt=Options(),logger=Logger(),
                kwargs...) where {MT <: AbstractMatrix}

    set_options!(opt,option_dict,kwargs...)
    fact = CuMatrix{Float64}(undef,size(dense))
    rhs = CuVector{Float64}(undef,size(dense,1))
    work  = CuVector{Float64}(undef, 1)
    lwork = Int32[1]
    work_host  = Vector{Float64}(undef, 1)
    lwork_host = Int32[1]
    info = CuVector{Int32}(undef,1)
    etc = Dict{Symbol,Any}()

    return Solver(dense,fact,rhs,work,lwork,work_host,lwork_host,info,etc,opt,logger)
end

function factorize!(M::Solver)
    if M.opt.lapackgpu_algorithm == BUNCHKAUFMAN
        factorize_bunchkaufman!(M)
    elseif M.opt.lapackgpu_algorithm == LU
        factorize_lu!(M)
    elseif M.opt.lapackgpu_algorithm == QR
        factorize_qr!(M)
    elseif M.opt.lapackgpu_algorithm == CHOLESKY
        factorize_cholesky!(M)
    else
        error(LOGGER,"Invalid lapackgpu_algorithm")
    end
end
function solve!(M::Solver,x)
    if M.opt.lapackgpu_algorithm == BUNCHKAUFMAN
        solve_bunchkaufman!(M,x)
    elseif M.opt.lapackgpu_algorithm == LU
        solve_lu!(M,x)
    elseif M.opt.lapackgpu_algorithm == QR
        solve_qr!(M,x)
    elseif M.opt.lapackgpu_algorithm == CHOLESKY
        solve_cholesky!(M,x)
    else
        error(LOGGER,"Invalid lapackgpu_algorithm")
    end
end

improve!(M::Solver) = false
introduce(M::Solver) = "Lapack-GPU ($(M.opt.lapackgpu_algorithm))"

if version() >= v"11.3.1"

    is_inertia(M::Solver) = M.opt.lapackgpu_algorithm == CHOLESKY  # TODO: implement inertia(M::Solver) for BUNCHKAUFMAN

    function factorize_bunchkaufman!(M::Solver)
        haskey(M.etc,:ipiv) || (M.etc[:ipiv] = CuVector{Int32}(undef,size(M.dense,1)))
        haskey(M.etc,:ipiv64) || (M.etc[:ipiv64] = CuVector{Int64}(undef,length(M.etc[:ipiv])))

        copyto!(M.fact,M.dense)
        cusolverDnDsytrf_bufferSize(
            dense_handle(),Int32(size(M.fact,1)),M.fact,Int32(size(M.fact,2)),M.lwork)
        length(M.work) < M.lwork[] && resize!(M.work,Int(M.lwork[]))
        cusolverDnDsytrf(
            dense_handle(),CUBLAS_FILL_MODE_LOWER,
            Int32(size(M.fact,1)),M.fact,Int32(size(M.fact,2)),
            M.etc[:ipiv],M.work,M.lwork[],M.info)
        return M
    end

    function solve_bunchkaufman!(M::Solver,x)

        copyto!(M.etc[:ipiv64],M.etc[:ipiv])
        copyto!(M.rhs,x)
        ccall((:cusolverDnXsytrs_bufferSize, libcusolver()), cusolverStatus_t,
              (cusolverDnHandle_t, cublasFillMode_t, Int64, Int64, cudaDataType,
               CuPtr{Cdouble}, Int64, CuPtr{Int64}, cudaDataType,
               CuPtr{Cdouble}, Int64, Ptr{Int64}, Ptr{Int64}),
              dense_handle(), CUBLAS_FILL_MODE_LOWER,
              size(M.fact,1),1,R_64F,M.fact,size(M.fact,2),
              M.etc[:ipiv64],R_64F,M.rhs,length(M.rhs),M.lwork,M.lwork_host)
        length(M.work) < M.lwork[] && resize!(M.work,Int(M.lwork[]))
        length(M.work_host) < M.lwork_host[] && resize!(work_host,Int(M.lwork_host[]))
        ccall((:cusolverDnXsytrs, libcusolver()), cusolverStatus_t,
              (cusolverDnHandle_t, cublasFillMode_t, Int64, Int64, cudaDataType,
               CuPtr{Cdouble}, Int64, CuPtr{Int64}, cudaDataType,
               CuPtr{Cdouble}, Int64, CuPtr{Cdouble}, Int64, Ptr{Cdouble}, Int64,
               CuPtr{Int64}),
              dense_handle(),CUBLAS_FILL_MODE_LOWER,
              size(M.fact,1),1,R_64F,M.fact,size(M.fact,2),
              M.etc[:ipiv64],R_64F,M.rhs,length(M.rhs),M.work,M.lwork[],M.work_host,M.lwork_host[],M.info)
        copyto!(x,M.rhs)

        return x
    end
else
    is_inertia(M::Solver) =
        M.opt.lapackgpu_algorithm == CHOLESKY || M.opt.lapackgpu_algorithm == CHOLESKY

    function factorize_bunchkaufman!(M::Solver)
        haskey(M.etc,:ipiv) || (M.etc[:ipiv] = CuVector{Int32}(undef,size(M.dense,1)))

        copyto!(M.fact,M.dense)
        CUSOLVER.cusolverDnDsytrf_bufferSize(
            CUSOLVER.dense_handle(),Int32(size(M.fact,1)),M.fact,Int32(size(M.fact,2)),M.lwork)
        length(M.work) < M.lwork[] && resize!(M.work,Int(M.lwork[]))
        CUSOLVER.cusolverDnDsytrf(
            CUSOLVER.dense_handle(),CUBLAS.CUBLAS_FILL_MODE_LOWER,
            Int32(size(M.fact,1)),M.fact,Int32(size(M.fact,2)),
            M.etc[:ipiv],M.work,M.lwork[],M.info)

        # need to send the factorization back to cpu to call mkl sytrs --------------
        haskey(M.etc,:fact_cpu) || (M.etc[:fact_cpu] = Matrix{Float64}(undef,size(M.dense)))
        haskey(M.etc,:ipiv_cpu) || (M.etc[:ipiv_cpu] = Vector{Int}(undef,length(M.etc[:ipiv])))
        haskey(M.etc,:info_cpu) || (M.etc[:info_cpu] = Vector{Int}(undef,length(M.info)))
        copyto!(M.etc[:fact_cpu],M.fact)
        copyto!(M.etc[:ipiv_cpu],M.etc[:ipiv])
        copyto!(M.etc[:info_cpu],M.info)
        # ---------------------------------------------------------------------------
        return M
    end

    function solve_bunchkaufman!(M::Solver,x)
        ccall(
            (:dsytrs_64_,"libopenblas64_"), # MKL doesn't work for some reason...
            Cvoid,
            (Ref{Cchar},Ref{Int},Ref{Int},Ptr{Cdouble},Ref{Int},Ptr{Int},Ptr{Cdouble},Ref{Int},Ptr{Int}),
            'L',size(M.fact,1),1,M.etc[:fact_cpu],size(M.fact,2),M.etc[:ipiv_cpu],x,length(x),[1])

        return x
    end
end

function factorize_lu!(M::Solver)
    haskey(M.etc,:ipiv) || (M.etc[:ipiv] = CuVector{Int32}(undef,size(M.dense,1)))
    tril_to_full!(M.dense)
    copyto!(M.fact,M.dense)
    cusolverDnDgetrf_bufferSize(
        dense_handle(),Int32(size(M.fact,1)),Int32(size(M.fact,2)),
        M.fact,Int32(size(M.fact,2)),M.lwork)
    length(M.work) < M.lwork[] && resize!(M.work,Int(M.lwork[]))
    cusolverDnDgetrf(
        dense_handle(),Int32(size(M.fact,1)),Int32(size(M.fact,2)),
        M.fact,Int32(size(M.fact,2)),M.work,M.etc[:ipiv],M.info)
    return M
end

function solve_lu!(M::Solver,x)
    copyto!(M.rhs,x)
    cusolverDnDgetrs(
        dense_handle(),CUBLAS_OP_N,
        Int32(size(M.fact,1)),Int32(1),M.fact,Int32(size(M.fact,2)),
        M.etc[:ipiv],M.rhs,Int32(length(M.rhs)),M.info)
    copyto!(x,M.rhs)
    return x
end

function factorize_qr!(M::Solver)
    haskey(M.etc,:tau) || (M.etc[:tau] = CuVector{Float64}(undef,size(M.dense,1)))
    haskey(M.etc,:one) || (M.etc[:one] = ones(1))
    tril_to_full!(M.dense)
    copyto!(M.fact,M.dense)
    cusolverDnDgeqrf_bufferSize(dense_handle(),Int32(size(M.fact,1)),Int32(size(M.fact,2)),M.fact,Int32(size(M.fact,2)),M.lwork)
    length(M.work) < M.lwork[] && resize!(M.work,Int(M.lwork[]))
    cusolverDnDgeqrf(dense_handle(),Int32(size(M.fact,1)),Int32(size(M.fact,2)),M.fact,Int32(size(M.fact,2)),M.etc[:tau],M.work,M.lwork[],M.info)
    return M
end

function solve_qr!(M::Solver,x)
    copyto!(M.rhs,x)
    cusolverDnDormqr_bufferSize(dense_handle(),CUBLAS_SIDE_LEFT,CUBLAS_OP_T,
                                Int32(size(M.fact,1)),Int32(1),Int32(length(M.etc[:tau])),M.fact,Int32(size(M.fact,2)),M.etc[:tau],M.rhs,Int32(length(M.rhs)),M.lwork)
    length(M.work) < M.lwork[] && resize!(M.work,Int(M.lwork[]))
    cusolverDnDormqr(dense_handle(),CUBLAS_SIDE_LEFT,CUBLAS_OP_T,
                     Int32(size(M.fact,1)),Int32(1),Int32(length(M.etc[:tau])),M.fact,Int32(size(M.fact,2)),M.etc[:tau],M.rhs,Int32(length(M.rhs)),M.work,M.lwork[],M.info)
    cublasDtrsm_v2(handle(),CUBLAS_SIDE_LEFT,CUBLAS_FILL_MODE_UPPER,CUBLAS_OP_N,CUBLAS_DIAG_NON_UNIT,
                   Int32(size(M.fact,1)),Int32(1),M.etc[:one],M.fact,Int32(size(M.fact,2)),M.rhs,Int32(length(M.rhs)))
    copyto!(x,M.rhs)
    return x
end

function factorize_cholesky!(M::Solver)
    copyto!(M.fact,M.dense)
    cusolverDnDpotrf_bufferSize(
        dense_handle(),CUBLAS_FILL_MODE_LOWER,
        Int32(size(M.fact,1)),M.fact,Int32(size(M.fact,2)),M.lwork)
    length(M.work) < M.lwork[] && resize!(M.work,Int(M.lwork[]))
    cusolverDnDpotrf(
        dense_handle(),CUBLAS_FILL_MODE_LOWER,
        Int32(size(M.fact,1)),M.fact,Int32(size(M.fact,2)),
        M.work,M.lwork[],M.info)
    return M
end

function solve_cholesky!(M::Solver,x)
    copyto!(M.rhs,x)
    cusolverDnDpotrs(
        dense_handle(),CUBLAS_FILL_MODE_LOWER,
        Int32(size(M.fact,1)),Int32(1),M.fact,Int32(size(M.fact,2)),
        M.rhs,Int32(length(M.rhs)),M.info)
    copyto!(x,M.rhs)
    return x
end

function inertia(M::Solver)
    if M.opt.lapackgpu_algorithm == BUNCHKAUFMAN
        inertia(M.etc[:fact_cpu],M.etc[:ipiv_cpu],M.etc[:info_cpu][])
    elseif M.opt.lapackgpu_algorithm == CHOLESKY
        sum(M.info) == 0 ? (size(M.fact,1),0,0) : (0,size(M.fact,1),0)
    else
        error(LOGGER,"Invalid lapackcpu_algorithm")
    end
end

end # module


text:
const SYSTEM_UUIDS = Dict{Int,UUID}()

function system_uuid_fallback()
    get!(SYSTEM_UUIDS, myid()) do
        username = get(ENV, "USER", "__global")
        uuid_file = joinpath(tempdir(), "dagger-system-uuid-$username")
        if !isfile(uuid_file)
            temp_uuid_file, temp_io = mktemp(; cleanup=false)
            uuid = uuid4()
            write(temp_io, string(uuid))
            flush(temp_io)
            close(temp_io)
            try
                # Try to make this the UUID file
                mv(temp_uuid_file, uuid_file; force=false)
                return uuid
            catch err
                err isa ArgumentError || rethrow(err)
                # Failed, clean up temp file, and read existing UUID file
                rm(temp_uuid_file)
            end
        end
        return parse(UUID, read(uuid_file, String))
    end
    SYSTEM_UUIDS[myid()]
end

if Sys.islinux()
    system_uuid() = system_uuid_fallback()
    # TODO: dmidecode
else
    system_uuid() = system_uuid_fallback()
end

"Gets the unique UUID for the server hosting worker `wid`."
function system_uuid(wid::Integer)
    get!(SYSTEM_UUIDS, wid) do
        remotecall_fetch(system_uuid, wid)
    end
end


text:
using Cropbox
using SimpleCrop
using Test

using CSV
using DataFrames
using TimeZones

loaddata(f) = CSV.File(joinpath(@__DIR__, "data/simplecrop", f)) |> DataFrame

config = @config (
    :Clock => :step => 1u"d",
    :Calendar => :init => ZonedDateTime(1987, 1, 1, tz"UTC"),
    :Weather => :weather_data => loaddata("weather.csv"),
    :SoilWater => :irrigation_data => loaddata("irrigation.csv"),
)

@testset "simplecrop" begin
    r = simulate(SimpleCrop.Model; config, stop = :endsim)
    visualize(r, :DATE, :LAI; kind = :line) |> display
    visualize(r, :DATE, :(SWC/DP); yunit = u"mm^3/mm^3", kind = :line) |> display
end


text:
struct CG{T <: AbstractFloat}
    u::Vector{T}
    r::Vector{T}
    c::Vector{T}
    function CG{T}(dim_b::Integer, dim_x::Integer) where {T <: AbstractFloat}
        new{T}(zeros(T, dim_b), zeros(T, dim_b), zeros(T, dim_x))
    end
end
CG{T}(dims_A::Dims{2}) where {T} = CG{T}(dims_A...)

function (cg::CG{T})(x, A, b; tol::Real = sqrt(eps(real(eltype(b)))), maxiter::Int = size(A, 2), initiallyzero::Bool = false) where {T}
    if !(T == eltype(x) == eltype(A) == eltype(b))
        @warn "eltype mismatch: expected x, A, and b to have eltype $T. CG may run slower as a result" maxlog=1
        error()
    end

    d1, d2 = size(A)
    u = cg.u
    r = cg.r
    c = cg.c

    if !(length(r) == d1 && length(u) == length(c) == d2)
        resize!(u, d2)
        resize!(r, d1)
        resize!(c, d2)
    end

    fill!(u, zero(T))
    copyto!(r, b)
    # do not need to zero out c as e overwrite it below

    if initiallyzero
        # ignore x and treat x0 = 0 --> r0 = b
        residual = norm(b)
        reltol = residual * tol
    else
        # x0 = x --> r0 = b - A*x0
        mul!(c, A, x) # A*x0
        r .-= c # b - A*x0
        residual = norm(r)
        reltol = norm(b) * tol # TODO norm(b) or norm(r)?
    end

    prev_residual = one(residual)
    for _ = 1:maxiter
        residual <= reltol && break

        β = residual^2 / prev_residual^2
        u .= r .+ β .* u

        mul!(c, A, u)
        α = residual^2 / dot(u, c)

        # Improve solution and residual
        x .+= α .* u
        r .-= α .* c

        prev_residual = residual
        residual = norm(r)
    end
    x
end


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

@doc raw"""
    ObjsMessageAttachments(;
        fallback=nothing,
        id=nothing,
        image_bytes=nothing,
        image_height=nothing,
        image_url=nothing,
        image_width=nothing,
    )

    - fallback::String
    - id::Int32
    - image_bytes::Int32
    - image_height::Int32
    - image_url::String
    - image_width::Int32
"""
mutable struct ObjsMessageAttachments <: SwaggerModel
    fallback::Any # spec type: Union{ Nothing, String } # spec name: fallback
    id::Any # spec type: Union{ Nothing, Int32 } # spec name: id
    image_bytes::Any # spec type: Union{ Nothing, Int32 } # spec name: image_bytes
    image_height::Any # spec type: Union{ Nothing, Int32 } # spec name: image_height
    image_url::Any # spec type: Union{ Nothing, String } # spec name: image_url
    image_width::Any # spec type: Union{ Nothing, Int32 } # spec name: image_width

    function ObjsMessageAttachments(;fallback=nothing, id=nothing, image_bytes=nothing, image_height=nothing, image_url=nothing, image_width=nothing)
        o = new()
        validate_property(ObjsMessageAttachments, Symbol("fallback"), fallback)
        setfield!(o, Symbol("fallback"), fallback)
        validate_property(ObjsMessageAttachments, Symbol("id"), id)
        setfield!(o, Symbol("id"), id)
        validate_property(ObjsMessageAttachments, Symbol("image_bytes"), image_bytes)
        setfield!(o, Symbol("image_bytes"), image_bytes)
        validate_property(ObjsMessageAttachments, Symbol("image_height"), image_height)
        setfield!(o, Symbol("image_height"), image_height)
        validate_property(ObjsMessageAttachments, Symbol("image_url"), image_url)
        setfield!(o, Symbol("image_url"), image_url)
        validate_property(ObjsMessageAttachments, Symbol("image_width"), image_width)
        setfield!(o, Symbol("image_width"), image_width)
        o
    end
end # type ObjsMessageAttachments

const _property_map_ObjsMessageAttachments = Dict{Symbol,Symbol}(Symbol("fallback")=>Symbol("fallback"), Symbol("id")=>Symbol("id"), Symbol("image_bytes")=>Symbol("image_bytes"), Symbol("image_height")=>Symbol("image_height"), Symbol("image_url")=>Symbol("image_url"), Symbol("image_width")=>Symbol("image_width"))
const _property_types_ObjsMessageAttachments = Dict{Symbol,String}(Symbol("fallback")=>"String", Symbol("id")=>"Int32", Symbol("image_bytes")=>"Int32", Symbol("image_height")=>"Int32", Symbol("image_url")=>"String", Symbol("image_width")=>"Int32")
Base.propertynames(::Type{ ObjsMessageAttachments }) = collect(keys(_property_map_ObjsMessageAttachments))
Swagger.property_type(::Type{ ObjsMessageAttachments }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ObjsMessageAttachments[name]))}
Swagger.field_name(::Type{ ObjsMessageAttachments }, property_name::Symbol) =  _property_map_ObjsMessageAttachments[property_name]

function check_required(o::ObjsMessageAttachments)
    (getproperty(o, Symbol("id")) === nothing) && (return false)
    true
end

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


text:
using LinearAlgebra

"""
Symmetric Matrix which only stores upper triangle

data::Vector{T}
"""
struct SymDenseMat{T} <: AbstractMatrix{T}
    data::Vector{T}
    N::Int
end

"""
    function vec_tril(M::AbstractMatrix{T}) where T

Export Symmetri
"""
function SymDenseMat{T}(M::AbstractMatrix{T}) where T
    m, n = size(M)
    m == n || throw(error("not square"))
    nv = n*(n+1) ÷ 2
    v = SymDenseMat{T}(Vector{T}(undef,nv), n)
    k = 0
    for i in 1:n
        for j in i:n
            k += 1
            v.data[k] = M[j, i]
        end
    end
    return v
end

function Base.Matrix(A::SymDenseMat{T}) where {T}
    M = Matrix{T}(undef,A.N,A.N)
    k = 1
    for i in 1:A.N
        for j in i:A.N
            M[i,j] = A.data[k]
            M[j,i] = A.data[k]
            k += 1
        end
    end
    return M
end
Base.size(A::SymDenseMat{T}) where {T} = return (A.N, A.N)
function Base.getindex(A::SymDenseMat{T}, i::Integer, j::Integer) where {T}
    #i>=j || error("i<j")
    if i<j
        return A.data[A.N*(A.N+1)÷2 - (A.N-i+1)*(A.N-i+2)÷2 + j - i + 1]
    else
        return A.data[A.N*(A.N+1)÷2 - (A.N-j+1)*(A.N-j+2)÷2 + i - j + 1]
    end
end
function Base.display(A::SymDenseMat{T}) where {T}
    display(typeof(A))
    display(("N=",A.N))
    display(A.data)
end
function Base.:(*)(A::SymDenseMat{T}, v::AbstractVector{T}) where {T}
    A.N == length(v)  || throw(DimensionMismatch)
    out = zeros(T,A.N)
    for i in 1:A.N
        for j in 1:A.N
            out[i] += A[i,j]*v[j]
        end
    end
    return out
end

max_stars_repo_name: nmayhall-vt/FermiCG

max_issues_repo_name: nmayhall-vt/FermiCG

max_forks_repo_name: nmayhall-vt/FermiCG


text:
#****************************************************************************
#                     Molecular Dynamics Potentials (MDP)
#                            CESMIX-MIT Project
#
#  Contributing authors: Ngoc-Cuong Nguyen (cuongng@mit.edu, exapde@gmail.com)
#****************************************************************************

function writelattice(lat,filename)

    tmp = [lat.style size(lat.basis,1) lat.spaceflag lat.scale];
    basis = transpose(lat.basis);

    nsize = zeros(11,1);
    nsize[1] = length(tmp);
    nsize[2] = length(lat.origin[:]);
    nsize[3] = length(lat.orientx[:]);
    nsize[4] = length(lat.orienty[:]);
    nsize[5] = length(lat.orientz[:]);
    nsize[6] = length(lat.spacing[:]);
    nsize[7] = length(lat.a1[:]);
    nsize[8] = length(lat.a2[:]);
    nsize[9] = length(lat.a3[:]);
    nsize[10] = length(basis[:]);
    nsize[11] = length(lat.type[:]);

    fileID = open(filename,"w");
    write(fileID,Float64(length(nsize[:])));
    write(fileID,Float64.(nsize[:]));
    write(fileID,Float64.(tmp[:]));
    write(fileID,Float64.(lat.origin[:]));
    write(fileID,Float64.(lat.orientx[:]));
    write(fileID,Float64.(lat.orienty[:]));
    write(fileID,Float64.(lat.orientz[:]));
    write(fileID,Float64.(lat.spacing[:]));
    write(fileID,Float64.(lat.a1[:]));
    write(fileID,Float64.(lat.a2[:]));
    write(fileID,Float64.(lat.a3[:]));
    write(fileID,Float64.(basis[:]));
    write(fileID,Float64.(lat.type[:]));
    close(fileID);

end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
using MultiThreadedLU
using Test, Random

n = 1024
a = rand(n,n);
b = rand(n);

@testset "hpl_seq" begin
    x = hpl_seq(a, b)
    @test a*x ≈ b
end

@testset "hpl_mt" begin
    x = hpl_mt(a, b)
    @test a*x ≈ b
end


text:
"""
    skeleton(m::AbstractModel)

Return the model skeleton.
"""
skeleton(m::AbstractModel)

"""
    update_atoms!(m::AbstractModel)
Apply 1 Gibbs update to the cluster atoms.
"""
update_atoms!(m::AbstractModel)

"""
    kernel_pdf(m::AbstractModel, y0, x0, j::Int)

Return the kernel pdf at `y0`, given a predictor `x0` a cluster label `j`.
"""
kernel_pdf(m::AbstractModel, y0, x0, j::Int)
# ok: 2022-03-22


text:
using BenchmarkTools, Random
using Mill, Flux

using Plots, StatsPlots, LaTeXStrings; pgfplotsx()

experiments = [
    (
     label="DenseMatrix * (Full) OneHot",
     A_init = nc -> rand(Float32, 100, nc),
     B_init = (nc, nd) -> Flux.onehotbatch(rand(1:nc, nd), 1:nc)
    ), (
     label="DenseMatrix * Full MaybeHotMatrix",
     A_init = nc -> rand(Float32, 100, nc),
     B_init = (nc, nd) -> Mill.maybehotbatch(rand(1:nc, nd), 1:nc)
    ), (
     label="PostImputing * Full MaybeHotMatrix",
     A_init = nc -> rand(Float32, 100, nc) |> PostImputingMatrix,
     B_init = (nc, nd) -> Mill.maybehotbatch(rand(1:nc, nd), 1:nc)
    ), (
     label="PostImputing * Half Full MaybeHotMatrix",
     A_init = nc -> rand(Float32, 100, nc) |> PostImputingMatrix,
     B_init = (nc, nd) -> Mill.maybehotbatch([rand() < 0.5 ? missing : rand(1:nc) for _ in 1:nd], 1:nc)
    )
]

params = Base.Iterators.product([10, 50, 100], [100, 200, 300]) |> collect |> vec
results = [Float64[] for _ in eachindex(experiments)]
for (nd, nc) in params
    @show nd, nc
    for (i,e) in experiments |> enumerate
        Random.seed!(0)
        A = e.A_init(nc)
        B = e.B_init(nc, nd)
        b = @benchmark $A * $B
        push!(results[i], median(b.times))
    end
end

xticklabels = [latexstring("$a|$b") for (a,b) in params]
groupedbar(xticklabels, hcat(results...); ylabel="ns", xlabel="Number of samples | Number of classes",
           labels=getindex.(experiments, :label) |> permutedims, title="A * B")
savefig("maybehot.png")

params = Base.Iterators.product([10, 50, 100], [100, 200, 300]) |> collect |> vec
results = [Float64[] for _ in eachindex(experiments)]
for (nd, nc) in params
    @show nd, nc
    for (i,e) in experiments |> enumerate
        Random.seed!(0)
        A = e.A_init(nc)
        B = e.B_init(nc, nd)
        b = @benchmark gradient(A -> sum(A*$B), $A)
        push!(results[i], median(b.times))
    end
end

xticklabels = [latexstring("$a|$b") for (a,b) in params]
groupedbar(xticklabels, hcat(results...); ylabel="ns", xlabel="Number of samples | Number of classes",
           labels=getindex.(experiments, :label) |> permutedims, title="A * B (gradient)")
savefig("maybehot_grad.png")


text:
module Hooking

using Gallium
using Gallium: X86_64, Registers
using Gallium.Registers: ip

export hook, unhook

typealias MachTask Ptr{Void}
typealias KernReturn UInt32

immutable MemoryRegion
    @static if is_apple()
        task::MachTask
    end
    addr::Ptr{Void}
    size::UInt64
end

region_to_array(region::MemoryRegion) =
    unsafe_wrap(Array, convert(Ptr{UInt8}, region.addr), (region.size,), false)

# Windows Wrappers
@static if is_windows()
    to_page(addr, size) =
        (@assert size <= 4096;
        MemoryRegion(
            addr-(reinterpret(UInt, addr)%4096),1))

    function VirtualProtect(region, perms)
        oldPerm = Ref{UInt32}()
        succ = ccall(:VirtualProtect,stdcall,Cint,(Ptr{Void}, Csize_t, UInt32, Ptr{UInt32}), region.addr, 4096, perms, oldPerm)
        succ == 0 && error(Libc.GetLastError())
        nothing
    end

    immutable MEMORY_BASIC_INFORMATION64
        BaseAddress::UInt64
        AllocationBase::UInt64
        AllocationProtect::UInt32
        __alignment1::UInt32
        RegionSize::UInt64
        State::UInt32
        Protect::UInt32
        Type::UInt32
        __alignment2::UInt32
    end

    function mem_validate(addr, length)
        ret = Ref{MEMORY_BASIC_INFORMATION64}()
        succ = ccall(:VirtualQuery, stdcall, Cint, (Ptr{Void}, Ptr{MEMORY_BASIC_INFORMATION64}, Csize_t), Ptr{Void}(addr), ret, sizeof(MEMORY_BASIC_INFORMATION64))
        if succ == 0
            error(Libc.GetLastError())
        end
        (ret[].Protect & (PAGE_EXECUTE_READ | PAGE_EXECUTE_READWRITE | PAGE_READONLY | PAGE_READWRITE)) != 0
    end

    const PAGE_READWRITE         = 0x04
    const PAGE_READONLY          = 0x02
    const PAGE_EXECUTE_READ      = 0x20
    const PAGE_EXECUTE_READWRITE = 0x40
    const PAGE_EXECUTE_WRITECOPY = 0x80
    const PAGE_NOACCESS          = 0x01

    const MEM_COMMIT             = 0x1000
    const MEM_RESERVE            = 0x2000
end

# mach vm wrappers
@static if is_apple()

    mach_task_self() = ccall(:mach_task_self,MachTask,())

    to_page(addr, size) =
        (@assert size <= 4096;
        MemoryRegion(mach_task_self(),
            addr-(reinterpret(UInt, addr)%4096),1))

    const KERN_SUCCESS     = 0x0

    const VM_PROT_NONE     = 0x0
    const VM_PROT_READ     = 0x1
    const VM_PROT_WRITE    = 0x2
    const VM_PROT_EXECUTE  = 0x4
    const VM_PROT_COPY     = 0x8

    function mach_vm_protect(task::MachTask, addr::Ptr{Void}, size::UInt64,
        prots; set_maximum::Bool = false)

        ccall(:mach_vm_protect, KernReturn,
            (MachTask, Ptr{Void}, UInt64, Bool, UInt32),
            task, addr, size, set_maximum, prots)
    end
    mach_vm_protect(region, prots) = mach_vm_protect(region.task,
        region.addr, region.size, prots)

    const VM_FLAGS_FIXED    = 0x0000
    const VM_FLAGS_ANYWHERE = 0x0001
    const VM_FLAGS_PURGABLE = 0x0002
    const VM_FLAGS_NO_CACHE = 0x0010

    function mach_vm_allocate(task, size; addr = C_NULL)
        x = Ref{Ptr{Void}}()
        x[] = addr
        ret = ccall(:mach_vm_allocate, KernReturn,
            (MachTask, Ref{Ptr{Void}}, Csize_t, Cint), task, x, size,
            addr == C_NULL ? VM_FLAGS_ANYWHERE : VM_FLAGS_FIXED)
        (ret, MemoryRegion(task,x[],size))
    end
    mach_vm_allocate(size) = mach_vm_allocate(mach_task_self(),size)

    function mach_check(status,args...)
        if status != KERN_SUCCESS
            error("Mach system call failed (error code $status)")
        end
        @assert length(args) <= 1
        length(args) == 1 ? args[1] : nothing
    end

    const VM_REGION_BASIC_INFO_64 = 9
    const KERN_INVALID_ADDRESS = 1
    function mem_validate(addr, length)
        x = Array(UInt8,2)
        addr_info = Ref{Ptr{Void}}(addr)
        size_info = Ref{Csize_t}(length)
        flavor = VM_REGION_BASIC_INFO_64
        buf = Array(UInt,1000)
        buf_size = Ref{Csize_t}(1000)
        object_name = Ref{Ptr{Void}}(C_NULL)
        ret = ccall(:mach_vm_region, KernReturn, (Ptr{Void}, Ptr{Ptr{Void}}, Ptr{Csize_t}, Csize_t, Ptr{Void}, Ptr{Csize_t}, Ptr{Void}),
            mach_task_self(), addr_info, size_info, flavor, buf, buf_size, object_name)
        ret == 0 &&
            addr_info[] <= Ptr{Void}(addr) &&
            Ptr{Void}(addr + length) < (addr_info[] + size_info[])
    end

end

@static if is_linux()
    to_page(addr, size) = (@assert size <= 4096;
        MemoryRegion(addr-(reinterpret(UInt, addr)%4096),1))

    mprotect(region, flags) = ccall(:mprotect, Cint,
        (Ptr{Void}, Csize_t, Cint), region.addr, region.size, flags)

    function mem_validate(addr, length)
        x = Array(UInt8,2)
        ret = ccall(:mincore, Cint, (Ptr{Void}, Csize_t, Ptr{UInt8}), to_page(addr, length).addr, length, x)
        ret == 0
    end

    const PROT_READ	 =  0x1
    const PROT_WRITE =  0x2
    const PROT_EXEC	 =  0x4
    const PROT_NONE	 =  0x0
end

include("backtraces.jl")

immutable Hook
    addr::Ptr{Void}
    orig_data::Vector{UInt8}
    callback::Function
end

using Base.llvmcall
hooks = Dict{Ptr{Void},Hook}()

immutable Deopt
    addr::Ptr{Void}
end

# The text section of jumpto-x86_64-macho.o minus one byte
const resume_length = 0x6c

function hook_asm_template(hook_addr, target_addr; call = true)
    diff = (target_addr - (hook_addr + 5))%Int32
    @show diff
    @assert typemin(Int32) <= diff <= typemax(Int32)
    [ call ? 0xe8 : 0xe9; reinterpret(UInt8, [Int32(diff)]); ]
end

function instrument_jmp_template(alt_stack, entry_func, exit_func)
    instrument_exit = exit_func != 0
    UInt8[
        0x41; 0x53; 	   # pushq	%r11
        0x49; 0x89; 0xe3   # movq %rsp, %r11
        # movq alt_stack, %rsp
        0x48; 0xbc; reinterpret(UInt8, [alt_stack]);
        0x41; 0x53; 	   # pushq	%r11
        0x41; 0x53; 	   # pushq	%r11
        # movq entry_func, %r11
        0x49; 0xbb; reinterpret(UInt8, [entry_func]);
        0x41; 0xff; 0xd3;       # callq *%r11
        0x41; 0x5b;             # popq	%r11
        (instrument_exit ? [
            # movq	8(%r11), %r11
            0x4d; 0x8b; 0x5b; 0x08;
            0x41; 0x53; 	   # pushq	%r11
            0x41; 0x5b;        # popq	%r11
        ] : []);
        0x5c;                   # popq %rsp
        instrument_exit ? [
            # movq exit_func, %r11
            0x49; 0xbb; reinterpret(UInt8, [exit_func]);
            0x4c; 0x89; 0x5c; 0x24; 0x08; # movq	%r11, 8(%rsp)
        ] : [];
        0x41; 0x5b;             # popq	%r11
    ]
end

function return_hook_template(alt_stack, func)
    @assert func != 0
    [
        0x41; 0x53;         # pushq	%r11
        0x41; 0x53;         # pushq	%r11
        0x49; 0x89; 0xe3    # movq %rsp, %r11
        # movq alt_stack-0x10, %rsp
        0x48; 0xbc; reinterpret(UInt8, [alt_stack-0x10]);
        0x41; 0x54;        # pushq  %r12
        0x41; 0x53; 	   # pushq	%r11
        0x49; 0xbb; reinterpret(UInt8, [func]);
        0x41; 0xff; 0xd3;  # callq *%r11
        0x41; 0x5b;        # popq	%r11
        # movq	8(%rsp), %r12
        0x4c; 0x8b; 0x64; 0x24; 0x08;
        0x4d; 0x89; 0x63; 0x08; # movq	%r12, 8(%r11)
        0x41; 0x5c;             # popq	%r12
        0x4c; 0x89; 0xdc;       # movq	%r11, %rsp
        0x41; 0x5b;             # popq	%r11
        0xc3;                   # retq
    ]
end

function hook_asm_template(addr)
    [
        0x90; # 0xcc
        0x50; #pushq   %rax
        # movq $hookto, %rax
        0x48; 0xb8; reinterpret(UInt8, [addr]);
        0xff; 0xd0; #callq *%rax
    ]
end
hook_asm_template() = (global thehook; hook_asm_template(thehook))
const hook_length = length(hook_asm_template(UInt(0)))

function hook_tail_template(extra_instructions, ret_addr)
    addr_bytes = reinterpret(UInt8,[ret_addr])
    [
    # Is this a good idea? Probably not
    extra_instructions...,
    0x66, 0x68, addr_bytes[7:8]...,
    0x66, 0x68, addr_bytes[5:6]...,
    0x66, 0x68, addr_bytes[3:4]...,
    0x66, 0x68, addr_bytes[1:2]...,
    0xc3
    ]
end

function aligned_xsave_RC()
    RCnew = Array(UInt8,length(X86_64.basic_regs)*sizeof(UInt64) +
        sizeof(fieldtype(X86_64.ExtendedRegs,:xsave_state)) + 64)

    rcptr = pointer(RCnew)
    # Align to 56 bytes (so the xsave state ends up at 64 byte alignment)
    npad = (64 - (UInt64(rcptr) % 64)) - 8
    @assert npad > 0
    rcptr += npad
    @assert (UInt64(rcptr) % 64) == 56

    RCnew, rcptr
end

# Split this out to avoid constructing a gc frame in the callback directly
@noinline function _callback(x::Ptr{Void})
    RC = X86_64.BasicRegs()
    regs = unsafe_wrap(Array, Ptr{UInt64}(x), (length(X86_64.basic_regs),), false)
    for i in X86_64.basic_regs
        set_dwarf!(RC, i, RegisterValue{UInt64}(regs[i+1], (-1%UInt64)))
    end
    xsave_ptr = Ptr{UInt8}(x+sizeof(Ptr{Void})*length(X86_64.basic_regs))
    RC = X86_64.ExtendedRegs(RC,
        unsafe_load(Ptr{fieldtype(X86_64.ExtendedRegs,:xsave_state)}(xsave_ptr)))
    hook_addr = UInt(ip(RC))-hook_length
    hook = hooks[reinterpret(Ptr{Void},hook_addr)]
    cb_RC = copy(RC)
    set_ip!(cb_RC, hook_addr+1)
    ret = hook.callback(hook, cb_RC)
    if isa(ret, Deopt)
        ret_addr = ret.addr
        extra_instructions = []
    else
        ret_addr = hook_addr+length(hook.orig_data)
        extra_instructions = hook.orig_data
    end
    resume_data = [
        #0xcc,
        resume_instructions;
        # Counteract the pushq %rip in the resume code
        0x48; 0x83; 0xc4; 0x8; # addq $8, %rsp
        hook_tail_template(extra_instructions, ret_addr)
    ]
    global callback_rwx
    callback_rwx[1:length(resume_data)] = resume_data

    # invalidate instruction cache here if ever ported to other
    # architectures

    RCnew, rcptr = aligned_xsave_RC()
    npad = UInt(rcptr - pointer(RCnew))

    # For alignment purposes
    RCnew[1:npad] = 0
    RCnew[npad+1:end] = UInt8[reinterpret(UInt8,
        UInt64[
            [get_dwarf(RC, i)[] for i in X86_64.basic_regs];
            reinterpret(UInt64,[RC.xsave_state])]);
        UInt8[0 for i = 1:(64-npad)]]

    ptr = convert(Ptr{Void},pointer(callback_rwx))::Ptr{Void}
    ptr, Ptr{UInt64}(rcptr)::Ptr{UInt64}
end
function callback(x::Ptr{Void})
    enabled = gc_enable(false)
    ptr, data = _callback(x)
    gc_enable(enabled)
    # jump to resume code
    ccall(ptr,Void,(Ptr{Void},),data)
    nothing
end

const hooking_lib = joinpath(dirname(@__FILE__),"hooking")

function __init__()
    global resume
    global thehook
    global callback_rwx
    global resume_instructions
    here = dirname(@__FILE__)
    function resume(RC)
        ccall((:hooking_jl_jumpto, hooking_lib),Void,(Ptr{UInt8},),pointer(RC.data))
    end
    theresume = cglobal((:hooking_jl_jumpto, hooking_lib), Ptr{UInt8})
    resume_instructions = unsafe_wrap(Array, convert(Ptr{UInt8}, theresume),
        (resume_length,), false)
    # Allocate an RWX page for the callback return
    callback_rwx = @static if is_apple()
        region = mach_check(mach_vm_allocate(4096)...)
        mach_check(mach_vm_protect(region,
            VM_PROT_READ | VM_PROT_WRITE | VM_PROT_EXECUTE))
        region_to_array(region)
    elseif is_linux()
        region = MemoryRegion(ccall(:mmap, Ptr{Void},
            (Ptr{Void}, Csize_t, Cint, Cint, Cint, Csize_t),
            C_NULL, 4096, PROT_EXEC | PROT_READ | PROT_WRITE,
            Base.Mmap.MAP_ANONYMOUS | Base.Mmap.MAP_PRIVATE,
            -1, 0), 4096)
        Base.systemerror("mmap", reinterpret(Int, region.addr) == -1)
        region_to_array(region)
    elseif is_windows()
        region = MemoryRegion(ccall(:VirtualAlloc, Ptr{Void},
            (Ptr{Void}, Csize_t, UInt32, UInt32),
            C_NULL, 4096, MEM_RESERVE | MEM_COMMIT, PAGE_EXECUTE_READWRITE), 4096)
        Base.systemerror("VirtualAlloc", reinterpret(Int, region.addr) == 0)
        region_to_array(region)
    end
    thecallback = Base.cfunction(callback,Void,Tuple{Ptr{Void}})::Ptr{Void}
    ccall((:hooking_jl_set_callback, hooking_lib), Void, (Ptr{Void},),
        thecallback)
    thehook = cglobal((:hooking_jl_savecontext, hooking_lib), Ptr{UInt8})
end

# High Level Implementation

# Temporarily allow writing to an executable page.
# It would be nice to have a general version of this, but unfortunately, it
# seems the only reliable version to write something to a protected page on
# linux is to either parse the /proc mappings file or to use ptrace, neither
# of which sounds like a lot of fun. For now, just do this and assume the page
# is executable to begin with.
function allow_writing(f, region)
    # On OS X, make sure that the page is mapped as COW
    @static if is_apple()
        mach_check(mach_vm_protect(region, VM_PROT_READ | VM_PROT_WRITE | VM_PROT_EXECUTE))
    elseif is_linux()
        mprotect(region, PROT_READ | PROT_WRITE | PROT_EXEC)
    elseif is_windows()
        VirtualProtect(region, PAGE_EXECUTE_READWRITE)
    end
    f()
    @static if is_apple()
        mach_check(mach_vm_protect(region, VM_PROT_EXECUTE | VM_PROT_READ))
    elseif is_linux()
        mprotect(region, PROT_READ | PROT_EXEC)
    elseif is_windows()
        VirtualProtect(region, PAGE_EXECUTE_READ)
    end
end

function determine_nbytes_to_replace(bytes_required, addr::Ptr{Void})
    triple = "x86_64-apple-darwin15.0.0"
    DC = ccall(:jl_LLVMCreateDisasm, Ptr{Void},
        (Ptr{UInt8},Ptr{Void},Cint,Ptr{Void},Ptr{Void}),
        triple, C_NULL, 0, C_NULL, C_NULL)
    @assert DC != C_NULL

    nbytes = 0
    template = hook_asm_template()
    while nbytes < bytes_required
        outs = Ref{UInt8}()
        nbytes += ccall(:jl_LLVMDisasmInstruction, Csize_t,
            (Ptr{Void}, Ptr{UInt8}, Csize_t, UInt64, Ptr{UInt8}, Csize_t),
            DC,          # Disassembler
            addr+nbytes, # bytes
            30,          # Size
            addr+nbytes, # PC
            outs, 1      # OutString
            )
    end

    nbytes
end
function determine_nbytes_to_replace(bytes_required, orig_bytes)
    determine_nbytes_to_replace(bytes_required, Ptr{Void}(pointer(orig_bytes)))
end

function hook(callback::Function, addr)
    # Compute number of bytes by disassembly
    # Ideally we would also check for uses of rip and branches here and error
    # out if any are found, but for now we don't need to

    template = hook_asm_template()
    nbytes = determine_nbytes_to_replace(length(template), addr)
    # Record the instructions that were there originally
    dest = unsafe_wrap(Array, convert(Ptr{UInt8}, addr), (nbytes,), false)
    orig_data = copy(dest)

    hook_asm = [ template; fill(0xcc,nbytes-length(template)) ]# Pad to nbytes
    @assert length(hook_asm) == length(orig_data)

    allow_writing(to_page(addr,nbytes)) do
        dest[:] = hook_asm
    end

    hooks[addr] = Hook(addr,orig_data,callback)
end

function hook(thehook::Hook)
    nbytes = length(thehook.orig_data)
    template = hook_asm_template();
    hook_asm = [ template; fill(0x90,nbytes-length(template)) ]# Pad to nbytes

    dest = unsafe_wrap(Array, convert(Ptr{UInt8}, thehook.addr),
        (nbytes,), false)
    allow_writing(to_page(thehook.addr,nbytes)) do
        dest[:] = hook_asm
    end

    hooks[thehook.addr] = thehook
end

function get_function_addr(t)
    llvmf = ccall(:jl_get_llvmf, Ptr{Void}, (Any, Bool, Bool), t, false, true)
    @assert llvmf != C_NULL
    fptr = ccall(:jl_get_llvm_fptr, UInt64, (Ptr{Void},), llvmf)
    @assert fptr != 0
    reinterpret(Ptr{Void},fptr)
end
get_function_addr(f, t) = get_function_addr(Tuple{typeof(f), Base.to_tuple_type(t).parameters...})

hook(callback, f, t) = hook(callback, get_function_addr(f, t))

function unhook(addr::Union{Ptr{Void},UInt64})
    hook = pop!(hooks, addr)

    nbytes = length(hook.orig_data)
    dest = unsafe_wrap(Array, convert(Ptr{UInt8}, addr),
        (nbytes,), false)

    allow_writing(to_page(addr,nbytes)) do
        dest[:] = hook.orig_data
    end
end

unhook(hook::Hook) = unhook(hook.addr)

@inline function getcontext()
    RCnew, rcptr = aligned_xsave_RC()
    ccall((:hooking_jl_simple_savecontext, hooking_lib),Void,(Ptr{UInt64},),rcptr)
    RC = X86_64.BasicRegs()
    for i in X86_64.basic_regs
        set_dwarf!(RC, i, RegisterValue{UInt64}(unsafe_load(Ptr{UInt64}(rcptr),i+1), (-1%UInt64)))
    end
    xsave_ptr = Ptr{UInt8}(rcptr+sizeof(Ptr{Void})*length(X86_64.basic_regs))
    RC = X86_64.ExtendedRegs(RC,
        unsafe_load(Ptr{fieldtype(X86_64.ExtendedRegs,:xsave_state)}(xsave_ptr)))
    RC
end

end # module


text:
using Infinity
"""
    One-NN(x, dist, data, labels)
    Devuelve la clase a la que clasifica x
    a partir del vecino más cercano dentro del conjunto de datos.
"""
function OneNN(x, dist, data, labels)
    distanciaMinima = ∞
    clase = nothing

    for (index,element) in enumerate(eachrow(data))
        d = dist(element, x)
        if d < distanciaMinima
            distanciaMinima = d
            clase = labels[index]
        end
    end
    return clase
end


text:
const type_to_label = Dict{String, String}([
    "article"         => "journal",
    "book"            => "book",
    "booklet"         => "booklet",
    "eprint"          => "eprint",
    "inbook"          => "book chapter",
    "incollection"    => "book section",
    "inproceedings"   => "conference",
    "manual"          => "manual",
    "masterthesis"    => "master thesis",
    "misc"            => "other",
    "phdthesis"       => "doctoral thesis",
    "proceedings"     => "proceedings",
    "techreport"      => "report",
    "unpublished"     => "other"
])

function xtitle(entry::T) where T <: BibInternal.Entry
    return :title ∈ fieldnames(typeof(entry)) ? entry.title : get(entry.fields, "title", "")
end

function xnames(
    entry::BibInternal.Entry,
    editors::Bool=false;
    names::Symbol=:full # Current options: :last, :full
    )
    # forces the names to be editors' name if the entry are Proceedings
    if !editors && entry.type ∈ ["proceedings"]
        return xnames(entry, true)
    end
    entry_names = editors ? entry.editors : entry.authors
    str = ""

    start = true
    for s in entry_names
        str *= start ? "" : ", "
        start = false
        if names == :last
            str *= s.particle * " " * s.last * " " * s.junior
        else
            str *= s.first * " " * s.middle * " " * s.particle * " " * s.last * " " * s.junior
        end
    end
    # str *= editors ? ", editors" : ""
    return replace(str, r"[\n\r ]+" => " ") # TODO: make it cleaner (is replace still necessary)
end

function xin(entry::BibInternal.Entry)
    str = ""
    if entry.type == "article"
        str *= entry.in.journal * ", " * entry.in.volume
        str *= entry.in.number != "" ? "($(entry.in.number))" : ""
        str *= entry.in.pages != "" ? ", $(entry.in.pages)" : ""
        str *= ", " * entry.date.year
    elseif entry.type == "book"
        str *= entry.in.publisher
        str *= entry.in.address != "" ? ", $(entry.in.address)" : ""
        str *= ", " * entry.date.year
    elseif entry.type == "booklet"
        str *= entry.access.howpublished * ", " * entry.date.year
    elseif entry.type == "eprint"
        if entry.eprint.archive_prefix == ""
            str *= entry.eprint.eprint
        else
            str *= "$(entry.eprint.archive_prefix):$(entry.eprint.eprint) [$(entry.eprint.primary_class)]"
        end
    elseif entry.type == "inbook"
        aux = entry.in.chapter == "" ? entry.in.pages : entry.in.chapter
        str *= entry.booktitle * ", " * aux * ", " * entry.in.publisher
        str *= entry.in.address != "" ? ", " * entry.in.address : ""
    elseif entry.type == "incollection"
        str = "In " * xnames(entry, true) * ", editors, " * entry.booktitle * ", " * entry.in.pages * "."
        str *= " " * entry.in.publisher
        str *= entry.in.address != "" ? ", $(entry.in.address)" : ""
        str *= ", " * entry.date.year
    elseif entry.type == "inproceedings"
        str *= " In " * entry.booktitle
        str *= entry.in.series != "" ? ", " * entry.in.series : ""
        str *= entry.in.pages != "" ? ", $(entry.in.pages)" : ""
        str *= entry.in.address != "" ? ", $(entry.in.address)" : ""
        str *= ", " * entry.date.year
        str *= entry.in.publisher != "" ? ". $(entry.in.publisher)" : ""
    elseif entry.type == "manual"
        str *= entry.in.organization
        str *= entry.in.address != "" ? ", $(entry.in.address)" : ""
        str *= ", " * entry.date.year
    elseif entry.type ∈ ["masterthesis", "phdthesis"]
        str *= entry.type == "masterthesis" ? "Master's" : "PhD"
        str *= " thesis, " * entry.in.school
        str *= entry.in.address != "" ? ", $(entry.in.address)" : ""
        str *= ", " * entry.date.year
    elseif entry.type == "misc"
        aux = entry.access.howpublished != "" != entry.date.year ? ", " : ""
        str *= entry.access.howpublished * aux * entry.date.year
        str *= aux != "" && get(entry.fields, "note", "") != "" ? ". " : ""
        str *= get(entry.fields, "note", "")
    elseif entry.type == "proceedings"
        str *= entry.in.volume != "" ? "Volume " * entry.in.volume * " of " : ""
        str *= entry.in.series
        str *= entry.in.address != "" ? ", $(entry.in.address)" : ""
        str *= ", " * entry.date.year
        str *= entry.in.publisher != "" ? ". $(entry.in.address)" : ""
    elseif entry.type == "techreport"
        str *= entry.in.number != "" ? "Technical Report " * entry.in.number * ", " : ""
        str *= entry.in.institution
        str *= entry.in.address != "" ? ", $(entry.in.address)" : ""
        str *= ", " * entry.date.year
    elseif entry.type == "unpublished"
        aux = get(entry.fields, "note", "")
        str *= aux != "" != entry.date.year ? aux * ", " : ""
        str *= entry.date.year
    end
    str *= str == "" ? "" : "."
    return str
end

function xyear(entry::BibInternal.Entry)
    return entry.date.year
end

function xlink(entry::BibInternal.Entry)
    if entry.access.doi != ""
        return "https://doi.org/" * entry.access.doi
    elseif entry.access.url != ""
        return entry.access.url
    end
    return ""
end

function xfile(entry::BibInternal.Entry)
    return "files/$(entry.id).pdf"
end

function xcite(entry::BibInternal.Entry)
    string(entry)
end

function xlabels(entry::BibInternal.Entry)
    str = get(entry.fields, "swp-labels", "")
    str = str == "" ? get(entry.fields, "labels", "") : str
    return str == "" ? [entry.type] : split(str, r"[\n\r ]*,[\n\r ]*")
end

struct Publication
    id::String
    type::String
    title::String
    names::String
    in::String
    year::String
    link::String
    file::String
    cite::String
    labels::Vector{String}
end

function Publication(entry::BibInternal.Entry)
    id = entry.id
    type = entry.type
    title = entry.title
    names = xnames(entry)
    in_ = xin(entry)
    year = xyear(entry)
    link = xlink(entry)
    file = xfile(entry)
    cite = export_bibtex(entry)
    labels = xlabels(entry)
    return Publication(id, type, title, names, in_, year, link, file, cite, labels)
end

"""
    export_web(bibliography::DataStructures.OrderedDict{String,BibInternal.Entry})
Export a biblography in internal format to the web format of the [StaticWebPages.jl](https://github.com/Azzaare/StaticWebPages.jl) pakcage.
"""
function export_web(
    bibliography::DataStructures.OrderedDict{String,BibInternal.Entry}
    )
    entries = Vector{Publication}()
    for entry in values(bibliography)
        p = Publication(entry)
        push!(entries, p)
    end
    return entries
end


text:
using Random
using LinearAlgebra
using Statistics
using Optim
using DataFrames
using DataFramesMeta
using CSV
using HTTP
using GLM

# read in function to create state transitions for dynamic model
include("create_grids.jl")

function ps6()
#1)
url = "https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2020/master/ProblemSets/PS5-ddc/busdata.csv"
df = CSV.read(HTTP.get(url).body,DataFrame)

# create bus id variable
df = @transform(df, bus_id = 1:size(df,1))

#---------------------------------------------------
# reshape from wide to long (must do this twice be-
# cause DataFrames.stack() requires doing it one
# variable at a time)
#---------------------------------------------------
# first reshape the decision variable
dfy = @select(df, :bus_id,:Y1,:Y2,:Y3,:Y4,:Y5,:Y6,:Y7,:Y8,:Y9,:Y10,:Y11,:Y12,:Y13,:Y14,:Y15,:Y16,:Y17,:Y18,:Y19,:Y20,:RouteUsage,:Branded)
dfy_long = DataFrames.stack(dfy, Not([:bus_id,:RouteUsage,:Branded]))
rename!(dfy_long, :value => :Y)
dfy_long = @transform(dfy_long, time = kron(collect([1:20]...),ones(size(df,1))))
select!(dfy_long, Not(:variable))

# next reshape the odometer variable
dfx = @select(df, :bus_id,:Odo1,:Odo2,:Odo3,:Odo4,:Odo5,:Odo6,:Odo7,:Odo8,:Odo9,:Odo10,:Odo11,:Odo12,:Odo13,:Odo14,:Odo15,:Odo16,:Odo17,:Odo18,:Odo19,:Odo20)
dfx_long = DataFrames.stack(dfx, Not([:bus_id]))
rename!(dfx_long, :value => :Odometer)
dfx_long = @transform(dfx_long, time = kron(collect([1:20]...),ones(size(df,1))))
select!(dfx_long, Not(:variable))

# join reshaped df's back together
df_long = leftjoin(dfy_long, dfx_long, on = [:bus_id,:time])
sort!(df_long,[:bus_id,:time])

#2)
#[print(string(":", x, ", ")) for x in names(df_long)]

θ̂_glm = glm(@formula(Y ~ Odometer^2 * RouteUsage^2 * Branded * time^2 ), df_long, Binomial(), LogitLink())
println(θ̂_glm)

#3)

# state transition matrices
zval,zbin,xval,xbin,xtran = create_grids()

df = CSV.read(HTTP.get(url).body,DataFrame)
Y = Matrix(df[:,[:Y1,:Y2,:Y3,:Y4,:Y5,:Y6,:Y7,:Y8,:Y9,:Y10,:Y11,:Y12,:Y13,:Y14,:Y15,:Y16,:Y17,:Y18,:Y19,:Y20]])
X = Matrix(df[:,[:Odo1,:Odo2,:Odo3,:Odo4,:Odo5,:Odo6,:Odo7,:Odo8,:Odo9,:Odo10,:Odo11,:Odo12,:Odo13,:Odo14,:Odo15,:Odo16,:Odo17,:Odo18,:Odo19,:Odo20]])
Z = Vector(df[:,:RouteUsage])
B = Vector(df[:,:Branded])
N = size(Y,1)
T = size(Y,2)
Xstate = Matrix(df[:,[:Xst1,:Xst2,:Xst3,:Xst4,:Xst5,:Xst6,:Xst7,:Xst8,:Xst9,:Xst10,:Xst11,:Xst12,:Xst13,:Xst14,:Xst15,:Xst16,:Xst17,:Xst18,:Xst19,:Xst20]])
Zstate = Vector(df[:,:Zst])

df1 = DataFrame(
Odometer = kron(ones(zbin),xval),
RouteUsage = kron(ones(xbin),zval),
Branded = 0,
time = 0
)

@views @inbounds function likebus2(θ̂_glm,T,Zstate,Xstate,xtran,zbin,xbin,xval,df1,df)
    FV = zeros(size(xtran,1),2,T+1)
    FVT1 = zeros(size(df,1), T)

    for t=2:T
        for b=0:1
            df1[t,"time"] = t
            df1[t,"Branded"] = b
            p0 = predict(θ̂_glm, df1)
            FV[:, b+1, t] = - 0.9 .*log.(p0)
        end
    end
    for i=1:size(df,1)
        row0 = (Zstate[i]-1)*xbin+1 # this is the same argument as the index of xtran in v0 above, but we use the actual Z
        for t=1:T
            row1  = Xstate[i,t] + (Zstate[i]-1)*xbin                                                                # this is the same as row in the first loop, except we use the actual X and Z
            FVT1[i,t] = (xtran[row1,:].-xtran[row0,:])'*FV[row0:row0+xbin-1,df[i,"Branded"]+1,t+1]
        end
    end
    return FVT1'[:]
end

fvt1 = likebus2(θ̂_glm,T,Zstate,Xstate,xtran,zbin,xbin,xval,df1,df)
df_long = @transform(df_long,fv = fvt1)

theta_hat_ccp_glm = glm(@formula(Y ~ Odometer + Branded),df_long, Binomial(), LogitLink(),offset=df_long.fv)
println(theta_hat_ccp_glm)

end

ps6()


text:
using KernelAbstractions: NoneEvent

using Oceananigans.Utils: arch_array
using Oceananigans.Grids: AbstractGrid

abstract type AbstractTimeDiscretization end

"""
    struct ExplicitTimeDiscretization <: AbstractTimeDiscretization

Represents fully-explicit time-discretization of a `TurbulenceClosure`.
"""
struct ExplicitTimeDiscretization <: AbstractTimeDiscretization end

Base.summary(::ExplicitTimeDiscretization) = "ExplicitTimeDiscretization"

"""
    struct VerticallyImplicitTimeDiscretization <: AbstractTimeDiscretization

Represents vertically-implicit time-discretization of a `TurbulenceClosure`.

This imples that a flux divergence such as ``∇ ⋅ q`` at the n-th timestep is
time-discretized as

```julia
[∇ ⋅ q]ⁿ = [explicit_flux_divergence]ⁿ + [∂z (κ ∂z c)]ⁿ⁺¹
```
"""
struct VerticallyImplicitTimeDiscretization <: AbstractTimeDiscretization end

Base.summary(::VerticallyImplicitTimeDiscretization) = "VerticallyImplicitTimeDiscretization"

@inline time_discretization(::AbstractTurbulenceClosure{TimeDiscretization}) where TimeDiscretization = TimeDiscretization()
@inline time_discretization(::Nothing) = ExplicitTimeDiscretization() # placeholder for closure::Nothing

#####
##### Explicit: move along, nothing to worry about here (use fallbacks).
#####

const ATD = AbstractTimeDiscretization

@inline diffusive_flux_x(i, j, k, grid, ::ATD, args...) = diffusive_flux_x(i, j, k, grid, args...)
@inline diffusive_flux_y(i, j, k, grid, ::ATD, args...) = diffusive_flux_y(i, j, k, grid, args...)
@inline diffusive_flux_z(i, j, k, grid, ::ATD, args...) = diffusive_flux_z(i, j, k, grid, args...)

@inline viscous_flux_ux(i, j, k, grid, ::ATD, args...) = viscous_flux_ux(i, j, k, grid, args...)
@inline viscous_flux_uy(i, j, k, grid, ::ATD, args...) = viscous_flux_uy(i, j, k, grid, args...)
@inline viscous_flux_uz(i, j, k, grid, ::ATD, args...) = viscous_flux_uz(i, j, k, grid, args...)

@inline viscous_flux_vx(i, j, k, grid, ::ATD, args...) = viscous_flux_vx(i, j, k, grid, args...)
@inline viscous_flux_vy(i, j, k, grid, ::ATD, args...) = viscous_flux_vy(i, j, k, grid, args...)
@inline viscous_flux_vz(i, j, k, grid, ::ATD, args...) = viscous_flux_vz(i, j, k, grid, args...)

@inline viscous_flux_wx(i, j, k, grid, ::ATD, args...) = viscous_flux_wx(i, j, k, grid, args...)
@inline viscous_flux_wy(i, j, k, grid, ::ATD, args...) = viscous_flux_wy(i, j, k, grid, args...)
@inline viscous_flux_wz(i, j, k, grid, ::ATD, args...) = viscous_flux_wz(i, j, k, grid, args...)


text:
abstract type Optimizer end
export Optimizer

# optimizers with L1 and L2 weight decay
include("./1-SGDL1L2.jl")
include("./2-MomentumL1L2.jl")
include("./3-AdamL1L2.jl")
include("./4-AdaGradL1L2.jl")
include("./5-RMSPropL1L2.jl")
export SGD
export Momentum
export Adam
export AdaGrad
export RMSProp

# auto gradient clipping
include("./auto-grad-cliper.jl")
export AutoGradCliper
export AutoGradNormCliper

export decay
export normclip
export LpNormClip
export L2NormClip
export L1NormClip
export L0NormClip
export LPInfNormClip
export LNInfNormClip
export setNanInfZero
export lrarray
export clip!

"""
    clip!(::Vector{XVariable}, kind='u'; L1decay=0.0, L2decay=0.0, clipvalue=1.0)

Limit the amplitude of parameters. `kind` has four options:\n
`'u'` for recurrent params\n
`'b'` for bias params\n
`'w'` for projection params\n
`'a'` for `'u'`, `'b'` and `'w'` params\n
as show in `yᵗ = f(w*xᵗ + u*hᵗ⁻¹ + b)` or other similar formulas
"""
function clip!(xparams::Vector{XVariable}, kind='u'; L1decay=0.0, L2decay=0.0, clipvalue=1.0)
    @assert clipvalue>0 "clipvalue is positive, but got $clipvalue"
    if !(kind=='u' || kind=='b' || kind=='w' || kind=='a')
        @error "type of XVariable not among u/w/b/a, but got $kind"
    end

    λ₁ = -L1decay
    λ₂ = -L2decay
    for (c, θ) in xparams
        if c == kind || kind=='a'
            𝒗 = ᵛ(θ)
            i = abs.(𝒗) .> clipvalue
            if sum(i) == 0 continue end
            if λ₁==0 && λ₂==0                     # Hard truncation
                @. 𝒗[i] = clipvalue * sign(𝒗[i])
            elseif λ₁==0 && λ₂!=0                 # Soft truncation (L2)
                @. 𝒗[i] += λ₂ * 𝒗[i]
            elseif λ₁!=0 && λ₂==0                 # Gradual truncation (L1)
                @. 𝒗[i] += λ₁ * sign(𝒗[i])
            else  # λ₁!=0 && λ₂!=0
                @. 𝒗[i] += λ₁ * sign(𝒗[i]) + λ₂ * 𝒗[i]
            end
        end
    end
end

function decay(params::Vector{Variable}; ratio=0.999)
    for p in params
        p.value .*= ratio
    end
end

"""
    setNanInfZero(x)
```julia
x = randn(1,4)
x[1] = Inf;
x[2] =-Inf;
x[3] = NaN;
```
```
julia> x
1×4 Array{Float64,2}:
 Inf  -Inf  NaN  0.602655

julia> x = setNanInfZero!(x)
1×4 Array{Float64,2}:
 0.0  0.0  0.0  0.602655
```
 """
function setNanInfZero(x)
    x[ isnan.(x) .⊻ isinf.(x) ] .= 0.0
    return x
end

function L2NormClip(x::AbstractArray, clipvalue)
    pnorm = sqrt(sum(x.^2) / length(x))
    scale = clipvalue / pnorm
    if pnorm > clipvalue
        x .*= scale
    end
    return x
end

function L1NormClip(x::AbstractArray, clipvalue)
    pnorm = sum(abs.(x)) / length(x)
    scale = clipvalue / pnorm
    if pnorm > clipvalue
        x .*= scale
    end
    return x
end

function L0NormClip(x::AbstractArray, clipvalue)
    pnorm = sum(x .!= 0.0) / length(x)
    scale = clipvalue / pnorm
    if pnorm > clipvalue
        x .*= scale
    end
    return x
end

function LPInfNormClip(x::AbstractArray, clipvalue)
    pnorm = maximum(abs.(x))
    scale = clipvalue / pnorm
    if pnorm > clipvalue
        x .*= scale
    end
    return x
end

function LNInfNormClip(x::AbstractArray, clipvalue)
    pnorm = minimum(abs.(x))
    scale = clipvalue / pnorm
    if pnorm > clipvalue
        x .*= scale
    end
    return x
end

function LpNormClip(x::AbstractArray, clipvalue; order::Union{Int,String}=2)
    order==2 && return L2NormClip(x, clipvalue)
    order==1 && return L1NormClip(x, clipvalue)
    order==0 && return L0NormClip(x, clipvalue)
    order=="inf"  && return LPInfNormClip(x, clipvalue)
    order=="-inf" && return LNInfNormClip(x, clipvalue)
    pnorm = (sum( abs.(x).^order ) / length(x)) ^ (1/order)
    scale = clipvalue / pnorm
    if pnorm > clipvalue
        x .*= scale
    end
    return x
end

"""
    lrarray(init,final,steps;func="exp") -> Array{Float32,1}
# Arguments
- ` init`: initial learning rate
- `final`: final learning rate
- `steps`: steps to change from `init` to `final`
- ` func`: function used to change the learning rate. (e.g. exp/cos/linear)
"""
function lrarray(init,final,steps;func="exp")
    func=="exp"    && return lrarrayexp(init,final,steps)
    func=="cos"    && return lrarraycos(init,final,steps)
    func=="linear" && return lrarraylinear(init,final,steps)
end

function lrarrayexp(i,f,n)
    if n == 1
        return [i]
    else
        lr = zeros(Float32,n)
        α  = log(f/i)/(n-1)
        for x = 1:n
            lr[x] = i * exp(α * (x-1))
        end
        return lr
    end
end

function lrarraylinear(i,f,n)
    if n == 1
        return [i]
    else
        lr = zeros(Float32,n)
        α  = (f-i)/(n-1)
        for x = 1:n
            lr[x] = α*(x-1) + i
        end
        return lr
    end
end

function lrarraycos(i,f,n)
    if n == 1
        return [i]
    else
        lr = zeros(Float32,n)
        α  = pi/2/(n-1)
        Δ  = (i-f)
        for x = 1:n
            lr[x] = cos(α*(x-1)) * Δ + f
        end
        return lr
    end
end

max_stars_repo_name: sonosole/Delta

max_issues_repo_name: sonosole/Delta

max_forks_repo_name: sonosole/Delta


text:
function cyclic_paths(n)
    n = (n - 3) % 28960854
    c = BigInt(8)
    d = BigInt(8)

    for i in 1:n
        c = (3 * c)^3 % 13^8
        d = (3 * d)^3 % (6 * 13^8)
    end

    return d, c
end

function main()
    n = 10000
    local c

    for i in 1:3
        n, c = cyclic_paths(n)
    end

    println(c)
end

main()

max_stars_repo_name: HarrisonGreen/Project-Euler-Solutions

max_issues_repo_name: HarrisonGreen/Project-Euler-Solutions

max_forks_repo_name: HarrisonGreen/Project-Euler-Solutions


text:
mutable struct IndexMap{T<:Number}
    m::Dict{Any,T}
    c::Int64
end

function sanitycheck()
    return "sanity"
    println("sanity")
end

IndexMap() = IndexMap(Dict(), 0)

function newindex!(oldindex, i::IndexMap)
    if haskey(i.m, oldindex)
        return i.m[oldindex]
    else
        i.c += 1
        i.m[oldindex] = i.c
        return i.c
    end
end

function normalize(sdp::SparseSDP{Float64})
    newsdp = SparseSDP(T, maximize=ismaximizationproblem(sdp), normalized=true)

    cm = IndexMap()
    bm = IndexMap()
    ems = Dict{Any,IndexMap}()

    for (ri, value) in rhs(sdp)
        newri = newindex!(ri, cm)
        setrhs!(newsdp, newri, value)
    end

    for (bi, matrix) in blocks(obj(sdp))
        newbi = newindex!(bi, bm)
        if !haskey(ems, newbi)
            ems[newbi] = IndexMap()
        end
        for ((i, j), v) in entries(matrix)
            newi = newindex!(i, ems[newbi])
            newj = newindex!(j, ems[newbi])
            setobj!(newsdp, newbi, newi, newj, v)
        end
    end

    for (ri, blockmatrix) in cons(sdp)
        newri = newindex!(ri, cm)
        for (bi, matrix) in blocks(blockmatrix)
            newbi = newindex!(bi, bm)
            if !haskey(ems, newbi)
                ems[newbi] = IndexMap()
            end
            for ((i, j), v) in entries(matrix)
                newi = newindex!(i, ems[newbi])
                newj = newindex!(j, ems[newbi])
                setcon!(newsdp, newri, newbi, newi, newj, v)
            end
        end
    end
    nems = Dict()
    for (k,v) in ems
        nems[k] = ems[k].m
    end
    newsdp, cm.m, bm.m, nems
end

max_stars_repo_name: JAndriesJ/SemidefiniteProgramming

max_issues_repo_name: JAndriesJ/SemidefiniteProgramming

max_forks_repo_name: JAndriesJ/SemidefiniteProgramming


text:
using SimJulia
using Base.Dates

type TestException <: Exception end

function test_cb(sim::Simulation, ev::Event)
  println(value(ev))
  println("$(now()): Hi, it's now $(now(sim))")
end

function test_another_cb(sim::Simulation, ev::Event)
  println("$(now()): Hi, I am a second callback at $(now(sim))")
end

function and_cb(sim::Simulation, ev::Event)
  if isa(value(ev), Exception)
    println("Exception has been thrown!")
  else
    println("$(now()): Both events are triggered at $(now(sim))!")
    println(value(ev))
  end
end

function or_cb(sim::Simulation, ev::Event)
  println("$(now()): One of both events is triggered at $(now(sim))!")
  println(value(ev))
end

sim = Simulation(now())
ev = Event()
append_callback(ev, test_cb)
schedule(sim, ev, Day(1))
another_ev = timeout(sim, 3600000*24*2)
schedule!(sim, ev, Hour(23), value=TestException())
append_callback(another_ev, test_cb)
and_event = ev & another_ev
append_callback(and_event, and_cb)
append_callback(and_event, test_another_cb)
run(sim, Month(1))

sim = Simulation(now())
ev = Event()
append_callback(ev, test_cb)
schedule(sim, ev, Day(1))
another_ev = timeout(sim, 3600000*24*2, value=TestException())
and_event = ev & another_ev
append_callback(and_event, and_cb)
append_callback(and_event, test_another_cb)
run(sim, Month(1))

sim = Simulation(now())
ev = Event()
append_callback(ev, test_cb)
schedule(sim, ev, Day(1))
another_ev = timeout(sim, 3600000*24*2)
and_event = ev & another_ev
append_callback(and_event, and_cb)
append_callback(and_event, test_another_cb)
run(sim, Month(1))

sim = Simulation()
ev = Event()
append_callback(ev, test_cb)
schedule(sim, ev, 1)
another_ev = timeout(sim, 3, value="π-day0314")
append_callback(another_ev, test_cb)
or_event = ev | another_ev
append_callback(or_event, or_cb)
append_callback(or_event, test_another_cb)
run(sim)

sim = Simulation()
ev = Event()
append_callback(ev, test_cb)
schedule(sim, ev, 1)
or_event = ev | timeout(sim, 1)
append_callback(or_event, or_cb)
run(sim)

sim = Simulation()
ev = Event()
append_callback(ev, test_cb)
schedule(sim, ev, 1)
or_event = timeout(sim, 1, value=TestException()) | ev
append_callback(or_event, or_cb)
run(sim)

sim = Simulation(today())
ev = Event()
append_callback(ev, test_cb)
schedule!(sim, ev, 1)
another_ev = timeout(sim, Month(2))
append_callback(another_ev, test_cb)
and_event = ev & another_ev
append_callback(and_event, and_cb)
append_callback(and_event, test_another_cb)
run(sim, today() + Month(1))

function print_cb(sim::Simulation, ev::Event, i::Int)
  println("At time $(now(sim)) event $i is processed")
end

function append_cb(sim::Simulation, ev::Event)
  try
    append_callback(ev, append_cb)
  catch exc
    println(exc)
  end
  try
    schedule(sim, ev)
  catch exc
    println(exc)
  end
  schedule!(sim, ev, 3.5)
end

sim = Simulation(2)
for i = 1:10
  ev = timeout(sim, rand())
  append_callback(ev, print_cb, i)
end
timeout(sim, 3)
append_callback(ev, append_cb)
try
  run(sim, 4)
catch exc
  println(exc)
end

sim = Simulation()
try
  run(sim, Event())
catch(exc)
  println(exc)
end


text:
# The Bio.Var module
# ==================
#
# Types and methods for analysing biological variation.
#
# Part of the Bio.Var module.
#
# This file is a part of BioJulia. License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

module Var
using Bio.Seq

export

    # Mutation types
    MutationType,
    AnyMutation,
    TransitionMutation,
    TransversionMutation,

    # Identifying and counting mutations
    count_mutations,
    is_mutation,
    flagmutations,

    # Genetic and Evolutionary distances
    EvolutionaryDistance,
    Count,
    Proportion,
    JukesCantor69,
    Kimura80,

    distance

include("mutation_counting.jl")
include("distances.jl")

end # module Var


text:
@ghdef mutable struct Reference
    ref::Union{String, Nothing}
    url::Union{HTTP.URI, Nothing}
    object::Union{Dict, Nothing}
end

name(ref::Reference) = String(split(ref.ref, "refs/")[2])

@api_default function reference(api::GitHubAPI, repo, ref_obj; options...)
    result = gh_get_json(api, "/repos/$(name(repo))/git/refs/$(name(ref_obj))"; options...)
    return Reference(result)
end

@api_default function references(api::GitHubAPI, repo; options...)
    results, page_data = gh_get_paged_json(api, "/repos/$(name(repo))/git/refs"; options...)
    return Reference.(results), page_data
end

@api_default function create_reference(api::GitHubAPI, repo; options...)
    result = gh_post_json(api, "/repos/$(name(repo))/git/refs"; options...)
    return Reference(result)
end

@api_default function update_reference(api::GitHubAPI, repo, ref_obj; options...)
    result = gh_patch_json(api, "/repos/$(name(repo))/git/refs/$(name(ref_obj))"; options...)
    return Reference(result)
end

@api_default function delete_reference(api::GitHubAPI, repo, ref_obj; options...)
    return gh_delete(api, "/repos/$(name(repo))/git/refs/$(name(ref_obj))"; options...)
end


text:
using Random
using Distributions
using MultiAgentPathFinding
using MultiAgentAllocationTransit
using BenchmarkTools
using Statistics
using JSON
using Logging
global_logger(SimpleLogger(stderr, Logging.Warn))

rng = MersenneTwister(2345)

# Script arguments relating to transit files etc.
const city_params_file = ARGS[1]
const stop_coords_file = ARGS[2]
const trips_file = ARGS[3]
const drone_params_file = ARGS[4]
const out_file = ARGS[5]

const TRANSIT_CAP_RANGE = (3,5)

tg = load_transit_graph_latlong(stop_coords_file, trips_file, TRANSIT_CAP_RANGE, rng)
drone_params = parse_drone_params(drone_params_file)

t = @benchmark transit_graph_preprocessing($tg, distance_lat_lon_euclidean, $drone_params)
pp_dict = Dict("city"=>city_params_file,
               "mean"=>mean(t.times),
               "median"=>median(t.times))

open(out_file, "w") do f
    JSON.print(f, pp_dict, 2)
end


text:
export nn_constitutive_law
"""
    nn_constitutive_law(input::Array{Float64,2}, θ::Array{Float64,1}, grad_input::Int64, grad_θ::Int64)

Computes the stress with the given neural network weights θ.
# Inputs
`input` is a ``n\\times 9`` array, θ is a `p` dimensional array
`grad_input = 0`: do not compute gradients with respect to input; otherwise compute the gradients, in this case, `g` must have the same size as σ
`grad_θ = 0`    : do not compute gradients with respect to θ; otherwise compute the gradients

output:
- σ
- ``\\frac{\\partial \\sigma}{\\partial ipt}``
- ``\\frac{\\partial \\sigma\\cdot g}{\\partial \\theta}``
"""
function nn_constitutive_law(input::Array{Float64,2}, θ::Array{Float64,1},
    g::Union{Array{Float64,2}, Nothing}=nothing, grad_input::Bool=false, grad_θ::Bool=false)
    nn_constitutive_law(input, θ, [9,20,20,4], g, grad_input, grad_θ)
end

function nn_constitutive_law(input::Array{Float64,2}, θ::Array{Float64,1},
        config::Array{Int64},
        g::Union{Array{Float64,2}, Nothing}=nothing, grad_input::Bool=false, grad_θ::Bool=false)

    @assert config[1]==9 && config[end]==4
    n_layer = length(config)
    nθ = 0
    for i = 1:n_layer-1
        nθ += config[i]*config[i+1]
        nθ += config[i+1]
    end
    @assert length(θ)==nθ

    n = size(input,1)
    @assert size(input,2)==9

    σ = zeros(n*3)
    input = input'[:]

    dinput = zeros()
    if(grad_input)
        dinput = zeros(n*9*3)
        grad_input = 1
    else
        grad_input = 0
    end

    dθ = zeros()
    if(grad_θ)
        @assert size(g,1)==n && size(g,2)==3
        g = g'[:]
        dθ = zeros(length(θ))
        grad_θ = 1
    else
        g = zeros()
        grad_θ = 0
    end

    lib = joinpath(@__DIR__, "../../deps/CustomOp/ADLaw/build/libnnlaw")
    @eval begin
        ccall((:constitutive_law_generic, $lib), Cvoid,
            (Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},
            Ptr{Cint}, Cint,
            Cint, Cint, Cint), $σ, $input, $θ, $g, $dinput, $dθ,
            $(Int32.(config)), $n_layer,
            $n, $grad_input, $grad_θ)
    end

    σ = reshape(σ, 3, n)'|>Array
    if grad_input>0
        new_input = zeros(n, 9, 3)
        for i = 1:3
            o = dinput[(i-1)*9*n+1:i*9*n]
            new_input[:,:,i] = reshape(o, 9, n)'
        end
        dinput = new_input
    end

    return σ, dinput, dθ
end


text:
export transport_mcmc

"""
    traces, lml_est = transport_mcmc(
      trace,
      kernel,
      family::TransportFamily,
      init_param::AbstractArray{<:AbstractArray{<:Real}},
      num_samples::Int,
      adapt_delay::Int
    )

"""
function transport_mcmc(
  trace,
  kernel,
  family::TransportFamily,
  init_param::AbstractArray{<:AbstractArray{<:Real}},
  num_samples::Int,
  adapt_delay::Int
)
end

max_stars_repo_name: mvarble/gen-transport-mcmc

max_issues_repo_name: mvarble/gen-transport-mcmc

max_forks_repo_name: mvarble/gen-transport-mcmc


text:
nz=43; nx=35
ny=randn(nz, nx);

my1=J.Fdtd.get_rhovxI(ny)
my2=J.Fdtd.get_rhovzI(ny)

myp1 = randn(size(my1));
myp2 = randn(size(my2));

# spray
nyp=zeros(nz, nx);

J.Fdtd.grad_modrr_sprayrr!(nyp, myp1, myp2)
J.Fdtd.grad_modrr_sprayrrvx!(nyp, myp1)
J.Fdtd.grad_modrr_sprayrrvz!(nyp, myp2)

# dot product test
@test LinearAlgebra.dot(vcat(my1,my2), vcat(myp1,myp2)) ≈ LinearAlgebra.dot(ny, nyp)


text:
using LinearAlgebra

"""
A kernel function takes two vectors and calculates
the distance between them.
"""
abstract type Kernel end

"""
    SquaredExponentialKernel(l, σf, σy)

Squared exponential kernel for noisy observations.
The parameter names follow (15.20) in MLAPP.

# Parameters
- `l`: length scale, which controls the horizontal scale.
      It can be in either of the forms
   1. `Float64`: each dimension will have the same length scale
   2. `Vector{Float64}`: length scale for each dimension.
- `σf`: vertical scale
- `σy`: observation noise
"""
mutable struct SquaredExponentialKernel <: Kernel
    l::Float64
    σf::Float64
    σy::Float64
end

function (kernel::SquaredExponentialKernel)(xp, xq, δpq)
    σf = kernel.σf
    σy = kernel.σy
    M = lengthscalematrix(kernel.l, length(xp))
    return σf^2 * exp(-0.5 * (xp - xq)' * M * (xp - xq)) + σy^2*δpq
end

lengthscalematrix(l::Float64, dimensionsize) = Diagonal(ones(dimensionsize)) / l^2

"""
    MixtureKernel(X, p, q)

A kernel mixture of Gaussian kernel and dot product.
The definition follows (6.63) in PRML.
"""
mutable struct MixtureKernel <: Kernel
    sekernel::SquaredExponentialKernel
    θ0::Float64
    θ1::Float64

    function MixtureKernel(l, σf, σy, θ0, θ1)
        sekernel = SquaredExponentialKernel(l, σf, σy)
        new(sekernel, θ0, θ1)
    end

end

function (kernel::MixtureKernel)(xp, xq, δpq)
    sekernel = kernel.sekernel
    θ0 = kernel.θ0
    θ1 = kernel.θ1
    return sekernel(xp, xq, δpq) + θ0 + θ1 * (xp' * xq)
end

function (kernel::Kernel)(inputtuple::Tuple)
    kernel(inputtuple...)
end

"""
    covarianceblock(kernel, X, Xs)

Construct covariance matrix between `X` and `Xs` based on
kernel function.

# Parameters

- `kernel`: kernel function.
- `X`: input matrix ``X``.
- `Xs`: another input matrix ``X_*`` (X-star).
- `addnoise`: if `true`, when building a covariance matrix, we add
      observation noise to the diagonal term.
"""
function covarianceblock(kernel, X, Xs; addnoise=false)
    N = size(X, 1)
    Ns = size(Xs, 1)
    K = zeros(N, Ns)
    for p in 1:N, q in 1:Ns
        xp = X[p, :]
        xq = Xs[q, :]
        if addnoise
            K[p, q] = kernel(xp, xq, p==q)
        else
            K[p, q] = kernel(xp, xq, false)
        end
    end
    K
end

function covariancematrix(kernel, X; addnoise=true)
    covarianceblock(kernel, X, X, addnoise=addnoise)
end

"""
Because inside squared exponential kernel we need Euclidean distance,
we can use the shortcut for it following (14.30). This can significantly
reduce the computational cost.
"""
function covariancematrix(kernel::SquaredExponentialKernel, X)
    l, σf, σy = kernel.l, kernel.σf, kernel.σy
    N = size(X, 1)
    euclideandist = X.^2 .+ (X.^2)' .- 2*X*X'
    return σf^2 * exp.(-0.5/l^2 .* euclideandist) + Diagonal(fill(σy^2, N))
end

max_stars_repo_name: Jetafull/AlgorithmsInJulia

max_issues_repo_name: Jetafull/algorithms-in-julia

max_forks_repo_name: Jetafull/algorithms-in-julia


text:
include("util.jl")
using ImageFiltering

function fuck(; debug_level=verbose)
       printd("fuckyou", verbose, debug_level)
       printd("fuckyou", concise, debug_level)
       printd("fuckyou", quiet, debug_level)

end
"""
Replaces by value e (by default machine precision eps())
"""
function zero_protection(a::Array{Float64,2}, e::Float64 = eps())
	(h,w) = size(a)
    #println(" e = $e")
	for j in 1:w
		for i in 1:h
			if @inbounds a[i,j] < e
				@inbounds a[i,j] = e
			end
		end
	end
end

function richardson_lucy(image::Array{Float64,2}, psf::AbstractArray;
                         iterations::Int=10,
                         view_window::Int=1,
                         fom_window::Int=64,
                         debug_level=mute)

    # init
    image_ref  = copy(image)
    psf_mirror = flip2d(psf)
    im_deconv  = 0.5 * ones(size(image))

    printd("view of image center = $(center_view(image_ref, view_window))",
           verbose, debug_level)
    printd("view of image deconv = $(center_view(im_deconv, view_window))",
           verbose, debug_level)
    printd("view of psf center   = $(center_view(collect(psf), view_window))",
           verbose, debug_level)
    printd("view of psf* deconv  = $(center_view(collect(psf_mirror), view_window))",
           verbose, debug_level)

    dff0 = 1.
    DFR = zeros(iterations) * 1.  # Relative differences

    for i in 1:iterations
        # Using imfilter with reflect to produce a convolution (rather than a correlation),
        # psf must be centered so that final image has the same size than original image
        x = imfilter(im_deconv, reflect(centered(psf)), Fill(0) )

        printd(" iteration = $i x ($i) view = $(center_view(x, view_window)) ",
               verbose, debug_level)

        # protect against division by zero
        limit_boundaries(x, eps())

        # relative blur
        relative_blur = image ./ x

        printd(" relBlur($i) view = $(center_view(relative_blur, view_window)) ",
               verbose, debug_level)

        # deconvolution
        im_deconv = im_deconv .* imfilter(relative_blur, reflect(centered(psf_mirror)), Fill(0))

        printd("im_deconv($i) view = $(center_view(relative_blur, view_window)) ",
               verbose, debug_level)

        # Relative differences
        dff = image_fom(collect(center_view(image_ref, 64)), collect(center_view(im_deconv, 64)))
        delta = (dff0 - dff) / dff0

        printd(" iteration $i ", concise, debug_level)
        printd(" diff($i)  = $dff ", concise, debug_level)
        printd(" delta($i)  = $delta ", concise, debug_level)

        DFR[i] = delta

        dff0 = dff
        image_ref = copy(im_deconv)
    end

    # clip image ( 0 < x <= 1)
    limit_boundaries(im_deconv, 0.,1.)
    return im_deconv, DFR

end

function frichardson_lucy(image::Array{Float64,2}, psf::AbstractArray; iterations::Int=10)

    println("Running RL with $iterations iterations")
    psf_mirror = flip2d(psf)
    im_deconv  = 0.5 * ones(size(image))

    for i in 1:iterations
        # Using imfilter with reflect to produce a convolution (rather than a correlation),
        # psf must be centered so that final image has the same size than original image
        x = imfilter(im_deconv, reflect(centered(psf)), Fill(0))

        # protect against division by zero
        #zero_protection(x)

        # relative blur
        relative_blur = image ./ x

        # deconvolution
        im_deconv = im_deconv .* imfilter(relative_blur,
                                          reflect(centered(psf_mirror)), Fill(0))
    end

    # clip image ( 0 < x <= 1)
    limit_boundaries(im_deconv, 0.,1.)
    return im_deconv

end

max_stars_repo_name: jjgomezcadenas/JImages

max_issues_repo_name: jjgomezcadenas/JImages

max_forks_repo_name: jjgomezcadenas/JImages