text:
using JuMP, Gurobi, Distributions, StatsFuns, StatsBase

include("normopt_noniid_initialization_template.jl")

const gurobi_env = Gurobi.Env()

const numVariables = 100
const numJCC = 100
const conRHS = (numVariables*1.0)^2

xi_mean = zeros(Float64,numVariables)
xi_covariance = 0.5*ones(Float64,numJCC,numJCC)

for j = 1:numVariables
	xi_mean[j] = j*1.0/numVariables
end

for i = 1:numJCC
	xi_covariance[i,i] = 1.0
end

xi_cov_chol_tmp = cholfact(xi_covariance)
xi_cov_chol = xi_cov_chol_tmp[:L]

xi_mean_mat = zeros(Float64,numJCC,numVariables)
for j = 1:numVariables
	xi_mean_mat[:,j] = xi_mean[j]*ones(numJCC)
end

# tailored implementation of multivariate normal distribution
function myMvNormal(numSamples::Int64)

	xi = zeros(Float64,numJCC,numVariables,numSamples)

	xi_tmp = rand(Normal(0.0,1.0),numJCC,numVariables,numSamples)
	for samp = 1:numSamples
		xi[:,:,samp] = xi_cov_chol*xi_tmp[:,:,samp] + xi_mean_mat
	end

	return xi
end

# generate random samples from multivariate normal distribution
function generateRandomSamples(numSamples::Int64)

	xi = myMvNormal(numSamples)

	return xi
end

srand(1234)

# generate samples to estimate actual risk level of solution
const numAnalyticalSamples = 20000
const xi_analytical = generateRandomSamples(numAnalyticalSamples)

srand()

# compute negative probability of all constraints being satisfied
# for given scenarios of demands d and a given decision vector x
function computeRiskLevel(x::Array{Float64},reliabilityLevel::Float64=1E-06)

	const numSamples = size(xi_analytical,3)
	simpleEst::Bool = true

	numViolated::Int64 = checkScenarioConstraints(x,xi_analytical)

	riskLevel::Float64 = Inf

	if(simpleEst)
		riskLevel = numViolated*1.0/numSamples
		return riskLevel
	else
		if(numViolated == 0)
			riskLevel = -log(reliabilityLevel)/numSamples
			return riskLevel
		end

		 # termination criterion for bisection
		gamma_tolerance::Float64 = 0.1/numSamples

		# use tail bounds on binomial distribution to estimate lower and upper bounds
		gamma_low::Float64 = (numViolated - sqrt(-numSamples*log(reliabilityLevel)/2.0))/numSamples
		if(gamma_low < 0.0)
			gamma_low = 0.0
		end
		gamma_up::Float64 = (numViolated + sqrt(-numSamples*log(reliabilityLevel)/2.0))/numSamples
		if(gamma_up > 1.0)
			gamma_up = 1.0
		end

		gamma::Float64 = 0.0
		func_value::Float64 = 0.0

		# use bisection to solve the nonlinear equation
		while (gamma_up - gamma_low > gamma_tolerance)

			gamma = (gamma_low + gamma_up)/2.0
			func_value = -reliabilityLevel
			for i = 1:numViolated
				tmp_func = i*log(gamma) + (numSamples-i)*log(1-gamma) + lfact(numSamples) - lfact(numSamples-i) - lfact(i)
				func_value += exp(tmp_func)
				if(func_value > 0)
					break
				end
			end
			if(func_value > 0)
				gamma_low = (gamma_low + gamma_up)/2.0
			else
				gamma_up = (gamma_low + gamma_up)/2.0
			end

		end

		riskLevel = gamma_up
		return riskLevel
	end

end

# determines if all constraints are satisfied for a given single scenario
# returns one if scenario constraints are all satisfied, zero otherwise
function checkScenarioConstraints(x::Array{Float64},xi::Array{Float64})

	numSamples::Int64 = size(xi,3)

	numViolated::Int64 = 0
	for samp = 1:numSamples
		for i = 1:numJCC
			con::Float64 = evaluateConstraint(i,x,xi[i,:,samp])
			if(con > 0)
				numViolated += 1
				break
			end
		end
	end

	return numViolated
end

