text:
"""Read REISE input matfiles, return parsed relevant data in a Dict."""
function read_case(filepath)
    println("Reading from folder: " * filepath)

    println("...loading case.mat")
    case_mat_file = MAT.matopen(joinpath(filepath, "case.mat"))
    mpc = read(case_mat_file, "mpc")

    # New case.mat analog
    case = Dict()

    # AC branches
    # dropdims() will remove extraneous dimension
    case["branchid"] = dropdims(mpc["branchid"]; dims=2)
    # convert() will convert float array to int array
    case["branch_from"] = convert(Array{Int,1}, mpc["branch"][:, 1])
    case["branch_to"] = convert(Array{Int,1}, mpc["branch"][:, 2])
    case["branch_reactance"] = mpc["branch"][:, 4]
    case["branch_rating"] = mpc["branch"][:, 6]

    # DC branches
    if "dcline" in keys(mpc)
        if isa(mpc["dclineid"], Int)
            case["dclineid"] = Int64[mpc["dclineid"]]
        else
            case["dclineid"] = dropdims(mpc["dclineid"]; dims=2)
        end
        case["dcline_from"] = convert(Array{Int,1}, mpc["dcline"][:, 1])
        case["dcline_to"] = convert(Array{Int,1}, mpc["dcline"][:, 2])
        case["dcline_pmin"] = mpc["dcline"][:, 10]
        case["dcline_pmax"] = mpc["dcline"][:, 11]
    else
        case["dclineid"] = Int64[]
        case["dcline_from"] = Int64[]
        case["dcline_to"] = Int64[]
        case["dcline_pmin"] = Float64[]
        case["dcline_pmax"] = Float64[]
    end

    # Buses
    case["busid"] = convert(Array{Int,1}, mpc["bus"][:, 1])
    case["bus_demand"] = mpc["bus"][:, 3]
    case["bus_zone"] = convert(Array{Int,1}, mpc["bus"][:, 7])

    # Generators
    case["genid"] = dropdims(mpc["genid"]; dims=2)
    genfuel = dropdims(mpc["genfuel"]; dims=2)
    case["genfuel"] = convert(Array{String,1}, genfuel)
    case["gen_bus"] = convert(Array{Int,1}, mpc["gen"][:, 1])
    case["gen_status"] = mpc["gen"][:, 8]
    case["gen_pmax"] = mpc["gen"][:, 9]
    case["gen_pmin"] = mpc["gen"][:, 10]
    case["gen_ramp30"] = mpc["gen"][:, 19]

    # Generator costs
    case["gencost"] = mpc["gencost"]

    # Load all relevant profile data from CSV files
    println("...loading demand.csv")
    case["demand"] = DataFrames.DataFrame(CSV.File(joinpath(filepath, "demand.csv")))

    println("...loading hydro.csv")
    case["hydro"] = DataFrames.DataFrame(CSV.File(joinpath(filepath, "hydro.csv")))

    println("...loading wind.csv")
    case["wind"] = DataFrames.DataFrame(CSV.File(joinpath(filepath, "wind.csv")))

    println("...loading solar.csv")
    case["solar"] = DataFrames.DataFrame(CSV.File(joinpath(filepath, "solar.csv")))

    return case
end

"""Read input matfile (if present), return parsed data in a Storage struct."""
function read_storage(filepath)::Storage
    # Fallback dataframe, in case there's no case_storage.mat file
    storage = Dict(
        "enabled" => false, "gen" => zeros(0, 21), "sd_table" => DataFrames.DataFrame()
    )
    try
        case_storage_file = MAT.matopen(joinpath(filepath, "case_storage.mat"))
        storage_mat_data = read(case_storage_file, "storage")
        println("...loading case_storage.mat")
        # Convert N x 1 array of strings into 1D array of Symbols (length N)
        column_symbols = Symbol.(vec(storage_mat_data["sd_table"]["colnames"]))
        storage = Dict(
            "enabled" => true,
            "gen" => storage_mat_data["gen"],
            "sd_table" =>
                DataFrames.DataFrame(storage_mat_data["sd_table"]["data"], column_symbols),
        )
    catch e
        println("File case_storage.mat not found in " * filepath)
    end

    # Convert Dict to NamedTuple
    storage = (; (Symbol(k) => v for (k, v) in storage)...)
    # Convert NamedTuple to Storage
    storage = Storage(; storage...)

    return storage
end

"""
    read_demand_flexibility(filepath, interval)

Load demand flexibility profiles and parameters from .csv files and return them in a
DemandFlexibility struct.
"""
function read_demand_flexibility(filepath, interval)::DemandFlexibility
    # Initialize demand flexibility
    demand_flexibility = Dict(
        "duration" => interval,
        "enabled" => "not_specified",
        "interval_balance" => true,
        "rolling_balance" => true,
        "input_granularity" => "AREA",
    )

    # Try loading the demand flexibility parameters
    demand_flexibility_parameters = DataFrames.DataFrame()
    try
        demand_flexibility_parameters = DataFrames.DataFrame(
            CSV.File(joinpath(filepath, "demand_flexibility_parameters.csv"))
        )
        println("...loading demand flexibility parameters")

        # Create a dictionary to hold the warning messages relevant to loading the
        # demand flexibility parameters
        demand_flexibility_params_warns = Dict(
            "duration" => (
                "The demand flexibility duration parameter is not defined. Will " *
                "default to being the size of the interval."
            ),
            "enabled" => (
                "The parameter that indicates if demand flexibility is enabled is not " *
                "defined. Will default to being enabled."
            ),
            "interval_balance" => (
                "The parameter that indicates if the interval load balance constraint " *
                "is enabled is not defined. Will default to being enabled."
            ),
            "rolling_balance" => (
                "The parameter that indicates if the rolling load balance constraint " *
                "is enabled is not defined. Will default to being enabled."
            ),
            "input_granularity" => (
                "The parameter that indicates the way demand flexibility " *
                "and demand flexibility cost are provided is not defined. " *
                "Will default to being area-based."
            ),
        )

        # Try assigning the different demand flexibility parameters from the file
        for k in keys(demand_flexibility_params_warns)
            try
                demand_flexibility[k] = demand_flexibility_parameters[1, k]
            catch e
                println(demand_flexibility_params_warns[k])
            end
        end

    catch e
        println("Demand flexibility parameters not found in " * filepath)
        println(
            "Demand flexibility parameters will default to allowing demand flexibility " *
            "to occur.",
        )
    end

    # Check the feasibility of the duration parameter
    if demand_flexibility["duration"] > interval
        @warn (
            "Demand flexibility durations greater than the interval length are " *
            "set equal to the interval length."
        )
        demand_flexibility["duration"] = interval
    end

    # Prevent the rolling_balance constraint according to the duration parameter
    demand_flexibility["rolling_balance"] &= !(demand_flexibility["duration"] == interval)

    # Try loading the demand flexibility and demand flexibility cost profiles
    for s in ["up", "dn"]
        # Pre-specify the demand flexibility and demand flexibility cost profiles
        demand_flexibility["flex_amt_" * s] = nothing
        demand_flexibility["cost_" * s] = nothing

        # Only try loading the profiles if demand flexibility is enabled
        if demand_flexibility["enabled"] == "not_specified" ||
            (demand_flexibility["enabled"])
            # Try loading the demand flexibility profiles
            try
                demand_flexibility["flex_amt_" * s] = DataFrames.DataFrame(
                    CSV.File(joinpath(filepath, "demand_flexibility_" * s * ".csv"))
                )
                println("...loading demand flexibility " * s * " profiles")
            catch e
                println("Demand flexibility " * s * " profile not found in " * filepath)
            end

            # Try loading the demand flexibility cost profiles
            try
                demand_flexibility["cost_" * s] = DataFrames.DataFrame(
                    CSV.File(joinpath(filepath, "demand_flexibility_cost_" * s * ".csv"))
                )
                println("...loading demand flexibility " * s * "-shift cost profiles")
            catch e
                println(
                    "Demand flexibility " *
                    s *
                    "-shift cost profiles not found in " *
                    filepath *
                    ". Will default to no cost for " *
                    s *
                    "-shifting demand.",
                )
            end
        end
    end

    # If demand flexibility is enabled but at least one demand flexibility profile is nothing
    if demand_flexibility["enabled"] == true && (
        isnothing(demand_flexibility["flex_amt_up"]) ||
        (isnothing(demand_flexibility["flex_amt_dn"]))
    )
        @error(
            "Demand flexibility was specified to be enabled, however at least one " *
                "demand flexibility profile is missing. Please make sure both demand " *
                "flexibility profiles are included in " *
                filepath
        )
        throw(ErrorException("See above."))
    elseif demand_flexibility["enabled"] == "not_specified"
        if !isnothing(demand_flexibility["flex_amt_up"]) &&
            (!isnothing(demand_flexibility["flex_amt_dn"]))
            demand_flexibility["enabled"] = true
        else
            if !isnothing(demand_flexibility["flex_amt_up"]) ||
                (!isnothing(demand_flexibility["flex_amt_dn"]))
                @warn (
                    "The exclusion of one of the demand flexibility profiles has resulted " *
                    "in demand flexibility not being enabled."
                )
            end
            demand_flexibility["enabled"] = false
        end
    end

    # Set the demand flexibility constraints to false if enabled is false
    if !demand_flexibility["enabled"]
        demand_flexibility["interval_balance"] = false
        demand_flexibility["rolling_balance"] = false
    end

    # Convert Dict to NamedTuple
    demand_flexibility = (; (Symbol(k) => v for (k, v) in demand_flexibility)...)

    # Convert NamedTuple to DemandFlexibility object
    demand_flexibility = DemandFlexibility(; demand_flexibility...)

    return demand_flexibility
end


text:
#!/usr/bin/env julia

variants =   ARGS[1]
bed = ARGS[2]

function read(file)
    if ismatch(r".gz$", file)   # Gzipped .gz files
		f = GZip.open(file)
		lines = readlines(f)
		close(f)
    else
		f = open(file)
		lines = readlines(f)
		close(f)
	end
	return lines
end

function expand(locci)
    list = []
    pos = split(locci, ':')
        chr = pos[1]
        region = split(pos[2], '-')
        n = parse(Int64, region[1])
        while n <= parse(Int64, region[end])
            out = "$chr" * ":" * "$n"
            push!(list, out)
            n +=  1
        end
	return list
end

function getchr(variant)
    chr = split(variant, ':')[1]
	return chr
end

function readBED(bed_file, chr)
    pos = []
    bed = read(bed_file)
    for line in bed
        bed = split(line, '\t')
        chrm = replace(bed[1], "chr", "")
        if chrm == chr
			start, last = bed[2], bed[3]
            chroms = "$chrm" * ":" * "$start" * "-" * "$last"
            region = expand(chroms)
            append!(pos, region)
        end
    end
    pos = unique(pos)
    return pos
end

function checkChromatin(variant, chromatin)
	parts  = split(chomp(variant), ':')
	chromatinPos = readBED(chromatin, parts[1])

    out = 0
	for pos in chromatinPos
       #println(typeof(variant), '\t', variant, typeof(pos), '\t', pos)
       cols = split(pos, ':')
        if (parts[1] == cols[1]) && (parts[2] == cols[2])
          out += 1
        end
    end
	return out
end

function filterChromatin(variants, bed)
	lines = read(variants)
	for line in lines
	    variant = chomp(line)
		openChrommatin = checkChromatin(variant, bed)
		if openChrommatin >= 1
            print(line)
        else
			continue
		end
	end
end

filterChromatin(variants, bed)

ext: jl


text:
# Parser
function Parser(pp::Preprocessor, sema::Sema, skip_func_body::Bool=false)
    @check_ptrs pp sema
    status = Ref{CXInit_Error}(CXInit_NoError)
    p = clang_Parser_create(pp, sema, skip_func_body, status)
    @assert status[] == CXInit_NoError
    return Parser(p)
end

dispose(x::Parser) = clang_Parser_dispose(x)

function Initialize(x::Parser)
    @check_ptrs x
    return clang_Parser_Initialize(x)
end

function getLangOpts(x::Parser)
    @check_ptrs x
    return clang_Parser_getLangOpts(x)
end

function getTargetInfo(x::Parser)
    @check_ptrs x
    return TargetInfo(clang_Parser_getTargetInfo(x))
end

function getPreprocessor(x::Parser)
    @check_ptrs x
    return Preprocessor(clang_Parser_getPreprocessor(x))
end

function getActions(x::Parser)
    @check_ptrs x
    return Sema(clang_Parser_getActions(x))
end

function parseOneTopLevelDecl(x::Parser, is_first_decl::Bool=false)
    @check_ptrs x
    return DeclGroupRef(clang_Parser_parseOneTopLevelDecl(x, is_first_decl))
end

function getCurToken(x::Parser)
    @check_ptrs x
    return Token(clang_Parser_getCurToken(x))
