text:
module VarCoeffLassoTest

const  NUMBER_REPEAT = 1

using CoordinateDescent
using Test
using ProximalBase
using Random
using LinearAlgebra
using SparseArrays

Random.seed!(1)

@testset "kernels" begin
  x = 0.3
  y = 0.4
  k = GaussianKernel(1.)
  @test evaluate(k, x, y) ≈ exp(-0.01)

  function evaluate!(g::GaussianKernel{T}, buffer::Array{T}, x::Array{T}, y::Float64) where {T <: AbstractFloat}
    size(buffer) == size(x) || throw(DimensionMismatch())

    @inbounds @simd for i=1:length(x)
       buffer[i] = evaluate(g, x[i], y)
    end
    buffer
  end

  x = rand(100, 1000)
  y = rand()
  k = GaussianKernel(0.5)

  o1 = zero(x)
  o2 = zero(x)
  o1 = evaluate.(Ref(k), x, Ref(y))
  evaluate!(k, o2, x, y)

  @test o1 ≈ o2
end

@testset "expand_X" begin
  X = reshape(collect(1.:6.), 2, 3)
  z = [0.2, 0.4]
  z0 = 0.3

  tX = similar(X)
  @test CoordinateDescent._expand_X!(tX, X, z, z0, 0) == X

  tX = zeros(2, 6)
  tX1 = zeros(2, 6)
  Q = [1. -0.1; 1. 0.1]
  for i=1:2
    tX1[i, :] = kron(X[i,:], Q[i, :])
  end
  @test CoordinateDescent._expand_X!(tX, X, z, z0, 1) ≈ tX1

  tX = zeros(2, 9)
  tX1 = zeros(2, 9)
  Q = [1. -0.1 0.01; 1. 0.1 0.01]
  for i=1:2
    tX1[i, :] = kron(X[i,:], Q[i, :])
  end
  @test CoordinateDescent._expand_X!(tX, X, z, z0, 2) ≈ tX1
end