# evaluate the constraint values for given decision vector x,
# realization of the random variables xi, and constraint scalings
function evaluateConstraint(x::Array{Float64},xi::Array{Float64},scalings::Array{Float64})

	con = (sum((xi.*x').^2,2) - conRHS*ones(numJCC))./scalings

	return con
end

# evaluate the constraint values for given decision vector x,
# realization of the random variables xi, and constraint scalings
function evaluateConstraint(x::Array{Float64},xi::Array{Float64})

	con = evaluateConstraint(x,xi,ones(numJCC))

	return con
end

# evaluate the constraint values for given decision vector x,
# realization of the random variables xi, and constraint scalings
function evaluateConstraint(index::Int64,x::Array{Float64},xi::Array{Float64})

	con::Float64 = norm(xi.*x)^2 - conRHS

	return con
end

# evaluate constraint gradients
function evaluateConstraintGradients(x::Array{Float64},xi::Array{Float64})

	return evaluateConstraintGradients(x,xi,ones(numJCC))
end

# evaluate constraint gradients
function evaluateConstraintGradients(x::Array{Float64},xi::Array{Float64},scalings::Array{Float64})

	grad_con = zeros(Float64,numJCC,numVariables)

	for i = 1:numJCC
		grad_con[i,:] = (2*(xi[i,:]).^2.*x)/scalings[i]
	end

	return grad_con
end

# evaluate constraint gradients
function evaluateConstraintGradients(x::Array{Float64},xi::Array{Float64},scalings::Float64)

	grad_con = (2*(xi.^2).*x)/scalings

	return grad_con
end

# project the current iterate onto the feasible set
function projectOntoFeasibleSet(x_curr::Array{Float64},objectiveBound::Float64,proj_params::Array{Float64})

	mod = Model(solver=GurobiSolver(gurobi_env,Presolve=0,OutputFlag=0))

	@variable(mod, y[j=1:numVariables] >= 0)

	@objective(mod, Min, sum((y[i]-x_curr[i])^2 for i = 1:numVariables))
	@constraint(mod, -sum(y[j] for j=1:numVariables) <= objectiveBound)

	solve(mod)
	x_proj = getvalue(y)

	return x_proj, proj_params
end

# get the stochastic gradient of the approximation
function getStochasticGradient(x::Array{Float64},numSamples::Int64,mu::Float64,tau::Float64,scalings::Array{Float64})

	xi = generateRandomSamples(numSamples)

	stochgrad = zeros(Float64,numVariables)

	for samp = 1:numSamples
		con = evaluateConstraint(x,xi[:,:,samp],scalings)

		max_info = findmax(con)
		max_con::Float64 = max_info[1]
		max_index::Int64 = max_info[2]

		baseFactor1::Float64 = tau/(sqrt(mu)*exp(0.5*tau*max_con) + sqrt(1.0/mu)*exp(-0.5*tau*max_con))^2

		grad_max_con = evaluateConstraintGradients(x,xi[max_index,:,samp],scalings[max_index])

		stochgrad += baseFactor1*grad_max_con/numSamples
	end

	return stochgrad

end

# get the stochastic gradient of the approximation for fixed random sample
function getStochasticGradientForFixedSample(x::Array{Float64},xi::SubArray{Float64,3,Array{Float64,3},Tuple{Base.Slice{Base.OneTo{Int64}},Base.Slice{Base.OneTo{Int64}},UnitRange{Int64}},true},mu::Float64,tau::Float64,scalings::Array{Float64})

	numSamples::Int64 = size(xi,3)

	stochgrad = zeros(Float64,numJCC,numVariables)

	for samp = 1:numSamples

		con = evaluateConstraint(x,xi[:,:,samp],scalings)

		grad_con = evaluateConstraintGradients(x,xi[:,:,samp],scalings)

		for i = 1:numJCC
			baseFactor1::Float64 = tau/(sqrt(mu)*exp(0.5*tau*con[i]) + sqrt(1.0/mu)*exp(-0.5*tau*con[i]))^2

			stochgrad[i,:] += baseFactor1*grad_con[i,:]/numSamples
		end
	end

	return stochgrad

end

# estimate constraint scalings
function estimateScalings(x::Array{Float64},numSamples::Int64)

	scaling_tol::Float64 = 1E-06

	xi = generateRandomSamples(numSamples)

	constraint_values = zeros(Float64,numJCC,numSamples)

	for samp = 1:numSamples
		constraint_values[:,samp] = evaluateConstraint(x,xi[:,:,samp],ones(numJCC))
	end

	scalings = zeros(Float64,numJCC)
	for index = 1:numJCC
		scalings[index] = max(median(abs.(constraint_values[index,:])),scaling_tol)
	end

	return scalings
end

# estimate Lipschitz constant of the gradient of the approximation
function estimateCompositeLipschitzConstant(x_ref::Array{Float64},numSamplesForLipVarEst::Int64,numGradientSamples::Int64,batchSize::Int64,mu::Float64,tau::Float64,scalings::Array{Float64},objectiveBound::Float64)

	proj_params = zeros(Float64,3)

	sample_radius::Float64 = norm(x_ref)/10.0

	xi = generateRandomSamples(batchSize)

	L_avg::Float64 = 0.0
	for batch = 1:batchSize
		xi_start::Int64 = batch
		xi_end::Int64 = batch

		L_max::Float64 = 0.0
		for samp = 1:numSamplesForLipVarEst
			x_1 = pickPointOnSphere(x_ref,sample_radius)
			x_1,~ = projectOntoFeasibleSet(x_1,objectiveBound,proj_params)

			x_2 = pickPointOnSphere(x_ref,sample_radius)
			x_2,~ = projectOntoFeasibleSet(x_2,objectiveBound,proj_params)

			@views grad_1 = getStochasticGradientForFixedSample(x_1,xi[:,:,xi_start:xi_end],mu,tau,scalings)
			@views grad_2 = getStochasticGradientForFixedSample(x_2,xi[:,:,xi_start:xi_end],mu,tau,scalings)

			L_est::Float64 = norm(grad_2 - grad_1)/norm(x_2 - x_1)
			if(L_est > L_max)
				L_max = L_est
			end
		end

		L_avg += L_max/batchSize
	end

	return L_avg
end

# estimate step length factor
function estimateStepLength(x_ref::Array{Float64},numGradientSamples::Int64,
							numSamplesForLipVarEst::Int64,batchSize::Int64,maxNumIterations::Int64,
							minNumReplicates::Int64,mu::Float64,tau::Float64,
							scalings::Array{Float64},objectiveBound::Float64)

	proj_params = zeros(Float64,3)

	sample_radius::Float64 = norm(x_ref)/10.0

	rho_est::Float64 = estimateCompositeLipschitzConstant(x_ref,numSamplesForLipVarEst,numGradientSamples,batchSize,mu,tau,scalings,objectiveBound)

	variance_est::Float64 = 0.0
	for samp = 1:numSamplesForLipVarEst
		x_1 = pickPointOnSphere(x_ref,sample_radius)
		x_1,~ = projectOntoFeasibleSet(x_1,objectiveBound,proj_params)

		variance_1::Float64 = 0.0
		for batch = 1:batchSize
			grad_1 = getStochasticGradient(x_1,numGradientSamples,mu,tau,scalings)

			variance_1 += norm(grad_1)^2/batchSize
		end

		if(variance_1 > variance_est)
			variance_est = variance_1
		end
	end

	gamma::Float64 = 1/sqrt(rho_est*variance_est)

	stepLength::Float64 = gamma/sqrt((maxNumIterations+1.0)*minNumReplicates)

	return stepLength
end

# given a reference point, pick a random point within a given radius
function pickPointOnSphere(x_orig::Array{Float64},radius::Float64)

	if(radius < 1E-09)
		error("Sampling radius in pickPointOnSphere is almost zero! Set manually.")
	end

	n::Int64 = size(x_orig,1)
	unifRand = rand(1)
	factor1::Float64 = radius*(unifRand[1])^(1.0/n)
	normRand = rand(Normal(0,1), n)
	factor2::Float64 = norm(normRand)

	x_samp = x_orig + factor1*normRand/factor2

	return x_samp
end

max_stars_repo_name: rohitkannan/SA-for-CCP

max_issues_repo_name: rohitkannan/SA-for-CCP

max_forks_repo_name: rohitkannan/SA-for-CCP


text:
export projectGradient, line_search, projected_newton

# Write projected Newton method, according to
# Bertsekas, Dimitri P. "Projected Newton methods for optimization problems with simple constraints." SIAM Journal on control and Optimization 20.2 (1982): 221-246.

# Code inspired from Ricardo Baptista, Youssef Marzouk MIT
function projectGradient(x::Array{Float64,1}, G::Array{Float64,1})
    PG = copy(G);
    for i=1:size(PG,1)
        @inbounds begin
        if x[i]==0.0 && G[i] >=0
        PG[i] = 0.0
        end
        end
    end
    return PG
end

function linesearch(x::Array{Float64,1}, G::Array{Float64,1},
         p::Array{Float64,1}, J::Float64,Idx::Array{Int64,1}, Lhd::LHD)

    # Fixed parameters
    itmax = 15
    σ = 1e-4
    β = 2

    it = 0
    Jx_α_lin = copy(J)+1 # To enter the while loop
    α = copy(β)# The first iteration inside the while loop will be for α = 1
    α_p = copy(p)
    x_α = zero(x)
    while (J<Jx_α_lin) && (it<itmax)
        α /= β
        α_p .= copy(p)
        rmul!(α_p, copy(α))
        x_α .= max.(0.0, x - α_p)
        Jx_α = copy(Lhd(x_α, false, noutput = 1))
        α_p[Idx] .= x[Idx] - x_α[Idx]
        Jx_α_lin = Jx_α + σ*dot(G, α_p)
        it += 1
    end

    if (J<Jx_α_lin)
        print("Line search method reached maximum number of iterations")
        xh = copy(x)
        xh[Idx] .= 0.0
        index = Int64[]

        for i=1:size(x,1)
            @inbounds begin
            if (xh[i]>0.0) && (p[i]>0.0)
                push!(index, copy(i))
            end
            end
        end
        if isempty(index)
            α = 1.0
        else
            α = min(x[index]./p[index])
        end
    end
return α
end

function projected_newton(x0::Array{Float64,1}, Lhd::LHD, type::String)
# Default parameters
rtol_J = 1e-6;
rtol_G = 1e-6;
itmax = 30;
ϵ = 0.01;

n = size(x0,1)

for i=1:n
    @assert sign(x0[i])==1.0 "Initial conditions are not feasible"
end

x = copy(x0)
Lhd(x, true)
@get Lhd (G,H)

Jold =copy(Lhd.J[1])
J = copy(Lhd.J[1])

norm_PG0 = norm(projectGradient(x, G))
tol_G = norm_PG0*rtol_G

norm_PG = copy((norm_PG0))
rδJ = rtol_J + 1
it = 0

# Pre-allocate space
p =zero(x)
# h = zero(x)
z = zero(x)

## Iteration
while (rδJ > rtol_J) && (norm_PG > tol_G) && (it<itmax)
    # Define search direction
    wk = norm(x - max.(0.0, x-G))
    ϵk = min(ϵ, wk)
    Idx= Int64[]
    @inbounds for i=1:n
        if x[i] <= ϵk && G[i] >0.0
        push!(Idx, copy(i))
        end
    end

    if isempty(Idx)==false
        z .= zeros(n)
        @inbounds for i in Idx
            z[i] = copy(H[i,i])
        end
        # h .= deepcopy(diag(H))
        # z .= zeros(n)
        # z[Idx] .= h[Idx]
        H[Idx,:] .= 0.0
        H[:,Idx] .= 0.0
        H .+= diagm(z)
    end

    # Switch according to search direction
    if type =="TrueHessian"
    p .= copy(G)
    # Perform in-place Cholesky decomposition and store H\G into p
    LAPACK.posv!('U', H, p)
    elseif type =="ModifHessian"
        itmax = 100
        # Try if the TrueHessian technique works
        if isposdef(H)
            p .= copy(G)
            # Perform in-place Cholesky decomposition and store H\G into p
            LAPACK.posv!('U', H, p)
        else
        # else Strategy 1- flip negative eigenvalues
            if !issymmetric(H)
                H .+= H'
                rmul!(H, 0.5)
            end

            print("Hessian H is not positive-definite")
            # EH = eigen(Symmetric(H))
            EH = eigen(H)
            for i=1:size(EH.values,1)
                @inbounds EH.values[i] = copy( 1/(max(1e-8, abs(EH.values[i]))) )
            end

            # p .= deepcopy(G)
            mul!(p, EH.vectors', G)
            p .*= EH.values
            mul!(p, EH.vectors, copy(p))

            # p .= deepcopy(EH.vectors*( EH.values .* (EH.vectors'*G)))
        end
    elseif type =="Gradient"
        itmax = 1000
        p = copy(G)

    else
        error("Method not implemented yet")
    end

    # Perform line search
    α = linesearch(x,G,p,J,Idx,Lhd)
    # @show α

    # Update the estimate
    x .= copy(max.(0.0, x - α*p))
    # @show x
    Lhd(x, true)
    @get Lhd (G,H)
    J = copy(Lhd.J[1])

    # Convergence criteria
    rδJ = abs(J-Jold)/abs(Jold)
    # @show rδJ
    Jold = copy(J)
    norm_PG = copy(norm(projectGradient(x,G)))
    it +=1
end

if it>=itmax
    print("Max number of iterations is reached during the optimization")
end

return x, rδJ, norm_PG, it

end

ext: jl


text:
import DataStructures: OrderedDict

"""
```
RealBondMkup{T} <: AbstractModel{T}
```

### Fields

#### Parameters and Steady-States
* `parameters::Vector{AbstractParameter}`: Vector of all time-invariant model
  parameters.

* `steady_state::Vector{AbstractParameter}`: Model steady-state values, computed
  as a function of elements of `parameters`.

* `keys::OrderedDict{Symbol,Int}`: Maps human-readable names for all model
  parameters and steady-states to their indices in `parameters` and
  `steady_state`.

#### Inputs to Measurement and Equilibrium Condition Equations

The following fields are dictionaries that map human-readable names to row and
column indices in the matrix representations of of the measurement equation and
equilibrium conditions.

* `endogenous_states::OrderedDict{Symbol,Int}`: Maps each state to a column in
  the measurement and equilibrium condition matrices.

* `exogenous_shocks::OrderedDict{Symbol,Int}`: Maps each shock to a column in
  the measurement and equilibrium condition matrices.

* `expected_shocks::OrderedDict{Symbol,Int}`: Maps each expected shock to a
  column in the measurement and equilibrium condition matrices.

* `equilibrium_conditions::OrderedDict{Symbol,Int}`: Maps each equlibrium
  condition to a row in the model's equilibrium condition matrices.

* `observables::OrderedDict{Symbol,Int}`: Maps each observable to a row in the
  model's measurement equation matrices.

#### Model Specifications and Settings

* `spec::String`: The model specification identifier, \"an_schorfheide\", cached
  here for filepath computation.

* `subspec::String`: The model subspecification number, indicating that some
  parameters from the original model spec (\"ss0\") are initialized
  differently. Cached here for filepath computation.

* `settings::Dict{Symbol,Setting}`: Settings/flags that affect computation
  without changing the economic or mathematical setup of the model.

* `test_settings::Dict{Symbol,Setting}`: Settings/flags for testing mode

#### Other Fields

* `rng::MersenneTwister`: Random number generator. Can be is seeded to ensure
  reproducibility in algorithms that involve randomness (such as
  Metropolis-Hastings).

* `testing::Bool`: Indicates whether the model is in testing mode. If `true`,
  settings from `m.test_settings` are used in place of those in `m.settings`.

* `observable_mappings::OrderedDict{Symbol,Observable}`: A dictionary that
  stores data sources, series mnemonics, and transformations to/from model units.
  DSGE.jl will fetch data from the Federal Reserve Bank of
  St. Louis's FRED database; all other data must be downloaded by the
  user. See `load_data` and `Observable` for further details.

"""
mutable struct RealBondMkup{T} <: AbstractModel{T}
    parameters::ParameterVector{T}                         # vector of all time-invariant model parameters
    steady_state::ParameterVector{T}                       # model steady-state values

    # Temporary to get it to work. Need to
    # figure out a more flexible way to define
    # "grids" that are not necessarily quadrature
    # grids within the model
    grids::OrderedDict{Symbol,Union{Grid, Vector}}

    keys::OrderedDict{Symbol,Int}                          # human-readable names for all the model
                                                           # parameters and steady-states

    state_variables::Vector{Symbol}                        # Vector of symbols of the state variables
    jump_variables::Vector{Symbol}                         # Vector of symbols of the jump variables
    normalized_model_states::Vector{Symbol}                # All of the distributional model
                                                           # state variables that need to be normalized
    endogenous_states_unnormalized::OrderedDict{Symbol,UnitRange} # Vector of unnormalized
                                                           # ranges of indices
    endogenous_states::OrderedDict{Symbol,UnitRange}       # Vector of ranges corresponding
                                                           # to normalized (post Klein solution) indices
    exogenous_shocks::OrderedDict{Symbol,Int}              #
    expected_shocks::OrderedDict{Symbol,Int}               #
    equilibrium_conditions::OrderedDict{Symbol,UnitRange}  #
    model_states_augmented::OrderedDict{Symbol, Int}
    observables::OrderedDict{Symbol,Int}                   #

    spec::String                                           # Model specification number (eg "m990")
    subspec::String                                        # Model subspecification (eg "ss0")
    settings::Dict{Symbol,Setting}                         # Settings/flags for computation
    test_settings::Dict{Symbol,Setting}                    # Settings/flags for testing mode
    rng::MersenneTwister                                   # Random number generator
    testing::Bool                                          # Whether we are in testing mode or not

    observable_mappings::OrderedDict{Symbol, Observable}
end

description(m::RealBondMkup) = "RealBondMkup, $(m.subspec)"

"""
`init_model_indices!(m::RealBondMkup)`

Arguments:
`m:: RealBondMkup`: a model object

Description:
Initializes indices for all of `m`'s states, shocks, and equilibrium conditions.
"""
function init_model_indices!(m::RealBondMkup)

    # Endogenous states
    endogenous_states = collect([
    # These states corresp. to the following in the original notation
    #    MUP,   LIIP,    ZP,    MONP,    ELLP,  RRP
        :μ′_t, :l_i′_t, :z′_t, :mon′_t, :l′_t, :R′_t,
    #    IIP    WWP    PIP    TTP
        :i′_t, :w′_t, :π′_t, :t′_t])

    # Exogenous shocks
    exogenous_shocks = collect([:z_sh, :mon_sh, :mkp_sh])

    # Equilibrium conditions
    equilibrium_conditions = collect([
        :eq_euler, :eq_kolmogorov_fwd, :eq_market_clearing, :eq_TFP,
        :eq_phillips, :eq_taylor, :eq_fisher, :eq_transfers, :eq_monetary_policy,
        :eq_markup, :eq_lagged_nominal_rate])

    model_states_augmented = [:y_t1]

    # Observables
    observables = keys(m.observable_mappings)

    ########################################################################################
    # Setting indices of endogenous_states and equilibrium conditions manually for now
    nx = get_setting(m, :nx)
    ns = get_setting(m, :ns)
    nxns = nx*ns
    endo = m.endogenous_states_unnormalized
    eqconds = m.equilibrium_conditions

    # State variables
    endo[:μ′_t]   = 1:nxns
    endo[:l_i′_t] = nxns+1:nxns+1
    endo[:z′_t]   = nxns+2:nxns+2
    endo[:mon′_t] = nxns+3:nxns+3
    endo[:mkp′_t] = nxns+4:nxns+4

    # Jump variables
    endo[:l′_t]   = nxns+5:2*nxns+4
    endo[:R′_t]   = 2*nxns+5:2*nxns+5
    endo[:i′_t]   = 2*nxns+6:2*nxns+6
    endo[:w′_t]   = 2*nxns+7:2*nxns+7
    endo[:π′_t]   = 2*nxns+8:2*nxns+8
    endo[:t′_t]   = 2*nxns+9:2*nxns+9

    eqconds[:eq_euler]              = 1:nxns
    eqconds[:eq_kolmogorov_fwd]     = nxns+1:2*nxns
    eqconds[:eq_market_clearing]    = 2*nxns+1:2*nxns+1
    eqconds[:eq_TFP]                = 2*nxns+2:2*nxns+2
    eqconds[:eq_phillips]           = 2*nxns+3:2*nxns+3
    eqconds[:eq_taylor]             = 2*nxns+4:2*nxns+4
    eqconds[:eq_fisher]             = 2*nxns+5:2*nxns+5
    eqconds[:eq_transfers]          = 2*nxns+6:2*nxns+6
    eqconds[:eq_monetary_policy]    = 2*nxns+7:2*nxns+7
    eqconds[:eq_markup]             = 2*nxns+8:2*nxns+8
    eqconds[:eq_lagged_nominal_rate]= 2*nxns+9:2*nxns+9
    ########################################################################################

    m.normalized_model_states = [:μ′_t]
    m.endogenous_states = deepcopy(endo)
    m.state_variables = m.endogenous_states.keys[get_setting(m, :state_indices)]
    m.jump_variables = m.endogenous_states.keys[get_setting(m, :jump_indices)]

    for (i,k) in enumerate(exogenous_shocks);            m.exogenous_shocks[k]            = i end
    for (i,k) in enumerate(model_states_augmented); m.model_states_augmented[k]   = i+length(endogenous_states) end
    for (i,k) in enumerate(observables);                 m.observables[k]                 = i end
end

function RealBondMkup(subspec::String="ss0";
                   custom_settings::Dict{Symbol, Setting} = Dict{Symbol, Setting}(),
                   testing = false)

    # Model-specific specifications
    spec               = "real_bond_mkup"
    subspec            = subspec
    settings           = Dict{Symbol,Setting}()
    test_settings      = Dict{Symbol,Setting}()
    rng                = MersenneTwister(0)

    # initialize empty model
    m = RealBondMkup{Float64}(
            # model parameters and steady state values
            Vector{AbstractParameter{Float64}}(), Vector{Float64}(),
            # grids and keys
            OrderedDict{Symbol,Grid}(), OrderedDict{Symbol,Int}(),

            # normalized_model_states, state_inds, jump_inds
            Vector{Symbol}(), Vector{Symbol}(), Vector{Symbol}(),

            # model indices
            # endogenous states unnormalized, endogenous states normalized
            OrderedDict{Symbol,UnitRange}(), OrderedDict{Symbol,UnitRange}(),
            OrderedDict{Symbol,Int}(), OrderedDict{Symbol,Int}(),
            OrderedDict{Symbol,UnitRange}(), OrderedDict{Symbol,Int}(), OrderedDict{Symbol, Int}(),

            spec,
            subspec,
            settings,
            test_settings,
            rng,
            testing,
            OrderedDict{Symbol,Observable}())
    # Set settings
    model_settings!(m)

    # default_test_settings!(m)
    for custom_setting in values(custom_settings)
        m <= custom_setting
    end

    # Set observable transformations
    init_observable_mappings!(m)

    # Initialize parameters
    init_parameters!(m)

    # Initialize grids
    init_grids!(m)

    # Initialize model indices
    init_model_indices!(m)

    # # Initialize subspecs
    init_subspec!(m)

    # Solve for the steady state
    steadystate!(m)

    # So that the indices of m.endogenous_states reflect the normalization
    normalize_model_state_indices!(m)

    return m
end

"""
```
init_parameters!(m::RealBondMkup)
```

Initializes the model's parameters, as well as empty values for the steady-state
parameters (in preparation for `steadystate!(m)` being called to initialize
those).
"""
function init_parameters!(m::RealBondMkup)

    # Steady State Parameters
    m <= parameter(:R, 1.04, fixed = true,
                   description = "R: Steady-state gross real interest rate.", tex_label = "R")
    m <= parameter(:γ, 1.0, fixed = true,
                   description = "γ: CRRA Parameter.", tex_label = "\\gamma")
    m <= parameter(:ν, 1.0, fixed = true,
                   description = "Inverse Frisch elasticity of labor supply.", tex_label = "\\nu")
    m <= parameter(:abar, -0.5, fixed = true,
                   description = "Borrowing floor.", tex_label = "\\bar{a}")

    # AR for observables (Non-steadystate)
    m <= parameter(:ρ_z, 0.03081565350294113, (1e-3, 0.999), (1e-3, 0.999), SquareRoot(), BetaAlt(0.5, 0.2),
                   fixed=false,
                   description="ρ_z: AR(1) coefficient in the technology process.",
                   tex_label="\\rho_z")
    m <= parameter(:ρ_mon, 0.5, (1e-3, 0.999), (1e-3, 0.999), SquareRoot(), BetaAlt(0.5, 0.2),
                   fixed = false,
                   description = "ρ_mon: Persistence of monetary policy shock",
                   tex_label = "\\rho_{mon}")
    m <= parameter(:ρ_mkp, 0.5, (1e-3, 0.999), (1e-3, 0.999), SquareRoot(), BetaAlt(0.5, 0.2),
                   fixed = false,
                   tex_label = "\\rho_{mkp}", description = "ρ_mkp: Persistence of the markup shock")

    m <= parameter(:σ_z, sqrt(.007), (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10),
                   fixed=false,
                   description="σ_z: The standard deviation of the process describing the stationary component of productivity.",
                   tex_label="\\sigma_{z}")
    m <= parameter(:σ_mon, 0.2380, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10),
                   fixed = false,
                   description="σ_mon: The standard deviation of the monetary policy shock.",
                   tex_label="\\sigma_{mon}")
    m <= parameter(:σ_mkp, 0.1314, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10),
                   fixed=false,
                   description="σ_mkp: The mean of the process that generates the price elasticity of the composite good. Specifically, the elasticity is (1+λ_{f,t})/(λ_{f_t}).",
                   tex_label="\\sigma_{\\lambda_f}")

    # More non-steadystate parameters
    m <= parameter(:ρ_tay, 0.5, (1e-5, 0.999), (1e-5, 0.999), SquareRoot(), BetaAlt(0.5, 0.2),
                   fixed = false,
                   description = "ρ_tay: Persistence in the taylor rule",
                   tex_label = "rho_{tay}")
    m <= parameter(:κ, 0.5, (1e-5, 0.999), (1e-5, 0.999), SquareRoot(), Distributions.Uniform(0.0, 1.0),
                   fixed = false,
                   description = "κ: The slope of the Phillips curve",
                   tex_label = "\\kappa")
    m <= parameter(:phipi, 1.5, (1e-5, 10.), (1e-5, 10.00), ModelConstructors.Exponential(), Normal(1.5, 0.25),
                   fixed = false,
                   description = "phipi: The slope of the taylor rule",
                   tex_label = "\\phi_\\pi")

    # More steadystate parameters
    m <= parameter(:μ_s, 0., fixed = true, description = "μ_s: Mu of log normal in income")
    m <= parameter(:σ_s, 0.1, fixed = true,
                   description = "σ_s: Sigma of log normal in income")
    m <= parameter(:e_y, 0., fixed = true, description = "e_y: Measurement error on GDP",
                   tex_label = "e_y")

    # Setting steady-state parameters
    nx = get_setting(m, :nx)
    ns = get_setting(m, :ns)
    nxns = nx*ns

    m <= SteadyStateParameterGrid(:lstar, fill(NaN, nxns), description = "Steady-state expected discounted
                                  marginal utility of consumption", tex_label = "l_*")
    m <= SteadyStateParameterGrid(:cstar, fill(NaN, nxns), description = "Steady-state consumption",
                                  tex_label = "c_*")
    m <= SteadyStateParameterGrid(:ηstar, fill(NaN, nxns), description = "Steady-state
                                  level of labor supply",
                                  tex_label = "\\eta_*")
    m <= SteadyStateParameterGrid(:μstar, fill(NaN, nxns), description = "Steady-state cross-sectional
                                  density of cash on hand", tex_label = "\\mu_*")

    # Figure out a better description for this...
    m <= SteadyStateParameterGrid(:χstar, fill(NaN, nxns), description = "Steady-state
                                  solution for constrained consumption and labor supply",
                                  tex_label = "\\chi_*")

    m <= SteadyStateParameter(:βstar, NaN, description = "Steady-state discount factor",
                              tex_label = "\\beta_*")
end

"""
```
init_grids!(m::RealBondMkup)
```
"""
function init_grids!(m::RealBondMkup)
    xscale  = get_setting(m, :xscale)
    xlo     = get_setting(m, :xlo)
    xhi     = get_setting(m, :xhi)
    nx      = get_setting(m, :nx)

    ns      = get_setting(m, :ns)
    λ       = get_setting(m, :λ)

    ehi     = get_setting(m, :ehi)

    grids = OrderedDict()

    # Cash on hand grid
    grids[:xgrid] = Grid(uniform_quadrature(xscale), xlo, xhi, nx, scale = xscale)

    # Skill grid
    lsgrid, sprob, sscale = tauchen86(m[:μ_s].value, m[:σ_s].value, ns, λ)
    swts = (sscale/ns)*ones(ns)
    sgrid = exp.(lsgrid) .+ ehi
    grids[:sgrid] = Grid(sgrid, swts, sscale)

    # Density of skill across skill grid
    grids[:ggrid] = sprob./swts

    # Total grid vectorized across both dimensions
    grids[:sgrid_total] = kron(sgrid, ones(nx))
    grids[:xgrid_total] = kron(ones(ns), grids[:xgrid].points)
    grids[:weights_total] = kron(swts, grids[:xgrid].weights)

    m.grids = grids
end

function model_settings!(m::RealBondMkup)
    default_settings!(m)

    # Defaults
    # Data settings for released and conditional data. Default behavior is to set vintage
    # of data to today's date.
    vint = Dates.format(now(), DSGE_DATE_FORMAT)
    m <= Setting(:data_vintage, vint, true, "vint", "Data vintage")

    saveroot = normpath(joinpath(dirname(@__FILE__), "../../../../","save"))
    datapath = normpath(joinpath(dirname(@__FILE__), "../../../../","save","input_data"))

    m <= Setting(:saveroot, saveroot, "Root of data directory structure")
    m <= Setting(:dataroot, datapath, "Input data directory path")

    # Solution method
    m <= Setting(:solution_method, :klein)

    # Likelihood method
    m <= Setting(:use_chand_recursion, true)

    # Anticipated shocks
    m <= Setting(:n_anticipated_shocks, 0,
                 "Number of anticipated policy shocks")
    # m <= Setting(:n_anticipated_shocks_padding, 20,
                 # "Padding for anticipated policy shocks")

    # Whether one wishes to re-compute βstar
    m <= Setting(:use_last_βstar, false, "Flag to avoid recomputing steady-state.")

    # Number of states and jumps
    m <= Setting(:normalize_distr_variables, true, "Whether or not to perform the
                 normalization of the distributional states in the Klein solution step")
    m <= Setting(:n_predetermined_variables, 0, "Number of predetermined variables after
                 removing the densities. Calculated with the Jacobian normalization.
                 This set to 0 is just the default setting, since it will always be
                 overwritten once the Jacobian is calculated.")

    m <= Setting(:state_indices, 1:5, "Which indices of m.endogenous_states correspond to
                 backward looking state variables")
    m <= Setting(:jump_indices, 6:11, "Which indices of m.endogenous_states correspond to jump
                 variables")

    # Mollifier setting parameters
    m <= Setting(:In, 0.443993816237631, "Normalizing constant for the mollifier")
    m <= Setting(:elo, 0.0, "Lower bound on stochastic consumption commitments")
    m <= Setting(:ehi, 1.0, "Upper bound on stochastic consumption commitments")

    # x: Cash on Hand Grid Setup
    m <= Setting(:nx, 150, "Cash on hand distribution grid points")
    m <= Setting(:xlo, -0.5 - 1.0, "Lower bound on cash on hand")
    m <= Setting(:xhi, 4.0, "Upper Bound on cash on hand")
    m <= Setting(:xscale, get_setting(m, :xhi) - get_setting(m, :xlo), "Size of the xgrid")

    # s: Skill Distribution/ "Units of effective labor" Grid Setup
    m <= Setting(:ns, 2, "Skill distribution grid points")
    m <= Setting(:λ, 3.0, "The λ parameter in the Tauchen distribution calculation")

    # Total grid x*s
    m <= Setting(:n, get_setting(m, :nx) * get_setting(m, :ns), "Total grid size, multiplying
                 across grid dimensions.")

    # Function-valued variables include distributional variables
    m <= Setting(:n_function_valued_backward_looking_states, 1, "Number of function-valued
                 backward looking state variables")
    m <= Setting(:n_backward_looking_distributional_vars, 1, "Number of state variables that are
                 distributional variables.")
    m <= Setting(:n_function_valued_jumps, 1, "Number of function-valued jump variables")
    m <= Setting(:n_jump_distributional_vars, 0, "Number of jump variables that are distributional
                 variables.")

    # Note, these settings assume normalization.
    # The n degrees of freedom removed depends on the distributions/dimensions
    # of heterogeneity that we have discretized over, in this case,
    # cash on hand and the skill distribution. In general the rule of
    # thumb is, remove one degree of freedom for the first endogenous distribution (cash on
    # hand), then one additional degree of freedom for each exogenous distribution (skill
    # distribution). Multiple endogenous distributions only permit removing a single degree
    # of freedom since it is then non-trivial to obtain the marginal distributions.
    m <= Setting(:n_degrees_of_freedom_removed, 2, "Number of degrees of freedom from the
                 distributional variables to remove.")
    n_dof_removed = get_setting(m, :n_degrees_of_freedom_removed)

    ####################################################
    # Calculating the number of backward-looking states
    ####################################################
    n_backward_looking_distr_vars = get_setting(m, :n_backward_looking_distributional_vars)
    m <= Setting(:backward_looking_states_normalization_factor,
                 n_dof_removed*n_backward_looking_distr_vars, "The number of dimensions removed from the
                 backward looking state variables for the normalization.")

    n = get_setting(m, :n)
    n_backward_looking_vars = length(get_setting(m, :state_indices))
    n_backward_looking_function_valued_vars = get_setting(m, :n_function_valued_backward_looking_states)
    n_backward_looking_scalar_vars = n_backward_looking_vars - n_backward_looking_function_valued_vars

    m <= Setting(:n_backward_looking_states, n*n_backward_looking_distr_vars +
                 n_backward_looking_scalar_vars - get_setting(m, :backward_looking_states_normalization_factor),
                 "Number of state variables, in the true sense (fully
                  backward looking) accounting for the discretization across the grid
                  of function-valued variables and the normalization of
                  distributional variables.")

    ##################################
    # Calculating the number of jumps
    ##################################
    n_jump_distr_vars = get_setting(m, :n_jump_distributional_vars)
    m <= Setting(:jumps_normalization_factor,
                 n_dof_removed*n_jump_distr_vars, "The number of dimensions removed from the
                 jump variables for the normalization.")

    n_jump_vars = length(get_setting(m, :jump_indices))
    n_jump_function_valued_vars = get_setting(m, :n_function_valued_jumps)
    n_jump_scalar_vars = n_jump_vars - n_jump_function_valued_vars

    m <= Setting(:n_jumps, n*n_jump_function_valued_vars +
                 n_jump_scalar_vars - get_setting(m, :jumps_normalization_factor),
                 "Number of jump variables (forward looking) accounting for
                  the discretization across the grid of function-valued variables and the
                  normalization of distributional variables.")

    m <= Setting(:n_model_states, get_setting(m, :n_backward_looking_states) + get_setting(m, :n_jumps),
                 "Number of 'states' in the state space model. Because backward and forward
                 looking variables need to be explicitly tracked for the Klein solution
                 method, we have n_states and n_jumps")

    # Reduction
    m <= Setting(:reduce_state_vars, true, "Reduce state variables or not")
    m <= Setting(:reduce_v, true, "Reduce value function or not")
    m <= Setting(:krylov_dim, 5, "Krylov reduction dimension")
    m <= Setting(:F, identity, "Function applied during block arnoldi deflation")
    m <= Setting(:n_knots, 12, "Number of knot points")
    m <= Setting(:c_power, 7, "Amount of bending of knot point locations to make them nonuniform")

end


text:
# Autogenerated wrapper script for ConnectFourSolver_jll for x86_64-apple-darwin
export libconnect4jl

using libcxxwrap_julia_jll
JLLWrappers.@generate_wrapper_header("ConnectFourSolver")
JLLWrappers.@declare_library_product(libconnect4jl, "@rpath/libconnect4jl.dylib")
function __init__()
    JLLWrappers.@generate_init_header(libcxxwrap_julia_jll)
    JLLWrappers.@init_library_product(
        libconnect4jl,
        "lib/libconnect4jl.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
import StatsBase: countmap

function solve_a(data)
    count = 0

    for line in data
        if length(Set(line)) == length(line)
            count += 1
        end
    end
    count
end

function solve_b(data)
    count = 0

    for line in data
        if length(Set(countmap.(line))) == length(line)
            count += 1
        end
    end
    count
end

data = split.(readlines("input.txt"))

sol_a = solve_a(data)
println("sol a: $(sol_a)")

sol_b = solve_b(data)
println("sol b: $(sol_b)")


text:
#!/usr/bin/env julia

using Distributed
using Random
using LinearAlgebra
procs_id = addprocs(8)
using DatagenCopulaBased
@everywhere using Distributions
@everywhere using Cumulants
using SymmetricTensors
using CumulantsFeatures
using JLD2
using FileIO
using ArgParse
@everywhere import CumulantsFeatures: reduceband
@everywhere using DatagenCopulaBased
@everywhere using CumulantsFeatures

function main(args)
  s = ArgParseSettings("description")
  @add_arg_table s begin
    "--nu", "-u"
    default = 5
    help = "the number of degrees of freedom for the t-Student copula"
    arg_type = Int
  end
  parsed_args = parse_args(s)
  ν = parsed_args["nu"]

  println(ν)
  @everywhere t = 100_000
  @everywhere n = 50
  @everywhere malf_size = 10
  @everywhere a = 1_000
  data_dir = "."
  test_number = 3
  filename = "tstudent_$(ν)-t_size-$(n)_malfsize-$malf_size-t_$(t)_$a.jld2"

  data = Dict{String, Any}("variables_no" => n,
                         "sample_number" => t,
                         "ν" => ν,
                         "test_number" => test_number,
                         "malf_size" => malf_size,
                         "a" => a,
                         "data" => Dict{String, Dict{String,Any}}())

  known_data_size = 0
  if isfile("$data_dir/$filename")
   data["data"] = load("$data_dir/$filename")["data"]
   known_data_size += length(data["data"])
   println("Already have $known_data_size samples \n Will generate $(test_number-known_data_size) more")
  end

  #true calculations
  println("Calculation started")
  for m=(known_data_size+1):test_number
    @time begin
      println(" > $m ($ν)")
      malf = randperm(n)[1:malf_size]
      Σ = cormatgen_rand(n)
      samples_orig = rand(MvNormal(Σ), t)'

      versions = [(x->x, "original"),
                  (x->vcat(gcop2tstudent(x[1:a, :], malf, ν), x[a+1:end, :]), "malf")]

      cur_dict = Dict{String, Any}("malf" => malf,
                                   "cor_source" => Σ)

      data_dict = @distributed (merge) for (sampler, label)=versions
        println(label)
        samples = sampler(samples_orig)
        Σ_malf = cov(samples)
        Dict("cor_$label" => Σ_malf,
             "x_$label" => samples)
      end

      data["data"]["$m"] = merge(cur_dict, data_dict)
      save("$data_dir/$filename", data)
    end
  end

end

main(ARGS)


text:
# Copyright 2018 Martin Holters
# See accompanying license file.
ACMEdir = joinpath(@__DIR__, "..")
println("Checking copyright headers...")
for dirname in ("src", "examples", "test")
    dirname = joinpath(ACMEdir, dirname)
    for name in readdir(dirname)
        if endswith(name, ".jl")
            name = joinpath(dirname, name)
            local years
            try
                years = sort!(unique(parse.([Int], readlines(`git log --format=%cd --date=format:%Y -- $name`))))
            catch e
                @static if isdefined(Base, Symbol("@warn"))
                    @warn e
                else
                    warn(e)
                end
                continue
            end
            if isempty(years)
                continue
            end
            println(name)
            open(name, "r") do io
                l = readline(io)
                m = match(r"#\s*Copyright\s+(([0-9]+(,\s*)?)*)", l)
                if m === nothing
                    error("Missing copyright header in $name")
                end
                headyears = parse.([Int], strip.(split(m.captures[1], ',')))
                d = setdiff(years, headyears)
                if !isempty(d)
                    error("Missing years in copyright header of $name: $d")
                end
            end
        end
    end
end


text:
function hline(varargs...; kwargs...)
  simple_plot = SimplePlot(; kwargs...)
  hline!(simple_plot, varargs...; kwargs...)
end

function vline(varargs...; kwargs...)
  simple_plot = SimplePlot(; kwargs...)
  vline!(simple_plot, varargs...; kwargs...)
end

function hline!(simple_plot::SimplePlot, varargs...; kwargs...)
  _line!(simple_plot, "x", varargs...; kwargs...)
end

function vline!(simple_plot::SimplePlot, varargs...; kwargs...)
  _line!(simple_plot, "y", varargs...; kwargs...)
end

hline!(varargs...; kwargs...) = hline!(_plot, varargs...; kwargs...)
vline!(varargs...; kwargs...) = vline!(_plot, varargs...; kwargs...)

export hline, hline!
export vline, vline!


text:
using DBAPI

include("./failedinterface.jl")
include("./arrays.jl")


text:
import UUIDs

"""
    latest(name::AbstractString;
           cache = RegistryParsingCache(),
           registries = [],
           include_default_registries::Bool = true)

Get the latest registered version of a package, given the name of the package.

## Examples

```julia
julia> Pkg.add(latest("Example"))
```
"""
@inline function latest(name::AbstractString;
                        cache = RegistryParsingCache(),
                        registries = [],
                        include_default_registries::Bool = true)
    package = PartialUnversionedPackage(name)
    result = _latest(
        package;
        cache = cache,
        registries = registries,
        include_default_registries = include_default_registries,
    )
    return result
end

"""
    latest(name::AbstractString,
           uuid::UUID;
           cache = RegistryParsingCache(),
           registries = [],
           include_default_registries::Bool = true)

Get the latest registered version of a package, given the name and UUID of the package.

## Examples

```julia
julia> Pkg.add(latest("Example", UUID("7876af07-990d-54b4-ab0e-23690620f79a")))
```
"""
@inline function latest(name::AbstractString,
                        uuid::UUIDs.UUID;
                        cache = RegistryParsingCache(),
                        registries = [],
                        include_default_registries::Bool = true)
    package = PartialUnversionedPackage(name, uuid)
    result = _latest(
        package;
        cache = cache,
        registries = registries,
        include_default_registries = include_default_registries,
    )
    return result
end

@inline function _latest(package::PartialUnversionedPackage,;
                         cache,
                         registries,
                         include_default_registries::Bool)
    registry_list = make_registry_list(registries;
                                       include_default_registries = include_default_registries)
    parse_registries!(cache, registry_list)
    full_package = identify_package(package, cache)
    package.name == full_package.name || throw(PackageNotFoundException("An unexpected error occured when trying to find the package"))
    if package.uuid !== nothing
        package.uuid == full_package.uuid || throw(PackageNotFoundException("An unexpected error occured when trying to find the package"))
    end
    versions = Vector{VersionAndTree}(undef, 0)
    for registry in keys(cache)
        if haskey(cache[registry], full_package)
            push!(versions, cache[registry][full_package])
        end
    end
    num_versions = length(versions)
    version_numbers = [VersionNumber(versions[j].version) for j in 1:num_versions]
    max_version_number = maximum(version_numbers)
    indices_of_max_version_numbers = findall(version_numbers .== max_version_number)
    versions_at_max = versions[indices_of_max_version_numbers]
    first_version_at_max = versions_at_max[1]
    for v in versions_at_max
        assert_no_disagreement(first_version_at_max, v)
    end
    result = VersionedPackage(
        full_package.name,
        full_package.uuid,
        first_version_at_max.version,
        first_version_at_max.tree,
    )
    return result
end


text:
module BipartiteGraphs

export BipartiteEdge, BipartiteGraph, DiCMOBiGraph, Unassigned, unassigned,
        Matching, ResidualCMOGraph, InducedCondensationGraph, maximal_matching,
        construct_augmenting_path!, MatchedCondensationGraph

export 𝑠vertices, 𝑑vertices, has_𝑠vertex, has_𝑑vertex, 𝑠neighbors, 𝑑neighbors,
       𝑠edges, 𝑑edges, nsrcs, ndsts, SRC, DST, set_neighbors!, invview,
       complete

using DocStringExtensions
using UnPack
using SparseArrays
using Graphs
using Setfield

### Matching
struct Unassigned
    global unassigned
    const unassigned = Unassigned.instance
end
# Behaves as a scalar
Base.length(u::Unassigned) = 1
Base.size(u::Unassigned) = ()
Base.iterate(u::Unassigned) = (unassigned, nothing)
Base.iterate(u::Unassigned, state) = nothing

Base.show(io::IO, ::Unassigned) =
    printstyled(io, "u"; color=:light_black)

struct Matching{U #=> :Unassigned =#, V<:AbstractVector} <: AbstractVector{Union{U, Int}}
    match::V
    inv_match::Union{Nothing, V}
end
# These constructors work around https://github.com/JuliaLang/julia/issues/41948
function Matching{V}(m::Matching) where {V}
    eltype(m) === Union{V, Int} && return M
    VUT = typeof(similar(m.match, Union{V, Int}))
    Matching{V}(convert(VUT, m.match),
        m.inv_match === nothing ? nothing : convert(VUT, m.inv_match))
end
Matching(m::Matching) = m
Matching{U}(v::V) where {U, V<:AbstractVector} = Matching{U, V}(v, nothing)
Matching{U}(v::V, iv::Union{V, Nothing}) where {U, V<:AbstractVector} = Matching{U, V}(v, iv)
Matching(v::V) where {U, V<:AbstractVector{Union{U, Int}}} =
    Matching{@isdefined(U) ? U : Unassigned, V}(v, nothing)
Matching(m::Int) = Matching{Unassigned}(Union{Int, Unassigned}[unassigned for _ = 1:m], nothing)

Base.size(m::Matching) = Base.size(m.match)
Base.getindex(m::Matching, i::Integer) = m.match[i]
Base.iterate(m::Matching, state...) = iterate(m.match, state...)
Base.copy(m::Matching) = Matching(copy(m.match), m.inv_match === nothing ? nothing : copy(m.inv_match))
function Base.setindex!(m::Matching{U}, v::Union{Integer, U}, i::Integer) where {U}
    if m.inv_match !== nothing
        oldv = m.match[i]
        isa(oldv, Int) && (m.inv_match[oldv] = unassigned)
        isa(v, Int) && (m.inv_match[v] = i)
    end
    return m.match[i] = v
end

function Base.push!(m::Matching{U}, v::Union{Integer, U}) where {U}
    push!(m.match, v)
    if v !== unassigned && m.inv_match !== nothing
        m.inv_match[v] = length(m.match)
    end
end

function complete(m::Matching{U}, N=length(m.match)) where {U}
    m.inv_match !== nothing && return m
    inv_match = Union{U, Int}[unassigned for _ = 1:N]
    for (i, eq) in enumerate(m.match)
        isa(eq, Int) || continue
        inv_match[eq] = i
    end
    return Matching{U}(collect(m.match), inv_match)
end

@noinline require_complete(m::Matching) =
    m.inv_match === nothing && throw(ArgumentError("Backwards matching not defined. `complete` the matching first."))

function invview(m::Matching{U, V}) where {U, V}
    require_complete(m)
    return Matching{U, V}(m.inv_match, m.match)
end

###
### Edges & Vertex
###
@enum VertType SRC DST

struct BipartiteEdge{I<:Integer} <: Graphs.AbstractEdge{I}
    src::I
    dst::I
    function BipartiteEdge(src::I, dst::V) where {I,V}
        T = promote_type(I, V)
        new{T}(T(src), T(dst))
    end
end

Graphs.src(edge::BipartiteEdge) = edge.src
Graphs.dst(edge::BipartiteEdge) = edge.dst

function Base.show(io::IO, edge::BipartiteEdge)
    @unpack src, dst = edge
    print(io, "[src: ", src, "] => [dst: ", dst, "]")
end

Base.:(==)(a::BipartiteEdge, b::BipartiteEdge) = src(a) == src(b) && dst(a) == dst(b)

###
### Graph
###
"""
$(TYPEDEF)

A bipartite graph representation between two, possibly distinct, sets of vertices
(source and dependencies). Maps source vertices, labelled `1:N₁`, to vertices
on which they depend (labelled `1:N₂`).

# Fields
$(FIELDS)

# Example
```julia
using ModelingToolkit

ne = 4
srcverts = 1:4
depverts = 1:2

# six source vertices
fadjlist = [[1],[1],[2],[2],[1],[1,2]]

# two vertices they depend on
badjlist = [[1,2,5,6],[3,4,6]]

bg = BipartiteGraph(7, fadjlist, badjlist)
```
"""
mutable struct BipartiteGraph{I<:Integer, M} <: Graphs.AbstractGraph{I}
    ne::Int
    fadjlist::Vector{Vector{I}} # `fadjlist[src] => dsts`
    badjlist::Union{Vector{Vector{I}},I} # `badjlist[dst] => srcs` or `ndsts`
    metadata::M
end
BipartiteGraph(ne::Integer, fadj::AbstractVector, badj::Union{AbstractVector,Integer}=maximum(maximum, fadj); metadata=nothing) = BipartiteGraph(ne, fadj, badj, metadata)
BipartiteGraph(fadj::AbstractVector, badj::Union{AbstractVector,Integer}=maximum(maximum, fadj); metadata=nothing) =
    BipartiteGraph(mapreduce(length, +, fadj; init=0), fadj, badj, metadata)

@noinline require_complete(g::BipartiteGraph) = g.badjlist isa AbstractVector || throw(ArgumentError("The graph has no back edges. Use `complete`."))

function invview(g::BipartiteGraph)
    require_complete(g)
    BipartiteGraph(g.ne, g.badjlist, g.fadjlist)
end

function complete(g::BipartiteGraph{I}) where {I}
    isa(g.badjlist, AbstractVector) && return g
    badjlist = Vector{I}[Vector{I}() for _ in 1:g.badjlist]
    for (s, l) in enumerate(g.fadjlist)
        for d in l
            push!(badjlist[d], s)
        end
    end
    BipartiteGraph(g.ne, g.fadjlist, badjlist)
end

"""
```julia
Base.isequal(bg1::BipartiteGraph{T}, bg2::BipartiteGraph{T}) where {T<:Integer}
```

Test whether two [`BipartiteGraph`](@ref)s are equal.
"""
function Base.isequal(bg1::BipartiteGraph{T}, bg2::BipartiteGraph{T}) where {T<:Integer}
    iseq = (bg1.ne == bg2.ne)
    iseq &= (bg1.fadjlist == bg2.fadjlist)
    iseq &= (bg1.badjlist == bg2.badjlist)
    iseq
end

"""
$(SIGNATURES)

Build an empty `BipartiteGraph` with `nsrcs` sources and `ndsts` destinations.
"""
function BipartiteGraph(nsrcs::T, ndsts::T, backedge::Val{B}=Val(true); metadata=nothing) where {T,B}
    fadjlist = map(_->T[], 1:nsrcs)
    badjlist = B ? map(_->T[], 1:ndsts) : ndsts
    BipartiteGraph(0, fadjlist, badjlist, metadata)
end

Base.copy(bg::BipartiteGraph) = BipartiteGraph(bg.ne, copy(bg.fadjlist), copy(bg.badjlist), deepcopy(bg.metadata))
Base.eltype(::Type{<:BipartiteGraph{I}}) where I = I
function Base.empty!(g::BipartiteGraph)
    foreach(empty!, g.fadjlist)
    g.badjlist isa AbstractVector && foreach(empty!, g.badjlist)
    g.ne = 0
    if g.metadata !== nothing
        foreach(empty!, g.metadata)
    end
    g
end
Base.length(::BipartiteGraph) = error("length is not well defined! Use `ne` or `nv`.")

if isdefined(Graphs, :has_contiguous_vertices)
    Graphs.has_contiguous_vertices(::Type{<:BipartiteGraph}) = false
end
Graphs.is_directed(::Type{<:BipartiteGraph}) = false
Graphs.vertices(g::BipartiteGraph) = (𝑠vertices(g), 𝑑vertices(g))
𝑠vertices(g::BipartiteGraph) = axes(g.fadjlist, 1)
𝑑vertices(g::BipartiteGraph) = g.badjlist isa AbstractVector ? axes(g.badjlist, 1) : Base.OneTo(g.badjlist)
has_𝑠vertex(g::BipartiteGraph, v::Integer) = v in 𝑠vertices(g)
has_𝑑vertex(g::BipartiteGraph, v::Integer) = v in 𝑑vertices(g)
𝑠neighbors(g::BipartiteGraph, i::Integer, with_metadata::Val{M}=Val(false)) where M = M ? zip(g.fadjlist[i], g.metadata[i]) : g.fadjlist[i]
function 𝑑neighbors(g::BipartiteGraph, j::Integer, with_metadata::Val{M}=Val(false)) where M
    require_complete(g)
    M ? zip(g.badjlist[j], (g.metadata[i][j] for i in g.badjlist[j])) : g.badjlist[j]
end
Graphs.ne(g::BipartiteGraph) = g.ne
Graphs.nv(g::BipartiteGraph) = sum(length, vertices(g))
Graphs.edgetype(g::BipartiteGraph{I}) where I = BipartiteEdge{I}

nsrcs(g::BipartiteGraph) = length(𝑠vertices(g))
ndsts(g::BipartiteGraph) = length(𝑑vertices(g))

function Graphs.has_edge(g::BipartiteGraph, edge::BipartiteEdge)
    @unpack src, dst = edge
    (src in 𝑠vertices(g) && dst in 𝑑vertices(g)) || return false  # edge out of bounds
    insorted(dst, 𝑠neighbors(g, src))
end
Base.in(edge::BipartiteEdge, g::BipartiteGraph) = Graphs.has_edge(g, edge)

### Maximal matching
"""
    construct_augmenting_path!(m::Matching, g::BipartiteGraph, vsrc, dstfilter, vcolor=falses(ndsts(g)), ecolor=falses(nsrcs(g))) -> path_found::Bool

Try to construct an augmenting path in matching and if such a path is found,
update the matching accordingly.
"""
function construct_augmenting_path!(matching::Matching, g::BipartiteGraph, vsrc, dstfilter, dcolor=falses(ndsts(g)), scolor=falses(nsrcs(g)))
    scolor[vsrc] = true

    # if a `vdst` is unassigned and the edge `vsrc <=> vdst` exists
    for vdst in 𝑠neighbors(g, vsrc)
        if dstfilter(vdst) && matching[vdst] === unassigned
            matching[vdst] = vsrc
            return true
        end
    end

    # for every `vsrc` such that edge `vsrc <=> vdst` exists and `vdst` is uncolored
    for vdst in 𝑠neighbors(g, vsrc)
        (dstfilter(vdst) && !dcolor[vdst]) || continue
        dcolor[vdst] = true
        if construct_augmenting_path!(matching, g, matching[vdst], dstfilter, dcolor, scolor)
            matching[vdst] = vsrc
            return true
        end
    end
    return false
end

"""
    maximal_matching(g::BipartiteGraph, [srcfilter], [dstfilter])

For a bipartite graph `g`, construct a maximal matching of destination to source
vertices, subject to the constraint that vertices for which `srcfilter` or `dstfilter`,
return `false` may not be matched.
"""
function maximal_matching(g::BipartiteGraph, srcfilter=vsrc->true, dstfilter=vdst->true)
    matching = Matching(ndsts(g))
    foreach(Iterators.filter(srcfilter, 𝑠vertices(g))) do vsrc
        construct_augmenting_path!(matching, g, vsrc, dstfilter)
    end
    return matching
end

###
### Populate
###
struct NoMetadata
end
const NO_METADATA = NoMetadata()

Graphs.add_edge!(g::BipartiteGraph, i::Integer, j::Integer, md=NO_METADATA) = add_edge!(g, BipartiteEdge(i, j), md)
function Graphs.add_edge!(g::BipartiteGraph, edge::BipartiteEdge, md=NO_METADATA)
    @unpack fadjlist, badjlist = g
    s, d = src(edge), dst(edge)
    (has_𝑠vertex(g, s) && has_𝑑vertex(g, d)) || error("edge ($edge) out of range.")
    @inbounds list = fadjlist[s]
    index = searchsortedfirst(list, d)
    @inbounds (index <= length(list) && list[index] == d) && return false  # edge already in graph
    insert!(list, index, d)
    if md !== NO_METADATA
        insert!(g.metadata[s], index, md)
    end

    g.ne += 1
    if badjlist isa AbstractVector
        @inbounds list = badjlist[d]
        index = searchsortedfirst(list, s)
        insert!(list, index, s)
    end
    return true  # edge successfully added
end

Graphs.rem_edge!(g::BipartiteGraph, i::Integer, j::Integer) =
    Graphs.rem_edge!(g, BipartiteEdge(i, j))
function Graphs.rem_edge!(g::BipartiteGraph, edge::BipartiteEdge)
    @unpack fadjlist, badjlist = g
    s, d = src(edge), dst(edge)
    (has_𝑠vertex(g, s) && has_𝑑vertex(g, d)) || error("edge ($edge) out of range.")
    @inbounds list = fadjlist[s]
    index = searchsortedfirst(list, d)
    @inbounds (index <= length(list) && list[index] == d) || error("graph does not have edge $edge")
    deleteat!(list, index)
    g.ne -= 1
    if badjlist isa AbstractVector
        @inbounds list = badjlist[d]
        index = searchsortedfirst(list, s)
        deleteat!(list, index)
    end
    return true  # edge successfully deleted
end

function Graphs.add_vertex!(g::BipartiteGraph{T}, type::VertType) where T
    if type === DST
        if g.badjlist isa AbstractVector
            push!(g.badjlist, T[])
        else
            g.badjlist += 1
        end
    elseif type === SRC
        push!(g.fadjlist, T[])
    else
        error("type ($type) must be either `DST` or `SRC`")
    end
    return true  # vertex successfully added
end

function set_neighbors!(g::BipartiteGraph, i::Integer, new_neighbors::AbstractVector)
    old_neighbors = g.fadjlist[i]
    old_nneighbors = length(old_neighbors)
    new_nneighbors = length(new_neighbors)
    g.fadjlist[i] = new_neighbors
    g.ne += new_nneighbors - old_nneighbors
    if isa(g.badjlist, AbstractVector)
        for n in old_neighbors
            @inbounds list = g.badjlist[n]
            index = searchsortedfirst(list, i)
            deleteat!(list, index)
        end
        for n in new_neighbors
            @inbounds list = g.badjlist[n]
            index = searchsortedfirst(list, i)
            insert!(list, index, i)
        end
    end
end

###
### Edges iteration
###
Graphs.edges(g::BipartiteGraph) = BipartiteEdgeIter(g, Val(SRC))
𝑠edges(g::BipartiteGraph) = BipartiteEdgeIter(g, Val(SRC))
𝑑edges(g::BipartiteGraph) = BipartiteEdgeIter(g, Val(DST))

struct BipartiteEdgeIter{T,G} <: Graphs.AbstractEdgeIter
    g::G
    type::Val{T}
end

Base.length(it::BipartiteEdgeIter) = ne(it.g)
Base.eltype(it::BipartiteEdgeIter) = edgetype(it.g)

function Base.iterate(it::BipartiteEdgeIter{SRC,<:BipartiteGraph{T}}, state=(1, 1, SRC)) where T
    @unpack g = it
    neqs = nsrcs(g)
    neqs == 0 && return nothing
    eq, jvar = state

    while eq <= neqs
        eq′ = eq
        vars = 𝑠neighbors(g, eq′)
        if jvar > length(vars)
            eq += 1
            jvar = 1
            continue
        end
        edge = BipartiteEdge(eq′, vars[jvar])
        state = (eq, jvar + 1, SRC)
        return edge, state
    end
    return nothing
end

function Base.iterate(it::BipartiteEdgeIter{DST,<:BipartiteGraph{T}}, state=(1, 1, DST)) where T
    @unpack g = it
    nvars = ndsts(g)
    nvars == 0 && return nothing
    ieq, jvar = state

    while jvar <= nvars
        eqs = 𝑑neighbors(g, jvar)
        if ieq > length(eqs)
            ieq = 1
            jvar += 1
            continue
        end
        edge = BipartiteEdge(eqs[ieq], jvar)
        state = (ieq + 1, jvar, DST)
        return edge, state
    end
    return nothing
end

###
### Utils
###
function Graphs.incidence_matrix(g::BipartiteGraph, val=true)
    I = Int[]
    J = Int[]
    for i in 𝑠vertices(g), n in 𝑠neighbors(g, i)
        push!(I, i)
        push!(J, n)
    end
    S = sparse(I, J, val, nsrcs(g), ndsts(g))
end

"""
    struct DiCMOBiGraph

This data structure implements a "directed, contracted, matching-oriented" view of an
original (undirected) bipartite graph. It has two modes, depending on the `Transposed`
flag, which switches the direction of the induced matching.

Essentially the graph adapter performs two largely orthogonal functions
[`Transposed == true` differences are indicated in square brackets]:

1. It pairs an undirected bipartite graph with a matching of the destination vertex.

    This matching is used to induce an orientation on the otherwise undirected graph:
    Matched edges pass from destination to source [source to desination], all other edges
    pass in the opposite direction.

2. It exposes the graph view obtained by contracting the destination [source] vertices
   along the matched edges.

The result of this operation is an induced, directed graph on the source [destination] vertices.
The resulting graph has a few desirable properties. In particular, this graph
is acyclic if and only if the induced directed graph on the original bipartite
graph is acyclic.

# Hypergraph interpretation

Consider the bipartite graph `B` as the incidence graph of some hypergraph `H`.
Note that a maching `M` on `B` in the above sense is equivalent to determining
an (1,n)-orientation on the hypergraph (i.e. each directed hyperedge has exactly
one head, but any arbitrary number of tails). In this setting, this is simply
the graph formed by expanding each directed hyperedge into `n` ordinary edges
between the same vertices.
"""
mutable struct DiCMOBiGraph{Transposed, I, G<:BipartiteGraph{I}, M <: Matching} <: Graphs.AbstractGraph{I}
    graph::G
    ne::Union{Missing, Int}
    matching::M
    DiCMOBiGraph{Transposed}(g::G, ne::Union{Missing, Int}, m::M) where {Transposed, I, G<:BipartiteGraph{I}, M} =
        new{Transposed, I, G, M}(g, ne, m)
end
function DiCMOBiGraph{Transposed}(g::BipartiteGraph) where {Transposed}
    DiCMOBiGraph{Transposed}(g, 0, Matching(ndsts(g)))
end
function DiCMOBiGraph{Transposed}(g::BipartiteGraph, m::M) where {Transposed, M}
    DiCMOBiGraph{Transposed}(g, missing, m)
end

invview(g::DiCMOBiGraph{Transposed}) where {Transposed} =
    DiCMOBiGraph{!Transposed}(invview(g.graph), g.ne, invview(g.matching))

Graphs.is_directed(::Type{<:DiCMOBiGraph}) = true
Graphs.nv(g::DiCMOBiGraph{Transposed}) where {Transposed} = Transposed ? ndsts(g.graph) : nsrcs(g.graph)
Graphs.vertices(g::DiCMOBiGraph{Transposed}) where {Transposed} = Transposed ? 𝑑vertices(g.graph) : 𝑠vertices(g.graph)

struct CMONeighbors{Transposed, V}
    g::DiCMOBiGraph{Transposed}
    v::V
    CMONeighbors{Transposed}(g::DiCMOBiGraph{Transposed}, v::V) where {Transposed, V} =
        new{Transposed, V}(g, v)
end

Graphs.outneighbors(g::DiCMOBiGraph{false}, v) = CMONeighbors{false}(g, v)
Graphs.inneighbors(g::DiCMOBiGraph{false}, v) = inneighbors(invview(g), v)
Base.iterate(c::CMONeighbors{false}) = iterate(c, (c.g.graph.fadjlist[c.v],))
function Base.iterate(c::CMONeighbors{false}, (l, state...))
    while true
        r = iterate(l, state...)
        r === nothing && return nothing
        # If this is a matched edge, skip it, it's reversed in the induced
        # directed graph. Otherwise, if there is no matching for this destination
        # edge, also skip it, since it got delted in the contraction.
        vsrc = c.g.matching[r[1]]
        if vsrc === c.v || vsrc === unassigned
            state = (r[2],)
            continue
        end
        return vsrc, (l, r[2])
    end
end
Base.length(c::CMONeighbors{false}) = count(_->true, c)

liftint(f, x) = (!isa(x, Int)) ? nothing : f(x)
liftnothing(f, x) = x === nothing ? nothing : f(x)

_vsrc(c::CMONeighbors{true}) = c.g.matching[c.v]
_neighbors(c::CMONeighbors{true}) = liftint(vsrc->c.g.graph.fadjlist[vsrc], _vsrc(c))
Base.length(c::CMONeighbors{true}) = something(liftnothing(length, _neighbors(c)), 1) - 1
Graphs.inneighbors(g::DiCMOBiGraph{true}, v) = CMONeighbors{true}(g, v)
Graphs.outneighbors(g::DiCMOBiGraph{true}, v) = outneighbors(invview(g), v)
Base.iterate(c::CMONeighbors{true}) = liftnothing(ns->iterate(c, (ns,)), _neighbors(c))
function Base.iterate(c::CMONeighbors{true}, (l, state...))
    while true
        r = iterate(l, state...)
        r === nothing && return nothing
        if r[1] === c.v
            state = (r[2],)
            continue
        end
        return r[1], (l, r[2])
    end
end

_edges(g::DiCMOBiGraph{Transposed}) where Transposed = Transposed ?
    ((w=>v for w in inneighbors(g, v)) for v in vertices(g)) :
    ((v=>w for w in outneighbors(g, v)) for v in vertices(g))

Graphs.edges(g::DiCMOBiGraph) = (Graphs.SimpleEdge(p) for p in Iterators.flatten(_edges(g)))
function Graphs.ne(g::DiCMOBiGraph)
    if g.ne === missing
        g.ne = mapreduce(x->length(x.iter), +, _edges(g))
    end
    return g.ne
end

Graphs.has_edge(g::DiCMOBiGraph{true}, a, b) = a in inneighbors(g, b)
Graphs.has_edge(g::DiCMOBiGraph{false}, a, b) = b in outneighbors(g, a)

# Condensation Graphs
abstract type AbstractCondensationGraph <: AbstractGraph{Int}; end
function (T::Type{<:AbstractCondensationGraph})(g, sccs::Vector{Union{Int, Vector{Int}}})
    scc_assignment = Vector{Int}(undef, isa(g, BipartiteGraph) ? ndsts(g) : nv(g))
    for (i, c) in enumerate(sccs)
        for v in c
            scc_assignment[v] = i
        end
    end
    T(g, sccs, scc_assignment)
end
(T::Type{<:AbstractCondensationGraph})(g, sccs::Vector{Vector{Int}}) =
    T(g, Vector{Union{Int, Vector{Int}}}(sccs))

Graphs.is_directed(::Type{<:AbstractCondensationGraph}) = true
Graphs.nv(icg::AbstractCondensationGraph) = length(icg.sccs)
Graphs.vertices(icg::AbstractCondensationGraph) = Base.OneTo(nv(icg))

"""
    struct MatchedCondensationGraph

For some bipartite-graph and an orientation induced on its destination contraction,
records the condensation DAG of the digraph formed by the orientation. I.e. this
is a DAG of connected components formed by the destination vertices of some
underlying bipartite graph.
N.B.: This graph does not store explicit neighbor relations of the sccs.
Therefor, the edge multiplicity is derived from the underlying bipartite graph,
i.e. this graph is not strict.
"""
struct MatchedCondensationGraph{G <: DiCMOBiGraph} <: AbstractCondensationGraph
    graph::G
    # Records the members of a strongly connected component. For efficiency,
    # trivial sccs (with one vertex member) are stored inline. Note: the sccs
    # here need not be stored in topological order.
    sccs::Vector{Union{Int, Vector{Int}}}
    # Maps the vertices back to the scc of which they are a part
    scc_assignment::Vector{Int}
end

Graphs.outneighbors(mcg::MatchedCondensationGraph, cc::Integer) =
    Iterators.flatten((mcg.scc_assignment[v′] for v′ in outneighbors(mcg.graph, v) if mcg.scc_assignment[v′] != cc) for v in mcg.sccs[cc])

Graphs.inneighbors(mcg::MatchedCondensationGraph, cc::Integer) =
    Iterators.flatten((mcg.scc_assignment[v′] for v′ in inneighbors(mcg.graph, v) if mcg.scc_assignment[v′] != cc) for v in mcg.sccs[cc])

"""
    struct InducedCondensationGraph

For some bipartite-graph and a topologicall sorted list of connected components,
represents the condensation DAG of the digraph formed by the orientation. I.e. this
is a DAG of connected components formed by the destination vertices of some
underlying bipartite graph.
N.B.: This graph does not store explicit neighbor relations of the sccs.
Therefor, the edge multiplicity is derived from the underlying bipartite graph,
i.e. this graph is not strict.
"""
struct InducedCondensationGraph{G <: BipartiteGraph} <: AbstractCondensationGraph
    graph::G
    # Records the members of a strongly connected component. For efficiency,
    # trivial sccs (with one vertex member) are stored inline. Note: the sccs
    # here are stored in topological order.
    sccs::Vector{Union{Int, Vector{Int}}}
    # Maps the vertices back to the scc of which they are a part
    scc_assignment::Vector{Int}
end

_neighbors(icg::InducedCondensationGraph, cc::Integer) =
    Iterators.flatten(Iterators.flatten(icg.graph.fadjlist[vsrc] for vsrc in icg.graph.badjlist[v]) for v in icg.sccs[cc])

Graphs.outneighbors(icg::InducedCondensationGraph, v::Integer) =
    (icg.scc_assignment[n] for n in _neighbors(icg, v) if icg.scc_assignment[n] > v)

Graphs.inneighbors(icg::InducedCondensationGraph, v::Integer) =
    (icg.scc_assignment[n] for n in _neighbors(icg, v) if icg.scc_assignment[n] < v)

end # module


text:
using Knet, ArgParse, HDF5, MAT, JSON
include("lib/convnet.jl")
include("lib/imgproc.jl")
include("lib/data.jl")
SPLITS = ["train", "restval", "val", "test"]

function main(args)
    s = ArgParseSettings()
    s.description = "Extract CNN features of images (now only just for VGG-19)"

    @add_arg_table s begin
        ("--images"; help="images file in JLD format")
        ("--captions"; help="captions archieve (karpathy)")
        ("--savefile"; help="save file in JLD format")
        ("--cnnfile"; help="CNN model file")
        ("--convnet"; default="vgg16")
        ("--lastlayer"; default="relu7"; help="layer for feature extraction")
        ("--seed"; arg_type=Int; default=0; help="random seed")
        ("--fc6drop"; arg_type=Float32; default=Float32(0.0))
        ("--feedback"; arg_type=Int; default=0; help="feedback in every N image")
        ("--extradata"; action=:store_true)
    end

    # parse args
    isa(args, AbstractString) && (args=split(args))
    o = parse_args(args, s; as_symbols=true); println(o); flush(STDOUT)
    o[:seed] > 0 && srand(o[:seed])

    # load CNN
    @printf("CNN loading...\n"); flush(STDOUT)
    lastlayer = o[:lastlayer]
    convnet = eval(parse(o[:convnet]))
    CNNmat = matread(o[:cnnfile])
    weights = get_vgg_weights(CNNmat; last_layer=lastlayer)
    avgimg = CNNmat["meta"]["normalization"]["averageImage"]
    avgimg = KnetArray(make_average_image(avgimg))

    # FIXME: this only works for vgg19 right now
    featuremaps = lastlayer == "conv5_4"
    if featuremaps
        @printf("featuremaps from %s layer.\n", lastlayer)
    end
    @printf("Done.\n"); flush(STDOUT)

    # feature extraction
    h5open(abspath(o[:savefile]), "w") do outfile
        entries = get_entries(o[:captions], SPLITS)
        for (k,splitname) in enumerate(SPLITS)
            @printf("Feature extraction for %s split...\n", splitname)
            flush(STDOUT)
            dropouts = Dict()
            if splitname == "train" || (splitname == "restval" && o[:extradata])
                dropouts = Dict(:fc6drop => o[:fc6drop])
            end
            for i = 1:length(entries[k])
                filename = entries[k][i]["filename"]
                image = h5open(o[:images], "r") do f
                    read(f, filename)
                end
                dict = Dict(:featuremaps => featuremaps)
                feats = convnet(weights, KnetArray(image) .- avgimg; o=dict)
                feats = convert(Array{Float32}, feats)
                if !featuremaps
                    feats = reshape(feats, length(feats), 1)
                end
                write(outfile, filename, feats)

                if o[:feedback] > 0 && i % o[:feedback] == 0
                    @printf("[%%%.1f] %d images processed so far [%s]\n",
                            round(100*i/length(entries[k]), 1), i, now())
                    flush(STDOUT)
                end
            end
        end
    end

    @printf("Save extracted features to %s\n", abspath(o[:savefile]))
end

!isinteractive() && !isdefined(Core.Main, :load_only) && main(ARGS)

max_stars_repo_name: ilkerkesen/kusrp

max_issues_repo_name: ilkerkesen/kusrp

max_forks_repo_name: ilkerkesen/kusrp


text:
# using Revise
using Test

# brittle hack b/s of https://github.com/dmlc/XGBoost.jl/issues/58:
# using Pkg
# Pkg.add(PackageSpec(url="https://github.com/dmlc/XGBoost.jl"))

using MLJRegistry
using MLJBase

#@test !(isempty(MLJRegistry.metadata()))

d = Dict()
d[:test] = Tuple{Union{Continuous,Missing},Discrete}
d["junk"] = Dict{Any,Any}("H" => Missing, :cross => "lemon", :t => :w, "r" => "r")
d["a"] = "b"
d[:f] = true
d["j"] = :post

@test MLJRegistry.decode_dic(MLJRegistry.encode_dic(d)) == d


text:
# Compatibility shims until users upgrade to new NNlib format
function conv!(y::CuArray{T}, x::CuArray{T}, w::CuArray{T}; pad=0, stride=1, flipkernel=0, dilation=1, kwargs...) where {T<:CUDNNFloat}
    cdims = DenseConvDims(x, w; padding=pad, stride=stride, flipkernel=flipkernel, dilation=dilation)
    return conv!(y, x, w, cdims; kwargs...)
end

function ∇conv_filter!(dw::CuArray{T}, dy::CuArray{T}, x::CuArray{T}; pad=0, stride=1, flipkernel=0, dilation=1, kwargs...) where {T<:CUDNNFloat}
    cdims = DenseConvDims(x, dw; padding=pad, stride=stride, flipkernel=flipkernel, dilation=dilation)
    # NOTE!!!  This compat shim re-arranges the argument order!
    return ∇conv_filter!(dw, x, dy, cdims; kwargs...)
end

function maxpool!(y::CuArray{T}, x::CuArray{T}, k; pad=map(_->0,k), stride=k) where {T<:CUDNNFloat}
    pdims = PoolDims(x, k; padding=pad, stride=stride)
    return maxpool!(y, x, pdims)
end

function meanpool!(y::CuArray{T}, x::CuArray{T}, k; pad=map(_->0,k), stride=k) where {T<:CUDNNFloat}
    pdims = PoolDims(x, k; padding=pad, stride=stride)
    return meanpool!(y, x, pdims)
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Test, Libdl, libblastrampoline_jll

@testset "libblastrampoline_jll" begin
    @test isa(Libdl.dlsym(Libdl.dlopen(:libblastrampoline), :dgemm_64_), Ptr{Nothing})
end

max_stars_repo_name: jonas-schulze/julia

max_issues_repo_name: jonas-schulze/julia

max_forks_repo_name: jonas-schulze/julia


text:
#
#--------------------------------------#
export mask
#--------------------------------------#
"""
 Mu = (R'*R)*u

 masks dirichlet boundary points
"""
function mask(u::Array
             ,M::Array)

if(length(M)==0) Mu = copy(u)
else             Mu = @. M*u
end

return Mu
end
#--------------------------------------#


text:
function solve()
    readline()
    as = [parse(BigInt, x) for x in split(readline())]
    p = prod(as)
    if p > 10^18
        -1
    else
        p
    end
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
using LinearAlgebra
using UnicodePlots
include("gemm.jl")

function plot_data(ns, data)
  scaling = 2 * ns.^3

  minv = floatmax()
  maxv = floatmin()
  for (s, v) = pairs(data)
    if s != :ref
    minv = min(minv, minimum(scaling ./ v.min))
    maxv = max(maxv, maximum(scaling ./ v.min))
  end
  end

  plt = lineplot(ns, scaling ./ data.ref.min, xlim = (ns[1], ns[end]),
                 ylim = (minv, maxv), title="max GFLOPS", xlabel="n",
                 ylabel="time", name = "ref")

  for (s, v) = pairs(data)
    if s != :ref
      lineplot!(plt, ns, scaling ./ v.min, name = string(s))
    end
  end
  plt
end

using BenchmarkTools

function run(ns, nsamples = 10)

  ref = (min=zeros(length(ns)), avg = zeros(length(ns)))
  jip_pji = (min=zeros(length(ns)), avg = zeros(length(ns)))
  pji = (min=zeros(length(ns)), avg = zeros(length(ns)))
  jpi = (min=zeros(length(ns)), avg = zeros(length(ns)))
  ji_4x4 = (min=zeros(length(ns)), avg = zeros(length(ns)))
  ji = (min=zeros(length(ns)), avg = zeros(length(ns)))

  for (i, n) = enumerate(ns)
    println()
    @show i,n
    m = k = n
    C = rand(m, n)
    A = rand(m, k)
    B = rand(k, n)
    C2 = copy(C)
    C1 = copy(C)

    t = @elapsed for n = 1:nsamples
      mul!(C2, A, B, 1, 1)
    end
    ref.min[i] = ref.avg[i] = t / nsamples
    @show ref.min[i], 2*m^3 / ref.min[i]

    #=
    C .= C1
    t = @elapsed for n = 1:nsamples
      mygemm_jpi(C, A, B)
    end
    jpi.min[i] = jpi.avg[i] = t / nsamples
    @show jpi.min[i], 2m^3/jpi.min[i]
    @show extrema(C - C2)
    =#

    #=
    t = @elapsed for n = 1:nsamples
      mygemm_pji(C, A, B)
    end
    pji.min[i] = pji.avg[i] = t / nsamples
    @show pji.min[i], 2m^3/pji.min[i]
    =#

    #=
    t = @elapsed for n = 1:nsamples
      mygemm_jip_pji(C, A, B, 4, 4)
    end
    jip_pji.min[i] = jip_pji.avg[i] = t / nsamples
    @show jip_pji.min[i], 2m^3/jip_pji.min[i]
    =#

    #=
    C .= C1
    t = @elapsed for n = 1:nsamples
      mygemm_ji_4x4(C, A, B)
    end
    ji_4x4.min[i] = ji_4x4.avg[i] = t / nsamples
    @show ji_4x4.min[i], 2m^3/ji_4x4.min[i]
    @show extrema(C - C2)
    =#

    C .= C1
    t = @elapsed for n = 1:nsamples
      mygemm_ji(C, A, B, Val(8), Val(4))
    end
    ji.min[i] = ji.avg[i] = t / nsamples
    @show ji.min[i], 2m^3/ji.min[i]
    @show extrema(C - C2)
  end

  data = (ref = ref,
          # jpi = jpi,
          # pji = pji,
          # jip_pji = jip_pji,
          # ji_4x4 = ji_4x4,
          ji = ji,
         )

  return data
end

ns = 48:48:1500
data = run(ns)
plt = plot_data(ns, data)

max_stars_repo_name: jkozdon/Julia-LAFF-On-PfHP

max_issues_repo_name: jkozdon/Julia-LAFF-On-PfHP

max_forks_repo_name: jkozdon/Julia-LAFF-On-PfHP


text:
"""
    struct CompatiblePaths
        S::States
        C::Vector{Node}
        Z::DecisionStrategy
        fixed::FixedPath
    end

CompatiblePaths type.
"""
struct CompatiblePaths
    S::States
    C::Vector{Node}
    Z::DecisionStrategy
    fixed::FixedPath
    function CompatiblePaths(S, C, Z, fixed)
        if !all(k∈Set(C) for k in keys(fixed))
            throw(DomainError("You can only fix chance states."))
        end
        new(S, C, Z, fixed)
    end
end

"""
    CompatiblePaths(diagram::InfluenceDiagram, Z::DecisionStrategy, fixed::FixedPath=Dict{Node, State}())

CompatiblePaths outer construction function. Interface for iterating over paths that are compatible and active given influence diagram and decision strategy.

1) Initialize path `s` of length `n`
2) Fill chance states `s[C]` by generating subpaths `paths(C)`
3) Fill decision states `s[D]` by decision strategy `Z` and path `s`

# Examples
```julia
for s in CompatiblePaths(diagram, Z)
    ...
end
```
"""
function CompatiblePaths(diagram::InfluenceDiagram, Z::DecisionStrategy, fixed::FixedPath=Dict{Node, State}())
    CompatiblePaths(diagram.S, diagram.C, Z, fixed)
end

function compatible_path(S::States, C::Vector{Node}, Z::DecisionStrategy, s_C::Path)
    s = Array{State}(undef, length(S))
    for (c, s_C_j) in zip(C, s_C)
        s[c] = s_C_j
    end
    for (d, I_d, Z_d) in zip(Z.D, Z.I_d, Z.Z_d)
        s[d] = Z_d((s[I_d]...,))
    end
    return (s...,)
end

function Base.iterate(S_Z::CompatiblePaths)
    if isempty(S_Z.fixed)
        iter = paths(S_Z.S[S_Z.C])
    else
        ks = sort(collect(keys(S_Z.fixed)))
        fixed = Dict{Node, State}(Node(i) => S_Z.fixed[k] for (i, k) in enumerate(S_Z.C) if k in ks)
        iter = paths(S_Z.S[S_Z.C], fixed)
    end
    next = iterate(iter)
    if next !== nothing
        s_C, state = next
        return (compatible_path(S_Z.S, S_Z.C, S_Z.Z, s_C), (iter, state))
    end
end

function Base.iterate(S_Z::CompatiblePaths, gen)
    iter, state = gen
    next = iterate(iter, state)
    if next !== nothing
        s_C, state = next
        return (compatible_path(S_Z.S, S_Z.C, S_Z.Z, s_C), (iter, state))
    end
end

Base.eltype(::Type{CompatiblePaths}) = Path
Base.length(S_Z::CompatiblePaths) = prod(S_Z.S[c] for c in S_Z.C)

"""
    struct UtilityDistribution
        u::Vector{Float64}
        p::Vector{Float64}
    end

UtilityDistribution type.

"""
struct UtilityDistribution
    u::Vector{Float64}
    p::Vector{Float64}
end

"""
    UtilityDistribution(diagram::InfluenceDiagram, Z::DecisionStrategy)

Construct the probability mass function for path utilities on paths that are compatible with given decision strategy.

# Examples
```julia
UtilityDistribution(diagram, Z)
```
"""
function UtilityDistribution(diagram::InfluenceDiagram, Z::DecisionStrategy)
    # Extract utilities and probabilities of active paths
    S_Z = CompatiblePaths(diagram, Z)
    utilities = Vector{Float64}(undef, length(S_Z))
    probabilities = Vector{Float64}(undef, length(S_Z))
    for (i, s) in enumerate(S_Z)
        utilities[i] = diagram.U(s)
        probabilities[i] = diagram.P(s)
    end

    # Filter zero probabilities
    nonzero = @. (!)(iszero(probabilities))
    utilities = utilities[nonzero]
    probabilities = probabilities[nonzero]

    # Sort by utilities
    perm = sortperm(utilities)
    u = utilities[perm]
    p = probabilities[perm]

    # Compute the probability mass function
    u2 = unique(u)
    p2 = similar(u2)
    j = 1
    p2[j] = p[1]
    for k in 2:length(u)
        if u[k] == u2[j]
            p2[j] += p[k]
        else
            j += 1
            p2[j] = p[k]
        end
    end

    UtilityDistribution(u2, p2)
end

"""
    struct StateProbabilities
        probs::Dict{Node, Vector{Float64}}
        fixed::FixedPath
    end

StateProbabilities type.
"""
struct StateProbabilities
    probs::Dict{Node, Vector{Float64}}
    fixed::FixedPath
end

"""
    StateProbabilities(diagram::InfluenceDiagram, Z::DecisionStrategy, node::Node, state::State, prior_probabilities::StateProbabilities)

Associate each node with array of conditional probabilities for each of its states occuring in compatible paths given
    fixed states and prior probability. Fix node and state using their indices.

# Examples
```julia
# Prior probabilities
prior_probabilities = StateProbabilities(diagram, Z)
StateProbabilities(diagram, Z, Node(2), State(1), prior_probabilities)
```
"""
function StateProbabilities(diagram::InfluenceDiagram, Z::DecisionStrategy, node::Node, state::State, prior_probabilities::StateProbabilities)
    prior = prior_probabilities.probs[node][state]
    fixed = deepcopy(prior_probabilities.fixed)

    push!(fixed, node => state)
    probs = Dict(i => zeros(diagram.S[i]) for i in 1:length(diagram.S))
    for s in CompatiblePaths(diagram, Z, fixed), i in 1:length(diagram.S)
        probs[i][s[i]] += diagram.P(s) / prior
    end
    StateProbabilities(probs, fixed)
end

"""
    StateProbabilities(diagram::InfluenceDiagram, Z::DecisionStrategy, node::Name, state::Name, prior_probabilities::StateProbabilities)

Associate each node with array of conditional probabilities for each of its states occuring in compatible paths given
    fixed states and prior probability. Fix node and state using their names.

# Examples
```julia
# Prior probabilities
prior_probabilities = StateProbabilities(diagram, Z)

# Select node and fix its state
node = "R"
state = "no test"
StateProbabilities(diagram, Z, node, state, prior_probabilities)
```
"""
function StateProbabilities(diagram::InfluenceDiagram, Z::DecisionStrategy, node::Name, state::Name, prior_probabilities::StateProbabilities)
    node_index = findfirst(j -> j ==node, diagram.Names)
    state_index = findfirst(j -> j == state, diagram.States[node_index])

    return StateProbabilities(diagram, Z, Node(node_index), State(state_index), prior_probabilities)
end

"""
    StateProbabilities(diagram::InfluenceDiagram, Z::DecisionStrategy)

Associate each node with array of probabilities for each of its states occuring in compatible paths.

# Examples
```julia
StateProbabilities(diagram, Z)
```
"""
function StateProbabilities(diagram::InfluenceDiagram, Z::DecisionStrategy)
    probs = Dict(i => zeros(diagram.S[i]) for i in 1:length(diagram.S))
    for s in CompatiblePaths(diagram, Z), i in 1:length(diagram.S)
        probs[i][s[i]] += diagram.P(s)
    end
    StateProbabilities(probs, Dict{Node, State}())
end

"""
    value_at_risk(U_distribution::UtilityDistribution, α::Float64)

Calculate value-at-risk.
"""
function value_at_risk(U_distribution::UtilityDistribution, α::Float64)
    @assert 0 ≤ α ≤ 1 "We should have 0 ≤ α ≤ 1."
    perm = sortperm(U_distribution.u)
    u, p = U_distribution.u[perm], U_distribution.p[perm]
    index = findfirst(x -> x≥α, cumsum(p))
    return if index === nothing; u[end] else u[index] end
end

"""
    conditional_value_at_risk(U_distribution::UtilityDistribution, α::Float64)

Calculate conditional value-at-risk.
"""
function conditional_value_at_risk(U_distribution::UtilityDistribution, α::Float64)
    x_α = value_at_risk(U_distribution, α)
    if iszero(α)
        return x_α
    else
        tail = U_distribution.u .≤ x_α
        return (sum(U_distribution.u[tail] .* U_distribution.p[tail]) - (sum(U_distribution.p[tail]) - α) * x_α) / α
    end
end


text:
__precompile__()

module AlphaStructures
	using LinearAlgebra
	using Delaunay # used in 3D triangulation --->> warning Test possible BUG in this package
	using Triangulate # used in 2D triangulation
	using Combinatorics, DataStructures
	using Distributed
	using SharedArrays

	#using MATLAB

	include("alpha_complex.jl")
	include("geometry.jl")
	include("delaunayTriangulation.jl")

	export alphaFilter, alphaSimplex, delaunayTriangulation
end


text:
using GraphMatFun
using Test
@testset "GraphMatFun.jl" begin
    include("bbc_sid_exp.jl")
    include("bbcs_cheb_exp.jl")
    include("eval_runerr.jl")
    include("code_gen_test.jl")
    include("compress.jl")
    include("degopt_poly.jl")
    include("degopt.jl")
    include("denman_beavers.jl")
    include("exp.jl")
    include("gauss_newton.jl")
    include("graph_ops.jl")
    include("import_export.jl")
    include("jac.jl")
    include("linear_fit.jl")
    include("merge_graphs.jl")
    include("newton_schulz.jl")
    include("opt_common.jl")
    include("polynomials_degopt.jl")
    include("polynomials.jl")
    include("rational.jl")
    include("sastre.jl")
    include("topo_order.jl")
end


text:
#---------------------Parse utils-----------------------------------------------------#
"""
Searches a directory for all files containing the key.
"""
searchdir(path::String, key) = filter(x -> occursin(key, x), readdir(path))

"""
Parse an array of strings into an array of a type.
"""
parse_string_array(T::Type, array) = map(x -> (v = tryparse(T, x); v==nothing ? 0.0 : v), array)

"""
Parse a line for occurrences of type T.
"""
parse_line(T::Type, line::String) = parse_string_array(T, split(line))

"""
Splits a line using arguments, then strips spaces from the splits.
"""
strip_split(line, args...) = strip.(split(line, args...))

function replace_multiple(str, replacements::Pair{String, String}...)
    tstr = deepcopy(str)
    for r in replacements
        tstr = replace(tstr, r)
    end
    return tstr
end

function cut_after(line, c)
    t = findfirst(isequal(c), line)
    if t == nothing
        return line
    elseif t == 1
        return ""
    else
        return line[1:t - 1]
    end
end

#--------------------------------------------------------------------------------------#
"""
Mutatatively applies the fermi level to all eigvals in the band. If fermi is a quantum espresso scf output file it will try to find it in there.
"""
function apply_fermi_level!(band::Band, fermi::Union{String,AbstractFloat})
    if typeof(fermi) == String
        fermi = qe_read_fermi_from_output(fermi)
    end
    for i = 1:size(band.eigvals)[1]
        band.eigvals[i] -= fermi
    end
end

#TODO: there is probably a slightly more optimal way for the frozen window, by checking max and min
#      such that every k-point has nbands inbetween.
function Emax(Emin, nbnd, bands)
    nbndfound = 0
    max = 0
    for b in bands
        if maximum(b.eigvals) >= Emin && nbndfound <= nbnd
            nbndfound += 1
            #maximum of allowed frozen window is the minimum of the first band>nbnd
            max = minimum(b.eigvals)-0.005
        end
    end

    nbndfound <= nbnd && error("Number of needed bands for the projections ($nbnd) exceeds the amount of bands starting from \nEmin=$Emin ($nbndfound).")
    return max
end

function wanenergyranges(Emin, nbnd, bands, Epad=5)
    max = Emax(Emin, nbnd, bands)
    (Emin - Epad, Emin, max, max + Epad)
end

"""
Applies the fermi level to all eigvals in the band. If fermi is a quantum espresso scf output file it will try to find it in there.
"""
function apply_fermi_level(band::Band, fermi)
    T = typeof(band.eigvals[1])
    if typeof(fermi) == String
        fermi = qe_read_fermi_from_output(fermi)
    end
    out = deepcopy(band)
    for i1 = 1:size(band.eigvals)[1]
        out.eigvals[i1] = band.eigvals[i1] - T(fermi)
    end
    return out
end

"""
    kgrid(na, nb, nc, input)

Returns an array of k-grid points that are equally spaced, input can be either `:wan` or `:nscf`, the returned grids are appropriate as inputs for wannier90 or an nscf calculation respectively.
"""
kgrid(na, nb, nc, ::Type{QE}) = reshape([[a, b, c, 1 / (na * nb * nc)] for a in collect(range(0, stop=1, length=na + 1))[1:end - 1], b in collect(range(0, stop=1, length=nb + 1))[1:end - 1], c in collect(range(0, stop=1, length=nc + 1))[1:end - 1]], (na * nb * nc))
kgrid(na, nb, nc, ::Type{Wannier90}) = reshape([[a, b, c] for a in collect(range(0, stop=1, length=na + 1))[1:end - 1], b in collect(range(0, stop=1, length=nb + 1))[1:end - 1], c in collect(range(0, stop=1, length=nc + 1))[1:end - 1]],(na * nb * nc))
kgrid(na, nb, nc, input::Symbol) = input==:wan ? kgrid(na, nb, nc, Wannier90) : kgrid(na, nb, nc, QE)
kgrid(na, nb, nc, input::DFInput{T}) where T = kgrid(na, nb, nc, T)

kakbkc(kgrid) = length.(unique.([[n[i] for n in kgrid] for i=1:3]))

function fort2julia(f_type)
    f_type = lowercase(f_type)
    if f_type == "real"
        return Float32
    elseif f_type == "real(kind=dp)"
        return Float64
    elseif f_type == "complex(kind=dp)"
        return Complex{Float64}
    elseif occursin("character", f_type)
        return String
    elseif f_type == "string"
        return String
    elseif f_type == "integer"
        return Int
    elseif f_type == "logical"
        return Bool
    elseif occursin(".D", f_type)
        return replace(f_type, "D" => "e")
    else
        return Nothing
    end
end

"""
It's like filter()[1].
"""
function getfirst(f::Function, A)
    for el in A
        if f(el)
            return el
        end
    end
end

"""
    parse_block(f, types...; to_strip=',')

Takes the specified types and parses each line into the types.
When it finds a line where it cannot match all the types, it stops and returns  the parsed values.
The split and strip keywords let the user specify how to first split the line, then strip the splits from the strip char.
"""
function parse_block(f, types...; to_strip=',')
    output = []
    len_typ = length(types)
    while !eof(f)
        line = strip.(split(readline(f)), to_strip)
        len_lin = length(line)
        if isempty(line)
            continue
        end
        i,j = 1,1
        tmp = []
        while i <= len_typ && j <= len_lin
            typ = types[i]
            l   = line[j]
            try
                t   = Meta.parse(l)
                if typeof(t) == typ
                    push!(tmp, t)
                    i+=1
                    j+=1
                else
                    j+=1
                end
            catch
                j+=1
            end
        end
        if length(tmp) < length(types)
            return output
        end
        push!(output, Tuple{types...}(tmp))
    end
    return output
end

yesterday() = today() - Day(1)
lastweek()  = today() - Week(1)
lastmonth() = today() - Month(1)


text:
module TestSpaces

using Test
using TreeParzen

space = Dict(
    :a => HP.Choice(:a, [:a, :b, :c]),
    :b => HP.Choice(:b, ["a", "b", "c"])
)

# just asking, isn't good enough
trial = ask(space)
sample = trial.hyperparams

@test sample[:a] in (:a, :b, :c)
@test sample[:b] in ("a", "b", "c")

# gotta ask both ways
trials = TreeParzen.Trials.Trial[]
tell!(trials, trial, 2.)
trial = ask(space)
tell!(trials, trial, 2.6)

tp_trial = ask(space, trials, TreeParzen.Config())
tp_sample = tp_trial.hyperparams

# same holds true but we have to make sure we can sample from strings literals via TP asks
@test tp_sample[:a] in (:a, :b, :c)
@test tp_sample[:b] in ("a", "b", "c")

# This is to test support for using inputs other than Dicts as space definitions
direct_delayed_space = HP.Choice(:mychoice,
    [
        HP.Uniform(:a1, 0., 1.),
        HP.Uniform(:a2, 1., 2.),
        HP.Uniform(:a3, 2., 3.),
    ]
)

config = Config()
trials = [ask(direct_delayed_space) for i in 1:config.random_trials]
# fill them in
[tell!(t, 1.) for t in trials]
# actually, we just wanna see that this stuff doesn't vom
newtrial = ask(direct_delayed_space, trials, config)
# Well and check that the hyperparams is a float ....
@test newtrial.hyperparams isa Float64

direct_dict_space = HP.Choice(:mychoice,
    [
        Dict(:a => HP.Uniform(:a1, 0., 1.), :b => HP.Uniform(:b1, 0., 1.),),
        Dict(:a => HP.Uniform(:a2, 0., 1.), :b => HP.Uniform(:b2, 0., 1.),),
        Dict(:a => HP.Uniform(:a3, 0., 1.), :b => HP.Uniform(:b3, 0., 1.),),
    ]
)

trials = [ask(direct_dict_space) for i in 1:config.random_trials]
# fill them in
[tell!(t, 1.) for t in trials]
# actually, we just wanna see that this stuff doesn't vom
newtrial = ask(direct_dict_space, trials, config)
# Well and check that the hyperparams is a dict ....
@test newtrial.hyperparams isa Dict
@test haskey(newtrial.hyperparams, :a)
@test haskey(newtrial.hyperparams, :b)

direct_array_space = [
    HP.Uniform(:a1, 0., 1.),
    HP.Uniform(:a2, 1., 2.),
    HP.Uniform(:a3, 2., 3.),
]

trials = [ask(direct_array_space) for i in 1:config.random_trials]
# fill them in
[tell!(t, 1.) for t in trials]
# actually, we just wanna see that this stuff doesn't vom
newtrial = ask(direct_array_space, trials, config)
# Well and check that the hyperparams is an array ....
@test newtrial.hyperparams isa Vector
@test length(newtrial.hyperparams) == 3
@test 0 <= newtrial.hyperparams[1] <= 1
@test 1 <= newtrial.hyperparams[2] <= 2
@test 2 <= newtrial.hyperparams[3] <= 3

end

true


text:
using Relief, Test

# Run tests
include("./relief_test.jl")
include("./relieff_test.jl")
include("./reliefseq_test.jl")
include("./reliefmss_test.jl")
include("./surfstar_test.jl")
include("./surf_test.jl")
include("./multisurfstar_test.jl")
include("./multisurf_test.jl")
include("./swrfstar_test.jl")
include("./boostedsurf_test.jl")
include("./iterative_relief_test.jl")
include("./irelief_test.jl")
include("./ecrelieff_test.jl")
include("./vlsrelief_test.jl")
include("./turf_test.jl")