text:
function test_z_boundary_condition_simple(arch, T, field_name, bctype, bc, Nx, Ny)
    Nz = 16
    model = Model(N=(Nx, Ny, Nz), L=(0.1, 0.2, 0.3), arch=arch, float_type=T)

    bc = BoundaryCondition(bctype, bc)
    bcs = getfield(model.boundary_conditions, field_name)
    bcs.z.top = bc

    time_step!(model, 1, 1e-16)

    typeof(model) <: Model
end

function test_z_boundary_condition_top_bottom_alias(arch, TF, field_name)
    N = 16
    model = Model(N=(N, N, N), L=(0.1, 0.2, 0.3), arch=arch, float_type=TF)

    bcval = 1.0
    bcs = getfield(model.boundary_conditions, field_name)
    bcs.z.top = BoundaryCondition(Value, bcval)
    bcs.z.bottom = BoundaryCondition(Value, -bcval)

    time_step!(model, 1, 1e-16)

    getbc(bcs.z.top) == bcval && getbc(bcs.z.bottom) == -bcval
end

function test_z_boundary_condition_array(arch, T, field_name)
    Nx = Ny = Nz = 16
    model = Model(N=(Nx, Ny, Nz), L=(0.1, 0.2, 0.3), arch=arch, float_type=T)

    bcarray = rand(T, Nx, Ny)
    if arch == GPU()
        bcarray = CuArray(bcarray)
    end

    bcs = getfield(model.boundary_conditions, field_name)
    bcs.z.top = BoundaryCondition(Value, bcarray)

    time_step!(model, 1, 1e-16)

    bcs.z.top[1, 2] == bcarray[1, 2]
end

function test_flux_budget(arch, TF, field_name)
    N, κ, Lz = 16, 1, 0.7
    model = Model(N=(N, N, N), L=(1, 1, Lz), ν=κ, κ=κ,
        arch=arch, float_type=TF, eos=LinearEquationOfState(βS=0, βT=0))

    if field_name ∈ (:u, :v, :w)
        field = getfield(model.velocities, field_name)
    else
        field = getfield(model.tracers, field_name)
    end

    @. field.data = 0

    bottom_flux = TF(0.3)
    flux_bc = BoundaryCondition(Flux, bottom_flux)
    bcs = getfield(model.boundary_conditions, field_name)
    bcs.z.bottom = flux_bc

    mean_init = mean(field.data)

    τκ = Lz^2 / κ # diffusion time-scale
    Δt = 1e-6 * τκ # time-step much less than diffusion time-scale
    Nt = 100

    time_step!(model, Nt, Δt)

    # budget is Lz * ∂<ϕ>/∂t = -Δflux = -top_flux/Lz (left) + bottom_flux/Lz (right);
    # therefore <ϕ> = bottom_flux * t / Lz
    isapprox(mean(field.data) - mean_init, bottom_flux * model.clock.time / Lz)