@testset "expand_X multiplications" begin

  p = 10
  X = randn(100, p)
  Y = randn(100)
  z = rand(100)
  w = zeros(100)
  k = GaussianKernel(0.2)
  w = evaluate.(Ref(k), z, Ref(0.5))

  for degree=0:2
    cp = p*(degree+1)
    eX = zeros(100, cp)
    ewX = zeros(100, cp)

    CoordinateDescent._expand_X!(eX, X, z, 0.5, degree)
    @test CoordinateDescent._expand_wX!(ewX, w, X, z, 0.5, degree) ≈ Diagonal(w)*eX

    Xt_w_Y = zeros(cp)
    @test CoordinateDescent._expand_Xt_w_Y!(Xt_w_Y, w, X, z, Y, 0.5, degree) ≈ eX' * Diagonal(w) * Y

    Xt_w_X = zeros(cp, cp)
    @test CoordinateDescent._expand_Xt_w_X!(Xt_w_X, w, X, z, 0.5, degree) ≈ (eX'*Diagonal(w))*eX
  end
end
#
# @testset "locpoly" begin
#
#   n, p = 500, 2
#   Y, X, Z, betaMat = genData(n, p)
#   zgrid = collect(0.01:0.2:0.99)
#
#   gk = GaussianKernel(0.4)
#
#   w = zeros(n)
#
#   for degree=0:2
#     @test locpoly(X, Z, Y, zgrid, degree, gk) ≈ locpoly_alt(X, Z, Y, zgrid, degree, gk)
#
#     z0 = 0.5
#     cp = p*(degree+1)
#
#     eX = zeros(n, cp)
#     _expand_X!(eX, X, Z, z0, degree)
#
#     @. w = evaluate(gk, Z, z0)
#
#     @test locpoly(X, Z, Y, z0, degree, gk) ≈ (eX' * diagm(w) * eX)\(eX' * diagm(w) * Y)
#   end
# end
#
#
#
# @testset "locpolyl1" begin
#
#   n, s, p = 500, 10, 50
#   gk = GaussianKernel(0.1)
#   zgrid = collect(0.01:0.1:0.99)
#   opt = CDOptions(;randomize=false)
#
#   for i=1:NUMBER_REPEAT
#     Y, X, Z, betaT = genData(n, s)
#     X = [X zeros(n, p-s)]
#
#     λ0 = rand() / 10
#
#     for degree=0:2
#       o1 = locpolyl1(X,Z,Y,zgrid,degree,gk,λ0, opt)
#       o2 = locpolyl1_alt(X,Z,Y,zgrid,degree,gk,λ0, opt)
#
#       @test maximum( maximum(abs.(o1[i] - o2[i])) for i=1:length(zgrid) ) ≈ 0. atol=1e-4
#     end
#   end
#
# end

end


text:
# This is file provides the ability for the user to override the default prior.

# This code is `include()`'ed inside of `/scripts/venus.jl`, so it will **literarly replace** the
# function lnprior inside `src/model.jl`. Be careful to make sure it accepts the same arguments!

using DiskJockey.model
import DiskJockey.model.lnprior_base
import DiskJockey.constants.ModelException

# for example, a replacement for the standard model might look like
function lnprior(pars::ParametersStandard, dpc_mu::Float64, dpc_sig::Float64, grid::Grid)
    lnp = lnprior_base(pars, dpc_mu, dpc_sig)

    if pars.mu_RA > -20
        println("Stupid RA")
        throw(ModelException("User prior violated."))
    end

    r_out = grid.Rs[end]/AU # [AU]
    # A somewhat arbitrary cutoff regarding the gridsize to prevent the disk from being too large
    # to fit on the model grid.
    if ((3 * pars.r_c) > r_out) || (pars.q > 0.75)
        throw(ModelException("User prior violated."))
    end

    return lnp
end

ext: jl

max_stars_repo_name: elnjensen/DiskJockey

max_issues_repo_name: elnjensen/DiskJockey

max_forks_repo_name: elnjensen/DiskJockey


text:
module SddLibrary

@static if Sys.isunix()
    @static if Sys.islinux()
        const LIBSDD = "$(@__DIR__)/../deps/sdd-2.0/lib/Linux/libsdd"
    elseif Sys.isapple()
        const LIBSDD = "$(@__DIR__)/../deps/sdd-2.0/lib/Darwin/libsdd"
    else
        LoadError("sddapi.jl", 0, "Sdd library only available on Linux and Darwin")
    end
else
    LoadError("sddapi.jl", 0, "Sdd library only available on Linux and Darwin")
end

const SddSize = Csize_t
const SddNodeSize = Cuint
const SddRefCount = Cuint
const SddModelCount = Culonglong
const SddWmc = Cdouble
const SddLiteral = Clong

const SddID = SddSize

const BoolOp = Cushort
const CONJOIN = convert(BoolOp, 0)
const DISJOIN = convert(BoolOp, 1)

struct VTree_c end
struct SddNode_c end
struct SddManager_c end
struct WmcManager_c end

# SDD MANAGER FUNCTIONS
function sdd_manager_new(vtree::Ptr{VTree_c})::Ptr{SddManager_c}
    return ccall((:sdd_manager_new, LIBSDD), Ptr{SddManager_c}, (Ptr{VTree_c},), vtree)
end
function sdd_manager_create(var_count::SddLiteral, auto_gc_and_minimize::Cint)::Ptr{SddManager_c}
    return ccall((:sdd_manager_create, LIBSDD),Ptr{SddManager_c}, (SddLiteral,Cint), var_count, auto_gc_and_minimize)
end
# function sdd_manager_new(size::SddSize, nodes::Array{Ptr{SddNode_c}}, from_manager::Ptr{SddManager_c})::Ptr{SddManager_c}
#     return ccall((:sdd_manager_copy, LIBSDD), Ptr{SddManager_c}, (SddSize,Array{Ptr{SddNode_c}},Ptr{SddManager_c}), size, nodes, from_manager)
# end
function sdd_manager_free(manager::Ptr{SddManager_c})
    ccall((:sdd_manager_free, LIBSDD), Cvoid, (Ptr{SddManager_c},), manager)
end
function sdd_manager_print(manager::Ptr{SddManager_c})
    ccall((:sdd_manager_print, LIBSDD), Cvoid, (Ptr{SddManager_c},), manager)
end
function sdd_manager_auto_gc_and_minimize_on(manager::Ptr{SddManager_c})
    ccall((:sdd_manager_auto_gc_and_minimize_on, LIBSDD), Cvoid, (Ptr{SddManager_c},), manager)
end
function sdd_manager_auto_gc_and_minimize_off(manager::Ptr{SddManager_c})
    ccall((:sdd_manager_auto_gc_and_minimize_off, LIBSDD), Cvoid, (Ptr{SddManager_c},), manager)
end
function sdd_manager_is_auto_gc_and_minimize_on(manager::Ptr{SddManager_c})::Cint
    return ccall((:sdd_manager_is_auto_gc_and_minimize_on, LIBSDD), Cint, (Ptr{SddManager_c},), manager)
end
# TODO void sdd_manager_set_minimize_function(SddVtreeSearchFunc func, SddManager* manager);
function sdd_manager_unset_minimize_function(manager::Ptr{SddManager_c})
    ccall((:sdd_manager_unset_minimize_function, LIBSDD), Cvoid, (Ptr{SddManager_c},), manager)
end
# function sdd_manager_options(manager::Ptr{SddManager_c})
#     ccall((:sdd_manager_options, LIBSDD), Ptr{Cvoid}, (Ptr{SddManager_c},), manager)
# end
# void sdd_manager_set_options(void* options, SddManager* manager);
function sdd_manager_is_var_used(var::SddLiteral, manager::Ptr{SddManager_c})::Cint
    return ccall((:sdd_manager_is_var_used, LIBSDD), Cint, (SddLiteral, Ptr{SddManager_c}), var, manager)
end
function sdd_manager_vtree_of_var(var::SddLiteral, manager::Ptr{SddManager_c})::Ptr{VTree_c}
    return ccall((:sdd_manager_vtree_of_var, LIBSDD), Ptr{VTree_c}, (SddLiteral, Ptr{SddManager_c}), var, manager)
end
function sdd_manager_lca_of_literals(count::SddLiteral, literals::Array{SddLiteral,1}, manager::Ptr{SddManager_c})::Ptr{VTree_c}
    return ccall((:sdd_manager_lca_of_literals, LIBSDD), Ptr{VTree_c}, (Int32, Ptr{SddLiteral}, Ptr{SddManager_c}), count, literals, manager)
end
function sdd_manager_var_count(manager::Ptr{SddManager_c})::SddLiteral
    return ccall((:sdd_manager_var_count, LIBSDD), SddLiteral, (Ptr{SddManager_c},), manager)
end
# TODO void sdd_manager_var_order(SddLiteral* var_order, SddManager *manager);
function sdd_manager_add_var_before_first(manager::Ptr{SddManager_c})
    ccall((:sdd_manager_add_var_before_first, LIBSDD), Cvoid, (Ptr{SddManager_c},), manager)
end
function sdd_manager_add_var_after_last(manager::Ptr{SddManager_c})
    ccall((:sdd_manager_add_var_after_last, LIBSDD), Cvoid, (Ptr{SddManager_c},), manager)
end
function sdd_manager_add_var_before(manager::Ptr{SddManager_c})
    ccall((:sdd_manager_add_var_before, LIBSDD), Cvoid, (Ptr{SddManager_c},), manager)
end
function sdd_manager_add_var_after(manager::Ptr{SddManager_c})
    ccall((:sdd_manager_add_var_after, LIBSDD), Cvoid, (Ptr{SddManager_c},), manager)
end

# TERMINAL SDDS
function sdd_manager_true(manager::Ptr{SddManager_c})::Ptr{SddNode_c}
    return ccall((:sdd_manager_true, LIBSDD), Ptr{SddNode_c}, (Ptr{SddManager_c},), manager)
end
function sdd_manager_false(manager::Ptr{SddManager_c})::Ptr{SddNode_c}
    return ccall((:sdd_manager_false, LIBSDD), Ptr{SddNode_c}, (Ptr{SddManager_c},), manager)
end
function sdd_manager_literal(literal::SddLiteral, manager::Ptr{SddManager_c})::Ptr{SddNode_c}
    return ccall((:sdd_manager_literal, LIBSDD), Ptr{SddNode_c}, (SddLiteral, Ptr{SddManager_c}), literal, manager)
end

# SDD QUERIES AND TRANSFORMATIONS
function sdd_apply(node1::Ptr{SddNode_c}, node2::Ptr{SddNode_c}, op::BoolOp ,manager::Ptr{SddManager_c})::Ptr{SddNode_c}
    return ccall((:sdd_apply, LIBSDD), Ptr{SddNode_c}, (Ptr{SddNode_c}, Ptr{SddNode_c}, BoolOp, Ptr{SddManager_c}), node1, node2, op, manager)
end
function sdd_conjoin(node1::Ptr{SddNode_c}, node2::Ptr{SddNode_c}, manager::Ptr{SddManager_c})::Ptr{SddNode_c}
    return ccall((:sdd_conjoin, LIBSDD), Ptr{SddNode_c}, (Ptr{SddNode_c}, Ptr{SddNode_c}, Ptr{SddManager_c}), node1, node2, manager)
end
function sdd_disjoin(node1::Ptr{SddNode_c}, node2::Ptr{SddNode_c}, manager::Ptr{SddManager_c})::Ptr{SddNode_c}
    return ccall((:sdd_disjoin, LIBSDD), Ptr{SddNode_c}, (Ptr{SddNode_c}, Ptr{SddNode_c}, Ptr{SddManager_c}), node1, node2, manager)
end
function sdd_negate(node::Ptr{SddNode_c}, manager::Ptr{SddManager_c})::Ptr{SddNode_c}
    return ccall((:sdd_negate, LIBSDD), Ptr{SddNode_c}, (Ptr{SddNode_c}, Ptr{SddManager_c}), node, manager)
end
function sdd_condition(lit::SddLiteral, node::Ptr{SddNode_c}, manager::Ptr{SddManager_c})::Ptr{SddNode_c}
    return ccall((:sdd_condition, LIBSDD), Ptr{SddNode_c}, (SddLiteral, Ptr{SddNode_c}, Ptr{SddManager_c}), lit, node, manager)
end
function sdd_exists(lit::SddLiteral, node::Ptr{SddNode_c}, manager::Ptr{SddManager_c})::Ptr{SddNode_c}
    return ccall((:sdd_exists, LIBSDD), Ptr{SddNode_c}, (SddLiteral, Ptr{SddNode_c}, Ptr{SddManager_c}), lit, node, manager)
end
function sdd_exists_multiple(exists_map::Array{Cint,1}, node::Ptr{SddNode_c}, manager::Ptr{SddManager_c})::Ptr{SddNode_c}
    return ccall((:sdd_exists_multiple, LIBSDD), Ptr{SddNode_c}, (Array{Cint,1}, Ptr{SddNode_c}, Ptr{SddManager_c}), exists_map, node, manager)
end
function sdd_exists_multiple_static(exists_map::Array{Cint,1}, node::Ptr{SddNode_c}, manager::Ptr{SddManager_c})::Ptr{SddNode_c}
    return ccall((:sdd_exists_multiple_static, LIBSDD), Ptr{SddNode_c}, (Array{Cint,1}, Ptr{SddNode_c}, Ptr{SddManager_c}), exists_map, node, manager)
end
function sdd_forall(lit::SddLiteral, node::Ptr{SddNode_c}, manager::Ptr{SddManager_c})::Ptr{SddNode_c}
    return ccall((:sdd_forall, LIBSDD), Ptr{SddNode_c}, (SddLiteral, Ptr{SddNode_c}, Ptr{SddManager_c}), lit, node, manager)
end
function sdd_minimize_cardinality(node::Ptr{SddNode_c}, manager::Ptr{SddManager_c})::Ptr{SddNode_c}
    return ccall((:sdd_minimize_cardinality, LIBSDD), Ptr{SddNode_c}, (Ptr{SddNode_c}, Ptr{SddManager_c}), node, manager)
end
function sdd_global_minimize_cardinality(node::Ptr{SddNode_c}, manager::Ptr{SddManager_c})::Ptr{SddNode_c}
    return ccall((:sdd_global_minimize_cardinality, LIBSDD), Ptr{SddNode_c}, (Ptr{SddNode_c}, Ptr{SddManager_c}), node, manager)
end
function sdd_minimum_cardinality(node::Ptr{SddNode_c})::SddLiteral
    return ccall((:sdd_minimum_cardinality, LIBSDD), SddLiteral, (Ptr{SddNode_c}, ), node)
end
function sdd_model_count(node::Ptr{SddNode_c}, manager::Ptr{SddManager_c})::SddModelCount
    return ccall((:sdd_model_count, LIBSDD), SddModelCount, (Ptr{SddNode_c}, Ptr{SddManager_c}), node, manager)
end
function sdd_global_model_count(node::Ptr{SddNode_c}, manager::Ptr{SddManager_c})::SddModelCount
    return ccall((:sdd_global_model_count, LIBSDD), SddModelCount, (Ptr{SddNode_c}, Ptr{SddManager_c}), node, manager)
end

# // SDD NAVIGATION
function sdd_node_is_true(node::Ptr{SddNode_c})::Int32
    return ccall((:sdd_node_is_true, LIBSDD), Cint, (Ptr{SddNode_c}, ), node)
end
function sdd_node_is_false(node::Ptr{SddNode_c})::Int32
    return ccall((:sdd_node_is_false, LIBSDD), Cint, (Ptr{SddNode_c}, ), node)
end
function sdd_node_is_literal(node::Ptr{SddNode_c})::Int32
    return ccall((:sdd_node_is_literal, LIBSDD), Cint, (Ptr{SddNode_c}, ), node)
end
function sdd_node_is_decision(node::Ptr{SddNode_c})::Int32
    return ccall((:sdd_node_is_decision, LIBSDD), Cint, (Ptr{SddNode_c}, ), node)
end
function sdd_node_size(node::Ptr{SddNode_c})::SddNodeSize
    return ccall((:sdd_node_size, LIBSDD), SddNodeSize, (Ptr{SddNode_c}, ), node)
end
function sdd_node_literal(node::Ptr{SddNode_c})::SddLiteral
    return ccall((:sdd_node_literal, LIBSDD), SddLiteral, (Ptr{SddNode_c}, ), node)
end
function sdd_node_elements(node::Ptr{SddNode_c})::Ptr{Ptr{SddNode_c}}
    return ccall((:sdd_node_elements, LIBSDD), Ptr{Ptr{SddNode_c}}, (Ptr{SddNode_c}, ), node)
end
function sdd_node_set_bit(bit::Int32, node::Ptr{SddNode_c})
    ccall((:sdd_node_set_bit, LIBSDD), Cvoid, (Cint, Ptr{SddNode_c}), bit, node)
end
function sdd_node_bit(node::Ptr{SddNode_c})::Int32
    return ccall((:sdd_node_bit, LIBSDD), Int32, (Ptr{SddNode_c}, ), node)
end

# # SDD FUNCTIONS
#
# SDD FILE I/O
function sdd_read(filename::Ptr{UInt8}, manager::Ptr{SddManager_c})::Ptr{SddNode_c}
    return ccall((:sdd_read, LIBSDD), Ptr{SddNode_c}, (Ptr{UInt8}, Ptr{SddManager_c}), filename, manager)
end
function sdd_save(filename::Ptr{UInt8}, node::Ptr{SddNode_c})
    ccall((:sdd_save, LIBSDD), Cvoid, (Ptr{UInt8}, Ptr{SddNode_c}), filename, node)
end
function sdd_save_as_dot(filename::Ptr{UInt8}, node::Ptr{SddNode_c})
    ccall((:sdd_save_as_dot, LIBSDD), Cvoid, (Ptr{UInt8}, Ptr{SddNode_c}), filename, node)
end
function sdd_shared_save_as_dot(filename::Ptr{UInt8}, manager::Ptr{SddManager_c})
    ccall((:sdd_shared_save_as_dot, LIBSDD), Cvoid, (Ptr{UInt8}, Ptr{SddManager_c}), filename, manager)
end

# // SDD SIZE AND NODE COUNT
# //SDD
function sdd_count(node::Ptr{SddNode_c})::SddSize
    return ccall((:sdd_count, LIBSDD), SddSize, (Ptr{SddNode_c},), node)
end
function sdd_size(node::Ptr{SddNode_c})::SddSize
    return ccall((:sdd_size, LIBSDD), SddSize, (Ptr{SddNode_c},), node)
end
# TODO SddSize sdd_shared_size(SddNode** nodes, SddSize count);
# //SDD OF MANAGER
manager_size_fnames_c = [
    "sdd_manager_size", "sdd_manager_live_size", "sdd_manager_dead_size",
    "sdd_manager_count", "sdd_manager_live_count", "sdd_manager_dead_count"
]
for fnc in manager_size_fnames_c
    @eval begin
        function $(Symbol(fnc))(manager::Ptr{SddManager_c})::SddSize
            return ccall(($(:($(fnc))), LIBSDD), SddSize, (Ptr{SddManager_c},), manager)
        end
    end
end

# SDD SIZE OF VTREE
vtree_size_fnames_c = [
    "sdd_vtree_size", "sdd_vtree_live_size", "sdd_vtree_dead_size",
    "sdd_vtree_size_at", "sdd_vtree_live_size_at", "sdd_vtree_dead_size_at",
    "sdd_vtree_size_above", "sdd_vtree_live_size_above", "sdd_vtree_dead_size_above",
    "sdd_vtree_count", "sdd_vtree_live_count", "sdd_vtree_dead_count",
    "sdd_vtree_count_at", "sdd_vtree_live_count_at", "sdd_vtree_dead_count_at",
    "sdd_vtree_count_above", "sdd_vtree_live_count_above", "sdd_vtree_dead_count_above"
]
for fnc in vtree_size_fnames_c
    @eval begin
        function $(Symbol(fnc))(vtree::Ptr{VTree_c})::SddSize
            return ccall(($(:($(fnc))), LIBSDD), SddSize, (Ptr{VTree_c},), vtree)
        end
    end
end

# CREATING VTREES
function sdd_vtree_new(var_count::SddLiteral, vtree_type::Ptr{UInt8})::Ptr{VTree_c}
    return ccall((:sdd_vtree_new, LIBSDD), Ptr{VTree_c}, (SddLiteral, Ptr{UInt8}), var_count, vtree_type)
end
function sdd_vtree_new_with_var_order(var_count::SddLiteral, var_order::Array{SddLiteral,1}, vtree_type::Ptr{UInt8})::Ptr{VTree_c}
    return ccall((:sdd_vtree_new_with_var_order, LIBSDD), Ptr{VTree_c}, (SddLiteral, Ptr{SddLiteral}, Ptr{UInt8}), var_count, var_order, vtree_type)
end
function sdd_vtree_new_X_constrained(var_count::SddLiteral, is_X_var::Array{SddLiteral,1}, vtree_type::Ptr{UInt8})::Ptr{VTree_c}
    return ccall((:sdd_vtree_new_X_constrained, LIBSDD), Ptr{VTree_c}, (SddLiteral, Ptr{SddLiteral}, Ptr{UInt8}), var_count, is_X_var, vtree_type)
end
function sdd_vtree_free(vtree::Ptr{VTree_c})
    ccall((:sdd_vtree_free, LIBSDD), Cvoid, (Ptr{VTree_c},), vtree)
end

# VTREE FILE I/O
function sdd_vtree_read(filename::Ptr{UInt8})::Ptr{VTree_c}
    return ccall((:sdd_vtree_read, LIBSDD), Ptr{VTree_c}, (Ptr{UInt8},), filename)
end
function sdd_vtree_save(filename::Ptr{UInt8}, vtree::Ptr{VTree_c})
    ccall((:sdd_vtree_save, LIBSDD), Cvoid, (Ptr{UInt8}, Ptr{VTree_c}), filename, vtree)
end
function sdd_vtree_save_as_dot(filename::Ptr{UInt8}, vtree::Ptr{VTree_c})
    ccall((:sdd_vtree_save_as_dot, LIBSDD), Cvoid, (Ptr{UInt8}, Ptr{VTree_c}), filename, vtree)
end

# // SDD MANAGER VTREE
function sdd_manager_vtree(manager::Ptr{SddManager_c})::Ptr{VTree_c}
    return ccall((:sdd_manager_vtree, LIBSDD), Ptr{VTree_c}, (Ptr{SddManager_c},), manager)
end
function sdd_manager_vtree_copy(manager::Ptr{SddManager_c})::Ptr{VTree_c}
    return ccall((:sdd_manager_vtree_copy, LIBSDD), Ptr{VTree_c}, (Ptr{SddManager_c},), manager)
end

# // VTREE NAVIGATION
function sdd_vtree_left(vtree::Ptr{VTree_c})::Ptr{VTree_c}
    return ccall((:sdd_vtree_left, LIBSDD), Ptr{VTree_c}, (Ptr{VTree_c},), vtree)
end
function sdd_vtree_right(vtree::Ptr{VTree_c})::Ptr{VTree_c}
    return ccall((:sdd_vtree_right, LIBSDD), Ptr{VTree_c}, (Ptr{VTree_c},), vtree)
end
function sdd_vtree_parent(vtree::Ptr{VTree_c})::Ptr{VTree_c}
    return ccall((:sdd_vtree_parent, LIBSDD), Ptr{VTree_c}, (Ptr{VTree_c},), vtree)
end

# VTREE FUNCTIONS
function sdd_vtree_is_leaf(vtree::Ptr{VTree_c})::Cint
    return ccall((:sdd_vtree_is_leaf, LIBSDD), Cint, (Ptr{VTree_c}, ), vtree)
end
function sdd_vtree_is_sub(vtree1::Ptr{VTree_c}, vtree2::Ptr{VTree_c})::Cint
    return ccall((:sdd_vtree_is_sub, LIBSDD), Cint, (Ptr{VTree_c}, Ptr{VTree_c}), vtree1, vtree2)
end
function sdd_vtree_lca(vtree1::Ptr{VTree_c}, vtree2::Ptr{VTree_c}, root::Ptr{VTree_c})::Ptr{VTree_c}
    return ccall((:sdd_vtree_lca, LIBSDD), Ptr{VTree_c}, (Ptr{VTree_c}, Ptr{VTree_c}, Ptr{VTree_c}), vtree1, vtree2, root)
end
function sdd_vtree_var_count(vtree::Ptr{VTree_c})::SddLiteral
    return ccall((:sdd_vtree_var_count, LIBSDD), SddLiteral, (Ptr{VTree_c}, ), vtree)
end
function sdd_vtree_var(vtree::Ptr{VTree_c})::SddLiteral
    return ccall((:sdd_vtree_var, LIBSDD), SddLiteral, (Ptr{VTree_c}, ), vtree)
end
function sdd_vtree_position(vtree::Ptr{VTree_c})::SddLiteral
    return ccall((:sdd_vtree_position, LIBSDD), SddLiteral, (Ptr{VTree_c}, ), vtree)
end
# Vtree** sdd_vtree_location(Vtree* vtree, SddManager* manager);

# VTREE/SDD EDIT OPERATIONS
function sdd_vtree_rotate_left(vtree::Ptr{VTree_c}, manager::Ptr{SddManager_c}, limited::Cint)::Cint
    return ccall((:sdd_vtree_rotate_left, LIBSDD), Cint, (Ptr{VTree_c},Ptr{SddManager_c}, Cint), vtree, manager, limited)
end
function sdd_vtree_rotate_right(vtree::Ptr{VTree_c}, manager::Ptr{SddManager_c}, limited::Cint)::Cint
    return ccall((:sdd_vtree_rotate_right, LIBSDD), Cint, (Ptr{VTree_c},Ptr{SddManager_c}, Cint), vtree, manager, limited)
end
function sdd_vtree_swap(vtree::Ptr{VTree_c}, manager::Ptr{SddManager_c}, limited::Cint)::Cint
    return ccall((:sdd_vtree_swap, LIBSDD), Cint, (Ptr{VTree_c}, Ptr{SddManager_c}, Cint), vtree, manager, limited)
end

# LIMITS FOR VTREE/SDD EDIT OPERATIONS
function sdd_manager_init_vtree_size_limit(vtree::Ptr{VTree_c}, manager::Ptr{SddManager_c})
    ccall((:sdd_manager_init_vtree_size_limit, LIBSDD), Cvoid, (Ptr{VTree_c}, Ptr{SddManager_c}), vtree, manager)
end
function sdd_manager_update_vtree_size_limit(manager::Ptr{SddManager_c})
    ccall((:sdd_manager_update_vtree_size_limit, LIBSDD), Cvoid, (Ptr{SddManager_c},), manager)
end

# # VTREE STATE

# GARBAGE COLLECTION
function sdd_ref_count(node::Ptr{SddNode_c})::SddRefCount
    return ccall((:sdd_ref_count, LIBSDD), SddRefCount, (Ptr{SddNode_c},), node)
end
function sdd_ref(node::Ptr{SddNode_c}, manager::Ptr{SddManager_c})::Ptr{SddNode_c}
    return ccall((:sdd_ref, LIBSDD), Ptr{SddNode_c}, (Ptr{SddNode_c},Ptr{SddManager_c}), node, manager)
end
function sdd_deref(node::Ptr{SddNode_c}, manager::Ptr{SddManager_c})::Ptr{SddNode_c}
    return ccall((:sdd_deref, LIBSDD), Ptr{SddNode_c}, (Ptr{SddNode_c},Ptr{SddManager_c}), node, manager)
end
function sdd_manager_garbage_collect(manager::Ptr{SddManager_c})
    ccall((:sdd_manager_garbage_collect, LIBSDD), Cvoid, (Ptr{SddManager_c},), manager)
end
function sdd_vtree_garbage_collect(vtree::Ptr{VTree_c}, manager::Ptr{SddManager_c})
    ccall((:sdd_vtree_garbage_collect, LIBSDD), Cvoid, (Ptr{VTree_c}, Ptr{SddManager_c}), vtree, manager)
end
function sdd_manager_garbage_collect_if(dead_node_threshold::Float32, manager::Ptr{SddManager_c})::Int32
    return ccall((:sdd_manager_garbage_collect_if, LIBSDD), Int32, (Float32, Ptr{SddManager_c}), dead_node_threshold, manager)
end
function sdd_vtree_garbage_collect_if(dead_node_threshold::Float32, vtree::Ptr{VTree_c}, manager::Ptr{SddManager_c})::Int32
    return ccall((:sdd_manager_garbage_collect_if, LIBSDD), Int32, (Float32, Ptr{VTree_c}, Ptr{SddManager_c}), dead_node_threshold, vtree, manager)
end

# MINIMIZATION
function sdd_manager_minimize(manager::Ptr{SddManager_c})
    ccall((:sdd_manager_minimize, LIBSDD), Cvoid, (Ptr{SddManager_c},), manager)
end
function sdd_vtree_minimize(vtree::Ptr{VTree_c}, manager::Ptr{SddManager_c})::Ptr{VTree_c}
    return ccall((:sdd_manager_minimize, LIBSDD), Ptr{VTree_c}, (Ptr{VTree_c}, Ptr{SddManager_c}), vtree, manager)
end
function sdd_manager_minimize_limited(manager::Ptr{SddManager_c})
    ccall((:sdd_manager_minimize_limited, LIBSDD), Cvoid, (Ptr{SddManager_c},), manager)
end
function sdd_vtree_minimize_limited(vtree::Ptr{VTree_c}, manager::Ptr{SddManager_c})::Ptr{VTree_c}
    return ccall((:sdd_vtree_minimize_limited, LIBSDD), Ptr{VTree_c}, (Ptr{VTree_c}, Ptr{SddManager_c}), vtree, manager)
end

function sdd_manager_set_vtree_search_convergence_threshold(threshold::Float32, manager::Ptr{SddManager_c})
    ccall((:sdd_manager_set_vtree_search_convergence_threshold, LIBSDD), Cvoid, (Float32, Ptr{SddManager_c}), threshold, manager)
end

function sdd_manager_set_vtree_search_time_limit(time_limit::Float32, manager::Ptr{SddManager_c})
    ccall((:sdd_manager_set_vtree_search_time_limit, LIBSDD), Cvoid, (Float32, Ptr{SddManager_c}), time_limit, manager)
end
function sdd_manager_set_vtree_fragment_time_limit(time_limit::Float32, manager::Ptr{SddManager_c})
    ccall((:sdd_manager_set_vtree_fragment_time_limit, LIBSDD), Cvoid, (Float32, Ptr{SddManager_c}), time_limit, manager)
end
function sdd_manager_set_vtree_operation_time_limit(time_limit::Float32, manager::Ptr{SddManager_c})
    ccall((:sdd_manager_set_vtree_operation_time_limit, LIBSDD), Cvoid, (Float32, Ptr{SddManager_c}), time_limit, manager)
end
function sdd_manager_set_vtree_apply_time_limit(time_limit::Float32, manager::Ptr{SddManager_c})
    ccall((:sdd_manager_set_vtree_apply_time_limit, LIBSDD), Cvoid, (Float32, Ptr{SddManager_c}), time_limit, manager)
end
function sdd_manager_set_vtree_operation_memory_limit(memory_limit::Float32, manager::Ptr{SddManager_c})
    ccall((:sdd_manager_set_vtree_operation_memory_limit, LIBSDD), Cvoid, (Float32, Ptr{SddManager_c}), memory_limit, manager)
end
function sdd_manager_set_vtree_operation_size_limit(size_limit::Float32, manager::Ptr{SddManager_c})
    ccall((:sdd_manager_set_vtree_operation_size_limit, LIBSDD), Cvoid, (Float32, Ptr{SddManager_c}), size_limit, manager)
end
function sdd_manager_set_vtree_cartesian_product_limit(size_limit::Float32, manager::Ptr{SddManager_c})
    ccall((:sdd_manager_set_vtree_cartesian_product_limit, LIBSDD), Cvoid, (Float32, Ptr{SddManager_c}), size_limit, manager)
end

# WMC
function wmc_manager_new(node::Ptr{SddNode_c}, log_mode::Cint, manager::Ptr{SddManager_c})::Ptr{WmcManager_c}
    return ccall((:wmc_manager_new, LIBSDD), Ptr{WmcManager_c}, (Ptr{SddNode_c}, Cint, Ptr{SddManager_c}), node, log_mode, manager)
end
function wmc_manager_free(wmc_manager::Ptr{WmcManager_c})
    ccall((:wmc_manager_free, LIBSDD), Cvoid, (Ptr{WmcManager_c},), wmc_manager)
end
function wmc_set_literal_weight(literal::SddLiteral, weight::SddWmc, wmc_manager::Ptr{WmcManager_c})
    ccall((:wmc_set_literal_weight, LIBSDD), Cvoid, (SddLiteral, SddWmc, Ptr{WmcManager_c}), literal, weight, wmc_manager)
end
function wmc_propagate(wmc_manager::Ptr{WmcManager_c})::SddWmc
    return ccall((:wmc_propagate, LIBSDD), SddWmc, (Ptr{WmcManager_c},), wmc_manager)
end
function wmc_zero_weight(wmc_manager::Ptr{WmcManager_c})::SddWmc
    return ccall((:wmc_zero_weight, LIBSDD), SddWmc, (Ptr{WmcManager_c},), wmc_manager)
end
function wmc_one_weight(wmc_manager::Ptr{WmcManager_c})::SddWmc
    return ccall((:wmc_one_weight, LIBSDD), SddWmc, (Ptr{WmcManager_c},), wmc_manager)
end
function wmc_literal_weight(literal::SddLiteral, wmc_manager::Ptr{WmcManager_c})::SddWmc
    return ccall((:wmc_literal_weight, LIBSDD), SddWmc, (SddLiteral, Ptr{WmcManager_c}), literal, wmc_manager)
end
function wmc_literal_derivative(literal::SddLiteral, wmc_manager::Ptr{WmcManager_c})::SddWmc
    return ccall((:wmc_literal_derivative, LIBSDD), SddWmc, (SddLiteral, Ptr{WmcManager_c}), literal, wmc_manager)
end
function wmc_literal_pr(literal::SddLiteral, wmc_manager::Ptr{WmcManager_c})::SddWmc
    return ccall((:wmc_literal_pr, LIBSDD), SddWmc, (SddLiteral, Ptr{WmcManager_c}), literal, wmc_manager)
end

end

# TO BE WRAPPED
# // SDD FUNCTIONS
# SddSize sdd_id(SddNode* node);
# int sdd_garbage_collected(SddNode* node, SddSize id);
# Vtree* sdd_vtree_of(SddNode* node);
# SddNode* sdd_copy(SddNode* node, SddManager* dest_manager);
# SddNode* sdd_rename_variables(SddNode* node, SddLiteral* variable_map, SddManager* manager);
# int* sdd_variables(SddNode* node, SddManager* manager);

# // VTREE STATE
# int sdd_vtree_bit(const Vtree* vtree);
# void sdd_vtree_set_bit(int bit, Vtree* vtree);
# void* sdd_vtree_data(Vtree* vtree);
# void sdd_vtree_set_data(void* data, Vtree* vtree);
# void* sdd_vtree_search_state(const Vtree* vtree);
# void sdd_vtree_set_search_state(void* search_state, Vtree* vtree);


text:
"""
`read_formula(pkg::Union{AbstractString,BrewPkg})`

Returns the string contents of a package's formula.
"""
function read_formula(pkg::StringOrPkg) end

function read_formula(name::AbstractString)
    return read(formula_path(name), String)
end

function read_formula(pkg::BrewPkg)
    return read(formula_path(pkg), String)
end

"""
`write_formula(name::AbstractString, formula::AbstractString)`

Write out fully-qualified formula `name` with contents `formula` to disk. Note
that writing out without a tap name is not allowed; we won't write new formulae
out to `Homebrew/core`, only to taps.
"""
function write_formula(name::AbstractString, formula::AbstractString)
    path, tap_path = formula_tap(name)

    if isempty(tap_path)
        error("Cannot write a formula out to Homebrew/core!")
    end

    # Insert the "homebrew-" that exists in all taps
    tap_path = "$(dirname(tap_path))/homebrew-$(basename(tap_path))"

    # Open this file, and write it out!
    path = joinpath(brew_prefix, "Library", "Taps", tap_path, "$path.rb")
    open(path, "w") do f
        write(f, formula)
    end

    # We done!
    return
end

"""
`delete_translated_formula(name::AbstractString; verbose::Bool=false)`

Delete a translated formula from the `staticfloat/juliatranslated` tap.
"""
function delete_translated_formula(name::AbstractString; verbose::Bool=false)
    # Throw out the tap_path part of any formula that someone passed into here.
    path, tap_path = formula_tap(name)

    # Override tap_path with our auto_tapname
    tap_path = "$(dirname(auto_tapname))/homebrew-$(basename(auto_tapname))"

    try
        del_path = joinpath(brew_prefix, "Library", "Taps", tap_path, "$path.rb")
        Base.rm(del_path)
        if verbose
            println("Deleting $(basename(del_path))...")
        end
    catch
    end
end

"""
`delete_all_translated_formulae(;verbose::Bool=false)`

Delete all translated formulae from the `staticfloat/juliatranslated` tap. This
is useful for debugging misbehaving formulae during translation.
"""
function delete_all_translated_formulae(;verbose::Bool=false)
    for f in readdir(auto_tappath)
        if endswith(f,".rb")
            if verbose
                println("Deleting $f...")
            end
            Base.rm(joinpath(auto_tappath,f))
        end
    end
end

"""
`translate_formula(pkg::Union{AbstractString,BrewPkg}; verbose::Bool=false)`

Given a formula `name`, return the fully-qualified name of a translated formula
if it is translatable.  Translation copies a `Homebrew/core` formula to
`$auto_tapname`, adding appropriate `cellar :any` and `root_url` lines to
any bottle stanzas.  This allows us to transparently install non-cellar-any
formulae from `Homebrew/core`.

This function is fairly strict, bailing out at every possible opportunity, and
returning the original name.  If a formula is non-translatable, it's possible
it needs manual intervention, check out the $tapname tap for examples.
"""
function translate_formula(pkg::StringOrPkg; verbose::Bool=false) end

function translate_formula(name::AbstractString; verbose::Bool=false)
    if verbose
        println("translation: beginning for $name")
    end

    path, tap_path = formula_tap(name)

    # We maintain a list of formulae that we WILL NOT translate. As an example,
    # 'xz' is automatically added as a dependency by Homebrew when a formula
    # has a `.tar.xz` source tarball.  This cannot be redirected to the
    # `staticfloat/juliatranslated` tap, causing diamonds of death where we have
    # the potential for both `xz` and `staticfloat/juliatranslated/xz` to
    # be in the dependency tree for a formula.
    translation_blacklist = ["xz"]
    if basename(name) in translation_blacklist
        if verbose
            println("translation: bailing because $name is in the translation blacklist")
        end
        return name
    end

    # Did we ask for any old name, or did we explicitly request a tap?
    if isempty(tap_path)
        # Bail if there is an overriding formula in our manually-curated tap
        if isfile(joinpath(tappath, "$(name).rb"))
            if verbose
                println("translation: using $tapname/$name for the source of our translation")
            end
            tap_path = tapname
            name = "$(tap_path)/$(path)"
        end
    else
        # If we explicitly asked for a tap, then make sure it's here!
        tap(tap_path)
    end

    # Bail if we have no bottles
    if !has_bottle(name)
        if verbose
            println("translation: bailing because $name has no bottles")
        end
        return name
    end

    # Delete any older translated formula if they exist
    auto_path = joinpath(auto_tappath, "$(path).rb")
    override_name = joinpath(auto_tapname, path)
    if isfile(auto_path)
        src_path = formula_path(name)
        if stat(auto_path).mtime < stat(src_path).mtime
            if verbose
                println("translation: deleting stale translation for $name")
            end
            delete_translated_formula(override_name; verbose=verbose)
        else
            if verbose
                println("translation: bailing becuase $name is already available in tap $auto_tapname")
            end
            return joinpath(auto_tapname, path)
        end
    end

    # Read formula source in, and also get a JSON representation of it
    obj = json(name)
    formula = read_formula(name)

    # Find bottle section. We allow 1 to 8 lines of code in a bottle stanza:
    # a root_url, a prefix, a cellar, a revision, and four OSX version bottles.

    ex = r"(?:\r\n|\r|\n)\s*bottle\s+do(?:\r\n|\r|\n)(?:[^\n]*(?:\r\n|\r|\n)){1,8}\s*end\s*(?:\r\n|\r|\n)"
    m = match(ex, formula)
    if m === nothing
        # This shouldn't happen, because we passed `has_bottle()` above
        @warn("Couldn't find bottle stanza in $name")
        return name
    end

    # We know there is no `cellar :any` or `cellar :any_skip_relocation` since
    # we made it past the `has_relocatable_bottle()` check.  Eliminate any
    # `cellar` lines still within the match, then replace that section with a
    # section that contains the `cellar :any` we rely so heavily upon.
    bottle_lines = split(m.match, "\n")

    # Eliminate any lines that start with "cellar" or "root_url"
    bottle_lines = filter(line -> !startswith(lstrip(line), "cellar"), bottle_lines)
    bottle_lines = filter(line -> !startswith(lstrip(line), "root_url"), bottle_lines)

    # Find at which line the "bottle do" actually starts
    bottle_idx = findfirst(line -> match(r"bottle\s+do", line) !== nothing, bottle_lines)

    # Add a "cellar :any" and "root_url" line to this formula just after the
    # `bottle do`.  Note that since `match()` returns  SubString's, we need to
    # explicitly convert our string to SubString; this should be fixed in 0.6

    # We should, however, preserve :any_skip_relocation if that is what this
    # bottle is marked as, which includes important bottles such as `cctools`.
    if occursin(":any_skip_relocation", m.match)
        insert!(bottle_lines, bottle_idx+1, SubString("    cellar :any_skip_relocation",1))
    else
        insert!(bottle_lines, bottle_idx+1, SubString("    cellar :any",1))
    end
    insert!(bottle_lines, bottle_idx+1, SubString("    root_url \"$(obj["bottle"]["stable"]["root_url"])\"",1))

    # Resynthesize the bottle stanza and embed it into `formula` once more
    bottle_stanza = join(bottle_lines, "\n")
    formula = formula[1:m.offset-1] * bottle_stanza * formula[m.offset+length(m.match):end]

    # Find any depends_on lines, substitute in any translated formulae
    adjustment = 0
    for m in eachmatch(r"depends_on\s+\"([^ ]+)\"", formula)
        # This is the path that this dependency would have if it has been translated
        dep_name = m.captures[1]
        auto_dep_path = joinpath(auto_tappath, "$(dep_name).rb")

        # If this dependency has been translated, then prepend "staticfloat/juliatranslated"
        # to it in the formula, so there is no confusion inside of Homebrew
        if isfile(auto_dep_path)
            if verbose
                println("translation: replacing dependency $dep_name because it's been translated before")
            end
            new_name = "$(auto_tapname)/$(dep_name)"
            offset = m.offsets[1]

            start_idx = offset-1+adjustment
            stop_idx = offset+length(dep_name)+adjustment
            formula = formula[1:start_idx] * new_name * formula[stop_idx:end]
            adjustment += length(auto_tapname) + 1
        end
    end

    # Write our patched formula out to our override tap
    write_formula(override_name, formula)

    # Read that formula in again as a JSON object, compare with the original
    new_obj = json(override_name)
    if !has_relocatable_bottle(override_name)
        @warn("New formula $override_name doesn't have a relocatable bottle despite our meddling")
        return name
    end

    # Wow.  We actually did it.
    if verbose
        println("translation: successfully finished for $name")
    end
    return override_name
end

function translate_formula(pkg::BrewPkg; verbose::Bool=false)
    return translate_formula(fullname(pkg); verbose=verbose)
end


text:
if Sys.isunix()
    "R eventloop hook on Unix system"
    function polled_events()
        event_callback()
        nothing
    end
end

"Event Callback: allows R to process Julia events when R is busy.
For example, writing output to stdout while running an expensive R command."
function event_callback()
    # dump output buffer to stdout when available
    handle_eval_stdout()
    nothing
end

# there is no use now, maybe useful for the future.
function interrupts_pending(s::Bool=true)
    @static if Sys.iswindows()
        unsafe_store!(cglobal((:UserBreak,libR),Cint), s ? 1 : 0)
    else
        unsafe_store!(cglobal((:R_interrupts_pending,libR),Cint), s ? 1 : 0)
    end
    nothing
end

# this shouldn't exist if we could hook into Julia eventloop.
function process_events()
    ##FIXME: a dirty fix to prevent segfault right after a sigint
    if unsafe_load(cglobal((:R_interrupts_pending,libR),Cint)) == 0
        @static if Sys.iswindows() || Sys.isapple()
            ccall((:R_ProcessEvents, libR), Nothing, ())
        end
        @static if Sys.isunix()
            what = ccall((:R_checkActivity,libR),Ptr{Cvoid},(Cint,Cint),0,1)
            if what != C_NULL
                R_InputHandlers = unsafe_load(cglobal((:R_InputHandlers,libR),Ptr{Cvoid}))
                ccall((:R_runHandlers,libR),Nothing,(Ptr{Cvoid},Ptr{Cvoid}),R_InputHandlers,what)
            end
        end
    end
    nothing
end

global timeout = nothing

function rgui_start(silent=false)
    global timeout
    if timeout == nothing
        timeout = Base.Timer(x -> process_events(), 0.05, interval = 0.05)
        return true
    else
        silent || error("eventloop is already running.")
        return false
    end
end

function rgui_stop(silent=false)
    global timeout
    if timeout != nothing
        close(timeout)
        timeout = nothing
        return true
    else
        silent || error("eventloop is not running.")
        return false
    end
end

function set_hook(hookname, value)
    l = rparse("""setHook(hookname, function(...) foo)""")
    l[1][2] = hookname
    l[1][3][3] = value
    reval(l)
end

function rgui_init()
    f = rlang(rgui_start, true)
    set_hook("plot.new", f)
    set_hook("persp", f)
    set_hook("grid.newpage", f)
    set_hook(rlang(:packageEvent, "rgl", "onLoad"), f)

    # inject rgui_start(TRUE) to utils::help
    help_type = rcopy(rcall(:options, "help_type")[1])
    if  help_type == "html"
        # need to hack both as.environment('package:utils') and  getNamespace("utils")
        # to make ?foo and help("foo") to work
        l = rparse("help <- function(...) { foo(); bar(...) }")
        l[1][3][3][2] = f
        l[1][3][3][3][1] = reval("utils:::help")
        for env in (reval("as.environment('package:utils')"), getNamespace("utils"))
            rcall(:unlockBinding, "help", env)
            reval(l, env)
            rcall(:lockBinding, "help", env)
        end
    end
end


text:
using Documenter, LargeVis
push!(LOAD_PATH,"../src/")

makedocs(sitename="LargeVis.jl Documentation")
deploydocs(
    repo = "https://github.com/hexie1995/LargeVis.jl.git",
    devbranch = "main"
)


text:
# TODO: possibly add GenericSchur to deps and test T == Float64 etc.

@testset "simple eigenvalues $T" for T in (Float32, ComplexF32)
    DT = widen(T)
    maxit = 5
    tol = 20.0
    e = eps(real(T))
    dmin = 1e3 * e
    for n in [8,32]
        A = mkmat_simple(n,dmin,T)
        Ad = convert.(DT,A)
        ewd = eigvals(Ad)
        ew, ev = eigen(A)

        if verbose
            ewerrs = [minimum(abs.(ew[j] .- ewd)) for j in 1:n]
            println("initial errors ", ewerrs / (e * n))
        end
        newvecs = similar(ev)
        newews = similar(ew)
        for j=1:n
            λ, x = rfeigen(A, ev[:,j], ew[j], real(DT), maxiter=maxit)
            newvecs[:,j] .= x
            newews[j] = λ
        end
        # allow for different orderings because of roundoff
        ewerrs = [minimum(abs.(newews[j] .- ewd)) for j in 1:n]
        if verbose
            println("final errors ", ewerrs / (e * n))
        end
        @test maximum(abs.(ewerrs))  < tol * e * n * norm(A)
        # TODO: check newvecs against DP version
    end
end

@testset "(nearly) defective eigenvalues $T" for T in (Float32, ComplexF32)
    DT = widen(T)
    etarget = T(2)
    maxit = 5
    tol = 20.0
    for n in [5,10,32]
        for k in [2,3]
            @label retry
            A = mkmat_defective(n,k,etarget,T)
            # we need a true Schur here
            S = schur(A .+ 0im)
            Ad = convert.(DT,A)
            ew = eigvals(Ad)
            idx = findall(abs.(S.values .- etarget) .< 0.2)
            # don't try if A is so nonnormal that initial estimates are bad
            if length(idx) != k
                @goto retry
            end
            e = eps(real(T))
            if verbose
                ewerrs = [minimum(abs.(S.values[j] .- ew)) for j in idx]
                println("initial errors ", ewerrs / (e * n))
            end
            newew, newvecs = rfeigen(A, S, idx, DT, maxit)
            ewerrs = [minimum(abs.(newew[j] .- ew)) for j in 1:k]
            if verbose
                println("final errors ", ewerrs / (e * n))
            end
            @test maximum(ewerrs) / abs(etarget) < tol * e * n
        end
    end
end

@testset "multiple eigenvalues $T" for T in (Float32, ComplexF32)
    DT = widen(T)
    etarget = T(2)
    dmin = 1e3 * eps(real(T))
    maxit = 5
    tol = 20.0
    for n in [5,10,32]
        for k in [2,3]
            @label retry
            A = mkmat_multiple(n,k,etarget,dmin,T)
            # we need a true Schur here
            S = schur(A .+ 0im)
            Ad = convert.(DT,A)
            ew = eigvals(Ad)
            idx = findall(abs.(S.values .- etarget) .< 0.2)
            # don't try if A is so nonnormal that initial estimates are bad
            if length(idx) != k
                @goto retry
            end
            e = eps(real(T))
            if verbose
                ewerrs = [minimum(abs.(S.values[j] .- ew)) for j in idx]
                println("initial errors ", ewerrs / (e * n))
            end
            newew, newvecs = rfeigen(A, S, idx, DT, maxit)
            ewerrs = [minimum(abs.(newew[j] .- ew)) for j in 1:k]
            if verbose
                println("final errors ", ewerrs / (e * n))
            end
            @test maximum(ewerrs) / abs(etarget) < tol * e * n
        end
    end
end


text:
const BOUND_LOWER = UInt8(1)
const BOUND_UPPER = UInt8(2)
const BOUND_EXACT = UInt8(3)

import Base.hash

# Below we override the base hashing method as we already generated a hash!
# Improves performance for larger depth searches where TT calls are numerous.
Base.hash(x::UInt32) = x
Base.hash(x::UInt16) = x

mutable struct TT_Entry
    eval::Int16
    move::Move
    depth::UInt8
    bound::UInt8
    hash16::UInt16
end
const NO_ENTRY = TT_Entry(0, Move(), 0, 0, 0)

mutable struct TT_Table
    table::Dict{UInt32, TT_Entry}
    hashmask::UInt64
end
TT_Table() = TT_Table(16)

#=
Chess engines usually face memory restrictions, as we can't possibly hash and store all positions.
The below function aims to cap the hash memory size, by restricting the keys in the dict.
=#
function TT_Table(size_MB::Int)
    bytes = size_MB << 20

    # work out size of slots in TT
    #slotsize = Base.summarysize(TT_Entry())
    slotsize = 14 # At present...
    maxslots = bytes / slotsize

    size = 1
    while size <= maxslots
        size <<= 1
    end
    # We shift right again as we "overshifted" by one in the while loop.
    size >>= 1
    size -= 1

    TT_Table(Dict{UInt32, TT_Entry}(), UInt64(size))
end

function getTTentry(tt::TT_Table, hash::ZobristHash)
    res = get(tt.table, UInt32(hash.hash & tt.hashmask), NO_ENTRY)
    if (res !== NO_ENTRY) && (res.hash16 === UInt16(hash.hash >> 48))
        return res
    else
        return NO_ENTRY
    end
end

function setTTentry!(tt::TT_Table, hash::ZobristHash, eval::Integer, move::Move, depth::Integer, bound::UInt8)
    hash16 = UInt16(hash.hash >> 48)
    entry = TT_Entry(eval, move, depth, bound, hash16)
    tt.table[UInt32(hash.hash & tt.hashmask)] = entry
end

function ttvalue(tt_entry::TT_Entry, ply::Int)
    if tt_entry.eval >= (MATE - MAX_PLY)
        return tt_entry.eval + ply
    elseif tt_entry.eval <= (MAX_PLY - MATE)
        return tt_entry.eval - ply
    else
        return tt_entry.eval
    end
end

function hashfull(tt::TT_Table)
    full = 0
    for i in 1:3000
        if haskey(tt.table, UInt32(i))
            full += 1
        end
    end
    fld(full,3)
end

function getPKentry(pktable::PawnKingTable, pkhash::ZobristHash)
    key16 = UInt16(pkhash.hash >> 48)
    if (res = get(pktable, key16, PKT_BLANK)) !== PKT_BLANK
        if res.pkhash === pkhash
            return res
        end
    end
    return PKT_BLANK
end

const PKT_BLANK = PKT_Entry(ZobristHash(0), EMPTY, 0)

function storePKentry!(pktable::PawnKingTable, pkhash::ZobristHash, passed::Bitboard, eval::Int)
    key16 = UInt16(pkhash.hash >> 48)
    pktable[key16] = PKT_Entry(pkhash, passed, eval)
end

max_stars_repo_name: J-Revell/Astellarn

max_issues_repo_name: J-Revell/Astellarn

max_forks_repo_name: J-Revell/Astellarn


text:
using Pivot
using Base.Test

# write your own tests here
@test GET == 1
@test PUT == 2
@test POST == 4
@test DELETE == 8
@test PATCH == 16

@test GET + POST == GET | POST == 5

@testset "Endpoint" begin
end


text:
module TestCompositeParameters

using Mimi
using Test

import Mimi: model_params

@defcomp A begin
    p1 = Parameter(unit = "\$", default=3)
    p2 = Parameter()

    v1 = Variable()
    v2 = Variable(index=[time])
end

@defcomp B begin
    p1 = Parameter(unit = "thous \$", default=3)
    p3 = Parameter(description="B p3")
    p4 = Parameter(index=[time])
    p5 = Parameter(default = 1)
end

#------------------------------------------------------------------------------
# Test a failure to resolve namespace collisions

fail_expr1 = :(
    @defcomposite TestFailComposite begin
        Component(A)
        Component(B)
        # Will fail because you need to resolve the namespace collision of the p1's
    end
)

err1 = try eval(fail_expr1) catch err err end
@test occursin("unresolved parameter name collisions from subcomponents", sprint(showerror, err1))

#------------------------------------------------------------------------------
# Test a failure to resolve conflicting "description" fields

fail_expr2 = :(
    @defcomposite TestFailComposite begin
        Component(A)
        Component(B)

        p1 = Parameter(A.p1, B.p1, unit="")
        p2 = Parameter(A.p2, B.p3)  # Will fail because the conflicting descriptions aren't resolved
    end
)

err2 = try eval(fail_expr2) catch err err end
@test occursin("subcomponents have conflicting values for the \"description\" field", sprint(showerror, err2))

#------------------------------------------------------------------------------
# Test trying to join parameters with different dimensions

fail_expr3 = :(
    @defcomposite TestFailComposite begin
        Component(A)
        Component(B)

        p1 = Parameter(A.p1, B.p1, unit="")
        p2 = Parameter(A.p2, B.p4)  # Will fail because different dimensions
    end
)

err3 = try eval(fail_expr3) catch err err end
@test occursin("subcomponents have conflicting values for the \"dim_names\" field", sprint(showerror, err3))

#------------------------------------------------------------------------------
# Test a failure to auto-import a parameter because it's name has already been used

fail_expr4 = :(
    @defcomposite TestFailComposite begin
        Component(A)
        Component(B)

        p3 = Parameter(A.p1, B.p1, unit="")    # should fail on auto import of other p3 parameter
    end
)

err4 = try eval(fail_expr4) catch err err end
@test occursin("this name has already been defined in the composite component's namespace", sprint(showerror, err4))

#------------------------------------------------------------------------------
# Test an attempt to import a parameter twice

fail_expr5 = :(
    @defcomposite TestFailComposite begin
        Component(A)
        Component(B)

        p1 = Parameter(A.p1, B.p1, unit="")
        p1_repeat = Parameter(B.p1) # should not allow a second import of B.p1 (already connected)
    end
)

err5 = try eval(fail_expr5) catch err err end
@test occursin("Duplicate import", sprint(showerror, err5))

#------------------------------------------------------------------------------
# Test set_param! with unit collision

function get_model()
    m = Model()
    set_dimension!(m, :time, 10)
    add_comp!(m, A)
    add_comp!(m, B)
    return m
end

m1 = get_model()
err6 = try set_param!(m1, :p1, 5) catch err err end
@test occursin("components have conflicting values for the :unit field of this parameter", sprint(showerror, err6))

# use ignoreunits flag
set_param!(m1, :p1, 5, ignoreunits=true)

err7 = try run(m1) catch err err end
@test occursin("Cannot build model; the following parameters still have values of `nothing` and need to be updated:", sprint(showerror, err7))

# Set separate values for p1 in A and B
m2 = get_model()
set_param!(m2, :A, :p1, 2)  # Set the value only for component A

# test that the proper connection has been made for :p1 in :A
@test Mimi.model_param(m2, :p1).value == 2
@test Mimi.model_param(m2, :p1).is_shared
# and that B.p1 is still the default value and unshared
sym = Mimi.get_model_param_name(m2, :B, :p1)
@test Mimi.model_param(m2, sym).value == 3
@test Mimi.model_param(m2, :B, :p1).value == 3
@test !Mimi.model_param(m2, :B, :p1).is_shared

# test defaults
m3 = get_model()
set_param!(m3, :p1, 1, ignoreunits=true)    # Need to set parameter values for all except :p5, which has a default
set_param!(m3, :p2, 2)
set_param!(m3, :p3, 3)
set_param!(m3, :p4, 1:10)
run(m3)

err8 = try set_param!(m3, :B, :p1, 2) catch err err end
@test occursin("the model already has a parameter with this name", sprint(showerror, err8))

set_param!(m3, :B, :p1, :B_p1, 2)   # Use a unique name to set B.p1
@test Mimi.model_param(m3, :B_p1).value == 2
@test Mimi.model_param(m3, :B_p1).is_shared
@test issubset(Set([:p1, :B_p1]), Set(keys(m3.md.model_params)))

#------------------------------------------------------------------------------
# Test update_param! with unit collision

function get_model()
    m = Model()
    set_dimension!(m, :time, 10)
    add_comp!(m, A)
    add_comp!(m, B)
    return m
end

m1 = get_model()
add_shared_param!(m1, :p1, 5)
connect_param!(m1, :A, :p1, :p1) # no conflict
err9 = try connect_param!(m1, :B, :p1, :p1) catch err err end
@test occursin("Cannot connect B:p1 to shared model parameter", sprint(showerror, err9))

# use ignoreunits flag
connect_param!(m1, :B, :p1, :p1, ignoreunits=true)

err10 = try run(m1) catch err err end
@test occursin("Cannot build model; the following parameters still have values of `nothing` and need to be updated:", sprint(showerror, err10))

# Set separate values for p1 in A and B
m2 = get_model()
add_shared_param!(m2, :p1, 2)
connect_param!(m2, :A, :p1, :p1) # Set the value only for component A

# test that the proper connection has been made for :p1 in :A
@test Mimi.model_param(m2.md, :p1).value == 2
@test Mimi.model_param(m2.md, :p1).is_shared
# and that B.p1 is still the default value and unshared
sym = Mimi.get_model_param_name(m2, :B, :p1)
@test Mimi.model_param(m2, sym).value == 3
@test Mimi.model_param(m2, :B, :p1).value == 3
@test !Mimi.model_param(m2, :B, :p1).is_shared

# test defaults - # Need to set parameter values for all except :p5, which has a default
m3 = get_model()
add_shared_param!(m3, :p1, 1)
connect_param!(m3, :A, :p1, :p1, ignoreunits = true)
connect_param!(m3, :B, :p1, :p1, ignoreunits = true)
update_param!(m3, :A, :p2, 2)
update_param!(m3, :B, :p3, 3)
update_param!(m3, :B, :p4, 1:10)
run(m3)

err11 = try add_shared_param!(m3, :p1, 2) catch err err end
@test occursin("the model already has a shared parameter with this name", sprint(showerror, err11))

add_shared_param!(m3, :B_p1, 2) # Use a unique name to set B.p1
connect_param!(m3, :B, :p1, :B_p1)
@test Mimi.model_param(m3, :B_p1).value == 2
@test Mimi.model_param(m3, :B_p1).is_shared
@test issubset(Set([:p1, :B_p1]), Set(keys(m3.md.model_params)))

#------------------------------------------------------------------------------
# Unit tests on default behavior

# different default and override
@defcomp A begin
    p1 = Parameter(default=3)
end
@defcomp B begin
    p1 = Parameter(default=2)
end

@defcomposite top begin
    Component(A)
    Component(B)
    superp1 = Parameter(A.p1, B.p1, default = nothing) # override default collision with nothing
end

m = Model()
set_dimension!(m, :time, 2000:2005)
add_comp!(m, top);
@test length(model_params(m)) == 1
model_param_name = Mimi.get_model_param_name(m.md, :top, :superp1)
@test Mimi.is_nothing_param(model_params(m)[model_param_name])
@test !model_params(m)[model_param_name].is_shared

@defcomposite top begin
    Component(A)
    Component(B)
    superp1 = Parameter(A.p1, B.p1, default = 8.0) # override default collision with value
end

m = Model()
set_dimension!(m, :time, 2000:2005)
add_comp!(m, top);
@test length(model_params(m)) == 1
model_param_name = Mimi.get_model_param_name(m.md, :top, :superp1)
@test model_params(m)[model_param_name].value == 8.0
@test !model_params(m)[model_param_name].is_shared

# same default and no override
@defcomp A begin
    p1 = Parameter(default=2)
end
@defcomp B begin
    p1 = Parameter(default=2)
end

@defcomposite top begin
    Component(A)
    Component(B)
    superp1 = Parameter(A.p1, B.p1)
end

m = Model()
set_dimension!(m, :time, 2000:2005)
add_comp!(m, top);
@test length(model_params(m)) == 1
model_param_name = Mimi.get_model_param_name(m.md, :top, :superp1)
@test model_params(m)[model_param_name].value == 2
@test !model_params(m)[model_param_name].is_shared

# simple case with no super parameter
@defcomp A begin
    p1 = Parameter(default=2)
end
@defcomp B begin
    p2 = Parameter(default=3)
end

@defcomposite top begin
    Component(A)
    Component(B)
end

m = Model()
set_dimension!(m, :time, 2000:2005)
add_comp!(m, top);
@test length(model_params(m)) == 2
model_param_name = Mimi.get_model_param_name(m.md, :top, :p1)
@test model_params(m)[model_param_name].value == 2
@test !model_params(m)[model_param_name].is_shared
model_param_name = Mimi.get_model_param_name(m.md, :top, :p2)
@test model_params(m)[model_param_name].value == 3
@test !model_params(m)[model_param_name].is_shared

#------------------------------------------------------------------------------
# Test set_param! for parameter that exists in neither model definition nor any subcomponent

m1 = get_model()
err12 = try set_param!(m1, :pDNE, 42) catch err err end
@test occursin("not found in ModelDef or children", sprint(showerror, err12))

# Test update_param! for parameter that exists in neither model definition nor any subcomponent
err13 = try update_param!(m1, :pDNE, 42) catch err err end
@test occursin("not found in composite's model parameters", sprint(showerror, err13))

end #module


text:
using Test, YaoBase
using YaoBase: @interface, handle

if VERSION < v"1.1.0"
    isnothing(x) = x === nothing
end

ex = [
    :(function foo end),
    :(foo(x) = x),
    :(foo(x, y::T) where T),
    :(foo(x, y::T) where T = x)]

function check_handle(ex, name, args, body)
    r = handle(ex)
    r_name, r_args, r_body = r

    flag = name == r_name
    if args == true
        flag = flag && !isnothing(r_args)
    else
        flag = flag && (args == r_args)
    end

    if body == true
        return flag && !isnothing(r_body)
    else
        return flag && isnothing(r_body)
    end
end

@test check_handle(ex[1], :foo, nothing, false)
@test check_handle(ex[2], :foo, [:x], true)
@test check_handle(ex[3], :foo, [:x, :y], false)
@test check_handle(ex[4], :foo, [:x, :y], true)

@test_throws ErrorException @interface x + 1

@test check_handle(:(cache_type(::Type{<:AbstractBlock})), :cache_type, true, false)


text:
function mc(O::Option; N::Int64 = 1000, P::Int64 = 10000)
    Δt = O.t / N
    R = exp(O.r * Δt)
    T = typeof(O.s * exp(-O.σ^2 * Δt / 2 + O.σ * √Δt * 0.1) / R)
    X = Array{T}(N+1, P)

    for p = 1:P
        X[1, p] = x = O.s
        for n = 1:N
            x *= R * exp(-O.σ^2 * Δt / 2 + O.σ * √Δt * randn())
            X[n+1, p] = x
        end
    end

    if O.call_or_put == "Put"
        V = [max(O.k - x, 0) / R for x in X[N+1, :]]
    else
        V = [max(x - O.k, 0) / R for x in X[N+1, :]]
    end

    for n = N-1:-1:1
        I = V .!= 0
        A = [x^d for d = 0:3, x in X[n+1, :]]
        β = A[:, I]' \ V[I]
        cV = A' * β
        for p = 1:P
            ev = max(O.k - X[n+1, p], 0)
            if I[p] && cV[p] < ev
                V[p] = ev / R
            else
                V[p] /= R
            end
        end
    end

    if O.call_or_put == "Put"
        return max(mean(V), O.k - O.s)
    else
        return max(mean(V), O.s - O.k)
    end
end


text:
function jacobian(m::HetDSGESimpleTaylor)

    # Load in endogenous state and eq cond indices
    endo = m.endogenous_states #=augment_model_states(m.endogenous_states_unnormalized,
                                n_model_states_unnormalized(m)) =#
    eq   = m.equilibrium_conditions

    # Load in parameters, steady-state parameters, and grids
    r::Float64     = m[:r].value
    α::Float64     = m[:α].value
    H::Float64     = m[:H].value
    δ::Float64     = m[:δ].value
    μ_sp::Float64  = m[:μ_sp].value
    ρ_sp::Float64  = m[:ρ_sp].value
    σ_sp::Float64  = m[:σ_sp].value
    γ::Float64     = m[:γ].value
    g::Float64     = m[:g].value
    η::Float64     = m[:η].value
    ρB::Float64    = m[:ρB].value
    ρG::Float64    = m[:ρG].value
    ρZ::Float64    = m[:ρZ].value
    ρμ::Float64    = m[:ρμ].value
    ρlamw::Float64 = m[:ρlamw].value
    ρlamf::Float64 = m[:ρlamf].value
    ρmon::Float64  = m[:ρmon].value
    spp::Float64   = m[:spp].value
    lamw::Float64  = m[:lamw].value
    ϕh::Float64    = m[:ϕh].value
    Φw::Float64    = m[:Φw].value
    lamf::Float64  = m[:lamf].value
    Φp::Float64    = m[:Φp].value
    ρR::Float64    = m[:ρR].value
    ψπ::Float64    = m[:ψπ].value
    ψy::Float64    = m[:ψy].value

    R = 1 + r

    ell::Vector{Float64}  = m[:lstar].value
    c::Vector{Float64}    = m[:cstar].value
    μ::Vector{Float64}    = m[:μstar].value
    β::Float64            = m[:βstar].value

    T::Float64     = m[:Tstar].value
    ω::Float64     = m[:ωstar].value
    xstar::Float64 = m[:xstar].value
    ystar::Float64 = m[:ystar].value
    Rk::Float64    = m[:Rkstar].value
    kstar::Float64 = m[:kstar].value

    xgrid::Vector{Float64} = m.grids[:xgrid].points
    xwts::Vector{Float64}  = m.grids[:xgrid].weights
    sgrid::Vector{Float64} = m.grids[:sgrid].points
    swts::Vector{Float64}  = m.grids[:sgrid].weights
    fgrid::Matrix{Float64} = m.grids[:fgrid]

    xswts = kron(swts,xwts)

    zlo::Float64 = get_setting(m, :zlo)
    zhi::Float64 = get_setting(m, :zhi)

    nx::Int = get_setting(m, :nx)
    ns::Int = get_setting(m, :ns)

    xlo = get_setting(m, :xlo)
