text:
@node MvNormalMeanCovariance Stochastic [ out, (μ, aliases = [ mean ]), (Σ, aliases = [ cov ]) ]

conjugate_type(::Type{ <: MvNormalMeanCovariance }, ::Type{ Val{ :out } }) = MvNormalMeanCovariance
conjugate_type(::Type{ <: MvNormalMeanCovariance }, ::Type{ Val{ :μ } })   = MvNormalMeanCovariance
conjugate_type(::Type{ <: MvNormalMeanCovariance }, ::Type{ Val{ :Σ } })   = InverseWishart

@average_energy MvNormalMeanCovariance (q_out::Any, q_μ::Any, q_Σ::Any) = begin
    (m_mean, v_mean) = mean(q_μ), cov(q_μ)
    (m_out, v_out)   = mean(q_out), cov(q_out)
    0.5 * (ndims(q_out) * log2π + logdet(mean(q_Σ)) + tr(cholinv(mean(q_Σ))*(v_out + v_mean + (m_out - m_mean)*(m_out - m_mean)')))
end

@average_energy MvNormalMeanCovariance (q_out_μ::Any, q_Σ::Any) = begin
    (m, V) = mean(q_out_μ), cov(q_out_μ)
    d = Int64(ndims(q_out_μ)/2)
    @views 0.5*(d*log2π + logdet(mean(q_Σ)) + tr(cholinv(mean(q_Σ))*( V[1:d,1:d] - V[1:d,d+1:end] - V[d+1:end,1:d] + V[d+1:end,d+1:end] + (m[1:d] - m[d+1:end])*(m[1:d] - m[d+1:end])' ) ))
end


text:
using CSV, DataFrames, DecisionTree, Distributions, Gurobi, JuMP, LightGraphs, Parameters, SparseArrays, Statistics, ArgParse
include("../oracles/shortest_path_oracle.jl")
include("../solver/util.jl")
include("../solver/sgd.jl")
include("../solver/reformulation.jl")
include("../solver/random_forests_po.jl")
include("../solver/validation_set.jl")
include("../experiments/replication_functions.jl")

# Fixed parameter settings (these never change)
p_features = 5
grid_dim = 5
num_trials = 50
n_test = 10000

num_lambda = 1
lambda_max = 0
lambda_min_ratio = 1
holdout_percent = 0.25
regularization = :lasso
different_validation_losses = false
include_rf = true

# Fixed parameter sets (these are also the same for all experiments)
n_train_vec = [5000]
polykernel_degree_vec = [1; 2; 4; 6; 8]
polykernel_noise_half_width_vec = [0; 0.5]

# Set this to get reproducible results
rng_seed = 64234

# Run experiment and get results
# Note that Gurobi enviornments are set within this function call
expt_results = shortest_path_multiple_replications(rng_seed, num_trials, grid_dim,
    n_train_vec, n_test,
    p_features, polykernel_degree_vec, polykernel_noise_half_width_vec;
    num_lambda = num_lambda, lambda_max = lambda_max, lambda_min_ratio = lambda_min_ratio,
    holdout_percent = holdout_percent, regularization = regularization,
    different_validation_losses = different_validation_losses,
    include_rf = include_rf)

csv_string = "shortest_path_5000.csv"
CSV.write(csv_string, expt_results)

ext: jl

max_stars_repo_name: paulgrigas/SmartPredictThenOptimize

max_issues_repo_name: paulgrigas/SmartPredictThenOptimize

max_forks_repo_name: paulgrigas/SmartPredictThenOptimize


text:
export SimPosteriorSampler, SimIexInsertDelete, SimIexInsertDeleteProportional
export SimIexInsertDeleteDependent, SimIexInsertDeleteWithKick
export SimIexSplitMerge

abstract type SimPosteriorSampler end

struct SimIexInsertDelete <: SimPosteriorSampler
    ν_ed::Int  # Maximum number of edit operations
    ν_td::Int  # Maximum increase/decrease in dimension
    β::Real  # Probability of trans-dimensional move
    α::Float64  # How much vertex proposal is informed by data (0.0 max and ∞ is uniform over vertice)
    len_dist::DiscreteUnivariateDistribution  # Distribution used to introduce new interactions
    ε::Float64
    aux_mcmc::SimMcmcSampler
    K::Int # Max number of interactions (used to determined how many pointers to store interactions)
    desired_samples::Int  # Final three set default values for MCMC samplers
    burn_in::Int
    lag::Int
    par_info::Dict
    curr_pointers::InteractionSequence{Int} # Storage for prev value in MCMC
    prop_pointers::InteractionSequence{Int} # Storage for curr value in MCMC
    ind_del::Vector{Int} # Storage for indexing of deletions of interactions
    ind_add::Vector{Int} # Storage for indexing of additions of interactions
    vals::Vector{Int} # Storage for values where are deleted from interactions
    ind_update::Vector{Int} # Storage of which values have been updated
    ind_td_add::Vector{Int} # Storage of where to insert/delete
    ind_td_del::Vector{Int} # Storage of where to insert/delete
    function SimIexInsertDelete(
        aux_mcmc::SimMcmcSampler;
        K=100,
        ν_ed=2, ν_td=1 , β=0.7, α=0.0, len_dist=TrGeometric(0.9,1,K),
        ε=0.1,
        desired_samples=1000, lag=1, burn_in=0
        )
        curr_pointers = [Int[] for i in 1:K]
        prop_pointers = [Int[] for i in 1:K]
        ind_del = zeros(Int, ν_ed)
        ind_add = zeros(Int, ν_ed)
        vals = zeros(Int, ν_ed)
        ind_update = zeros(Int, ν_ed)
        ind_td_add = zeros(Int, ν_td)
        ind_td_del = zeros(Int, ν_td)
        par_info = Dict()
        par_info[:ν_ed] = "(maximum number of edit operations)"
        par_info[:ν_td] = "(maximum change in dimension)"
        par_info[:len_dist] = "(distribution to sample length of path insertions)"
        par_info[:α] = "(controls how much data informs entry insertion proposals)"
        par_info[:aux_mcmc] = "(mcmc sampler from auxiliary data)"
        par_info[:β] = "(probability of update move)"
        par_info[:K] = "(maximum number of interactions, used to initialise storage)"
        par_info[:ε] = "(neighborhood for γ proposals)"

        new(
            ν_ed, ν_td, β, α,
            len_dist,
            ε,
            aux_mcmc, K,
            desired_samples, burn_in, lag,
            par_info,
            curr_pointers, prop_pointers, ind_del, ind_add, vals,
            ind_update,
            ind_td_add, ind_td_del
        )
    end
end

struct SimIexSplitMerge <: SimPosteriorSampler
    ν_ed::Int  # Maximum number of edit operations
    ν_td::Int  # Maximum increase/decrease in dimension
    β::Real  # Probability of trans-dimensional move
    α::Float64  # How much vertex proposal is informed by data (0.0 max and ∞ is uniform over vertice)
    len_dist::DiscreteUnivariateDistribution  # Distribution used to introduce new interactions
    ε::Float64
    aux_mcmc::SimMcmcSampler
    K::Int # Max number of interactions (used to determined how many pointers to store interactions)
    desired_samples::Int  # Final three set default values for MCMC samplers
    burn_in::Int
    lag::Int
    par_info::Dict
    curr_pointers::InteractionSequence{Int} # Storage for prev value in MCMC
    prop_pointers::InteractionSequence{Int} # Storage for curr value in MCMC
    store::InteractionSequence{Int} # Temporary storage between for use in move
    ind_del::Vector{Int} # Storage for indexing of deletions of interactions
    ind_add::Vector{Int} # Storage for indexing of additions of interactions
    vals::Vector{Int} # Storage for values where are deleted from interactions
    ind_update::Vector{Int} # Storage of which values have been updated
    ind_td_split::Vector{Int} # Storage for split/merge locations
    ind_td_merge::Vector{Int} # Storage for split/merge locations
    function SimIexSplitMerge(
        aux_mcmc::SimMcmcSampler;
        K=100,
        ν_ed=2, ν_td=1 , β=0.7, α=0.0, len_dist=TrGeometric(0.9,1,K),
        ε=0.1,
        desired_samples=1000, lag=1, burn_in=0
        )
        curr_pointers = [Int[] for i in 1:K]
        prop_pointers = [Int[] for i in 1:K]
        store = Vector{Int}[]
        ind_del = zeros(Int, ν_ed)
        ind_add = zeros(Int, ν_ed)
        vals = zeros(Int, ν_ed)
        ind_update = zeros(Int, ν_ed)
        ind_td_split = zeros(Int, ν_td)
        ind_td_merge = zeros(Int, 2*ν_td) # Need to index double
        par_info = Dict()
        par_info[:ν_ed] = "(maximum number of edit operations)"
        par_info[:ν_td] = "(maximum change in dimension)"
        par_info[:len_dist] = "(distribution to sample length of path insertions)"
        par_info[:α] = "(controls how much data informs entry insertion proposals)"
        par_info[:aux_mcmc] = "(mcmc sampler from auxiliary data)"
        par_info[:β] = "(probability of update move)"
        par_info[:K] = "(maximum number of interactions, used to initialise storage)"
        par_info[:ε] = "(neighborhood for γ proposals)"

        new(
            ν_ed, ν_td, β, α,
            len_dist,
            ε,
            aux_mcmc, K,
            desired_samples, burn_in, lag,
            par_info,
            curr_pointers, prop_pointers, store,
            ind_del, ind_add, vals,
            ind_update,
            ind_td_split, ind_td_merge
        )
    end
end

struct SimIexInsertDeleteGibbs <: SimPosteriorSampler
    ν_ed::Int  # Maximum number of edit operations
    ν_td::Int  # Maximum increase/decrease in dimension
    β::Real  # Probability of trans-dimensional move
    α::Float64  # How much vertex proposal is informed by data (0.0 max and ∞ is uniform over vertice)
    len_dist::DiscreteUnivariateDistribution  # Distribution used to introduce new interactions
    ε::Float64
    aux_mcmc::SimMcmcSampler
    K::Int # Max number of interactions (used to determined how many pointers to store interactions)
    desired_samples::Int  # Final three set default values for MCMC samplers
    burn_in::Int
    lag::Int
    par_info::Dict
    count::Vector{Int}      # Counts for attempted moves
    acc_count::Vector{Int}  # Counts for accepted moves
    curr_pointers::InteractionSequence{Int} # Storage for prev value in MCMC
    prop_pointers::InteractionSequence{Int} # Storage for curr value in MCMC
    ind_del::Vector{Int} # Storage for indexing of deletions of interactions
    ind_add::Vector{Int} # Storage for indexing of additions of interactions
    vals::Vector{Int} # Storage for values where are deleted from interactions
    ind_update::Vector{Int} # Storage of which values have been updated
    ind_td_add::Vector{Int} # Storage of where to insert/delete
    ind_td_del::Vector{Int} # Storage of where to insert/delete
    function SimIexInsertDeleteGibbs(
        aux_mcmc::SimMcmcSampler;
        K=100,
        ν_ed=2, ν_td=1 , β=0.7, α=0.0, len_dist=TrGeometric(0.9,1,K),
        ε=0.1,
        desired_samples=1000, lag=1, burn_in=0
        )
        count, acc_count = (zeros(Int, 3), zeros(Int, 3))
        curr_pointers = [Int[] for i in 1:K]
        prop_pointers = [Int[] for i in 1:K]
        ind_del = zeros(Int, ν_ed)
        ind_add = zeros(Int, ν_ed)
        vals = zeros(Int, ν_ed)
        ind_update = zeros(Int, ν_ed)
        ind_td_add = zeros(Int, ν_td)
        ind_td_del = zeros(Int, ν_td)
        par_info = Dict()
        par_info[:ν_ed] = "(maximum number of edit operations)"
        par_info[:ν_td] = "(maximum change in dimension)"
        par_info[:len_dist] = "(distribution to sample length of path insertions)"
        par_info[:α] = "(controls how much data informs entry insertion proposals)"
        par_info[:aux_mcmc] = "(mcmc sampler from auxiliary data)"
        par_info[:β] = "(probability of update move)"
        par_info[:K] = "(maximum number of interactions, used to initialise storage)"
        par_info[:ε] = "(neighborhood for γ proposals)"

        new(
            ν_ed, ν_td, β, α,
            len_dist,
            ε,
            aux_mcmc, K,
            desired_samples, burn_in, lag,
            par_info,
            count, acc_count,
            curr_pointers, prop_pointers,
            ind_del, ind_add, vals,
            ind_update,
            ind_td_add, ind_td_del
        )
    end
end

get_summary_measures(mcmc::SimMcmcInsertDeleteGibbs) = Dict(
    "Insert/Delete Acc. Prob." => mcmc.acc_count[1]/mcmc.count[1],
    "Flip Acc. Prob." => mcmc.acc_count[2]/mcmc.count[2],
    "Trans-dim Acc. Prob." => mcmc.acc_count[3]/mcmc.count[3]
)

struct SimIexInsertDeleteProportional <: SimPosteriorSampler
    ν_ed::Int  # Maximum number of edit operations
    ν_td::Int  # Maximum increase/decrease in dimension
    β::Real  # Probability of trans-dimensional move
    α::Float64  # How much vertex proposal is informed by data (0.0 max and ∞ is uniform over vertice)
    len_dist::DiscreteUnivariateDistribution  # Distribution used to introduce new interactions
    ε::Float64
    aux_mcmc::SimMcmcSampler
    K::Int # Max number of interactions (used to determined how many pointers to store interactions)
    desired_samples::Int  # Final three set default values for MCMC samplers
    burn_in::Int
    lag::Int
    par_info::Dict
    curr_pointers::InteractionSequence{Int} # Storage for prev value in MCMC
    prop_pointers::InteractionSequence{Int} # Storage for curr value in MCMC
    ind_del::Vector{Int} # Storage for indexing of deletions of interactions
    ind_add::Vector{Int} # Storage for indexing of additions of interactions
    vals::Vector{Int} # Storage for values where are deleted from interactions
    ind_update::Vector{Int} # Storage of which values have been updated
    ind_td_add::Vector{Int} # Storage of where to insert/delete
    ind_td_del::Vector{Int} # Storage of where to insert/delete
    function SimIexInsertDeleteProportional(
        aux_mcmc::SimMcmcSampler;
        K=100,
        ν_ed=2, ν_td=1 , β=0.7, α=0.0, len_dist=TrGeometric(0.9,1,K),
        ε=0.1,
        desired_samples=1000, lag=1, burn_in=0
        )
        curr_pointers = [Int[] for i in 1:K]
        prop_pointers = [Int[] for i in 1:K]
        ind_del = zeros(Int, ν_ed)
        ind_add = zeros(Int, ν_ed)
        vals = zeros(Int, ν_ed)
        ind_update = zeros(Int, ν_ed)
        ind_td_add = zeros(Int, ν_td)
        ind_td_del = zeros(Int, ν_td)
        par_info = Dict()
        par_info[:ν_ed] = "(maximum number of edit operations)"
        par_info[:ν_td] = "(maximum change in dimension)"
        par_info[:len_dist] = "(distribution to sample length of path insertions)"
        par_info[:α] = "(controls how much data informs entry insertion proposals)"
        par_info[:aux_mcmc] = "(mcmc sampler from auxiliary data)"
        par_info[:β] = "(probability of update move)"
        par_info[:K] = "(maximum number of interactions, used to initialise storage)"
        par_info[:ε] = "(neighborhood for γ proposals)"

        new(
            ν_ed, ν_td, β, α,
            len_dist,
            ε,
            aux_mcmc, K,
            desired_samples, burn_in, lag,
            par_info,
            curr_pointers, prop_pointers, ind_del, ind_add, vals,
            ind_update,
            ind_td_add, ind_td_del
            )
    end
end

struct SimIexInsertDeleteDependent <: SimPosteriorSampler
    ν_ed::Int  # Maximum number of edit operations
    ν_td::Int  # Maximum increase/decrease in dimension
    β::Real  # Probability of trans-dimensional move
    α::Float64  # How much vertex proposal is informed by data (0.0 max and ∞ is uniform over vertice)
    len_dist::DiscreteUnivariateDistribution  # Distribution used to introduce new interactions
    ε_ed::Float64
    ε_td::Float64
    aux_mcmc::SimMcmcSampler
    K::Int # Max number of interactions (used to determined how many pointers to store interactions)
    desired_samples::Int  # Final three set default values for MCMC samplers
    burn_in::Int
    lag::Int
    par_info::Dict
    curr_pointers::InteractionSequence{Int} # Storage for prev value in MCMC
    prop_pointers::InteractionSequence{Int} # Storage for curr value in MCMC
    ind_del::Vector{Int} # Storage for indexing of deletions of interactions
    ind_add::Vector{Int} # Storage for indexing of additions of interactions
    vals::Vector{Int} # Storage for values where are deleted from interactions
    ind_update::Vector{Int} # Storage of which values have been updated
    ind_td::Vector{Int} # Storage of where to insert/delete
    function SimIexInsertDeleteDependent(
        aux_mcmc::SimMcmcSampler;
        K=100,
        ν_ed=2, ν_td=1 , β=0.7, α=0.0, len_dist=TrGeometric(0.9,1,K),
        ε_ed=0.1, ε_td=2.0,
        desired_samples=1000, lag=1, burn_in=0
        )
        curr_pointers = [Int[] for i in 1:K]
        prop_pointers = [Int[] for i in 1:K]
        ind_del = zeros(Int, ν_ed)
        ind_add = zeros(Int, ν_ed)
        vals = zeros(Int, ν_ed)
        ind_update = zeros(Int, ν_ed)
        ind_td = zeros(Int, ν_td)
        par_info = Dict()
        par_info[:ν_ed] = "(maximum number of edit operations)"
        par_info[:ν_td] = "(maximum change in dimension)"
        par_info[:len_dist] = "(distribution to sample length of path insertions)"
        par_info[:α] = "(controls how much data informs entry insertion proposals)"
        par_info[:aux_mcmc] = "(mcmc sampler from auxiliary data)"
        par_info[:β] = "(probability of update move)"
        par_info[:K] = "(maximum number of interactions, used to initialise storage)"
        par_info[:ε] = "(neighborhood for γ proposals)"

        new(
            ν_ed, ν_td, β, α,
            len_dist,
            ε_ed, ε_td,
            aux_mcmc, K,
            desired_samples, burn_in, lag,
            par_info,
            curr_pointers, prop_pointers, ind_del, ind_add, vals,
            ind_update, ind_td
            )
    end
end

struct SimIexInsertDeleteWithKick <: SimPosteriorSampler
    ν_ed::Int  # Maximum number of edit operations
    ν_td::Int  # Maximum increase/decrease in dimension
    β::Real  # Probability of trans-dimensional move
    α::Float64  # How much vertex proposal is informed by data (0.0 max and ∞ is uniform over vertice)
    len_dist::DiscreteUnivariateDistribution  # Distribution used to introduce new interactions
    ε_ed::Float64
    ε_td::Float64
    ε_kick::Float64
    aux_mcmc::SimMcmcSampler
    K::Int # Max number of interactions (used to determined how many pointers to store interactions)
    desired_samples::Int  # Final three set default values for MCMC samplers
    burn_in::Int
    lag::Int
    par_info::Dict
    curr_pointers::InteractionSequence{Int} # Storage for prev value in MCMC
    prop_pointers::InteractionSequence{Int} # Storage for curr value in MCMC
    ind_del::Vector{Int} # Storage for indexing of deletions of interactions
    ind_add::Vector{Int} # Storage for indexing of additions of interactions
    vals::Vector{Int} # Storage for values where are deleted from interactions
    ind_update::Vector{Int} # Storage of which values have been updated
    ind_td::Vector{Int} # Storage of where to insert/delete
    function SimIexInsertDeleteWithKick(
        aux_mcmc::SimMcmcSampler;
        K=100,
        ν_ed=2, ν_td=1 , β=0.7, α=0.0, len_dist=TrGeometric(0.9,1,K),
        ε_ed=0.1, ε_td=1.0, ε_kick=ε_td/4,
        desired_samples=1000, lag=1, burn_in=0
        )
        curr_pointers = [Int[] for i in 1:K]
        prop_pointers = [Int[] for i in 1:K]
        ind_del = zeros(Int, ν_ed)
        ind_add = zeros(Int, ν_ed)
        vals = zeros(Int, ν_ed)
        ind_update = zeros(Int, ν_ed)
        ind_td = zeros(Int, ν_td)
        par_info = Dict()
        par_info[:ν_ed] = "(maximum number of edit operations)"
        par_info[:ν_td] = "(maximum change in dimension)"
        par_info[:len_dist] = "(distribution to sample length of path insertions)"
        par_info[:α] = "(controls how much data informs entry insertion proposals)"
        par_info[:aux_mcmc] = "(mcmc sampler from auxiliary data)"
        par_info[:β] = "(probability of update move)"
        par_info[:K] = "(maximum number of interactions, used to initialise storage)"
        par_info[:ε] = "(neighborhood for γ proposals)"

        new(
            ν_ed, ν_td, β, α,
            len_dist,
            ε_ed, ε_td, ε_kick,
            aux_mcmc, K,
            desired_samples, burn_in, lag,
            par_info,
            curr_pointers, prop_pointers, ind_del, ind_add, vals,
            ind_update, ind_td
            )
    end
end


text:
mutable struct SDEIntegrator{algType,IIP,uType,uEltype,tType,tdirType,P2,eigenType,tTypeNoUnits,uEltypeNoUnits,randType,randType2,rateType,solType,cacheType,F4,F5,F6,OType,noiseType,EventErrorType,CallbackCacheType,RCs} <: AbstractSDEIntegrator{algType,IIP,uType,tType}
  f::F4
  g::F5
  c::F6
  noise::noiseType
  uprev::uType
  tprev::tType
  t::tType
  u::uType
  p::P2
  dt::tType
  dtnew::tType
  dtpropose::tType
  dtcache::tType
  T::tType
  tdir::tdirType
  just_hit_tstop::Bool
  isout::Bool
  event_last_time::Int
  vector_event_last_time::Int
  last_event_error::EventErrorType
  accept_step::Bool
  last_stepfail::Bool
  force_stepfail::Bool
  dtchangeable::Bool
  u_modified::Bool
  saveiter::Int
  alg::algType
  sol::solType
  cache::cacheType
  callback_cache::CallbackCacheType
  sqdt::tType
  W::randType
  P::randType2
  rate_constants::RCs
  opts::OType
  iter::Int
  success_iter::Int
  eigen_est::eigenType
  EEst::tTypeNoUnits
  q::tTypeNoUnits
  qold::tTypeNoUnits
  q11::tTypeNoUnits
  destats::DiffEqBase.DEStats
end


text:
export IterativeSolver,getIterativeSolver,solveLinearSystem!

"""
mutable struct jInv.LinearSolvers.IterativeSolver <: AbstractSolver

Fields:

	IterMethod - iterative method to apply
	PC      - symbol, (:ssor, :jac,...)
	maxIter - maximum number of iterations
	tol     - tolerance
	Ainv    - preconditioner
	out     - flag for output
	doClear - flag for deleting preconditioner
	nthreads - number of threads for spmatvecs
	sym      - 0=unsymmetric, 1=symm. pos def, 2=general symmetric
	isTranspose - if true, transpose(A) is provided to solver, else A is provided to solver
		      default=false, use isTranspose=true for efficiency with caution
		      note that A_mul_B! is slower than Ac_mul_B for SparseMatrixCSC

Example:

	getIterativeSolver(cg)

"""
mutable struct IterativeSolver<: AbstractSolver
	IterMethod::Function
	PC::Symbol
	maxIter::Int
	tol::Real
	Ainv
	out::Int
	doClear::Bool
	nthreads::Int
	sym::Int
	isTranspose::Bool
	nIter::Int
	nBuildPC::Int
	timePC::Real
	timeSolve::Real
	timeMV::Real
end

"""
function jInv.LinearSolvers.getIterativeSolver

constructs IterativeSolver

Required Input:

	IterMethod::Function   - function handle for linear solvers
		Inputs are: (A,b,M), A is matrix, b is right hand side, M is preconditioner
			Examples: IterMethod = KrylovMethods.cg   #KrylovMethods.cg already has required API
				  IterMethod(A,b;M=M,tol=1e-1,maxIter=10,out=-1) =
				                  bicgstb(A,b,M1=M,tol=tol,maxIter=maxIter,out=out)
				  IterMethod(A,b;M=M,tol=1e-1,maxIter=10,out=-1)  =
				                  gmres(A,b,5,M1=M,tol=tol,maxIter=maxIter,out=out)
			The keyword arguments of IterMethod for bicgstb and gmres
			will be initialized with the fields in the IterativeSolver type.
		Outputs are: (x,flag,err,iter), x is approximate solution

Optional Inputs:

	PC::Symbol     - specifies preconditioner, default:ssor
	maxIter        - maximum number of iterations, default:500
	tol            - tolerance on relative residual, default=1e-5
	Ainv           - preconditioner, default=identity
	out            - flag for output, default=-1 (no output)
	doClear        - flag for clearing the preconditioner, default=true
	nthreads       - number of threads to use for matvecs (requires ParSpMatVec.jl), default=4
	sym            - 0=unsymmetric, 1=symm. pos def, 2=general symmetric
    isTranspose    - if true, transpose(A) is provided to solver, else A is proved to solver
		              default=false, use isTranspose=true for efficiency with caution
		              note that A_mul_B! is slower than Ac_mul_B for SparseMatrixCSC
"""
function getIterativeSolver(IterMethod::Function;PC=:ssor,maxIter=500,tol=1e-5,
					Ainv=identity,out=-1,doClear::Bool=true,nthreads::Int=4,sym=0,isTranspose=false)
 	return IterativeSolver(IterMethod,PC,maxIter,tol,Ainv,out,doClear,nthreads,sym,isTranspose,0,0,.0,.0,.0)
end

function solveLinearSystem!(A,B,X,param::IterativeSolver,doTranspose=0)
	if param.doClear
		# clear preconditioner
		clear!(param)
		param.doClear=false
	end

	# build preconditioner
	if param.Ainv == []
		if param.PC==:ssor
			OmInvD = 1 ./Vector(diag(A));
			x      = zeros(eltype(A),size(B,1))
			M = r -> (x[:].=0.0; param.timePC += @elapsed x=ssorPrecTrans!(A,x,r,OmInvD); return x);
			param.Ainv= M
		elseif param.PC==:jac
			OmInvD = 1 ./Vector(diag(A))
			M = r -> (param.timePC += @elapsed x=r.*OmInvD; return x);
			param.Ainv= M
		else
			error("Iterativesolver: preconditioner $(param.PC) not implemented.")
		end
		param.nBuildPC+=1
	end

	# solve systems
	y     = zeros(eltype(A),size(X,1))
	doTranspose = (param.isTranspose) ? mod(doTranspose+1,2) : doTranspose
	if hasParSpMatVec
		if (param.sym==1) ||  ((param.sym != 1) && (doTranspose == 1))
			Af = x -> (y[:].=0.0; param.timeMV+=@elapsed ParSpMatVec.Ac_mul_B!(one(eltype(A)),A,x,zero(eltype(A)),y,param.nthreads);  return y)
		elseif (param.sym != 1) && (doTranspose == 0)
			Af = x -> (y[:].=0.0; param.timeMV+=@elapsed ParSpMatVec.A_mul_B!(one(eltype(A)),A,x,zero(eltype(A)),y,param.nthreads);  return y)
		end

	else
		if (param.sym==1) ||  ((param.sym != 1) && (doTranspose == 1))
			Af = x -> (y[:].=0.0; param.timeMV+=@elapsed mul!(y,adjoint(A),x,one(eltype(A)),zero(eltype(A))); return y)
		elseif (param.sym != 1) && (doTranspose == 0)
			Af = x -> (y[:].=0.0; param.timeMV+=@elapsed mul!(y,A,x,one(eltype(A)),zero(eltype(A)));  return y)
		end
	end

	t = time_ns();
	for i=1:size(X,2)
		bi      = Vector(B[:,i]);
		X[:,i],flag,err,iter = param.IterMethod(Af,bi,M = param.Ainv,tol=param.tol,maxIter=param.maxIter,out=param.out)
		param.nIter+=iter
	end
	param.timeSolve+=(time_ns()-t)/1e+9;
	return X, param
end # function solveLinearSystem PCGsolver


text:
import Pkg
Pkg.activate(@__DIR__)

module StaticRBD

using RigidBodyDynamics
using DelimitedFiles
using LinearAlgebra

Base.@ccallable function create_mechanism(c_urdf::Cstring, c_floating::Cint, c_scalar_type::Cint)::Any
    T = Dict(Cint(1) => Float64, Cint(2) => Float32)[c_scalar_type]
    urdf = unsafe_string(c_urdf)
    floating = Bool(c_floating)
    root_joint_type = floating ? QuaternionFloating{T}() : Fixed{T}()
    mechanism = parse_urdf(urdf; scalar_type=T, root_joint_type=root_joint_type)
    remove_fixed_tree_joints!(mechanism)
    mechanism.gravitational_acceleration = zero(mechanism.gravitational_acceleration)
    mechanism
end

Base.@ccallable function create_state(mechanism::Any)::Any
    MechanismState(mechanism)
end

Base.@ccallable function create_dynamics_result(mechanism::Any)::Any
    DynamicsResult(mechanism)
end

Base.@ccallable function inverse_dynamics(τ::Any, jointwrenches::Any, accelerations::Any, state::Any, v̇::Any)::Cvoid
    setdirty!(state)
    inverse_dynamics!(τ, jointwrenches, accelerations, state, v̇)
    nothing
end

Base.@ccallable function mass_matrix(M::Any, state::Any)::Cvoid
    setdirty!(state)
    mass_matrix!(M, state)
    nothing
end

Base.@ccallable function dynamics(result::Any, state::Any, τ::Any)::Cvoid
    setdirty!(state)
    dynamics!(result, state, τ)
    nothing
end

end

max_stars_repo_name: CobbledSteel/rbd-benchmarks-riscv

max_issues_repo_name: CobbledSteel/rbd-benchmarks-riscv

max_forks_repo_name: CobbledSteel/rbd-benchmarks-riscv


text:
# struct ParentsBeforeChildrenIterator
#     ops::Vector{Operation}
#     visited::Vector{Bool}
# end
# function iterate(pbci::ParentsBeforeChildrenIterator)
#     for (i,op) ∈ enumerate(pbci.ops)
#         if iszero(length(parents(op)))
#             pbci.visited[i] = true
#             return op, pbci
#         end
#     end
#     nothing
# end
# function iterate()

# end

# function dependent_outer_reducts(ls::LoopSet, op)
#     for i ∈ ls.outer_reductions
#         search_tree(parents(operations(ls)[i]), name(op)) && return true
#     end
#     false
# end

function isnopidentity(ls::LoopSet, op::Operation, u₁loop::Symbol, u₂loop::Symbol, vectorized::Symbol, u₂max::Int)
  parents_op = parents(op)
  if iscompute(op) && instruction(op).instr === :identity
    # loopistiled = u₂max ≠ -1
    # parents_u₁syms, parents_u₂syms = parent_unroll_status(op, u₁loop, u₂loop, u₂max)
    # if (u₁unrolledsym == first(parents_u₁syms)) && (isu₂unrolled(op) == parents_u₂syms[1])
    oppstate = Base.iterate(parents_op)
    oppstate === nothing && return false
    opp, state = oppstate
    Base.iterate(parents_op, state) === nothing || return false
    name(opp) === name(op) || return false
    # @show op opp isu₁unrolled(op), isu₁unrolled(opp), isu₂unrolled(op), isu₂unrolled(opp)
    (isu₁unrolled(op) == isu₁unrolled(opp)) & (isu₂unrolled(op) == isu₂unrolled(opp))
  else
    false
  end
end

function set_upstream_family!(adal::Vector{T}, op::Operation, val::T, ld::Vector{Symbol}, id::Int) where {T}
    adal[identifier(op)] == val && return # must already have been set
    if ld != loopdependencies(op) || id == identifier(op)
        (adal[identifier(op)] = val)
    end
    for opp ∈ parents(op)
        identifier(opp) == identifier(op) && continue
        set_upstream_family!(adal, opp, val, ld, id)
    end
end
function search_for_reductinit!(op::Operation, opswap::Operation, var::Symbol, loopdeps::Vector{Symbol})
  for (i,opp) ∈ enumerate(parents(op))
    if (name(opp) === var) && (length(reduceddependencies(opp)) == 0) && (length(loopdependencies(opp)) == length(loopdeps)) && (length(children(opp)) == 1)
      if all(in(loopdeps), loopdependencies(opp))
        parents(op)[i] = opswap
        return opp
      end
    end
    opcheck = search_for_reductinit!(opp, opswap, var, loopdeps)
    opcheck === opp || return opcheck
  end
  return op
end
function addoptoorder!(
    ls::LoopSet, included_vars::Vector{Bool}, place_after_loop::Vector{Bool}, op::Operation,
    loopsym::Symbol, _n::Int, u₁loop::Symbol, u₂loop::Symbol, vectorized::Symbol, u₂max::Int
)
  lo = ls.loop_order
  id = identifier(op)
  included_vars[id] || return nothing
  loopsym ∈ loopdependencies(op) || return nothing
  for opp ∈ parents(op) # ensure parents are added first
    addoptoorder!(ls, included_vars, place_after_loop, opp, loopsym, _n, u₁loop, u₂loop, vectorized, u₂max)
  end
  included_vars[id] || return nothing
  included_vars[id] = false
  isunrolled = (isu₁unrolled(op)) + 1
  istiled = isu₂unrolled(op) + 1
  # optype = Int(op.node_type) + 1
  after_loop = place_after_loop[id] + 1
  if !isloopvalue(op)
    isnopidentity(ls, op, u₁loop, u₂loop, vectorized, u₂max) || push!(lo[isunrolled,istiled,after_loop,_n], op)
    # if istiled
    #     isnopidentity(ls, op, u₁loop, u₂loop, vectorized, u₂max, u₂max) || push!(lo[isunrolled,2,after_loop,_n], op)
    # else
    #     isnopidentity(ls, op, u₁loop, u₂loop, vectorized, u₂max, nothing) || push!(lo[isunrolled,1,after_loop,_n], op)
    # end
  end
  # @show op, after_loop
  # isloopvalue(op) || push!(lo[isunrolled,istiled,after_loop,_n], op)
  # all(opp -> iszero(length(reduceddependencies(opp))), parents(op)) &&
  set_upstream_family!(place_after_loop, op, false, loopdependencies(op), identifier(op)) # parents that have already been included are not moved, so no need to check included_vars to filter
  nothing
end
function replace_reduct_init!(ls::LoopSet, op::Operation, opsub::Operation, opcheck::Operation)
  deleteat!(parents(op), 2)
  op.variable = opcheck.variable
  opsub.variable = opcheck.variable
  op.mangledvariable = opcheck.mangledvariable
  opsub.mangledvariable = opcheck.mangledvariable
  op.instruction = instruction(:identity)
  fill_children!(ls)
end
function nounrollreduction(op::Operation, u₁loop::Symbol, u₂loop::Symbol, vectorized::Symbol)
  reduceddeps = reduceddependencies(op)
  (vectorized ∉ reduceddeps) &&
    (u₁loop ∉ reduceddeps) &&
    (u₂loop ∉ reduceddeps)
end
function load_short_static_reduction_first!(ls::LoopSet, u₁loop::Symbol, u₂loop::Symbol, vectorized::Symbol)
  for op ∈ operations(ls)
    iscompute(op) || continue
    length(reduceddependencies(op)) == 0 && continue
    parents_op = parents(op)
    length(parents_op) == 2 || continue
    found = false
    parent₁deps = loopdependencies(parents_op[1])
    parent₂deps = loopdependencies(parents_op[2])
    for reduced_dep ∈ reduceddependencies(op)
      if (reduced_dep ∈ parent₁deps) || (reduced_dep ∈ parent₂deps)
        found = true
        break
      end
    end
    found || continue
    if (instruction(op).instr === :reduced_add)
      vecloop = getloop(ls, vectorized)
      if isstaticloop(vecloop) && (length(vecloop) ≤ 16) && nounrollreduction(op, u₁loop, u₂loop, vectorized)
        opsub = parents(op)[2]
        length(children(opsub)) == 1 || continue
        opsearch = parents(op)[1]
        opcheck = search_for_reductinit!(opsearch, opsub, name(opsearch), loopdependencies(op))
        opcheck === opsearch || replace_reduct_init!(ls, op, opsub, opcheck)
      end
    elseif (instruction(op).instr === :add_fast) && (instruction(first(parents(op))).instr === :identity)
      vecloop = getloop(ls, vectorized)
      if isstaticloop(vecloop) && (length(vecloop) ≤ 16) && nounrollreduction(op, u₁loop, u₂loop, vectorized)
        opsub = parents(op)[2]
        ((length(reduceddependencies(opsub)) == 0) & (length(children(opsub)) == 1)) || continue
        opsearch = parents(op)[1]
        opcheck = search_for_reductinit!(opsearch, opsub, name(opsearch), loopdependencies(op))
        opcheck === opsearch || replace_reduct_init!(ls, op, opsub, opcheck)
      end
    end
  end
end

function fillorder!(ls::LoopSet, order::Vector{Symbol}, u₁loop::Symbol, u₂loop::Symbol, u₂max::Int, vectorized::Symbol)
    load_short_static_reduction_first!(ls, u₁loop, u₂loop, vectorized)
    lo = ls.loop_order
    resize!(lo, length(ls.loopsymbols))
    ro = lo.loopnames # reverse order; will have same order as lo
    nloops = length(order)
    ops = operations(ls)
    nops = length(ops)
    included_vars = getroots!(resize!(ls.included_vars, nops), ls)
    #for i ∈ eachindex(included_vars)
    #    included_vars[i] = !included_vars[i]
    #end
    #included_vars = fill!(resize!(ls.included_vars, nops), false)
    place_after_loop = fill!(resize!(ls.place_after_loop, nops), true)
    # to go inside out, we just have to include all those not-yet included depending on the current sym
    empty!(lo)
    for _n ∈ 1:nloops
        n = 1 + nloops - _n
        ro[_n] = loopsym = order[n]
        #loopsym = order[n]
        for op ∈ ops
          addoptoorder!( ls, included_vars, place_after_loop, op, loopsym, _n, u₁loop, u₂loop, vectorized, u₂max )
        end
    end
end


text:
################################################################################
# Abstract Types for Placement
################################################################################

##########
# SANode #
##########

"""
Abstract super types for the `SA` representation of `TaskgraphNodes`.

API
---
* [`location`](@ref)
* [`assign`](@ref)
* [`getclass`](@ref)
* [`setclass!`](@ref)

Implementations
---------------
* [`BasicNode`](@ref)
"""
abstract type SANode end

# API
"""
    location(node::SANode{T}) :: T

Return the location of `node`. Must be parameterized by `T`.
"""
function location end

"""
    assign(node::SANode{T}, location::T)

Set the location of `node` to `location`.
"""
function assign end

"""
    getclass(node)

Return the class of `node`.
"""
function getclass end

"""
    setclass!(node, class::Integer)

Set the class of `node` to `class`.
"""
function setclass! end

#############
# SAChannel #
#############
"""
[`SAStruct`](@ref) representation of a [`TaskgraphEdge`](@ref). Comes in two
varieties: [`TwoChannel`](@ref) and [`MultiChannel`](@ref)
"""
abstract type SAChannel end

"""
Abstract supertype for channels with only one source and sink.

Required Fields
---------------
* `source::Int64`
* `sink::Int64`

Implementations
---------------
* [`BasicChannel`](@ref)
"""
abstract type TwoChannel <: SAChannel end

"""
Abstract supertype for channels with multiple sources/sinks.

Required Fields
---------------
* `sources::Vector{Int}`
* `sinks::Vector{Int}`

Implementations
---------------
* [`BasicMultiChannel`](@ref)
"""
abstract type MultiChannel <: SAChannel end

###############
# AddressData #
###############

"""
Supertype for containers of data for address specific placement. There is no
API for this type since the specific needs of address data vary between
applications. If a custom type is used, extend [`address_cost`] to get the
desired behavior.

Implementations
---------------
* [`EmptyAddressData`](@ref)
* [`DefaultAddressData`](@ref)
"""
abstract type AddressData end

"""
Null representation of [`AddressData`](@ref). Used when there is no address
data to be used during placement.
"""
struct EmptyAddressData <: AddressData end

"""
Default implementation of address data when it is to be used. In its normal
state, it is just a wrapper for a `Dict` mapping addresses to a cost. Look
at the implementation of [`address_cost`](@ref) to see how this is used. This
function may be exteneded on `node` to provide different behavior.

To use this type, the method [`address_data`](@ref) must be defined to
encode the values for the dict.
"""
struct DefaultAddressData{U,T} <: AddressData
    dict::Dict{U,T}
end

getindex(A::DefaultAddressData{U}, i::U) where {U} = getindex(A.dict, i)

################################################################################
# SA Struct
################################################################################

"""
Datastructure for simulated annealing placement.

Important parameters:

* `A` - The concrete Architecture type.

Constructor
-----------
Arguments:
* `m`: The `Map` to translate into an `SAStruct`.

Keyword Arguments:
* `distance`: The distance type to use. Defaults: [`BasicDistance`](@ref)

* `enable_flattness :: Bool`: Enable the flat architecture optimization if
    it is applicable. Default: `true`.

* `enable_address :: Bool`: Enable address-specific data to be incorporated
    into the struct. Default: `false`.

* `aux`: Auxiliary data struct to provide any extra information that may be
    needed for specializations of placement. Default: `nothing`.
"""
struct SAStruct{
    T<:RuleSet,U<:SADistance,D,N<:SANode,L<:SAChannel,M<:AbstractMapTable,A<:AddressData,Q
}
    ruleset::T

    "`Vector{N}`: Container of nodes."
    nodes::Vector{N}
    "`Vector{L}`: Container of edges."
    channels::Vector{L}
    maptable::M
    distance::U
    grid::Array{Int64,D}
    address_data::A
    aux::Q
    # Map back fields
    pathtable::PathTable{D,Path{Component}}
    tasktable::Dict{String,Int64}
end

# Convenience decoding methods
dimension(::SAStruct{T,U,D}) where {T,U,D} = D
MapperCore.rules(sa_struct::SAStruct) = sa_struct.ruleset
nodetype(s::SAStruct) = typeof(s.nodes)
channeltype(s::SAStruct) = typeof(s.channels)
distancetype(::SAStruct{T,U}) where {T,U} = U

isflat(x) = false
isflat(::SAStruct{A,U,D,D}) where {A,U,D} = true
Base.eltype(sa_struct::SAStruct) = location_type(sa_struct.maptable)

################################################################################
# Basis SA Node
################################################################################

"""
The standard implementation of [`SANode`](@ref).
"""
mutable struct BasicNode{T} <: SANode
    "Location this node is assigned in the architecture. Must be parametric."
    location::T
    "The class of this node."
    class::Int
    "Adjacency list of outgoing channels."
    outchannels::Vector{Int64}
    "Adjacency list of incoming channels."
    inchannels::Vector{Int64}
end

# Node Interface
@inline location(n::SANode) = n.location
@inline location(x::CartesianIndex) = x
@inline assign(n::SANode, l) = (n.location = l)
@inline getclass(n::SANode) = n.class
@inline setclass!(n::SANode, class) = n.class = class

isnormal(node::SANode) = isnormal(class(node))
isnormal(class::Int64) = class > 0

# -------------
# Construction
# -------------
function buildnode(::RuleSet, n::TaskgraphNode, x)
    return BasicNode(x, 0, Int64[], Int64[])
end

"""
    setup_node_build(ruleset, taskgraph, location_type)
"""
function setup_node_build(ruleset::RuleSet, t::Taskgraph, ::PathTable{D}) where {D}
    return [buildnode(ruleset, n, zero(Address{D})) for n in getnodes(t)]
end

################################################################################
# Basis SA Edge
################################################################################

"Basic Implementation of [`TwoChannel`](@ref)"
struct BasicChannel <: TwoChannel
    source::Int64
    sink::Int64
end

"Basic Implementation of [`MultiChannel`](@ref)"
struct BasicMultiChannel <: MultiChannel
    sources::Vector{Int64}
    sinks::Vector{Int64}
end

function setup_channel_build(ruleset::RuleSet, taskgraph)
    edges = getedges(taskgraph)
    nodes = getnodes(taskgraph)
    node_dict = Dict(n.name => i for (i, n) in enumerate(nodes))
    # Make source and sink vectors
    sources = map(edges) do edge
        [node_dict[i] for i in edge.sources]
    end
    sinks = map(edges) do edge
        [node_dict[i] for i in edge.sinks]
    end
    # Pass this to the `buildchannels` function
    return buildchannels(ruleset, edges, sources, sinks)
end

buildchannel(ruleset::RuleSet, edge, sources, sinks) = BasicMultiChannel(sources, sinks)

function buildchannels(ruleset::RuleSet, edges, sources, sinks)
    # Get the maximum length of sources and sinks. Use this to determine
    # which type of channels to build.
    channels = [
        buildchannel(ruleset, edge, source, sink) for
        (edge, source, sink) in zip(edges, sources, sinks)
    ]
    max_length = maximum(Iterators.flatten((length.(sources), length.(sinks))))

    # Optimization - if the maximum length is 1 and the element type of `channels` is the
    # default multichannel, change `channels` to be just Basic Channels
    if max_length == 1 && eltype(channels) == BasicMultiChannel
        @info "Optimizing for single source-sink connections"
        channels = [BasicChannel(first(i), first(j)) for (i, j) in zip(sources, sinks)]
    end

    return channels
end

################################################################################
# Address Data
################################################################################

# Document the "address_data" that users must provide to opt-in to address
# specific data generation.

"""
    address_data(ruleset::RuleSet, component::Component) :: T where T

Return some token representing address specific data for `component` under
`ruleset`.
"""
function address_data end

# Default to assigning zero(Float64) to all addresses
address_data(ruleset::RuleSet, component::Component) = zero(Float64)

function build_address_data(
    ruleset::RuleSet, toplevel::TopLevel{D}, pathtable; isflat = false
) where {D}
    data(address) = [address_data(ruleset, toplevel[path]) for path in pathtable[address]]
    address_data = Dict(
        addr => data(addr) for
        addr in CartesianIndices(pathtable) if length(pathtable[addr]) > 0
    )
    # If the flat optimization is turned on - remove the vectors from the values
    # in the dictionary.
    if isflat
        return Dict(k => first(v) for (k, v) in address_data)
    else
        return address_data
    end
end

################################################################################
# Constructor for the SA Structure
################################################################################
function SAStruct(
    m::Map{D};
    # Enable address-specific data.
    enable_address = false,
    aux = nothing,
    kwargs...,
) where {D}

    #@debug "Building SA Placement Structure\n"

    # Unpack some data structures for easier reference.
    toplevel = m.toplevel
    taskgraph = m.taskgraph
    ruleset = rules(m)

    pathtable = build_pathtable(toplevel, ruleset)

    # TODO: Make extendable.
    distance = BasicDistance(toplevel, pathtable)

    # Build SA Node types and make a record mapping node name to the index of
    # that nodes representation in this data structure.
    nodes = setup_node_build(rules(m), taskgraph, pathtable)
    tasktable = Dict(n.name => i for (i, n) in enumerate(getnodes(taskgraph)))

    # Build Channel
    channels = setup_channel_build(rules(m), taskgraph)
    @debug """
        Node Type: $(typeof(nodes))
        Edge Type: $(typeof(channels))
        """

    # Assign adjacency information to nodes.
    record_channels!(nodes, channels)

    #----------------------------------------------------#
    # Obtain task equivalence classes from the taskgraph #
    #----------------------------------------------------#
    classes = equivalence_classes(rules(m), getnodes(taskgraph))

    # Assign each node with its given classification.
    for (index, class) in enumerate(classes.classes)
        setclass!(nodes[index], class)
    end

    maptable = MapTable(toplevel, rules(m), classes, pathtable)
    grid = zeros(size(pathtable)...)

    if enable_address
        address_data = build_address_data(rules(m), toplevel, pathtable; isflat = isflat)
    else
        address_data = EmptyAddressData()
    end

    sa = SAStruct{
        typeof(rules(m)),     # RuleSet Type
        typeof(distance),     # Encoding of Tile Distance
        ndims(grid),          # Architecture Dimensionality + 1
        eltype(nodes),        # Type of the taskgraph nodes
        eltype(channels),     # Type of the taskgraph channels
        typeof(maptable),
        typeof(address_data),  # Type of address data
        typeof(aux),
    }(
        rules(m),
        nodes,
        channels,
        maptable,
        distance,
        grid,
        address_data,
        aux,
        pathtable,
        tasktable,
    )

    # Run initial placement and verify result.
    initial_placement!(sa)
    verify_placement(m, sa)
    return sa
end

cleargrid(sa::SAStruct) = clear(sa.grid)
clear(x::Array{T}) where {T} = x .= zero(T)

function preplace(m::Map, sa::SAStruct)
    cleargrid(sa)

    offset = getoffset(m.toplevel)
    for (taskname, path) in m.mapping.nodes
        address = getaddress(m.toplevel, path) + offset
        location = findaddress(sa.pathtable, address, path)
        # Get the index to assign
        index = sa.tasktable[taskname]
        # Assign the nodes
        assign(sa, index, location)
    end
    return nothing
end

################################################################################
# Write back method
################################################################################
function record(m::Map, sa::SAStruct)
    verify_placement(m, sa)
    mapping = m.mapping
    tasktable_rev = rev_dict(sa.tasktable)

    for (index, node) in enumerate(sa.nodes)
        # Get the mapping for the node
        address = location(node)
        # Get the component name in the original toplevel
        path = sa.pathtable[address]
        task_node_name = tasktable_rev[index]
        # Create an entry in the "Mapping" data structure
        mapping.nodes[task_node_name] = path
    end
    return nothing
end

function record(nodes, i, channel::TwoChannel)
    push!(nodes[channel.source].outchannels, i)
    push!(nodes[channel.sink].inchannels, i)
    return nothing
end

function record(nodes, i, channel::MultiChannel)
    for j in channel.sources
        push!(nodes[j].outchannels, i)
    end
    for j in channel.sinks
        push!(nodes[j].inchannels, i)
    end
    return nothing
end

function record_channels!(nodes, channels)
    # Reverse populate the nodes so they track their channels.
    for (i, channel) in enumerate(channels)
        record(nodes, i, channel)
    end
    return nothing
end

function equivalence_classes(ruleset::RuleSet, nodes)
    classes = zeros(Int, length(nodes))

    # Representitives for each class
    reps = eltype(nodes)[]

    for (index, node) in enumerate(nodes)
        # See if this node class has already been discovered.
        i = findfirst(x -> isequivalent(ruleset, x, node), reps)
        # If not, add this node to the vector of representatives.
        if i == nothing
            push!(reps, node)
            i = length(reps)
        end
        # Record class
        classes[index] = i
    end
    return (classes = classes, reps = reps)
end

################################################################################
# Verification routine for SA Placement
################################################################################
function verify_placement(m::Map, sa::SAStruct)
    # Assert that the SAStruct belongs to the same toplevel
    @debug "Verifying Placement"

    bad_nodes = check_grid_population(sa)
    append!(bad_nodes, check_consistency(sa))
    append!(bad_nodes, check_mapability(m, sa))
    # Gather all the unique bad nodes and sort the final list.
    bad_nodes = sort(unique(bad_nodes))
    # Routine passes check if length of bad_nodes is 0
    passed = length(bad_nodes) == 0
    if passed
        @debug "Placement Verified"
    else
        @error begin
            bad_dict = Dict{Int64,String}()
            for (i, name) in enumerate(keys(m.taskgraph.nodes))
                if i in bad_nodes
                    bad_dict[i] = name
                end
            end
            """
            Placement Failed.
            Offending Node Names:
            $bad_dict
            """
        end
    end
    return passed
end

function check_grid_population(sa::SAStruct)
    # Iterate through all entries in the grid. Record the indices encountered
    # along the way. When an index is discovered, mark it as discovered.
    #
    # If an index is found twice - this is a problem. Print an error and mark
    # the test as failed.
    #
    # After this routine - make sure that all nodes are accounted for.

    # Use this to return the indices of tasks that are troublesome.
    bad_nodes = Int64[]
    found = fill(false, length(sa.nodes))
    for g in sa.grid
        g == 0 && continue
        if found[g]
            @warn "Found node $g more than once"
            push!(bad_nodes, g)
        else
            found[g] = true
        end
    end
    # Make sure all nodes have been found
    for i in 1:length(found)
        if found[i] == false
            @warn "Node $i not placed."
            push!(bad_nodes, i)
        end
    end
    return bad_nodes
end

function check_consistency(sa::SAStruct)
    bad_nodes = Int64[]
    # Verify that addresses for the nodes match the grid
    for (index, node) in enumerate(sa.nodes)
        node_assigned = sa.grid[location(node)]
        if index != node_assigned
            push!(bad_nodes, index)
            push!(bad_nodes, node_assigned)
            @warn """
                Data structure inconsistency for node $index.
                Node assigned to location: $(location(node)).

                Node assigned in the grid at this location: $node_assigned.
                """
        end
    end
    return bad_nodes
end

function check_mapability(m::Map, sa::SAStruct)
    bad_nodes = Int64[]
    toplevel = m.toplevel
    # Iterate through each node in the SA
    for (index, (m_node_name, m_node)) in enumerate(m.taskgraph.nodes)
        sa_node = sa.nodes[index]
        # Get the mapping for the node
        address = location(sa_node)
        #pathindex = getindex(sa_node)
        # Get the component name in the original toplevel
        path = sa.pathtable[address]
        # Get the component from the toplevel
        component = toplevel[path]
        if !canmap(rules(m), m_node, component)
            push!(bad_nodes, index)
            @warn """
                Node index $m_node_name incorrectly assigned to toplevel
                node $(path).
                """
        end
    end
    return bad_nodes
end


text:
module Exceptions

import Revise
import Genie
import HTTP

export ExceptionalResponse, RuntimeException, InternalServerException, NotFoundException, FileExistsException

struct ExceptionalResponse <: Exception
  response::HTTP.Response
end

Base.show(io::IO, ex::ExceptionalResponse) = print(io, "ExceptionalResponseException: $(ex.response.status) - $(Dict(ex.response.headers))")

###

struct RuntimeException <: Exception
  message::String
  info::String
  code::Int
  ex::Union{Nothing,Exception}
end

RuntimeException(message::String, code::Int) = RuntimeException(message, "", code, nothing)
RuntimeException(message::String, info::String, code::Int) = RuntimeException(message, info, code, nothing)

Base.show(io::IO, ex::RuntimeException) = print(io, "RuntimeException: $(ex.code) - $(ex.info) - $(ex.message)")

###

struct InternalServerException <: Exception
  message::String
  info::String
  code::Int
end

InternalServerException(message::String) = InternalServerException(message, "", 500)
InternalServerException() = InternalServerException("Internal Server Error")

###

struct NotFoundException <: Exception
  message::String
  info::String
  code::Int
  resource::String
end

NotFoundException(resource::String) = NotFoundException("$resource can not be found", "", 404, resource)
NotFoundException() = NotFoundException("")

###

struct FileExistsException <: Exception
  path::String
end

Base.show(io::IO, ex::FileExistsException) = print(io, "FileExistsException: $(ex.path)")

end


text:
# This file contains a function to reset
# an instance to default values

"""
    reset!(input::InputStruct)

Reset an instance of [`InputStruct`](@ref) to default values.

# Usage

```jldoctest; output = false
using Scats
s = Scats.API()
s.Input.reset!()

# output

```
"""
function reset!(input::InputStruct)

    # Reset all elements
    input.N = 0
    input.Δt = input.q = 0.0
    input.t = input.x = []

    return nothing

end


text:
using Gadfly

x1 = rand(40)
y1 = 4.*x1 .+ 2 .+randn(40)
x2 = rand(40)
y2 = -6.*x2 .+ 3 .+ randn(40)
x  = [x1;x2]
y  = [y1;y2]
col = [fill("Slope 4",40); fill("Slope -6",40)]
plot(x=x,y=y,colour=col, Geom.point, Geom.smooth(method=:lm))


text:
#=================================================================================
  NTR is a trust region sub-problem solver that solves the quadratic programming
  problem subject to the solution being in a Euclidean Ball. It is appropriate
  for all types of Hessians. It solves the problem by solving the secular equ-
  ation by a safeguarded Newton's method. The secular equation is simply the
  inverse of the difference between the step length and the trust region radius.
  It accepts either the Newton step if this is interior and the Hessian is po-
  sitive definite, or steps to the boundary (with some slack) in a clever way.

  NTR accepts all Hessians and is therefore well-suited for Newton's method for
  any vexity. It is expensive for very large problems, but uses Cholesky fac-
  torizations in clever ways, so it should be less expensive than the version
  that calculates the full eigensolution.
=================================================================================#

# TODO allow passing in a lambda and previous cholesky if the
# solution was not accepted
# As described in Algorith, 7.3.4 in [CGTBOOK]
struct NTR <: TRSPSolver end
summary(::NTR) = "Trust Region (Newton, cholesky)"

function initial_λs(∇f, H, Δ)
    n = length(∇f)
    T = eltype(∇f)
    Hfrob = norm(H, 2)
    Hinf = norm(H, Inf)
    max_h = maximum(diag(H))

    norm_to_Δ = norm(∇f) / Δ

    sumabs_j = sum(abs, view(H, 2:n, 1))
    max_i_p = sumabs_j + H[1, 1]
    max_i_m = sumabs_j - H[1, 1]
    for j = 2:n
        sumabs_j = T(0)
        for i = 1:n-1
            i == j && continue
            sumabs_j += abs(H[i, j])
        end
        max_i_p = max(max_i_p, sumabs_j + H[j, j])
        max_i_m = max(max_i_m, sumabs_j - H[j, j])
    end
    λL = max(T(0), norm_to_Δ - min(max_i_p, Hfrob, Hinf), -max_h)
    λU = max(T(0), norm_to_Δ + min(max_i_m, Hfrob, Hinf))
    return λL, λU
end
λ⁺_newton(λ, w, Δ) = λ + (s₂^2 / dot(w, w)) * (s₂ - Δ) / Δ
function (ms::NTR)(
    ∇f,
    H,
    Δ::T,
    s,
    scheme,
    λ0 = 0;
    abstol = 1e-10,
    maxiter = 50,
    κeasy = T(1) / 10,
    κhard = T(2) / 10,
) where {T}
    # λ0 might not be 0 if we come from a failed TRS solve
    λ = T(λ0)
    θ = T(1) / 2
    n = length(∇f)
    h = H isa UniformScaling ? copy(∇f) .* 0 .+ 1 : diag(H)
    H = H isa UniformScaling ? Diagonal(copy(∇f) .* 0 .+ 1) : H

    # Check for interior convergence
    if λ == T(0)
        F = cholesky(Symmetric(H); check = false)
        s .= -∇f
        s .= F \ s
        s₂ = norm(s, 2)
        if issuccess(F) && s₂ < Δ
            H = update_H!(H, h)
            return tr_return(;
                λ = λ,
                ∇f = ∇f,
                H = H,
                s = s,
                interior = true,
                solved = true,
                hard_case = false,
                Δ = Δ,
            )
        end
    end

    # If the solution was not internal, start Newton's method on the
    # secular equation.
    isg = initial_safeguards(H, h, ∇f, Δ)
    λ = safeguard_λ(λ, isg)
    λL, λU = isg.L, isg.U

    for iter = 1:maxiter
        H = update_H!(H, h, λ)
        F = cholesky(Symmetric(H); check = false)
        in𝓖, linpack = false, false
        #===========================================================================
         If F is successful, then H is positive definite, and we can safely look
         at the Newton step. If this is interior, we're done, if not, we're either
         in L or G. In L, the Newton step on the secular equation  stays in L and
         is safe to take. G requires more work.
        ===========================================================================#
        if issuccess(F)
            # H(λ) is PD, so we're in 𝓕

            # Algorithm 7.3.1 on p. 185 in [ConnGouldTointBook]
            # Step 1 was factorizing
            # Step 2
            s .= -∇f
            s .= F \ s

            # Check if step is approximately equal to the radius
            s₂ = norm(s, 2)
            if s₂ ≈ Δ
                H = update_H!(H, h)
                return tr_return(;
                    λ = λ,
                    ∇f = ∇f,
                    H = H,
                    s = s,
                    interior = false,
                    solved = true,
                    hard_case = false,
                    Δ = Δ,
                )
            end
            if s₂ < Δ # in 𝓖 because we're in 𝓕, but curve below Δ
                # we're in 𝓖 so λ is a new upper bound; λᴹ < λ
                in𝓖 = true
                λU = λ
            else # λ ∈ 𝓛
                # in 𝓛 λ is a *lower* bound instead
                λL = λ
            end

            # Step 3
            w = F.U' \ s

            # Step 4
            # Newton trial step
            λ⁺ = λ + ((s₂ - Δ) / Δ) * (s₂^2 / dot(w, w))
            if in𝓖
                linpack = true
                w, u = λL_with_linpack(F)
                λL = max(λL, λ - dot(u, H * u))

                α, s_g, m_g = 𝓖_root(u, s, Δ, ∇f, H)
                s .= s_g

                s₂ = norm(s)
                # check hard case convergnce
                if α^2 * dot(u, H * u) ≤ κhard * (dot(s, H * s) + λ * Δ^2)
                    H = update_H!(H, h)
                    return tr_return(;
                        λ = λ,
                        ∇f = ∇f,
                        H = H,
                        s = s,
                        interior = false,
                        solved = true,
                        hard_case = true,
                        Δ = Δ,
                        m = m_g,
                    )
                end
                # If not the hard case solution, try to factorize H(λ⁺)
                H = update_H!(H, h, λ⁺)
                F = cholesky(H; check = false)
                if issuccess(F) # Then we're in L, great! lemma 7.3.2
                    λ = λ⁺
                else # we landed in N, this is bad, so use bounds to approach L
                    λ = max(sqrt(λL * λU), λL + θ * (λU - λL))
                end
            else # in L, we can safely step
                λ = λ⁺
            end

            # check for convergence
            if in𝓖 && abs(s₂ - Δ) ≤ κeasy * Δ
                H = update_H!(H, h)
                return tr_return(;
                    λ = λ,
                    ∇f = ∇f,
                    H = H,
                    s = s,
                    interior = false,
                    solved = true,
                    hard_case = false,
                    Δ = Δ,
                )
            elseif abs(s₂ - Δ) ≤ κeasy * Δ # implicitly "if in 𝓕" since we're in that branch
                # u and α comes from linpack
                if linpack
                    if α^2 * dot(u, H * u) ≤ κhard * (dot(sλ, H * sλ) * Δ^2)
                        s .= s .+ α * u
                        H = update_H!(H, h)
                        return tr_return(;
                            λ = λ,
                            ∇f = ∇f,
                            H = H,
                            s = s,
                            interior = false,
                            solved = true,
                            hard_case = false,
                            Δ = Δ,
                        )
                    end
                end
            end
        else # λ ∈ 𝓝, because the factorization failed (typo in CGT)
            # Use partial factorization to find δ and v such that
            # H(λ) + δ*e*e' = 0. All we can do here is to find a better
            # lower bound, we cannot apply the Newton step here.
            δ, v = λL_in_𝓝(H, F)
            λL = max(λL, λ + δ / dot(v, v)) # update lower bound
            λ = max(sqrt(λL * λU), λL + θ * (λU - λL)) # no converence possible, so step in bracket
        end
    end
    H = update_H!(H, h)
    tr_return(;
        λ = λ,
        ∇f = ∇f,
        H = H,
        s = s,
        interior = true,
        solved = false,
        hard_case = false,
        Δ = Δ,
    )
end

function λL_in_𝓝(H, F)
    T = eltype(F)
    n = first(size(F))
    δ = sum(abs2, view(F.factors, 1:(F.info-1), F.info)) - H[F.info, F.info]
    v = zeros(T, n)
    v[F.info] = 1
    for j = (F.info-1):-1:1
        vj = zero(T)
        for i = (j+1):F.info
            vj += F.factors[j, i] * v[i]
        end
        v[j] = -vj / F.factors[j, j]
    end
    return δ, v
end

function λL_with_linpack(F)
    T = eltype(F)
    n = first(size(F))
    w = zeros(T, n)
    num_p1 = inv(F.factors[end, end])
    num_m1 = -inv(F.factors[end, end])
    w[end] = max(num_p1, num_m1)
    for k = n-1:-1:1
        ukk = F.factors[k, k]
        num = sum(abs2, view(F.factors, 1:(k-1), k))
        w[k] = max((1 - num) / ukk, (-1 - num) / ukk)
    end
    sol = F.factors \ w
    w, sol ./ norm(sol)
end

function 𝓖_root(u, s, Δ, ∇f, H)
    pa = sum(abs2, u)
    pb = 2 * dot(u, s)
    pd = sqrt(4 * pb^2 - pa * (sum(abs2, s) - Δ^2))
    α₁ = (-pb + pd) / 2pa
    α₂ = (-pb - pd) / 2pa

    s₁ = s + α₁ * u
    m₁ = dot(∇f, s₁) + dot(s₁, H * s₁) / 2
    s₂ = s + α₂ * u
    m₂ = dot(∇f, s₂) + dot(s₂, H * s₂) / 2
    α, s, m = m₁ ≤ m₂ ? (α₁, s₁, m₁) : (α₂, s₂, m₂)
    α, s, m
end

max_stars_repo_name: JuliaNLSolvers/NLSolvers

max_issues_repo_name: JuliaNLSolvers/NLSolvers

max_forks_repo_name: JuliaNLSolvers/NLSolvers


text:
################################################################################
## Problem Definition -- frictionless pendulum

struct Pendulum{P,S,T,K}

    prob::P
    solver::S
    sensealg::T
    kwargs::K

    function Pendulum(; solver = Tsit5(), sensalg = ForwardDiffSensitivity(), kwargs...)
        # Parameters and initial conditions only
        # used to initialize the ODE problem
        u₀ = Float32[1.0, 1.0]
        p = Float32[1.]
        tspan = (0.f0, 1.f0)

        # Define differential equations
        function f!(du, u, p, t)
                x, y = u
                G = 10.0f0
                L = p[1]

                du[1] = y
                du[2] =  -G/L*sin(x)
        end

        # Build ODE Problem
        _prob = ODEProblem(f!, u₀, tspan, p)

        @info "Optimizing ODE Problem"
        # Make ODE Problem more performant using ModelingToolkit
        # not necessary for this simple pendulum but can very usefull
        # for other more complex systems
        sys = modelingtoolkitize(_prob)
        ODEFunc = ODEFunction(sys, tgrad=true, jac = true, sparse = false, simplify = false)
        prob = ODEProblem(ODEFunc, u₀, tspan, p)

        P = typeof(prob)
        S = typeof(solver)
        T = typeof(sensalg)
        K = typeof(kwargs)
        new{P,S,T,K}(prob, solver, sensalg, kwargs)
    end

end

################################################################################
## Problem Definition -- Pendulum with friction

struct Pendulum_friction{P,S,T,K}

    prob::P
    solver::S
    sensealg::T
    kwargs::K

    function Pendulum_friction(; solver = Tsit5(), sensalg = ForwardDiffSensitivity(), kwargs...)
        # Default parameters and initial conditions
        u₀ = Float32[1.0, 1.0]
        p = Float32[1.]
        tspan = (0.f0, 1.f0)

        # Define differential equations
        function f!(du, u, p, t)
                x, y = u
                G = 10.0f0
                m = 1.0f0
                b = 0.7f0
                L = p[1]

                du[1] = y
                du[2] =  -G/L*sin(x) - (b/m) * y
        end

        # Build ODE Problem
        _prob = ODEProblem(f!, u₀, tspan, p)

        @info "Optimizing ODE Problem"
        sys = modelingtoolkitize(_prob)
        ODEFunc = ODEFunction(sys, tgrad=true, jac = true, sparse = false, simplify = false)
        prob = ODEProblem(ODEFunc, u₀, tspan, p)

        P = typeof(prob)
        S = typeof(solver)
        T = typeof(sensalg)
        K = typeof(kwargs)
        new{P,S,T,K}(prob, solver, sensalg, kwargs)
    end

end

###############################################################################
# Problem Definition -- Stochastic Pendulum

struct SPendulum{P,S,T,K}

    prob::P
    solver::S
    sensealg::T
    kwargs::K

    function SPendulum(;solver = SOSRI(), sensalg = ForwardDiffSensitivity(), kwargs...)
        # Default parameters and initial conditions
        u₀ = Float32[1.0, 1.0]
        p = Float32[1.]
        tspan = (0.f0, 1.f0)

        # Define differential equations
        function f!(du, u, p, t)
                x, y = u
                G = 10.0f0
                L = p[1]

                du[1] = y
                du[2] =  -G/L*sin(x)
        end

        # Build ODE Problem
        _prob = ODEProblem(f!, u₀, tspan, p)

        @info "Optimizing ODE Problem"
        sys = modelingtoolkitize(_prob)
        ODEFunc = ODEFunction(sys, tgrad=true, jac = true, sparse = false, simplify = false)
        prob = ODEProblem(ODEFunc, u₀, tspan, p)

        function σ(du,u,p,t)
            du .= 0.01f0
        end

        prob_sde = SDEProblem(prob.f.f, σ, prob.u0, prob.tspan, prob.p, jac = prob.f.jac, tgrad = prob.f.tgrad)

        P = typeof(prob_sde)
        S = typeof(solver)
        T = typeof(sensalg)
        K = typeof(kwargs)
        new{P,S,T,K}(prob_sde, solver, sensalg, kwargs)
    end

end

max_issues_repo_name: gabrevaya/LatentDiffEq

max_forks_repo_name: gabrevaya/LatentDiffEq


text:
"""
Makes sure all sage benchmarks import properly.
For now, just doing first 5 out of 25, since polynomial examples are not in R+.
"""
function test_sagemark_to_GlobalModel()
    idxs = 1:5; # 25
    # max_min = [26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38]
    for idx in idxs
        ex = sagemark_to_GlobalModel(idx);
    end
    idx = 1;
    gm = sagemark_to_GlobalModel(idx, false);
    gm_lse = sagemark_to_GlobalModel(idx, true);
    classify_constraints.([gm, gm_lse])

    # Testing proper importing using sagemark #1
    gm = sagemark_to_GlobalModel(1, false);
    gm_lse = sagemark_to_GlobalModel(1, true);

    bounds = Dict(all_variables(gm) .=> [[0.1, 1], [5., 10.], [8., 15.], [0.01, 1.], [-Inf, Inf]])
    @test all(get_bounds(gm)[var] ≈ bounds[var] for var in all_variables(gm))
    inp = Dict(all_variables(gm) .=> [1,1.9,3,3.9, 10.]);
    log_inp = Dict(vk => log(val) for (vk, val) in inp);

    @test gm_lse.bbls[1](log_inp)[1] ≈ gm.bbls[1](inp)[1] ≈ [inp[gm.vars[3]] ^ 0.8 * inp[gm.vars[4]] ^ 1.2][1]

    # Checking OCTHaGOnException for sampling unbounded model
    uniform_sample_and_eval!(gm)
    @test true
end

""" Testing that problems are correctly imported, with some random checking. """
function test_gams_to_GlobalModel()
    filenums = [2.15, 2.16, 2.17, 2.18] # [3.2, "3.10", 3.13, 3.15, 3.16, 3.18, 3.25]
    filenames = ["problem" * string(filenum) * ".gms" for filenum in filenums]
    gms = Dict()
    for filename in filenames
        try
            gm = GAMS_to_GlobalModel(OCTHaGOn.GAMS_DIR, filename)
#             println("    Problem NL constraints: " * string(length(gm.bbls)))
            types = JuMP.list_of_constraint_types(gm.model)
            if !isempty(types)
                total_constraints = sum(length(all_constraints(gm.model, type[1], type[2])) for type in types)
#                 println("    Problem total constraints: " * string(total_constraints))
            end
        catch
            throw(OCTHaGOnException(filename * " has an import issue."))
        end
    end

    # Testing problem import.
    filename = "problem3.13.gms"
    gm =  GAMS_to_GlobalModel(OCTHaGOn.GAMS_DIR, filename)
    x = gm.model[:x]
    @test length(gm.vars) == 8
    @test all(bound == [0,100] for bound in values(get_bounds(Array{JuMP.VariableRef}(flat(gm.model[:x])))))
    @test length(gm.bbls) == 1
end

# function test_load_fits()
#     gm = sagemark_to_GlobalModel(3, false)
#     set_optimizer(gm, SOLVER_SILENT)
#     # Testing finding bounds of bounded model
#     @test isnothing(get_unbounds(gm.bbls))
#     @test isnothing(find_bounds!(gm))
#     # Test reloading
#     load_fit(gm)
#     @test all([get_param(bbl, :reloaded) == true for bbl in gm.bbls])
#     add_tree_constraints!(gm)
#     optimize!(gm)
#     @test gm.cost[end] <= -120
# end

test_sagemark_to_GlobalModel()

test_gams_to_GlobalModel()

# test_load_fits()


text:
using Documenter
using NighttimeLights

DocMeta.setdocmeta!(NighttimeLights, :DocTestSetup, :(using NighttimeLights); recursive=true)
makedocs(;
    modules=[NighttimeLights],
    authors="Ayush Patnaik, Ajay Shah, Anshul Tayal, Susan Thomas",
    repo="https://github.com/xKDR/NighttimeLights.jl/blob/{commit}{path}#{line}",
    sitename="NighttimeLights.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://github.com/xKDR/NighttimeLights.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Basic concepts" => "concepts.md",
        "Polygons and shapefiles"=> "polygons.md",
        "Data Cleaning"=> "data_cleaning.md",
        "Tutorial" =>"tutorial.md",
    ],
)

deploydocs(;
    repo="github.com/xKDR/NighttimeLights.jl",
    target = "build",
    devbranch="main"
)


text:
using Flux, Flux.Data.MNIST
using Flux: onehotbatch, argmax, crossentropy, throttle
using Base.Iterators: repeated

# Classify MNIST digits with a simple multi-layer-perceptron

imgs = MNIST.images()
# Stack images into one large batch
X = hcat(float.(reshape.(imgs, :))...)

labels = MNIST.labels()
# One-hot-encode the labels
Y = onehotbatch(labels, 0:9)

m = Chain(
  Dense(28^2, 32, relu),
  Dense(32, 10),
  softmax)

# using CuArrays
# x, y = cu(x), cu(y)
# m = mapleaves(cu, m)

loss(x, y) = crossentropy(m(x), y)

accuracy(x, y) = mean(argmax(m(x)) .== argmax(y))

dataset = repeated((X, Y), 200)
evalcb = () -> @show(loss(X, Y))
opt = ADAM(params(m))

Flux.train!(loss, dataset, opt, cb = throttle(evalcb, 10))

accuracy(X, Y)

# Test set accuracy
tX = hcat(float.(reshape.(MNIST.images(:test), :))...)
tY = onehotbatch(MNIST.labels(:test), 0:9)
accuracy(tX, tY)

max_stars_repo_name: domluna/model-zoo

max_issues_repo_name: domluna/model-zoo

max_forks_repo_name: domluna/model-zoo


text:
# In addition to illustrating transformations, we'll start loading some
# shaders from a file.  We'll also demo using Reactive for animations.

using ModernGL, GeometryTypes, GLAbstraction, GLFW, Images, Reactive
const GLA = GLAbstraction

# Load our texture. See "downloads.jl" to get the images.
kitten = load(GLAbstraction.dir("tutorials", "images", "kitten.png"))
puppy  = load(GLAbstraction.dir("tutorials", "images", "puppy.png"))

# Create the window. This sets all the hints and makes the context current.
window = GLFW.Window(name="Transformations 2", resolution=(800,600))
GLFW.MakeContextCurrent(window)
GLA.set_context!(window)

# The positions of the vertices in our rectangle
vertex_positions = Point{2,Float32}[(-0.5,  0.5),     # top-left
                                    ( 0.5,  0.5),     # top-right
                                    ( 0.5, -0.5),     # bottom-right
                                    (-0.5, -0.5)]     # bottom-left

# The colors assigned to each vertex
vertex_colors = Vec3f0[(1, 0, 0),                     # top-left
                       (0, 1, 0),                     # top-right
                       (0, 0, 1),                     # bottom-right
                       (1, 1, 1)]                     # bottom-left

# The texture coordinates of each vertex
vertex_texcoords = Vec2f0[(0, 0),
                          (1, 0),
                          (1, 1),
                          (0, 1)]

# Specify how vertices are arranged into faces
# Face{N,T,O} type specifies a face with N vertices, with index type
# T (you should choose UInt32), and index-offset O. If you're
# specifying faces in terms of julia's 1-based indexing, you should set
# O=0. (If you instead number the vertices starting with 0, set
# O=-1.)
elements = Face{3,UInt32}[(0,1,2),          # the first triangle
                          (2,3,0)]          # the second triangle

# The vertex shader---note the `vert` in front of """
vertex_shader = GLA.vert"""
#version 150

in vec2 position;
in vec3 color;
in vec2 texcoord;

out vec3 Color;
out vec2 Texcoord;

uniform mat4 trans;

void main()
{
    Color = color;
    Texcoord = texcoord;
    gl_Position = trans*vec4(position, 0.0, 1.0);
}
"""
fragment_shader = load(joinpath(dirname(@__FILE__), "shaders", "puppykitten.frag"))

prog = GLA.Program(GLA.Shader(vertex_shader), fragment_shader)

buffers = GLA.generate_buffers(prog, GLA.GEOMETRY_DIVISOR, position = vertex_positions,
                                                           color = vertex_colors,
                                                           texcoord=vertex_texcoords)

vao  = GLA.VertexArray(buffers, elements)
tex_kitten = GLA.Texture(collect(kitten'))
tex_puppy  = GLA.Texture(collect(puppy'))

# Define the rotation matrix (could also use rotationmatrix_z)
function rotmat_z(angle::T) where T
    T0, T1 = zero(T), one(T)
    Mat{4}(
        cos(angle), sin(angle), T0, T0,
        -sin(angle), cos(angle),  T0, T0,
        T0, T0, T1, T0,
        T0, T0, T0, T1
    )
end

# By wrapping it in a Signal, we can easily update it.
trans = Signal(rotmat_z(0f0))

# Draw until we receive a close event
while !GLFW.WindowShouldClose(window)

    push!(trans, rotmat_z(time()*deg2rad(180)))
    Reactive.run_till_now()

    glClear(GL_COLOR_BUFFER_BIT)
    GLA.bind(prog)
    GLA.gluniform(prog, :texKitten, 0, tex_kitten) #first texture sampler
    GLA.gluniform(prog, :texPuppy, 1, tex_puppy) # second texture sampler
    GLA.gluniform(prog, :trans, value(trans))
    GLA.bind(vao)
    GLA.draw(vao)
    GLFW.SwapBuffers(window)
    GLFW.PollEvents()
    if GLFW.GetKey(window, GLFW.KEY_ESCAPE) == GLFW.PRESS
        GLFW.SetWindowShouldClose(window, true)
    end
end

GLFW.DestroyWindow(window)  # needed if you're running this from the REPL


text:
# This file contains functions for the Revenue General Direct Approach DEA model
"""
    RevenueGDADEAModel
An data structure representing a Revenue General Direct Approach DEA model.
"""
struct RevenueGDADEAModel <: AbstractProfitDEAModel
    n::Int64
    m::Int64
    s::Int64
    measure::Symbol
    rts::Symbol
    monetary::Bool
    dmunames::Union{Vector{String},Nothing}
    eff::Vector
    lambda::SparseMatrixCSC{Float64, Int64}
    techeff::Vector
    alloceff::Vector
    normalization::Vector
    Xtarget::Matrix
    Ytarget::Matrix
end

"""
    dearevenuegda(X, Y, P, measure)
Compute revenue efficiency using data envelopment analysis General Direct Approach model for
inputs `X`, outputs `Y`, price of outputs `P`, and efficiency `measure`.

# Measure specification:
- `:ERG`: Enhanced Russell Graph (or Slack Based Measure (SBM)).

# Optional Arguments
- `rts=:VRS`: choose between constant returns to scale `:CRS` or variable returns to scale `:VRS`.
- `atol=1e-6`: tolerance for DMU to be considered efficient.
- `names`: a vector of strings with the names of the decision making units.

# Examples
```jldoctest
julia> X = [1; 1; 1; 1; 1; 1; 1; 1];

julia> Y = [7 7; 4 8; 8 4; 3 5; 3 3; 8 2; 6 4; 1.5 5];

julia> P = [1 1; 1 1; 1 1; 1 1; 1 1; 1 1; 1 1; 1 1];

julia> dearevenuegda(X, Y, P, :ERG)
General Direct Approach Revenue DEA Model
DMUs = 8; Inputs = 1; Outputs = 2
Returns to Scale = VRS
Returns to Scale = VRS
───────────────────────────────────
    Revenue  Technical   Allocative
───────────────────────────────────
1  0.0        0.0       0.0
2  0.25       0.0       0.25
3  0.25       0.0       0.25
4  0.866667   0.866667  0.0
5  1.33333    1.33333   2.96059e-16
6  1.0        0.5       0.5
7  0.458333   0.458333  0.0
8  2.5        2.05556   0.444444
───────────────────────────────────
```
"""
function dearevenuegda(X::Union{Matrix,Vector}, Y::Union{Matrix,Vector},
    P::Union{Matrix,Vector},
    measure::Symbol;
    rts::Symbol = :VRS, monetary::Bool = false, atol::Float64 = 1e-6,
    names::Union{Vector{String},Nothing} = nothing,
    optimizer::Union{DEAOptimizer,Nothing} = nothing)::RevenueGDADEAModel

    # Check parameters
    nx, m = size(X, 1), size(X, 2)
    ny, s = size(Y, 1), size(Y, 2)

    np, sp = size(P, 1), size(P, 2)

    if nx != ny
        throw(DimensionMismatch("number of rows in X and Y ($nx, $ny) are not equal"));
    end
    if np != ny
        throw(DimensionMismatch("number of rows in P and Y ($np, $ny) are not equal"));
    end
    if sp != s
        throw(DimensionMismatch("number of columns in P and Y ($sp, $s) are not equal"));
    end

    n = nx

    # Get efficiency measure model
    if measure == :ERG
        measuremodel = dearevenuerussell(X, Y, P, rts = rts, optimizer = optimizer)
        measuretechnicalmodel = dearussell(X, Y, orient = :Output, rts = rts, slack = false, optimizer = optimizer)
    else
        throw(ArgumentError("Invalid efficiency `measure`"));
    end

    normmeasure = normfactor(measuremodel)
    rimeasure = efficiency(measuremodel)  .* normmeasure
    timeasure = efficiency(measuremodel, :Technical)
    rlambdaeff = peersmatrix(measuremodel)
    Xtarget = targets(measuremodel, :X)
    Ytarget = targets(measuremodel, :Y)
    Xtargetti = targets(measuretechnicalmodel, :X)
    Ytargetti = targets(measuretechnicalmodel, :Y)

    # Calculate slack as the difference between target and observed
    slackYmeasure = Ytargetti .- Y

    # Get profit inefficiency of projection (target)
    if measure == :ERG
        targetmodel = dearevenuerussell(Xtargetti, Ytargetti, P, rts = rts, optimizer = optimizer)
    end

    ritarget = efficiency(targetmodel) .* normfactor(targetmodel)

    monetaryti = sum(P .* slackYmeasure, dims = 2)

    # Normalization factor (equation 13.13)
    normalization = zeros(n)
    for i=1:n
        if abs(timeasure[i] ) <= atol
            # For technical efficient DMU's
            normalization[i] = normmeasure[i]
        else
            # For inefficient DMU's
            normalization[i] = monetaryti[i] ./ timeasure[i]
        end
    end

    if monetary
        refficiency = rimeasure
        techefficiency = timeasure .* normalization
        allocefficiency = ritarget
    else
        refficiency = rimeasure ./ normalization
        techefficiency = timeasure
        allocefficiency = ritarget ./ normalization
    end

    return RevenueGDADEAModel(n, m, s, measure, rts, monetary, names, refficiency, rlambdaeff, techefficiency, allocefficiency, normalization, Xtarget, Ytarget)
end

function Base.show(io::IO, x::RevenueGDADEAModel)
    compact = get(io, :compact, false)

    n = nobs(x)
    m = ninputs(x)
    s = noutputs(x)
    dmunames = names(x)

    eff = efficiency(x)
    techeff = efficiency(x, :Technical)
    alloceff = efficiency(x, :Allocative)

    if !compact
        print(io, "General Direct Approach Revenue DEA Model \n")
        print(io, "DMUs = ", n)
        print(io, "; Inputs = ", m)
        print(io, "; Outputs = ", s)
        print(io, "\n")
        print(io, "Returns to Scale = ", string(x.rts))
        print(io, "\n")
        print(io, "Associated efficiency measure = ", string(x.measure))
        print(io, "\n")
        show(io, CoefTable(hcat(eff, techeff, alloceff), ["Revenue", "Technical", "Allocative"], dmunames))
    end

end


text:
# dbn_distributions.jl
# David Merrell
# 2019-11
#
# We define some probability distribution objects:
#     * prior distributions for parent sets
#     * a marginal likelihood function for Gaussian DBNs

using LRUCache
using Combinatorics
using LinearAlgebra
import Distributions: Bernoulli
import Gen: Distribution, random, logpdf

###########################################
# PARENT SET PRIOR DISTRIBUTIONS
###########################################

# This distribution assumes prior knowledge in the form of
# "reference edges" -- roughly speaking, it assigns uniform probability
# to all subsets of the reference parent set.
struct RefParentPrior <: Distribution{Vector{Int64}} end
const refparentprior = RefParentPrior()

refparentprior(V::Int64,
               ref_parents::Vector{Int64},
               lambda::Float64) = random(refparentprior,
                                         V, ref_parents, lambda)

function random(pp::RefParentPrior, V::Int64,
                                    ref_parents::Vector{Int64},
                                    lambda::Float64)
    enl = exp(-lambda)
    p = enl /(1.0+enl)
    parents = Vector{Int64}()
    for v=1:V
        if v in ref_parents
            if Bool(rand(Bernoulli(0.5)))
                push!(parents, v)
            end
        else
            if Bool(rand(Bernoulli(p)))
                push!(parents, v)
            end
        end
    end

    return parents
end

function logpdf(pp::RefParentPrior, x_parents, V, ref_parents, lambda)

    lp = -log(1.0 + exp(-lambda))
    lhalf = -log(2)
    n_ref = length(ref_parents)
    total = n_ref*lhalf + (V - n_ref)*lp

    n_extra = length(setdiff(x_parents, ref_parents))
    total += n_extra*(-lambda + lp - lhalf)

    return total
end

# This distribution assumes prior knowledge in the form of
# edge confidences -- it's a "smooth" variant of the
# RefParentPrior distribution.
struct ConfParentPrior <: Distribution{Vector{Int64}} end
const confparentprior = ConfParentPrior()

confparentprior(idx, V, parent_confs::Dict,
                lambda::Float64) = random(confparentprior, idx, V, parent_confs, lambda)

function random(pp::ConfParentPrior, idx, V, parent_confs::Dict, lambda::Float64)

    enl = exp(-lambda)

    parents = Vector{Int64}()
    for i=1:V
        conf = get(parent_confs, i, 0.0)
        p = enl/(enl + exp(-conf * lambda))
        if Bool(rand(Bernoulli(p)))
            push!(parents, i)
        end
    end

    return parents
end

# We use LRU caching to reduce redundant computation.
TENMB = 10000000
global prior_cache = LRU{Vector{Int64}, Float64}(maxsize=10*TENMB, by=Base.summarysize)

function logpdf(pp::ConfParentPrior, x_parents::Vector{Int64},
				     idx::Int64,
				     V::Int64,
                                     parent_confs::Dict,
                                     lambda::Float64)

    lp = get!(prior_cache, pushfirst!(copy(x_parents), idx)) do

        enl = exp(-lambda)
        ref_nonzero = keys(parent_confs)
        lp = 0.0

        # If the number of nonzero confidences is large,
        # then we use array computations
        if length(ref_nonzero) > 50

            # All zero-confidence vertices
            lp += -log(1.0 + enl) * (V - length(parent_confs))

            # All nonzero-confidence vertices
            npl = map(x-> -lambda * x, values(parent_confs))
            lp += sum(map(x-> x - log(enl + exp(x)), npl))

            # parents
            for parent in x_parents
                # nonzero-confidence parents
                if haskey(parent_confs, parent)
                    lp += -lambda * (1.0 - parent_confs[parent])
                # zero-confidence parents
                else
                    lp += -lambda
                end
            end

        # Otherwise, we operate on the dictionaries
        else
            # parents
            for parent in x_parents
                conf = get(parent_confs, parent, 0.0)
                lp += -lambda - log(enl + exp(-conf*lambda))
            end

            # nonparents AND ref_nonzero
            for nonparent in setdiff(ref_nonzero, x_parents)
                npl = -parent_confs[nonparent]*lambda
                lp += npl - log(enl + exp(npl))
            end

            # nonparents AND NOT ref_nonzero
            zero_n_lp = - log(1.0 + enl)
            lp += zero_n_lp*(V - length(union(ref_nonzero, x_parents)))
        end

        return lp
    end

    return lp
end

###########################################
# DBN MARGINAL LIKELIHOOD HELPERS
###########################################

# We use LRU caching to reduce redundant computation.
global lml_cache = LRU{Vector{Int64}, Float64}(maxsize=10*TENMB, by=Base.summarysize)

function clear_caches()
    empty!(prior_cache)
    empty!(lml_cache)
end

function compute_B_col(inds::Vector{Int64}, Xminus)

    # Multiply the columns together
    col = prod(Xminus[:,inds], dims=2)[:,1]

    # The result ought to be standardized.
    col = (col .- mean(col)) ./std(col)
    return col
end

function construct_B(parent_inds::Vector{Int64}, Xminus::Array{Float64,2}, regression_deg::Int64)

    n_parents = length(parent_inds)
    n_cols = sum([binomial(n_parents, k) for k=1:min(regression_deg, n_parents)])
    B = zeros(size(Xminus)[1], n_cols)

    col = 1

    for deg=1:regression_deg
        for comb in Combinatorics.combinations(parent_inds, deg)
            B[:,col] = compute_B_col(comb, Xminus)
            col += 1
        end
    end
    return B
end

"""
Wrap Julia's linear solver so that singular matrices have
small "ridge" terms added to them until they become nonsingular.
"""
function safe_solve(A, y; step=1e-6, max_steps=5)
    try
        return A\y
    catch e
        if isa(e, LinearAlgebra.SingularException)
            AA = transpose(A) * A + step*I
            Ay = transpose(A) * y
            steps = 1
            while steps <= max_steps
                try
                    solved = AA\Ay
                    return solved
                catch ee
                    if isa(ee, LinearAlgebra.SingularException)
                        AA += step*I
                        steps += 1
                        continue
                    else
                        throw(ee)
                    end
                end
            end
        else
            throw(e)
        end

    end
    throw(LinearAlgebra.SingularException)
end

"""

"""
function compute_lml(parent_inds::Vector{Int64},
		     Xminus::Array{Float64,2}, Xp::Vector{Float64},
		     regression_deg::Int64)

    B = construct_B(parent_inds, Xminus, regression_deg)
    m = length(parent_inds)
    Bwidth = size(B,2)

    n = size(Xp)[1]

    ip = dot(Xp, Xp - (n/(n+1.0))*B*safe_solve(B, Xp))
    return -0.5*Bwidth*log(1.0 + n) - 0.5*n*log(ip)

end

"""
This log marginal likelihood is used by the "vertex-specific" model
formulation. `Xp` is a single column of forward-time data values.
"""
function log_marg_lik(ind::Int64, parent_inds::Vector{Int64},
		      Xminus::Array{Float64,2}, Xp::Vector{Float64},
		      regression_deg::Int64)

    lml = get!(lml_cache, pushfirst!(copy(parent_inds), ind)) do
        compute_lml(parent_inds, Xminus, Xp, regression_deg)
    end

    return lml
end

#######################
# END HELPER FUNCTIONS
#######################

##################################################
# Marginal likelihood
##################################################

"""
    CPDMarginal

Implementation of the regression marginal likelihood
described in Hill et al. 2012:

P(X|G) ~ (1.0 + n)^(-Bwidth/2) * ( X+^T X+ - n/(n+1) * X+^T B2invconj X+)^(n/2)

This is expensive to compute. In order to reduce redundant
comptuation, we use multiple levels of caching.
"""
struct CPDMarginal <: Distribution{Vector{Float64}} end
const cpdmarginal = CPDMarginal()

cpdmarginal(Xminus, ind, parents, regression_deg) = random(cpdmarginal,
							   Xminus,
							   ind, parents,
							   regression_deg)

"""
The random function does nothing -- we will always observe these values.
"""
function random(cpdm::CPDMarginal, Xminus::Array{Float64,2},
		ind::Int64, parents::Vector{Int64}, regression_deg::Int64)
    return zeros(size(Xminus)[1])
end

function logpdf(cpdm::CPDMarginal, Xp::Vector{Float64}, Xminus::Array{Float64,2},
		ind::Int64, parents::Vector{Int64}, regression_deg::Int64)
    return log_marg_lik(ind, parents, Xminus, Xp, regression_deg)
end

max_stars_repo_name: gitter-lab/ssps

max_issues_repo_name: gitter-lab/ssps

max_forks_repo_name: gitter-lab/ssps


text:
export
    TransformFilter

import .Ranger.Nodes.calc_transform!

# By default the TransformFilter will exclude or "block" rotations and scales from
# propagating to the children, but allow translations.
#
mutable struct TransformFilter <: Ranger.AbstractNode
    base::Nodes.NodeData
    transform::Nodes.TransformProperties{Float64}

    # Collection of nodes.
    children::Array{Ranger.AbstractNode,1}

    # This node's parent translation component
    components::Math.AffineTransform{Float64}

    # What to exclude from the immediate parent
    exclude_translation::Bool
    exclude_rotation::Bool
    exclude_scale::Bool

    function TransformFilter(world::Ranger.World, name::String, parent::Ranger.AbstractNode)
        o = new()

        o.base = Nodes.NodeData(name, parent, world)
        o.transform = Nodes.TransformProperties{Float64}()
        o.children = Array{Ranger.AbstractNode,1}[]

        o.components = Math.AffineTransform{Float64}()

        o.exclude_translation = false
        o.exclude_rotation = true
        o.exclude_scale = true

        o
    end
end

# Filters are special in that they overload the visit() method
function Nodes.visit(node::TransformFilter, context::Rendering.RenderContext, interpolation::Float64)
    if !Nodes.is_visible(node)
        return;
    end

    Rendering.save!(context)

    children = Nodes.get_children(node)
    if children ≠ nothing
        for child in children
            Rendering.save!(context)

            if Nodes.has_parent(node)
                parent = node.base.parent

                # Re-introduce only the parent's components as defined by
                # exclusion flags.
                Nodes.calc_filtered_transform!(parent.transform,
                    node.exclude_translation, node.exclude_rotation, node.exclude_scale,
                    node.components)

                # Combine using pre-multiply
                # "parent.transform.inverse" removes the immediate parent's transform effects
                Math.multiply!(node.components, parent.transform.inverse, node.transform.aft)

                # Merge them with the current context.
                Rendering.apply!(context, node.transform.aft)
            else
                println("TransformFilter::visit: ", node, " has NO parent")
                return;
            end

            # Now visit the child with the modified context
            Nodes.visit(child, context, interpolation)

            Rendering.restore!(context)
        end
    end

    Rendering.restore!(context)
end

# Filter nodes calc their own transform on behalf of their children
function Nodes.calc_transform!(node::TransformFilter)
    node.transform.aft
end

# --------------------------------------------------------
# Grouping
# --------------------------------------------------------
function Nodes.get_children(node::TransformFilter)
    node.children
end

# Older version that doesn't set aft
# function Nodes.visit(node::TransformFilter, context::Rendering.RenderContext, interpolation::Float64)
#     if !Nodes.is_visible(node)
#         return;
#     end

#     Rendering.save!(context)

#     children = Nodes.get_children(node)
#     if children ≠ nothing
#         for child in children
#             Rendering.save!(context)

#             if Nodes.has_parent(node)
#                 parent = node.base.parent

#                 # This removes the immediate parent's transform effects
#                 Rendering.apply!(context, parent.transform.inverse)

#                 # Re-introduce only the parent's components as defined by
#                 # exclusion flags.
#                 Nodes.calc_filtered_transform!(parent.transform,
#                     node.exclude_translation, node.exclude_rotation, node.exclude_scale,
#                     node.components)

#                 # Merge them with the current context.
#                 Rendering.apply!(context, node.components)
#             else
#                 println("TransformFilter::visit: ", node, " has NO parent")
#                 return;
#             end

#             # Now visit the child with the modified context
#             Nodes.visit(child, context, interpolation)

#             Rendering.restore!(context)
#         end
#     end

#     Rendering.restore!(context)
# end

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


text:
using Base.Test

include("./List.jl")
include("./Adwin.jl")

#-- start test counter
a = ListIt()
b = ListIt()
c = ListIt()
@test myCounter(a) == 1
a.next = b
b.previous = a
@test myCounter(a) == 2
a.next.next = c
c.previous = b
@test myCounter(a) == 3
#-- end test counter

L = List()
for i = 1:4
  addToHead!(L)
end

removeFromHead!(L)
@test myCounter(L.head) == 4
removeFromHead!(L)
@test myCounter(L.head) == 3
#
removeFromTail!(L)
# @test myCounter(L.head) == 2   ----------- todo FIX test case
run_adwin()

# addToTail!(L)
# @test myCounter(L.head) == 3

#################### end
# @test myCounter(L.head) == 3
#removeFromHead!(L)

#println("L.count = ", L.count)
#clearLst!(L)
#println("L.count = ", L.count)
#addToTail!(L)
#addToTail!(L)
#println("L.count = ", L.count)
#removeFromTail!(L)
#println("L.count = ", L.count)
#removeFromTail!(L)
#println("L.count = ", L.count)

#addToHead!(L, item)
#println("L.count = ", L.count)
#println(isnull(item.next))

#mutable struct ListItem
#    id::Int
#    next::Nullable{ListItem}
#    ListItem(idnum) = (x = new();
#                       x.id=idnum;
#                       x.next = Nullable{ListItem}();
#                       x)
#end
#
#a = ListItem(1)
#b = ListItem(2)
#c = ListItem(3)
#a.next = b
#b.next = c
#
#println(isnull(a.next))
#println(isnull(b.next))
#println(isnull(c.next))
#

max_stars_repo_name: av-maslov/Recurrency

max_issues_repo_name: av-maslov/Recurrency

max_forks_repo_name: av-maslov/Recurrency


text:
Base.@kwdef mutable struct GroupLassoRegressor{G,P,T<:Number} <: AbstractGroupRegressor
    decomposition::Symbol = :default
    groups::G
    groups_multiplier::P = sqrt.(groups.ps) ./ sqrt(groups.p)
    λ::T = 1.0
    center::Bool = true
    scale::Bool = true
    maxiter::Int = 100
    η_reg::T = 1e-7
    η_threshold::T = 1e-2
    abs_tol::T = 1e-4
    truncate_to_zero::Bool = true
end

_main_hyperparameter(::GroupLassoRegressor) = :λ

function _default_hyperparameter_maximum(model::GroupLassoRegressor, fitted_machine)
    @unpack groups, groups_multiplier = model
    _norms = group_summary(groups, _workspace(fitted_machine.cache).XtY, norm)
    maximum(_norms ./ groups_multiplier)
end

_default_param_min_ratio(::GroupLassoRegressor, fitted_machine) = 1e-5

function _glasso_fit!(workspace, glasso::GroupLassoRegressor)
    @unpack η_reg, η_threshold, abs_tol, groups, maxiter, λ, groups_multiplier, truncate_to_zero = glasso

    tmp_λs = copy(workspace.λs)
    ηs_new = group_summary(groups, StatsBase.coef(workspace), norm)
    ηs_old = copy(ηs_new)

    converged = false
    iter_cnt = 0
    for i = 1:maxiter
        tmp_λs .= λ .* groups_multiplier ./ sqrt.(abs2.(ηs_new) .+ η_reg)
        fit!(workspace, tmp_λs)
        ηs_new .= group_summary(groups, StatsBase.coef(workspace), norm)
        #converged = norm(ηs_new .- ηs_old, Inf) < abs_tol
        #@show (ηs_new .- ηs_old) ./ sqrt.(abs2.(ηs_old) .+ η_reg)
        #@show (ηs_new .- ηs_old)
        #@show sqrt.( abs2.(ηs_old) .+ η_reg)
        converged =
            (norm((ηs_new .- ηs_old) ./ sqrt.(abs2.(ηs_old) .+ η_reg), Inf) < abs_tol) ||
            (norm(ηs_new .- ηs_old, Inf) < abs_tol)
        ηs_old .= ηs_new
        iter_cnt += 1
        converged && break
    end
    #@show "conv"
    ηs = group_summary(groups, StatsBase.coef(workspace), norm)
    final_λs = deepcopy(workspace.λs)
    #zero_groups = group_summary(groups, StatsBase.coef(workspace), norm) .< η_threshold .* groups_multiplier
    #final_λs[zero_groups] .= Inf
    #fit!(workspace, final_λs)

    (workspace = workspace, converged = converged, iter_count = iter_cnt)
end

function MMI.fit(m::GroupLassoRegressor, verb::Int, X, y)
    @unpack decomposition, center, scale = m
    Xmatrix = MMI.matrix(X)
    p = size(Xmatrix, 2)
    m_tmp = MultiGroupRidgeRegressor(;
        groups = m.groups,
        decomposition = decomposition,
        scale = scale,
        center = center,
    )
    multiridge_machine = MLJ.machine(m_tmp, X, y)
    fit!(multiridge_machine)
    workspace = _workspace(multiridge_machine.cache)
    glasso_workspace = _glasso_fit!(workspace, m)
    βs = StatsBase.coef(glasso_workspace.workspace)
    x_transform = multiridge_machine.fitresult.x_transform
    y_transform = multiridge_machine.fitresult.y_transform
    fitresult = (coef = βs, x_transform = x_transform, y_transform = y_transform)

    # return
    return fitresult, glasso_workspace, NamedTuple{}()
end

function MMI.update(
    model::GroupLassoRegressor,
    verbosity::Int,
    old_fitresult,
    old_cache,
    X,
    y,
)
    glasso_workspace = _glasso_fit!(old_cache.workspace, model)
    βs = StatsBase.coef(glasso_workspace.workspace)

    x_transform = old_fitresult.x_transform
    y_transform = old_fitresult.y_transform
    fitresult = (coef = βs, x_transform = x_transform, y_transform = y_transform)

    return fitresult, glasso_workspace, NamedTuple{}()
end


text:
using PyCall
using CSV
using CodecZlib, Mmap
using JLD2

const DEFAULT_FEATURE_COLS = [:ensembleid, :genesymbol, :type]
const DEFAULT_BARCODE_COLS = [:barcode]
const FEATURE_COLS = [:featurekey, :featurename, :featuretype, :chromosome, :featurestart, :featureend, :isgene, :genus_species]

function read_mtx(filename::String)
    py"""
    from scipy.io import mmread
    mtx = mmread($filename)
    """
    sparse(py"mtx.row" .+ 1, py"mtx.col" .+ 1, py"mtx.data")
end

function read_tsv_gz(filename::String, header::Vector{Symbol}=Symbol[])
    gz_unzipper = transcode(GzipDecompressor, Mmap.mmap(filename))
    if isempty(header)
        return CSV.File(gz_unzipper; delim='\t') |> DataFrame
    else
        return CSV.File(gz_unzipper; delim='\t', header=header) |> DataFrame
    end
end

read_features(filename::String, names=DEFAULT_FEATURE_COLS) = read_tsv_gz(filename, names)
read_barcodes(filename::String, names=DEFAULT_BARCODE_COLS) = read_tsv_gz(filename, names)

## Human Cell Atlas (HCA)
read_genes(filename::String) = read_tsv_gz(filename)
read_cells(filename::String) = read_tsv_gz(filename)