text:
##############################################################################
# Some utilities
##############################################################################

function add_to_index_maps!(forward_map::Dict{Int, Int},
                            backward_map::Dict{Int, Int},
                            inds::Array{Int},
                            offset::Int)
    for i in inds
        index = length(forward_map) + offset
        forward_map[i] = index
        backward_map[index] = i
    end
end

function add_to_index_maps!(forward_map::Dict{Int, Int},
                            inds::Array{Int},
                            offset::Int)
    for i in inds
        index = length(forward_map) + offset
        forward_map[i] = index
    end
end

##############################################################################
# The main structure definition
##############################################################################

mutable struct ReSHOPMathProgBaseModel <: MPB.AbstractMathProgModel
    options::Dict{String, Any}

    solver_name::String

    x_l::Vector{Float64}
    x_u::Vector{Float64}
    g_l::Vector{Float64}
    g_u::Vector{Float64}

    nvar::Int
    ncon::Int

    obj
    constrs::Vector{Any}

    lin_constrs::Vector{Dict{Int32, Float64}}
    lin_obj::Dict{Int32, Float64}

    r_codes::Vector{Int}
    j_counts::Vector{Int}

    vartypes::Vector{Symbol}
    varlinearities_con::Vector{Symbol}
    varlinearities_obj::Vector{Symbol}
    conlinearities::Vector{Symbol}
    objlinearity::Symbol

    v_index_map::Dict{Int, Int}
    v_index_map_rev::Dict{Int, Int}

    # was:
    # c_index_map::Dict{Int, Int}
    # c_index_map_rev::Dict{Int, Int}

    nonquad_idx::Dict{Int, Int}
    quad_idx::Dict{Int, Int}

    sense::Symbol

    x_0::Vector{Float64}

    objval::Float64
    solution::Vector{Float64}

    status::Symbol
    solve_exitcode::Int
    solve_result_num::Int
    solve_result::String
    model_result_num::Int
    model_result::String
    solve_message::String
    solve_time::Float64

    model_type::MODEL_TYPE

    quad_equs::Vector{Any}
#    quad_obj::Tuple{Vector{Int}, Vector{Int}, Vector{Float64}}
    quad_obj::Tuple
    offset::Int
    emp::Union{Nothing,Function}

    d::Union{Nothing,MPB.AbstractNLPEvaluator}

    reshop_ctx::Ptr{context}
    reshop_ctx_dest::Ptr{context}
    reshop_mdl::Ptr{reshop_model}
    reshop_mdl_solver::Ptr{reshop_model}
    reshop_options::Ptr{reshop_options}
    gams_dir::String

    function ReSHOPMathProgBaseModel(solver_name::String,
                                options::Dict{String,Any},
                                model_type::MODEL_TYPE,
                                emp)
        o = new(options,
            solver_name,
            zeros(0),
            zeros(0),
            zeros(0),
            zeros(0),
            0,
            0,
            :(0),
            [],
            Dict{Int, Float64}[],
            Dict{Int, Float64}(),
            Int[],
            Int[],
            Symbol[],
            Symbol[],
            Symbol[],
            Symbol[],
            :Lin,
            Dict{Int, Int}(),
            Dict{Int, Int}(),
            Dict{Int, Int}(),
            Dict{Int, Int}(),
            :Min,
            zeros(0),
            NaN,
            zeros(0),
            :NotSolved,
            -1,
            -1,
            "?",
            -1,
            "?",
            "",
            NaN,
            model_type,
            Vector{Any}(),
            (),
            0,
            emp,
            nothing,
            Ptr{context}(C_NULL),
            Ptr{context}(C_NULL),
            Ptr{reshop_model}(C_NULL),
            Ptr{reshop_model}(C_NULL),
            Ptr{reshop_options}(C_NULL),
            "")
        @compat finalizer(reshop_cleanup, o)
        o
    end
end

include(joinpath("MBP_wrapper", "utils.jl"))

struct ReSHOPLinearQuadraticModel <: MPB.AbstractLinearQuadraticModel
    inner::ReSHOPMathProgBaseModel
end

struct ReSHOPNonlinearModel <: MPB.AbstractNonlinearModel
    inner::ReSHOPMathProgBaseModel
end

struct ReSHOPConicModel <: MPB.AbstractNonlinearModel
    inner::ReSHOPMathProgBaseModel
end

MPB.NonlinearModel(s::ReSHOPSolver) = ReSHOPNonlinearModel(
    ReSHOPMathProgBaseModel(s.solver_name, s.options, nlp, s.emp)
)

MPB.LinearQuadraticModel(s::ReSHOPSolver) = ReSHOPLinearQuadraticModel(
    ReSHOPMathProgBaseModel(s.solver_name, s.options, qcp, s.emp)
)

function MPB.ConicModel(s::ReSHOPSolver)
    error("ConicModel is not yet supported")
end

function MPB.loadproblem!(outer::ReSHOPNonlinearModel, nvar::Integer, ncon::Integer,
                      x_l, x_u, g_l, g_u, sense::Symbol,
                      d::MPB.AbstractNLPEvaluator)
    m = outer.inner

    m.nvar, m.ncon = nvar, ncon
    loadcommon!(m, x_l, x_u, g_l, g_u, sense)

    m.d = d

    MPB.initialize(m.d.value, [:ExprGraph])

    # Process constraints
    m.constrs = map(1:m.ncon) do i
        c = MPB.constr_expr(m.d.value, i)

        # Remove relations and bounds from constraint expressions
        if length(c.args) == 3
            if VERSION < v"0.5-"
                expected_head = :comparison
                expr_index = 1
                rel_index = 2
            else
                expected_head = :call
                expr_index = 2
                rel_index = 1
            end

            @assert c.head == expected_head
            # Single relation constraint: expr rel bound
            rel = c.args[rel_index]
            m.r_codes[i] = relation_to_reshop[rel]
            if rel == [:<=, :(==)]
                m.g_u[i] = c.args[3]
            end
            if rel in [:>=, :(==)]
                m.g_l[i] = c.args[3]
            end
            c = c.args[expr_index]
        else
            # Double relation constraint: bound <= expr <= bound
            @assert c.head == :comparison
            m.r_codes[i] = relation_to_reshop[:multiple]
            m.g_u[i] = c.args[5]
            m.g_l[i] = c.args[1]
            c = c.args[3]
        end

        # Convert non-linear expression to non-linear, linear and constant
        c, constant, m.conlinearities[i] = process_expression!(
            c, m.lin_constrs[i], m.varlinearities_con)

        # Update bounds on constraint
        m.g_l[i] -= constant
        m.g_u[i] -= constant

        # Update jacobian counts using the linear constraint variables
        for j in keys(m.lin_constrs[i])
            m.j_counts[j] += 1
        end
        c
    end

    # Process objective
    m.obj = MPB.obj_expr(m.d.value)

    if length(m.obj.args) < 2
        m.obj = 0
    else
        # Convert non-linear expression to non-linear, linear and constant
        m.obj, constant, m.objlinearity = process_expression!(
            m.obj, m.lin_obj, m.varlinearities_obj)

        # Add constant back into non-linear expression
        if constant != 0
            m.obj = add_constant(m.obj, constant)
        end
    end
    m
end

function MPB.loadproblem!(outer::ReSHOPLinearQuadraticModel, A::AbstractMatrix,
                      x_l, x_u, c, g_l, g_u, sense)
    m = outer.inner
    m.ncon, m.nvar = size(A)

    loadcommon!(m, x_l, x_u, g_l, g_u, sense)

    # Load A into the linear constraints
    @assert (m.ncon, m.nvar) == size(A)
    load_A!(m, A)
    m.constrs = zeros(m.ncon)  # Dummy constraint expression trees

    # Load c
    for (index, val) in enumerate(c)
        m.lin_obj[index] = val
    end
    # TODO(xhub) see if we can get rid of that
    m.obj = 0  # Dummy objective expression tree

    # Process variables bounds
    for j = 1:m.ncon
        lower = m.g_l[j]
        upper = m.g_u[j]
        if lower == -Inf
            if upper == Inf
                error("Neither lower nor upper bound on constraint $j")
            else # <=
                m.r_codes[j] = 2
            end
        else
            if lower == upper  # ==
                m.r_codes[j] = 4
            elseif upper == Inf # >=
                m.r_codes[j] = 1
            else # lb <= expr <= ub
                m.r_codes[j] = -1
            end
        end
    end
    m
end

function load_A!(m::ReSHOPMathProgBaseModel, A::SparseMatrixCSC{Float64})
    for var = 1:A.n, k = A.colptr[var] : (A.colptr[var + 1] - 1)
        m.lin_constrs[A.rowval[k]][var] = A.nzval[k]
        m.j_counts[var] += 1
    end
end

function load_A!(m::ReSHOPMathProgBaseModel, A::Matrix{Float64})
    for con = 1:m.ncon, var = 1:m.nvar
        val = A[con, var]
        if val != 0
            m.lin_constrs[A.rowval[k]][var] = A.nzval[k]
            m.j_counts[var] += 1
        end
    end
end

function loadcommon!(m::ReSHOPMathProgBaseModel, x_l, x_u, g_l, g_u, sense)
    m.x_l, m.x_u = x_l, x_u
    m.g_l, m.g_u = g_l, g_u
    setsense!(m, sense)

    m.lin_constrs = [Dict{Int, Float64}() for _ in 1:m.ncon]
    m.j_counts = zeros(Int, m.nvar)

    m.r_codes = Vector{Int}(undef, m.ncon)

    m.varlinearities_con = fill(:Lin, m.nvar)
    m.varlinearities_obj = fill(:Lin, m.nvar)
    m.conlinearities = fill(:Lin, m.ncon)
    m.objlinearity = :Lin

    m.vartypes = fill(:Cont, m.nvar)
    m.x_0 = zeros(m.nvar)
end

getvartype(m::ReSHOPMathProgBaseModel) = copy(m.vartypes)
function setvartype!(m::ReSHOPMathProgBaseModel, cat::Vector{Symbol})
    @assert all(x-> (x in [:Cont,:Bin,:Int,:external]), cat)
    m.vartypes = copy(cat)
end

getsense(m::ReSHOPMathProgBaseModel) = m.sense
function setsense!(m::ReSHOPMathProgBaseModel, sense::Symbol)
    @assert sense == :Min || sense == :Max
    m.sense = sense
end

setwarmstart!(m::ReSHOPMathProgBaseModel, v::Vector{Float64}) = m.x_0 = v

function MPB.addquadconstr!(m::ReSHOPLinearQuadraticModel, linearidx, linearval, quadrowidx, quadcolidx, quadval, sense, rhs)
   # we have to do a little translation of the sense here ...
   push!(m.inner.quad_equs, tuple(linearidx,linearval,quadrowidx,quadcolidx,quadval, quad_relation_sense[sense], rhs))
end

function MPB.setquadobj!(m::ReSHOPLinearQuadraticModel, rowidx, colidx, quadval)
   m.inner.quad_obj = tuple(rowidx, colidx, quadval)
end

function make_var_index!(m::ReSHOPMathProgBaseModel)
    m.v_index_map = Dict(zip(1:m.nvar, 0:(m.nvar-1)))
    m.v_index_map_rev = Dict(zip(0:(m.nvar-1), 1:m.nvar))
end

#function make_var_index!(m::ReSHOPMathProgBaseModel)
#    nonlin_cont = Int[]
#    nonlin_int = Int[]
#    lin_cont = Int[]
#    lin_int = Int[]
#    lin_bin = Int[]
#
#    # TODO(xhub) we do that multiple times in the EMP context
#    for i in 1:length(m.vartypes)
#        if m.varlinearities_obj[i] == :Nonlin ||
#           m.varlinearities_con[i] == :Nonlin
#            if m.vartypes[i] == :Cont || m.vartypes[i] == :external
#                push!(nonlin_cont, i)
#            else
#                push!(nonlin_int, i)
#            end
#        else
#           if m.vartypes[i] == :Cont || m.vartypes[i] == :external
#                push!(lin_cont, i)
#            elseif m.vartypes[i] == :Int
#                push!(lin_int, i)
#            else
#                push!(lin_bin, i)
#            end
#        end
#    end
#
#    # Index variables in required order
#    for var_list in (nonlin_cont, nonlin_int, lin_cont, lin_bin, lin_int)
#        add_to_index_maps!(m.v_index_map, m.v_index_map_rev, var_list, 0)
#    end
#    CONFIG[:debug] && println("DEBUG: $(m.v_index_map)")
#end

function make_con_index!(m::ReSHOPMathProgBaseModel)
    nonlin_cons = Int[]
    lin_cons = Int[]

    for i in 1:m.ncon
        if m.conlinearities[i] == :Nonlin
            push!(nonlin_cons, i)
        else
            push!(lin_cons, i)
        end
    end
    for con_list in (nonlin_cons, lin_cons)
        add_to_index_maps!(m.nonquad_idx, con_list, 1)
    end

    if length(m.quad_idx) == 0
       m.quad_idx = Dict(enumerate(Vector{Int}(Compat.range(1+m.ncon, length=length(m.quad_equs)))))
    end

    CONFIG[:debug] && println("DEBUG: make_con_index: nonquad_idx: $(m.nonquad_idx)\n quad_idx: $(m.quad_idx)")
end