#####
    mindens = 1e-8
    rescale_xwts = true
    truncate = true
    if truncate == true
        oldnx = nx
        nx = maximum(find(μ[1:nx]+μ[nx+1:2*nx].>mindens)) # used to be 1e-8
        μ = μ[[1:nx;oldnx+1:oldnx+nx]]
        ell = ell[[1:nx;oldnx+1:oldnx+nx]]
        # experiment with this:
        if rescale_xwts == true
            xgrid   = xgrid[1:nx] #Evenly spaced grid
            xhi = xgrid[nx]
            xscale = xhi-xlo
            xwts     = (xscale/nx)*ones(nx)          #quadrature weights
            xswts = kron(swts,xwts)
        end
    end

    m <= Setting(:nx, nx)
    m <= Setting(:xhi, xhi)
    m <= Setting(:xscale, xscale)
    m.grids[:xgrid] =  Grid(uniform_quadrature(xscale), xlo, xhi, nx, scale = xscale)
    xswts = kron(swts,xwts)

    nxns = nx*ns
    nvars = 4*nxns+35
    nscalars = 35 # number of eqs which output scalars
    nyscalars = 13 # number of scalar jumps
    nxscalars = nscalars - nyscalars # number of scalar states
    m <= Setting(:nxns, nxns)
    m <= Setting(:nvars, nvars)
    m <= Setting(:nscalars, nscalars)
    m <= Setting(:nyscalars, nyscalars)
    m <= Setting(:nxscalars, nxscalars)

    endo[:l′_t1]    = 1:nxns            # lagged ell function
    endo[:μ′_t1]    = nxns+1:2*nxns     # lagged wealth distribution
    #endogenous scalar-valued states
    endo[:k′_t]  = 2*nxns+1:2*nxns+1             # capital –dont get confused with steadystate object K
    endo[:R′_t1] = 2*nxns+2:2*nxns+2              # lagged real interest rate
    endo[:i′_t1] = 2*nxns+3:2*nxns+3             # lagged nominal interest rate
    endo[:π′_t1] = 2*nxns+4:2*nxns+4             # lagged inflation
    endo[:π′_t2] = 2*nxns+5:2*nxns+5             # double-lagged inflation
    endo[:π′_t3] = 2*nxns+6:2*nxns+6             # triple-lagged inflation
    endo[:y′_t1] = 2*nxns+7:2*nxns+7             # lagged gdp
    endo[:y′_t2] = 2*nxns+8:2*nxns+8             # double-lagged gdp
    endo[:y′_t3] = 2*nxns+9:2*nxns+9             # triple-lagged gdp
    endo[:y′_t4] = 2*nxns+10:2*nxns+10           # quad-lagged gdp
    endo[:z′_t1] = 2*nxns+11:2*nxns+11           # lag tfp growth
    endo[:z′_t2] = 2*nxns+12:2*nxns+12           # double-lag tfp growth
    endo[:z′_t3] = 2*nxns+13:2*nxns+13           # triple-lag tfp growth
    endo[:w′_t1] = 2*nxns+14:2*nxns+14           # lag real wages
    endo[:I′_t1] = 2*nxns+15:2*nxns+15           # lag investment–don't get this confused with i, the nominal interest rate
    # exogenous scalar-valued states:
    endo[:BP]    = 2*nxns+16:2*nxns+16        # discount factor shock
    endo[:GP]    = 2*nxns+17:2*nxns+17        # govt spending
    endo[:ZP]    = 2*nxns+18:2*nxns+18        # tfp growth
    endo[:MUP]   = 2*nxns+19:2*nxns+19        # investment shock
    endo[:LAMWP] = 2*nxns+20:2*nxns+20        # wage markup
    endo[:LAMFP] = 2*nxns+21:2*nxns+21        # price markup
    endo[:MONP]  = 2*nxns+22:2*nxns+22        # monetary policy shock

 # function-valued jumps
    endo[:l′_t]  = 2*nxns+23:3*nxns+22 # ell function
    endo[:μ′_t]  = 3*nxns+23:4*nxns+22 # wealth distribution
    #scalar-valued jumps
    endo[:R′_t]  = 4*nxns+23:4*nxns+23        # real interest rate
    endo[:i′_t]  = 4*nxns+24:4*nxns+24        # nominal interest rate
    endo[:t′_t]  = 4*nxns+25:4*nxns+25        # transfers + dividends
    endo[:w′_t]  = 4*nxns+26:4*nxns+26        # real wage
    endo[:L′_t]  = 4*nxns+27:4*nxns+27        # hours worked
    endo[:π′_t]  = 4*nxns+28:4*nxns+28        # inflation
    endo[:wageinflation′_t] = 4*nxns+29:4*nxns+29        # nominal wage inflation
    endo[:mu′_t] =  4*nxns+30:4*nxns+30        # average marginal utility
    endo[:y′_t]  = 4*nxns+31:4*nxns+31       # gdp
    endo[:I′_t]  = 4*nxns+32:4*nxns+32        # investment
    endo[:mc′_t] = 4*nxns+33:4*nxns+33        # marginal cost - this is ζ in HetDSGEₖd.pdf
    endo[:Q′_t]  = 4*nxns+34:4*nxns+34        # Tobin's qfunction
    endo[:capreturn′_t] = 4*nxns+35:4*nxns+35        # return on capital