end

function NextToken(x::Parser)
    @check_ptrs x
    return Token(clang_Parser_NextToken(x))
end

function ConsumeToken(x::Parser)
    @check_ptrs x
    return SourceLocation(clang_Parser_ConsumeToken(x))
end

function ConsumeAnyToken(x::Parser)
    @check_ptrs x
    return SourceLocation(clang_Parser_ConsumeAnyToken(x))
end

function getCurScope(x::Parser)
    @check_ptrs x
    return Scope(clang_Parser_getCurScope(x))
end

# should be sync to Clang's implementation
function getDeclSpecContextFromDeclaratorContext(ctx::CXDeclaratorContext)
    if ctx == CXDeclaratorContext_Member
        return CXDeclSpecContext_DSC_class
    elseif ctx == CXDeclaratorContext_File
        return CXDeclSpecContext_DSC_top_level
    elseif ctx == CXDeclaratorContext_TemplateParam
        return CXDeclSpecContext_DSC_template_param
    elseif ctx == CXDeclaratorContext_TemplateArg ||
           ctx == CXDeclaratorContext_TemplateTypeArg
        return CXDeclSpecContext_DSC_template_type_arg
    elseif ctx == CXDeclaratorContext_TrailingReturn ||
           ctx == CXDeclaratorContext_TrailingReturnVar
        return CXDeclSpecContext_DSC_trailing
    elseif ctx == CXDeclaratorContext_AliasDecl || ctx == CXDeclaratorContext_AliasTemplate
        return CXDeclSpecContext_DSC_alias_declaration
    else
        return CXDeclSpecContext_DSC_normal
    end
end

function shouldEnterContext(ctx::CXDeclSpecContext)
    return ctx == CXDeclSpecContext_DSC_top_level || ctx == CXDeclSpecContext_DSC_class
end

"""
    TryAnnotateCXXScopeToken(x::Parser, entering_context)
Return true if there was an error.
"""
function TryAnnotateCXXScopeToken(x::Parser, entering_context::Bool)
    @check_ptrs x
    return clang_Parser_TryAnnotateCXXScopeToken(x, entering_context)
end

function TryAnnotateCXXScopeToken(x::Parser, ctx::CXDeclSpecContext)
    return TryAnnotateCXXScopeToken(x, shouldEnterContext(ctx))
end

function TryAnnotateCXXScopeToken(x::Parser, dlctx::CXDeclaratorContext)
    ctx = getDeclSpecContextFromDeclaratorContext(dlctx)
    return TryAnnotateCXXScopeToken(x, ctx)
end

function TryAnnotateCXXScopeToken(x::Parser)
    return TryAnnotateCXXScopeToken(x, CXDeclSpecContext_DSC_top_level)
end

"""
    TryAnnotateTypeOrScopeTokenAfterScopeSpec(x::Parser, ss::CXXScopeSpec, is_new_scope::Bool=false)
Return true if there was an error.
"""
function TryAnnotateTypeOrScopeTokenAfterScopeSpec(x::Parser, ss::CXXScopeSpec,
                                                   is_new_scope::Bool=false)
    @check_ptrs x ss
    return clang_Parser_TryAnnotateTypeOrScopeTokenAfterScopeSpec(x, ss, is_new_scope)
end


text:
using Documenter, TaylorModels

makedocs(
    modules = [TaylorModels],
    format = Documenter.HTML(),
    sitename = "TaylorModels.jl",
    authors  = "Luis Benet and David P. Sanders",
    pages = [
        "Home" => "index.md",
        "API" => "api.md"
    ]
)

Documenter.deploydocs(
    repo = "github.com/JuliaIntervals/TaylorModels.jl.git",
    target = "build",
    deps = nothing,
    make = nothing
)


text:
using AdversarialDriving
using POMDPs, POMDPPolicies, POMDPSimulators
using AutomotiveVisualization

sut_agent = BlinkerVehicleAgent(get_ped_vehicle(id=1, s=5., v=15.), TIDM(ped_TIDM_template, noisy_observations = true))
adv_ped = NoisyPedestrianAgent(get_pedestrian(id=2, s=7., v=2.0), AdversarialPedestrian())
mdp = AdversarialDrivingMDP(sut_agent, [adv_ped], ped_roadway, 0.2)

# Render a single frame
s = rand(initialstate(mdp))
c = render([mdp.roadway,
                crosswalk,
                VelocityArrow(entity = s[1], color = colorant"black"),
                VelocityArrow(entity = s[2], color = colorant"black"),
                FancyCar(car = s[2], color = colorant"blue"),
                FancyPedestrian(ped = s[1])],
                surface=AutomotiveVisualization.CairoPDFSurface(IOBuffer(), DEFAULT_CANVAS_WIDTH, DEFAULT_CANVAS_HEIGHT))
write("pedestrian_crosswalk.pdf", c)

null_action = Disturbance[PedestrianControl()]
noisy_action = Disturbance[PedestrianControl(noise = Noise((-10.,0.), -2))]

# Nominal Behavior
hist = POMDPSimulators.simulate(HistoryRecorder(), mdp, FunctionPolicy((s) -> null_action))
scenes_to_gif([s.s for s in hist], mdp.roadway, "ped_crosswalk_nominal.gif", others = [crosswalk])

# Behavior with noise
hist = POMDPSimulators.simulate(HistoryRecorder(), mdp, FunctionPolicy((s) -> noisy_action))
scenes_to_gif([s.s for s in hist], mdp.roadway, "ped_crosswalk_failure.gif", others = [crosswalk])

max_stars_repo_name: ancorso/MAST

max_issues_repo_name: ancorso/MAST


text:
using AcuteML

# Type Definition
@aml mutable struct Person "person", courseCheck
    age::UInt64, "~"
    field, "study-field"
    GPA::Float64 = 4.5, "~", GPAcheck
    courses::Vector{String}, "taken-courses"
    id::Int64, a"~"
end

@aml mutable struct University "university"
    name, a"university-name"
    people::Vector{Person}, "person"
end

@aml mutable struct Doc xd""
    university::University, "~"
end

# Value Checking Functions
GPAcheck(x) = x <= 4.5 && x >= 0

function courseCheck(age, field, GPA, courses, id)

    if field == "Mechanical Engineering"
        relevant = ["Artificial Intelligence", "Robotics", "Machine Design"]
    elseif field == "Computer Engineering"
        relevant = ["Julia", "Algorithms"]
    else
        error("study field is not known")
    end

    return any(in.(courses, Ref(relevant)))
end

# Extractor
xml = parsexml("""
<?xml version="1.0" encoding="UTF-8" standalone="yes"?><!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.0 Transitional//EN" "http://www.w3.org/TR/REC-html40/loose.dtd">
<university university-name="Julia University">
  <person id="1">
    <age>24</age>
    <study-field>Mechanical Engineering</study-field>
    <GPA>4.5</GPA>
    <taken-courses>Artificial Intelligence</taken-courses>
    <taken-courses>Robotics</taken-courses>
  </person>
  <person id="2">
    <age>18</age>
    <study-field>Computer Engineering</study-field>
    <GPA>4</GPA>
    <taken-courses>Julia</taken-courses>
  </person>
</university>
""")

# extract Doc
D = Doc(xml) # StructName(xml) like Doc(xml) extracts the data and stores them in proper format

# Now you can access all of the data by calling the fieldnames

# extract University
U = D.university

U.name # Julia University

# extract Person

P1 = U.people[1]

P1.age # 24
P1.field # Mechanical Engineering
P1.GPA # 4.5
P1.courses # ["Artificial Intelligence", "Robotics"]
P1.id # 1


text:
# This example starts a coroutine for websocket handling,
# which again splits input and output on different tasks.
# Note that if
using WebSockets
inbox = Channel{String}(10)
outbox = Channel{String}(10)

ws_task = @async WebSockets.with_logger(WebSocketLogger(stderr, Base.CoreLogging.Debug)) do
    WebSockets.open("wss://echo.websocket.org") do ws
        @sync begin
            inbox_task = @async try
                while true #!eof(ws)
                    in_data, success = readguarded(ws)
                    success || break
                    in_msg = String(in_data)
                    @wslog in_msg
                    put!(inbox, in_msg)
                end
            finally
                @debug ws " Now closing outbox, to be sure"
                close(outbox)
            end
            outbox_task = @async try
                for outmsg in outbox
                    isopen(ws) && writeguarded(ws, outmsg) || break
                end
            finally
                @debug "Closing " ws
                close(ws)
            end
        end
    end
end

put!(outbox, "Hello")
put!(outbox, "World!")

@show take!(inbox)
@show take!(inbox)

close(outbox) # close(outbox) causes outbox_task to call close(ws)
wait(ws_task)

@show istaskdone(ws_task)
@show Base.istaskfailed(ws_task)
nothing


text:
push!(LOAD_PATH,"../src/") #eso es por si Julia no encuentra el ./src
using Pkg
Pkg.activate(".")

using Documenter, Covid

makedocs(
    sitename = "Covid",
    format = Documenter.HTML(),
    modules = [Covid],
    pages=[
           "Home" => "index.md"
          ])
#aqui esta la informacion general de la pagina
deploydocs(
    repo = "github.com/mucinoab/Covid.git"
)

max_stars_repo_name: MarleenMs/Covid

max_issues_repo_name: MarleenMs/Covid

max_forks_repo_name: MarleenMs/Covid


text:
function key_index(nt::NamedTuple, k::Symbol)
    @assert haskey(nt, k)
    ntkeys = keys(nt)
    for i in eachindex(ntkeys)
        if ntkeys[i] == k
            return i
        end
    end
end


text:
# Dynamically aware constraints manager

include("constraintManager.jl")
include("../objective/LQR_objective.jl")

"""
    constraintManager_Dynamics

Contains a list of constraints and the dual variables (lambda) that correspond
with each constraint.

Unlike [`constraintManager_Base`](@ref), this constraint manager also has the
dynamics constraints embedded such that the nearest dynamically feasible
trajectory can be determined at any given time.
"""
mutable struct constraintManager_Dynamics <: constraintManager
    cList::Array{constraint, 1}
    lambdaList
    dynamicsConstraint::AL_AffineEquality
    affineLambdaList
end

#############################################################
###            Evaluate the Affine Constraints            ###
#############################################################

function evalAffineEq(cM::constraintManager_Dynamics, yCurr)
    return getRaw(cM.dynamicsConstraint, yCurr)
end

function affineBlock(cM::constraintManager_Dynamics)
    return cM.dynamicsConstraint.A
end

max_stars_repo_name: danineamati/TrajOptSOCPs


text:
function dq_ri(δ::T) where {T <: Number}
    ## Uses the referenceframe of the Kundur page 852 of dq to RI
    return T[
        sin(δ) cos(δ)
        -cos(δ) sin(δ)
    ]
end

function ri_dq(δ::T) where {T <: Number}
    #Uses the reference frame of the Kundur page 852 of RI to dq
    return T[
        sin(δ) -cos(δ)
        cos(δ) sin(δ)
    ]
end

function dq_ri(δ::Float64)
    ## Uses the referenceframe of the Kundur page 852 of dq to RI
    return Float64[
        sin(δ) cos(δ)
        -cos(δ) sin(δ)
    ]
end

function ri_dq(δ::Float64)
    #Uses the reference frame of the Kundur page 852 of RI to dq
    return Float64[
        sin(δ) -cos(δ)
        cos(δ) sin(δ)
    ]
end


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

abstract type AbstractMsg end

## Wire format description
#
# Each message has three parts, which are written in order to the worker's stream.
#  1) A header of type MsgHeader is serialized to the stream (via `serialize`).
#  2) A message of type AbstractMsg is then serialized.
#  3) Finally, a fixed boundary of 10 bytes is written.

# Message header stored separately from body to be able to send back errors if
# a deserialization error occurs when reading the message body.
struct MsgHeader
    response_oid::RRID
    notify_oid::RRID
    MsgHeader(respond_oid=RRID(0,0), notify_oid=RRID(0,0)) =
        new(respond_oid, notify_oid)
end

# Special oid (0,0) uses to indicate a null ID.
# Used instead of Union{Int, Nothing} to decrease wire size of header.
null_id(id) =  id == RRID(0, 0)

struct CallMsg{Mode} <: AbstractMsg
    f::Any
    args::Tuple
    kwargs
end
struct CallWaitMsg <: AbstractMsg
    f::Any
    args::Tuple
    kwargs
end
struct RemoteDoMsg <: AbstractMsg
    f::Any
    args::Tuple
    kwargs
end
struct ResultMsg <: AbstractMsg
    value::Any
end

# Worker initialization messages
struct IdentifySocketMsg <: AbstractMsg
    from_pid::Int
end

struct IdentifySocketAckMsg <: AbstractMsg
end