function optimize!(m::ReSHOPMathProgBaseModel)
    m.status = :NotSolved
    m.solve_exitcode = -1
    m.solve_result_num = -1
    m.solve_result = "?"
    m.model_result_num = -1
    m.model_result = "?"
    m.solve_message = ""

    # There is no non-linear binary type, only non-linear discrete, so make
    # sure binary vars have bounds in [0, 1]
    for i in 1:m.nvar
        if m.vartypes[i] == :Bin
            if m.x_l[i] < 0
                m.x_l[i] = 0
            end
            if m.x_u[i] > 1
                m.x_u[i] = 1
            end
        end
    end

    m.reshop_options = reshop_options_set(m.options)

    # TODO(Xhub) this hack has to go.
    if m.emp.hasvalue
       return m.emp.value()
    end

    make_var_index!(m)
    make_con_index!(m)

    # Run solver and save exitcode
    t = time()
    m.reshop_ctx = create_reshop_ctx(m)
    reshop_set_modeltype(m)
    # Solve via gams for now
    m.reshop_ctx_dest, m.gams_dir = reshop_setup_gams()

    m.reshop_mdl = reshop_alloc(m.reshop_ctx)
    m.reshop_mdl_solver = reshop_alloc(m.reshop_ctx_dest)
    m.solve_exitcode = reshop_solve(m.reshop_mdl, m.reshop_mdl_solver, m.reshop_ctx_dest, m.solver_name)

#    ccall((:print_model, libreshop), Cint, (Ptr{context},), m.reshop_ctx)
    m.solve_time = time() - t

    if m.solve_exitcode == 0
        report_results(m)
    else
        println("ReSHOP: solver failed with status $(m.solve_exitcode)")
        m.status = :Error
        m.solution = fill(NaN, m.nvar)
        m.solve_result = "failure"
        m.solve_result_num = 999
    end
end

function getconstrduals(m::ReSHOPMathProgBaseModel)

    ctx = m.reshop_ctx

    x = fill(NaN, numconstr(m))

   if has_objective(m)
        offset = m.offset
    else
        offset = m.offset - 1
    end

    for idx in 1:m.ncon
       eidx = m.nonquad_idx[idx] + offset
       x[idx] = ctx_getequmult(ctx, eidx)
    end

    return x
end

function getquadconstrduals(quadm::ReSHOPLinearQuadraticModel)

    m = quadm.inner
    ctx = m.reshop_ctx

    x = fill(NaN, numconstr(m))

    if has_objective(m)
        offset = m.offset
    else
        offset = m.offset - 1
    end

    for (idx, equ) in enumerate(m.quad_equs)
       eidx = m.quad_idx[idx] + offset
       x[idx] = ctx_getequmult(ctx, eidx)
    end

    return x
end

function getreducedcosts(m::ReSHOPMathProgBaseModel)

    ctx = m.reshop_ctx

    x = fill(NaN, numvar(m))

    for idx in 1:numvar(m)
       x[idx] = ctx_getvarmult(ctx, idx-1)
    end

    return x
end

function report_results_common(m::ReSHOPMathProgBaseModel)
    x = fill(NaN, m.nvar)
    m.objval = NaN

    for index in 0:(m.nvar - 1)
        i = m.v_index_map_rev[index]
        x[i] = ctx_getvarval(m.reshop_ctx, index)
    end

    m.solution = x

    ###########################################################################
    # Convert solve_result
    #
    # GAMS return two information:
    #  - the solve status
    #  - the model status
    #
    # - :Optimal
    # - :Infeasible
    # - :Unbounded
    # - :UserLimit (iteration limit or timeout)
    # - :Error (and maybe others)
    ###########################################################################

    m.model_result_num = rhp_get_modelstat(m.reshop_ctx)
    m.solve_result_num = rhp_get_solvestat(m.reshop_ctx)

    m.model_result = unsafe_string(ccall((:ctx_getmodelstattxt, libreshop), Cstring, (Ptr{context}, Cint), m.reshop_ctx_dest, m.model_result_num))
    m.solve_result = unsafe_string(ccall((:ctx_getsolvestattxt, libreshop), Cstring, (Ptr{context}, Cint), m.reshop_ctx, m.solve_result_num))

    # GAMS already uses an 1-indices
    solver_code = solver_stat[m.solve_result_num]
    model_code = model_stat[m.model_result_num]

    CONFIG[:debug] && println("solver stat $(m.solve_result) ($(m.solve_result_num)); model stat $(m.model_result) ($(m.model_result_num))")

    if solver_code == :Optimal
       if model_code == :OptimalGlobal || model_code == :OptimalLocal || model_code == :Integer
          m.status = :Optimal
       elseif model_code == :Unbounded || model_code == :UnboundedNoSolution
          m.status = :Unbounded
       elseif model_code == :InfeasibleGlobal || model_code == :InfeasibleLocal || model_code == :InfeasibleIntermed || model_code == :InfeasibleNoSolution
          m.status = :Infeasible
       elseif model_code == :Feasible
          # TODO investigate that. Baron is weird
          m.status = :Optimal
       elseif model_code == :NoSolutionReturned
          gams_solver = ctx_get_solvername(m.reshop_ctx_dest)
          m.status = :Optimal
#          if gams_solver == "jams" || gams_solver == "JAMS"
             # This is fine, we have a kludge in the code
#          else
#             println("ReSHOP: Solve successed, but no solution was returned by solver $(gams_solver)!")
#             m.status = :Error
#          end
       else
          println("ReSHOP: unhandle case: solver stat $(m.solve_result); model stat $(m.model_result)")
          m.status = :Error
       end
    elseif solver_code == :Iteration || solver_code == :Resource
       m.status == :UserLimit
    elseif solver_code == :License || model_code == :LicenseError
       println("ReSHOP: License error. Check that you have a valid license")
       m.status == :Error
    elseif solver_code == :Capability
       if m.solver_name == ""
          sname = default
       else
          sname = m.solver_name
       end

       println("ReSHOP: solver $(sname) cannot solve the specific problem")
       m.status = :Error

    else
       println("ReSHOP: solver stat is $(m.solve_result) and model stat is $(m.model_result)")
       m.status = :Error
    end

    CONFIG[:debug] && println("status is $(m.status)")
 end

getsolvername(s::ReSHOPSolver) = basename(s.solver_name)

MPB.getreducedcosts(nlpm::ReSHOPNonlinearModel) = getreducedcosts(nlpm.inner)
MPB.getreducedcosts(quadm::ReSHOPLinearQuadraticModel) = getreducedcosts(quadm.inner)
MPB.getconstrduals(nlpm::ReSHOPNonlinearModel) = getconstrduals(nlpm.inner)
MPB.getconstrduals(quadm::ReSHOPLinearQuadraticModel) = getconstrduals(quadm.inner)

# Wrapper functions
status(m::ReSHOPMathProgBaseModel) = m.status
getsolution(m::ReSHOPMathProgBaseModel) = copy(m.solution)
getobjval(m::ReSHOPMathProgBaseModel) = m.objval
numvar(m::ReSHOPMathProgBaseModel) = m.nvar
numconstr(m::ReSHOPMathProgBaseModel) = m.ncon + length(m.quad_equs)
getsolvetime(m::ReSHOPMathProgBaseModel) = m.solve_time

# Access to solve results
get_solve_result(m::ReSHOPMathProgBaseModel) = m.solve_result
get_solve_result_num(m::ReSHOPMathProgBaseModel) = m.solve_result_num
get_model_result(m::ReSHOPMathProgBaseModel) = m.model_result
get_model_result_num(m::ReSHOPMathProgBaseModel) = m.model_result_num
get_solve_message(m::ReSHOPMathProgBaseModel) = m.solve_message
get_solve_exitcode(m::ReSHOPMathProgBaseModel) = m.solve_exitcode

for f in [:getvartype,:getsense,:optimize!,:status,:getsolution,:getobjval,:numvar,:numconstr,:getsolvetime]
    @eval MPB.$f(m::ReSHOPNonlinearModel) = $f(m.inner)
    @eval MPB.$f(m::ReSHOPLinearQuadraticModel) = $f(m.inner)
end

for f in [:get_solve_result,:get_solve_result_num,:get_solve_message,:get_solve_exitcode]
    @eval $f(m::ReSHOPNonlinearModel) = $f(m.inner)
    @eval $f(m::ReSHOPLinearQuadraticModel) = $f(m.inner)
end
for f in [:setvartype!,:setsense!,:setwarmstart!]
    @eval MPB.$f(m::ReSHOPNonlinearModel, x) = $f(m.inner, x)
    @eval MPB.$f(m::ReSHOPLinearQuadraticModel, x) = $f(m.inner, x)
end

# Deallocate the data
function reshop_cleanup(o::ReSHOPMathProgBaseModel)
    ctx_dealloc(o.reshop_ctx)
    ctx_dealloc(o.reshop_ctx_dest)
    reshop_options_dealloc(o.reshop_options)
    reshop_free(o.reshop_mdl)
    reshop_free(o.reshop_mdl_solver)
    if (!isempty(o.gams_dir))
        try
            rm(o.gams_dir, recursive=true, force=true)
        catch
            iswin && run(`cmd /C RMDIR /s /q $(o.gams_dir)`)
        end
    end
end


text:
using MiscUtilities
using Test

@testset "AST" begin include("ast.jl") end
@testset "ilog2" begin include("ilog2.jl") end

ext: jl


text:
using ITensors
using ITensorParallel
using Random

include(joinpath(pkgdir(ITensors), "examples", "src", "electronk.jl"))
include(joinpath(pkgdir(ITensors), "examples", "src", "hubbard.jl"))

function main(;
  Nx::Int=6,
  Ny::Int=3,
  U::Float64=4.0,
  t::Float64=1.0,
  maxdim::Int=3000,
  conserve_ky=true,
  use_splitblocks=true,
  seed=1234,
  in_partition=ITensorParallel.default_in_partition,
)
  Random.seed!(seed)
  @show Threads.nthreads()
  @show ITensors.using_threaded_blocksparse()

  N = Nx * Ny

  sweeps = Sweeps(10)
  maxdims = min.([100, 200, 400, 800, 2000, 3000, maxdim], maxdim)
  setmaxdim!(sweeps, maxdims...)
  setcutoff!(sweeps, 1e-6)
  setnoise!(sweeps, 1e-6, 1e-7, 1e-8, 0.0)
  @show sweeps

  sites = siteinds("ElecK", N; conserve_qns=true, conserve_ky=conserve_ky, modulus_ky=Ny)

  ℋ = hubbard(; Nx=Nx, Ny=Ny, t=t, U=U, ky=true)
  ℋs = opsum_sum(ℋ, Threads.nthreads(); in_partition=in_partition)

  H = Vector{MPO}(undef, Threads.nthreads())
  Threads.@threads for n in 1:length(ℋs)
    Hn = MPO(ℋs[n], sites)
    Hn = use_splitblocks ? splitblocks(linkinds, Hn) : Hn
    H[n] = Hn
  end

  @show maxlinkdim.(H)

  # Number of structural nonzero elements in a bulk
  # Hamiltonian MPO tensor
  @show nnz(H[1][end ÷ 2])
  @show nnzblocks(H[1][end ÷ 2])

  # Create start state
  state = Vector{String}(undef, N)
  for i in 1:N
    x = (i - 1) ÷ Ny
    y = (i - 1) % Ny
    if x % 2 == 0
      if y % 2 == 0
        state[i] = "Up"
      else
        state[i] = "Dn"
      end
    else
      if y % 2 == 0
        state[i] = "Dn"
      else
        state[i] = "Up"
      end
    end
  end

  psi0 = randomMPS(sites, state, 10)

  energy, psi = @time dmrg(ThreadedProjMPOSum(H), psi0, sweeps; svd_alg="divide_and_conquer")
  @show Nx, Ny
  @show t, U
  @show flux(psi)
  @show maxlinkdim(psi)
  @show energy
  return energy, H, psi
end

function custom_in_partition(sites::Tuple, p, nparts)
  return p == mod1(sites[1], nparts)
end

Random.seed!(1234)

ITensors.BLAS.set_num_threads(1)
ITensors.Strided.disable_threads()
ITensors.disable_threaded_blocksparse()

# A function that specifies which partition/processor that an
# MPO term will be on.
in_partition = ITensorParallel.default_in_partition # or: custom_in_partition
main(; in_partition=in_partition)

max_stars_repo_name: mtfishman/ITensorParallel

max_issues_repo_name: mtfishman/ITensorParallel

max_forks_repo_name: mtfishman/ITensorParallel


text:
using ImageFiltering, OffsetArrays, Base.Test