endo[:l_t1] = endo[:l′_t1] .+ nvars #LELL
endo[:μ_t1] = endo[:μ′_t1] .+ nvars #LM
endo[:k_t] = endo[:k′_t] .+ nvars #KK
endo[:R_t1] = endo[:R′_t1] .+ nvars #LRR
endo[:i_t1] = endo[:i′_t1] .+ nvars #LII
endo[:π_t1] = endo[:π′_t1]  .+ nvars #LPI
endo[:π_t2] = endo[:π′_t2]  .+ nvars #L2PI
endo[:π_t3] = endo[:π′_t3]  .+ nvars #L3PI
endo[:y_t1] = endo[:y′_t1] .+ nvars #LY
endo[:y_t2] = endo[:y′_t2] .+ nvars #L2Y
endo[:y_t3] = endo[:y′_t3] .+ nvars #l3Y
endo[:y_t4] = endo[:y′_t4] .+ nvars #L4Y
endo[:z_t1] = endo[:z′_t1] .+ nvars #LZ
endo[:z_t2] = endo[:z′_t2] .+ nvars #L2Z
endo[:z_t3] = endo[:z′_t3] .+ nvars #L3Z
endo[:w_t1] = endo[:w′_t1] .+ nvars #LW
endo[:I_t1] = endo[:I′_t1] .+ nvars #LX
endo[:B] = endo[:BP] .+ nvars #B
endo[:G] = endo[:GP] .+ nvars #G
endo[:Z] = endo[:ZP] .+ nvars #Z
endo[:MU] = endo[:MUP] .+ nvars #MU
endo[:LAMW] = endo[:LAMWP] .+ nvars #LAMW
endo[:LAMF] = endo[:LAMFP] .+ nvars #LAMF
endo[:MON] = endo[:MONP] .+ nvars #MON
endo[:ELL] = endo[:l′_t] .+ nvars #ELL
endo[:μ_t] = endo[:μ′_t] .+ nvars #M
endo[:R_t] = endo[:R′_t] .+ nvars #RR
endo[:i_t] = endo[:i′_t] .+ nvars #II
endo[:t_t] = endo[:t′_t]  .+ nvars #TT
endo[:w_t] = endo[:w′_t] .+ nvars #W
endo[:L_t] = endo[:L′_t] .+ nvars #HH
endo[:π_t] = endo[:π′_t] .+ nvars #PI
endo[:wageinflation_t] = endo[:wageinflation′_t] .+ nvars #PIW
endo[:mu_t] = endo[:mu′_t] .+ nvars #LAM
endo[:y_t] = endo[:y′_t] .+ nvars #Y
endo[:I_t] = endo[:I′_t]  .+ nvars #X
endo[:mc_t] = endo[:mc′_t]  .+ nvars #MC
endo[:Q_t] = endo[:Q′_t]  .+ nvars #Q
endo[:capreturn_t] = endo[:capreturn′_t] .+ nvars #RK