struct JoinPGRPMsg <: AbstractMsg
    self_pid::Int
    other_workers::Array
    topology::Symbol
    enable_threaded_blas::Bool
    lazy::Bool
end
struct JoinCompleteMsg <: AbstractMsg
    cpu_threads::Int
    ospid::Int
end

# Avoiding serializing AbstractMsg containers results in a speedup
# of approximately 10%. Can be removed once module Serialization
# has been suitably improved.

const msgtypes = Any[CallWaitMsg, IdentifySocketAckMsg, IdentifySocketMsg,
                     JoinCompleteMsg, JoinPGRPMsg, RemoteDoMsg, ResultMsg,
                     CallMsg{:call}, CallMsg{:call_fetch}]

for (idx, tname) in enumerate(msgtypes)
    exprs = Any[ :(serialize(s, o.$fld)) for fld in fieldnames(tname) ]
    @eval function serialize_msg(s::AbstractSerializer, o::$tname)
        write(s.io, UInt8($idx))
        $(exprs...)
        return nothing
    end
end

let msg_cases = :(@assert false "Message type index ($idx) expected to be between 1:$($(length(msgtypes)))")
    for i = length(msgtypes):-1:1
        mti = msgtypes[i]
        msg_cases = :(if idx == $i
                          $(Expr(:call, QuoteNode(mti), fill(:(deserialize(s)), fieldcount(mti))...))
                      else
                          $msg_cases
                      end)
    end
    @eval function deserialize_msg(s::AbstractSerializer)
        idx = read(s.io, UInt8)
        return $msg_cases
    end
end

function send_msg_unknown(s::IO, header, msg)
    error("attempt to send to unknown socket")
end

function send_msg(s::IO, header, msg)
    id = worker_id_from_socket(s)
    if id > -1
        return send_msg(worker_from_id(id), header, msg)
    end
    send_msg_unknown(s, header, msg)
end

function send_msg_now(s::IO, header, msg::AbstractMsg)
    id = worker_id_from_socket(s)
    if id > -1
        return send_msg_now(worker_from_id(id), header, msg)
    end
    send_msg_unknown(s, header, msg)
end
function send_msg_now(w::Worker, header, msg)
    send_msg_(w, header, msg, true)
end

function send_msg(w::Worker, header, msg)
    send_msg_(w, header, msg, false)
end

function flush_gc_msgs(w::Worker)
    if !isdefined(w, :w_stream)
        return
    end
    lock(w.msg_lock) do
        w.gcflag || return # early exit if someone else got to this
        w.gcflag = false
        msgs = w.add_msgs
        w.add_msgs = Any[]
        if !isempty(msgs)
            remote_do(add_clients, w, msgs)
        end

        msgs = w.del_msgs
        w.del_msgs = Any[]
        if !isempty(msgs)
            remote_do(del_clients, w, msgs)
        end
    end
end

# Boundary inserted between messages on the wire, used for recovering
# from deserialization errors. Picked arbitrarily.
# A size of 10 bytes indicates ~ ~1e24 possible boundaries, so chance of collision
# with message contents is negligible.
const MSG_BOUNDARY = UInt8[0x79, 0x8e, 0x8e, 0xf5, 0x6e, 0x9b, 0x2e, 0x97, 0xd5, 0x7d]

# Faster serialization/deserialization of MsgHeader and RRID
function serialize_hdr_raw(io, hdr)
    write(io, hdr.response_oid.whence, hdr.response_oid.id, hdr.notify_oid.whence, hdr.notify_oid.id)
end

function deserialize_hdr_raw(io)
    data = read!(io, Ref{NTuple{4,Int}}())[]
    return MsgHeader(RRID(data[1], data[2]), RRID(data[3], data[4]))
end

function send_msg_(w::Worker, header, msg, now::Bool)
    check_worker_state(w)
    if myid() != 1 && !isa(msg, IdentifySocketMsg) && !isa(msg, IdentifySocketAckMsg)
        wait(w.initialized)
    end
    io = w.w_stream
    lock(io)
    try
        reset_state(w.w_serializer)
        serialize_hdr_raw(io, header)
        invokelatest(serialize_msg, w.w_serializer, msg)  # io is wrapped in w_serializer
        write(io, MSG_BOUNDARY)

        if !now && w.gcflag
            flush_gc_msgs(w)
        else
            flush(io)
        end
    finally
        unlock(io)
    end
end

function flush_gc_msgs()
    try
        for w in (PGRP::ProcessGroup).workers
            if isa(w,Worker) && w.gcflag && (w.state == W_CONNECTED)
                flush_gc_msgs(w)
            end
        end
    catch e
        bt = catch_backtrace()
        @async showerror(stderr, e, bt)
    end
end

function send_connection_hdr(w::Worker, cookie=true)
    # For a connection initiated from the remote side to us, we only send the version,
    # else when we initiate a connection we first send the cookie followed by our version.
    # The remote side validates the cookie.
    if cookie
        write(w.w_stream, LPROC.cookie)
    end
    write(w.w_stream, rpad(VERSION_STRING, HDR_VERSION_LEN)[1:HDR_VERSION_LEN])
end

max_stars_repo_name: barucden/julia

max_issues_repo_name: barucden/julia

max_forks_repo_name: barucden/julia


text:
import Pkg
Pkg.instantiate() # downloads all dependencies for the current project
Pkg.add("JSON")

using JSON

# dictionary to write
dictionary1 = Dict(
                "James" => 9001,
                "Jo" => 3474,
                "Jess" =>11926,
        )

# pass data as a json string (how it shall be displayed in a file)
stringdata = JSON.json(dictionary1)

# write the file with the stringdata variable information
open("write_read-dictionary1.json", "w") do f
        write(f, stringdata)
     end

println(stringdata)
println(dictionary1)

dictionary2 = Dict(
                "James" => 9001,
                "Jo" => 3474,
                "Jess" =>11926,
                Dict(
                        "Darko" => 11,
                        "Moljko" => 42,

                )
        )

stringdata = JSON.json(dictionary2)

open("write_read-dictionary2.json", "w") do f
        write(f, stringdata)
end


text:
immutable Beta <: ContinuousUnivariateDistribution
    alpha::Float64
    beta::Float64
    function Beta(a::Real, b::Real)
        if a > 0.0 && b > 0.0
            new(float64(a), float64(b))
        else
            error("Both alpha and beta must be positive")
        end
    end
end

Beta(a::Real) = Beta(a, a) # symmetric in [0, 1]
Beta() = Beta(1.0) # uniform

@_jl_dist_2p Beta beta

function entropy(d::Beta)
    o = lbeta(d.alpha, d.beta)
    o -= (d.alpha - 1.0) * digamma(d.alpha)
    o -= (d.beta - 1.0) * digamma(d.beta)
    o += (d.alpha + d.beta - 2.0) * digamma(d.alpha + d.beta)
    return o
end

insupport(d::Beta, x::Number) = isreal(x) && 0.0 < x < 1.0

function kurtosis(d::Beta)
    a, b = d.alpha, d.beta
    d = 6.0 * ((a - b)^2 * (a + b + 1.0) - a * b * (a + b + 2.0))
    n = a * b * (a + b + 2.0) * (a + b + 3.0)
    return d / n
end

mean(d::Beta) = d.alpha / (d.alpha + d.beta)

function modes(d::Beta)
    if d.alpha > 1.0 && d.beta > 1.0
        return [(d.alpha - 1.0) / (d.alpha + d.beta - 2.0)]
    else
        error("Beta distribution with a <= 1 || b <= 1 has no modes")
    end
end

function rand(d::Beta)
    u = rand(Gamma(d.alpha))
    return u / (u + rand(Gamma(d.beta)))
end

# TODO: Don't create temporaries here
function rand(d::Beta, dims::Dims)
    u = rand(Gamma(d.alpha), dims)
    return u ./ (u + rand(Gamma(d.beta), dims))
end

function rand!(d::Beta, A::Array{Float64})
    for i in 1:length(A)
        A[i] = rand(d)
    end
    return A
end

function skewness(d::Beta)
    d = 2.0 * (d.beta - d.alpha) * sqrt(d.alpha + d.beta + 1.0)
    n = (d.alpha + d.beta + 2.0) * sqrt(d.alpha * d.beta)
    return d / n
end

function var(d::Beta)
    ab = d.alpha + d.beta
    return d.alpha * d.beta / (ab * ab * (ab + 1.0))
end


text:
"""
    GaussianLikelihood(σ²)

Gaussian likelihood with `σ²` variance. This is to be used if we assume that the
uncertainity associated with the data follows a Gaussian distribution.

```math
    p(y|f) = Normal(y | f, σ²)
```
On calling, this would return a normal distribution with mean `f` and variance σ².
"""
struct GaussianLikelihood{T<:Real}
    σ²::Vector{T}
end

GaussianLikelihood() = GaussianLikelihood(1e-6)

GaussianLikelihood(σ²::Real) = GaussianLikelihood([σ²])

@functor GaussianLikelihood

(l::GaussianLikelihood)(f::Real) = Normal(f, sqrt(first(l.σ²)))

(l::GaussianLikelihood)(fs::AbstractVector{<:Real}) = MvNormal(fs, first(l.σ²) * I)

"""
    HeteroscedasticGaussianLikelihood(l::AbstractLink=ExpLink())

Heteroscedastic Gaussian likelihood.
This is a Gaussian likelihood whose mean and variance are functions of
latent processes.

```math
    p(y|[f, g]) = Normal(y | f, sqrt(l(g)))
```
On calling, this would return a normal distribution with mean `f` and variance `l(g)`.
Where `l` is link going from R to R^+
"""
struct HeteroscedasticGaussianLikelihood{Tl<:AbstractLink}
    invlink::Tl
end

HeteroscedasticGaussianLikelihood() = HeteroscedasticGaussianLikelihood(ExpLink())

function (l::HeteroscedasticGaussianLikelihood)(f::AbstractVector{<:Real})
    return Normal(f[1], sqrt(l.invlink(f[2])))
end

function (l::HeteroscedasticGaussianLikelihood)(fs::AbstractVector)
    return MvNormal(first.(fs), Diagonal(l.invlink.(last.(fs))))
end


text:
include("src/BioEnergeticFoodWebs.jl")
using BioEnergeticFoodWebs
using Plots

S = 15
A = nichemodel(S, 0.2)
parameters = model_parameters(A)
bm = rand(Float64, S).*(1/S)

@time sim = simulate(parameters, bm; stop=500, use=:stiff)

minimum(sim[:B])

pl2 = vec(sum(sim[:B].>0.0, dims = 2)) |> x -> plot(sim[:t], x, leg=false)
pl1 = plot(sim[:t], sim[:B], leg=false, c=:teal)
plot(pl1, pl2)

max_stars_repo_name: SimmonsBI/BioEnergeticFoodWebs

max_issues_repo_name: SimmonsBI/BioEnergeticFoodWebs

max_forks_repo_name: SimmonsBI/BioEnergeticFoodWebs


text:
using Plots
import Base.:+
import Base.:-
import Base.:*
import Base.:/
import Base.:^

"""
     Discontinuous

 A type that contains a discontinuous function and an index of any of its discontinuities. This allows the construction of step parameter profiles that retain discontinuity information. You can also call the Discontinuous type like a normal function, if you want to.
"""
struct Discontinuous <: Function
    f::Function
    d::Set
end

arithmetics = (:+, :-, :*, :/, :^)
for a ∈ arithmetics
    eval(quote
        ($a)(c::Real, D1::Discontinuous) = Discontinuous(x -> ($a)(c, D1.f(x)), D1.d)
        ($a)(D1::Discontinuous, c::Real, ) = Discontinuous(x -> ($a)(D1.f(x), c), D1.d)
        ($a)(D1::Discontinuous, D2::Discontinuous) = Discontinuous((x...) -> ($a)(D1.f(x...), D2.f(x...)), union(D1.d, D2.d))
        ($a)(f::Function, D::Discontinuous) = Discontinuous((x...) -> ($a)(f(x...), D.f(x...)), D.d)
        ($a)(D::Discontinuous, f::Function) = Discontinuous((x...) -> ($a)(D.f(x...), f(x...)), D.d)
    end)
end

# Overload call so that you can use a Discontinuous like a normal (vectorised) function, if you want
(D::Discontinuous)(x::Real) = D.f(x)
(D::Discontinuous)(x::Union{Array, StepRange, StepRangeLen, UnitRange}) = D.f.(x)

function Plots.plot(xr, D::Discontinuous; kwargs...)
    Plots.plot(xr, D(xr); kwargs...)
end
function Plots.plot(D::Discontinuous; xguide="t", yguide="p(t)", kwargs...) # Actually this can be a recipe?
    xr = extrema(collect(D.d))
    scale = abs(-(xr...))
    if scale == 0.0
        xr = collect(D.d)[1]-1:0.001:collect(D.d)[1]+1
    else
        xr = (xr[1] - 0.2*scale):(scale/1000):(xr[2] + 0.2*scale)
    end
    Plots.plot(xr, D; xguide, yguide, kwargs...)