@testset "basic" begin
    @test eltype(KernelFactors.IIRGaussian(3)) == Float64
    @test eltype(KernelFactors.IIRGaussian(Float32, 3)) == Float32
    kern = KernelFactors.IIRGaussian([1,2.0f0])
    @test isa(kern, Tuple{KernelFactors.ReshapedOneD,KernelFactors.ReshapedOneD})
    @test isa(ImageFiltering.filter_algorithm([1],[1],kern), Algorithm.IIR)
    @test isa(ImageFiltering.filter_algorithm([1],[1],(kern...,Kernel.Laplacian())), Algorithm.Mixed)

    @test KernelFactors.kernelfactors(([0,3], [1,7]))  == (reshape([0,3], 1:2, 0:0), reshape([1,7], 0:0, 1:2))
    @test KernelFactors.kernelfactors(([0,3], [1,7]')) == (reshape([0,3], 2, 1), reshape([1,7], 1, 2))

    # Warnings
    const OLDERR = STDERR
    fname = tempname()
    open(fname, "w") do f
        redirect_stderr(f)
        try
            KernelFactors.IIRGaussian(0.5, emit_warning=false)
        finally
            redirect_stderr(OLDERR)
        end
    end
    @test isempty(chomp(readstring(fname)))
    open(fname, "w") do f
        redirect_stderr(f)
        try
            KernelFactors.IIRGaussian(0.5)
        finally
            redirect_stderr(OLDERR)
        end
    end
    @test contains(readstring(fname), "too small for accuracy")
end

nothing


text:
function smallF(anF, bigPhi, nn)

    (fRows, fCols) = size(anF)
    lilFL = anF[(nn + 1):fRows, (nn + 1):fRows]
    uu = nullspace(lilFL)
    theNull = size(uu, 2)
    #eOpts.disp = 0;
    lilFU = anF[1:nn, (nn + 1):fRows]
    onLeft = lilFU
    onRight = bigPhi[(nn + 1):fRows, :]
    inMiddle = lilFL

    return onLeft, inMiddle, onRight

end # function


text:
import LinearAlgebra.axpy!

function axpy!(a::DoubleFloat{T}, xv::StridedVector{DoubleFloat{T}}, yv::StridedVector{DoubleFloat{T}})  where {N, T <: AbstractFloat}
    _axpy!(a,xv,yv,Vec{Npref,T})
end
function axpy!(a::Complex{DoubleFloat{T}}, xv::StridedVector{Complex{DoubleFloat{T}}}, yv::StridedVector{Complex{DoubleFloat{T}}})  where {N, T <: AbstractFloat}
     n = length(xv)
    (length(yv) == n) || throw(ArgumentError("arguments must have equal lengths"))
    _axpy!(n,a,xv,1,yv,1,Vec{Npref,T})
end

function _axpy!(a::DoubleFloat{T}, xv::StridedVector{DoubleFloat{T}}, yv::StridedVector{DoubleFloat{T}},::Type{Vec{N,T}}) where {N, T <: AbstractFloat}
     n = length(xv)
    (length(yv) == n) || throw(ArgumentError("arguments must have equal lengths"))
    nd,nr = divrem(n, N)
    ahi = HI(a)
    alo = LO(a)

    shi = Vec{N,T}(ahi)    # this lets us avoid @generated (Thanks, SIMD.jl!)
    slo = Vec{N,T}(alo)
    @inbounds begin
        for i in 1:nd
            i0=(i-1)*N

            xhi = vgethi(xv,i0,Vec{N,T})
            xlo = vgetlo(xv,i0,Vec{N,T})
            zhi, zlo = dfvmul(xhi, xlo, shi, slo)

            yhi = vgethi(yv,i0,Vec{N,T})
            ylo = vgetlo(yv,i0,Vec{N,T})

            zhi, zlo = dfvadd(yhi, ylo, zhi, zlo)
            vputhilo!(yv,i0,zhi,zlo)

        end
    end
    (nr == 0) && return yv
    @inbounds begin
        @simd for i in (nd*N)+1:n
            yv[i] += a * xv[i]
        end
    end
    yv
end

# warning: unsafe indexing
function _axpy!(n, a::DoubleFloat{T},
                xv::StridedVecOrMat{DoubleFloat{T}}, ix1::Int,
                yv::StridedVecOrMat{DoubleFloat{T}}, iy1::Int,
                ::Type{Vec{N,T}}) where {N, T <: AbstractFloat}
    nd,nr = divrem(n, N)
    ahi = HI(a)
    alo = LO(a)
    shi = Vec{N,T}(ahi)
    slo = Vec{N,T}(alo)
    ixoff = ix1-1
    iyoff = iy1-1
    @inbounds begin
        for i in 1:nd
            i0=(i-1)*N
            ix0 = ixoff + i0
            iy0 = iyoff + i0

            xhi = vgethi(xv,ix0,Vec{N,T})
            xlo = vgetlo(xv,ix0,Vec{N,T})
            yhi = vgethi(yv,iy0,Vec{N,T})
            ylo = vgetlo(yv,iy0,Vec{N,T})

            zhi, zlo = dfvmul(xhi, xlo, shi, slo)
            zhi, zlo = dfvadd(yhi, ylo, zhi, zlo)
            vputhilo!(yv,iy0,zhi,zlo)

        end
    end
    (nr == 0) && return yv
    @inbounds begin
        @simd for i in (nd*N)+1:n
            yv[iyoff+i] += a * xv[ixoff+i]
        end
    end
    yv
end

function _axpy!(n, a::Complex{DoubleFloat{T}},
                xv::StridedVecOrMat{Complex{DoubleFloat{T}}}, ix1::Integer,
                yv::StridedVecOrMat{Complex{DoubleFloat{T}}}, iy1::Integer,
                ::Type{Vec{N,T}}) where {N, T <: AbstractFloat}
    nd,nr = divrem(n, N)
    arhi = HI(real(a))
    arlo = LO(real(a))
    srhi = Vec{N,T}(arhi)
    srlo = Vec{N,T}(arlo)
    aihi = HI(imag(a))
    ailo = LO(imag(a))
    sihi = Vec{N,T}(aihi)
    silo = Vec{N,T}(ailo)
    ixoff = ix1-1
    iyoff = iy1-1
    @inbounds begin
        for i in 1:nd
            i0=(i-1)*N
            ix0 = ixoff + i0
            iy0 = iyoff + i0

            xrhi = vgethire(xv,ix0,Vec{N,T})
            xrlo = vgetlore(xv,ix0,Vec{N,T})
            yrhi = vgethire(yv,iy0,Vec{N,T})
            yrlo = vgetlore(yv,iy0,Vec{N,T})
            xihi = vgethiim(xv,ix0,Vec{N,T})
            xilo = vgetloim(xv,ix0,Vec{N,T})
            yihi = vgethiim(yv,iy0,Vec{N,T})
            yilo = vgetloim(yv,iy0,Vec{N,T})

            z1hi, z1lo = dfvmul(xrhi, xrlo, srhi, srlo)
            z2hi, z2lo = dfvmul(xrhi, xrlo, sihi, silo)
            z3hi, z3lo = dfvmul(xihi, xilo, srhi, srlo)
            z4hi, z4lo = dfvmul(xihi, xilo, sihi, silo)

            zrhi, zrlo = dfvsub(z1hi, z1lo, z4hi, z4lo)
            zihi, zilo = dfvadd(z2hi, z2lo, z3hi, z3lo)

            zrhi, zrlo = dfvadd(yrhi, yrlo, zrhi, zrlo)
            zihi, zilo = dfvadd(yihi, yilo, zihi, zilo)
            vputhilo!(yv,iy0,zrhi,zrlo,zihi,zilo)

        end
    end
    (nr == 0) && return yv
    @inbounds begin
        @simd for i in (nd*N)+1:n
            yv[iyoff+i] += a * xv[ixoff+i]
        end
    end
    yv
end

function rmul!(xv::StridedVecOrMat{DT},a::DT) where {DT <: Union{DoubleFloat{T}, Complex{DoubleFloat{T}}}} where T
    _rmul!(xv,1,a,Vec{Npref,T})
end

# Pretend multiplication of these types is commutative.
function lmul!(a::DT, xv::StridedVecOrMat{DT}) where {DT <: Union{DoubleFloat{T}, Complex{DoubleFloat{T}}}} where T
    _rmul!(xv,1,a,Vec{Npref,T})
end

function _rmul!(xv::StridedVecOrMat{DoubleFloat{T}}, ix1::Integer,
                a::DoubleFloat{T},
                ::Type{Vec{N,T}}) where {N, T <: AbstractFloat}
    n = length(xv)
    ixoff = ix1-1
    nd,nr = divrem(n, N)
    ahi = HI(a)
    alo = LO(a)

    shi = Vec{N,T}(ahi)    # this lets us avoid @generated (Thanks, SIMD.jl!)
    slo = Vec{N,T}(alo)
    @inbounds begin
        for i in 1:nd
            i0=(i-1)*N
            ix0 = i0 + ixoff
            xhi = vgethi(xv,ix0,Vec{N,T})
            xlo = vgetlo(xv,ix0,Vec{N,T})
            zhi, zlo = dfvmul(xhi, xlo, shi, slo)
            vputhilo!(xv,ix0,zhi,zlo)
        end
    end
    (nr == 0) && return xv
    @inbounds begin
        @simd for i in (nd*N)+1:n
            xv[ixoff+i] = a * xv[ixoff+i]
        end
    end
    xv
end

function _rmul!(xv::StridedVecOrMat{Complex{DoubleFloat{T}}}, ix1::Integer,
                a::Complex{DoubleFloat{T}},
                ::Type{Vec{N,T}}) where {N, T <: AbstractFloat}
    n = length(xv)
    ixoff = ix1-1
    nd,nr = divrem(n, N)
    arhi = HI(real(a))
    arlo = LO(real(a))
    srhi = Vec{N,T}(arhi)
    srlo = Vec{N,T}(arlo)
    aihi = HI(imag(a))
    ailo = LO(imag(a))
    sihi = Vec{N,T}(aihi)
    silo = Vec{N,T}(ailo)
    @inbounds begin
        for i in 1:nd
            i0=(i-1)*N
            ix0 = ixoff + i0

            xrhi = vgethire(xv,ix0,Vec{N,T})
            xrlo = vgetlore(xv,ix0,Vec{N,T})
            xihi = vgethiim(xv,ix0,Vec{N,T})
            xilo = vgetloim(xv,ix0,Vec{N,T})

            z1hi, z1lo = dfvmul(xrhi, xrlo, srhi, srlo)
            z2hi, z2lo = dfvmul(xrhi, xrlo, sihi, silo)
            z3hi, z3lo = dfvmul(xihi, xilo, srhi, srlo)
            z4hi, z4lo = dfvmul(xihi, xilo, sihi, silo)

            zrhi, zrlo = dfvsub(z1hi, z1lo, z4hi, z4lo)
            zihi, zilo = dfvadd(z2hi, z2lo, z3hi, z3lo)
            vputhilo!(xv,ix0,zrhi,zrlo,zihi,zilo)

        end
    end
    (nr == 0) && return xv
    @inbounds begin
        @simd for i in (nd*N)+1:n
            xv[ixoff+i] = a * xv[ixoff+i]
        end
    end
    xv
end


text:
immutable BoxplotGeometry <: Gadfly.GeometryElement
end

const boxplot = BoxplotGeometry

element_aesthetics(::BoxplotGeometry) = [:x, :y, :color]

default_statistic(::BoxplotGeometry) = Gadfly.Stat.boxplot()

function render(geom::BoxplotGeometry, theme::Gadfly.Theme,
                aes::Gadfly.Aesthetics)
    Gadfly.assert_aesthetics_defined("Geom.bar", aes,
                                     :lower_fence, :lower_hinge, :middle,
                                     :upper_hinge, :upper_fence, :outliers)
    Gadfly.assert_aesthetics_equal_length("Geom.bar", aes,
                                     :lower_fence, :lower_hinge, :middle,
                                     :upper_hinge, :upper_fence, :outliers)

    default_aes = Gadfly.Aesthetics()
    default_aes.color = PooledDataArray(ColorValue[theme.default_color])
    default_aes.x = Float64[0.5]
    aes = inherit(aes, default_aes)

    aes_iter = zip(aes.lower_fence,
                   aes.lower_hinge,
                   aes.middle,
                   aes.upper_hinge,
                   aes.upper_fence,
                   aes.outliers,
                   cycle(aes.x),
                   cycle(aes.color.refs))

    forms = Compose.Form[]
    middle_forms = Compose.Form[]

    r = theme.default_point_size
    bw = 1.0cx - theme.boxplot_spacing

    # TODO: handle color non-nothing color

    for (lf, lh, mid, uh, uf, outliers, x, cref) in aes_iter
        c = aes.color.pool[cref]
        sc = theme.highlight_color(c) # stroke color
        mc = theme.middle_color(c) # middle color

        # Middle
        push!(middle_forms, compose(lines((x - 1/6, mid), (x + 1/6, mid)),
                                    linewidth(theme.middle_width), stroke(mc)))

        # Box
        push!(forms, compose(rectangle(x*cx - bw/2, lh, bw, uh - lh),
                            fill(c), stroke(sc),
                            linewidth(theme.highlight_width)))

        # Whiskers
        push!(forms, compose(lines((x, lh), (x, lf)),
                            linewidth(theme.line_width), stroke(sc)))

        push!(forms, compose(lines((x, uh), (x, uf)),
                            linewidth(theme.line_width), stroke(sc)))

        # Fences
        push!(forms, compose(lines((x - 1/6, lf), (x + 1/6, lf)),
                            linewidth(theme.line_width), stroke(sc)))

        push!(forms, compose(lines((x - 1/6, uf), (x + 1/6, uf)),
                            linewidth(theme.line_width), stroke(sc)))

        # Outliers
        if !isempty(outliers)
            push!(forms, compose(combine([circle(x, y, r) for y in outliers]...),
                                fill(c), stroke(sc)))
        end
    end

    compose(canvas(units_inherited=true),
            (canvas(units_inherited=true), combine(forms...)),
            (canvas(units_inherited=true, order=1), combine(middle_forms...)),
            svgclass("geometry"))
end


text:
using DataStructures
using Test
using Aqua
using ParallelUtilities
using Documenter
using OffsetArrays
import ParallelUtilities: pval, value, BinaryTreeNode, BranchChannel,
ProductSplit, SegmentedOrderedBinaryTree
import ParallelUtilities.ClusterQueryUtils: chooseworkers

@testset "Project quality" begin
    if VERSION < v"1.6.0"
        Aqua.test_all(ParallelUtilities, ambiguities=false)
    else
        Aqua.test_all(ParallelUtilities)
    end
end

DocMeta.setdocmeta!(ParallelUtilities, :DocTestSetup, :(using ParallelUtilities); recursive=true)

@testset "doctest" begin
    doctest(ParallelUtilities, manual = false)
end

@testset "pval" begin
    p1 = pval{Float64}(1, false, 2.0)
    p2 = pval{Int}(1, false, 2)

    @test pval{Float64}(p1) === p1
    @test pval{Int}(p1) === p2

    @test value(p1) === 2.0
    @test value(p2) === 2
    @test value(2) === 2
    @test value(nothing) === nothing
end

@testset "chooseworkers" begin
    workers = 1:8
    workers_on_hosts = OrderedDict("host1" => 1:4, "host2" => 5:8)
    @test chooseworkers(workers, 3, workers_on_hosts) == 1:3
    @test chooseworkers(workers, 5, workers_on_hosts) == 1:5

    workers_on_hosts = OrderedDict(Libc.gethostname() => 1:4, "host2" => 5:8)
    @test chooseworkers(workers, 3, workers_on_hosts) == 1:3
    @test chooseworkers(workers, 5, workers_on_hosts) == 1:5

    workers_on_hosts = OrderedDict("host1" => 1:4, Libc.gethostname() => 5:8)
    @test chooseworkers(workers, 3, workers_on_hosts) == 5:7
    @test chooseworkers(workers, 5, workers_on_hosts) == [5:8; 1]
end

@testset "Reduction functions" begin
    # BroadcastStack with OffsetArrays
    @testset "BroadcastStack" begin
        arr = ParallelUtilities.BroadcastStack(+, 1)(ones(2:4), ones(3:5))
        @test arr == OffsetArray([1, 2, 2, 1], 2:5)

        arr = ParallelUtilities.BroadcastStack(+, 1:2)(ones(1:2, 2:4), ones(2:3, 3:5))
        arr_exp = OffsetArray([1.0  1.0  1.0  0.0
                                1.0  2.0  2.0  1.0
                                0.0  1.0  1.0  1.0], 1:3, 2:5)
        @test arr == arr_exp
    end

    @testset "BroadcastFunction" begin
        x = ones(3); y = ones(3);
        b = ParallelUtilities.BroadcastFunction{1}(+)
        @test b(x, y) == ones(3) * 2
        @test x == ones(3) * 2
        @test y == ones(3)

        b = ParallelUtilities.BroadcastFunction{2}(+)
        x = ones(3); y = ones(3);
        @test b(x, y) == ones(3) * 2
        @test x == ones(3)
        @test y == ones(3) * 2
    end

    @testset "Flip" begin
        x = ones(3); y = ones(3);
        f = ParallelUtilities.Flip(ParallelUtilities.elementwisesum!)
        @test f(x,y) == ones(3) * 2
        @test x == ones(3)
        @test y == ones(3) * 2

        x = ones(3); y = ones(3);
        f = ParallelUtilities.Flip(ParallelUtilities.broadcastinplace(+, Val(2)))
        @test f(x,y) == ones(3) * 2
        @test x == ones(3) * 2
        @test y == ones(3)
    end
end

@testset "show" begin

    @testset "ProductSplit" begin
        io = IOBuffer()
        ps = ProductSplit((1:20, 1:30), 4, 1)
        show(io, ps)
        showstr = String(take!(io))
        startstr = string(length(ps))*"-element ProductSplit"
        @test startswith(showstr, startstr)
    end

    @testset "error" begin
        io = IOBuffer()

        showerror(io, ParallelUtilities.TaskNotPresentError((1:4,), (5,)))
        strexp = "could not find the task $((5,)) in the list $((1:4,))"
        @test String(take!(io)) == strexp
    end;

    @testset "BranchChannel" begin
        io = IOBuffer()

        b = BranchChannel(1, 0)
        show(io, b)
        strexp = "Leaf  : 1 ← 1"
        @test String(take!(io)) == strexp

        b = BranchChannel(1, 1)
        show(io, b)
        strexp = "Branch: 1 ← 1 ← 1 child"
        @test String(take!(io)) == strexp

        b = BranchChannel(1, 2)
        show(io, b)
        strexp = "Branch: 1 ← 1 ⇇ 2 children"
        @test String(take!(io)) == strexp
    end;

    @testset "BinaryTreeNode" begin
        io = IOBuffer()
        b = BinaryTreeNode(2, 3, 1)
        show(io, b)
        strexp = "BinaryTreeNode(p = 2, parent = 3, nchildren = 1)"
        @test String(take!(io)) == strexp
    end;

    @testset "BinaryTree" begin
        # check that show is working
        io = IOBuffer()
        tree = SegmentedOrderedBinaryTree(1:8, OrderedDict("host1" => 1:4, "host2" => 5:8))
        show(io, tree)
        show(io, ParallelUtilities.toptree(tree))
        show(io, ParallelUtilities.toptree(tree).tree)
    end
end;


text:
# what is the path?
_BASE_PATH = pwd()
_PATH_TO_SRC = joinpath(_BASE_PATH, "src")
_PATH_TO_DATA = joinpath(_BASE_PATH, "data")
_PATH_TO_MODEL = joinpath(_BASE_PATH, "model")
_PATH_TO_TMP = joinpath(_BASE_PATH, "tmp")

# load external packages -
using CSV
using DataFrames
using LinearAlgebra
using Plots
using Distributions
using DifferentialEquations
using Tables
using NumericalIntegration
using Optim

# load my codes -
include(joinpath(_PATH_TO_SRC, "Balances.jl"))
include(joinpath(_PATH_TO_SRC, "Kinetics.jl"))
include(joinpath(_PATH_TO_SRC, "Factory.jl"))
include(joinpath(_PATH_TO_SRC, "Utility.jl"))
include(joinpath(_PATH_TO_SRC, "Simulation.jl"))
include(joinpath(_PATH_TO_SRC, "Compute.jl"))
include(joinpath(_PATH_TO_SRC, "Learn.jl"))

max_stars_repo_name: varnerlab/BST-S-Sytem-Model-Coagulation

max_issues_repo_name: varnerlab/BST-S-Sytem-Model-Coagulation

max_forks_repo_name: varnerlab/BST-S-Sytem-Model-Coagulation


text:
# ---
# title: White Additive Gaussian Noise
# id: wagn-demo
# cover: assets/awgn_cover.png
# description: Add white additive gaussian noise to an image
# date: 2019-12-22
# author: Johnny Chen
# ---

# Mathematically, adding a white additive gaussian noise to an image is as simple
# as $Y = X + N$ where $N \sim \mathcal{N}(\mu,\,\sigma^{2})\,.$

using ImageNoise
using TestImages, ImageShow, ImageCore, ImageQualityIndexes, ImageTransformations
using FileIO, Random #src

# Let's add noise to gray image first

gray_img = float.(imresize(testimage("cameraman"), ratio=0.5))
noise = AdditiveWhiteGaussianNoise(0.1)
Random.seed!(0) #src
gray_noisy_img = apply_noise(gray_img, noise)

# For RGB image, noise aren't added to it channel by channel, instead, we generate
# a $3\times M\times N$ gaussian noise and directly added to its channelview result.
# Adding noise channel by channel would indeed get a different distribution (but
# still a gaussian noise).

rgb_img = float.(imresize(testimage("mandrill"), ratio=0.5))
Random.seed!(0) #src
rgb_noisy_img = apply_noise(rgb_img, noise)

# Colorful images of other formats are converted to RGB first since RGB color
# space is considered "linear".

hsv_img = HSV.(rgb_img)
Random.seed!(0) #src
hsv_noisy_img = apply_noise(hsv_img, noise)
incorrect_hsv_noisy_img = colorview(HSV, apply_noise(channelview(hsv_img), noise))

## If we compare their psnr, the second one is significantly lower than the first one
[assess_psnr(rgb_noisy_img, RGB.(hsv_noisy_img))
assess_psnr(rgb_noisy_img, RGB.(incorrect_hsv_noisy_img))]

mkpath("assets") #src
save("assets/awgn_cover.png", rgb_noisy_img) #src


text:
using ValueSymbols
using Base.Test

syms = [:car, :boat, :plane]
@test !isbits(syms[1])

@test syms[1] === :car
@test syms[2] !== :car
@test syms[3] !== syms[2]
@test syms[1] == :car
@test syms[2] != :car
@test syms[3] != syms[2]

@test syms[1] <= syms[1]
@test syms[1] > syms[2]
@test syms[1] < syms[3]

vsyms = map(ValueSymbol, syms)
@test isbits(vsyms[1])

io = IOBuffer()
print(io, vsyms[1])
@test takebuf_string(io) == ":car"

@test vsyms[1] === ValueSymbol(:car)
@test vsyms[2] !== ValueSymbol(:car)
@test vsyms[3] !== vsyms[2]
@test vsyms[1] == ValueSymbol(:car)
@test vsyms[2] != ValueSymbol(:car)
@test vsyms[3] != vsyms[2]

# ValueSymbol and Symbol are different types, so comparing them via
# `===` is always false
@test vsyms[1] !== :car
@test vsyms[2] !== :car
@test vsyms[3] !== vsyms[2]

@test vsyms[1] == :car
@test vsyms[2] != :car
@test vsyms[3] != vsyms[2]

@test vsyms[1] <= vsyms[1]
@test vsyms[1] > vsyms[2]
@test vsyms[1] < vsyms[3]

@test vsyms[1] <= syms[1]
@test vsyms[1] > syms[2]
@test vsyms[1] <vsyms[3]

io = IOBuffer()
serialize(io, vsyms)
vsyms2 = deserialize(IOBuffer(takebuf_array(io)))
for n in 1:length(vsyms2)
    @test vsyms2[n] === vsyms[n]
end


text:
# StridedView
struct StridedView{T,N,A<:DenseArray,F<:Union{FN,FC,FA,FT}} <: AbstractStridedView{T,N,F}
    parent::A
    size::NTuple{N,Int}
    strides::NTuple{N,Int}
    offset::Int
    op::F
end
function StridedView(parent::Array{S},
                        size::NTuple{N,Int} = size(parent),
                        strides::NTuple{N,Int} = strides(parent),
                        offset::Int = 0,
                        op::F = identity) where {S, N, F}

    T = Base.promote_op(op, S)
    # reshape array to vector in order to reduce number of element types
    StridedView{T,N,Vector{S},F}(reshape(parent, length(parent)), size,
                                    _normalizestrides(size, strides), offset, op)
end
function StridedView(parent::A,
                        size::NTuple{N,Int} = size(parent),
                        strides::NTuple{N,Int} = strides(parent),
                        offset::Int = 0,
                        op::F = identity) where {A<:DenseArray, N, F}

    T = Base.promote_op(op, eltype(parent))
    StridedView{T,N,A,F}(parent, size, strides, offset, op)
end
StridedView(a::StridedView) = a

StridedView(a::Adjoint) = StridedView(a')'
StridedView(a::Transpose) = transpose(StridedView(transpose(a)))
StridedView(a::Base.SubArray) = sview(StridedView(a.parent), a.indices...)
StridedView(a::Base.ReshapedArray) = sreshape(StridedView(a.parent), a.dims)

# Methods for StridedView
Base.size(a::StridedView) = a.size
Base.strides(a::StridedView) = a.strides
Base.stride(a::StridedView{<:Any, 0}, n::Int) = 1
Base.stride(a::StridedView{<:Any, N}, n::Int) where N =
    (n <= N) ? a.strides[n] : a.strides[N]*a.size[N]
offset(a::StridedView) = a.offset
Base.parent(a::StridedView) = a.parent

function blasstrides(a::StridedView{<:Any,2})
    # canonialize strides to make compatible with gemm
    if size(a, 2) == 1 && stride(a, 1) == 1
        return StridedView(a.parent, a.size, (1, size(a,1)), a.offset, a.op)
    else
        return a
    end
end

Base.similar(a::StridedView, ::Type{T}, dims::NTuple{N,Int}) where {N,T} =
    StridedView(similar(a.parent, T, dims))
Base.copy(a::StridedView) = copyto!(similar(a), a)

Base.unsafe_convert(::Type{Ptr{T}}, a::StridedView{T}) where T =
    pointer(a.parent, a.offset+1)

Base.dataids(a::StridedView) = Base.dataids(a.parent)

Base.IndexStyle(::Type{<:StridedView}) = Base.IndexCartesian()

# Indexing with N integer arguments
@inline function Base.getindex(a::StridedView{<:Any,N}, I::Vararg{Int,N}) where N
    @boundscheck checkbounds(a, I...)
    @inbounds r = a.op(a.parent[a.offset+_computeind(I, a.strides)])
    return r
end
@inline function Base.setindex!(a::StridedView{<:Any,N}, v, I::Vararg{Int,N}) where N
    @boundscheck checkbounds(a, I...)
    @inbounds a.parent[a.offset+_computeind(I, a.strides)] = a.op(v)
    return a
end

# force inlining so that view typically does not need to be created
@inline function Base.getindex(a::StridedView{<:Any,N},
        I::Vararg{SliceIndex,N}) where N
    StridedView(a.parent, _computeviewsize(a.size, I), _computeviewstrides(a.strides, I),
                a.offset + _computeviewoffset(a.strides, I), a.op)
end

@propagate_inbounds Base.getindex(a::StridedView, I::ParentIndex) =
    a.op(getindex(a.parent, I.i))
@propagate_inbounds Base.setindex!(a::StridedView, v, I::ParentIndex) =
    (setindex!(a.parent, a.op(v), I.i); return a)

# Specialized methods for `StridedView` which produce views/share data
Base.conj(a::StridedView{<:Real}) = a
Base.conj(a::StridedView) = StridedView(a.parent, a.size, a.strides, a.offset, _conj(a.op))

@inline function Base.permutedims(a::StridedView{<:Any,N}, p) where {N}
    _isperm(N, p) || throw(ArgumentError("Invalid permutation of length $N: $p"))
    newsize = TupleTools._permute(a.size, p)
    newstrides = TupleTools._permute(a.strides, p)
    return StridedView(a.parent, newsize, newstrides, a.offset, a.op)
end
_isperm(N::Integer, p::AbstractVector) = (length(p) == N && isperm(p))
_isperm(N::Integer, p::NTuple{M,Integer}) where M = (M == N && TupleTools.isperm(p))
_isperm(N::Integer, p) = false

LinearAlgebra.transpose(a::StridedView{<:Number,2}) = permutedims(a, (2,1))
LinearAlgebra.adjoint(a::StridedView{<:Number,2}) = permutedims(conj(a), (2,1))
LinearAlgebra.adjoint(a::StridedView{<:Any,2}) = # act recursively, like Base
    permutedims(StridedView(a.parent, a.size, a.strides, a.offset, _adjoint(a.op)), (2,1))
LinearAlgebra.transpose(a::StridedView{<:Any,2}) = # act recursively, like Base
    permutedims(StridedView(a.parent, a.size, a.strides, a.offset, _transpose(a.op)), (2,1))

Base.map(::FC, a::StridedView{<:Real}) = a
Base.map(::FT, a::StridedView{<:Number}) = a
Base.map(::FA, a::StridedView{<:Number}) = conj(a)
Base.map(::FC, a::StridedView) =
    StridedView(a.parent, a.size, a.strides, a.offset, _conj(a.op))
Base.map(::FT, a::StridedView) =
    StridedView(a.parent, a.size, a.strides, a.offset, _transpose(a.op))
Base.map(::FA, a::StridedView) =
    StridedView(a.parent, a.size, a.strides, a.offset, _adjoint(a.op))

@inline function sreshape(a::StridedView, newsize::Dims)
    if any(isequal(0), newsize)
        any(isequal(0), size(a)) || throw(DimensionMismatch())
        newstrides = one.(newsize)
    else
        newstrides = _computereshapestrides(newsize, _simplify(size(a), strides(a))...)
    end
    StridedView(a.parent, newsize, newstrides, a.offset, a.op)
end


text:
"""
Power nonlinearity, with energy C ∫ρ^α.
"""
struct PowerNonlinearity
    C::Real
    α::Real
end
(P::PowerNonlinearity)(basis) = TermPowerNonlinearity(basis, P.C, P.α)

struct TermPowerNonlinearity <: Term
    basis::PlaneWaveBasis
    C::Real
    α::Real
end

function ene_ops(term::TermPowerNonlinearity, ψ, occ; ρ, kwargs...)
    basis = term.basis
    dVol = basis.model.unit_cell_volume / prod(basis.fft_size)

    E = term.C * sum(ρ.real .^ term.α) * dVol
    potential = @. term.C * term.α * ρ.real^(term.α-1)

    ops = [RealSpaceMultiplication(basis, kpoint, potential) for kpoint in basis.kpoints]
    (E=E, ops=ops)
end
_kernel(C, α, ρ) = @. C * α * (α-1) * ρ^(α-2)
function compute_kernel(term::TermPowerNonlinearity; ρ::RealFourierArray, kwargs...)
    k = _kernel(term.C, term.α, ρ.real)
    Diagonal(vec(k))
end
function apply_kernel(term::TermPowerNonlinearity, dρ; ρ::RealFourierArray, kwargs...)
    k = _kernel(term.C, term.α, ρ.real)
    k .* dρ
end


text:
# ==============================================================================
# This file fixes a some issues found in Julia.
# ==============================================================================

# ==============================================================================
# https://github.com/JuliaLang/julia/issues/33332
# Fix bug in Sparse Matrix in Julia. They used vals1[j1] != 0 instead of using
# iszero(vals1[j1]) idem for vals2[j2]

function Base.:(==)(A1::SparseMatrixCSC, A2::SparseMatrixCSC)
    size(A1) != size(A2) && return false
    vals1, vals2 = nonzeros(A1), nonzeros(A2)
    rows1, rows2 = rowvals(A1), rowvals(A2)
    m, n = size(A1)
    @inbounds for i = 1:n
        nz1,nz2 = nzrange(A1,i), nzrange(A2,i)
        j1,j2 = first(nz1), first(nz2)
        # step through the rows of both matrices at once:
        while j1 <= last(nz1) && j2 <= last(nz2)
            r1,r2 = rows1[j1], rows2[j2]
            if r1==r2
                vals1[j1]!=vals2[j2] && return false
                j1+=1
                j2+=1
            else
                if r1<r2
                    !iszero(vals1[j1]) && return false
                    j1+=1
                else
                    !iszero(vals2[j2]) && return false
                    j2+=1
                end
            end
        end
        # finish off any left-overs:
        for j = j1:last(nz1)
            !iszero(vals1[j]) && return false
        end
        for j = j2:last(nz2)
            !iszero(vals2[j]) && return false
        end
    end
    return true
end

# ==============================================================================
# https://github.com/JuliaLang/julia/issues/33332 fixed in Julia 1.3
# Fix bug of sparse[i,j] = mp0 confusing with sparse[i,j] = 0

getcolptr(S::SparseMatrixCSC) = getfield(S, :colptr)

function _insert!(v::Vector, pos::Integer, item, nz::Integer)
    if nz > length(v)
        insert!(v, pos, item)
    else # nz < length(v)
        Base.unsafe_copyto!(v, pos+1, v, pos, nz - pos)
        v[pos] = item
        v
    end
end

function _setindex_scalar!(A::SparseMatrixCSC{MP{Tv},Ti}, _v, _i::Integer, _j::Integer) where {Tv,Ti<:Integer}
    v = convert(MP{Tv}, _v)
    i = convert(Ti, _i)
    j = convert(Ti, _j)
    if !((1 <= i <= size(A, 1)) & (1 <= j <= size(A, 2)))
        throw(BoundsError(A, (i,j)))
    end
    coljfirstk = Int(getcolptr(A)[j])
    coljlastk = Int(getcolptr(A)[j+1] - 1)
    searchk = searchsortedfirst(rowvals(A), i, coljfirstk, coljlastk, Base.Order.Forward)
    if searchk <= coljlastk && rowvals(A)[searchk] == i
        nonzeros(A)[searchk] = v
        return A
    end
    if !iszero(v)
        nz = getcolptr(A)[size(A, 2)+1]
        !isbitstype(Ti) || nz < typemax(Ti) ||
            throw(ArgumentError("nnz(A) going to exceed typemax(Ti) = $(typemax(Ti))"))
        _insert!(rowvals(A), searchk, i, nz)
        _insert!(nonzeros(A), searchk, v, nz)
        @simd for m in (j + 1):(size(A, 2) + 1)
            @inbounds getcolptr(A)[m] += Ti(1)
        end
    end
    return A
end

@inline Base.:setindex!(A::SparseMatrixCSC{MP{Tv},Ti}, _v, _i::Integer, _j::Integer) where {Tv,Ti<:Integer} =
    _setindex_scalar!(A, _v, _i, _j)

# Used by SimpleWeightedDiGraph
function SparseMatrixCSC{MP{Tv},Ti}(M::StridedMatrix) where {Tv,Ti}
    nz = count(!iszero, M)
    colptr = zeros(Ti, size(M, 2) + 1)
    nzval = Vector{MP{Tv}}(undef, nz)
    rowval = Vector{Ti}(undef, nz)
    colptr[1] = 1
    cnt = 1
    @inbounds for j in 1:size(M, 2)
        for i in 1:size(M, 1)
            v = M[i, j]
            if !iszero(v)
                rowval[cnt] = i
                nzval[cnt] = v
                cnt += 1
            end
        end
        colptr[j+1] = cnt
    end
    return SparseMatrixCSC(size(M, 1), size(M, 2), colptr, rowval, nzval)
end

# ==============================================================================
# Bug https://github.com/JuliaLang/julia/issues/33036

"""
    mpI

Fix an algebra conception in Julia official LinearAlgebra (uniformscaling.jl)
`I` representes an identity matrix of any size is defined with a booleen instead
of the function `one()`. As consequence, in Julia 0.4 the `eye(T,m,n)` could
created a max-plus identity matrix. Since Julia 0.7 `eye()` has been deprecated
and replaced by the buggy `Matrix{T}(I, m, n)`. This function uses `zero()` but
not `one()` and as consequence the max-plus identity matrix is not well formed.

This const allows to be more algebra compliant by calling `one()` and fixing the
fucntion `Matrix{T}(I, m, n)`.
"""
const global mpI = UniformScaling(one(MP{Float64}).λ)

# ==============================================================================
# Because Julia will create the ill-formed identity matrix mixing zero() and true
# instead of zero() and one()

@inline Base.literal_pow(::typeof(^), A::ArrMP{T}, ::Val{0}) where T =
    mpeye(T, size(A,1), size(A,2))


text:
module RawMutex

const shim = joinpath(dirname(@__FILE__), "..", "deps", "libshimuv.so")
Libdl.dlopen(shim, Libdl.RTLD_GLOBAL)

const UV_BARRIER_SIZE = ccall(:sizeof_uv_barrier_t, Cint, ())
const UV_MUTEX_SIZE = ccall(:jl_sizeof_uv_mutex, Cint, ())

typealias Mutex Ptr{Void}
typealias Barrier Ptr{Void}

function create_mutex()
  m = Libc.malloc(UV_MUTEX_SIZE)
  ccall(:uv_mutex_init, Void, (Ptr{Void},), m)
  return m
end

function create_barrier(count)
  b = Libc.malloc(UV_BARRIER_SIZE)
  err = ccall(:uv_barrier_init, Cint, (Ptr{Void}, Cuint), b, count)
  err != 0 && error("uv error code", err)
  return b
end

function close_mutex(m::Mutex)
    if m != C_NULL
        ccall(:uv_mutex_destroy, Void, (Ptr{Void},), m)
        Libc.free(m)
        nothing
    end
end

function close_barrier(b::Barrier)
  if b != C_NULL
    ccall(:uv_barrier_destroy, Void, (Ptr{Void},), b)
    Libc.free(b)
    nothing
  end
end

function trylock(m::Mutex)
    r = ccall(:uv_mutex_trylock, Cint, (Ptr{Void},), m)
    return r == 0
end

function unlock(m::Mutex)
    ccall(:uv_mutex_unlock, Void, (Ptr{Void},), m)
    return
end

function lock(m::Mutex)
  ccall(:uv_mutex_lock, Void, (Ptr{Void},), m)
end

function wait(b::Barrier)
  ccall(:uv_barrier_wait, Cint, (Ptr{Void},), b)
  nothing
end

end # module


text:
@add_problem lp function lp_dual_abs_atom(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    x = Variable()
    p = minimize(abs(x), x <= -1; numeric_type = T)

    if test
        @test vexity(p) == ConvexVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 1 atol = atol rtol = rtol
        @test evaluate(abs(x)) ≈ 1 atol = atol rtol = rtol
        @test p.constraints[1].dual ≈ 1 atol = atol rtol = rtol
    end

    x = Variable(2, 2)
    p = minimize(
        sum(abs(x)),
        x[2, 2] >= 1,
        x[1, 1] >= 1,
        x >= 0;
        numeric_type = T,
    )

    if test
        @test vexity(p) == ConvexVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 2 atol = atol rtol = rtol
        @test evaluate(sum(abs(x))) ≈ 2 atol = atol rtol = rtol
        @test p.constraints[1].dual ≈ 1 atol = atol rtol = rtol
        @test p.constraints[2].dual ≈ 1 atol = atol rtol = rtol
        @test p.constraints[3].dual[1, 1] ≈ 0 atol = atol rtol = rtol
        @test p.constraints[3].dual[2, 2] ≈ 0 atol = atol rtol = rtol
        @test p.constraints[3].dual[1, 2] <= 1 + atol
        @test p.constraints[3].dual[2, 1] <= 1 + atol
        @test p.constraints[3].dual[1, 2] >= -atol
        @test p.constraints[3].dual[2, 1] >= -atol
    end
end

@add_problem lp function lp_maximum_atom(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    x = Variable(10)
    a = shuffle(collect(0.1:0.1:1.0))
    p = minimize(maximum(x), x >= a; numeric_type = T)

    if test
        @test vexity(p) == ConvexVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ maximum(a) atol = atol rtol = rtol
        @test evaluate(maximum(x)) ≈ maximum(a) atol = atol rtol = rtol
    end
end

@add_problem lp function lp_minimum_atom(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    x = Variable(1)
    a = reshape(shuffle(collect(0.01:0.01:1.0)), (10, 10))
    p = maximize(minimum(x), x <= a; numeric_type = T)

    if test
        @test vexity(p) == ConvexVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ minimum(a) atol = atol rtol = rtol
        @test evaluate(minimum(x)) ≈ minimum(a) atol = atol rtol = rtol
    end

    x = Variable(4, 4)
    y = Variable(4, 6)
    z = Variable(1)
    c = ones(4, 1)
    d = fill(2.0, (6, 1))
    constraints = [[x y] <= 2, z <= 0, z <= x, 2z >= -1]
    objective = sum(x + z) + minimum(y) + c' * y * d
    p = maximize(objective, constraints; numeric_type = T)

    if test
        @test vexity(p) == ConvexVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 130 atol = atol rtol = rtol
        @test (evaluate(objective))[1] ≈ 130 atol = atol rtol = rtol
    end
end

@add_problem lp function lp_max_atom(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    x = Variable(10, 10)
    y = Variable(10, 10)
    a = reshape(shuffle(collect(0.01:0.01:1.0)), (10, 10))
    b = reshape(shuffle(collect(0.01:0.01:1.0)), (10, 10))
    p = minimize(maximum(max(x, y)), [x >= a, y >= b]; numeric_type = T)

    if test
        @test vexity(p) == ConvexVexity()
    end
    handle_problem!(p)
    max_a = maximum(a)
    max_b = maximum(b)
    if test
        @test p.optval ≈ max(max_a, max_b) atol = 10atol atol = atol rtol = rtol
        @test evaluate(maximum(max(x, y))) ≈ max(max_a, max_b) atol = 10atol atol =
            atol rtol = rtol
    end
end

@add_problem lp function lp_min_atom(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    x = Variable(10, 10)
    y = Variable(10, 10)
    a = reshape(shuffle(collect(0.01:0.01:1.0)), (10, 10))
    b = reshape(shuffle(collect(0.01:0.01:1.0)), (10, 10))
    p = maximize(minimum(min(x, y)), [x <= a, y <= b]; numeric_type = T)

    if test
        @test vexity(p) == ConvexVexity()
    end
    handle_problem!(p)
    min_a = minimum(a)
    min_b = minimum(b)
    if test
        @test p.optval ≈ min(min_a, min_b) atol = 10atol atol = atol rtol = rtol
        @test evaluate(minimum(min(x, y))) ≈ min(min_a, min_b) atol = 10atol atol =
            atol rtol = rtol
    end
end

@add_problem lp function lp_pos_atom(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    x = Variable(3)
    a = [-2; 1; 2]
    p = minimize(sum(pos(x)), [x >= a, x <= 2]; numeric_type = T)

    if test
        @test vexity(p) == ConvexVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 3 atol = atol rtol = rtol
        @test evaluate(sum(pos(x))) ≈ 3 atol = atol rtol = rtol
    end
end

@add_problem lp function lp_neg_atom(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    x = Variable(3)
    p = minimize(1, [x >= -2, x <= -2, neg(x) <= 3]; numeric_type = T)

    if test
        @test vexity(p) == ConvexVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 1 atol = atol rtol = rtol
        @test evaluate(sum(neg(x))) ≈ 6 atol = atol rtol = rtol
    end
end

@add_problem lp function lp_sumlargest_atom(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    x = Variable(2)
    p = minimize(
        sumlargest(x, 2) + sumlargest(x, 0),
        x >= [1; 1];
        numeric_type = T,
    )

    if test
        @test vexity(p) == ConvexVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 2 atol = atol rtol = rtol
        @test evaluate(sumlargest(x, 2)) ≈ 2 atol = atol rtol = rtol
    end

    x = Variable(4, 4)
    p = minimize(
        sumlargest(x, 3),
        x >= eye(4),
        x[1, 1] >= 1.5,
        x[2, 3] >= 2.1;
        numeric_type = T,
    )

    if test
        @test vexity(p) == ConvexVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 4.6 atol = atol rtol = rtol
        @test evaluate(sumlargest(x, 2)) ≈ 3.6 atol = atol rtol = rtol
    end
end

@add_problem lp function lp_sumsmallest_atom(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    x = Variable(4, 4)
    p = minimize(
        sumlargest(x, 2) + sumsmallest(x, 0),
        sumsmallest(x, 4) >= 1;
        numeric_type = T,
    )

    if test
        @test vexity(p) == ConvexVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 0.5 atol = atol rtol = rtol
        @test evaluate(sumsmallest(x, 4)) ≈ 1 atol = atol rtol = rtol
    end

    x = Variable(3, 2)
    p = maximize(
        sumsmallest(x, 3),
        x >= 2,
        x <= 5,
        sumlargest(x, 3) <= 12;
        numeric_type = T,
    )

    if test
        @test vexity(p) == ConvexVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 12 atol = atol rtol = rtol
        @test evaluate(sumsmallest(x, 3)) ≈ 12 atol = atol rtol = rtol
    end
end

@add_problem lp function lp_dotsort_atom(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    x = Variable(4, 1)
    p = minimize(
        dotsort(x, [1, 2, 3, 4]),
        sum(x) >= 7,
        x >= 0,
        x <= 2,
        x[4] <= 1;
        numeric_type = T,
    )

    if test
        @test vexity(p) == ConvexVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 19 atol = atol rtol = rtol
        @test vec(evaluate(x)) ≈ [2; 2; 2; 1] atol = atol rtol = rtol
        @test evaluate(dotsort(x, [1, 2, 3, 4])) ≈ 19 atol = atol rtol = rtol
    end

    x = Variable(2, 2)
    p = minimize(
        dotsort(x, [1 2; 3 4]),
        sum(x) >= 7,
        x >= 0,
        x <= 2,
        x[2, 2] <= 1;
        numeric_type = T,
    )

    if test
        @test vexity(p) == ConvexVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 19 atol = atol rtol = rtol
        @test evaluate(dotsort(x, [1, 2, 3, 4])) ≈ 19 atol = atol rtol = rtol
    end
end

@add_problem lp function lp_hinge_loss_atom(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    # TODO: @davidlizeng. We should finish this someday.
end

@add_problem lp function lp_dual_norm_inf_atom(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    x = Variable(3)
    p = minimize(norm_inf(x), [-2 <= x, x <= 1]; numeric_type = T)

    if test
        @test vexity(p) == ConvexVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 0 atol = atol rtol = rtol
        @test evaluate(norm_inf(x)) ≈ 0 atol = atol rtol = rtol
        @test norm(p.constraints[1].dual) ≈ 0 atol = atol rtol = rtol
        @test norm(p.constraints[2].dual) ≈ 0 atol = atol rtol = rtol
    end
end

@add_problem lp function lp_dual_norm_1_atom(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    x = Variable(3)
    p = minimize(norm_1(x), [-2 <= x, x <= 1]; numeric_type = T)

    if test
        @test vexity(p) == ConvexVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 0 atol = atol rtol = rtol
        @test evaluate(norm_1(x)) ≈ 0 atol = atol rtol = rtol
        @test norm(p.constraints[1].dual) ≈ 0 atol = atol rtol = rtol
        @test norm(p.constraints[2].dual) ≈ 0 atol = atol rtol = rtol
    end
end


text:
# tests for the usual functions: optBL, extractQuartet, moves
# for network with n=6 (1_astral.out, add hybridization; true network)
# messy code intended to find specific bugs, does not follow an order
# Claudia April 2015

include("../src/types.jl")
include("../src/functions.jl")

quartets = readListQuartets("allQuartets.txt");

truenet = readTopologyUpdate("truenetwork.txt");
printEdges(truenet)
printNodes(truenet)

tree = readTopologyUpdate("1_astral.out");
printEdges(tree)
writeTopologyLevel1(tree)
net = deepcopy(tree);
srand(1234)
addHybridizationUpdate!(net);
printEdges(net)
writeTopologyLevel1(net)

extractQuartet!(truenet, quartets) # no errors
extractQuartet!(net,quartets) # errors: after delete leaf!: fixed!

printNodes(net)
nettest = deepcopy(net);
extractQuartet!(nettest,quartets[1]) # error here

qnet = QuartetNetwork(net);
deleteLeaf!(qnet,qnet.node[3])
printEdges(qnet)

## qnet2 = deepcopy(qnet); # for testing
## printNodes(qnet2)

## deleteLeaf!(qnet2,qnet2.node[4]) # no error!!!

deleteLeaf!(qnet,qnet.node[4])
updateHasEdge!(qnet,net)
parameters!(qnet,net) # error here: fixed

# now, error in minor edge 12 set as identifiable AND node.k=-1 even after updateInCycle
tree = readTopologyUpdate("1_astral.out");
printEdges(tree)
writeTopologyLevel1(tree)
net = deepcopy(tree);
srand(1234)
hybrid = addHybridization!(net);
printNodes(net)
printEdges(net)
flag, nocycle, edgesInCycle, nodesInCycle = updateInCycle!(net,hybrid);
net.node[11].k #4, so correct!
updateMajorHybrid!(net,hybrid);
flag2, edgesGammaz = updateGammaz!(net,hybrid,allow);

success,hybrid,flag,nocycle,flag2,flag3 = addHybridizationUpdate!(net);

# ----------------
include("../src/types.jl")
include("../src/functions.jl")

quartets = readListQuartets("allQuartets.txt");
df = readtable("HGT_truenet_expCF.csv")
d = readTableCF(df); #expCF

tree = readTopologyUpdate("1_astral.out");
printEdges(tree)
net = deepcopy(tree);
srand(1234)
addHybridizationUpdate!(net);
printEdges(net)

currT = deepcopy(net);
optBL!(currT,d)
newT = deepcopy(currT);
count = 0
move = whichMove(currT)
move = :CHdir
move = :MVorigin
move = :MVtarget
move = :nni

movescount = zeros(Int,18)
movesfail = zeros(Int,6)
flag = proposedTop!(move,newT,true,count,1, movescount, movesfail)
printEdges(newT)
printNodes(newT)
optBL!(newT,d)

newT0 = deepcopy(newT);

qnet = QuartetNetwork(newT);
q=quartets[3];
q.taxon
newT.names
extractQuartet!(newT,quartets) # error

qnet = QuartetNetwork(newT);
printNodes(qnet)
deleteLeaf!(qnet,qnet.node[4])
printNodes(qnet)
printEdges(qnet)
qnet0 = deepcopy(qnet);
qnet = deepcopy(qnet0);

identifyQuartet!(qnet)
qnet.which != 1 ? error("qnet which not correctly assigned") : nothing
qnet.hybrid[1].k != 2 ? error("qnet.hybrid[1].k not correctly assigned") : nothing
qnet.hybrid[1].typeHyb != 1 ? error("qnet.hybrid[1].typeHyb not correctly assigned") : nothing
qnet.hybrid[1].prev.number != -5 ? error("qnet.hybrid[1].prev not correctly assigned") : nothing

eliminateHybridization!(qnet) #error here
size(qnet.hybrid,1) != 0 || qnet.numHybrids != 0 ? error("qnet should not have hybrid nodes anymore") : nothing
qnet.t1 != 0.2-log(1-0.1*(1-exp(-1.1))) ? error("internal edge length not correctly updated") : nothing

#---
node = qnet.hybrid[1];

updateSplit!(qnet)
qnet.split != [1,1,2,2] ? error("qnet.split not correctly assigned") : nothing

updateFormula!(qnet)
qnet.formula != [2,1,2] ? error("qnet.formula not correctly assigned") :nothing

calculateExpCF!(qnet)
qnet.expCF != [1/3*exp(-qnet.t1),1-2/3*exp(-qnet.t1),1/3*exp(-qnet.t1)] ? error("qnet.expCF wrongly calculated") : nothing


text:
#=
Multi-file related type definitions
=#

mutable struct MFAttributes{T} <: BaseAttributes where T <: BaseAttributes
    as::Vector{T}
end

function Base.getindex(a::MFAttributes,name::AbstractString)
    return a.as[1][name]
end

function Base.setindex!(a::MFAttributes,data,name::AbstractString)
    for a in a.as
        a[name] = data
    end
    return data
end

Base.keys(a::MFAttributes) = keys(a.as[1])

mutable struct MFDimensions{T} <: AbstractDimensions where T <: AbstractDimensions
    as::Vector{T}
    aggdim::String
end

mutable struct MFGroups{T} <: AbstractGroups where T <: AbstractGroups
    as::Vector{T}
    aggdim::String
end

function Base.getindex(a::MFDimensions,name::AbstractString)
    if name == a.aggdim
        return sum(d[name] for d in a.as)
    else
        return a.as[1][name]
    end
end

function Base.setindex!(a::MFDimensions,data,name::AbstractString)
    for a in a.as
        a[name] = data
    end
    return data
end

Base.keys(a::Union{MFDimensions,MFGroups}) = keys(a.as[1])

unlimited(a::MFDimensions) = unique(reduce(hcat,unlimited.(a.as)))

function Base.getindex(a::MFGroups,name::AbstractString)
    ds = getindex.(a.as,name)
    attrib = MFAttributes([d.attrib for d in ds])
    dim = MFDimensions([d.dim for d in ds],a.aggdim)
    group = MFGroups([d.group for d in ds],a.aggdim)

    return MFDataset(ds,a.aggdim,attrib,dim,group)
end

#---
mutable struct MFDataset{T,N,S<:AbstractString,TA,TD,TG} <: AbstractDataset where T <: AbstractDataset
    ds::Array{T,N}
    aggdim::S
    attrib::MFAttributes{TA}
    dim::MFDimensions{TD}
    group::MFGroups{TG}
    _boundsmap::Dict{String,String}
end

mutable struct MFVariable{T,N,M,TA,A} <: AbstractVariable{T,N}
    var::CatArrays.CatArray{T,N,M,TA}
    attrib::MFAttributes{A}
    dimnames::NTuple{N,String}
    varname::String
end

Base.Array(v::MFVariable) = Array(v.var)

iswritable(mfds::MFDataset) = iswritable(mfds.ds[1])

function MFDataset(ds,aggdim,attrib,dim,group)
    _boundsmap = Dict{String,String}()
    mfds = MFDataset(ds,aggdim,attrib,dim,group,_boundsmap)
    if !iswritable(mfds)
        initboundsmap!(mfds)
    end
    return mfds
end

"""
    mfds = NCDataset(fnames, mode = "r"; aggdim = nothing, deferopen = true)

Opens a multi-file dataset in read-only "r" or append mode "a". `fnames` is a
vector of file names. You can use [Glob.jl](https://github.com/vtjnash/Glob.jl)
to make `fnames`, e.g.
```julia
using NCDatasets, Glob
ds = NCDataset(glob("ERA5_monthly3D_reanalysis_*.nc"))
```

Variables are aggregated over the first unlimited dimension or over
the dimension `aggdim` if specified. The append mode is only implemented when
`deferopen` is `false`.

All variables containing the dimension `aggdim` are aggregated. The variable who
do not contain the dimension `aggdim` are assumed constant.

If deferopen is `false`, all files are opened at the same time.
However the operating system might limit the number of open files. In Linux,
the limit can be controled with the [command `ulimit`](https://stackoverflow.com/questions/34588/how-do-i-change-the-number-of-open-files-limit-in-linux).

All metadata (attributes and dimension length are assumed to be the same for all
NetCDF files. Otherwise reading the attribute of a multi-file dataset would be
ambiguous. An exception to this rule is the length of the dimension over which
the data is aggregated. This aggregation dimension can varify from file to file.

Setting the experimental flag `_aggdimconstant` to `true` means that the
length of the aggregation dimension is constant. This speeds up the creating of
a multi-file dataset as only the metadata of the first file has to be loaded.
"""
function NCDataset(fnames::AbstractArray{TS,N},mode = "r"; aggdim = nothing, deferopen = true,
                   _aggdimconstant = false,
                   ) where N where TS <: AbstractString
    if !(mode == "r" || mode == "a")
        throw(NetCDFError(-1,"""Unsupported mode for multi-file dataset (mode = $(mode)). Mode must be "r" or "a". """))
    end

    if deferopen
        @assert mode == "r"

        if _aggdimconstant
            # load only metadata from master
            master_index = 1
            ds_master = NCDataset(fnames[master_index],mode);
            data_master = metadata(ds_master)
            ds = Vector{Union{NCDataset,DeferDataset}}(undef,length(fnames))
            #ds[master_index] = ds_master
            for i = 1:length(fnames)
                #if i !== master_index
                ds[i] = DeferDataset(fnames[i],mode,data_master)
                #end
            end
        else
            ds = DeferDataset.(fnames,mode)
        end
    else
        ds = NCDataset.(fnames,mode);
    end

    if aggdim == nothing
        # first unlimited dimensions
        aggdim = NCDatasets.unlimited(ds[1].dim)[1]
    end

     attrib = MFAttributes([d.attrib for d in ds])
    dim = MFDimensions([d.dim for d in ds],aggdim)
    group = MFGroups([d.group for d in ds],aggdim)

    mfds = MFDataset(ds,aggdim,attrib,dim,group)
    return mfds
end

function close(mfds::MFDataset)
    close.(mfds.ds)
    return nothing
end

function sync(mfds::MFDataset)
    sync.(mfds.ds)
    return nothing
end

function path(mfds::MFDataset)
    path(mfds.ds[1]) * "…" * path(mfds.ds[end])
end
groupname(mfds::MFDataset) = groupname(mfds.ds[1])

function Base.keys(mfds::MFDataset)
    if mfds.aggdim == ""
        return unique(Iterators.flatten(keys.(mfds.ds)))
    else
        keys(mfds.ds[1])
    end
end

Base.getindex(v::MFVariable,indexes::Union{Int,Colon,UnitRange{Int},StepRange{Int,Int}}...) = getindex(v.var,indexes...)
Base.setindex!(v::MFVariable,data,indexes::Union{Int,Colon,UnitRange{Int},StepRange{Int,Int}}...) = setindex!(v.var,data,indexes...)
Base.size(v::MFVariable) = size(v.var)
dimnames(v::MFVariable) = v.dimnames
name(v::MFVariable) = v.varname

function variable(mfds::MFDataset,varname::SymbolOrString)
    if mfds.aggdim == ""
        # merge all variables

        # the latest dataset should be used if a variable name is present multiple times
        for ds in reverse(mfds.ds)
            if haskey(ds,varname)
                return variable(ds,varname)
            end
        end
    else
        # aggregated along a given dimension
        vars = variable.(mfds.ds,varname)

        dim = findfirst(dimnames(vars[1]) .== mfds.aggdim)
        @debug "dim $dim"

        if (dim != nothing)
            v = CatArrays.CatArray(dim,vars...)
            return MFVariable(v,MFAttributes([var.attrib for var in vars]),
                          dimnames(vars[1]),varname)
        else
            return vars[1]
        end
    end
end


text:
include("derivative.jl")

function quadratic(f)
    # Compute the first derivative of f
	f1 = derivative(f)

    # Get the y intercept (explictly passing a floating point number)
	c = f(0.0)

	# Get the y intercept of the first derivative
	b = f1(0.0)

	a = f(1.0) - b - c

    # Our old friend the quadratic formula
    # Notice how Julia lets you return multiple values
    # Multiple values should be separated by a comma and
    # are returned as a tuple
    # return (-b + sqrt(b^2 -4a*c)) / 2a, (-b - sqrt(b^2 -4a*c)) / 2a
	return (-b + sqrt(b^2 -4a*c + 0im)) / 2a, (-b - sqrt(b^2 -4a*c + 0im)) / 2a
end


text:
data_union = Union{Nothing, AbstractString, Dict, Array{T, 1} where T}
init_union = Union{Nothing, StanBase.Init, AbstractString, Dict, Array{T, 1} where T}

"""

Default `output_base` data files, in tmpdir.

$(SIGNATURES)

# Extended help

Internal, not exported.
"""
data_file_path(output_base::AbstractString, id::Int) = output_base * "_data_$(id).R"

"""

Default `output_base` init files, in tmpdir.

$(SIGNATURES)

# Extended help

Internal, not exported.
"""
init_file_path(output_base::AbstractString, id::Int) =
  output_base * "_init_$(id).R"

"""

Default `output_base` chain files, in tmpdir.

$(SIGNATURES)

# Extended help

Internal, not exported.
"""
sample_file_path(output_base::AbstractString, id::Int) =
  output_base * "_chain_$(id).csv"

"""

Default `output_base` for generated_quatities files, in tmpdir.

$(SIGNATURES)

# Extended help

Internal, not exported.
"""
generated_quantities_file_path(output_base::AbstractString, id::Int) =
  output_base * "_generated_quantities_$(id).csv"

"""

Default `output_base` log files, in tmpdir.

$(SIGNATURES)

# Extended help

Internal, not exported.
"""
log_file_path(output_base::AbstractString, id::Int) =
  output_base * "_log_$(id).log"

"""

Default `output_base` diagnostic files, in tmpdir.

$(SIGNATURES)

# Extended help

Internal, not exported.
"""
diagnostic_file_path(output_base::AbstractString, id::Int) =
  output_base * "_diagnostic_$(id).csv"

"""

Execute the method contained in a CmdStanModel.

$(SIGNATURES)

# Extended help

### Required arguments
```julia
* `model <: CmdStanModels`             : CmdStanModel subtype, e.g. SampleModel
```

### Keyword arguments
```julia
* `init`                               : Init dict
* `data`                               : Data dict
* `n_chains=4`                         : Update number of chains.
* `seed=-1`                            : Set seed value (::Int)
```

### Returns
```julia
* `rc`                                 : Return code, 0 is success
```
"""
function stan_sample(model::T; kwargs...) where {T <: CmdStanModels}

    # How many chains?
    if :n_chains in keys(kwargs)
        model.n_chains = [kwargs[:n_chains]]
    end
    n_chains = get_n_chains(model)

    # Diagnostics files requested?
    diagnostics = false
    if :diagnostics in keys(kwargs)
        diagnostics = kwargs[:diagnostics]
        setup_diagnostics(model, n_chains)
    end

    if :seed in keys(kwargs)
        model.seed = RandomSeed(kwargs[:seed])
    end

    # Remove existing sample files
    for id in 1:n_chains
        sfile = sample_file_path(model.output_base, id)
        isfile(sfile) && rm(sfile)
    end

    :init in keys(kwargs) && update_R_files(model, kwargs[:init], n_chains, "init")
    :data in keys(kwargs) && update_R_files(model, kwargs[:data], n_chains, "data")

    model.cmds = [stan_cmds(model, id; kwargs...) for id in 1:n_chains]

    #println(typeof(model.cmds))
    #println()
    #println(model.cmds)

    run(pipeline(par(model.cmds), stdout=model.log_file[1]))
end

"""

Generate a cmdstan command line (a run `cmd`).

$(SIGNATURES)

Internal, not exported.
"""
function stan_cmds(model::T, id::Integer; kwargs...) where {T <: CmdStanModels}
    append!(model.sample_file, [sample_file_path(model.output_base, id)])
    append!(model.log_file, [log_file_path(model.output_base, id)])
    if length(model.diagnostic_file) > 0
      append!(model.diagnostic_file, [diagnostic_file_path(model.output_base, id)])
    end
    cmdline(model, id)
end

"""

Update data or init R files.

$(SIGNATURES)

# Extended help

### Required arguments
```julia
* `model`                         : CmdStanModels object
* `input`                         : Input data or init values
* `n_chains`                      : Number of chains in model
* `fname_part="data"`             : Data or init R files to be created
```

Not exported.
"""
function update_R_files(model, input, n_chains, fname_part="data")

  model_field = fname_part == "data" ? model.data_file : model.init_file
  if typeof(input) <: NamedTuple || typeof(input) <: Dict
    for i in 1:n_chains
      stan_dump(model.output_base*"_$(fname_part)_$i.R", input, force=true)
      append!(model_field, [model.output_base*"_$(fname_part)_$i.R"])
    end
  elseif  typeof(input) <: Array
    if length(input) == n_chains
      for (i, d) in enumerate(input)
        stan_dump(model.output_base*"_$(fname_part)_$i.R", d, force=true)
        append!(model_field, [model.output_base*"_$(fname_part)_$i.R"])
      end
    else
      @info "Data vector length does not match number of chains,"
      @info "only first element in data vector will be used,"
      for i in 1:n_chains
        stan_dump(model.output_base*"_$(fname_part)_$i.R", input[1], force=true)
        append!(model_field, [model.output_base*"_$(fname_part)_$i.R"])
      end
    end
  elseif typeof(input) <: AbstractString && length(input) > 0
    for i in 1:n_chains
      cp(input, "$(model.output_base)_$(fname_part)_$i.R", force=true)
      append!(model_field, [model.output_base*"_$(fname_part)_$i.R"])
    end
  else
    error("\nUnrecognized input argument: $(typeof(input))\n")
  end

end

"""
Helper function for the (deprecated) diagnostics file generation.

$(SIGNATURES)

# Extended help

I am not aware the diagnostic files contain other info then the regular .csv files.
Currently I have not disabled this functionality. Please let me know if this
feature should be included/enabled.
"""
function setup_diagnostics(model, n_chains)
  for i in 1:n_chains
    append!(model.diagnostic_file, [model.output_base*"_diagnostic_$i.log"])
  end
end


text:
module MeshSampling

using StaticArrays
using GeometryTypes
using FileIO

export samplePointsFromMesh

function getFaceArea(points, face::Face)
    a = points[face[1]]
    b = points[face[2]]
    c = points[face[3]]
    return norm(cross(b - a, c - a)) / 2.
end

function samplePointsFromMesh(N_points, points, faces)
    # Get normalized cumulative area of faces
    cum_area = Array{Float64}(size(faces, 1))
    for i = 1:size(faces, 1)
        cum_area[i] = getFaceArea(points, faces[i])
    end
    cum_area = cumsum(cum_area)
    cum_area /= cum_area[end]

    sampled_pts = Array{Float64}(3, N_points)
    sampled_faces = zeros(N_points, size(faces, 1))
    for i=1:N_points
        # Pick the face we'll sample from
        sample = rand()
        k = findfirst(cum_area .>= sample)
        a = points[faces[k][1]]
        b = points[faces[k][2]]
        c = points[faces[k][3]]

        sampled_faces[i, k] = 1

        s1 = rand()
        s2 = rand()
        while (s1 + s2 > 1.0)
            s1 = rand()
            s2 = rand()
        end
        pt = a + s1*(b - a) + s2*(c - a)
        for k=1:3
            sampled_pts[k, i] = pt[k]
        end
    end
    return (sampled_pts, sampled_faces)
end

end # module

max_stars_repo_name: gizatt/PoseEstimationSandbox

max_issues_repo_name: gizatt/PoseEstimationSandbox

max_forks_repo_name: gizatt/PoseEstimationSandbox


text:
splat(f) = args -> f(args...)
push1(x) = push(x, 1.0)

function calibrationmatrices(check, extrinsic)
    path = @__DIR__()
    mat"""addpath(path,genpath($path));
    warning('off','all')"""
    mat"[$ttform, $errors] = justextrinsic($check, $extrinsic)"
    tform = ttform'
    itform = inv(tform)
    return (tform, itform, errors)
end

"""
    buildcalibration(check, extrinsic)

Build a calibration object from just a single checkerboard image (i.e. extrinsic). `check`
is the width of the checker, and `extrinsic` is the file name of the image.
"""
function buildcalibration(check, extrinsic::AbstractString)
    _tform, _itform, ϵ = calibrationmatrices(check, extrinsic)
    M = LinearMap(SMatrix{3,3, Float64}(_tform))
    tform = PerspectiveMap() ∘ M ∘ push1
    iM = LinearMap(SMatrix{3,3, Float64}(_itform))
    itform = PerspectiveMap() ∘ iM ∘ push1
    (; tform, itform, ϵ)
end

"""
    buildcalibration(check, extrinsic, ::Missing)

A convinience function for building a calibration object from just a single checkerboard
image (i.e. extrinsic).
"""
buildcalibration(check, extrinsic, ::Missing) = buildcalibration(check, extrinsic)

function calibrationmatrices(check, extrinsic, intrinsic)
    path = @__DIR__()
    mat"""addpath(path,genpath($path));
    warning('off','all')"""

    mat"[$sz, cameraParams, R, t, $errors] = bothcalibrations($check, $extrinsic, $intrinsic)"

    h, w = Int.(sz)
    rows = 1:h
    cols = 1:w
    points = vcat(map.(hcat, rows, cols')...)

    mat"$world = pointsToWorld(cameraParams, R, t, $points);"
    tform = reshape(V2.(eachrow(world)), h, w)'
    _tform = (; tform, cols, rows)

    mx, Mx = extrema(world[:,1])
    my, My = extrema(world[:,2])
    xs = UnitRange(mx, Mx)
    ys = UnitRange(my, My)
    worldPoints = vcat(map.(hcat, xs, ys', 0.0)...)
    mat"$projectedPoints = worldToImage(cameraParams,R,t,$worldPoints);"
    itform = reshape(V2.(eachrow(projectedPoints)), length(xs), length(ys))
    _itform = (; itform, xs, ys)

    return (_tform, _itform, errors)
end

"""
    buildcalibration(check, extrinsic, intrinsic)

Build a calibration object from extrinsic and intrinsic images. `check`
is the width of the checker, `extrinsic` is the file name of the image, `intrinsic` is a
collection of the intrinsic image files.
"""
function buildcalibration(check, extrinsic::AbstractString, intrinsic::AbstractVector)
    tf, itf, ϵ = calibrationmatrices(check, extrinsic, intrinsic)
    _tform = interpolate(tf.tform, BSpline(Linear()))
    tform = splat(extrapolate(scale(_tform, tf.cols, tf.rows), Flat()))
    _itform = interpolate(itf.itform, BSpline(Linear()))
    itform = splat(extrapolate(scale(_itform, itf.xs, itf.ys), Flat()))
    (; tform, itform, ϵ)
end


text:
@defcomp nice_recycle begin

    # --------------------
    # Model Indices
    # --------------------

    regions                  = Index()                                     # Index for RICE regions.
    quintiles                = Index()                                     # Index for regional income quintiles.

    # --------------------
    # Model Parameters
    # --------------------

    min_study_gdp            = Parameter()                                 # Minimum observed per capita GDP value found in elasticity studies ($/person).
    max_study_gdp            = Parameter()                                 # Maximum observed per capita GDP value found in elasticity studies ($/person).
    elasticity_intercept     = Parameter()                                 # Intercept term to estimate time-varying income elasticity.
    elasticity_slope         = Parameter()                                 # Slope term to estimate time-varying income elasticity.
    damage_elasticity        = Parameter()                                 # Income elasticity of climate damages (1 = proportional to income)
    lost_revenue_share       = Parameter()                                 # Share of carbon tax revenue that is lost and cannot be recycled (1 = 100% of revenue lost, 0 = nothing lost)
    global_carbon_tax        = Parameter(index=[time])                     # Carbon tax ($/ton CO₂).
    global_recycle_share     = Parameter(index=[regions])                  # Shares of regional revenues that are recycled globally as international transfers (1 = 100% of revenue recycled globally).
    regional_population      = Parameter(index=[time, regions])            # Regional population levels (millions of people).
    ABATEFRAC                = Parameter(index=[time, regions])            # Cost of CO₂ emission reductions as share of gross economic output.
    DAMFRAC                  = Parameter(index=[time, regions])            # Climate damages as share of gross output.
    Y                        = Parameter(index=[time, regions])            # Gross world product net of abatement and damages (trillions 2005 USD yr⁻¹).
    CPC                      = Parameter(index=[time, regions])            # Regional per capita consumption (thousands 2005 USD yr⁻¹)
    industrial_emissions     = Parameter(index=[time, regions])            # Industrial carbon emissions (GtC yr⁻¹).
    recycle_share            = Parameter(index=[regions, quintiles])       # Share of carbon tax revenue recycled back to each quintile.
    quintile_income_shares   = Parameter(index=[time, regions, quintiles]) # Quintile share of regional income (can vary over time).

    # --------------------
    # Model Variables
    # --------------------

    global_pc_revenue        = Variable(index=[time])                      # Per capita carbon tax revenue from globally recycled regional revenues ($1000s/person).
    pc_gdp                   = Variable(index=[time, regions])             # Per capita output net of abatement and damages (2005 USD / person).
    CO₂_income_elasticity    = Variable(index=[time, regions])             # Elasticity of CO₂ price exposure with respect to income.
    tax_revenue              = Variable(index=[time, regions])             # Total carbon tax revenue (2005 USD).
    regional_pc_revenue      = Variable(index=[time, regions])             # Total per capita carbon tax revenue, including any international transfers ($1000s/person).
    damage_dist              = Variable(index=[time, regions, quintiles])  # Quintile share of regional climate damages (can vary over time).
    abatement_cost_dist      = Variable(index=[time, regions, quintiles])  # Time-varying regional CO₂ tax distribution share by quintile.
    carbon_tax_dist          = Variable(index=[time, regions, quintiles])  # Time-varying regional CO₂ abatement cost distribution share by quintile.
    qc_base                  = Variable(index=[time, regions, quintiles])  # Pre-damage, pre-abatement cost, pre-tax quintile consumption (thousands 2005 USD yr⁻¹).
    qc_post_damage_abatement = Variable(index=[time, regions, quintiles])  # Post-damage, post-abatement cost per capita quintile consumption (thousands 2005 USD/person yr⁻¹).
    qc_post_tax              = Variable(index=[time, regions, quintiles])  # Quintile per capita consumption after subtracting out carbon tax (thousands 2005 USD/person yr⁻¹).
    qc_post_recycle          = Variable(index=[time, regions, quintiles])  # Quintile per capita consumption after recycling tax back to quintiles (thousands 2005 USD/person yr⁻¹).

    function run_timestep(p, v, d, t)

        for r in d.regions

            # Calculate net per capita income ($/person).
            # Note, Y in $trillions and population in millions, so scale by 1e6.
            v.pc_gdp[t,r] = p.Y[t,r] / p.regional_population[t,r] * 1e6

            # Calculate time-varying income elasticity of CO₂ price exposure (requires pc_gdp units in $/person).
            # Note, hold elasticity constant at boundary value if GDP falls outside the study support range.

            if v.pc_gdp[t,r] < p.min_study_gdp
                # GDP below observed study values.
                v.CO₂_income_elasticity[t,r] = p.elasticity_intercept + p.elasticity_slope * log(p.min_study_gdp)
            elseif v.pc_gdp[t,r] > p.max_study_gdp
                # GDP above observed study values.
                v.CO₂_income_elasticity[t,r] = p.elasticity_intercept + p.elasticity_slope * log(p.max_study_gdp)
            else
                # GDP within observed study values.
                v.CO₂_income_elasticity[t,r] = p.elasticity_intercept + p.elasticity_slope * log(v.pc_gdp[t,r])
            end

            # Calculate total carbon tax revenue for each region (dollars).
            # Note, emissions in GtC and tax in dollars, so scale by 1e9.
            v.tax_revenue[t,r] = (p.industrial_emissions[t,r] * p.global_carbon_tax[t] * 1e9) * (1.0 - p.lost_revenue_share)
        end

        # Calculate per capita tax revenue from globally recycled revenue (convert to $1000s/person to match pc consumption units).
        # Note, tax in dollars and population in millions, so scale by 1e9.
        v.global_pc_revenue[t] = sum(v.tax_revenue[t,:] .* p.global_recycle_share[:]) / sum(p.regional_population[t,:]) / 1e9

        for r in d.regions

            # Calculate total recycled per capita tax revenue for each region (this also includes the globally recycled revenue).
            v.regional_pc_revenue[t,r] = (v.tax_revenue[t,r] * (1-p.global_recycle_share[r])) / p.regional_population[t,r] / 1e9 + v.global_pc_revenue[t]

            # Calculate quintile distribution shares of CO₂ tax burden and mitigation costs (assume both distributions are equal) and cliamte damages.
            v.abatement_cost_dist[t,r,:] = regional_quintile_distribution(v.CO₂_income_elasticity[t,r], p.quintile_income_shares[t,r,:])
            v.carbon_tax_dist[t,r,:]     = regional_quintile_distribution(v.CO₂_income_elasticity[t,r], p.quintile_income_shares[t,r,:])
            v.damage_dist[t,r,:]         = regional_quintile_distribution(p.damage_elasticity, p.quintile_income_shares[t,r,:])

            # Create a temporary variable used to calculate NICE baseline quintile consumption (just for convenience).
            temp_C = 5.0 * p.CPC[t,r] * (1.0 + p.DAMFRAC[t,r]) / (1.0 - p.ABATEFRAC[t,r])

            for q in d.quintiles

                # Calculate pre-damage, pre-abatement cost quintile consumption.
                v.qc_base[t,r,q] = temp_C * p.quintile_income_shares[t,r,q]

                # Calculate post-damage, post-abatement cost per capita quintile consumption (bounded below to ensure consumptions don't collapse to zero or go negative).
                # Note, this differs from standard NICE equation because quintile CO₂ abatement cost and climate damage shares can now vary over time.
                v.qc_post_damage_abatement[t,r,q] = max(v.qc_base[t,r,q] - (5.0 * p.CPC[t,r] * p.DAMFRAC[t,r] * v.damage_dist[t,r,q]) - (temp_C * p.ABATEFRAC[t,r] * v.abatement_cost_dist[t,r,q]), 1e-8)

                # Subtract tax revenue from each quintile based on quintile CO₂ tax burden distributions.
                # Note, per capita tax revenue and consumption should both be in 1000s dollars/person.
                v.qc_post_tax[t,r,q] = v.qc_post_damage_abatement[t,r,q] - (5 * v.regional_pc_revenue[t,r] * v.carbon_tax_dist[t,r,q])

                # Recycle tax revenue by adding shares back to all quintiles (assume recycling shares constant over time).
                v.qc_post_recycle[t,r,q] = v.qc_post_tax[t,r,q] + (5 * v.regional_pc_revenue[t,r] * p.recycle_share[r,q])
            end
        end
    end
end


text:
type BlackConstantVol{BC <: BusinessCalendar, DC <: DayCount} <: BlackVolTermStructure
  referenceDate::Date
  settlementDays::Int
  calendar::BC
  volatility::Quote
  dc::DC
end

BlackConstantVol(refDate::Date, cal::BusinessCalendar, volatility::Float64, dc::DayCount) = BlackConstantVol(refDate, 0, cal, Quote(volatility), dc)

black_vol_impl(bts::BlackConstantVol, ::Float64, ::Float64) = bts.volatility.value

function black_variance_impl(bts::BlackConstantVol, t::Float64, strike::Float64)
  vol = black_vol_impl(bts, t, strike)
  return vol * vol * t
end

function black_vol(bts::BlackVolTermStructure, t::Float64, strike::Float64)
  #TODO check stuff
  return black_vol_impl(bts, t, strike)
end

black_vol(bts::BlackVolTermStructure, d::Date, strike::Float64) = black_vol(bts, time_from_reference(bts, d), strike)

function black_variance(bts::BlackVolTermStructure, d::Date, strike::Float64)
  #TODO check stuff
  t = time_from_reference(bts, d)
  return black_variance_impl(bts, t, strike)
end

function black_variance(bts::BlackVolTermStructure, t::Float64, strike::Float64)
  #TODO check stuff
  return black_variance_impl(bts, t, strike)
end


text:
#!/usr/bin/env julia
using CSV, DataFrames, PGFPlotsX
push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\usepackage{fontspec}")
push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\setmainfont{Asana-Math}")
seflu1617 = CSV.File("data/Rpt_LabSurveillanceDataAnyWeekByCtry_Chart_se_201640-201720.tsv") |> DataFrame
seflu21 = CSV.File("data/Rpt_LabSurveillanceDataAnyWeekByCtry_Chart_se_202120-202149.tsv") |> DataFrame
seflu1617[!,:AH3_NS] = seflu1617[!,:AH3] + seflu1617[!,:ANOTSUBTYPED]
seflu21[!,:AH3_NS] = seflu21[!,:AH3] + seflu21[!,:ANOTSUBTYPED]
p = @pgf Axis({ymode="log",xlabel="vecka", ylabel="antal labbfall", xmajorgrids,
    title="Fall influensa A(H3N2) och A ej subtypad Sverige",
    xtick="{1,10,13,33}",xticklabels="{40,49,52,20}"})
@pgf push!(p, PlotInc({"mark=x", "color=blue"}, Table([1:33,seflu1617[!,:AH3_NS]])), LegendEntry("2016/17"))
@pgf push!(p, PlotInc({"mark=*", "color=red"}, Table([1:10,seflu21[21:30,:AH3_NS]])), LegendEntry("2021/22"))
inflpl = "se_flu_1617_4020_21_4049"
inflplpath = "../../images/$(inflpl).svg"
print(*("![Labbfall influensa A(H3N2) och A ej subtypad Sverige 2016/17 och 2021/22.]",
    "($(inflplpath)){#fig:$(inflpl) width=100%}\n\n"))
pgfsave(inflplpath, p)

max_stars_repo_name: klpn/static-dust

max_issues_repo_name: klpn/static-dust

max_forks_repo_name: klpn/static-dust


text:
module DemoProject

export calcHypotenuse
include("pythagoras.jl")
end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export 'feffery'_antdimage

"""
    'feffery'_antdimage(;kwargs...)

An AntdImage component.

Keyword arguments:
- `id` (String; optional)
- `alt` (String; optional)
- `className` (String; optional)
- `fallback` (String; optional)
- `height` (Real; optional)
- `loading_state` (optional): . loading_state has the following type: lists containing elements 'is_loading', 'prop_name', 'component_name'.
Those elements have the following types:
  - `is_loading` (Bool; optional): Determines if the component is loading or not
  - `prop_name` (String; optional): Holds which property is loading
  - `component_name` (String; optional): Holds the name of the component that is loading
- `locale` (a value equal to: 'zh-cn', 'en-us'; optional)
- `multiImageMode` (a value equal to: 'fold', 'unfold'; optional)
- `preview` (Bool; optional)
- `src` (String | Array of Strings; optional)
- `style` (Dict; optional)
- `width` (Real; optional)
"""
function 'feffery'_antdimage(; kwargs...)
        available_props = Symbol[:id, :alt, :className, :fallback, :height, :loading_state, :locale, :multiImageMode, :preview, :src, :style, :width]
        wild_props = Symbol[]
        return Component("'feffery'_antdimage", "AntdImage", "feffery_antd_components", available_props, wild_props; kwargs...)
end

max_stars_repo_name: ReilYoung/feffery-antd-components

max_issues_repo_name: ReilYoung/feffery-antd-components

max_forks_repo_name: ReilYoung/feffery-antd-components


text:
using LinearAlgebra

@doc """
Gaussian kernel function.

# Examples
```jldoctest
y = k(x)
```
where the 'x' is Vector{Float64} sample and 'y' is density in point 'x'.
""" ->
k(x::Vector{Float64}) = √(2*π)*ℯ^(-(x'*x)/2)

@doc """
Function to return Parzen window estimation.

# Examples
```jldoctest
fun = Σ(h, X, f, x)
```
where the 'fun' is parzen window estimation, 'h' is window size,
'X' is set of train samples (size(dim, N) ~ (dimension of data, data count),
'f' is kernel function and 'x' is given Vector{Float64} sample.
""" ->
create_parzen_window(h::Float64, X::Matrix{Float64},
                        f, x::Vector{Float64}) = (1/(h*size(X)[1]))*mapreduce(a->f((x-a)/h), +, eachcol(X))

max_stars_repo_name: kozvojtex/AnomalyDetection

max_issues_repo_name: kozvojtex/AnomalyDetection

max_forks_repo_name: kozvojtex/AnomalyDetection