eqconds = m.equilibrium_conditions
    # function blocks which output a function
    eqconds[:eq_euler]              = 1:nxns
    eqconds[:eq_kolmogorov_fwd]     = nx*ns+1:2*nx*ns
    eqconds[:eq_lag_ell]            = 2*nxns+1:3*nxns
    eqconds[:eq_lag_wealth]         = 3*nxns+1:4*nxns
    # function blocks which map functions to scalars
    eqconds[:eq_market_clearing]    = 4*nxns+1:4*nxns+1
    eqconds[:eq_lambda]             = 4*nxns+2:4*nxns+2
    #scalar blocks involving endogenous variables
    eqconds[:eq_transfers]            = 4*nxns+3:4*nxns+3 # transfers
    eqconds[:eq_investment]          = 4*nxns+4:4*nxns+4 # investment
    eqconds[:eq_tobin_q]              = 4*nxns+5:4*nxns+5 # tobin's q
    eqconds[:eq_capital_accumulation] = 4*nxns+6:4*nxns+6 # capital accumulation
    eqconds[:eq_wage_phillips] = 4*nxns+7:4*nxns+7 # wage phillips curve
    eqconds[:eq_price_phillips] = 4*nxns+8:4*nxns+8 # price phillips curve
    eqconds[:eq_marginal_cost]  = 4*nxns+9:4*nxns+9 # marginal cost
    eqconds[:eq_gdp]  = 4*nxns+10:4*nxns+10 # gdp
    eqconds[:eq_optimal_kl] = 4*nxns+11:4*nxns+11 # optimal K/L ratio
    eqconds[:eq_taylor] = 4*nxns+12:4*nxns+12 # taylor rule
    eqconds[:eq_fisher] = 4*nxns+13:4*nxns+13 # fisher eqn
    eqconds[:eq_nominal_wage_inflation] = 4*nxns+14:4*nxns+14 # nominal wage inflation
    # lagged variables
    eqconds[:LR] = 4*nxns+15:4*nxns+15 # LR
    eqconds[:LI] = 4*nxns+16:4*nxns+16 # LI
    eqconds[:LPI] = 4*nxns+17:4*nxns+17 # LPI
    eqconds[:L2PI] = 4*nxns+18:4*nxns+18 # L2PI
    eqconds[:L3PI] = 4*nxns+19:4*nxns+19 # L3PI
    eqconds[:LY] = 4*nxns+20:4*nxns+20 # LY
    eqconds[:L2Y] = 4*nxns+21:4*nxns+21 # L2Y
    eqconds[:L3Y] = 4*nxns+22:4*nxns+22 # L3Y
    eqconds[:L4Y] = 4*nxns+23:4*nxns+23 # L4Y
    eqconds[:LZ]  = 4*nxns+24:4*nxns+24 # LZ
    eqconds[:L2Z] = 4*nxns+25:4*nxns+25 # L2Z
    eqconds[:L3Z] = 4*nxns+26:4*nxns+26 # L3Z
    eqconds[:LW]  = 4*nxns+27:4*nxns+27 # LW
    eqconds[:LX] = 4*nxns+28:4*nxns+28 # LX
    # shocks
    eqconds[:F33] = 4*nxns+29:4*nxns+29 # discount factor B
    eqconds[:F34] = 4*nxns+30:4*nxns+30 # govt spending G
    eqconds[:F35] = 4*nxns+31:4*nxns+31 # tfp growth Z
    eqconds[:F36] = 4*nxns+32:4*nxns+32 # investment MU
    eqconds[:F37] = 4*nxns+33:4*nxns+33 # wage mkup LAMW
    eqconds[:F38] = 4*nxns+34:4*nxns+34 # price mkup LAMF
    eqconds[:F39] = 4*nxns+35:4*nxns+35 # monetary policy MON