end

# ! To use stiff solvers that error about instability with discontinuous profiles, try e.g. Rosenbrock23(autodiff = false) and turning off adaptive timestepping


text:
using NiLang, NiLang.AD
using KernelAbstractions, CUDA

@i @kernel function kernel_f(A, B::AbstractVector{TB}) where TB
    # turng off reversibility check, since GPU can not handle errors
    @invcheckoff begin
        # allocate
        batch ← @index(Global)
        s ← zero(TB)
        # computing
        for i in axes(A, 1)
            s += A[i, i, batch]
        end
        B[batch] += s
        # deallocate safely
        s → zero(TB)
        batch → @index(Global)
    end
end

@i function batched_tr!(A::CuArray{T, 3}, B::CuVector{T}) where T
    @launchkernel CUDADevice() 256 length(B) kernel_f(A, B)
end

A = CuArray(randn(ComplexF32, 10, 10, 100))
B = CUDA.zeros(ComplexF32, 100)
A_out, B_out = batched_tr!(A, B)
# put random values in the gradient field of B
grad_B = CuArray(randn(ComplexF32, 100))
A_with_g, B_with_g = (~batched_tr!)(GVar(A_out), GVar(B_out, grad_B))
# will see nonzero gradients in complex diagonal parts of A
grad_A = grad(A_with_g |> Array)


text:
# SPGL

export spgl1, project, SpotFunForward, snr

"""
# INFO
    Use:
        x, r, g, info = spgl1(A::AbstractArray, b::AbstractVector{ETb};
                            x::AbstractVector{ETx}=Array{ETb,1}(),
                            tau::AbstractFloat=NaN,
                            sigma::AbstractFloat=NaN,
                            options::spgOptions = spgOptions(),
                            params::Dict{String,Any} = Dict{String,Any}())\n
        Solve regularized composite programs, including:
            a) basis pursuit, basis pursuit denoise, and lasso
            b) non-linear versions of above problems, where forward model is non-linear

# Inputs
    A:
        - An explicit 'm' x 'n' matrix
        - An implicit JOLI opperator (in development)
        - A non-linear function handle (in development)\n
    b:
        - An m-vector\n
    tau:
        - A non-negative scalar\n
    sigma:
        - If sigma != NaN then GenSPGL will launch into a root-finding mode to find the
            tau above that solves. In this case it is STRONGLY recommended that tau = 0.\n
    x:
        - An n-vector estimate of the solution (possibly all zeros). If empty, then
            GenSPGL determines the legnth n via n = length(A'b) and sets x0 = zeros(n).\n
    options:
        - An instance of the composite tpe spgOptions. If not specified, default values
            are used. See the spgOptions documents for default values.

# Outputs
    x:
        - A solution of the problem\n
    r:
        - The residual, r = b - f(x)\n
    g:
        - The gradient, g = \u2207h(b - f(x))\n
    info:
        - An instance of the spgInfo composite type. See the spgInfo documents for
            more information.
# Author
    Keegan Lensink
        Seismic Laboratory for Imaging and Modeling
        The University of British Columbia
        keeganlensink@gmail.com

    This code is an adaptation of Michael P. Friedlander, Ewout van den Berg,
    and Aleksandr Aravkin's MATLAB program SPGL1.

"""
function spgl1(A::TA,
               b::AbstractVector{ETb};
               x::AbstractVector{ETx}=Array{ETb,1}(),
               tau::AbstractFloat=NaN,
               sigma::AbstractFloat=NaN,
               options::spgOptions = spgOptions(),
               params::Dict{String,Any} = Dict{String,Any}()) where
                 {TA<:Union{joAbstractLinearOperator,AbstractArray}, ETx<:Number, ETb<:Number}

    REVISION = "0.1"
    DATE = "June, 2017"

    #DEVNOTE# Could make Tau and Sigma nullable types? However, as long as
    # Tau and Sigma are always Float64 this wont be a problem
    (options.verbosity > 1) && println("Script made it to spgl1 for A::AbstractArray")

    # removed for 0.7+ #tic()
    m = length(b)

    #Add options proxy to params dict
    params["proxy"] = options.proxy

    # Check Tau and Sigma
    if isnan(tau) & isnan(sigma)
        tau = 0.
        sigma = 0.
        singleTau = false
    elseif isnan(sigma)
        singleTau = true
    else
        if isnan(tau)
            tau = 0.
        end
        singleTau = false
    end

    # Definitely dont do subspacemin in the non LS case
    #DEVNOTE# Make sure name matches once funLS is written
    string(options.funPenalty)=="GenSPGL.funLS" || (options.subspaceMin = 0)

    # Threshold for signifigant Newton step
    pivTol = 1e-12

    # Account for precision of data format in line search
    if options.stepMin < 10*eps(real(ETx))
        options.stepMin = 10 * eps(real(ETx))
        println("WARNING: options.stepMin is below the precision of the data type. Setting to 10*eps(DT)")
    end

    ##--------------------------------------------------------------------------------
    # Initialize Local Variables
    ##--------------------------------------------------------------------------------
    iter = 0; itnTotLSQR = 0#Total SPGL1 and LSQR iterations.
    nProdA = 0; nProdAt = 0
    lastFv = [-Inf for i=1:options.nPrevVals] # Last m functions values
    nLineTot = 0            # Total number of linesearch steps
    printTau = false
    nNewton = 0;
    bNorm, b_normalized = options.funPenalty(b, params)
    stat = false
    timeProject = zero(Float64)
    timeMatProd = zero(Float64)
    nnzIter = zero(Int64) # No. of Its with fixed pattern
    nnzIdx = BitArray{1}()  # Active set indicator
    subspace = false        # Flag if did subspace min in current itn
    stepG = 1.0               # Step length for projected gradient
    testUpdateTau = false

    ##-------------------------------------------------------------------------------
    # End Init
    ##-------------------------------------------------------------------------------

    #DEVNOTE# This could be a splitting point for multiple dispatch

    # Determine Initial x, vector length n, and check if complex
    # Explicit Method
    #DEVNOTE# Change name to JOLI  once things are working
    funForward = SpotFunForward
    options.linear = true

    if isempty(x)
        if eltype(A)<:Number
            n = size(A,2)
            realx = isreal(A) & isreal(b)
        else
            x = funForward(A, x, -b) #DEVNOTE# Check to make sure this is tmp
            n = length(x)
            realx = isreal(x) & isreal(b)
        end
        x = zeros(ETb, n)
    else
        n = length(x)
        realx = isreal(x) & isreal(A)
    end

    (eltype(A)<:Number) && (realx = realx & isreal(A))

    # Override if complex flag was used
    isnull(options.iscomplex) || (realx = ~get(options.iscomplex))

    # Check if all weights (if any) are strictly positive. In previous
    # versions we also checked if the number of weights was equal to
    # n. In the case of multiple measurement vectors, this no longer
    # needs to apply, so the check was removed.
    any(isinf.(options.weights)) && error("Weights must be finite")
    any(options.weights .> 0) || error("Weights must be strictly positive")

    # Quick exit if sigma >= ||b||. Set Tau = 0 to short circuit the loop
    if bNorm <= sigma
        println("W: sigma >= ||b||.  Exact solution is x = 0.")
        tau = zero(real(ETb))
        singleTau = true
    end

    # Don't do subspaceMin if x is complex
    if (~realx & options.subspaceMin)
        println("W: Subspace minimization disabled when variables are complex")
        subspaceMin = false
    end

    #DEVNOTE# Once these are being updated check to see if type should be inferred
            # from promotion rules
    # Pre-Allocate iteration info vectors
    xNorm1 = zeros(Float64, min(options.iterations,10000))
    rNorm2 = zeros(Float64, min(options.iterations,10000))
    lambda = zeros(Float64, min(options.iterations,10000))

    # Create ExitCondition with null trigger
    exit_status = spgExitCondition()

    #Prepare Log Header
    logheader_1 = """
    ================================================================================
    GenSPGL.jl, Rev $(REVISION), $(DATE)
    ================================================================================\n

    No. Rows                :$(m)
    No. Columns             :$(n)
    Initial Tau             :$(tau)
    Penalty                 :$(options.funPenalty)
    Regularizer             :$(options.primal_norm)
    Penalty(b)              :$(bNorm)
    Optimality tol          :$(options.optTol)
    Basis Pursuit tol       :$(options.bpTol)
    Maximum Iterations      :$(options.iterations)
    """

    options.verbosity > 1 && println("singleTau: ", singleTau)

    if singleTau
        (logheader_2 = """
    Target reg. Norm of x   :$(tau)\n

    Iter    Objective      Relative_Error  gNorm         stepG
    ------------------------------------------------------------------------
    """)

    else
        (logheader_2 = """
    Target Objective        :$(sigma)\n

    Iter   Objective      Relative_Error  RelError  gNorm        stepG  tau
    ----------------------------------------------------------------------------------------
    """)
    end

    logheader = logheader_1*logheader_2

    # Decide what to do with log header based on verbosity setting
    (options.verbosity > 0) && println(logheader)

    # Project the stating point and evaluate function and gradient
    r = typeof(b)()

    if isempty(x)#DEVNOTE# matlab Legacy, this can never be invoked. see 90-99

        #DEVNOTE# Why copy? Waste of mem and time. Check if b is used again
        r = deepcopy(b)
        f,g,g2 = options.funCompositeR(A, x, r, funForward, options.funPenalty, timeMatProd, nProdAt)
    else
        x,itn = project(x,tau, timeProject, options, params)
        r = b - funForward(A, x, [], params)
        nProdA += 1
        f,g,g2 = options.funCompositeR(A, x,r, funForward, options.funPenalty, nProdAt, params)
        dx_tmp, itn_tmp = project(x-g, tau, timeProject, options, params)
        dx = dx_tmp - x
        itn += itn_tmp
    end

    dxNorm = norm(dx,Inf)
    if dxNorm < (1/options.stepMax)
        gStep = options.stepMax
    else
        gStep = min(options.stepMax, max(options.stepMin, 1/dxNorm))
    end

    # Required for non-monotone strategy
    lastFv[1] = copy(f)
    fBest = copy(f)
    xBest = copy(x)
    fOld = copy(f)

    (options.verbosity > 1) && println("Init Finished")

    # Wrap up initialized variables
    init = spgInit(A,
                    b,
                    x,
                    tau,
                    sigma,
                    g,
                    g2,
                    f,
                    r,
                    nnzIdx,
                    nnzIter,
                    options,
                    params,
                    timeProject,
                    timeMatProd,
                    exit_status,
                    singleTau,
                    bNorm,
                    fOld,
                    testUpdateTau,
                    iter,
                    nNewton,
                    printTau,
                    subspace,
                    stepG,
                    xNorm1,
                    rNorm2,
                    lambda,
                    lastFv,
                    gStep,
                    funForward,
                    nLineTot,
                    nProdAt,
                    nProdA,
                    fBest,
                    xBest)

    # Wrap main loop in a function to ease type stability
    init, rNorm, gNorm, rErr  = spglcore(init)

    # Prepare output
    info = spgInfo( init.tau,
                    rNorm,
                    gNorm,
                    rErr,
                    init.exit_status,
                    init.iter,
                    init.nProdA,
                    init.nProdAt,
                    init.nNewton,
                    init.timeProject,
                    init.timeMatProd,
                    init.options,
                    init.xNorm1[1:init.iter],
                    init.rNorm2[1:init.iter],
                    init.lambda[1:init.iter])

    if (options.verbosity > 0)
        println("-----------------------------------------------------------------------\n")
        print(info.exit_status)
    end

    return init.x, init.r, init.g, info
end #func

"""
Use:    x = project(x::AbstractArray, tau::Number, timeProject::Float64
                    options::spgOptions, params::Dict)
"""
function project(x::Tx, tau::Number, timeProject::Float64,
                    options::spgOptions, params::Dict) where
                      {ETx<:Number, Tx<:AbstractVector{ETx}}

    (options.verbosity > 1) && println("Begin Project")

    x_out::Tx, itn::Int64 = options.project(x, tau, options.weights, params)

    #DEVNOTE# Replace with @elapsed at call #timeProject[1] += (toc() - tStart)

    (options.verbosity > 1) && println("Finish Project")

    return x_out, itn

end