end


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*--------------------------------------------------------------------
#* Author: Gerson Balbueno Bicca
#* $Id: DoublePipeIncr.mso								$
#*------------------------------------------------------------------
type Physical_Properties
	Physical_Properties()=begin
		new(
			Properties_In_Out(Dict{Symbol,Any}(
				:Brief=>"Properties at Inlet Stream",
				:Symbol=>"^{in}"
			)),
			Properties_Average(Dict{Symbol,Any}(
				:Brief=>"Properties at Average Temperature",
				:Symbol=>"^{avg}"
			)),
			Properties_In_Out(Dict{Symbol,Any}(
				:Brief=>"Properties at Outlet Stream",
				:Symbol=>"^{out}"
			)),
			Properties_Wall(Dict{Symbol,Any}(
				:Brief=>"Properties at Wall Temperature",
				:Symbol=>"^{wall}"
			)),
			[:Inlet,:Average,:Outlet,:Wall,]
		)
	end
	Inlet::Properties_In_Out
	Average::Properties_Average
	Outlet::Properties_In_Out
	Wall::Properties_Wall
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export Physical_Properties
function atributes(in::Physical_Properties,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=false
	fields[:Brief]="to be documented"
	fields[:Info]="to be documented"
	drive!(fields,_)
	return fields
end
Physical_Properties(_::Dict{Symbol,Any})=begin
	newModel=Physical_Properties()
	newModel.attributes=atributes(newModel,_)
	newModel
end

ext: jl


text:
# ridge regression as example of Bayesian estimation
using Plots, LinearAlgebra
function main()
n = 100
a = 0.95 # controls collinearity: close to 0 indep close to 1 highly collinear
bs = zeros(1000,6)
truebeta = ones(3,1)
for i = 1:1000
	# generate regressors that are collinear
	x1 = rand(n,1)
	x2 = a*x1 + (1-a)*randn(n,1)
	x = [ones(n,1) x1 x2]
	# true params and dep var
	e = randn(n,1)
	y = x*truebeta+e;
	# OLS
	bols = x\y
	# ridge regression
	k = 1.0
	yy = [y; k*zeros(3,1)]
	xx = [x; k*Matrix{Float64}(I, 3, 3)]
	bridge = xx\yy
	# store results
	bs[i,:] = [bols' bridge']

end
# histograms to compare
bs = bs .- [truebeta' truebeta'] # take away true, so centered around 0
p1 = histogram(bs[:,2], title="OLS")
p2 = histogram(bs[:,5],title="Ridge")
plot(p1,p2,layout=(1,2), legend=false)
#savefig("RidgeExample.svg")
gui()
return
end
main()

max_stars_repo_name: nilshg/Econometrics

max_issues_repo_name: nilshg/Econometrics

max_forks_repo_name: nilshg/Econometrics


text:
using Documenter, QueryOperators

makedocs(
	modules = [QueryOperators],
	sitename = "QueryOperators.jl",
	analytics="UA-132838790-1",
	pages = [
        "Introduction" => "index.md"
    ]
)

deploydocs(
    repo = "github.com/queryverse/QueryOperators.jl.git"
)


text:
type Template{T<:Function}
    width::Int
    render::T
end

Template{T}(w, f::T) = Template{T}(w, f)

const templates = Dict(
    :text => Template(10, function(buf, data::DataSet)
        error("TODO")
    end),

    :list => Template(10, function(buf, data::DataSet)
        error("TODO")
    end),

    :title => Template(10, function(buf, data::DataSet)
        error("TODO")
    end),

    :scatter => Template(10, function(buf, data::DataSet)
        id = uid()
        println(buffer, """

""")
        println(buffer, """""")
    end),

    :line => Template(10, function(buf, data::DataSet)
        error("TODO")
    end),

    :pie => Template(10, function(buf, data::DataSet)
        error("TODO")
    end),

    :bar => Template(10, function(buf, data::DataSet)
        error("TODO")
    end),

    :box => Template(10, function(buf, data::DataSet)
        error("TODO")
    end),

    :heat => Template(10, function(buf, data::DataSet)
        error("TODO")
    end)
)


text:
function DSO_min_LoadShift(hours,
	nodes::Nodes,
	av_el_price,
	av_rd_price_dict,
	x_ag,
	x_re,
	x_ind,
	s_ag,
	s_re,
	s_ind,
	isworkhour
	)

# create subset of nodes
T = hours
J = nodes.id
Sbase = nodes.bmva

# calculate average nodal load
loadmean = Dict(i => sum(nodes.load[i][hours]/Sbase)/(length(hours)) for i in keys(nodes.load))

DSO_mod = JuMP.Model(with_optimizer(Gurobi.Optimizer))

	### VARIABLE DECLARATION
		### Variables for sectoral shift
			@variable(DSO_mod, ΔD_up_ag[T, J] >= 0);
			@variable(DSO_mod, ΔD_up_re[T, J] >= 0);
			@variable(DSO_mod, ΔD_up_ind[T, J] >= 0);
			@variable(DSO_mod, ΔD_dn_ag[T, J] >= 0);
			@variable(DSO_mod, ΔD_dn_re[T, J] >= 0);
			@variable(DSO_mod, ΔD_dn_ind[T, J] >= 0);
		### Variables for total shift
			@variable(DSO_mod, ΔD_up[T, J] >= 0);
			@variable(DSO_mod, ΔD_dn[T, J] >= 0);
		### Variable for new Loadprofile
			@variable(DSO_mod, Load[T, J] >= 0);

	### CONSTRAINT DECLARATION
		### Constraints for total shift
			@constraint(DSO_mod, ShiftUpTotal[j = J, t = T], ΔD_up[t, j] == ΔD_up_ag[t, j] + ΔD_up_re[t, j] + ΔD_up_ind[t, j]);
			@constraint(DSO_mod, ShiftDnTotal[j = J, t = T], ΔD_dn[t, j] == ΔD_dn_ag[t, j] + ΔD_dn_re[t, j] + ΔD_dn_ind[t, j]);
		### Constraints for agricultural, residential shift
			@constraint(DSO_mod, ShiftAgUpMax[j = J, t = T], ΔD_up_ag[t, j] <= x_ag*s_ag*loadmean[j]/Sbase);
			@constraint(DSO_mod, ShiftAgDnMax[j = J, t = T], ΔD_dn_ag[t, j] <= x_ag*s_ag*loadmean[j]/Sbase);
			@constraint(DSO_mod, ShiftReUpMax[j = J, t = T], ΔD_up_re[t, j] <= x_re*s_re*loadmean[j]/Sbase);
			@constraint(DSO_mod, ShiftReDnMax[j = J, t = T], ΔD_dn_re[t, j] <= x_re*s_re*loadmean[j]/Sbase);
		### Constraints for industrial shift
			# if mod(hours[24],168) < 120
			@constraint(DSO_mod, ShiftIndUpMax[j = J, t = T], ΔD_up_ind[t, j] <= x_ind*s_ind*loadmean[j]/Sbase*isworkhour[mod(t,24)]);
			@constraint(DSO_mod, ShiftIndDnMax[j = J, t = T], ΔD_dn_ind[t, j] <= x_ind*s_ind*loadmean[j]/Sbase*isworkhour[mod(t,24)]);
			#	else
			#@constraint(DSO_mod, ShiftIndUpMax[j = J, t = T], ΔD_up_ind[t, j] == 0);
			#@constraint(DSO_mod, ShiftIndDnMax[j = J, t = T], ΔD_dn_ind[t, j] == 0);
			# end
		### Constraints for ShiftBalances
			@constraint(DSO_mod, ShiftBalAg[j = J],
			sum(ΔD_up_ag[t, j] - ΔD_dn_ag[t, j] for t in hours) == 0);
			@constraint(DSO_mod, ShiftBalRe[j = J],
			sum(ΔD_up_re[t, j] - ΔD_dn_re[t, j] for t in hours) == 0);
			@constraint(DSO_mod, ShiftBalInd[j = J],
			sum(ΔD_up_ind[t, j] - ΔD_dn_ind[t, j] for t in hours) == 0);

		### Constraint for new Loadprofile
			@constraint(DSO_mod, ShiftedLoad[j = J, t = T],
			Load[t, j] == nodes.load[j][t]/Sbase + ΔD_up[t, j] - ΔD_dn[t, j]
			);

	### Objective Function
			@objective(DSO_mod, Min,
			sum(sum(av_el_price*Load[t, j] + av_rd_price_dict[j][t]*(ΔD_up[t, j]-ΔD_dn[t, j]) for t in T) for j in J)
				)

		# Initiate optimisation process
		JuMP.optimize!(DSO_mod)

ShiftUp = JuMP.value.(ΔD_up) * Sbase
ShiftDn = JuMP.value.(ΔD_dn) * Sbase
ShiftAgUp = JuMP.value.(ΔD_up_ag) * Sbase
ShiftAgDn = JuMP.value.(ΔD_dn_ag) * Sbase
ShiftReUp = JuMP.value.(ΔD_up_re) * Sbase
ShiftReDn = JuMP.value.(ΔD_dn_re) * Sbase
ShiftIndUp = JuMP.value.(ΔD_up_ind) * Sbase
ShiftIndDn = JuMP.value.(ΔD_dn_ind) * Sbase
Load_shifted_opt = JuMP.value.(Load) * Sbase

return(Load_shifted_opt,
	ShiftUp,
	ShiftDn,
)


text:
@testset "Gabor" begin
    v1 = rand(3)
    v2 = rand(3)
    ell = rand()
    p = rand()
    k = gaborkernel(;
        sqexponential_transform=ScaleTransform(inv(ell)),
        cosine_transform=ScaleTransform(inv(p)),
    )
    @test k isa KernelProduct{
        <:Tuple{
            TransformedKernel{<:SqExponentialKernel,<:ScaleTransform},
            TransformedKernel{<:CosineKernel,<:ScaleTransform},
        },
    }
    @test k.kernels[1].transform.s[1] == inv(ell)
    @test k.kernels[2].transform.s[1] == inv(p)

    k_manual = exp(-sqeuclidean(v1, v2) / (2 * ell^2)) * cospi(euclidean(v1, v2) / p)
    @test k_manual ≈ k(v1, v2) atol = 1e-5

    lhs_manual = (SqExponentialKernel() ∘ ScaleTransform(1 / ell))(v1, v2)
    rhs_manual = (CosineKernel() ∘ ScaleTransform(1 / p))(v1, v2)
    @test lhs_manual * rhs_manual ≈ k(v1, v2) atol = 1e-5

    @test gaborkernel() isa KernelProduct{<:Tuple{<:SqExponentialKernel,<:CosineKernel}}

    test_ADs(
        x -> gaborkernel(;
            sqexponential_transform=ScaleTransform(x[1]),
            cosine_transform=ScaleTransform(x[2]),
        ),
        [ell, p],
    )
end


text:
module OSQPMathProgBaseInterface

using LinearAlgebra
using SparseArrays
using OSQP: Model, Results, setup!, solve!, update!, clean!, update_settings!, warm_start!
using MathProgBase

struct OSQPSolver <: MathProgBase.AbstractMathProgSolver
    settings::Dict{Symbol,Any}
end

OSQPSolver(; kwargs...) = OSQPSolver(Dict{Symbol, Any}(k => v for (k, v) in kwargs))

mutable struct OSQPMathProgModel <: MathProgBase.AbstractLinearQuadraticModel
    # Problem settings
    settings::Dict{Symbol, Any}

    # Inner OSQP model structure
    inner::Model

    # Setup needed
    perform_setup::Bool

    # Problem data in the form
    #
    # minimize    (1/2) x' P x + q' x
    # subject to  lb <= A x <= ub
    #             l <= x <= u
    sense::Symbol
    P::SparseMatrixCSC{Float64, Int}
    q::Vector{Float64}
    A::SparseMatrixCSC{Float64, Int}
    lb::Vector{Float64}
    ub::Vector{Float64}
    l::Vector{Float64}
    u::Vector{Float64}

    # Warm start auxiliary variables
    x_ws::Vector{Float64}  # Warm start vector
    # Call warm start function in optimize!.
    # NB. Warmstart is enabled internally by default when changing problem data
    run_warmstart::Bool

    # Solved
    solved::Bool

    # Results
    results::Results

    function OSQPMathProgModel(settings::AbstractDict{Symbol, Any})
        P = spzeros(Float64, 0, 0)
        q = Float64[]
        A = spzeros(Float64, 0, 0)
        lb = Float64[]
        ub = Float64[]
        l = Float64[]
        u = Float64[]
        perform_setup = true
        x_ws = Float64[]
        run_warmstart = false
        solved = false
        new(settings, Model(), perform_setup,
            :Min, P, q, A, lb, ub, l, u,
            x_ws, run_warmstart, solved) # leave other fields undefined
    end
end

function MathProgBase.numvar(model::OSQPMathProgModel)
    # Check problem dimensions
    if size(model.P, 2) == 0
        if size(model.A, 2) != 0
            n = size(model.A, 2)
        else
            error("the model has no variables")
        end
    else
        n = size(model.P, 1)
    end

    return n
end

MathProgBase.numconstr(model::OSQPMathProgModel) = size(model.A, 1)

# Initialize
function initialize_dimensions!(model::OSQPMathProgModel, n, m)
    # Initialize empty vectors and matrices with the right dimensions
    model.P = spzeros(n, n)
    model.q = zeros(n)
    model.A = spzeros(m, n)
    model.lb = zeros(m)
    model.ub = zeros(m)
    model.l = zeros(n)
    model.u = zeros(n)
    model.x_ws = zeros(n)

    model
end

# Maybe add resetproblem function to remove results and reset perform_setup?
checksolved(model::OSQPMathProgModel) = model.solved || error("Model has not been solved.")

# Reset problem when setup has to be performed again
resetproblem(model::OSQPMathProgModel) = (model.solved = false)  # Clear model.results!

# http://mathprogbasejl.readthedocs.io/en/latest/solverinterface.html
MathProgBase.LinearQuadraticModel(solver::OSQPSolver) = OSQPMathProgModel(solver.settings)
MathProgBase.getsolution(model::OSQPMathProgModel) = (checksolved(model); model.results.x)
function MathProgBase.getobjval(model::OSQPMathProgModel)
    checksolved(model)

    if model.sense == :Min
        return model.results.info.obj_val
    else
        return -model.results.info.obj_val
    end
end

function MathProgBase.optimize!(model::OSQPMathProgModel)
    # Perform setup only if necessary
    if model.perform_setup
        n = MathProgBase.numvar(model)
        setup!(model.inner; P = model.P, q = model.q,
               A = [model.A; sparse(I, n, n)],
               l = [model.lb; model.l],
               u = [model.ub; model.u],
               model.settings...)
        model.perform_setup = false  # No longer needed to perform setup
    end

    # Call internal warm start function if needed
    if model.run_warmstart
        warm_start!(model.inner, x=model.x_ws)
    end

    # Solve the problem and store the results
    model.results = solve!(model.inner)

    # Problem is now solved
    model.solved = true
end

function MathProgBase.status(model::OSQPMathProgModel)::Symbol
    checksolved(model)
    osqpstatus = model.results.info.status
    status = osqpstatus # if OSQP status can't be mapped to a standard return value, just return as is

    # map to standard return status values:
    osqpstatus == :Solved && (status = :Optimal)
    osqpstatus == :Max_iter_reached && (status = :UserLimit)
    osqpstatus == :Interrupted && (status = :UserLimit) # following Gurobi.jl
    osqpstatus == :Primal_infeasible && (status = :Infeasible)
    osqpstatus == :Primal_infeasible_inaccurate && (status = :Infeasible)
    osqpstatus == :Dual_infeasible && (status = :Unbounded)
    osqpstatus == :Dual_infeasible_inaccurate && (status = :Unbounded)
    osqpstatus == :Time_limit_reached && (status = :UserLimit)
    osqpstatus == :Non_convex && (status = :Error)

    return status
end

# TODO: getobjbound
# TODO: getobjgap

MathProgBase.getrawsolver(model::OSQPMathProgModel) = model.inner
MathProgBase.getsolvetime(model::OSQPMathProgModel) = (checksolved(model); model.results.info.run_time)

function MathProgBase.setsense!(model::OSQPMathProgModel, sense::Symbol)
    if sense in [:Min, :Max]
        if (sense == :Max) & (model.sense != :Max)
            model.q .= .-model.q
            model.P .= .-model.P

            # Update problem data
            if !model.perform_setup
                update!(model.inner, q=model.q)
                update!(model.inner, Px=triu(model.P).nzval)
            end
        end
        model.sense = sense
    else
        error("sense not recognized")
    end

    resetproblem(model)
end

MathProgBase.getsense(model::OSQPMathProgModel) = model.sense
MathProgBase.freemodel!(model::OSQPMathProgModel) = clean!(model.inner)

function Base.copy(model::OSQPMathProgModel)
    ret = OSQPMathProgModel(model.settings)
    initialize_dimensions!(ret, MathProgBase.numvar(model), MathProgBase.numconstr(model))
    ret.sense = model.sense
    copyto!(ret.P, model.P)
    copyto!(ret.q, model.q)
    copyto!(ret.A, model.A)
    copyto!(ret.lb, model.lb)
    copyto!(ret.ub, model.ub)
    copyto!(ret.l, model.l)
    copyto!(ret.u, model.u)
    ret
end

MathProgBase.setvartype!(model::OSQPMathProgModel, v::Vector{Symbol}) = any(x -> x != :Cont, v) && error("OSQP only supports continuous variables.")
MathProgBase.getvartype(model::OSQPMathProgModel) = fill(:Cont, MathProgBase.numvar(model))

# Set solver independent parameters: Silent is the only supported one for now
function MathProgBase.setparameters!(x::OSQPMathProgModel; Silent = nothing, TimeLimit = nothing)
    if Silent != nothing
        Silent::Bool
        x.settings[:verbose] = !Silent

        # Update silent setting if setup has already been performed
        if !x.perform_setup
            update_settings!(x.inner, verbose=Silent)
        end
    end

    if TimeLimit != nothing
        TimeLimit::Float
        x.settings[:time_limit] = !TimeLimit

        # Update time_limit setting if setup has already been performed
        if !x.perform_setup
            update_settings!(x.inner, time_limit=TimeLimit)
        end
    end

    x
end

# Do not update the problem instance settings using OSQP internal functions if the x is OSQPSolver and not OSQPMathProgModel
function MathProgBase.setparameters!(x::OSQPSolver; Silent = nothing, TimeLimit = nothing)
    if Silent != nothing
        Silent::Bool
        x.settings[:verbose] = !Silent
    end
    if TimeLimit != nothing
        TimeLimit::Float
        x.settings[:time_limit] = TimeLimit
    end

    x
end

"""
    setwarmstart!(model::OSQPMathProgModel, x)

Warm start the solver solution with the primal variable `x`.

NB. The `warm_start!` function supports setting also the dual variables but in MathProgBase only
the primal variable warm start is supported. Note that OSQP performs warm starting automatically
when parts of the problem data change and a new optimize! is called.
"""
function MathProgBase.setwarmstart!(model::OSQPMathProgModel, x)
    copyto!(model.x_ws, x)
    model.run_warmstart = true
end

# http://mathprogbasejl.readthedocs.io/en/latest/lpqcqp.html#linearquadratic-models
function MathProgBase.loadproblem!(model::OSQPMathProgModel, A, l, u, c, lb, ub, sense)

    (m, n) = size(A)  # Get problem dimensions
    initialize_dimensions!(model, n, m)  # Initialize problem dimensions

    if sense == :Min
        copyto!(model.q, c)
    elseif sense == :Max
        model.q .= .-c
    else
        error("Objective sense not recognized")
    end

    # Fix model sense
    model.sense = sense

    # Copy variables
    copyto!(model.lb, lb)
    copyto!(model.ub, ub)
    copyto!(model.l, l)
    copyto!(model.u, u)
    copyto!(model.A, sparse(A))  # Sparsify matrix A
    model
end

MathProgBase.getvarLB(model::OSQPMathProgModel) = model.l
function MathProgBase.setvarLB!(model::OSQPMathProgModel, l)
    # Copy new variable lower bounds
    copyto!(model.l, l)

    # Update lower bounds in the model
    if !model.perform_setup
        update!(model.inner, l=[model.lb; model.l])
    end

    # Reset problem solution
    resetproblem(model)

    return model
end

MathProgBase.getvarUB(model::OSQPMathProgModel) = model.u
function MathProgBase.setvarUB!(model::OSQPMathProgModel, u)
    # Copy new variable upper bounds
    copyto!(model.u, u)

    # Update lower bounds in the model
    if !model.perform_setup
        update!(model.inner, u=[model.ub; model.u])
    end

    # Reset problem solution
    resetproblem(model)

    return model
end

MathProgBase.getconstrLB(model::OSQPMathProgModel) = model.lb
function MathProgBase.setconstrLB!(model::OSQPMathProgModel, lb)
    # Copy new constraints lower bounds
    copyto!(model.lb, lb)

    # Update lower bounds in the model
    if !model.perform_setup
        update!(model.inner, l=[model.lb; model.l])
    end

    # Reset problem solution
    resetproblem(model)

    return model
end

MathProgBase.getconstrUB(model::OSQPMathProgModel) = model.ub
function MathProgBase.setconstrUB!(model::OSQPMathProgModel, ub)
    # Copy new constraints upper bounds
    copyto!(model.ub, ub)

    # Update upper bounds in the model
    if !model.perform_setup
        update!(model.inner, u=[model.ub; model.u])
    end

    # Reset problem solution
    resetproblem(model)

    return model
end

function MathProgBase.setobj!(model::OSQPMathProgModel, c)
    # Copy cost
    copyto!(model.q, c)

    # Negate cost if necessary
    if model.sense == :Max
        model.q .= .- model.q
    end

    if !model.perform_setup
        update!(model.inner, q=model.q)
    end

    # Reset problem solution
    resetproblem(model)

    return model

end

function MathProgBase.getobj(model::OSQPMathProgModel)
    if model.sense == :Min
        return model.q
    else
        return -model.q
    end
end

MathProgBase.getconstrmatrix(model::OSQPMathProgModel) = model.A

function MathProgBase.addvar!(model::OSQPMathProgModel, constridx, constrcoef, l, u, objcoef)

    # Get bounds if they are not set
    ((l == nothing) && (li = -Inf)) || (li = l)
    ((u == nothing) && (ui = Inf)) || (ui = u)

    # Change cost P, q
    model.q = [model.q; objcoef]
    if model.sense == :Max
        model.q[end] = -model.q[end]
    end
    model.P = blockdiag(model.P, spzeros(1,1))

    # Change constraints
    if !isempty(constrcoef) && !isempty(constridx)
        # Update A
        m = MathProgBase.numconstr(model)
        constr_xi = sparsevec(constridx, constrcoef, m)
        model.A = [model.A constr_xi]
    end

    # Update l and u
    model.l = [model.l; li]
    model.u = [model.u; ui]

    # Change ws
    model.x_ws = [model.x_ws; 0]

    # With new variable we need to perform setup
    model.perform_setup = true

    # Reset problem solution
    resetproblem(model)

end
MathProgBase.addvar!(model::OSQPMathProgModel, l, u, objcoef) = MathProgBase.addvar!(model, [], [], l, u, objcoef)

function MathProgBase.addconstr!(model::OSQPMathProgModel, varidx, coef, lb, ub)
    # Construct sparse vector with the new constraint
    n = MathProgBase.numvar(model)
    ai = sparsevec(varidx, coef, n)

    # Augment A, lb, ub
    if isempty(model.A)
        model.A = sparse(ai')
    else
        model.A = [model.A; ai']
    end
    model.lb = [model.lb; lb]
    model.ub = [model.ub; ub]

    # With the new constraint we must perform setup
    model.perform_setup = true

    # Reset problem solution
    resetproblem(model)

end

# TODO: delvars!
# TODO: delconstrs!
# TODO: changecoeffs!

MathProgBase.numlinconstr(model::OSQPMathProgModel) = MathProgBase.numconstr(model)
MathProgBase.getconstrsolution(model::OSQPMathProgModel) = model.A * MathProgBase.getsolution(model)
function MathProgBase.getreducedcosts(model::OSQPMathProgModel)
    checksolved(model)
    if model.sense == :Min
        return -model.results.y[end-(MathProgBase.numvar(model)-1):end]
    else
        return model.results.y[end-(MathProgBase.numvar(model)-1):end]
    end
end

function MathProgBase.getconstrduals(model::OSQPMathProgModel)
    checksolved(model)
    if model.sense == :Min
        return -model.results.y[1:MathProgBase.numconstr(model)]
    else
        return model.results.y[1:MathProgBase.numconstr(model)]
    end
end

function MathProgBase.getinfeasibilityray(model::OSQPMathProgModel)
    checksolved(model)

    if model.results.info.status in [:Primal_infeasible, :Primal_infeasible_inaccurate]
        m = MathProgBase.numconstr(model)
        # Returns infeasibility ray taking into account both bounds and constraints
        ray = -model.results.prim_inf_cert[1:MathProgBase.numconstr(model)]
    else
        error("Problem not infeasible")
    end
end

function MathProgBase.getunboundedray(model::OSQPMathProgModel)
    checksolved(model)

    if model.results.info.status in [:Dual_infeasible, :Dual_infeasible_inaccurate]
        return model.results.dual_inf_cert
    else
        error("Problem not unbounded")
    end
end

# TODO: getbasis

# http://mathprogbasejl.readthedocs.io/en/latest/lpqcqp.html#quadratic-programming
MathProgBase.numquadconstr(model::OSQPMathProgModel) = 0

# setquadobj!(model::OSQPMathProgModel, Q::Matrix) = ((Qi, Qj, Qx) = findnz(Q); setquadobj!(model, Qi, Qj, Qx))
function MathProgBase.setquadobj!(model::OSQPMathProgModel, rowidx::Vector, colidx::Vector, quadval::Vector)
    nterms = length(quadval)
    @boundscheck length(rowidx) == nterms || error()
    @boundscheck length(colidx) == nterms || error()

    # Check if only the values have changed
    Pi, Pj, Px = findnz(model.P)
    if (rowidx == Pi) & (colidx == Pj) & !model.perform_setup
        if model.sense == :Max
            # Update matrix in MathProgBase model
            copyto!(model.P.nzval, -Px)
        else
            # Update matrix in MathProgBase model
            copyto!(model.P.nzval, Px)
        end
        # Update only nonzeros of P
        update!(model.inner, Px=model.P.nzval)

        # Reset solution status
        resetproblem(model)

        return model
    end

    # Create sparse matrix from indices
    # zero out coeffs
    for i = 1 : nterms
        @inbounds row = rowidx[i]
        @inbounds col = colidx[i]
        model.P[row, col] = 0
    end

    # add new coeffs and symmetrize
    for i = 1 : nterms
        @inbounds row = rowidx[i]
        @inbounds col = colidx[i]
        @inbounds val = quadval[i]
        model.P[row, col] += val
        model.P[col, row] = model.P[row, col]
    end

    # Change sign if maximizing
    if model.sense == :Max
        model.P .= .-model.P
    end

    # Need new setup when we set a new P
    model.perform_setup = true

    # Reset problem solution
    resetproblem(model)

    model
end

# Note: skipping quadconstr methods

end # module


text:
@testset "1599.maximum-profit-of-operating-a-centennial-wheel.jl" begin
    @test min_operations_max_profit([8, 3], 5, 6) == 3
    @test min_operations_max_profit([4, 1, 1, 1], 5, 6) == 1
    @test min_operations_max_profit([4, 0, 0, 4], 5, 6) == 4
    @test min_operations_max_profit([3, 4, 0, 5, 1], 1, 92) == -1
    @test min_operations_max_profit([10,10,6,4,7], 3, 8) == 9
end


text:
using Documenter, AoC2021

for symbol in names(AoC2021; all=true)
    object = eval(:(AoC2021.$symbol))
    if isa(object, Module) && object != AoC2021
        qualifiedname = string(object)
        name = split(qualifiedname, '.')[2]
        content = """
        # $qualifiedname

        ```@index
        Modules = [$qualifiedname]
        Order = [:function, :type]
        ```

        ```@autodocs
        Modules = [$qualifiedname]
        Order = [:function, :type]
        ```
        """
        open(joinpath(@__DIR__, "src", lowercase(name) * ".md"), "w") do file
            write(file, content)
        end
    end
end

makedocs(
    sitename = "Advent of Code 2021",
    modules = [AoC2021],
    strict = true
)


text:
function circuitfunction(Circuit)
    for (f,t) in zip(["-","[",",","]"],["+","((",")^-1+(",")^-1)^-1"])
        Circuit=replace(Circuit,f=>t)
    end
    for I in 2:-1:1
        Es = eachmatch(Regex("([CLRPW])([0-9]){$(I)}"),Circuit)
            for e in Es
                match = e.match
                if match[1] == 'C'
                 Circuit = replace(Circuit,match=>"(1/(2im*π*f*"*"T"*"))")
                elseif match[1] == 'R'
                 Circuit = replace(Circuit,match=>"T")
                elseif match[1] == 'L'
                 Circuit = replace(Circuit,match=>"(2im*π*f*"*"T"*")")
             elseif match[1] == 'P'
                 Circuit = replace(Circuit,match=>"T*(2*π*f)^(-N)"*"*(cos((π*N)*0.5)-sin((π*T)*0.5)im)")
             elseif match[1] == 'W'
                 Circuit = replace(Circuit,match=>"T*(2*π*f)^(-0.5)"*"*(cos((π*0.5)*0.5)-sin((π*0.5)*0.5)im)")
            end
        end
end
new_circuit = ""
counter = 1
for i in Circuit
    if i == 'T'
        new_circuit = new_circuit*"T["*string(counter)*"]"
        counter += 1
    elseif i == 'N'
        new_circuit = new_circuit*"T["*string(counter)*"]"
    else
        new_circuit = new_circuit*i
    end
end

    Circuit_expression = Meta.parse(new_circuit)
    return mk_function([:T,:f],[],Circuit_expression)
end

function tree_to_function(tree_array)
    dictionary_of_calculations = Dict('R'=>"T",'C'=>"(1/(2im*π*f*"*"T"*"))",'L'=>"(2im*π*f*"*"T"*")",'P'=>"T[1]*(2*π*f)^(-T[2])"*"*(cos((π*T[2])*0.5)-sin((π*T[2])*0.5)im)")
    tree = copy(tree_array)
    nodecount = length(tree)
    essential_info = [[node.ParentIndex,node.Type,node.Index] for node in tree]
    for i in 1:((length(essential_info)-1)/2)
        parent1,type1,index1 = pop!(essential_info)
        parent2,type2,index2 = pop!(essential_info)
        type1 = haskey(dictionary_of_calculations,type1) ? replace(dictionary_of_calculations[type1],"T"=>"T["*string(index1)*']') : type1
        type2 = haskey(dictionary_of_calculations,type2) ? replace(dictionary_of_calculations[type2],"T"=>"T["*string(index2)*']') : type2
        operation = essential_info[parent1][2]
        if operation == '+'
            essential_info[parent1][2] = type1*'+'*type2
        elseif operation == '-'
            essential_info[parent1][2] = "((("*type1*")^-1)+("*type2*")^-1)^-1"
        end
    end
    CircuitExpression = Meta.parse(essential_info[1][2])
return mk_function([:T,:f],[],CircuitExpression)
end

karva_to_function(karva::String) = tree_to_function(karva_to_tree(karva))
karva_to_function(circuit::Circuit) = karva_to_function(circuit.karva)

function flatten(params)
    new_array = Float64[]
    for param in params
        if length(param) == 1
            push!(new_array,param)
        else
            push!(new_array,param[1])
            push!(new_array,param[2])
        end
    end
    return new_array
end

function get_target_impedance(circuit,circuit_parameters,frequency=1000)
    circfunc = circuitfunction(circuit)
    target_impedance = circfunc(flatten(circuit_parameters),frequency)
   return target_impedance
end


text:
using Hecke, Test

DEFAULT_NPROCS = 4

# Test if _adjust_path works on Windows
x = Hecke._adjust_path("GrpAb/Elem")
y = joinpath(Hecke.pkgdir, "test", "$x.jl")
@test isfile(y)

################################################################################
#
#  Analyze the arguments
#
################################################################################

# Is short?
short_test = false

if "short" in ARGS || get(ENV, "HECKE_TESTSHORT", "false") in ["1", "true"]
  global short_test = true
end

# Is long?
long_test = false

if "long" in ARGS || get(ENV, "HECKE_TESTLONG", "false") in ["1", "true"]
  global long_test = true
end

# Is GAP there?
with_gap = false

push!(Base.LOAD_PATH, "@v#.#")

try
  using GAP
  println("Found GAP. Add FieldFactory.jl to the long tests")
  global with_gap = true
catch e
  if !(isa(e, ArgumentError))
    rethrow(e)
  else
    println("using GAP failed.")
  end
end

# Parallel
isparallel = false
n_procs = 0

fl = get(ENV, "HECKE_TEST_PARALLEL", "false")
if fl != "false"
  isparallel = true
  n_procs = parse(Int, fl)
else
  for a in ARGS
    r = match(r"j([0-9])*", a)
    if r === nothing
      continue
    else
      global isparallel = true
      if r.captures[1] === nothing
        global n_procs = DEFAULT_NPROCS
      else
        global n_procs = parse(Int, r.captures[1])
      end
      @assert n_procs > 0 "Number of processes ($(n_procs)) must be > 0"
    end
  end
end

fl = get(ENV, "CI", "false")
if fl === "true"
  isparallel = true
  # CPU_THREADS reports number of logical cores (including hyperthreading)
  # So be pessimistic and divide by 2
  n_procs = div(Sys.CPU_THREADS, 1)
end

if VERSION < v"1.5.0"
  if isparallel
    @warn "Parallel testing might be hanging on julia < 1.5.0"
  end
end

if v"1.3" <= VERSION < v"1.4.0"
  @warn "Parallel testing disabled on julia 1.3"
  isparallel = false
  n_procs = 0
end

# Now collect the tests we want to run

const exclude = ["setup.jl", "runtests.jl", "parallel.jl", "testdefs.jl", "FieldFactory.jl"]

test_directory = joinpath(@__DIR__)

const long_tests = String[]

if with_gap
  push!(long_tests, "FieldFactory.jl")
end

tests = String[]

for t in readdir(test_directory)
  if !isfile(joinpath(test_directory, t))
    continue
  end

  if t in long_tests
    if long_test
      push!(tests, t)
    else
      continue
    end
  end

  if !(t in exclude)
    push!(tests, t)
  end
end

# Put FieldFactory.jl and QuadForm.jl at the beginning, because they take the
# longest

for s in ["QuadForm.jl", "FieldFactory.jl"]
  if s in tests
    i = findfirst(isequal(s), tests)
    deleteat!(tests, i)
    pushfirst!(tests, s)
  end
end

test_path(test) = joinpath(@__DIR__, test)

@info "Hecke test setup"
@info "long_test : $long_test"
@info "short_test: $short_test"
if isparallel
  @info "parallel  : $isparallel ($(n_procs))"
else
  @info "parallel  : $isparallel"
end
@info "tests     : $tests"

if short_test
  include("setup.jl")
  # Short tests are always running on one machine
  @info "Running short tests"
  include(joinpath("..", "system", "precompile.jl"))
else
  if !isparallel
    # We are not short
    k, a = quadratic_field(5)
    @test fmpz(1) - a == -(a - 1)
    @test 1 - a == -(a - 1)

    include("setup.jl")

    for t in tests
      @time include(test_path(t))
    end
  else
    # Now we are parallel
    include("parallel.jl")
  end

  #try
  #  using Polymake
  #  @time include("AlgAssRelOrd/Eichler.jl")
  #catch e
  #  if !(isa(e, ArgumentError))
  #    rethrow(e)
  #  else
  #    println("using Polymake failed. Not running sophisticated norm equation tests")
  #  end
  #end
end


text:
struct Triangle <: Shape
    id::Int64
    p1::Point3
    p2::Point3
    p3::Point3
    obj_to_world::Transformation
    world_to_obj::Transformation
end

Triangle(T::Transformation) = Triangle(1, Point3(0,0,0), Point3(1,0,0), Point3(0,1,0), T, inv(T))
Triangle(p1::Point3, p2::Point3, p3::Point3, T::Transformation) = Triangle(1, p1, p2, p3, T, inv(T))
Triangle(id::Integer, p1::Point3, p2::Point3, p3::Point3, T::Transformation) = Triangle(id, p1, p2, p3, T, inv(T))
Triangle(p1::Point3, p2::Point3, p3::Point3) = Triangle(p1, p2, p3, Transformation())

can_intersect(T::Triangle) = true
area(T::Triangle) = 0.5 * norm(cross(T.p2 - T.p1, T.p3 - T.p1))

function obj_bounds(T::Triangle)
    return BoundingBox(T.p1, T.p2, T.p3)
end

function world_bounds(T::Triangle)
    T.obj_to_world(obj_bounds(T))
end

function (T::Transformation)(C::Triangle)
    T2 = T * C.obj_to_world
    Triangle(C.id, C.p1, C.p2, C.p3, T2, inv(T2))
end

function shape_intersect(R::Ray, T::Triangle)
    ray = T.world_to_obj(R)

    e1 = T.p2 - T.p1
    e2 = T.p3 - T.p1

    s1 = cross(ray.direction, e2)
    divisor = dot(s1, e1)
    if divisor == 0.0
        return nothing, NaN, NaN
    end
    invd = 1.0 / divisor

    d = ray.origin - T.p1
    b1 = dot(d, s1) * invd
    if b1 < 0.0 || b1 > 1.0
        return nothing, NaN, NaN
    end

    s2 = cross(d, e1)
    b2 = dot(ray.direction, s2) * invd
    if b2 < 0.0 || b1+b2 > 1.0
        return nothing, NaN, NaN
    end

    t = dot(e2, s2) * invd
    if t < ray.tmin || t > ray.tmax
        return nothing, NaN, NaN
    end

    P = ray(t)

    n = normalize(Normal3(cross(e1, e2)))

    n = flip_normal(n, ray.direction)

    DG = DifferentialGeometry(T.obj_to_world(P), T.obj_to_world(n), T.obj_to_world(e1))

    return DG, t, 5e-4 * t
end

# Fast true/false intersection test between a ray and a triangle
function intersectP(R::Ray, T::Triangle)
    ray = T.world_to_obj(R)

    e1 = T.p2 - T.p1
    e2 = T.p3 - T.p1

    s1 = cross(ray.direction, e2)
    divisor = dot(s1, e1)
    if divisor == 0.0
        return false
    end
    invd = 1.0 / divisor

    d = ray.origin - T.p1
    b1 = dot(d, s1) * invd
    if b1 < 0.0 || b1 > 1.0
        return false
    end

    s2 = cross(d, e1)
    b2 = dot(ray.direction, s2) * invd
    if b2 < 0.0 || b1+b2 > 1.0
        return false
    end

    t = dot(e2, s2) * invd
    if t < ray.tmin || t > ray.tmax
        return false
    end

    return true
end

max_stars_repo_name: dronir/GRAYTR

max_issues_repo_name: dronir/GRAYTR

max_forks_repo_name: dronir/GRAYTR


text:
module FitMethodsTest
using OnlineStats, Base.Test

@testset "fit! methods" begin
    @testset "Single Observation" begin
        o = Mean()
        fit!(o, randn())
        fit!(o, randn(), rand())
        @test nobs(o) == 2

        o = Means(2)
        fit!(o, randn(2))
        fit!(o, randn(2), rand())
        @test nobs(o) == 2

        o = LinReg(2)
        fit!(o, randn(2), randn())
        fit!(o, randn(2), randn(), rand())
        @test nobs(o) == 2
    end
    @testset "Multiple Observations" begin
        o = Mean()
        fit!(o, randn(100))
        fit!(o, randn(100), rand(100))
        fit!(o, randn(100), .1)
        @test nobs(o) == 300

        o = Means(2)
        fit!(o, randn(100, 2))
        fit!(o, randn(100, 2), rand(100))
        fit!(o, randn(100, 2), .1)
        @test nobs(o) == 300

        o = LinReg(2)
        fit!(o, randn(100, 2), randn(100))
        fit!(o, randn(100, 2), randn(100), rand(100))
        fit!(o, randn(100, 2), randn(100), .1)
        @test nobs(o) == 300
    end
    @testset "maprows" begin
        o = Mean()
        o2 = Mean()
        y = randn(100)
        fit!(o, y)
        maprows(10, y) do yi
            fit!(o2, yi)
        end
        @test mean(o) == mean(o2)
    end
end

end  # module


text:
export
    DCPPowerModel, StandardDCPForm,
    DCPLLPowerModel, StandardDCPLLForm

""
abstract type AbstractDCPForm <: AbstractPowerFormulation end

""
abstract type StandardDCPForm <: AbstractDCPForm end

""
const DCPPowerModel = GenericPowerModel{StandardDCPForm}

"default DC constructor"
DCPPowerModel(data::Dict{String,Any}; kwargs...) =
    GenericPowerModel(data, StandardDCPForm; kwargs...)

""
function variable_voltage(pm::GenericPowerModel{T}, n::Int=pm.cnw; kwargs...) where T <: AbstractDCPForm
    variable_voltage_angle(pm, n; kwargs...)
end

"nothing to add, there are no voltage variables on branches"
function variable_voltage_ne(pm::GenericPowerModel{T}, n::Int=pm.cnw; kwargs...) where T <: AbstractDCPForm
end

"dc models ignore reactive power flows"
function variable_reactive_generation(pm::GenericPowerModel{T}, n::Int=pm.cnw; bounded = true) where T <: AbstractDCPForm
end

"dc models ignore reactive power flows"
function variable_reactive_branch_flow(pm::GenericPowerModel{T}, n::Int=pm.cnw; bounded = true) where T <: AbstractDCPForm
end

"dc models ignore reactive power flows"
function variable_reactive_branch_flow_ne(pm::GenericPowerModel{T}, n::Int=pm.cnw) where T <: AbstractDCPForm
end

""
function variable_active_branch_flow(pm::GenericPowerModel{T}, n::Int=pm.cnw; bounded = true) where T <: StandardDCPForm
    if bounded
        pm.var[:nw][n][:p] = @variable(pm.model,
            [(l,i,j) in pm.ref[:nw][n][:arcs_from]], basename="$(n)_p",
            lowerbound = -pm.ref[:nw][n][:branch][l]["rate_a"],
            upperbound =  pm.ref[:nw][n][:branch][l]["rate_a"],
            start = getstart(pm.ref[:nw][n][:branch], l, "p_start")
        )
    else
        pm.var[:nw][n][:p] = @variable(pm.model,
            [(l,i,j) in pm.ref[:nw][n][:arcs_from]], basename="$(n)_p",
            start = getstart(pm.ref[:nw][n][:branch], l, "p_start")
        )
    end

    # this explicit type erasure is necessary
    p_expr = Dict{Any,Any}([((l,i,j), pm.var[:nw][n][:p][(l,i,j)]) for (l,i,j) in pm.ref[:nw][n][:arcs_from]])
    p_expr = merge(p_expr, Dict([((l,j,i), -1.0*pm.var[:nw][n][:p][(l,i,j)]) for (l,i,j) in pm.ref[:nw][n][:arcs_from]]))
    pm.var[:nw][n][:p] = p_expr
end

""
function variable_active_branch_flow_ne(pm::GenericPowerModel{T}, n::Int=pm.cnw) where T <: StandardDCPForm
    pm.var[:nw][n][:p_ne] = @variable(pm.model,
        [(l,i,j) in pm.ref[:nw][n][:ne_arcs_from]], basename="$(n)_p_ne",
        lowerbound = -pm.ref[:nw][n][:ne_branch][l]["rate_a"],
        upperbound =  pm.ref[:nw][n][:ne_branch][l]["rate_a"],
        start = getstart(pm.ref[:nw][n][:ne_branch], l, "p_start")
    )

    # this explicit type erasure is necessary
    p_ne_expr = Dict{Any,Any}([((l,i,j), 1.0*pm.var[:nw][n][:p_ne][(l,i,j)]) for (l,i,j) in pm.ref[:nw][n][:ne_arcs_from]])
    p_ne_expr = merge(p_ne_expr, Dict([((l,j,i), -1.0*pm.var[:nw][n][:p_ne][(l,i,j)]) for (l,i,j) in pm.ref[:nw][n][:ne_arcs_from]]))
    pm.var[:nw][n][:p_ne] = p_ne_expr
end

"do nothing, this model does not have complex voltage variables"
function constraint_voltage(pm::GenericPowerModel{T}, n::Int) where T <: AbstractDCPForm
end

"do nothing, this model does not have complex voltage variables"
function constraint_voltage_ne(pm::GenericPowerModel{T}, n::Int) where T <: AbstractDCPForm
end

"do nothing, this model does not have voltage variables"
function constraint_voltage_magnitude_setpoint(pm::GenericPowerModel{T}, n::Int, i, vm) where T <: AbstractDCPForm
end

"do nothing, this model does not have reactive variables"
function constraint_reactive_gen_setpoint(pm::GenericPowerModel{T}, n::Int, i, qg) where T <: AbstractDCPForm
end

""
function constraint_kcl_shunt(pm::GenericPowerModel{T}, n::Int, i, bus_arcs, bus_arcs_dc, bus_gens, pd, qd, gs, bs) where T <: AbstractDCPForm
    pg = pm.var[:nw][n][:pg]
    p = pm.var[:nw][n][:p]
    p_dc = pm.var[:nw][n][:p_dc]

    pm.con[:nw][n][:kcl_p][i] = @constraint(pm.model, sum(p[a] for a in bus_arcs) + sum(p_dc[a_dc] for a_dc in bus_arcs_dc) == sum(pg[g] for g in bus_gens) - pd - gs*1.0^2)
    # omit reactive constraint
end

""
function constraint_kcl_shunt_ne(pm::GenericPowerModel{T}, n::Int, i, bus_arcs, bus_arcs_dc, bus_arcs_ne, bus_gens, pd, qd, gs, bs) where T <: AbstractDCPForm
    pg = pm.var[:nw][n][:pg]
    p = pm.var[:nw][n][:p]
    p_ne = pm.var[:nw][n][:p_ne]
    p_dc = pm.var[:nw][n][:p_dc]

    @constraint(pm.model, sum(p[a] for a in bus_arcs) + sum(p_ne[a] for a in bus_arcs_ne) + sum(p_dc[a_dc] for a_dc in bus_arcs_dc) == sum(pg[g] for g in bus_gens) - pd - gs*1.0^2)
end

"""
Creates Ohms constraints (yt post fix indicates that Y and T values are in rectangular form)

```
p[f_idx] == -b*(t[f_bus] - t[t_bus])
```
"""
function constraint_ohms_yt_from(pm::GenericPowerModel{T}, n::Int, f_bus, t_bus, f_idx, t_idx, g, b, c, tr, ti, tm) where T <: AbstractDCPForm
    p_fr = pm.var[:nw][n][:p][f_idx]
    va_fr = pm.var[:nw][n][:va][f_bus]
    va_to = pm.var[:nw][n][:va][t_bus]

    @constraint(pm.model, p_fr == -b*(va_fr - va_to))
    # omit reactive constraint
end

"Do nothing, this model is symmetric"
function constraint_ohms_yt_to(pm::GenericPowerModel{T}, n::Int, f_bus, t_bus, f_idx, t_idx, g, b, c, tr, ti, tm) where T <: AbstractDCPForm
end

function constraint_ohms_yt_from_ne(pm::GenericPowerModel{T}, n::Int, i, f_bus, t_bus, f_idx, t_idx, g, b, c, tr, ti, tm, vad_min, vad_max) where T <: AbstractDCPForm
    p_fr = pm.var[:nw][n][:p_ne][f_idx]
    va_fr = pm.var[:nw][n][:va][f_bus]
    va_to = pm.var[:nw][n][:va][t_bus]
    z = pm.var[:nw][n][:branch_ne][i]

    @constraint(pm.model, p_fr <= -b*(va_fr - va_to + vad_max*(1-z)) )
    @constraint(pm.model, p_fr >= -b*(va_fr - va_to + vad_min*(1-z)) )
end

"Do nothing, this model is symmetric"
function constraint_ohms_yt_to_ne(pm::GenericPowerModel{T}, n::Int, i, f_bus, t_bus, f_idx, t_idx, g, b, c, tr, ti, tm, vad_min, vad_max) where T <: AbstractDCPForm
end

"`-rate_a <= p[f_idx] <= rate_a`"
function constraint_thermal_limit_from(pm::GenericPowerModel{T}, n::Int, f_idx, rate_a) where T <: AbstractDCPForm
    p_fr = pm.con[:nw][n][:sm_fr][f_idx[1]] = pm.var[:nw][n][:p][f_idx]
    getlowerbound(p_fr) < -rate_a && setlowerbound(p_fr, -rate_a)
    getupperbound(p_fr) > rate_a && setupperbound(p_fr, rate_a)
end

"Do nothing, this model is symmetric"
function constraint_thermal_limit_to(pm::GenericPowerModel{T}, n::Int, t_idx, rate_a) where T <: AbstractDCPForm
end

""
function add_bus_voltage_setpoint(sol, pm::GenericPowerModel{T}) where T <: AbstractDCPForm
    add_setpoint(sol, pm, "bus", "vm", :vm; default_value = (item) -> 1)
    add_setpoint(sol, pm, "bus", "va", :va)
end

""
function variable_voltage_on_off(pm::GenericPowerModel{T}, n::Int=pm.cnw; kwargs...) where T <: AbstractDCPForm
    variable_voltage_angle(pm, n; kwargs...)
end

"do nothing, this model does not have complex voltage variables"
function constraint_voltage_on_off(pm::GenericPowerModel{T}, n::Int) where T <: AbstractDCPForm
end

"`-b*(t[f_bus] - t[t_bus] + vad_min*(1-branch_z[i])) <= p[f_idx] <= -b*(t[f_bus] - t[t_bus] + vad_max*(1-branch_z[i]))`"
function constraint_ohms_yt_from_on_off(pm::GenericPowerModel{T}, n::Int, i, f_bus, t_bus, f_idx, t_idx, g, b, c, tr, ti, tm, vad_min, vad_max) where T <: AbstractDCPForm
    p_fr = pm.var[:nw][n][:p][f_idx]
    va_fr = pm.var[:nw][n][:va][f_bus]
    va_to = pm.var[:nw][n][:va][t_bus]
    z = pm.var[:nw][n][:branch_z][i]

    @constraint(pm.model, p_fr <= -b*(va_fr - va_to + vad_max*(1-z)) )
    @constraint(pm.model, p_fr >= -b*(va_fr - va_to + vad_min*(1-z)) )
end

"Do nothing, this model is symmetric"
function constraint_ohms_yt_to_on_off(pm::GenericPowerModel{T}, n::Int, i, f_bus, t_bus, f_idx, t_idx, g, b, c, tr, ti, tm, vad_min, vad_max) where T <: AbstractDCPForm
end

"""
Generic on/off thermal limit constraint

```
-rate_a*branch_z[i] <= p[f_idx] <=  rate_a*branch_z[i]
```
"""
function constraint_thermal_limit_from_on_off(pm::GenericPowerModel{T}, n::Int, i, f_idx, rate_a) where T <: AbstractDCPForm
    p_fr = pm.var[:nw][n][:p][f_idx]
    z = pm.var[:nw][n][:branch_z][i]

    @constraint(pm.model, p_fr <=  rate_a*z)
    @constraint(pm.model, p_fr >= -rate_a*z)
end

"""
Generic on/off thermal limit constraint

```
-rate_a*branch_ne[i] <= p_ne[f_idx] <=  rate_a*branch_ne[i]
```
"""
function constraint_thermal_limit_from_ne(pm::GenericPowerModel{T}, n::Int, i, f_idx, rate_a) where T <: AbstractDCPForm
    p_fr = pm.var[:nw][n][:p_ne][f_idx]
    z = pm.var[:nw][n][:branch_ne][i]

    @constraint(pm.model, p_fr <=  rate_a*z)
    @constraint(pm.model, p_fr >= -rate_a*z)
end

"nothing to do, from handles both sides"
function constraint_thermal_limit_to_on_off(pm::GenericPowerModel{T}, n::Int, i, t_idx, rate_a) where T <: AbstractDCPForm
end

"nothing to do, from handles both sides"
function constraint_thermal_limit_to_ne(pm::GenericPowerModel{T}, n::Int, i, t_idx, rate_a) where T <: AbstractDCPForm
end

"`angmin*branch_z[i] + vad_min*(1-branch_z[i]) <= t[f_bus] - t[t_bus] <= angmax*branch_z[i] + vad_max*(1-branch_z[i])`"
function constraint_voltage_angle_difference_on_off(pm::GenericPowerModel{T}, n::Int, i, f_bus, t_bus, angmin, angmax, vad_min, vad_max) where T <: AbstractDCPForm
    va_fr = pm.var[:nw][n][:va][f_bus]
    va_to = pm.var[:nw][n][:va][t_bus]
    z = pm.var[:nw][n][:branch_z][i]

    @constraint(pm.model, va_fr - va_to <= angmax*z + vad_max*(1-z))
    @constraint(pm.model, va_fr - va_to >= angmin*z + vad_min*(1-z))
end

"`angmin*branch_ne[i] + vad_min*(1-branch_ne[i]) <= t[f_bus] - t[t_bus] <= angmax*branch_ne[i] + vad_max*(1-branch_ne[i])`"
function constraint_voltage_angle_difference_ne(pm::GenericPowerModel{T}, n::Int, i, f_bus, t_bus, angmin, angmax, vad_min, vad_max) where T <: AbstractDCPForm
    va_fr = pm.var[:nw][n][:va][f_bus]
    va_to = pm.var[:nw][n][:va][t_bus]
    z = pm.var[:nw][n][:branch_ne][i]

    @constraint(pm.model, va_fr - va_to <= angmax*z + vad_max*(1-z))
    @constraint(pm.model, va_fr - va_to >= angmin*z + vad_min*(1-z))
end

""
abstract type AbstractDCPLLForm <: AbstractDCPForm end

""
abstract type StandardDCPLLForm <: AbstractDCPLLForm end

""
const DCPLLPowerModel = GenericPowerModel{StandardDCPLLForm}

"default DC constructor"
DCPLLPowerModel(data::Dict{String,Any}; kwargs...) = GenericPowerModel(data, StandardDCPLLForm; kwargs...)

"`sum(p[a] for a in bus_arcs) + sum(p_dc[a_dc] for a_dc in bus_arcs_dc)== sum(pg[g] for g in bus_gens) - pd - gs*1.0^2`"
function constraint_kcl_shunt(pm::GenericPowerModel{T}, n::Int, i, bus_arcs, bus_arcs_dc, bus_gens, pd, qd, gs, bs) where T <: AbstractDCPLLForm
    pg = pm.var[:nw][n][:pg]
    p = pm.var[:nw][n][:p]
    p_dc = pm.var[:nw][n][:p_dc]

    @constraint(pm.model, sum(p[a] for a in bus_arcs) + sum(p_dc[a_dc] for a_dc in bus_arcs_dc) == sum(pg[g] for g in bus_gens) - pd - gs*1.0^2i)
end

"""
Creates Ohms constraints (yt post fix indicates that Y and T values are in rectangular form)
"""
function constraint_ohms_yt_from(pm::GenericPowerModel{T}, n::Int, f_bus, t_bus, f_idx, t_idx, g, b, c, tr, ti, tm) where T <: AbstractDCPLLForm
    p_fr = pm.var[:nw][n][:p][f_idx]
    p_to = pm.var[:nw][n][:p][t_idx]
    va_fr = pm.var[:nw][n][:va][f_bus]
    va_to = pm.var[:nw][n][:va][t_bus]

    @constraint(pm.model, p_fr == -b*(va_fr - va_to))
end

"""
Creates Ohms constraints (yt post fix indicates that Y and T values are in rectangular form)
"""
function constraint_ohms_yt_to(pm::GenericPowerModel{T}, n::Int, f_bus, t_bus, f_idx, t_idx, g, b, c, tr, ti, tm) where T <: AbstractDCPLLForm
    p_fr = pm.var[:nw][n][:p][f_idx]
    p_to = pm.var[:nw][n][:p][t_idx]
    va_fr = pm.var[:nw][n][:va][f_bus]
    va_to = pm.var[:nw][n][:va][t_bus]

    r = g/(g^2 + b^2)
    @constraint(pm.model, p_fr + p_to >= r*(p_fr^2))
end

"""
"""
function constraint_ohms_yt_from_on_off(pm::GenericPowerModel{T}, n::Int, i, f_bus, t_bus, f_idx, t_idx, g, b, c, tr, ti, tm, vad_min, vad_max) where T <: AbstractDCPLLForm
    p_fr = pm.var[:nw][n][:p][f_idx]
    p_to = pm.var[:nw][n][:p][t_idx]
    va_fr = pm.var[:nw][n][:va][f_bus]
    va_to = pm.var[:nw][n][:va][t_bus]
    z = pm.var[:nw][n][:branch_z][i]

    @constraint(pm.model, p_fr <= -b*(va_fr - va_to + vad_max*(1-z)) )
    @constraint(pm.model, p_fr >= -b*(va_fr - va_to + vad_min*(1-z)) )
end

"""
"""
function constraint_ohms_yt_to_on_off(pm::GenericPowerModel{T}, n::Int, i, f_bus, t_bus, f_idx, t_idx, g, b, c, tr, ti, tm, vad_min, vad_max) where T <: AbstractDCPLLForm
    p_fr = pm.var[:nw][n][:p][f_idx]
    p_to = pm.var[:nw][n][:p][t_idx]
    va_fr = pm.var[:nw][n][:va][f_bus]
    va_to = pm.var[:nw][n][:va][t_bus]
    z = pm.var[:nw][n][:branch_z][i]

    r = g/(g^2 + b^2)
    t_m = max(abs(vad_min),abs(vad_max))
    @constraint(pm.model, p_fr + p_to >= r*( (-b*(va_fr - va_to))^2 - (-b*(t_m))^2*(1-z) ) )
    @constraint(pm.model, p_fr + p_to >= 0)
end

"""
"""
function constraint_ohms_yt_from_ne(pm::GenericPowerModel{T}, n::Int, i, f_bus, t_bus, f_idx, t_idx, g, b, c, tr, ti, tm, vad_min, vad_max) where T <: AbstractDCPLLForm
    p_fr = pm.var[:nw][n][:p_ne][f_idx]
    p_to = pm.var[:nw][n][:p_ne][t_idx]
    va_fr = pm.var[:nw][n][:va][f_bus]
    va_to = pm.var[:nw][n][:va][t_bus]
    z = pm.var[:nw][n][:branch_ne][i]

    @constraint(pm.model, p_fr <= -b*(va_fr - va_to + vad_max*(1-z)) )
    @constraint(pm.model, p_fr >= -b*(va_fr - va_to + vad_min*(1-z)) )
end

"""
"""
function constraint_ohms_yt_to_ne(pm::GenericPowerModel{T}, n::Int, i, f_bus, t_bus, f_idx, t_idx, g, b, c, tr, ti, tm, vad_min, vad_max) where T <: AbstractDCPLLForm
    p_fr = pm.var[:nw][n][:p_ne][f_idx]
    p_to = pm.var[:nw][n][:p_ne][t_idx]
    va_fr = pm.var[:nw][n][:va][f_bus]
    va_to = pm.var[:nw][n][:va][t_bus]
    z = pm.var[:nw][n][:branch_ne][i]

    r = g/(g^2 + b^2)
    t_m = max(abs(vad_min),abs(vad_max))
    @constraint(pm.model, p_fr + p_to >= r*( (-b*(va_fr - va_to))^2 - (-b*(t_m))^2*(1-z) ) )
    @constraint(pm.model, p_fr + p_to >= 0)
end

"`-rate_a*branch_z[i] <= p[t_idx] <= rate_a*branch_z[i]`"
function constraint_thermal_limit_to_on_off(pm::GenericPowerModel{T}, n::Int, i, t_idx, rate_a) where T <: AbstractDCPLLForm
    p_to = pm.var[:nw][n][:p][t_idx]
    z = pm.var[:nw][n][:branch_z][i]

    @constraint(pm.model, p_to <=  rate_a*z)
    @constraint(pm.model, p_to >= -rate_a*z)
end


text:
path = dirname(@__FILE__)

# Set up arguments
global m = AnSchorfheide(testing = true)

global system = JLD2.jldopen("$path/../reference/forecast_args.jld2","r") do file
    read(file, "system")
end

# Run impulse responses
states, obs, pseudo = impulse_responses(m, system)

# Compare to expected output
exp_states, exp_obs, exp_pseudo =
    JLD2.jldopen("$path/../reference/impulse_responses_out.jld2", "r") do file
        read(file, "exp_states"), read(file, "exp_obs"), read(file, "exp_pseudo")
    end

@testset "Compare irfs to expected output" begin
    @test @test_matrix_approx_eq exp_states states
    @test @test_matrix_approx_eq exp_obs    obs
    @test @test_matrix_approx_eq exp_pseudo pseudo
end

# Test impulse response methods for:
# (1) Choosing a set of shocks and their shock values
# (2) Choosing a shock, a state, and the value the state should hit in period 1 based on
# an impulse from that shock at time 1.
# (3) Same as (2) but with an observable.

# Setup
endo   = m.endogenous_states
obsvs  = m.observables
exo    = m.exogenous_shocks
horizon = 10

# Testing both for exact output match and for certain values to match to be
# explicit about the expected outcome of the impulse response

# Shock
shock_names  = [:z_sh, :g_sh]
shock_values = [1., 1.]
states, obs, pseudo = impulse_responses(m, system, horizon, shock_names, shock_values)

# Explicit checks
@testset "Compare shock-set irfs to the states they affect 1-for-1" begin
    @test states[endo[:z_t], 1, exo[:z_sh]] ≈ 1.
    @test states[endo[:y_t], 1, exo[:g_sh]] ≈ 1.
    @test states[endo[:g_t], 1, exo[:g_sh]] ≈ 1.
end

exp_states_shockset, exp_obs_shockset, exp_pseudo_shockset =
    jldopen("$path/../reference/impulse_responses_out.jld2", "r") do file
        read(file, "exp_states_shockset"), read(file, "exp_obs_shockset"), read(file, "exp_pseudo_shockset")
    end

@testset "Checking shockset irfs" begin
    @test states ≈ exp_states_shockset
    @test obs    ≈ exp_obs_shockset
    @test pseudo ≈ exp_pseudo_shockset
end

# State
shock_name = :z_sh
state_name   = :y_t
state_value  = 1.
states, obs, pseudo = impulse_responses(m, system, horizon, shock_name, state_name, state_value)

# Explicit checks
@testset "Check state irf initial period is exactly the specified state value" begin
    @test states[endo[:y_t], 1, exo[:z_sh]] ≈ 1.
end

exp_states_shockstates, exp_obs_shockstates, exp_pseudo_shockstates =
    jldopen("$path/../reference/impulse_responses_out.jld2", "r") do file
        read(file, "exp_states_shockstate"), read(file, "exp_obs_shockstate"), read(file, "exp_pseudo_shockstate")
    end

@testset "Checking shockstate irfs" begin
    @test states ≈ exp_states_shockstates
    @test obs    ≈ exp_obs_shockstates
    @test pseudo ≈ exp_pseudo_shockstates
end

# Observable
shock_name = :z_sh
obs_name   = :obs_gdp
obs_value  = 1.
states, obs, pseudo = impulse_responses(m, system, horizon, shock_name, obs_name, obs_value)

# Explicit checks
@testset "Check obs irf initial period is exactly the specified obs value" begin
    @test obs[obsvs[:obs_gdp], 1, exo[:z_sh]] ≈ 1.
end

exp_states_shockobs, exp_obs_shockobs, exp_pseudo_shockobs =
    jldopen("$path/../reference/impulse_responses_out.jld2", "r") do file
        read(file, "exp_states_shockobs"), read(file, "exp_obs_shockobs"), read(file, "exp_pseudo_shockobs")
    end

@testset "Checking shockobs irfs" begin
    @test states ≈ exp_states_shockobs
    @test obs    ≈ exp_obs_shockobs
    @test pseudo ≈ exp_pseudo_shockobs
end

# Test impulse response methods for:
# (1) Choosing a set of shocks and their shock values
# (2) Choosing a shock, a state, and the value the state should hit in period 1 based on
# an impulse from that shock at time 1.
# (3) Same as (2) but with an observable.

# Setup
endo   = m.endogenous_states
obsvs  = m.observables
exo    = m.exogenous_shocks
horizon = 10

# Testing both for exact output match and for certain values to match to be
# explicit about the expected outcome of the impulse response

# Shock
shock_names  = [:z_sh, :g_sh]
shock_values = [1., 1.]
states, obs, pseudo = impulse_responses(m, system, horizon, shock_names, shock_values)

# Explicit checks
@testset "Compare shock-set irfs to the states they affect 1-for-1" begin
    @test states[endo[:z_t], 1, exo[:z_sh]] ≈ 1.
    @test states[endo[:y_t], 1, exo[:g_sh]] ≈ 1.
    @test states[endo[:g_t], 1, exo[:g_sh]] ≈ 1.
end

exp_states_shockset, exp_obs_shockset, exp_pseudo_shockset =
    jldopen("$path/../reference/impulse_responses_out.jld2", "r") do file
        read(file, "exp_states_shockset"), read(file, "exp_obs_shockset"), read(file, "exp_pseudo_shockset")
    end

@testset "Checking shockset irfs" begin
    @test states ≈ exp_states_shockset
    @test obs    ≈ exp_obs_shockset
    @test pseudo ≈ exp_pseudo_shockset
end

# State
shock_name = :z_sh
state_name   = :y_t
state_value  = 1.
states, obs, pseudo = impulse_responses(m, system, horizon, shock_name, state_name, state_value)

# Explicit checks
@testset "Check state irf initial period is exactly the specified state value" begin
    @test states[endo[:y_t], 1, exo[:z_sh]] ≈ 1.
end

exp_states_shockstates, exp_obs_shockstates, exp_pseudo_shockstates =
    jldopen("$path/../reference/impulse_responses_out.jld2", "r") do file
        read(file, "exp_states_shockstate"), read(file, "exp_obs_shockstate"), read(file, "exp_pseudo_shockstate")
    end

@testset "Checking shockstate irfs" begin
    @test states ≈ exp_states_shockstates
    @test obs    ≈ exp_obs_shockstates
    @test pseudo ≈ exp_pseudo_shockstates
end

# Observable
shock_name = :z_sh
obs_name   = :obs_gdp
obs_value  = 1.
states, obs, pseudo = impulse_responses(m, system, horizon, shock_name, obs_name, obs_value)

# Explicit checks
@testset "Check obs irf initial period is exactly the specified obs value" begin
    @test obs[obsvs[:obs_gdp], 1, exo[:z_sh]] ≈ 1.
end

exp_states_shockobs, exp_obs_shockobs, exp_pseudo_shockobs =
    jldopen("$path/../reference/impulse_responses_out.jld2", "r") do file
        read(file, "exp_states_shockobs"), read(file, "exp_obs_shockobs"), read(file, "exp_pseudo_shockobs")
    end

@testset "Checking shockobs irfs" begin
    @test states ≈ exp_states_shockobs
    @test obs    ≈ exp_obs_shockobs
    @test pseudo ≈ exp_pseudo_shockobs
end

# Test impulse response method for computing
# a short-run Cholesky-identified shock
obs_shock = zeros(n_observables(m))
obs_shock[1] = 1.
states_chol, obs_chol, pseudo_chol, struct_shock =
    impulse_responses(system, horizon, Matrix{Float64}(I, n_observables(m), n_observables(m)),
                      obs_shock, flip_shocks = false, get_shocks = true)
states_chol1, obs_chol1, pseudo_chol1 =
    impulse_responses(system, horizon, Matrix{Float64}(I, n_observables(m), n_observables(m)),
                      obs_shock, flip_shocks = false, get_shocks = false)
states_chol2, obs_chol2, pseudo_chol2 =
    impulse_responses(system, horizon, Matrix{Float64}(I, n_observables(m), n_observables(m)),
                      flip_shocks = true, get_shocks = false)
states_chol3, obs_chol3, pseudo_chol3, struct_shock3 =
    impulse_responses(system, horizon, Matrix{Float64}(I, n_observables(m), n_observables(m)),
                      obs_shock, flip_shocks = false, get_shocks = true,
                      restriction = :cholesky)

exp_states_chol, exp_obs_chol, exp_pseudo_chol, exp_struct_shock =
    JLD2.jldopen("$path/../reference/impulse_responses_out.jld2", "r") do file
        read(file, "exp_states_chol"), read(file, "exp_obs_chol"), read(file, "exp_pseudo_chol"),
        read(file, "exp_struct_shock")
    end

@testset "Compare irfs to expected output for a short-run Cholesky-identified shock" begin
    @test @test_matrix_approx_eq states_chol1 states_chol
    @test @test_matrix_approx_eq obs_chol1 obs_chol
    @test @test_matrix_approx_eq pseudo_chol1 pseudo_chol
    @test @test_matrix_approx_eq states_chol2 -states_chol
    @test @test_matrix_approx_eq obs_chol2 -obs_chol
    @test @test_matrix_approx_eq pseudo_chol2  -pseudo_chol
    @test @test_matrix_approx_eq states_chol  exp_states_chol
    @test @test_matrix_approx_eq obs_chol  exp_obs_chol
    @test @test_matrix_approx_eq pseudo_chol  exp_pseudo_chol
    @test @test_matrix_approx_eq struct_shock exp_struct_shock
    @test @test_matrix_approx_eq states_chol3  exp_states_chol
    @test @test_matrix_approx_eq obs_chol3  exp_obs_chol
    @test @test_matrix_approx_eq pseudo_chol3  exp_pseudo_chol
    @test @test_matrix_approx_eq struct_shock3 exp_struct_shock
    @test @test_matrix_approx_eq -system[:RRR]*struct_shock states_chol[:, 1]
    @test @test_matrix_approx_eq -system[:ZZ]*system[:RRR]*struct_shock obs_chol[:, 1]
end

# Test impulse response method for computing
# a long-run Cholesky-identified shock
if VERSION >= v"1.1"
    # A PosDefException is triggered in Julia 1.0
    # but it is not triggered in Julia 1.1. The matrix
    # that is problematic is constructed from the following code
    # nstates = size(system[:TTT], 1)
    # obs_std = system[:ZZ] * inv(Matrix{Float64}(I, nstates, nstates) - system[:TTT]) * system[:RRR] * sqrt.(system[:QQ])
    # obs_cov = obs_std * obs_std'
    # cholesky(obs_cov) # this triggers the exception in 1.0
    # cholesky((obs_cov + obs_cov') ./ 2) # this also does too
    obs_shock = zeros(n_observables(m))
    obs_shock[1] = 1.
    states_chol, obs_chol, pseudo_chol, struct_shock =
        impulse_responses(system, horizon, Matrix{Float64}(I, n_observables(m), n_observables(m)),
                          obs_shock, flip_shocks = false, get_shocks = true,
                          restriction = :long_run)
    states_chol1, obs_chol1, pseudo_chol1 =
        impulse_responses(system, horizon, Matrix{Float64}(I, n_observables(m), n_observables(m)),
                          obs_shock, flip_shocks = false, get_shocks = false,
                          restriction = :long_run)
    states_chol2, obs_chol2, pseudo_chol2 =
        impulse_responses(system, horizon, Matrix{Float64}(I, n_observables(m), n_observables(m)),
                          flip_shocks = true, get_shocks = false, restriction = :long_run)
    states_chol3, obs_chol3, pseudo_chol3, struct_shock3 =
        impulse_responses(system, horizon, Matrix{Float64}(I, n_observables(m), n_observables(m)),
                          obs_shock, flip_shocks = false, get_shocks = true,
                          restriction = :cholesky_long_run)
    states_chol4, obs_chol4, pseudo_chol4, struct_shock4 =
        impulse_responses(system, horizon, Matrix{Float64}(I, n_observables(m), n_observables(m)),
                          obs_shock, flip_shocks = false, get_shocks = true,
                          restriction = :choleskyLR)

    exp_states_chol, exp_obs_chol, exp_pseudo_chol, exp_struct_shock =
        JLD2.jldopen("$path/../reference/impulse_responses_out.jld2", "r") do file
            read(file, "exp_states_chollr"), read(file, "exp_obs_chollr"), read(file, "exp_pseudo_chollr"),
            read(file, "exp_struct_shocklr")
        end

    @testset "Compare irfs to expected output for a long-run Cholesky-identified shock" begin
        @test @test_matrix_approx_eq states_chol1 states_chol
        @test @test_matrix_approx_eq obs_chol1 obs_chol
        @test @test_matrix_approx_eq pseudo_chol1  pseudo_chol
        @test @test_matrix_approx_eq states_chol2 -states_chol
        @test @test_matrix_approx_eq obs_chol2 -obs_chol
        @test @test_matrix_approx_eq pseudo_chol2  -pseudo_chol
        @test @test_matrix_approx_eq states_chol  exp_states_chol
        @test @test_matrix_approx_eq obs_chol  exp_obs_chol
        @test @test_matrix_approx_eq pseudo_chol  exp_pseudo_chol
        @test @test_matrix_approx_eq struct_shock exp_struct_shock
        @test @test_matrix_approx_eq states_chol3  exp_states_chol
        @test @test_matrix_approx_eq obs_chol3  exp_obs_chol
        @test @test_matrix_approx_eq pseudo_chol3  exp_pseudo_chol
        @test @test_matrix_approx_eq struct_shock3 exp_struct_shock
        @test @test_matrix_approx_eq states_chol4  exp_states_chol
        @test @test_matrix_approx_eq obs_chol4  exp_obs_chol
        @test @test_matrix_approx_eq pseudo_chol4  exp_pseudo_chol
        @test @test_matrix_approx_eq struct_shock4 exp_struct_shock
        @test @test_matrix_approx_eq -system[:RRR]*struct_shock states_chol[:, 1]
        @test @test_matrix_approx_eq -system[:ZZ]*system[:RRR]*struct_shock obs_chol[:, 1]
    end
end

# Test impulse response method for computing
# a shock to maximizes business-cycle variance
states_bc1, obs_bc1, pseudo_bc1 = impulse_responses(system, 4, (2 * π / 32, 2 * π / 6), 2)
states_bc2, _ = impulse_responses(system, 4, (2 * π / 32, 2 * π / 6), 2,
                                                    flip_shocks = true)
exp_states_bc, exp_obs_bc, exp_pseudo_bc =
    JLD2.jldopen("$path/../reference/impulse_responses_out.jld2", "r") do file
        read(file, "exp_states_bc"), read(file, "exp_obs_bc"), read(file, "exp_pseudo_bc")
    end

@testset "Compare irfs to expected output for shock maximizing business cycle variance of an observable" begin
    @test @test_matrix_approx_eq states_bc1 -states_bc2
    @test @test_matrix_approx_eq states_bc1  exp_states_bc
    @test @test_matrix_approx_eq obs_bc1  exp_obs_bc
    @test @test_matrix_approx_eq pseudo_bc1  exp_pseudo_bc
end

@testset "Check impulse responses to pegged interest rate" begin
    global m = SmetsWouters(
        custom_settings = Dict{Symbol, Setting}(:n_anticipated_shocks => Setting(:n_anticipated_shocks, 8)))
    global system = compute_system(m)
    for hor=1:8
        DSGE.impulse_responses_peg(m, system, horizon, H = hor, peg = :all_periods, real_rate = false)
        DSGE.impulse_responses_peg(m, system, horizon, H = hor, peg = :all_periods, real_rate = true)
        DSGE.impulse_responses_peg(m, system, horizon, H = hor, peg = :some_periods, real_rate = false)
        DSGE.impulse_responses_peg(m, system, horizon, H = hor, peg = :some_periods, real_rate = true)
    end
end

nothing


text:
using Base.Test
using ThreePhotons

# Determination

params,state = rotation_search(Dict("reference_pdb_path"=>"$(ENV["DETERMINATION_DATA"])/structures/crambin.pdb","stepsizefactor"=>1.01, "initial_stepsize" => Float64(pi/100), "L"=>L, "K3_range" =>1:K, "N"=>N, "histograms"=>"/tmp/correlations.dat", "optimizer"=>rotate_all_at_once, "initial_temperature_factor"=>1.0e-6, "measure"=>"Bayes", "temperature_decay"=>0.96, "LMAX"=>LMAX, "K2_range"=>1:KMAX, "qmax"=>float(KMAX), "lambda"=>0.0, "include_negativity"=>false))


text:
function (loss::Loss)(par, model)
    F = zero(eltype(par))
    for i = 1:length(loss.functions)
        F += loss.functions[i](par, model)
        # all functions have to have those arguments??
    end
    return F
end

# function (loss::Loss)(par, model, E, G)
#     if E != nothing
#         F = zero(eltype(model.imply.imp_cov))
#         for i = 1:length(loss.functions)
#             F += loss.functions[i](par, model, E, G)
#         end
#         return F
#     end
#     if G != nothing
#         G .= zero(eltype(G))
#         for i = 1:length(loss.functions)
#             loss.functions[i](par, model, E, G)
#         end
#         for i = 1:length(loss.functions)
#             G .+= loss.functions[i].grad
#         end
#
#     end
# end

function (loss::Loss)(par, model, E, G)

    if E != nothing

        for i = 1:length(loss.functions)
            loss.functions[i](par, model, E, G)
        end

        objective = zero(eltype(model.imply.imp_cov))
        for i = 1:length(loss.functions)
            objective += loss.functions[i].objective[1]
        end

        return objective
    end

    if G != nothing

        for i = 1:length(loss.functions)
            loss.functions[i](par, model, E, G)
        end

        G .= zero(eltype(G))
        for i = 1:length(loss.functions)
            G .+= loss.functions[i].grad
        end

    end
end

# function (loss::Loss)(par, model, E, G)
#
#     #common computations with 2 arguments
#     #store computations in fields
#     for i = 1:length(loss.functions)
#         loss.functions[i](par, model)
#     end
#
#     for i = 1:length(loss.functions)
#         loss.functions[i](par, model, E, G)
#     end
#
#     if E != nothing
#         objective = zero(eltype(model.imply.imp_cov))
#         for i = 1:length(loss.functions)
#             objective += loss.functions[i].objective[1]
#         end
#         return objective
#     end
#
#     if G != nothing
#
#         G .= zero(eltype(G))
#         for i = 1:length(loss.functions)
#             G .+= loss.functions[i].grad
#         end
#
#     end
# end


text:
module FunctionSyntax

"f_1"
function f_1 end

"""
f_2
"""
function f_2 end

end


text:
using Test
using UnitTestDesign
#### Set-based

### tuples_in_trials
tt_cases = [
    [[[1,2,3], [1,2,1], [2,3,4]], 2, 8],
    [[[1,1,1], [1,1,1], [1,1,1]], 2, 3],
    [[[1,1,1], [1,1,1], [1,1,1]], 3, 1],
]
for tt_case in tt_cases
    res = UnitTestDesign.tuples_in_trials(tt_case[1], tt_case[2])
    total_covered = sum([length(aset) for aset in values(res)])
    @test total_covered == tt_case[3]
end

sc = UnitTestDesign.SetCoverage([2,3,2], 2)
@test UnitTestDesign.remaining(sc) == 0
@test eltype(sc) == Int64
UnitTestDesign.build_all_combinations!(sc, 2)
@test UnitTestDesign.remaining(sc) == 16
@test UnitTestDesign.add_coverage!(sc, [1, 1, 1]) == 3
@test UnitTestDesign.remaining(sc) == 13
@test UnitTestDesign.add_coverage!(sc, [1, 2, 1]) == 2
@test UnitTestDesign.remaining(sc) == 11

UnitTestDesign.build_all_combinations!(sc, 2)
multiple = [2 2; 2 3; 1 2]
@test UnitTestDesign.add_coverage!(sc, multiple) == 6
@test UnitTestDesign.remaining(sc) == 10

# Add two types of coverage
sc = UnitTestDesign.SetCoverage([2,3,2], 2)
UnitTestDesign.build_all_combinations!(sc, 2)
UnitTestDesign.build_all_combinations!(sc, 3)
@test UnitTestDesign.remaining(sc) == 28


text:
using Pkg.TOML, JSON

red(msg)    = println("\033[0;31m$msg\033[0m")
green(msg)  = println("\033[0;32m$msg\033[0m")
yellow(msg) = println("\033[0;33m$msg\033[0m")

function reg_url(pkg)
  return "https://raw.githubusercontent.com/JuliaRegistries/General/master/$(pkg[1])/$pkg/Deps.toml"
end

function jso_pkgs()
  json = JSON.parse(read(`curl https://api.github.com/orgs/JuliaSmoothOptimizers/repos`, String))
  pkgs = map(x -> x["name"], json)
  pkgs = filter(x -> x[end-1:end] == "jl", pkgs)
  pkgs = map(x -> x[1:end-3], pkgs)
  return sort(pkgs)
end

function jso_deps()
  pkgs = jso_pkgs()
  pkgs = setdiff(pkgs, ["qr_mumps", "QPSReader", "SuiteSparseMatrixCollection", "QuadraticModels", "JSOSuite"])
  #pkgs = ["LinearOperators", "Krylov", "NLPModels", "NLPModelsIpopt", "CUTEst"]
  deps = dependency_tree(pkgs, limitto=pkgs)
  return deps
end

dependency_tree(pkg :: String; limitto=[]) = dependency_tree([pkg], limitto=limitto)

function dependency_tree(pkgs :: Array{<: String}; limitto=[])
  todo = copy(pkgs)
  deps = Dict()
  versions = Dict()
  ignored = String[]
  fname = tempname()
  while length(todo) > 0
    pkg = pop!(todo)
    green("Processing $pkg")
    url = reg_url(pkg)
    try
      download(url, fname)
      download(replace(url, "Deps" => "Versions"), fname * "_v")
    catch ex
      yellow("$pkg not found in Registry, possible a stdlib package")
      deps[pkg] = []
      push!(ignored, pkg)
      continue
    end
    toml = TOML.parsefile(fname)
    dep_pkgs = String[]
    for p in keys.(values(toml))
      append!(dep_pkgs, p)
    end
    if length(limitto) > 0
      dep_pkgs = dep_pkgs ∩ limitto
    end
    deps[pkg] = dep_pkgs
    for p in dep_pkgs
      if p == "Homebrew"
        println(dep_pkgs)
        println(limitto)
      end
      if !(p in keys(deps)) && !(p in ignored)
        push!(todo, p)
      end
    end
    versions[pkg] = sort(VersionNumber.(keys(TOML.parsefile(fname * "_v"))))[end]
  end
  return deps, versions
end

function depth_computation(deps :: Dict)
  pkgs = sort(collect(keys(deps)))
  n = length(pkgs)
  depth = zeros(Int, n)

  function depth_update(i, v)
    if depth[i] ≥ v
      depth[i] = v
      J = indexin(deps[pkgs[i]], pkgs)
      for j in J
        depth_update(j, v - 1)
      end
      return true
    else
      return false
    end
  end

  # Find feasible depth
  for i = 1:n
    depth_update(i, 0)
  end
  # Increase for positivity
  depth .+= -minimum(depth)

  done = false
  while !done
    done = true
    for i = 1:n
      J = indexin(deps[pkgs[i]], pkgs)
      d = if length(J) == 0
        0
      else
        max(0, maximum(depth[J]) + 1)
      end
      if depth[i] > d
        depth[i] = d
        done = false
      end
    end
  end

  return depth
end

function weight_computation(deps, depth)
  pkgs = sort(collect(keys(deps)))
  n = length(pkgs)
  weight = ones(Int, n)
  for d = maximum(depth):-1:0
    for i = findall(depth .== d)
      J = indexin(deps[pkgs[i]], pkgs)
      weight[J] .+= weight[i]
    end
  end
  return weight
end

function tikz_draw(deps, depth, weight, versions)
  pkgs = sort(collect(keys(deps)))
  #S = sortperm(weight)
  #depth = depth[S]
  #pkgs  = pkgs[S]
  open("tree.tex", "w") do f
    println(f, "\\begin{tikzpicture}")
    M = maximum(depth)
    for d = 0:M
      J = findall(depth .== d)
      for (c,j) = enumerate(J)
        pkg = replace(pkgs[j], "_" => "\\_")
        v = versions[pkgs[j]]
        san = lowercase(replace(pkgs[j], "_" => ""))
        x = 6d
        y = -10 * (c - 1) / (length(J) - 1)
        println(f, "\\node[draw,fill=white] ($san) at ($x,$y) {$pkg - $v};")
      end
    end
    println(f, "\\begin{scope}[on background layer]")
    for (i,pkg) in enumerate(pkgs)
      J = indexin(deps[pkg], pkgs)
      san_dest = lowercase(replace(pkg, "_" => ""))
      for j in J
        san_orig = lowercase(replace(pkgs[j], "_" => ""))
        Δd = abs(depth[j] - depth[i])
        println(f, "\\draw[->,gray] ($san_orig) to [in=180,out=0,looseness=0.3] ($san_dest);")
      end
    end
    println(f, "\\end{scope}")
    println(f, "\\end{tikzpicture}")
  end
end

#deps = dependency_tree("NLPModels")
#pkgs = jso_pkgs()
deps, versions = jso_deps()
depth = depth_computation(deps)
weight = weight_computation(deps, depth)
tikz_draw(deps, depth, weight, versions)

max_stars_repo_name: abelsiqueira/dependency-tree-julia

max_issues_repo_name: abelsiqueira/dependency-tree-julia

max_forks_repo_name: abelsiqueira/dependency-tree-julia


text:
using Netpbm, ColorTypes, FixedPointNumbers, IndirectArrays, ImageCore, FileIO
using Test

@testset "IO" begin
    workdir = joinpath(tempdir(), "Images")
    isdir(workdir) && rm(workdir, recursive=true)
    mkdir(workdir)

    @testset "Gray pgm" begin
        af = rand(2, 3)
        for T in (N0f8, N4f12, N0f16,
                  Gray{N0f8}, Gray{N4f12}, Gray{N0f16})
            ac = convert(Array{T}, af)
            fn = File(format"PGMBinary", joinpath(workdir, "3by2.pgm"))
            Netpbm.save(fn, ac)
            b = Netpbm.load(fn)
            @test b == ac
        end
        a8 = convert(Array{N0f8}, af)
        for T in (Float32, Float64, Gray{Float32}, Gray{Float64})
            ac = convert(Array{T}, af)
            fn = File(format"PGMBinary", joinpath(workdir, "3by2.pgm"))
            Netpbm.save(fn, ac)
            b = Netpbm.load(fn)
            @test b == a8
        end

        # Issue #11
        @test_nowarn load("test_header_space.pgm")
    end

    @testset "Color ppm" begin
        af = rand(RGB{Float64}, 2, 3)
        for T in (RGB{N0f8}, RGB{N4f12}, RGB{N0f16})
            ac = convert(Array{T}, af)
            fn = File(format"PPMBinary", joinpath(workdir, "3by2.ppm"))
            Netpbm.save(fn, ac)
            b = Netpbm.load(fn)
            @test b == ac
        end
        a8 = convert(Array{RGB{N0f8}}, af)
        for T in (RGB{Float32}, RGB{Float64}, HSV{Float64})
            ac = convert(Array{T}, af)
            fn = File(format"PPMBinary", joinpath(workdir, "3by2.ppm"))
            Netpbm.save(fn, ac)
            b = Netpbm.load(fn)
            @test b == a8
        end
    end

    @testset "Colormap" begin
        datafloat = reshape(range(0.5, stop=1.5, length=6), 2, 3)
        dataint = map(x->round(UInt8, x), 254*(datafloat .- 0.5) .+ 1) # ranges 1 to 255
        # build our colormap
        b = RGB(0,0,1)
        w = RGB(1,1,1)
        r = RGB(1,0,0)
        cmaprgb = Array{RGB{N0f8}}(undef, 255)
        f = range(0, stop=1, length=128)
        cmaprgb[1:128] = [(1-x)*b + x*w for x in f]
        cmaprgb[129:end] = [(1-x)*w + x*r for x in f[2:end]]
        img = IndirectArray(dataint, cmaprgb)
        fn = File(format"PPMBinary", joinpath(workdir,"cmap.ppm"))
        Netpbm.save(fn, img)
        imgr = Netpbm.load(fn)
        @test imgr == img
        cmaprgb = Array{RGB}(undef, 255) # poorly-typed cmap, Images issue #336
        @test !isconcretetype(eltype(cmaprgb))
        cmaprgb[1:128] = RGB{N0f16}[(1-x)*b + x*w for x in f]
        cmaprgb[129:end] = RGB{N4f12}[(1-x)*w + x*r for x in f[2:end]]
        img = IndirectArray(dataint, cmaprgb)
        @test_throws ErrorException Netpbm.save(fn, img) # widens to unsupported type
        cmaprgb[129:end] = RGB{N0f8}[(1-x)*w + x*r for x in f[2:end]]
        img = IndirectArray(dataint, cmaprgb)
        Netpbm.save(fn, img)
        imgr = Netpbm.load(fn)
        @test imgr == img
    end

    # Images issue #256
    @testset "Clamping" begin
        A = rand(2,3)
        A[1,1] = -0.4
        fn = File(format"PGMBinary", joinpath(workdir, "2by3.pgm"))
        @test_throws InexactError Netpbm.save(fn, A)
        Netpbm.save(fn, A, mapf=clamp01nan)
        B = Netpbm.load(fn)
        A[1,1] = 0
        @test B == Gray{N0f8}.(A)
    end
end

nothing


text:
using Nextion
using Test

@testset "number" begin

    nexSerial = NexSerial("/dev/ttyUSB0")

    nexNumber = NexNumber(nexSerial, Name("n0"), cid=ComponentID(1))

    send(nexSerial, "page pg_num")

    #=
    sleep(1)

    nexNumber.value = 1

    sleep(1)

    nexNumber.value = 2
    nexNumber.backcolor = Colour.RED
    nexNumber.forecolor = Colour.WHITE

    sleep(1)

    @test nexNumber.value == 2

    sleep(1)

    nexNumber.value = 3
    nexNumber.backcolor = Colour.WHITE
    nexNumber.forecolor = Colour.RED

    sleep(1)
    =#

    #=
    n = typemax(Int32)
    nexNumber.value = n
    sleep(1)
    @test nexNumber.value == n
    =#

    #=
    n = typemin(Int32)
    nexNumber.value = n
    sleep(1)
    @test nexNumber.value == n
    =#

    n = -2
    nexNumber.value = n
    sleep(1)
    @test nexNumber.value == n

    close(nexSerial)

    @test true

end


text:
"""
    load_model(file_name::String)::MetabolicModel

Generic function for loading models that chooses a specific loader function
from the `file_name` extension, or throws an error.

Currently, these model types are supported:

- SBML models (`*.xml`, loaded with [`load_sbml_model`](@ref))
- JSON models (`*.json`, loaded with [`load_json_model`](@ref))
- MATLAB models (`*.mat`, loaded with [`load_mat_model`](@ref))
"""
function load_model(file_name::String)::MetabolicModel
    if endswith(file_name, ".json")
        return load_json_model(file_name)
    elseif endswith(file_name, ".xml")
        return load_sbml_model(file_name)
    elseif endswith(file_name, ".mat")
        return load_mat_model(file_name)
    else
        throw(DomainError(file_name, "Unknown file extension"))
    end
end

"""
    load_model(type::Type{T}, file_name::String)::T where T

Helper function tht loads the model using [`load_model`](@ref) and return it
converted to `type`.

# Example:

    load_model(CoreModel, "mySBMLModel.xml")
"""
function load_model(type::Type{T}, file_name::String)::T where {T<:MetabolicModel}
    convert(type, load_model(file_name))
end

"""
    save_model(model::MetabolicModel, file_name::String)

Generic function for saving models that chooses a specific writer function
from the `file_name` extension, or throws an error.

Currently, these model types are supported:

- JSON models (`*.json`, loaded with [`save_json_model`](@ref))
- MATLAB models (`*.mat`, loaded with [`save_mat_model`](@ref))
"""
function save_model(model::MetabolicModel, file_name::String)
    if endswith(file_name, ".json")
        return save_json_model(model, file_name)
    elseif endswith(file_name, ".mat")
        return save_mat_model(model, file_name)
    else
        throw(DomainError(file_name, "Unknown file extension"))
    end
end


text:
using Oceananigans.Grids

#####
##### Centered fourth-order advection scheme
#####

struct CenteredFourthOrder <: AbstractAdvectionScheme end

const C4 = CenteredFourthOrder

const centered_second_order = CenteredSecondOrder()

@inline ℑ³xᶠᵃᵃ(i, j, k, grid, u) = @inbounds u[i, j, k] - δxᶠᵃᵃ(i, j, k, grid, δxᶜᵃᵃ, u) / 6
@inline ℑ³xᶜᵃᵃ(i, j, k, grid, c) = @inbounds c[i, j, k] - δxᶜᵃᵃ(i, j, k, grid, δxᶠᵃᵃ, c) / 6

@inline ℑ³yᵃᶠᵃ(i, j, k, grid, v) = @inbounds v[i, j, k] - δyᵃᶠᵃ(i, j, k, grid, δyᵃᶜᵃ, v) / 6
@inline ℑ³yᵃᶜᵃ(i, j, k, grid, c) = @inbounds c[i, j, k] - δyᵃᶜᵃ(i, j, k, grid, δyᵃᶠᵃ, c) / 6

@inline ℑ³zᵃᵃᶠ(i, j, k, grid, w) = @inbounds w[i, j, k] - δzᵃᵃᶠ(i, j, k, grid, δzᵃᵃᶜ, w) / 6
@inline ℑ³zᵃᵃᶜ(i, j, k, grid, c) = @inbounds c[i, j, k] - δzᵃᵃᶜ(i, j, k, grid, δzᵃᵃᶠ, c) / 6

# Momentum

@inline fourth_order_momentum_flux_uu(i, j, k, grid, u)    = ℑxᶜᵃᵃ(i, j, k, grid, Ax_ψᵃᵃᶠ, u) * ℑxᶜᵃᵃ(i, j, k, grid, ℑ³xᶠᵃᵃ, u)
@inline fourth_order_momentum_flux_uv(i, j, k, grid, u, v) = ℑxᶠᵃᵃ(i, j, k, grid, Ay_ψᵃᵃᶠ, v) * ℑyᵃᶠᵃ(i, j, k, grid, ℑ³yᵃᶜᵃ, u)
@inline fourth_order_momentum_flux_uw(i, j, k, grid, u, w) = ℑxᶠᵃᵃ(i, j, k, grid, Az_ψᵃᵃᵃ, w) * ℑzᵃᵃᶠ(i, j, k, grid, ℑ³zᵃᵃᶜ, u)

@inline fourth_order_momentum_flux_vu(i, j, k, grid, u, v) = ℑyᵃᶠᵃ(i, j, k, grid, Ax_ψᵃᵃᶠ, u) * ℑxᶠᵃᵃ(i, j, k, grid, ℑ³xᶜᵃᵃ, v)
@inline fourth_order_momentum_flux_vv(i, j, k, grid, v)    = ℑyᵃᶜᵃ(i, j, k, grid, Ay_ψᵃᵃᶠ, v) * ℑyᵃᶜᵃ(i, j, k, grid, ℑ³yᵃᶠᵃ, v)
@inline fourth_order_momentum_flux_vw(i, j, k, grid, v, w) = ℑyᵃᶠᵃ(i, j, k, grid, Az_ψᵃᵃᵃ, w) * ℑzᵃᵃᶠ(i, j, k, grid, ℑ³zᵃᵃᶜ, v)

@inline fourth_order_momentum_flux_wu(i, j, k, grid, u, w) = ℑzᵃᵃᶠ(i, j, k, grid, Ax_ψᵃᵃᶠ, u) * ℑxᶠᵃᵃ(i, j, k, grid, ℑ³xᶜᵃᵃ, w)
@inline fourth_order_momentum_flux_wv(i, j, k, grid, v, w) = ℑzᵃᵃᶠ(i, j, k, grid, Ay_ψᵃᵃᶠ, v) * ℑyᵃᶠᵃ(i, j, k, grid, ℑ³yᵃᶜᵃ, w)
@inline fourth_order_momentum_flux_ww(i, j, k, grid, w)    = ℑzᵃᵃᶜ(i, j, k, grid, Az_ψᵃᵃᵃ, w) * ℑzᵃᵃᶜ(i, j, k, grid, ℑ³zᵃᵃᶠ, w)

# Periodic directions!

@inline momentum_flux_uu(i, j, k, grid, ::C4, u)    = fourth_order_momentum_flux_uu(i, j, k, grid, u)
@inline momentum_flux_uv(i, j, k, grid, ::C4, u, v) = fourth_order_momentum_flux_uv(i, j, k, grid, u, v)
@inline momentum_flux_uw(i, j, k, grid, ::C4, u, w) = fourth_order_momentum_flux_uw(i, j, k, grid, u, w)

@inline momentum_flux_vu(i, j, k, grid, ::C4, u, v) = fourth_order_momentum_flux_vu(i, j, k, grid, u, v)
@inline momentum_flux_vv(i, j, k, grid, ::C4, v)    = fourth_order_momentum_flux_vv(i, j, k, grid, v)
@inline momentum_flux_vw(i, j, k, grid, ::C4, v, w) = fourth_order_momentum_flux_vw(i, j, k, grid, v, w)

@inline momentum_flux_wu(i, j, k, grid, ::C4, u, w) = fourth_order_momentum_flux_wu(i, j, k, grid, u, w)
@inline momentum_flux_wv(i, j, k, grid, ::C4, v, w) = fourth_order_momentum_flux_wv(i, j, k, grid, v, w)
@inline momentum_flux_ww(i, j, k, grid, ::C4, w)    = fourth_order_momentum_flux_ww(i, j, k, grid, w)

# Bounded directions

@inline function momentum_flux_uu(i, j, k, grid::AbstractGrid{FT, <:Bounded}, ::C4, u) where FT
    if i > 1 && i < grid.Nx
        return fourth_order_momentum_flux_uu(i, j, k, grid, u)
    else
        return momentum_flux_uu(i, j, k, grid, centered_second_order, u)
    end
end

@inline function momentum_flux_uv(i, j, k, grid::AbstractGrid{FT, TX, <:Bounded}, ::C4, u, v) where {FT, TX}
    if j > 1 && j < grid.Ny
        return fourth_order_momentum_flux_uv(i, j, k, grid, u, v)
    else
        return momentum_flux_uv(i, j, k, grid, centered_second_order, u, v)
    end
end

@inline function momentum_flux_uw(i, j, k, grid::AbstractGrid{FT, TX, TY, <:Bounded}, ::C4, u, w) where {FT, TX, TY}
    if k > 1 && k < grid.Nz
        return fourth_order_momentum_flux_uw(i, j, k, grid, u, w)
    else
        return momentum_flux_uw(i, j, k, grid, centered_second_order, u, w)
    end
end

@inline function momentum_flux_vu(i, j, k, grid::AbstractGrid{FT, <:Bounded}, ::C4, u, v) where {FT}
    if i > 1 && i < grid.Nx
        return fourth_order_momentum_flux_vu(i, j, k, grid, u, v)
    else
        return momentum_flux_vu(i, j, k, grid, centered_second_order, u, v)
    end
end

@inline function momentum_flux_vv(i, j, k, grid::AbstractGrid{FT, TX, <:Bounded}, ::C4, v) where {FT, TX}
    if j > 1 && j < grid.Ny
        return fourth_order_momentum_flux_vv(i, j, k, grid, v)
    else
        return momentum_flux_vv(i, j, k, grid, centered_second_order, v)
    end
end

@inline function momentum_flux_vw(i, j, k, grid::AbstractGrid{FT, TX, TY, <:Bounded}, ::C4, v, w) where {FT, TX, TY}
    if k > 1 && k < grid.Nz
        return fourth_order_momentum_flux_vw(i, j, k, grid, v, w)
    else
        return momentum_flux_vw(i, j, k, grid, centered_second_order, v, w)
    end
end

@inline function momentum_flux_wu(i, j, k, grid::AbstractGrid{FT, <:Bounded}, ::C4, u, w) where FT
    if i > 1 && i < grid.Nx
        return fourth_order_momentum_flux_wu(i, j, k, grid, u, w)
    else
        return momentum_flux_wu(i, j, k, grid, centered_second_order, u, w)
    end
end

@inline function momentum_flux_wv(i, j, k, grid::AbstractGrid{FT, TX, <:Bounded}, ::C4, v, w) where {FT, TX}
    if j > 1 && j < grid.Ny
        return fourth_order_momentum_flux_wv(i, j, k, grid, v, w)
    else
        return momentum_flux_wv(i, j, k, grid, centered_second_order, v, w)
    end
end

@inline function momentum_flux_ww(i, j, k, grid::AbstractGrid{FT, TX, TY, <:Bounded}, ::C4, w) where {FT, TX, TY}
    if k > 1 && i < grid.Nz
        return fourth_order_momentum_flux_ww(i, j, k, grid, w)
    else
        return momentum_flux_ww(i, j, k, grid, centered_second_order, w)
    end
end

# Tracers

@inline fourth_order_advective_tracer_flux_x(i, j, k, grid, u, c) = Ax_ψᵃᵃᶠ(i, j, k, grid, u) * ℑxᶠᵃᵃ(i, j, k, grid, ℑ³xᶜᵃᵃ, c)
@inline fourth_order_advective_tracer_flux_y(i, j, k, grid, v, c) = Ay_ψᵃᵃᶠ(i, j, k, grid, v) * ℑyᵃᶠᵃ(i, j, k, grid, ℑ³yᵃᶜᵃ, c)
@inline fourth_order_advective_tracer_flux_z(i, j, k, grid, w, c) = Az_ψᵃᵃᵃ(i, j, k, grid, w) * ℑzᵃᵃᶠ(i, j, k, grid, ℑ³zᵃᵃᶜ, c)

@inline advective_tracer_flux_x(i, j, k, grid, ::C4, u, c) = fourth_order_advective_tracer_flux_x(i, j, k, grid, u, c)
@inline advective_tracer_flux_y(i, j, k, grid, ::C4, v, c) = fourth_order_advective_tracer_flux_y(i, j, k, grid, v, c)
@inline advective_tracer_flux_z(i, j, k, grid, ::C4, w, c) = fourth_order_advective_tracer_flux_z(i, j, k, grid, w, c)

@inline function advective_tracer_flux_x(i, j, k, grid::AbstractGrid{FT, <:Bounded}, ::C4, u, c) where FT
    if i > 1 && i < grid.Nx
        return fourth_order_advective_tracer_flux_x(i, j, k, grid, u, c)
    else
        return advective_tracer_flux_x(i, j, k, grid, centered_second_order, u, c)
    end
end

@inline function advective_tracer_flux_y(i, j, k, grid::AbstractGrid{FT, TX, <:Bounded}, ::C4, v, c) where {FT, TX}
    if j > 1 && i < grid.Ny
        return fourth_order_advective_tracer_flux_y(i, j, k, grid, v, c)
    else
        return advective_tracer_flux_y(i, j, k, grid, centered_second_order, v, c)
    end
end

@inline function advective_tracer_flux_z(i, j, k, grid::AbstractGrid{FT, TX, TY, <:Bounded}, ::C4, w, c) where {FT, TX, TY}
    if k > 1 && i < grid.Nz
        return fourth_order_advective_tracer_flux_z(i, j, k, grid, w, c)
    else
        return advective_tracer_flux_z(i, j, k, grid, centered_second_order, w, c)
    end
end


text:
using Documenter, DocumenterMarkdown, Turing, AdvancedHMC, Bijectors
using LibGit2: clone

# Include the utility functions.
include("make-utils.jl")

# Paths.
source_path = joinpath(@__DIR__, "src")
build_path = joinpath(@__DIR__, "_docs")
tutorial_path = joinpath(@__DIR__, "_tutorials")

with_clean_docs(source_path, build_path) do source, build
    makedocs(
        sitename = "Turing.jl",
        source = source,
        build = build,
        format = Markdown(),
        checkdocs = :all
    )
end

# You can skip this part if you are on a metered
# connection by calling `julia make.jl no-tutorials`
in("no-tutorials", ARGS) || copy_tutorial(tutorial_path)

# set default baseurl for the master branch
baseurl = "/dev"
# set baseurl for version tag when current head is tagged
vtag = get(ENV, "TRAVIS_TAG", "") # if on Travis
if isempty(vtag) # so we are using GitHub Actions
    vtag = replace(get(ENV, "GITHUB_REF", ""), r"/?refs/tags/?" => "")
end

version_match = match(r"^(v\d+\.\d+\.\d+)$", vtag)
if !isnothing(version_match)
    baseurl = "/" * version_match[1]
end

jekyll_build = joinpath(@__DIR__, "jekyll-build")
with_baseurl(() -> run(`$jekyll_build`), baseurl)

# deploy
devurl = "dev"
repo = "github.com:TuringLang/Turing.jl.git"
deploydocs(
    target = "_site",
    repo = repo,
    branch = "gh-pages",
    devbranch = "master",
    devurl = devurl,
    versions = ["stable" => "v^", "v#.#", devurl => devurl]
)


text:
@doc raw"""
Tableau of a Partitioned Runge-Kutta method
```math
\begin{aligned}
Q_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} a_{ij} \, v(t_{n} + c_j \Delta t, Q_{n,j}, P_{n,j}) , &
q_{n+1} &= q_{n} + h \sum \limits_{i=1}^{s} b_{i}  \, v(t_{n} + c_j \Delta t, Q_{n,i}, P_{n,i}) , \\
P_{n,i} &= p_{n} + h  \sum \limits_{i=1}^{s} \bar{a}_{ij} \, f(t_{n} + c_j \Delta t, Q_{n,j}, P_{n,j}) , &
p_{n+1} &= p_{n} + h \sum \limits_{i=1}^{s} \bar{b}_{i}   \, f(t_{n} + c_j \Delta t, Q_{n,i}, P_{n,i}) .
\end{aligned}
```

Parameters:
 * `T`: datatype of coefficient arrays

Fields:
 * `name`: symbolic name of the tableau
 * `o`: order of the method
 * `s`: number of stages
 * `q`: Tableau for `q`
 * `p`: Tableau for `p`

The actual tableaus are stored in `q` and `p`:
 * `a`: coefficients $a_{ij}$ with $ 1 \le i,j \le s$
 * `b`: weights $b_{i}$  with $ 1 \le i \le s$
 * `c`: nodes $c_{i}$  with $ 1 \le i \le s$

Constructors:
```julia
PartitionedTableau{T}(name, o, s, q, p)
PartitionedTableau{T}(name, q, p)
PartitionedTableau(name::Symbol, q::Tableau, p::Tableau)
PartitionedTableau(name::Symbol, q::Tableau)
```
"""
struct PartitionedTableau{T}
    name::Symbol
    o::Int
    s::Int

    q::Tableau{T}
    p::Tableau{T}

    function PartitionedTableau{T}(name, o, s, q, p) where {T}
        @assert s == q.s == p.s
        new(name, o, s, q, p)
    end

    function PartitionedTableau{T}(name, q, p) where {T}
        new(name, min(q.o, p.o), q.s, q, p)
    end
end

PartitionedTableau(name::Symbol, q::Tableau{T}, p::Tableau{T}) where {T} = PartitionedTableau{T}(name, q, p)
PartitionedTableau(name::Symbol, q::Tableau) = PartitionedTableau(name, q, q)

Base.hash(tab::PartitionedTableau, h::UInt) = hash(tab.o, hash(tab.s, hash(tab.q, hash(tab.p, hash(:PartitionedTableau, h)))))

Base.:(==)(tab1::PartitionedTableau, tab2::PartitionedTableau) = (
                                            tab1.o == tab2.o
                                         && tab1.s == tab2.s
                                         && tab1.q == tab2.q
                                         && tab1.p == tab2.p)

Base.isapprox(tab1::PartitionedTableau, tab2::PartitionedTableau; kwargs...) = (
                                            tab1.o == tab2.o
                                         && tab1.s == tab2.s
                                         && isapprox(tab1.q, tab2.q; kwargs...)
                                         && isapprox(tab1.p, tab2.p; kwargs...))

Base.isequal(tab1::PartitionedTableau{T1}, tab2::PartitionedTableau{T2}) where {T1,T2} = (tab1 == tab2 && T1 == T2 && tab1.name == tab2.name)

Base.eltype(::PartitionedTableau{T}) where {T} = T

name(tab::PartitionedTableau) = tab.name
order(tab::PartitionedTableau) = tab.o
nstages(tab::PartitionedTableau) = tab.s
eachstage(tab::PartitionedTableau) = 1:tab.s

isexplicit(tab::PartitionedTableau) = istril(tab.q.a) && istril(tab.p.a) && all([tab.q.a[i,i] == 0 || tab.p.a[i,i] == 0 for i in 1:tab.s]) && (tab.q.c[1] == 0 || tab.p.c[1] == 0)
isimplicit(tab::PartitionedTableau) = !isexplicit(tab)
isdiagnonallyimplicit(tab::PartitionedTableau) = isimplicit(tab) && tab.s != 1 && istril(tab.q.a) && istril(tab.p.a)
isfullyimplicit(tab::PartitionedTableau) = !isexplicit(tab) && !(isdiagnonallyimplicit(tab))