####

    zgrid = range(zlo,stop = zhi,length = nx)
    zwts = (zhi-zlo)/nx
    sumz=0.
    for i=1:nx
        sumz += mollifier_hetdsge(zgrid[i],zhi,zlo)*zwts
    end

    qp(z) = dmollifier_hetdsge(z, zhi, zlo)
    qfunction(x) = mollifier_hetdsge(x, zhi, zlo) #/sumz

    unc = 1 ./ ell .<= repeat(xgrid,ns) .+ η

    # Euler Equation
    dF1_dELL, dF1_dRZ, dF1_dELLP, dF1_dWHP, dF1_dTTP, ee =
        euler_equation_hetdsge(nx, ns, qp, qfunction, xgrid, sgrid, fgrid, unc, xswts,
                               R, γ, β, η, ell, T, ω, H)

    # KF Equation
    dF2_dWH, dF2_dRZ, dF2_dTT,dF2_dELL, bigΨ, dF2_dM =
        kolmogorov_fwd_hetdsge(nx, ns, qfunction, qp, xgrid, sgrid, fgrid, unc, xswts,
                               R, γ, ell, μ, η, T, ω, H, ee)

    # Market clearing, lambda function
    c = min.(1 ./ ell,repeat(xgrid,ns).+η)
    lam = (xswts.*μ)'*(1 ./ c) # average marginal utility which the union uses to set wages
    ϕ = lam*ω/(H^ϕh) # now that we know lam in steady state, choose disutility to target hours H

    nvars = get_setting(m, :nvars)
    # Make the Jacobian
    JJ = zeros(nvars, 2*nvars)

    # Euler equation
    JJ[eq[:eq_euler],endo[:l′_t]] = dF1_dELLP
    JJ[eq[:eq_euler],endo[:ZP]]   = -dF1_dRZ
    JJ[eq[:eq_euler],endo[:w′_t]]   = dF1_dWHP
    JJ[eq[:eq_euler],endo[:L′_t]]  = dF1_dWHP
    JJ[eq[:eq_euler],endo[:t′_t]]  = dF1_dTTP
    JJ[eq[:eq_euler],endo[:B]]    = ell
    JJ[eq[:eq_euler],endo[:ELL]]  = dF1_dELL
    JJ[eq[:eq_euler],endo[:R_t]]   = dF1_dRZ

    # KF eqn
    JJ[eq[:eq_kolmogorov_fwd],endo[:μ_t1]]   = dF2_dM
    JJ[eq[:eq_kolmogorov_fwd],endo[:R_t1]]  = dF2_dRZ
    JJ[eq[:eq_kolmogorov_fwd],endo[:l_t1]] = dF2_dELL
    JJ[eq[:eq_kolmogorov_fwd],endo[:Z]]    = -dF2_dRZ
    JJ[eq[:eq_kolmogorov_fwd],endo[:μ_t]]    = Matrix{Float64}(I, nxns, nxns)
    JJ[eq[:eq_kolmogorov_fwd],endo[:w_t]]    = dF2_dWH
    JJ[eq[:eq_kolmogorov_fwd],endo[:L_t]]   = dF2_dWH
    JJ[eq[:eq_kolmogorov_fwd],endo[:t_t]]   = dF2_dTT

    # update lagged ELL
    JJ[eq[:eq_lag_ell],endo[:l′_t1]] = Matrix{Float64}(I, nxns, nxns)
    JJ[eq[:eq_lag_ell],endo[:ELL]]   = -Matrix{Float64}(I, nxns, nxns)

    # update lagged M
    JJ[eq[:eq_lag_wealth],endo[:μ′_t1]] = Matrix{Float64}(I, nxns, nxns)
    JJ[eq[:eq_lag_wealth],endo[:μ_t]]   = -Matrix{Float64}(I, nxns, nxns)

    # mkt clearing
    JJ[first(eq[:eq_market_clearing]),first(endo[:y_t])]   = ystar/g
    JJ[first(eq[:eq_market_clearing]),first(endo[:G])]   = -ystar/g
    JJ[first(eq[:eq_market_clearing]),first(endo[:I_t])]   = -xstar
    JJ[eq[:eq_market_clearing],endo[:ELL]] = (μ.*unc.*xswts.*c)'
    JJ[eq[:eq_market_clearing],endo[:μ_t]]   = -(xswts.*c)' # note, now we linearize

    # lambda = average marginal utility
    JJ[first(eq[:eq_lambda]),first(endo[:mu_t])] = lam
    JJ[first(eq[:eq_lambda]),endo[:μ_t]]   = -(xswts./c)' # note, now we linearize
    JJ[first(eq[:eq_lambda]),endo[:ELL]] = -(xswts.*unc.*μ./c)'

    # transfer
    JJ[first(eq[:eq_transfers]),first(endo[:t_t])]  = T
    JJ[first(eq[:eq_transfers]),first(endo[:capreturn_t])]  = -Rk*kstar
    JJ[first(eq[:eq_transfers]),first(endo[:k_t])]  = -Rk*kstar
    JJ[first(eq[:eq_transfers]),first(endo[:Z])]  = Rk*kstar
    JJ[first(eq[:eq_transfers]),first(endo[:I_t])]   = xstar
    JJ[first(eq[:eq_transfers]),first(endo[:mc_t])]  = ystar
    JJ[first(eq[:eq_transfers]),first(endo[:y_t])]   = (1-1/g)
    JJ[first(eq[:eq_transfers]),first(endo[:G])]   = (ystar/g)

    # investment
    JJ[first(eq[:eq_investment]),first(endo[:Q_t])]  = 1.
    JJ[first(eq[:eq_investment]),first(endo[:MU])] = 1.
    JJ[first(eq[:eq_investment]),first(endo[:I′_t])] = spp*(exp(3*γ))/R
    JJ[first(eq[:eq_investment]),first(endo[:ZP])] = spp*(exp(3*γ))/R
    JJ[first(eq[:eq_investment]),first(endo[:I_t])]  = -spp*(exp(3*γ))/R - spp*exp(2*γ)
    JJ[first(eq[:eq_investment]),first(endo[:I_t1])] = spp*exp(2*γ)
    JJ[first(eq[:eq_investment]),first(endo[:Z])]  = -spp*exp(2*γ)

    # tobin's q
    JJ[first(eq[:eq_tobin_q]),first(endo[:R_t])]  = R
    JJ[first(eq[:eq_tobin_q]),first(endo[:Q_t])]   = R
    JJ[first(eq[:eq_tobin_q]),first(endo[:capreturn′_t])] = -Rk
    JJ[first(eq[:eq_tobin_q]),first(endo[:Q′_t])]  = -(1-δ)

    # capital accumulation
    JJ[first(eq[:eq_capital_accumulation]),first(endo[:k′_t])] = 1.
    JJ[first(eq[:eq_capital_accumulation]),first(endo[:k_t])]  = -(1-δ)
    JJ[first(eq[:eq_capital_accumulation]),first(endo[:Z])]   = (1-δ)
    JJ[first(eq[:eq_capital_accumulation]),first(endo[:MU])]  = -xstar/kstar
    JJ[first(eq[:eq_capital_accumulation]),first(endo[:I_t])]   = -xstar/kstar

    # wage phillips curve
    JJ[first(eq[:eq_wage_phillips]),first(endo[:wageinflation_t])]  = -1.
    JJ[first(eq[:eq_wage_phillips]),first(endo[:LAMW])] = (ϕ*H^ϕh)/Φw
    JJ[first(eq[:eq_wage_phillips]),first(endo[:L_t])]   = (ϕ*(H^ϕh)*(1+lamw)/lamw*Φw)*ϕh
    JJ[first(eq[:eq_wage_phillips]),first(endo[:mu_t])]  = -(ϕ*(H^ϕh)*(1+lamw)/lamw*Φw)
    JJ[first(eq[:eq_wage_phillips]),first(endo[:w_t])]    = -(ϕ*(H^ϕh)*(1+lamw)/lamw*Φw)
    JJ[first(eq[:eq_wage_phillips]),first(endo[:wageinflation′_t])]  = β

    # price phillips curve
    JJ[first(eq[:eq_price_phillips]),first(endo[:π_t])]   = -1.
    JJ[first(eq[:eq_price_phillips]),first(endo[:mc_t])]  = (1+lamf)/(lamf*Φp)
    JJ[first(eq[:eq_price_phillips]),first(endo[:LAMF])] = 1/Φp
    JJ[first(eq[:eq_price_phillips]),first(endo[:π′_t])]  = 1/R

    # marginal cost
    JJ[first(eq[:eq_marginal_cost]),first(endo[:mc_t])] = 1.
    JJ[first(eq[:eq_marginal_cost]),first(endo[:w_t])]  = -(1-α)
    JJ[first(eq[:eq_marginal_cost]),first(endo[:capreturn_t])] = -α

    # gdp
    JJ[first(eq[:eq_gdp]),first(endo[:y_t])]  = 1.
    JJ[first(eq[:eq_gdp]),first(endo[:Z])]  = α
    JJ[first(eq[:eq_gdp]),first(endo[:k_t])] = -α
    JJ[first(eq[:eq_gdp]),first(endo[:L_t])] = -(1-α)

    # optimal k/l ratio
    JJ[first(eq[:eq_optimal_kl]),first(endo[:capreturn_t])] = 1.
    JJ[first(eq[:eq_optimal_kl]),first(endo[:w_t])]  = -1.
    JJ[first(eq[:eq_optimal_kl]),first(endo[:L_t])] = -1.
    JJ[first(eq[:eq_optimal_kl]),first(endo[:k_t])] = 1.
    JJ[first(eq[:eq_optimal_kl]),first(endo[:Z])]  = -1.

    # taylor rule
    JJ[first(eq[:eq_taylor]),first(endo[:i_t])]   = -1.
    JJ[first(eq[:eq_taylor]),first(endo[:π_t])]   = ψπ
    JJ[first(eq[:eq_taylor]),first(endo[:MONP])] = 1.

    # fisher eqn
    JJ[first(eq[:eq_fisher]),first(endo[:R_t])]  = 1.
    JJ[first(eq[:eq_fisher]),first(endo[:π′_t])] = 1.
    JJ[first(eq[:eq_fisher]),first(endo[:i_t])]  = -1.

    # wage inflation
    JJ[first(eq[:eq_nominal_wage_inflation]),first(endo[:wageinflation_t])] = 1.
    JJ[first(eq[:eq_nominal_wage_inflation]),first(endo[:π_t])]  = -1.
    JJ[first(eq[:eq_nominal_wage_inflation]),first(endo[:Z])]   = -1.
    JJ[first(eq[:eq_nominal_wage_inflation]),first(endo[:w_t])]   = -1.
    JJ[first(eq[:eq_nominal_wage_inflation]),first(endo[:w_t1])]  = 1.

    # update lagged variables
    JJ[first(eq[:LR]),first(endo[:R′_t1])] = 1.
    JJ[first(eq[:LR]),first(endo[:R_t])]   = -1.

    JJ[first(eq[:LI]),first(endo[:i′_t1])] = 1.
    JJ[first(eq[:LI]),first(endo[:i_t])]   = -1.

    JJ[first(eq[:LPI]),first(endo[:π′_t1])] = 1.
    JJ[first(eq[:LPI]),first(endo[:π_t])]   = -1.

    JJ[first(eq[:L2PI]),first(endo[:π′_t2])] = 1.
    JJ[first(eq[:L2PI]),first(endo[:π_t1])]   = -1.

    JJ[first(eq[:L3PI]),first(endo[:π′_t3])] = 1.
    JJ[first(eq[:L3PI]),first(endo[:π_t2])]  = -1.

    JJ[first(eq[:LY]),first(endo[:y′_t1])] = 1.
    JJ[first(eq[:LY]),first(endo[:y_t])]   = -1.

    JJ[first(eq[:L2Y]),first(endo[:y′_t2])] = 1.
    JJ[first(eq[:L2Y]),first(endo[:y_t1])]   = -1.

    JJ[first(eq[:L3Y]),first(endo[:y′_t3])] = 1.
    JJ[first(eq[:L3Y]),first(endo[:y_t2])]  = -1.

    JJ[first(eq[:L4Y]),first(endo[:y′_t4])] = 1.
    JJ[first(eq[:L4Y]),first(endo[:y_t3])]  = -1.

    JJ[first(eq[:LZ]),first(endo[:z′_t1])] = 1.
    JJ[first(eq[:LZ]),first(endo[:Z])]   = -1.

    JJ[first(eq[:L2Z]),first(endo[:z′_t2])] = 1.
    JJ[first(eq[:L2Z]),first(endo[:z_t1])]   = -1.

    JJ[first(eq[:L3Z]),first(endo[:z′_t3])] = 1.
    JJ[first(eq[:L3Z]),first(endo[:z_t2])]  = -1.

    JJ[first(eq[:LW]),first(endo[:w′_t1])] = 1.
    JJ[first(eq[:LW]),first(endo[:w_t])]   = -1.

    JJ[first(eq[:LX]),first(endo[:I′_t1])] = 1.
    JJ[first(eq[:LX]),first(endo[:I_t])]   = -1.

    # discount factor shock
    JJ[first(eq[:F33]),first(endo[:BP])] = 1.
    JJ[first(eq[:F33]),first(endo[:B])]  = -ρB

    # g/y shock
    JJ[first(eq[:F34]),first(endo[:GP])] = 1.
    JJ[first(eq[:F34]),first(endo[:G])]  = -ρG

    # tfp growth shock
    JJ[first(eq[:F35]),first(endo[:ZP])] = 1.
    JJ[first(eq[:F35]),first(endo[:Z])]  = -ρZ

    # investment shock
    JJ[first(eq[:F36]),first(endo[:MUP])] = 1.
    JJ[first(eq[:F36]),first(endo[:MU])]  = -ρμ

    # wage mkup shock
    JJ[first(eq[:F37]),first(endo[:LAMWP])] = 1.
    JJ[first(eq[:F37]),first(endo[:LAMW])]  = -ρlamw

    # price mkup shock
    JJ[first(eq[:F38]),first(endo[:LAMFP])] = 1.
    JJ[first(eq[:F38]),first(endo[:LAMF])]  = -ρlamf

    # monetary policy shock
    JJ[first(eq[:F39]),first(endo[:MONP])] = 1.
    JJ[first(eq[:F39]),first(endo[:MON])]  = -ρmon

    if !m.testing && get_setting(m, :normalize_distr_variables)
        JJ = normalize(m, JJ)
    end

    return JJ