"""
Activated when an explicit or JOLI operator is passed in.
"""
function SpotFunForward(A::TA, x::AbstractArray, g::AbstractArray, params::Dict) where
                       {TA<:Union{joAbstractLinearOperator,AbstractArray}}
    #DEVNOTE# Double check type of g once in use

    isempty(g) && (f = A*x)
    isempty(g) || (f = A'*g)

    return f

end

"""
    Use:
        x, r, g, info = spgl1(A::Function, b::AbstractVector{ETb};
                            x::AbstractVector{ETx}=Array{ETb,1}(),
                            tau::AbstractFloat=NaN,
                            sigma::AbstractFloat=NaN,
                            options::spgOptions = spgOptions(),
                            params::Dict{String,Any} = Dict{String,Any}())\n
"""
function spgl1(A::Function, b::AbstractVector{ETb};
               x::AbstractVector{ETx}=Array{ETb,1}(),
               tau::AbstractFloat=NaN,
               sigma::AbstractFloat=NaN,
               options::spgOptions = spgOptions(),
               params::Dict{String,Any} = Dict{String,Any}()) where
                  {ETx<:Number, ETb<:Number}

    REVISION = "0.1"
    DATE = "June, 2017"

    #DEVNOTE# Could make Tau and Sigma nullable types? However, as long as
    # Tau and Sigma are always Float64 this wont be a problem
    (options.verbosity > 1) && println("Script made it to spgl1 for A::Function")

    # removed for 0.7+ #tic()
    m = length(b)

    #Add options proxy to params dict
    params["proxy"] = options.proxy

    # Check Tau and Sigma
    if isnan(tau) & isnan(sigma)
        tau = 0.
        sigma = 0.
        singleTau = false
    elseif isnan(sigma)
        singleTau = true
    else
        if isnan(tau)
            tau = 0.
        end
        singleTau = false
    end

    # Definitely dont do subspacemin in the non LS case
    #DEVNOTE# Make sure name matches once funLS is written
    string(options.funPenalty)=="GenSPGL.funLS" || (options.subspaceMin = 0)

    # Threshold for signifigant Newton step
    pivTol = 1e-12

    # Account for precision of data format in line search
    if options.stepMin < eps(real(ETx))
        options.stepMin = 10 * eps(real(ETx))
        println("WARNING: options.stepMin is below the precision of the data type. Setting to 10*eps(DT)")
    end

    ##--------------------------------------------------------------------------------
    # Initialize Local Variables
    ##--------------------------------------------------------------------------------
    iter = 0; itnTotLSQR = 0#Total SPGL1 and LSQR iterations.
    nProdA = 0; nProdAt = 0
    lastFv = [-Inf for i=1:options.nPrevVals] # Last m functions values
    nLineTot = 0            # Total number of linesearch steps
    printTau = false
    nNewton = 0;
    bNorm, b_normalized = options.funPenalty(b, params)
    stat = false
    timeProject = zero(Float64)
    timeMatProd = zero(Float64)
    nnzIter = zero(Int64) # No. of Its with fixed pattern
    nnzIdx = BitArray{1}()  # Active set indicator
    subspace = false        # Flag if did subspace min in current itn
    stepG = 1.0               # Step length for projected gradient
    testUpdateTau = false

    ##-------------------------------------------------------------------------------
    # End Init
    ##-------------------------------------------------------------------------------

    #DEVNOTE# This could be a splitting point for multiple dispatch

    # Determine Initial x, vector length n, and check if complex
    # Explicit Method
    funForward = A

    if isempty(x)
        x_tmp = funForward(A, x, -b, params)
        n = length(x_tmp)
        realx = isreal(x_tmp) & isreal(b)

        x = zeros(n)
    else
        n = length(x)
        realx = isreal(x) & isreal(b)
    end

    (eltype(A)<:Number) && (realx = realx & isreal(A))

    # Override if complex flag was used
    isnull(options.iscomplex) || (realx = ~get(options.iscomplex))

    # Check if all weights (if any) are strictly positive. In previous
    # versions we also checked if the number of weights was equal to
    # n. In the case of multiple measurement vectors, this no longer
    # needs to apply, so the check was removed.
    any(isinf.(options.weights)) && error("Weights must be finite")
    any(options.weights .> 0) || error("Weights must be strictly positive")

    # Quick exit if sigma >= ||b||. Set Tau = 0 to short circuit the loop
    if bNorm <= sigma
        println("W: sigma >= ||b||.  Exact solution is x = 0.")
        tau = 0
        singleTau = true
    end

    # Don't do subspaceMin if x is complex
    if (~realx & options.subspaceMin)
        println("W: Subspace minimization disabled when variables are complex")
        subspaceMin = false
    end

    #DEVNOTE# Once these are being updated check to see if type should be inferred
            # from promotion rules
    # Pre-Allocate iteration info vectors
    xNorm1 = zeros(Float64, min(options.iterations,10000))
    rNorm2 = zeros(Float64, min(options.iterations,10000))
    lambda = zeros(Float64, min(options.iterations,10000))

    # Create ExitCondition with null trigger
    exit_status = spgExitCondition()

    #Prepare Log Header
    logheader_1 = """
    ================================================================================
    GenSPGL.jl, Rev $(REVISION), $(DATE)
    ================================================================================\n

    No. Rows                :$(m)
    No. Columns             :$(n)
    Initial Tau             :$(tau)
    Penalty                 :$(options.funPenalty)
    Regularizer             :$(options.primal_norm)
    Penalty(b)              :$(bNorm)
    Optimality tol          :$(options.optTol)
    Basis Pursuit tol       :$(options.bpTol)
    Maximum Iterations      :$(options.iterations)
    """

    options.verbosity > 1 && println("singleTau: ", singleTau)

    if singleTau
        (logheader_2 = """
    Target reg. Norm of x   :$(tau)\n

    Iter    Objective      Relative_Error  gNorm         stepG
    -------------------------------------------------------------------
    """)

    else
        (logheader_2 = """
    Target Objective        :$(sigma)\n

    Iter   Objective      Relative_Error  RelError  gNorm         stepG  tau
    -------------------------------------------------------------------------
    """)
    end

    logheader = logheader_1*logheader_2

    # Decide what to do with log header based on verbosity setting
    (options.verbosity > 0) && println(logheader)

    # Project the stating point and evaluate function and gradient
    r = typeof(b)()

    if isempty(x)#DEVNOTE# matlab Legacy, this can never be invoked. see 90-99

        #DEVNOTE# Why copy? Waste of mem and time. Check if b is used again
        r = deepcopy(b)
        f,g,g2 = options.funCompositeR(A, x, r, funForward, options.funPenalty, timeMatProd, nProdAt)
    else
        x,itn = project(x,tau, timeProject, options, params)
        r = b - funForward(A, x, Array{ETx,1}(), params)[1]
        nProdA += 1
        f,g,g2 = options.funCompositeR(A, x,r, funForward, options.funPenalty, nProdAt, params)
        dx_tmp, itn_tmp = project(x-g, tau, timeProject, options, params)
        dx = dx_tmp - x
        itn += itn_tmp
    end

    dxNorm = norm(dx,Inf)
    if dxNorm < (1/options.stepMax)
        gStep = options.stepMax
    else
        gStep = min(options.stepMax, max(options.stepMin, 1/dxNorm))
    end

    # Required for non-monotone strategy
    lastFv[1] = copy(f)
    fBest = copy(f)
    xBest = copy(x)
    fOld = copy(f)

    (options.verbosity > 1) && println("Init Finished")

    # Wrap up initialized variables
    init = spgInit(A,
                    b,
                    x,
                    tau,
                    sigma,
                    g,
                    g2,
                    f,
                    r,
                    nnzIdx,
                    nnzIter,
                    options,
                    params,
                    timeProject,
                    timeMatProd,
                    exit_status,
                    singleTau,
                    bNorm,
                    fOld,
                    testUpdateTau,
                    iter,
                    nNewton,
                    printTau,
                    subspace,
                    stepG,
                    xNorm1,
                    rNorm2,
                    lambda,
                    lastFv,
                    gStep,
                    funForward,
                    nLineTot,
                    nProdAt,
                    nProdA,
                    fBest,
                    xBest)

    # Wrap main loop in a function to ease type stability
    init, rNorm, gNorm, rErr  = spglcore(init)
    # Prepare output
    info = spgInfo(  init.tau,
                    rNorm,
                    gNorm,
                    rErr,
                    init.exit_status,
                    init.iter,
                    init.nProdA,
                    init.nProdAt,
                    init.nNewton,
                    init.timeProject,
                    init.timeMatProd,
                    init.options,
                    init.xNorm1[1:init.iter],
                    init.rNorm2[1:init.iter],
                    init.lambda[1:init.iter])

    if (options.verbosity > 0)
        println("---------------------------------------------------------------------------\n")
        print(info.exit_status)
    end

    return init.x, init.r, init.g, info
end #func

snr(raw,interp) = -20log10(norm(interp-raw)/norm(raw))


text:
using Test
using OpenStreetMapX
using SmartTransitionSim
using Random

test_map = OpenStreetMapX.get_map_data("reno_east3.osm", use_cache = false)
Rect1 = [Rect((39.50,-119.70),(39.55,-119.74))]
Rect2 = [Rect((39.50,-119.80),(39.55,-119.76))]
KDict = Dict{Tuple{Int,Int},Array{Vector{Int}}}()
AgentsSet = generate_agents(test_map,10,Rect1,Rect2, 0.5, 3, 1.0, 1.0, KDict)

#generate_agents.jl
@testset "agents" begin

Random.seed!(0);
@test rand(Int) == -4635026124992869592

@test in(pick_random_node(test_map, Rect1), keys(test_map.nodes))
@test in(pick_random_node(test_map, [Rect1[1], Rect2[1]]), keys(test_map.nodes))
@test typeof(get_nodes_set(test_map, Rect1)) == Array{Int64,1}

@test all([in(x, keys(test_map.nodes)) for x in getfield.(AgentsSet,:start_node)])
@test sum(getfield.(AgentsSet,:smart)) == 5

end

#rerouting.jl
@testset "rerouting" begin

constantNode = AgentsSet[1].route[1]
speeds = OpenStreetMapX.get_velocities(test_map)
k_shortest_path_rerouting!(test_map, KDict, AgentsSet[1], speeds, speeds, 3, 1.0, 100)
@test AgentsSet[1].route[1] == constantNode

k_routes = yen_a_star(test_map.g,test_map.v[AgentsSet[1].route[1]],test_map.v[AgentsSet[1].route[end]],speeds,3)
@test length(k_routes.dists) == 3
end

#simulations.jl
@testset "simulations" begin

newAgents = generate_agents(test_map,10,Rect1,Rect2, 0.5, 3, 1.0,1.0, KDict)
output = simulation_run(:base,test_map, newAgents)
@test length(output) == 3
@test typeof(output) == NamedTuple{(:Steps, :Simtime, :TravelTimes),Tuple{Int64,Float64,Array{Float64,1}}}

ITSOutput = simulation_run(:smart,test_map, newAgents, 5.0, KDict, 100, 1.0, 3, true)

@test length(ITSOutput) == 3
@test typeof(ITSOutput) == NamedTuple{(:Steps, :Simtime, :TravelTimes),Tuple{Int64,Float64,Array{Float64,1}}}

output = simulation_run(:base,test_map, AgentsSet)
stats = gather_statistics(getfield.(AgentsSet,:smart),
                    output.TravelTimes,
                    ITSOutput.TravelTimes)

@test length(stats) == 7
@test typeof(stats) == NamedTuple{(:overall_time, :smart_time, :other_time, :avg_base, :avg_overall_V2I, :avg_smart_V2I, :avg_regular_V2I),
                        Tuple{Float64,Float64,Float64,Float64,Float64,Float64,Float64}}

Params = (reps = 1, α = 0.05, N = 8, U = 150, T = 1.0, k = 3, AvgT = 1.0)
ParametersGridProd = collect(Base.Iterators.product(Params.reps,Params.α, Params.N, Params.U, Params.T, Params.k, Params.AvgT))
ParametersGrid = [ParametersGridProd[i] for i in 1:length(ParametersGridProd)]

run_parameter_analysis(1,ParametersGrid,newAgents,5.0,KDict,test_map)
run_parameter_analysis(1,ParametersGrid,Rect1,Rect2,5.0,KDict,test_map)

end

#traffic_model.jl
@testset "traffic_model" begin

max_dens = get_max_densities(test_map, 5.0)
max_speeds = OpenStreetMapX.get_velocities(test_map)
AgentsSet[1].start_time = 0.0
AgentsSet[1].active = true
densities, speeds = init_traffic_variables(test_map, AgentsSet)
#get_max_densities tests
@test round(sum(max_dens);digits=3) == 206315.887

edge = rand(keys(densities))
update_weights!(speeds, densities, max_dens, max_speeds)
#update_weights! tests
@test speeds[edge[1],edge[2]] < max_speeds[edge[1],edge[2]] && speeds[487,1] == max_speeds[487,1]

#traffic_constants tests
max_d, max_s = traffic_constants(test_map, 5.0)
@test max_d == max_dens && max_s == max_speeds

#init_traffic_variables tests
i_densities, i_speeds = init_traffic_variables(test_map, AgentsSet)
@test sum(values(i_densities)) == 1
@test i_speeds == max_speeds

#next_edge tests
events = next_edge(AgentsSet, max_speeds, test_map.w)
@test typeof(events) == Vector{Float64}
event = next_edge(AgentsSet[1], max_speeds, test_map.w)
@test typeof(event) == Float64

#update_event_agent! tests
Agent1 = deepcopy(AgentsSet[1])
update_event_agent!(AgentsSet[1],event[1], i_densities, test_map.v)
@test AgentsSet[1].edge[1] == Agent1.edge[2]
@test length(AgentsSet[1].route) == length(Agent1.route) - 1

end

#utilities.jl
@testset "utilities" begin
    #Check if print_map_statistics() is not returning errors
    @test print_map_statistics("./","reno_east3.osm",5)
end


text:
module Dirac_operators

import Gaugefields:AbstractGaugefields,Abstractfields,CovNeuralnet,shift_U,clear_U!,set_wing_U!,add_U!
#import Gaugefields:Verbose_level,Verbose_3,Verbose_2,Verbose_1,println_verbose3
#import Gaugefields:Verbose_level,Verbose_3,Verbose_2,Verbose_1,println_verbose3
using LinearAlgebra
import Gaugefields.Verboseprint_mpi:Verbose_print,println_verbose_level1,println_verbose_level2,println_verbose_level3
using SparseArrays

include("./cgmethods.jl")
include("./SakuraiSugiura/SakuraiSugiuramethod.jl")

abstract type Operator#<: AbstractMatrix{ComplexF64}
end

abstract type Dirac_operator{Dim}  <: Operator
end

abstract type DdagD_operator  <: Operator
end

abstract type γ5D_operator  <: Operator #hermitian matrix
end

struct γ5D{Dirac} <: γ5D_operator
    dirac::Dirac
    function γ5D(D)
        return new{typeof(D)}(D)
    end
end

abstract type Adjoint_Dirac_operator <: Operator
end

function Base.size(D::Operator)
    error("Base.size(D::T)  is not implemented in type $(typeof(D))")
end

function Base.size(D::Dirac_operator)
    x = get_temporaryvectors_forCG(D)[1]
    NN = length(x)
    return (NN,NN)
end

function Base.size(D::DdagD_operator)
    return size(D.dirac)
end

function Base.size(D::γ5D_operator)
    return size(D.dirac)
end

function Base.adjoint(A::Dirac_operator{Dim} ) where Dim
    error("Base.adjoint(A::T)  is not implemented in type $(typeof(A))")
end

Base.adjoint(A::Adjoint_Dirac_operator) = A.parent

Base.adjoint(A::γ5D_operator) = A

function get_temporaryvectors_forCG(A::T) where T <: Dirac_operator
    return A._temporary_fermion_forCG
end

function get_temporaryvectors_forCG(A::T) where T <: Adjoint_Dirac_operator
    return A.parent._temporary_fermion_forCG
end

function get_temporaryvectors_forCG(A::T) where T <: γ5D_operator
    return A.dirac._temporary_fermion_forCG
end

const default_eps_CG = 1e-19
const default_MaxCGstep = 3000

include("./AbstractFermions.jl")
include("./StaggeredFermion/StaggeredFermion.jl")
include("./WilsonFermion/WilsonFermion.jl")
include("./DomainwallFermion/DomainwallFermion.jl")
include("./MobiusDomainwallFermion/MobiusDomainwallFermion.jl")
include("./action/FermiAction.jl")

function Dirac_operator(U::Array{<: AbstractGaugefields{NC,Dim},1},x,parameters) where {NC,Dim}
    @assert haskey(parameters,"Dirac_operator") "parameters should have Dirac_operator keyword!"
    if parameters["Dirac_operator"] == "staggered"
        Staggered_Dirac_operator(U,x,parameters)
    elseif parameters["Dirac_operator"] == "Wilson"
        fasterversion = check_parameters(parameters,"faster version",false)
        if fasterversion
            Wilson_Dirac_operator_faster(U,x,parameters)
        else
            Wilson_Dirac_operator(U,x,parameters)
        end
    elseif parameters["Dirac_operator"] == "Wilson_general"
        Wilson_GeneralDirac_operator(U,x,parameters)
    elseif parameters["Dirac_operator"] == "Domainwall"
        Domainwall_Dirac_operator(U,x,parameters)
    elseif parameters["Dirac_operator"] == "MobiusDomainwall"
        MobiusDomainwall_Dirac_operator(U,x,parameters)
    else
        error("$(parameters["Dirac_operator"]) is not supported")
    end
end

function DdagD_operator(U::Array{<: AbstractGaugefields{NC,Dim},1},x,parameters) where {NC,Dim}
    @assert haskey(parameters,"Dirac_operator") "parameters should have Dirac_operator keyword!"
    if parameters["Dirac_operator"] == "staggered"
        DdagD_Staggered_operator(U,x,parameters)
    elseif parameters["Dirac_operator"] == "Wilson"
        DdagD_Wilson_operator(U,x,parameters)
    elseif parameters["Dirac_operator"] == "Domainwall"
        DdagD_Domainwall_operator(U,x,parameters)
    else
        error("$(parameters["Dirac_operator"]) is not supported")
    end
end

function get_temporaryvectors_forCG(A::T) where T <: DdagD_operator
    return A.dirac._temporary_fermion_forCG
end

function solve_DinvX!(y::T1,A::T2,x::T3) where {T1 <: AbstractFermionfields,T2 <:  Operator, T3 <: AbstractFermionfields}
    error("solve_DinvX!(y,A,x) (y = A^{-1} x) is not implemented in type y:$(typeof(y)),A:$(typeof(A)) and x:$(typeof(x))")
end

function solve_DinvX!(y::T1,A::T2,x::T3) where {T1 <: AbstractFermionfields,T2 <:  γ5D_operator, T3 <: AbstractFermionfields}
    x2 = deepcopy(x)
    apply_γ5!(x2)
    solve_DinvX!(y,A.dirac,x2)
end

function solve_DinvX!(y::T1,A::T2,x::T3) where {T1 <: AbstractFermionfields,T2 <:  Dirac_operator, T3 <: AbstractFermionfields}
    if A.method_CG == "bicg"
        bicg(y,A,x;eps=A.eps_CG,maxsteps = A.MaxCGstep,verbose = A.verbose_print)#set_verbose(A.verbose_level))
        set_wing_fermion!(y,A.boundarycondition)
    elseif A.method_CG == "bicgstab"
        bicgstab(y,A,x;eps=A.eps_CG,maxsteps = A.MaxCGstep,verbose = A.verbose_print)
        set_wing_fermion!(y,A.boundarycondition)
    elseif A.method_CG == "preconditiond_bicgstab"
        #@assert A.Dirac_operator == "Wilson" "preconditiond_bicgstab is supported only in Wilson Dirac operator"
        WW = Wilson_Dirac_operator_evenodd(A)
        #b = A._temporary_fermi[6]
        #substitute_fermion!(beff,x)
        bout = A._temporary_fermi[7]
        calc_beff!(bout,A.U,x,A)
        iseven = true
        isodd = false

        #bout = x
        #bicgstab(y,WW,bout;eps=A.eps_CG,maxsteps = A.MaxCGstep,verbose = set_verbose(A.verbose_level))
        bicgstab_evenodd(y,WW,bout,iseven;eps=A.eps_CG,maxsteps = A.MaxCGstep,verbose = A.verbose_print)
        Tx = A._temporary_fermi[6]
        set_wing_fermion!(y,A.boundarycondition,iseven)

        Toex!(Tx,A.U,y,A,iseven)
        #set_wing_fermion!(Tx,A.boundarycondition)
        add_fermion!(y,1,x,1,Tx,isodd)
        set_wing_fermion!(y,A.boundarycondition,isodd)

        #Toex!(y,U,x,A,iseven)
        #xo = K Toe xe + b0
    else
        error("A.method_CG = $(A.method_CG) is not supported")
    end

end

function solve_DinvX!(y::T1,A::T2,x::T3) where {T1 <: AbstractFermionfields,T2 <:  Adjoint_Dirac_operator, T3 <: AbstractFermionfields}
    if A.parent.method_CG == "bicg"
        bicg(y,A,x;eps=A.parent.eps_CG,maxsteps = A.parent.MaxCGstep,verbose = A.parent.verbose_print)
    elseif A.parent.method_CG == "bicgstab"
        #println("Adag")
        bicgstab(y,A,x;eps=A.parent.eps_CG,maxsteps = A.parent.MaxCGstep,verbose = A.parent.verbose_print)
    elseif A.parent.method_CG == "preconditiond_bicgstab"
        bicgstab(y,A,x;eps=A.parent.eps_CG,maxsteps = A.parent.MaxCGstep,verbose = A.parent.verbose_print)
        #=
        #@assert A.Dirac_operator == "Wilson" "preconditiond_bicgstab is supported only in Wilson Dirac operator"
        WWdag = Wilson_Dirac_operator_evenodd(A.parent)'
        #b = A._temporary_fermi[6]
        #substitute_fermion!(beff,x)
        bout = A.parent._temporary_fermi[7]
        #calc_beff!(bout,A.parent.U,x,A.parent)
        calc_beff_dag!(bout,A.parent.U,x,A.parent)
        iseven = true

        #bout = x
        #bicgstab(y,WW,bout;eps=A.eps_CG,maxsteps = A.MaxCGstep,verbose = set_verbose(A.verbose_level))
        bicgstab_evenodd(y,WW',bout,iseven;eps=A.parent.eps_CG,maxsteps = A.parent.MaxCGstep,verbose = set_verbose(A.parent.verbose_level))
        Tx = A.parent._temporary_fermi[6]
        set_wing_fermion!(y,A.parent.boundarycondition)
        Tdagoex!(Tx,A.parent.U,y,A.parent,iseven)
        #set_wing_fermion!(Tx,A.boundarycondition)
        add_fermion!(y,1,x,1,Tx,isodd)
        =#

    else
        error("A.method_CG = $(A.method_CG) is not supported")
    end

    set_wing_fermion!(y,A.parent.boundarycondition)
end

using InteractiveUtils

function solve_DinvX!(y::T1,A::T2,x::T3) where {T1 <: AbstractFermionfields,T2 <:  DdagD_operator, T3 <: AbstractFermionfields}
    cg(y,A,x;eps=A.dirac.eps_CG,maxsteps = A.dirac.MaxCGstep,verbose = A.dirac.verbose_print)
    set_wing_fermion!(y,A.dirac.boundarycondition)
end

function LinearAlgebra.mul!(y::T1,A::T2,x::T3) where {T1 <: AbstractFermionfields,T2 <:  Operator, T3 <: AbstractFermionfields}
    error("LinearAlgebra.mul!(y,A,x) is not implemented in type y:$(typeof(y)),A:$(typeof(A)) and x:$(typeof(x))")
end

function LinearAlgebra.mul!(y::AbstractFermionfields{NC,Dim},A::T,x::AbstractFermionfields{NC,Dim})  where {T <: DdagD_operator,NC,Dim} #y = A*x
    temp = A.dirac._temporary_fermi[5]

    mul!(temp,A.dirac,x)
    mul!(y,A.dirac',temp)

    return
end

function LinearAlgebra.mul!(y::AbstractFermionfields{NC,Dim},A::T,x::AbstractFermionfields{NC,Dim})  where {T <: γ5D_operator,NC,Dim} #y = A*x
    mul!(y,A.dirac,x)
    apply_γ5!(y)
    return
end

function check_parameters(parameters,key,initial)
    if haskey(parameters,key)
        value = parameters[key]
    else
        value = initial
    end
    return value
end

function check_important_parameters(parameters,key,sample = nothing)
    if sample == nothing
        errstring = ""
    else
        errstring = "sample is $sample"
    end

    @assert haskey(parameters,key) "\"parameters\" should have the keyword $key. $errstring. Now \"parameters\" have $parameters"
    return parameters[key]
end

function set_verbose(verbose_level)
    if verbose_level == 1
        verbose = Verbose_1()
    elseif verbose_level == 2
        verbose = Verbose_2()
    elseif verbose_level == 3
        verbose = Verbose_3()
    else
        error("verbose_level = $verbose_level is not supported")
    end
    return verbose
end

function construct_sparsematrix(D::Operator) # D_ij = e_i D e_j
    NN,_ = size(D)
    mat_D = spzeros(ComplexF64,NN,NN)
    temp1 = get_temporaryvectors_forCG(D)[1]
    temp2 = get_temporaryvectors_forCG(D)[2]

    for j=1:NN
        clear_fermion!(temp1)
        temp1[j] = 1
        set_wing_fermion!(temp1)
        mul!(temp2,D,temp1)
        for i=1:NN
            mat_D[i,j] = temp2[i]
        end
    end
    return mat_D
end

end


text:
##### MEAN-VARIANCE OPTIMIZATION
#This is not supposed to produce great results, since mean and covariance
#is supposed to be very exact before good weights are found

#Load packages and support functions
using JuMP
using Gurobi
using StatsBase
using DataFrames
using CSV
include("SupportFunction.jl")
include("DataLoad.jl")
include("CEOSupportFunctions.jl")
println("Leeeeroooy Jenkins")
#Esben's path
#path = "$(homedir())/Documents/GitHub/Thesis/Data/IndexData"

#Set path for dataloading
path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Thesis/Data/IndexDataDiff/"
path = "/zhome/9f/d/88706/SpecialeCode/Thesis/Data/IndexDataDiff/"

#Set parameters
trainingSize = 240
possibilities = 5
industries = ["NoDur", "Durbl", "Manuf", "Enrgy", "HiTec", "Telcm", "Shops", "Hlth", "Utils", "Other"]
industriesTotal = length(industries)

#Define the input and transformations
modelMatrix = zeros(industriesTotal, possibilities)
noDurModel = [1 0 0 1 1]
testModel = [1 0 0 1 1]
modelMatrix[1, :] = noDurModel
for i=2:industriesTotal
    modelMatrix[i, :] = noDurModel
end

#Create and fill X and Y arrays
XArrays = Array{Array{Float64, 2}}(industriesTotal)
YArrays = Array{Array{Float64, 2}}(industriesTotal)
XArrays, YArrays = generateXandYs(industries, modelMatrix)

#Set
nGammas = 5
standY = YArrays[1]
SSTO = sum((standY[i]-mean(standY[:]))^2 for i=1:length(standY))
lambdaValues = [0.1, 0.5, 1, 2, 10]
#lambda4Values = [1e4,1e5,1e6,1e7,1e8]
nRows = size(standY)[1]
amountOfModels = nGammas^3

modelConfig = zeros(amountOfModels, 4)
counter = 1
for l2 = 1:nGammas
    for l3 = 1:nGammas
        for l4 = 1:nGammas
            modelConfig[counter,:] = [1 lambdaValues[l2] lambdaValues[l3] lambdaValues[l4]]
            counter += 1
        end
    end
end
modelConfig

#Initialization of parameters
w1N = repeat([1/10], outer = 10) #1/N weights
gamma = 2.4 #risk aversion
validationPeriod = 0 # Number of periods to do parameter training
PMatrix = zeros(nRows-trainingSize-2, amountOfModels)

bestModelAmount = nGammas^3
bestModelConfigs = zeros(bestModelAmount, 4)
bestModelIndexes = zeros(bestModelAmount)

return1NMatrix = zeros(nRows-trainingSize)
returnCEOMatrix = zeros(nRows-trainingSize, bestModelAmount)
returnPerfectMatrix = zeros(nRows-trainingSize)

weightsPerfect = zeros(nRows-trainingSize, 10)
weightsCEO     = zeros(nRows-trainingSize, 10, bestModelAmount)
expectedReturnMatrix = zeros(nRows-trainingSize, 10)
forecastErrors = zeros(nRows-trainingSize, 10)

#= comment out as we run with no risk free option
rfRates = loadRiskFreeRate("NoDur", path)
rfRates = rfRates[:,1]
=#

#path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Results/CEO/RFR/VIXTimeTA2.4/"
path = "/zhome/9f/d/88706/SpecialeCode/Results/CEO/"
weightsPerfect = Array{Float64}(CSV.read(path*"weightsPerfect.csv",header=false, datarow=1, nullable=false))
returnPerfectMatrix = Array{Float64}(CSV.read(path*"returnPerfectMatrix.csv", header=false, datarow=1, nullable=false))
#=
for i = 1:5
    for j = 0:24
        #path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Results/CEO/"
        #path = "/zhome/9f/d/88706/SpecialeCode/Results/CEORFR/"
        PMatrixTemp = Array{Float64}(CSV.read(path*string(i)*"_"*string(j)*"_PMatrix.csv",header=false, datarow=1, nullable=false, types=[Float64]))
        return1NMatrixTemp = Array{Float64}(CSV.read(path*string(i)*"_"*string(j)*"_return1NMatrix.csv",header=false, datarow=1, nullable=false, types=[Float64]))
        returnCEOMatrixTemp = Array{Float64}(CSV.read(path*string(i)*"_"*string(j)*"_returnCEOMatrix.csv",header=false, datarow=1, nullable=false, types=[Float64]))

        PMatrix[i,:] += PMatrixTemp
        return1NMatrix[i] += return1NMatrixTemp[i]
        returnCEOMatrix[i,:] += returnCEOMatrixTemp
    end
end
=#
#Load weightsPerfect and return perfect matrix
#Load PMatrix

#modelMeans = mean(PMatrix, 1)
counter = 1
for l2 = 1:nGammas
    for l3 = 1:nGammas
        for l4 = 1:nGammas
            bestModelConfigs[counter,:] = [1 lambdaValues[l2] lambdaValues[l3] lambdaValues[l4]]
            counter += 1
        end
    end
end

for i = 1:bestModelAmount
	bestModelIndexes[i]=i
end

inputArg1 = 0 #83 breaks the trainingSize
inputArg1 = parse(Int64,ARGS[1])
println("Starting CEO Validation loop")
t=1+inputArg1
println("Time $t/50")
trainingXArrays, trainingYArrays, validationXRows, validationY, OOSXArrays, OOSYArrays, OOSRow, OOSY = createDataSplits(XArrays, YArrays, t, trainingSize)

for m = 1:bestModelAmount
    expectedReturns = zeros(10)
    #CHANGES
    rfRatesVec = rfRates[t:(t+trainingSize-1)]
    betaArray, U = @time(runCEO(trainingXArrays, trainingYArrays, bestModelConfigs[m, :], gamma))

    #betaArrayCopy = betaArray
    ##PREVIOUS
    #betaArray, U = @time(runCEO(trainingXArrays, trainingYArrays, modelConfig[m, :], gamma))
    expectedReturns[1:10] = generateExpectedReturns(betaArray, trainingXArrays, trainingYArrays, validationXRows)
    #expectedReturns[11] = rfRates[t+trainingSize]
    if m == 1
        expectedReturnMatrix[t, 1:10] = (exp.(expectedReturns[1:10])-1)*100
        #expectedReturnMatrix[t, 11] = (exp.(expectedReturns[11])-1)
    end
    #Need to send OOSRow to mean-variance optimization to get "perfect information" since validationY is the values in OOSRow[1:10]
    valY = zeros(10)
    for i = 1:10
        valY[i] = validationY[i][1]
    end
    #valY[11] = rfRates[t+trainingSize]

    #rfRatesVec = rfRates[t:(t+trainingSize-1)]
    trainX = trainingXArrays[1][:,1:10]
    Sigma =  cov(trainX)
    F = lufact(Sigma)
    U = F[:U]  #Cholesky factorization of Sigma

    return1N, returnCEO, wStar, forecastRow = performMVOptimization(expectedReturns, U, gamma, valY, valY)
    weightsCEO[t, 1:10, Int64(m)] = wStar
    return1NMatrix[t]      = return1N
    returnCEOMatrix[t, Int64(bestModelIndexes[m])] = returnCEO
    returnPerfect = returnPerfectMatrix[t]
    if m == 1
        forecastErrors[t, 1:10]  = forecastRow-expectedReturnMatrix[t,:]
    end
    println("1N returns is $return1N, returnPerfect is $returnPerfect and returnCEO is $returnCEO")
    PMatrix[t, Int64(bestModelIndexes[m])] = calculatePvalue(return1N, returnPerfect, returnCEO)
end

combinedPortfolios = hcat(returnPerfectMatrix[1:nRows-trainingSize-2, 1], return1NMatrix[1:nRows-trainingSize-2, 1],
    returnCEOMatrix[1:nRows-trainingSize-2, Array{Int64}(bestModelIndexes)], PMatrix[1:nRows-trainingSize-2, Array{Int64}(bestModelIndexes)])
#path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Results/CEO/"
path = "/zhome/9f/d/88706/SpecialeCode/Results/CEOSearch/Iter0/"
writedlm(path*string(inputArg1)*"_returnPvalueOutcome1to200.csv", combinedPortfolios, ",")
for i = 1:bestModelAmount
    writedlm(path*"wightsCEO_Model"*string(i)*"_"*string(inputArg1)*".csv",weightsCEO[:,:,i], ",")
end

combinedPortfolios2 = hcat(expectedReturnMatrix, forecastErrors)
writedlm(path*string(inputArg1)*"_forecastsAndErrors1to200.csv", combinedPortfolios2, ",")
println("Finished Everything")

max_stars_repo_name: skipperafrosenborg/Thesis

max_issues_repo_name: skipperafrosenborg/Thesis

max_forks_repo_name: skipperafrosenborg/Thesis


text:
#=
# Dependency Graph of a Package
In this example we'll plot a dependency graph of a package using
[`PkgDeps.jl`](https://github.com/JuliaEcosystem/PkgDeps.jl) and
and a DAG layout from [`LayeredLayouts.jl`](https://github.com/oxinabox/LayeredLayouts.jl)
=#
using CairoMakie
CairoMakie.activate!(type="png") #hide
CairoMakie.inline!(true) #hide
set_theme!(resolution=(800, 600)) #hide
using GraphMakie
using Graphs
using LayeredLayouts
using PkgDeps
using Makie.GeometryBasics
using Makie.Colors

#=
First we need a small function which goes recursively through the dependencies of a package and
builds a `SimpleDiGraph` object.
=#
function depgraph(root)
    packages = [root]
    connections = Vector{Pair{Int,Int}}()

    for pkg in packages
        pkgidx = findfirst(isequal(pkg), packages)
        deps = direct_dependencies(pkg)

        for dep in keys(deps)
            idx = findfirst(isequal(dep), packages)
            if idx === nothing
                push!(packages, dep)
                idx = lastindex(packages)
            end
            push!(connections, idx => pkgidx)
        end
    end
    g = SimpleDiGraph(length(packages))
    for c in connections
        add_edge!(g, c)
    end
    return (packages, g)
end
nothing #hide

#=
As an example we'll plot the dependency Graph of [`Revise.jl`](https://github.com/timholy/Revise.jl)
because it is one of the most important packages in the Julia ecosystem but does not have a huge
dependency tree.
=#
(packages, g) = depgraph("Revise")
N = length(packages)
xs, ys, paths = solve_positions(Zarate(), g)

## we scale the y coordinates so the plot looks nice in `DataAspect()`
ys .= 0.3 .* ys
foreach(v -> v[2] .= 0.3 .* v[2], values(paths))
nothing #hide

#=
In `GraphMakie` the layout always needs to be function. So we're creating a dummy function...
We will use the [Edge waypoints](@ref) attribute to get the graph with the least crossings.
=#
lay = _ -> Point.(zip(xs,ys))
## create a vector of Point2f per edge
wp = [Point2f.(zip(paths[e]...)) for e in edges(g)]

## manually tweak some of the lable aligns
align = [(:right, :center) for i in 1:N]
align[1]  = (:left, :center)  # Revise
align[3]  = (:right, :top)    # LoweredCodeUtils
align[6]  = (:left, :bottom)  # CodeTracking
align[10] = (:left, :bottom)  # JuliaInterpreter
align[13] = (:left, :bottom)  # Requires

## shift "JuliaInterpreter" node in data space
offset = [Point2f(0,0) for i in 1:N]
offset[10] = Point(-0.1, 0.1)

f, ax, p = graphplot(g; layout=lay,
                     arrow_size=15,
                     edge_color=:gray,
                     nlabels=packages,
                     nlabels_align=align,
                     nlabels_distance=10,
                     nlabels_textsize=15,
                     nlabels_offset=offset,
                     node_size=[9.0 for i in 1:N],
                     edge_width=[3 for i in 1:ne(g)],
                     waypoints=wp,
                     waypoint_radius=0.5)
ax.title = "Dependency Graph of Revise.jl"
xlims!(ax, -0.6, 5.6)
hidedecorations!(ax); hidespines!(ax); ax.aspect = DataAspect()
f #hide

#=
If you run this example using `GLMakie` you can add this code to play
around with the interactive features.
```julia
deregister_interaction!(ax, :rectanglezoom)
register_interaction!(ax, :nodehover, NodeHoverHighlight(p))
register_interaction!(ax, :edgehover, EdgeHoverHighlight(p))
register_interaction!(ax, :edrag, EdgeDrag(p))
register_interaction!(ax, :ndrag, NodeDrag(p))
```
=#


text:
"""
$(SIGNATURES)

Initialize lid-driven cavity
"""
function ib_cavity(
    set::AbstractSetup,
    ps::AbstractPhysicalSpace,
    vs::Union{AbstractVelocitySpace,Nothing},
    gas::AbstractProperty,
    Um = 0.15,
    Vm = 0.0,
    Tm = 1.0,
)

    if set.nSpecies == 1

        prim = [1.0, 0.0, 0.0, 1.0]
        w = prim_conserve(prim, gas.γ)
        primU = [1.0, Um, Vm, Tm]

        p = (y1 = ps.y1, w = w, prim = prim, primU = primU)

        fw = function (args...)
            p = args[end]
            return p.w
        end

        bc = function (x, y, args...)
            p = args[end]

            if y == p.y1
                return p.primU
            else
                return p.prim
            end
        end

        if set.space[1:4] == "2d0f"
            return fw, bc, p
        elseif set.space == "2d1f2v"
            h = maxwellian(vs.u, vs.v, prim)
            p = (p..., h = h)
            ff = function (args...)
                p = args[end]
                return p.h
            end

            return fw, ff, bc, p
        elseif set.space == "2d2f2v"
            h = maxwellian(vs.u, vs.v, prim)
            b = h .* gas.K / 2.0 / prim[end]
            p = (p..., h = h, b = b)
            ff = function (args...)
                p = args[end]
                return p.h, p.b
            end

            return fw, ff, bc, p
        end

    end

    return nothing

end


text:
# Raby, C. R. and Alexis, D. M. and Dickinson, A. and Clayton, N. S., Planning
# for the future by western scrub-jays, Nature 445:919-921, 2007,
# http://dx.doi.org/10.1038/nature05575

function _summarize(::Experiment{:Raby07_planningforbreakfast}, data)
    combine(groupby(data, [:tray, :action]),
				df -> DataFrame(μ = mean(df.items), sem = sem(df.items)))
end
function statistical_tests(exp::Experiment{:Raby07_planningforbreakfast}, data)
    test = is_significantly_bigger(data, conditionkey = :tray,
                                   condbigger = "no-breakfast",
                                   condsmaller = "breakfast", valuekey = :items)
    TestCollection(Test("More pine nuts in no-breakfast compartment",
                     test,
                     exp.tests))
end
function run!(exp::Experiment{:Raby07_planningforbreakfast}, models; N = 8)
    results = DataFrame(tray = String[], action = String[], id = Int[], items = Float64[])
	for id in 1:N
		m = models[id]
		trays = [Tray(i) for i in 1:2]
		for day in 1:7
			add!(m, MaintenanceDiet) # start at 9am on first day
			wait!(m, 8u"hr")
			remove!(m, MaintenanceDiet)
			add!(m, trays)
			wait!(m, 1.5u"hr")
			if day == 7
				add!(m, Pinenut, 30)
			else
				add!(m, PowderedPinenut, 30)
			end
			wait!(m, 30u"minute")
			if day == 7
				push!(results, ["breakfast", "cache", id, countcache(trays[1])])
				push!(results, ["no-breakfast", "cache", id, countcache(trays[2])])
				break
			end
			remove!(m, Any)
			wait!(m, 12u"hr")
			if (id + day) % 2 == 0
				add!(m, PowderedPinenut, 30)
				add!(m, trays[1])
			else
				add!(m, trays[2])
			end
			wait!(m, 2u"hr")
			remove!(m, Any)
		end
	end
	results
end

function _summarize(::Experiment{:Raby07_breakfastchoice}, results)
	combine(groupby(results, [:tray, :action, :foodtype]),
				df -> DataFrame(μ = mean(df.items), sem = sem(df.items)))
end
function statistical_tests(exp::Experiment{:Raby07_breakfastchoice}, data)
    id = "More different food cache"
    test = splitplotanova(data,
                          exp.tests[id].locals,
                          betweenfactors = [],
                          withinfactors = [:tray, :foodtype],
                          id = :id, data = :items)
    TestCollection(Test(id, test, exp.tests))
end
function run!(exp::Experiment{:Raby07_breakfastchoice}, models; N = 9)
    results = DataFrame(tray = String[], action = String[], foodtype = String[],
						id = Int[], items = Float64[])
	for id in 1:N
		m = models[id]
		trays = [Tray(i) for i in 1:2]
		for day in 1:7
			add!(m, MaintenanceDiet) # start at 9am on first day
			wait!(m, 8u"hr")
			remove!(m, MaintenanceDiet)
			add!(m, trays)
			wait!(m, 1.5u"hr")
			if day == 7
				add!(m, Peanut, 30)
				add!(m, Kibble, 30)
			else
				add!(m, PowderedPeanut, 30)
				add!(m, PowderedKibble, 30)
			end
			wait!(m, 30u"minute")
			if day == 7
				push!(results, ["same", "cache", "peanut", id, countcache(trays[1], Peanut)])
				push!(results, ["same", "cache", "kibble", id, countcache(trays[2], Kibble)])
				push!(results, ["other", "cache", "peanut", id, countcache(trays[2], Peanut)])
				push!(results, ["other", "cache", "kibble", id, countcache(trays[1], Kibble)])
				break
			end
			remove!(m, Any)
			wait!(m, 12u"hr")
			if (id + day) % 2 == 0
				add!(m, PowderedPeanut, 30)
				add!(m, trays[1])
			else
				add!(m, PowderedKibble, 30)
				add!(m, trays[2])
			end
			wait!(m, 2u"hr")
			remove!(m, Any)
		end
	end
	results
end


text:
module Initial

include("./Node.jl")
using .Node

const collectionName = "animation"  # Collection name who will store data in DB

const tick = 0.1        # Single time tick (Seconds)
const totalTicks = 10   # Total frames to calculate
# framesPerSecond = 1

const onCollusion = "bounce"

# All the nodes
nodes = [
    node(1E20, 5, [0, 0], [0.1, 0]),
    node(1E20, 5, [1, 1], [0.1, 0]),
    node(1E20, 5, [3, 4], [0.1, 0]),
    node(1E20, 5, [0, 7], [0.1, 0])
]

export nodes, tick, totalTicks, collectionName, onCollusion
end

max_stars_repo_name: PseudoDevelopers/Physics-Simulations

max_issues_repo_name: PseudoDevelopers/Physics-Simulations

max_forks_repo_name: PseudoDevelopers/Physics-Simulations


text:
struct Nmix <: UnmarkedModel
  data::UmData
  opt::UnmarkedOpt
  submodels::NamedTuple
  K::Int
end

"""
    nmix(λ_formula::FormulaTerm, p_formula::FormulaTerm, data::UmData; K::Int)

Fit single-season N-mixture models. Covariates on abundance and detection
are specified with `λ_formula` and `p_formula`, respectively. The `UmData`
object should contain `site_covs` and `obs_covs` `DataFrames` containing
columns with names matching the formulas. The likelihood is marginalized over
possible latent abundance values `0:K` for each site. `K` should be set high
enough that the likelihood of true abundance `K` for any site is ≈ 0 (defaults
to maximum observed count + 20). Larger `K` will increase runtime.
"""
function nmix(λ_formula::FormulaTerm, p_formula::FormulaTerm,
              data::UmData; K::Union{Int,Nothing}=nothing)

  abun = UmDesign(:Abundance, λ_formula, LogLink(), data.site_covs)
  det = UmDesign(:Detection, p_formula, LogitLink(), data.obs_covs)
  add_idx!([abun, det])
  np = get_np([abun, det])

  y = data.y
  N, J = size(y)
  K = check_K(K, y)
  Kmin = get_Kmin(y)

  function loglik(β::Array)

    λ = transform(abun, β)
    p = transform(det, β)

    ll = zeros(eltype(β),N)

    idx = 0
    for i = 1:N

      if all(ismissing.(y[i,:]))
        idx += J
        continue
      end

      fg = zeros(eltype(β),K+1)
      for k = Kmin[i]:K
        fk = pdf(Poisson(λ[i]), k)
        gk = 1.0
        for j = 1:J
          idx += 1
          if ismissing(y[i,j])
            continue
          end
          gk *= pdf(Binomial(k, p[idx]), y[i,j])
        end
        fg[k+1] = fk * gk

        #Return p index to start for site if necessary
        if k != K idx -= J end
      end
      ll[i] = log(sum(fg))
    end

    return -sum(ll)
  end

  opt = optimize_loglik(loglik, np)
  Nmix(data, opt, (abun=UnmarkedSubmodel(abun,opt),
            det=UnmarkedSubmodel(det,opt)), K)
end

#Fit multiple models at once
"""
    nmix(λ_formula::Union{Array,FormulaTerm},
         p_formula::Union{Array,FormulaTerm}, data::UmData, K::Int)

Provide multiple formulas for λ and/or p and fit models for all formula
combinations. To fit all model subsets, wrap formula with `allsubs()`.
"""

function nmix(λ_form::Union{Array,FormulaTerm},
              p_form::Union{Array,FormulaTerm},
              data::UmData, K::Union{Int,Nothing}=nothing)
  cf = combine_formulas(λ_form, p_form)
  out = Array{UnmarkedModel}(undef, length(cf))
  p = Progress(length(cf), desc=string("Fitting ", length(cf), " models "))
  Threads.@threads for i in 1:length(cf)
    out[i] = nmix(cf[i][1],cf[i][2], data, K)
    next!(p)
  end
  return UnmarkedModels(out)
end

# Simulations------------------------------------------------------------------

"Simulate new dataset from a fitted model"
function simulate(fit::Nmix)

  ydims = collect(size(fit.data.y))
  λ = predict(abundance(fit))
  p = predict(detection(fit))

  y = sim_y(Nmix, λ, p)
  rep_missing!(y, fit.data.y)

  UmData(y, fit.data.site_covs, fit.data.obs_covs)
end

"Simulate new dataset from provided formulas"
function simulate(::Type{Nmix}, λ_form::FormulaTerm, p_form::FormulaTerm,
                  ydims::Tuple{Int,Int}, coef::Array)

  sc = gen_covs(λ_form, ydims[1])
  oc = gen_covs(p_form, prod(ydims))

  abun = UmDesign(:Abun, λ_form, LogLink(), sc)
  det = UmDesign(:Det, p_form, LogitLink(), oc)
  add_idx!([abun, det])
  np = get_np([abun, det])

  if np != length(coef) error(string("Coef array must be length ",np)) end

  λ = transform(abun, coef)
  p = transform(det, coef)

  y = sim_y(Nmix, λ, p)

  UmData(y, sc, oc)
end

function sim_y(::Type{Nmix}, λ::Array, p::Array)

  N = length(λ)
  J = Integer(length(p) / N)
  z = map(x -> rand(Poisson(x))[1], λ)
  y = Array{Union{Int,Missing}}(undef, N, J)
  idx = 0
  for i = 1:N
    for j = 1:J
      idx += 1
      y[i,j] = z[i] == 0 ? 0 : rand(Binomial(z[i], p[idx]))[1]
    end
  end

  return y
end

# Update ----------------------------------------------------------------------

function update(fit::Nmix, data::UmData)
  nmix(abundance(fit).formula, detection(fit).formula, data, K=fit.K)
end

# Goodness-of-fit -------------------------------------------------------------

function fitted(fit::Nmix)
  y = fit.data.y
  N,J = size(fit.data.y)
  λ = repeat(predict(abundance(fit)), inner=J)
  return λ .* predict(detection(fit))
end

function residuals(fit::Nmix)
  ft = fitted(fit)
  yl = reshape(transpose(fit.data.y), length(fit.data.y))
  return ft .- yl
end

function chisq(fit::Nmix)
  return sum(skipmissing(residuals(fit) .^ 2 ./ fitted(fit)))
end

"""
    gof(fit::Nmix, nsims::Int=30)

Compute an N-mixture goodness-of-fit test using Pearson residuals.
The distribution of the test statistic is generated using parametric
bootstrapping with `nsims` simulations.
"""
function gof(fit::Nmix, nsims::Int=30)
  test_name = "N-mixture Goodness-of-fit"
  t0 = chisq(fit)
  tstar = parboot(fit, nsims, chisq)
  return UmGOF(test_name,t0, tstar)
end

max_stars_repo_name: kenkellner/julia-hme

max_issues_repo_name: kenkellner/julia-hme

max_forks_repo_name: kenkellner/julia-hme


text:
export ZabbixTemplate, get_templates

struct ZabbixTemplate
  id::Int
  name::String
end
ZabbixTemplate(dict) =
  ZabbixTemplate(
    parse(Int, dict["templateid"]),
    dict["name"]
  )

showres(t::ZabbixTemplate; offset=0) =
  println(repeat(" ", offset), "ZabbixTemplate: ", t.id, " => ", t.name)

function get_templates(z)
  _templates = @checkerr make_request(z, "template.get")
  templates = ZabbixTemplate[]
  for t in _templates
    push!(templates, ZabbixTemplate(t))
  end
  templates
end