end

function euler_equation_hetdsge_simple_taylor(nx::Int, ns::Int,
                                qp::Function, qfunction::Function,
                                xgrid::Vector{Float64}, sgrid::Vector{Float64},
                                fgrid::Matrix{Float64},
                                unc::BitArray,
                                xswts::Vector{Float64},
                                R::Float64, γ::Float64, β::Float64,
                                η::Float64, ell::Vector{Float64}, T::Float64,
                                ω::Float64, H::Float64)
    nxns = nx*ns
    ee  = zeros(nxns,nxns) # ee[i,j] takes you from i to j
    ξ   = zeros(nxns,nxns)
    Ξ   = zeros(nxns,nxns)
    dF1_dELL = zeros(nxns,nxns)
    dF1_dRZ  = zeros(nxns)
    dF1_dELLP = zeros(nxns,nxns)
    dF1_dWHP = zeros(nxns)
    dF1_dTTP = zeros(nxns)
    ellRHS = zeros(nxns)
    for iss=1:ns
        for ia=1:nx
            i  = nx*(iss-1)+ia
            sumELL = 0.
            sumRZ  = 0.
            sumWH  = 0.
            sumTT  = 0.
            sum_ellRHS = 0.
            for isp=1:ns
                sum_ellRHS = 0.
                for iap=1:nx
                    ip = nx*(isp-1)+iap
                    ee[i,ip] = (xgrid[iap] - R*(exp(-γ))*max(xgrid[ia]-1/ell[i], -η) - T)/(ω*H*sgrid[isp])
                    ξ[i,ip] = ((β*R*xswts[i]*exp(-γ))/(ω*H*sgrid[isp])^2)*max(ell[ip],1/(xgrid[iap]+η))*qp(ee[i,ip])*fgrid[iss,isp]
                    Ξ[i,ip] = ((β*R*xswts[i]*exp(-γ))/(ω*H*sgrid[isp]))*max(ell[ip],1/(xgrid[iap]+η))*qfunction(ee[i,ip])*fgrid[iss,isp]
                    sumELL += ξ[i,ip]*R*(exp(-γ))*unc[i]/ell[i]
                    sumRZ  += ξ[i,ip]*R*(exp(-γ))*max(xgrid[ia] - 1/ell[i],-η)
                    dF1_dELLP[i,ip] = Ξ[i,ip]*unc[ip]
                    sumWH  += Ξ[i,ip] + ξ[i,ip]*ee[i,ip]/(ω*H*sgrid[isp])
                    sumTT  += ξ[i,ip]*T
                    sum_ellRHS += Ξ[i,ip]
                end
            end
            ellRHS[i] = sum_ellRHS
            dF1_dELL[i,i] = -ell[i] - sumELL
            dF1_dRZ[i] = ellRHS[i] - sumRZ
            dF1_dWHP[i] = -sumWH
            dF1_dTTP[i] = -sumTT
        end
    end
    return dF1_dELL, dF1_dRZ, dF1_dELLP, dF1_dWHP, dF1_dTTP, ee
end

function kolmogorov_fwd_hetdsge_simple_taylor(nx::Int, ns::Int,
                                qfunction::Function, qp::Function,
                                xgrid::Vector{Float64}, sgrid::Vector{Float64},
                                fgrid::Matrix{Float64}, unc::BitArray,
                                xswts::Vector{Float64},
                                R::Float64, γ::Float64,
                                ell::Vector{Float64}, μ::Vector{Float64},
                                η::Float64, T::Float64, ω::Float64, H::Float64, ee::Matrix{Float64})
    nxns = nx*ns
    bigΨ    = zeros(nxns,nxns) # this is also dF2_dM
    smallψ = zeros(nxns,nxns)
    dF2_dRZ = zeros(nxns)
    dF2_dELL = zeros(nxns,nxns)
    dF2_dM = zeros(nxns,nxns)
    dF2_dWH = zeros(nxns)
    dF2_dTT = zeros(nxns)

    for isp=1:ns
        for iap=1:nx
            ip  = nx*(isp-1)+iap
            sumWH = 0.
            sumRZ = 0.
            sumTT = 0.
            for iss=1:ns
                for ia=1:nx
                    i = nx*(iss-1)+ia
                    bigΨ[ip,i] = xswts[i]*μ[i]*qfunction(ee[i,ip])*fgrid[iss,isp]/(ω*H*sgrid[isp])
                    smallψ[ip,i] = xswts[i]*μ[i]*qp(ee[i,ip])*fgrid[iss,isp]/((ω*H*sgrid[isp])^2)
                    sumWH += bigΨ[ip,i] + smallψ[ip,i]*ee[i,ip]/(ω*H*sgrid[isp])
                    sumRZ += smallψ[ip,i]*(R*exp(-γ))*max(xgrid[ia] - 1/ell[i],-η)
                    dF2_dELL[ip,i] = smallψ[ip,i]*(R*exp(-γ))*(unc[i]/ell[i])
                    sumTT += smallψ[ip,i]*T
                    dF2_dM[ip,i] = xswts[i]*qfunction(ee[i,ip])*fgrid[iss,isp]/(ω*H*sgrid[isp]) # note, now we linearize
                end
            end
            dF2_dWH[ip] = -sumWH
            dF2_dRZ[ip] = -sumRZ
            dF2_dTT[ip] = -sumTT
        end
    end
    return dF2_dWH, dF2_dRZ, dF2_dTT, dF2_dELL, bigΨ, dF2_dM
end

function normalize(m::HetDSGESimpleTaylor, JJ::Matrix{Float64})

    Qx, _, Qleft, Qright = compose_normalization_matrices(m)

    m <= Setting(:n_predetermined_variables, size(Qx, 1))

	Jac1 = Qleft*sparse(JJ)*Qright

    return Jac1
end

function compose_normalization_matrices(m::HetDSGESimpleTaylor)
    nx::Int = get_setting(m, :nx)
    ns::Int = get_setting(m, :ns)
    nscalars::Int = get_setting(m, :nscalars)
    nyscalars::Int = get_setting(m, :nyscalars)
    nxscalars::Int = get_setting(m, :nxscalars)

    # Create PPP matrix
    P1 = kron(Matrix{Float64}(I, ns,ns),ones(nx,1))
    Ptemp = Matrix{Float64}(I, nx, nx)
    Ptemp = Ptemp[:, 2:end]
    P2 = kron(Matrix{Float64}(I, ns, ns), Ptemp)
    P  = hcat(P1, P2)

    Q,R = qr(P)
    Q = Array(Q)
    S         = Q[:, ns+1:end]'

    nxns = nx*ns

    Qleft     = sparse(cat(Matrix{Float64}(I, nxns, nxns),S,
                    Matrix{Float64}(I, nxns, nxns),S,
                    Matrix{Float64}(I, nscalars, nscalars), dims = [1 2]))
    Qx        = sparse(cat(Matrix{Float64}(I, nxns, nxns),S,
                    Matrix{Float64}(I, nxscalars, nxscalars), dims = [1 2]))
    Qy        = sparse(cat(Matrix{Float64}(I, nxns, nxns),S,
                    Matrix{Float64}(I, nyscalars, nyscalars), dims = [1 2]))
    Qright    = sparse(cat(Qx',Qy',Qx',Qy', dims = [1,2]))

    return Qx, Qy, Qleft, Qright
end


text:
# test.jl test file

using PkgTemplates;
t = Template(; user="jipq6175",
               authors=["Yen"; "Cornell University"],
               plugins=[License(; name="MIT"); TravisCI(); Codecov(); Coveralls(); AppVeyor()]);

t("PkgDev.jl")


text:
@testset "test_writer" begin
    @testset "write latitude" begin
        lat = IGCLatitude(51.29971666666667)
        s_expected_lat = "5117983N"
        @test string(lat) == s_expected_lat

        lat = IGCLatitude(-33.94611666666667)
        s_expected_lat = "3356767S"
        @test string(lat) == s_expected_lat
    end

    @testset "write longitude" begin
        long = IGCLongitude(6.956383333333333)
        s_expected_long = "00657383E"
        @test string(long) == s_expected_long

        long = IGCLongitude(-99.71176666666666)
        s_expected_long = "09942706W"
        @test string(long) == s_expected_long
    end

    @testset "write time" begin
        t = Time(16, 2, 45)
        s_expected_time = "160245"
        @test string(t, IGCTimeFormat) == s_expected_time
    end

    @testset "write date" begin
        d = Date(2019, 8, 20)
        s_expected_date = "200819"
        @test string(d, IGCDateFormat) == s_expected_date
    end

    @testset "write pressure altitude" begin
        alt = IGCPressureAltitude(288)
        s_expected_alt = "00288"
        @test string(alt) == s_expected_alt
    end

    @testset "write GPS altitude" begin
        alt = IGCGpsAltitude(429)
        s_expected_alt = "00429"
        @test string(alt) == s_expected_alt
    end

    @testset "write FixValidity" begin
        fv = FixValidity.Fix3D
        io = IOBuffer()
        FixValidity.write(io, fv)
        @test String(take!(io)) == "A"
        @test FixValidity.string(fv) == "A"

        fv = FixValidity.Fix2D
        io = IOBuffer()
        FixValidity.write(io, fv)
        @test String(take!(io)) == "V"
        @test FixValidity.string(fv) == "V"
    end

    @testset "write A record" begin
        rec = A_record(
            "XXX",
            "ABC",
            "FLIGHT:1"
        )
        s_expected_rec = "AXXXABC FLIGHT:1"
        @test string(rec) == s_expected_rec
    end

    @testset "write B record" begin
        rec = B_record(
            Time(16, 2, 45),
            IGCLatitude(51.118766666666666),
            IGCLongitude(-1.8216666666666668),
            FixValidity.parse('A'),
            IGCPressureAltitude(288),
            IGCGpsAltitude(429),
            36,
            "19509020"
        )
        s_expected_rec = "B1602455107126N00149300WA002880042919509020"
        @test string(rec) == s_expected_rec
    end

    @testset "write C record task info" begin
        rec = C_record_task_info(
            Date(2001, 7, 15),  # declaration_date
            Time(21, 38, 41),  # declaration_time
            Date(2001, 7, 16),  # flight_date
            "0001",  # number
            2,  # num_turnpoints
            "500K Tri"  # description
        )
        s_expected_rec = "C150701213841160701000102 500K Tri"
        @test string(rec) == s_expected_rec
    end

    @testset "write C record waypoint info" begin
        rec = C_record_waypoint_info(
            IGCLatitude(51.18931666666667),  # latitude
            IGCLongitude(-1.03165),  # longitude
            "Lasham Clubhouse"  # description
        )
        s_expected_rec = "C5111359N00101899W Lasham Clubhouse"
        @test string(rec) == s_expected_rec
    end

    @testset "write GpsQualifier" begin
        qual = GpsQualifier.GPS
        io = IOBuffer()
        GpsQualifier.write(io, qual)
        @test String(take!(io)) == "1"

        qual = GpsQualifier.DGPS
        io = IOBuffer()
        GpsQualifier.write(io, qual)
        @test String(take!(io)) == "2"
    end

    @testset "write D record" begin
        rec = D_record(
            GpsQualifier.DGPS,  # qualifier
            "0331"  # station_id
        )
        s_expected_rec = "D20331"
        @test string(rec) == s_expected_rec

        rec = D_record(
            GpsQualifier.GPS,  # qualifier
            "9999"  # station_id
        )
        s_expected_rec = "D19999"
        @test string(rec) == s_expected_rec
    end

    @testset "write E record" begin
        rec = E_record(
            Time(16, 2, 45),  # time
            "PEV",  # tlc
            ""  # extension_string
        )
        s_expected_rec = "E160245PEV"
        @test string(rec) == s_expected_rec
    end

    @testset "write F record" begin
        rec = F_record(
            Time(16, 2, 40),  # time
            ["04", "06", "09", "12", "36", "24", "22", "18", "21"]  # satellites
        )
        s_expected_rec = "F160240040609123624221821"
        @test string(rec) == s_expected_rec
    end

    @testset "write G record" begin
        rec = G_record(
            "REJNGJERJKNJKRE31895478537H43982FJN9248F942389T433T"
        )
        s_expected_rec = "GREJNGJERJKNJKRE31895478537H43982FJN9248F942389T433T"
        @test string(rec) == s_expected_rec
    end

    @testset "write I record" begin
        rec = I_record([
            IGCExtension(36:38, "FXA"),
            IGCExtension(39:40, "SIU"),
            IGCExtension(41:43, "ENL")
        ])
        s_expected_rec = "I033638FXA3940SIU4143ENL"
        @test string(rec) == s_expected_rec
    end

    @testset "write J record" begin
        rec = J_record([
            IGCExtension(8:12, "HDT")
        ])
        s_expected_rec = "J010812HDT"
        @test string(rec) == s_expected_rec
    end

    @testset "write K record" begin
        rec = K_record(
            Time(16, 2, 48),
            "00090",
            7  # tofix: is it really necessary?
        )
        s_expected_rec = "K16024800090"
        @test string(rec) == s_expected_rec
    end

    @testset "write L record" begin
        rec = L_record(
            "XXX",
            "RURITANIAN STANDARD NATIONALS DAY 1"
        )
        s_expected_rec = "LXXXRURITANIAN STANDARD NATIONALS DAY 1"
        @test string(rec) == s_expected_rec
    end

    @testset "write H record fix accuracy" begin
        rec = H_record_FiXAccuracy(35)
        s_expected_rec = "HFFXA035"
        @test string(rec) == s_expected_rec
    end

    @testset "write H record utc date" begin
        rec = H_record_DaTE(2001, 7, 16)
        s_expected_rec = "HFDTE160701"
        @test string(rec) == s_expected_rec
    end

    @testset "write H record pilot" begin
        rec = H_record_PiLoT("Bloggs Bill D")
        s_expected_rec = "HFPLTPILOTINCHARGE:Bloggs Bill D"
        @test string(rec) == s_expected_rec
    end

    @testset "write H record copilot" begin
        rec = H_record_Copilot("Smith-Barry John A")
        s_expected_rec = "HFCM2CREW2:Smith-Barry John A"
        @test string(rec) == s_expected_rec
    end

    @testset "write H record glider type" begin
        rec = H_record_GliderType("Schleicher ASH-25")
        s_expected_rec = "HFGTYGLIDERTYPE:Schleicher ASH-25"
        @test string(rec) == s_expected_rec
    end

    @testset "write H record glider registration" begin
        rec = H_record_GliderRegistration("ABCD-1234")
        s_expected_rec = "HFGIDGLIDERID:ABCD-1234"
        @test string(rec) == s_expected_rec
    end

    @testset "write H record gps datum" begin
        rec = H_record_GpsDatum("WGS-1984")
        s_expected_rec = "HFDTM100GPSDATUM:WGS-1984"
        @test string(rec) == s_expected_rec
    end

    @testset "write H record firmware revision" begin
        rec = H_record_FirmwareRevision("6.4")
        s_expected_rec = "HFRFWFIRMWAREVERSION:6.4"
        @test string(rec) == s_expected_rec
    end

    @testset "write H record hardware revision" begin
        rec = H_record_HardwareRevision("3.0")
        s_expected_rec = "HFRHWHARDWAREVERSION:3.0"
        @test string(rec) == s_expected_rec
    end

    @testset "write H record manufacturer model" begin
        rec = H_record_ManufacturerModel("Manufacturer", "Model")
        s_expected_rec = "HFFTYFRTYPE:Manufacturer,Model"
        @test string(rec) == s_expected_rec
    end

    @testset "write H record gps receiver" begin
        rec = H_record_GpsReceiver(
            "MarconiCanada",
            "Superstar",
            12,
            MaxAlt(10000, "m")
        )
        s_expected_rec = "HFGPS:MarconiCanada,Superstar,12ch,max10000m"
        @test string(rec) == s_expected_rec
    end

    @testset "write H record pressure sensor" begin
        rec = H_record_PressureAltitudeSensor(
            "Sensyn",
            "XYZ1111",
            MaxAlt(11000, "m")
        )
        s_expected_rec = "HFPRSPRESSALTSENSOR:Sensyn,XYZ1111,max11000m"
        @test string(rec) == s_expected_rec
    end

    @testset "write H record competition id" begin
        rec = H_record_CompetitionId("XYZ-78910")
        s_expected_rec = "HFCIDCOMPETITIONID:XYZ-78910"
        @test string(rec) == s_expected_rec
    end

    @testset "write H record competition class" begin
        rec = H_record_CompetitionClass("15m Motor Glider")
        s_expected_rec = "HFCCLCOMPETITIONCLASS:15m Motor Glider"
        @test string(rec) == s_expected_rec
    end

    @testset "write H record time zone offset" begin
        rec = H_record_TimeZoneOffset(3 * 60)
        s_expected_rec = "HFTZNTIMEZONE:3"
        @test string(rec) == s_expected_rec

        rec = H_record_TimeZoneOffset(-(3 * 60 + 30))
        s_expected_rec = "HFTZNTIMEZONE:-3.30"
        @test string(rec) == s_expected_rec
    end

    @testset "write H record mop sensor" begin
        rec = H_record_MeansOfPropulsion("MOP-(SN:1,ET=1375,0,1375,0,3.05V,p=0),Ver:0")
        s_expected_rec = "HFMOPSENSOR:MOP-(SN:1,ET=1375,0,1375,0,3.05V,p=0),Ver:0"
        @test string(rec) == s_expected_rec
    end

    @testset "write H record site" begin
        rec = H_record_Site("lk15comp")
        s_expected_rec = "HFSITSite:lk15comp"
        @test string(rec) == s_expected_rec
    end

    @testset "write H record units of measure" begin
        rec = H_record_UnitsOfMeasure(["km", "ft", "kt"])
        s_expected_rec = "HFUNTUnits:km,ft,kt"
        @test string(rec) == s_expected_rec
    end

    @testset "write IGCDocument" begin
        io = IOBuffer()  # this buffer is filled at each update! call

        igcdoc = IGCDocument(stream = io, parsing_mode = IGC.ParsingMode.STRICT)

        update!(igcdoc, A_record("XXX", "ABC", "FLIGHT:1"))
        s_expected = """AXXXABC FLIGHT:1"""

        update!(igcdoc, H_record_FiXAccuracy('F', 35))
        update!(igcdoc, H_record_DaTE('F', Date(2001, 07, 16)))
        update!(igcdoc, H_record_PiLoT('F', "Bloggs Bill D"))
        update!(igcdoc, H_record_Copilot('F', "Smith-Barry John A"))
        update!(igcdoc, H_record_GliderType('F', "Schleicher ASH-25"))
        update!(igcdoc, H_record_GliderRegistration('F', "ABCD-1234"))
        update!(igcdoc, H_record_GpsDatum('F', "WGS-1984"))
        update!(igcdoc, H_record_FirmwareRevision('F', "6.4"))
        update!(igcdoc, H_record_HardwareRevision('F', "3.0"))
        update!(igcdoc, H_record_ManufacturerModel('F', "Manufacturer", "Model"))
        update!(igcdoc, H_record_GpsReceiver('F', "MarconiCanada", "Superstar", 12, MaxAlt(10000, "m")))
        update!(igcdoc, H_record_PressureAltitudeSensor('F', "Sensyn", "XYZ1111", MaxAlt(11000, "m")))
        update!(igcdoc, H_record_CompetitionId('F', "XYZ-78910"))
        update!(igcdoc, H_record_CompetitionClass('F', "15m Motor Glider"))
        s_expected = s_expected * EOL_DEFAULT * """HFFXA035
HFDTE160701
HFPLTPILOTINCHARGE:Bloggs Bill D
HFCM2CREW2:Smith-Barry John A
HFGTYGLIDERTYPE:Schleicher ASH-25
HFGIDGLIDERID:ABCD-1234
HFDTM100GPSDATUM:WGS-1984
HFRFWFIRMWAREVERSION:6.4
HFRHWHARDWAREVERSION:3.0
HFFTYFRTYPE:Manufacturer,Model
HFGPS:MarconiCanada,Superstar,12ch,max10000m
HFPRSPRESSALTSENSOR:Sensyn,XYZ1111,max11000m
HFCIDCOMPETITIONID:XYZ-78910
HFCCLCOMPETITIONCLASS:15m Motor Glider"""

        rec = I_record([IGCExtension(36:38, "FXA"), IGCExtension(39:40, "SIU"), IGCExtension(41:43, "ENL")])
        update!(igcdoc, rec)
        s_expected = s_expected * EOL_DEFAULT * "I033638FXA3940SIU4143ENL"

        rec = J_record([IGCExtension(8:12, "HDT")])
        update!(igcdoc, rec)
        s_expected = s_expected * EOL_DEFAULT * "J010812HDT"

        rec = C_record_task_info(Date(2001, 7, 15), Time(21, 38, 41), Date(2001, 7, 16), "0001", 2, "500K Tri")
        update!(igcdoc, rec)
        s_expected = s_expected * EOL_DEFAULT * """C150701213841160701000102 500K Tri"""

        records = [
            C_record_waypoint_info(IGCLatitude(51.18931666666667), IGCLongitude(-1.03165), "Lasham Clubhouse"),
            C_record_waypoint_info(IGCLatitude(51.16965), IGCLongitude(-1.0440666666666667), "Lasham Start S, Start"),
            C_record_waypoint_info(IGCLatitude(52.15153333333333), IGCLongitude(-2.9204499999999998), "Sarnesfield, TP1"),
            C_record_waypoint_info(IGCLatitude(52.50245), IGCLongitude(-0.2935333333333333), "Norman Cross, TP2"),
            C_record_waypoint_info(IGCLatitude(51.16965), IGCLongitude(-1.0440666666666667), "Lasham Start S, Finish"),
            C_record_waypoint_info(IGCLatitude(51.18931666666667), IGCLongitude(-1.03165), "Lasham Clubhouse")
        ]
        update!(igcdoc, records)
        s_expected = s_expected * EOL_DEFAULT * """C5111359N00101899W Lasham Clubhouse
C5110179N00102644W Lasham Start S, Start
C5209092N00255227W Sarnesfield, TP1
C5230147N00017612W Norman Cross, TP2
C5110179N00102644W Lasham Start S, Finish
C5111359N00101899W Lasham Clubhouse"""

        rec = F_record(Time(16, 02, 40), ["04", "06", "09", "12", "36", "24", "22", "18", "21"])
        update!(igcdoc, rec)
        s_expected = s_expected * EOL_DEFAULT * "F160240040609123624221821"

        records = [
            B_record(Time(16, 02, 40), IGCLatitude(54.11868333333334), IGCLongitude(-2.8223666666666665), FixValidity.Fix3D, IGCPressureAltitude(280), IGCGpsAltitude(421), 36, "20509950"),
            D_record(GpsQualifier.DGPS, "0331"),
            E_record(Time(16, 02, 45), "PEV", ""),
            B_record(Time(16, 02, 45), IGCLatitude(51.118766666666666), IGCLongitude(-1.8216666666666668), FixValidity.Fix3D, IGCPressureAltitude(288), IGCGpsAltitude(429), 36, "19509020"),
            B_record(Time(16, 02, 50), IGCLatitude(51.1189), IGCLongitude(-1.8213833333333334), FixValidity.Fix3D, IGCPressureAltitude(290), IGCGpsAltitude(432), 36, "21009015"),
            B_record(Time(16, 02, 55), IGCLatitude(51.119), IGCLongitude(-1.82035), FixValidity.Fix3D, IGCPressureAltitude(290), IGCGpsAltitude(430), 36, "20009012"),
            F_record(Time(16, 03, 00), ["06", "09", "12", "36", "24", "22", "18", "21"]),
            B_record(Time(16, 03, 00), IGCLatitude(51.119166666666665), IGCLongitude(-1.8200333333333334), FixValidity.Fix3D, IGCPressureAltitude(291), IGCGpsAltitude(432), 36, "25608009"),
            E_record(Time(16, 03, 05), "PEV", ""),
            B_record(Time(16, 03, 05), IGCLatitude(51.11966666666667), IGCLongitude(-1.81975), FixValidity.Fix3D, IGCPressureAltitude(291), IGCGpsAltitude(435), 36, "21008015"),
            B_record(Time(16, 03, 10), IGCLatitude(51.1202), IGCLongitude(-1.8195666666666668), FixValidity.Fix3D, IGCPressureAltitude(293), IGCGpsAltitude(435), 36, "19608024"),
            K_record(Time(16, 02, 48), "00090", 7),
            B_record(Time(16, 02, 48), IGCLatitude(51.120333333333335), IGCLongitude(-1.8191666666666666), FixValidity.Fix3D, IGCPressureAltitude(494), IGCGpsAltitude(436), 36, "19008018"),
            B_record(Time(16, 02, 52), IGCLatitude(51.122166666666665), IGCLongitude(-1.8187833333333334), FixValidity.Fix3D, IGCPressureAltitude(496), IGCGpsAltitude(439), 36, "19508015")
        ]
        update!(igcdoc, records)
        s_expected = s_expected * EOL_DEFAULT * """B1602405407121N00249342WA002800042120509950
D20331
E160245PEV
B1602455107126N00149300WA002880042919509020
B1602505107134N00149283WA002900043221009015
B1602555107140N00149221WA002900043020009012
F1603000609123624221821
B1603005107150N00149202WA002910043225608009
E160305PEV
B1603055107180N00149185WA002910043521008015
B1603105107212N00149174WA002930043519608024
K16024800090
B1602485107220N00149150WA004940043619008018
B1602525107330N00149127WA004960043919508015"""

        records = [
            L_record("XXX", "RURITANIAN STANDARD NATIONALS DAY 1"),
            L_record("XXX", "FLIGHT TIME: 4:14:25, TASK SPEED:58.48KTS")
        ]
        update!(igcdoc, records)
        s_expected = s_expected * EOL_DEFAULT * """LXXXRURITANIAN STANDARD NATIONALS DAY 1
LXXXFLIGHT TIME: 4:14:25, TASK SPEED:58.48KTS"""

        records = [
            G_record("REJNGJERJKNJKRE31895478537H43982FJN9248F942389T433T"),
            G_record("JNJK2489IERGNV3089IVJE9GO398535J3894N358954983O0934"),
            G_record("SKTO5427FGTNUT5621WKTC6714FT8957FGMKJ134527FGTR6751"),
            G_record("K2489IERGNV3089IVJE39GO398535J3894N358954983FTGY546"),
            G_record("12560DJUWT28719GTAOL5628FGWNIST78154INWTOLP7815FITN")
        ]
        update!(igcdoc, records)
        s_expected = s_expected * EOL_DEFAULT * """GREJNGJERJKNJKRE31895478537H43982FJN9248F942389T433T
GJNJK2489IERGNV3089IVJE9GO398535J3894N358954983O0934
GSKTO5427FGTNUT5621WKTC6714FT8957FGMKJ134527FGTR6751
GK2489IERGNV3089IVJE39GO398535J3894N358954983FTGY546
G12560DJUWT28719GTAOL5628FGWNIST78154INWTOLP7815FITN"""
        @test String(take!(io)) == s_expected

        # this buffer is filled only at end of file
        io_at_end = IOBuffer()
        write(io_at_end, igcdoc)
        #@test String(take!(io_at_end)) == s_expected
    end

end