text:
## Functions for scanning for motifs that may be disrupted by variants

### build scanning sequence for a given motif and mutation described by indexes to the left and to the right exclusive of mutation
### this allows scanning of all positions in which motif intersects with at least one bp of mutated sequence
function motifscanseq(seq, leftind, rightind, motif)
    n = size(motif.pwm, 2)
    start = max(last(leftind) - n + 1 + 1, 1)
    stop  = min(first(rightind) + n - 1 - 1, length(seq))
    seq[start:stop], start
end

### find maximal forward or reverse match of a sequence
function scanmax(seq, motif)
    fs, rs = scanmotif(seq, motif.pbg)
    fm, fi = findmax(fs)
    rm, ri = findmax(rs)
    n = size(motif.pbg, 2)
    if fm > rm
        maxscore = fm
        start = fi
        stop = fi + n - 1
        strand = "+"

    else
        maxscore = rm
        start = ri
        stop  = ri + n - 1
        strand = "-"
    end
    maxscore, start, stop, strand
end

### scan all mutations with all motifs
function motifscannall(refseq, altseqs, distinctmutations, motifs)

    dfs = DataFrame[]
    for (as, row) in zip(altseqs, eachrow(distinctmutations))
        df = scanmots(refseq, as, row.LeftInd, row.RefRightInd, row.AltRightInd, motifs)
        df[!, :ID] .= row.ID
        push!(dfs, df)
    end
    reduce(vcat, dfs)
end

### scan a given mutation with all motifs
function scanmots(refseq, altseq, leftind, refrightind, altrightind, motifs)
    df = DataFrame(MotifName=String[], MotifID=String[], MotSeqStart=Int[], RefMaxScore=Float64[], RefStart=Int[], RefStop=Int[], RefStrand=String[], AltMaxScore=Float64[], AltStart=Int[], AltStop=Int[], AltStrand=String[],
        RefPrMax=Float64[], AltPrMax=Float64[], LR_RefAlt=Float64[], PR_RefAlt=Float64[])

    for m in motifs
        maxscore = sum(maximum(m.pbg, dims=1))
        refmseq, refstart = motifscanseq(refseq, leftind, refrightind, m)
        altmseq, altstart = motifscanseq(altseq, leftind, altrightind, m)
        refres = scanmax(refmseq, m)
        altres = scanmax(altmseq, m)

        ref_prmax = first(refres)/maxscore
        alt_prmax = first(altres)/maxscore
        lr_refalt = first(refres) - first(altres)
        pr_refalt = ref_prmax - alt_prmax

        push!(df, (m.name, m.id, refstart, refres..., altres..., ref_prmax, alt_prmax, lr_refalt, pr_refalt))
    end
    df.RefSeqStart = df.MotSeqStart + df.RefStart .- 1
    df.RefSeqStop  = df.MotSeqStart + df.RefStop .- 1

    df
end

### load motifs from JASPAR transfac file
function loadtransfac(file="JASPAR2020_CORE_vertebrates_non-redundant_pfms_transfac.txt.gz", w = 1e-2, background=MotifScanner.background_human_zero_markov(), projdir=getprojectdir())

    filepath = joinpath(projdir, "data", file)
    io = open(filepath) |> GzipDecompressorStream
    mode = 1

    motifname = ""
    motifid = ""
    currentmotif = Vector{Vector{Float64}}()
    pwms = Vector{Matrix{Float64}}()
    nsites = Float64[]

    motifdata = DataFrame(MotifName=String[], MotifID=String[], Property=String[], Value=String[])

    for line in eachline(io)

        fields = split(line)

        if mode == 1
            if fields[1] == "AC"
                mode = 2
                currentmotif = Vector{Vector{Float64}}()
            end
        elseif mode == 2
            if fields[1] == "DE"
                motifid = fields[2]
                motifname = fields[3]
                if startswith(motifname, motifid) ## JASPAR 2022 release appears to have altered the format of JASPAR 2020, remove the leading motifid
                    motifname = replace(motifname, string(motifid, ".") => "")
                end
            elseif fields[1] == "PO"
                @assert fields[2] == "A"
                @assert fields[3] == "C"
                @assert fields[4] == "G"
                @assert fields[5] == "T"
                mode = 3
            elseif fields[1] == "CC"
                ind = findfirst(':', line)

                cfields = [line[4:(ind-1)], line[(ind+1):end]]

                if length(cfields) != 2
                    @show line
                    @show cfields
                    error("")
                end
                push!(motifdata, (motifname, motifid, cfields...))
            elseif fields[1] == "//"
                mode = 1
            end
        elseif mode == 3
            if fields[1] == "XX"
                mode = 2
                pwm = reduce(hcat, currentmotif)
                v = sum(pwm, dims=1)
                push!(nsites, mean(v))
                pwm ./= sum(pwm, dims=1)
                pwm .+= w
                pwm ./= sum(pwm, dims=1)
                push!(pwms, pwm)
            else
                v = parse.(Float64, fields[2:end])
                push!(currentmotif, v)
            end
        end
    end
    close(io)
    motifmeta = unstack(motifdata, :Property, :Value)

    motifmeta[!, :NumSites] = nsites
    motifmeta[!, :Length] = size.(pwms, 2)

    @assert iszero(mapreduce(c -> sum(ismissing, c), +, eachcol(motifmeta)))
    dropmissing!(motifmeta)

    pbgs = [log2.(pwm) .- log2.(background) for pwm in pwms]

    motifdata = [(name=motifname, id=id, pwm=pwm, pbg=pbg) for (motifname, id, pwm, pbg) in zip(motifmeta.MotifName, motifmeta.MotifID, pwms, pbgs)]
    motifmeta.MotifFam = motiffamily.(motifmeta.MotifName);
    motifmeta, motifdata

end

### motif family from motif tf name, crude but does the job
function motiffamily(motifname)
    mf = uppercase.(motifname)
    mf = replace(mf, r"\([A-Za-z0-9.]*\)$" => "")

    mf = replace(mf, r"[0-9]*$" => "")
    mf = replace(mf, r"-$" => "")
    mf = replace(mf, r"FOX[A-Z]" => "FOX")
end

function moteq(mA, mB)
    (mA.name != mB.name) && return false
    (mA.id != mB.id) && return false

    (maximum(mA.pwm - mB.pwm) > 1e-5) && return false
    return true
end

overlap(grp, starts, stops) = overlappinglocations(fill(grp, length(starts)), [s:e for (s, e) in zip(starts, stops)])
function filtermatches(motscan; minscore=0.6, maxscore=Inf)

    motmatch = @chain motscan begin
        @subset(minscore .<= :RefPrMax .< maxscore)                                   # filter by norm motif score
        groupby(:ID)                                                                  # group by mutation id
        transform(:PR_RefAlt => (x -> sortperm(sortperm(x, rev=true))) => :Rank)      # Rank all motif matches (grouped by ids)
        sort([:MotifFamily, :ID, :Rank])                                              # sort
        transform([:MotifFamily, :RefSeqStart, :RefSeqStop] => overlap => :Group)     # Assign grouping label for overlapping motifs
        sort([:ID, :Rank])                                                            # sort
    end
    motmatch
end

function motif_rank1(motmatch)
    motsel = @chain motmatch begin
        @subset(:Rank .== 1)                                                          # rank of motif for each mutation
        groupby(:Group)                                                               # multiple rank 1 motifs overlap, group by those that overlap
        combine(df -> df[argmax(df.PR_RefAlt), :])                                    # select the motif with largest disruption score
        sort(:ID)
    end
    motsel
end

function motifmatchtable(motmatch_tier1, motmatch_tier2, region)

    motmatch_tier1[!, :Tier] .= 1
    motmatch_tier2[!, :Tier] .= 2

    motmatchall = [ motmatch_tier1 ; motmatch_tier2]

    motmatchall.MutSeqStart = motmatchall.MotSeqStart .+ motmatchall.AltStart .- 1
    motmatchall.MutSeqStop = motmatchall.MotSeqStart .+ motmatchall.AltStop .- 1
    motmatchall.RefSeqMotifStart = leftposition(region) .+ motmatchall.RefSeqStart .- 1
    motmatchall.RefSeqMotifStop  = leftposition(region) .+ motmatchall.RefSeqStop .- 1
    motmatchall.MutSeqMotifStart = leftposition(region) .+ motmatchall.MutSeqStart .- 1
    motmatchall.MutSeqMotifStop  = leftposition(region) .+ motmatchall.MutSeqStop .- 1
    motmatchall[!, :chrom] .= "chr10"
    motmatchall[!, :Mutation] = [string.("g.", f[1], ifelse(f[3] == "", string("del", f[2]), string(f[2], ">", f[3]))) for f in split.(motmatchall.ID, "_")]
    rename!(motmatchall, :RefStrand => :RefSeqStrand, :RefMaxScore => :RefSeqMaxScore, :AltStrand => :MutSeqStrand, :AltMaxScore => :MutSeqMaxScore,
                :RefPrMax => :RefSeqNormMaxScore, :AltPrMax => :AltSeqNormMaxScore, :LR_RefAlt => :LR_RefMut, :PR_RefAlt => :Norm_RefMut)
    mma = motmatchall[!, [:Mutation, :Tier, :MotifFamily, :MotifName, :MotifID, :chrom,
                        :RefSeqMotifStart, :RefSeqMotifStop, :RefSeqStrand, :RefSeqMaxScore,
                        :MutSeqMotifStart, :MutSeqMotifStop, :MutSeqStrand, :MutSeqMaxScore,
                        :RefSeqNormMaxScore, :AltSeqNormMaxScore, :LR_RefMut, :Norm_RefMut, :Rank, :Group]]

    mma
end


text:
# This contains the main slacker program

const INIT_ERROR_MSG = "Config file not found. Please create a settings file using `addConfig(name)`."
const KEY_NOT_FOUND_ERROR_MSG = "No configuration found for this name. Please add using `addConfig(name)`"

# A structure that contains a configuration for an incoming webhook on slack
mutable struct SlackConfig
    webhook_url:: String
    user:: String
    channel:: String
    icon_emoji:: String
    function SlackConfig()
        SlackConfig("http://www.example.org", "JuliaBot", "#general", ":ghost:")
    end

    function SlackConfig(dic)
        SlackConfig(dic["webhook_url"], dic["user"], dic["channel"], dic["icon_emoji"])
    end

    """
    ...
    # Arguments
    - `webhook_url::String`: the incoming webhook_url
    - `user::String`: the name of the user Slacker should have in slack.
    - `channel::String`: the channel Slacker should post to. Allows channel with # and direc messages using @.
    - `icon_emoji::String` : the icon the user should have when posting. By adding a custom emoji e.g. :juliabot: you can easily give julia its own avatar.
    ...
    """
    function SlackConfig(webhook_url, user, channel, icon_emoji)
        new(webhook_url, user, channel, icon_emoji)
    end

end

function getConfigFile()
    joinpath(homedir(), ".slacker_config.json")
end

"""
Load the registry and return all registered servers as a Dict.
"""
function readConfigFile()
    fn = getConfigFile()
    isfile(fn) || error(INIT_ERROR_MSG)
    dic = Dict{String, SlackConfig}()

    # read file
    str = open(fn) do file
        read(file, String)
    end

    # parse json to Dict
    res = JSON.parse(str)
    for entry in res
        val = SlackConfig(entry[2])
        dic[entry[1]] = val
    end
    dic
end

"""
Add a Slack server configuration to the registry. Shorthand to create a
 configuration that posts to the random channel (uses Juliabot and :ghost:
 as defaults).

# Examples
```julia-repl
julia> addConfig("server.url", "config_name")
```
"""
function addConfig(webhook_url::String, name::String)
    cfg = SlackConfig(webhook_url, "Juliabot", "#random",":ghost:")
    addConfig(cfg, name)
end

"""
Add a Slack server configuration to the registry. Multiple servers can be
defined by adjusting the `name` parameter.

# Examples
```julia-repl
julia> addConfig(SlackConfig("server.url", "username", "#channel", ":smile:"))
```
"""
function addConfig(config::SlackConfig, name::String = "default")
    fn = getConfigFile()
    dic = Dict(name => config)
    js = JSON.json(dic)

    if isfile(fn)
        @debug "File exists."
        configs = readConfigFile()
        if haskey(configs,name)
            @warn "Overwriting exisiting configuration."
            configs[name] = config
            open(fn, write = true) do file
                write(file, JSON.json(configs))
            end
        else
            configs[name] = config
            @debug "Append Key"
            open(fn, write = true) do file
                write(file, JSON.json(configs))
            end
        end
    else
        @debug "New File."
        open(fn, write = true) do file
            write(file, js)
        end
    end
end

"""
Deletes the whole registry from the computer.
"""
function removeConfigFile()
    fn = getConfigFile()
    if isfile(fn)
        rm(fn)
    end
end

"""
Load the default configuration from the registry. Other configurations may be loaded
by specifying the `name` parameter.
"""
function loadConfig(name="default")
    settings = readConfigFile()
    if haskey(settings, name)
        return settings[name]
    else
        error(KEY_NOT_FOUND_ERROR_MSG)
    end

end

"""
Sends a message `text`to Server defined in the configuration `cfg`.
"""
function sendMessage(text, config::SlackConfig)
    a = Dict(
        "channel" => config.channel,
        "username" => config.user,
        "text" => text,
        "icon_emoji" => config.icon_emoji)

    message = "payload=" * JSON.json(a)
    webhook = config.webhook_url
    try
        resp = HTTP.request("POST",
            webhook,
            ["Content-Type" => "application/x-www-form-urlencoded"], message;
            verbose = 0)
    catch e
        @warn "Could not send message $e"
    end
end

"""
Sends a message `text`to the configured default slack server or the one specified in `cfg`.
"""
function sendMessage(text, cfg::String = "default")
    config = loadConfig(cfg)
    sendMessage(text, config)
end

"""
Tests whether a configuration exists
"""
function hasConfig(name = "default")
    fn = getConfigFile()
    if isfile(fn)
        dic = readConfigFile()
        if haskey(dic, name)
            return true
        end
    end
    false
end

"""
removes an individual configuration
call with parameter "default" to delete the default configuration
"""
function removeConfig(name)
    #helps debugging
    #name = "default"
    fn = getConfigFile()
    dic = readConfigFile()
    if haskey(dic, name)
        cfg = pop!(dic, name)
    end
    # helps for debugging
    #dic[name] = cfg
    open(fn, write = true) do file
        write(file, JSON.json(dic))
    end
    "Configuration $name removed."
end

ext: jl


text:
# callbacks.jl
# Define callbacks for periodic dilution
# Separate definitions for different models
# July 2020, Nadanai Laohakunakorn (nadanai.laohakunakorn@ed.ac.uk)

###############################################
######### 1. Resource dependent model #########
###############################################

function diluteSR!(integrator)
    if integrator.t<SWITCHTIMES[1]*60
        # Stage 1
        INDEX_REFRESH = [idx_R,idx_dG,idx_pT,idx_dT]
        CONC_REFRESH = [R0,dG0,pT0,dT0]
        for j in 1:NSPECIES
            integrator.u[j] = integrator.u[j]*(1-DIL_FRAC)
        end
        for j in 1:size(INDEX_REFRESH)[1]
            integrator.u[INDEX_REFRESH[j]] = integrator.u[INDEX_REFRESH[j]] + DIL_FRAC*CONC_REFRESH[j]
        end

    elseif SWITCHTIMES[1]*60<=integrator.t<SWITCHTIMES[2]*60
        # Stage 2
        INDEX_REFRESH = [idx_R,idx_dG,idx_dT]
        CONC_REFRESH = [R0,dG0,dT0]
        for j in 1:NSPECIES
            integrator.u[j] = integrator.u[j]*(1-DIL_FRAC)
        end
        for j in 1:size(INDEX_REFRESH)[1]
            integrator.u[INDEX_REFRESH[j]] = integrator.u[INDEX_REFRESH[j]] + DIL_FRAC*CONC_REFRESH[j]
        end

    elseif SWITCHTIMES[2]*60<=integrator.t
        # Stage 3
        INDEX_REFRESH = [idx_R,idx_dG]
        CONC_REFRESH = [R0,dG0]
        for j in 1:NSPECIES
            integrator.u[j] = integrator.u[j]*(1-DIL_FRAC)
        end
        for j in 1:size(INDEX_REFRESH)[1]
            integrator.u[INDEX_REFRESH[j]] = integrator.u[INDEX_REFRESH[j]] + DIL_FRAC*CONC_REFRESH[j]
        end

    end
end

function dilutePC!(integrator)
    if integrator.t<SWITCHTIMES[1]*60
        # Stage 1
        INDEX_REFRESH = [idx_R,idx_dG,idx_pT]
        CONC_REFRESH = [R0,dG0,pT0]
        for j in 1:NSPECIES
            integrator.u[j] = integrator.u[j]*(1-DIL_FRAC)
        end
        for j in 1:size(INDEX_REFRESH)[1]
            integrator.u[INDEX_REFRESH[j]] = integrator.u[INDEX_REFRESH[j]] + DIL_FRAC*CONC_REFRESH[j]
        end

    elseif SWITCHTIMES[1]*60<=integrator.t<SWITCHTIMES[2]*60
        # Stage 2
        INDEX_REFRESH = [idx_R,idx_dG,idx_pT]
        CONC_REFRESH = [R0,dG0,pT0]
        for j in 1:NSPECIES
            integrator.u[j] = integrator.u[j]*(1-DIL_FRAC)
        end
        for j in 1:size(INDEX_REFRESH)[1]
            integrator.u[INDEX_REFRESH[j]] = integrator.u[INDEX_REFRESH[j]] + DIL_FRAC*CONC_REFRESH[j]
        end

    elseif SWITCHTIMES[2]*60<=integrator.t
        # Stage 3
        INDEX_REFRESH = [idx_R,idx_dG,idx_pT]
        CONC_REFRESH = [R0,dG0,pT0]
        for j in 1:NSPECIES
            integrator.u[j] = integrator.u[j]*(1-DIL_FRAC)
        end
        for j in 1:size(INDEX_REFRESH)[1]
            integrator.u[INDEX_REFRESH[j]] = integrator.u[INDEX_REFRESH[j]] + DIL_FRAC*CONC_REFRESH[j]
        end

    end
end

function diluteNC!(integrator)
    if integrator.t<SWITCHTIMES[1]*60
        # Stage 1
        INDEX_REFRESH = [idx_R,idx_dG,idx_pT]
        CONC_REFRESH = [R0,dG0,pT0]
        for j in 1:NSPECIES
            integrator.u[j] = integrator.u[j]*(1-DIL_FRAC)
        end
        for j in 1:size(INDEX_REFRESH)[1]
            integrator.u[INDEX_REFRESH[j]] = integrator.u[INDEX_REFRESH[j]] + DIL_FRAC*CONC_REFRESH[j]
        end

    elseif SWITCHTIMES[1]*60<=integrator.t<SWITCHTIMES[2]*60
        # Stage 2
        INDEX_REFRESH = [idx_R,idx_dG]
        CONC_REFRESH = [R0,dG0]
        for j in 1:NSPECIES
            integrator.u[j] = integrator.u[j]*(1-DIL_FRAC)
        end
        for j in 1:size(INDEX_REFRESH)[1]
            integrator.u[INDEX_REFRESH[j]] = integrator.u[INDEX_REFRESH[j]] + DIL_FRAC*CONC_REFRESH[j]
        end

    elseif SWITCHTIMES[2]*60<=integrator.t
        # Stage 3
        INDEX_REFRESH = [idx_R,idx_dG]
        CONC_REFRESH = [R0,dG0]
        for j in 1:NSPECIES
            integrator.u[j] = integrator.u[j]*(1-DIL_FRAC)
        end
        for j in 1:size(INDEX_REFRESH)[1]
            integrator.u[INDEX_REFRESH[j]] = integrator.u[INDEX_REFRESH[j]] + DIL_FRAC*CONC_REFRESH[j]
        end

    end
end

###############################################
######## 2. Resource independent model ########
###############################################

function diluteSR_RI!(integrator)
    if integrator.t<SWITCHTIMES[1]*60
        # Stage 1
        INDEX_REFRESH = [idx_TX,idx_TL,idx_dG,idx_pT,idx_dT]
        CONC_REFRESH = [TX0,TL0,dG0,pT0,dT0]
        for j in 1:NSPECIES
            integrator.u[j] = integrator.u[j]*(1-DIL_FRAC)
        end
        for j in 1:size(INDEX_REFRESH)[1]
            integrator.u[INDEX_REFRESH[j]] = integrator.u[INDEX_REFRESH[j]] + DIL_FRAC*CONC_REFRESH[j]
        end

    elseif SWITCHTIMES[1]*60<=integrator.t<SWITCHTIMES[2]*60
        # Stage 2
        INDEX_REFRESH = [idx_TX,idx_TL,idx_dG,idx_dT]
        CONC_REFRESH = [TX0,TL0,dG0,dT0]
        for j in 1:NSPECIES
            integrator.u[j] = integrator.u[j]*(1-DIL_FRAC)
        end
        for j in 1:size(INDEX_REFRESH)[1]
            integrator.u[INDEX_REFRESH[j]] = integrator.u[INDEX_REFRESH[j]] + DIL_FRAC*CONC_REFRESH[j]
        end

    elseif SWITCHTIMES[2]*60<=integrator.t
        # Stage 3
        INDEX_REFRESH = [idx_TX,idx_TL,idx_dG]
        CONC_REFRESH = [TX0,TL0,dG0]
        for j in 1:NSPECIES
            integrator.u[j] = integrator.u[j]*(1-DIL_FRAC)
        end
        for j in 1:size(INDEX_REFRESH)[1]
            integrator.u[INDEX_REFRESH[j]] = integrator.u[INDEX_REFRESH[j]] + DIL_FRAC*CONC_REFRESH[j]
        end

    end
end

function dilutePC_RI!(integrator)
    if integrator.t<SWITCHTIMES[1]*60
        # Stage 1
        INDEX_REFRESH = [idx_TX,idx_TL,idx_dG,idx_pT]
        CONC_REFRESH = [TX0,TL0,dG0,pT0]
        for j in 1:NSPECIES
            integrator.u[j] = integrator.u[j]*(1-DIL_FRAC)
        end
        for j in 1:size(INDEX_REFRESH)[1]
            integrator.u[INDEX_REFRESH[j]] = integrator.u[INDEX_REFRESH[j]] + DIL_FRAC*CONC_REFRESH[j]
        end

    elseif SWITCHTIMES[1]*60<=integrator.t<SWITCHTIMES[2]*60
        # Stage 2
        INDEX_REFRESH = [idx_TX,idx_TL,idx_dG,idx_pT]
        CONC_REFRESH = [TX0,TL0,dG0,pT0]
        for j in 1:NSPECIES
            integrator.u[j] = integrator.u[j]*(1-DIL_FRAC)
        end
        for j in 1:size(INDEX_REFRESH)[1]
            integrator.u[INDEX_REFRESH[j]] = integrator.u[INDEX_REFRESH[j]] + DIL_FRAC*CONC_REFRESH[j]
        end

    elseif SWITCHTIMES[2]*60<=integrator.t
        # Stage 3
        INDEX_REFRESH = [idx_TX,idx_TL,idx_dG,idx_pT]
        CONC_REFRESH = [TX0,TL0,dG0,pT0]
        for j in 1:NSPECIES
            integrator.u[j] = integrator.u[j]*(1-DIL_FRAC)
        end
        for j in 1:size(INDEX_REFRESH)[1]
            integrator.u[INDEX_REFRESH[j]] = integrator.u[INDEX_REFRESH[j]] + DIL_FRAC*CONC_REFRESH[j]
        end

    end
end

function diluteNC_RI!(integrator)
    if integrator.t<SWITCHTIMES[1]*60
        # Stage 1
        INDEX_REFRESH = [idx_TX,idx_TL,idx_dG,idx_pT]
        CONC_REFRESH = [TX0,TL0,dG0,pT0]
        for j in 1:NSPECIES
            integrator.u[j] = integrator.u[j]*(1-DIL_FRAC)
        end
        for j in 1:size(INDEX_REFRESH)[1]
            integrator.u[INDEX_REFRESH[j]] = integrator.u[INDEX_REFRESH[j]] + DIL_FRAC*CONC_REFRESH[j]
        end

    elseif SWITCHTIMES[1]*60<=integrator.t<SWITCHTIMES[2]*60
        # Stage 2
        INDEX_REFRESH = [idx_TX,idx_TL,idx_dG]
        CONC_REFRESH = [TX0,TL0,dG0]
        for j in 1:NSPECIES
            integrator.u[j] = integrator.u[j]*(1-DIL_FRAC)
        end
        for j in 1:size(INDEX_REFRESH)[1]
            integrator.u[INDEX_REFRESH[j]] = integrator.u[INDEX_REFRESH[j]] + DIL_FRAC*CONC_REFRESH[j]
        end

    elseif SWITCHTIMES[2]*60<=integrator.t
        # Stage 3
        INDEX_REFRESH = [idx_TX,idx_TL,idx_dG]
        CONC_REFRESH = [TX0,TL0,dG0]
        for j in 1:NSPECIES
            integrator.u[j] = integrator.u[j]*(1-DIL_FRAC)
        end
        for j in 1:size(INDEX_REFRESH)[1]
            integrator.u[INDEX_REFRESH[j]] = integrator.u[INDEX_REFRESH[j]] + DIL_FRAC*CONC_REFRESH[j]
        end

    end
end

###############################################
## 3. Resource independent model  no loading ##
###############################################

function diluteSR_RInl!(integrator)
    if integrator.t<SWITCHTIMES[1]*60
        # Stage 1
        INDEX_REFRESH = [idx_TX,idx_TL,idx_dG,idx_pT,idx_dT]
        CONC_REFRESH = [TX0,TL0,dG0,pT0,dT0]
        for j in 1:NSPECIES
            integrator.u[j] = integrator.u[j]*(1-DIL_FRAC)
        end
        for j in 1:size(INDEX_REFRESH)[1]
            integrator.u[INDEX_REFRESH[j]] = integrator.u[INDEX_REFRESH[j]] + DIL_FRAC*CONC_REFRESH[j]
        end

    elseif SWITCHTIMES[1]*60<=integrator.t<SWITCHTIMES[2]*60
        # Stage 2
        INDEX_REFRESH = [idx_TX,idx_TL,idx_dG,idx_dT]
        CONC_REFRESH = [TX0,TL0,dG0,dT0]
        for j in 1:NSPECIES
            integrator.u[j] = integrator.u[j]*(1-DIL_FRAC)
        end
        for j in 1:size(INDEX_REFRESH)[1]
            integrator.u[INDEX_REFRESH[j]] = integrator.u[INDEX_REFRESH[j]] + DIL_FRAC*CONC_REFRESH[j]
        end

    elseif SWITCHTIMES[2]*60<=integrator.t
        # Stage 3
        INDEX_REFRESH = [idx_TX,idx_TL,idx_dG]
        CONC_REFRESH = [TX0,TL0,dG0]
        for j in 1:NSPECIES
            integrator.u[j] = integrator.u[j]*(1-DIL_FRAC)
        end
        for j in 1:size(INDEX_REFRESH)[1]
            integrator.u[INDEX_REFRESH[j]] = integrator.u[INDEX_REFRESH[j]] + DIL_FRAC*CONC_REFRESH[j]
        end

    end
end

function dilutePC_RInl!(integrator)
    if integrator.t<SWITCHTIMES[1]*60
        # Stage 1
        INDEX_REFRESH = [idx_TX,idx_TL,idx_dG,idx_pT]
        CONC_REFRESH = [TX0,TL0,dG0,pT0]
        for j in 1:NSPECIES
            integrator.u[j] = integrator.u[j]*(1-DIL_FRAC)
        end
        for j in 1:size(INDEX_REFRESH)[1]
            integrator.u[INDEX_REFRESH[j]] = integrator.u[INDEX_REFRESH[j]] + DIL_FRAC*CONC_REFRESH[j]
        end

    elseif SWITCHTIMES[1]*60<=integrator.t<SWITCHTIMES[2]*60
        # Stage 2
        INDEX_REFRESH = [idx_TX,idx_TL,idx_dG,idx_pT]
        CONC_REFRESH = [TX0,TL0,dG0,pT0]
        for j in 1:NSPECIES
            integrator.u[j] = integrator.u[j]*(1-DIL_FRAC)
        end
        for j in 1:size(INDEX_REFRESH)[1]
            integrator.u[INDEX_REFRESH[j]] = integrator.u[INDEX_REFRESH[j]] + DIL_FRAC*CONC_REFRESH[j]
        end

    elseif SWITCHTIMES[2]*60<=integrator.t
        # Stage 3
        INDEX_REFRESH = [idx_TX,idx_TL,idx_dG,idx_pT]
        CONC_REFRESH = [TX0,TL0,dG0,pT0]
        for j in 1:NSPECIES
            integrator.u[j] = integrator.u[j]*(1-DIL_FRAC)
        end
        for j in 1:size(INDEX_REFRESH)[1]
            integrator.u[INDEX_REFRESH[j]] = integrator.u[INDEX_REFRESH[j]] + DIL_FRAC*CONC_REFRESH[j]
        end

    end
end

function diluteNC_RInl!(integrator)
    if integrator.t<SWITCHTIMES[1]*60
        # Stage 1
        INDEX_REFRESH = [idx_TX,idx_TL,idx_dG,idx_pT]
        CONC_REFRESH = [TX0,TL0,dG0,pT0]
        for j in 1:NSPECIES
            integrator.u[j] = integrator.u[j]*(1-DIL_FRAC)
        end
        for j in 1:size(INDEX_REFRESH)[1]
            integrator.u[INDEX_REFRESH[j]] = integrator.u[INDEX_REFRESH[j]] + DIL_FRAC*CONC_REFRESH[j]
        end

    elseif SWITCHTIMES[1]*60<=integrator.t<SWITCHTIMES[2]*60
        # Stage 2
        INDEX_REFRESH = [idx_TX,idx_TL,idx_dG]
        CONC_REFRESH = [TX0,TL0,dG0]
        for j in 1:NSPECIES
            integrator.u[j] = integrator.u[j]*(1-DIL_FRAC)
        end
        for j in 1:size(INDEX_REFRESH)[1]
            integrator.u[INDEX_REFRESH[j]] = integrator.u[INDEX_REFRESH[j]] + DIL_FRAC*CONC_REFRESH[j]
        end

    elseif SWITCHTIMES[2]*60<=integrator.t
        # Stage 3
        INDEX_REFRESH = [idx_TX,idx_TL,idx_dG]
        CONC_REFRESH = [TX0,TL0,dG0]
        for j in 1:NSPECIES
            integrator.u[j] = integrator.u[j]*(1-DIL_FRAC)
        end
        for j in 1:size(INDEX_REFRESH)[1]
            integrator.u[INDEX_REFRESH[j]] = integrator.u[INDEX_REFRESH[j]] + DIL_FRAC*CONC_REFRESH[j]
        end

    end
end

###############################################
######## 4. Coupled model TL saturation #######
###############################################

function diluteSR_1RI!(integrator)
    if integrator.t<SWITCHTIMES[1]*60
        # Stage 1
        INDEX_REFRESH = [idx_R,idx_TX,idx_TL,idx_dG,idx_pT,idx_dT]
        CONC_REFRESH = [R0,TX0,TL0,dG0,pT0,dT0]
        for j in 1:NSPECIES
            integrator.u[j] = integrator.u[j]*(1-DIL_FRAC)
        end
        for j in 1:size(INDEX_REFRESH)[1]
            integrator.u[INDEX_REFRESH[j]] = integrator.u[INDEX_REFRESH[j]] + DIL_FRAC*CONC_REFRESH[j]
        end

    elseif SWITCHTIMES[1]*60<=integrator.t<SWITCHTIMES[2]*60
        # Stage 2
        INDEX_REFRESH = [idx_R,idx_TX,idx_TL,idx_dG,idx_dT]
        CONC_REFRESH = [R0,TX0,TL0,dG0,dT0]
        for j in 1:NSPECIES
            integrator.u[j] = integrator.u[j]*(1-DIL_FRAC)
        end
        for j in 1:size(INDEX_REFRESH)[1]
            integrator.u[INDEX_REFRESH[j]] = integrator.u[INDEX_REFRESH[j]] + DIL_FRAC*CONC_REFRESH[j]
        end

    elseif SWITCHTIMES[2]*60<=integrator.t
        # Stage 3
        INDEX_REFRESH = [idx_R,idx_TX,idx_TL,idx_dG]
        CONC_REFRESH = [R0,TX0,TL0,dG0]
        for j in 1:NSPECIES
            integrator.u[j] = integrator.u[j]*(1-DIL_FRAC)
        end
        for j in 1:size(INDEX_REFRESH)[1]
            integrator.u[INDEX_REFRESH[j]] = integrator.u[INDEX_REFRESH[j]] + DIL_FRAC*CONC_REFRESH[j]
        end

    end
end

function dilutePC_1RI!(integrator)
    if integrator.t<SWITCHTIMES[1]*60
        # Stage 1
        INDEX_REFRESH = [idx_R,idx_TX,idx_TL,idx_dG,idx_pT]
        CONC_REFRESH = [R0,TX0,TL0,dG0,pT0]
        for j in 1:NSPECIES
            integrator.u[j] = integrator.u[j]*(1-DIL_FRAC)
        end
        for j in 1:size(INDEX_REFRESH)[1]
            integrator.u[INDEX_REFRESH[j]] = integrator.u[INDEX_REFRESH[j]] + DIL_FRAC*CONC_REFRESH[j]
        end

    elseif SWITCHTIMES[1]*60<=integrator.t<SWITCHTIMES[2]*60
        # Stage 2
        INDEX_REFRESH = [idx_R,idx_TX,idx_TL,idx_dG,idx_pT]
        CONC_REFRESH = [R0,TX0,TL0,dG0,pT0]
        for j in 1:NSPECIES
            integrator.u[j] = integrator.u[j]*(1-DIL_FRAC)
        end
        for j in 1:size(INDEX_REFRESH)[1]
            integrator.u[INDEX_REFRESH[j]] = integrator.u[INDEX_REFRESH[j]] + DIL_FRAC*CONC_REFRESH[j]
        end

    elseif SWITCHTIMES[2]*60<=integrator.t
        # Stage 3
        INDEX_REFRESH = [idx_R,idx_TX,idx_TL,idx_dG,idx_pT]
        CONC_REFRESH = [R0,TX0,TL0,dG0,pT0]
        for j in 1:NSPECIES
            integrator.u[j] = integrator.u[j]*(1-DIL_FRAC)
        end
        for j in 1:size(INDEX_REFRESH)[1]
            integrator.u[INDEX_REFRESH[j]] = integrator.u[INDEX_REFRESH[j]] + DIL_FRAC*CONC_REFRESH[j]
        end

    end
end

function diluteNC_1RI!(integrator)
    if integrator.t<SWITCHTIMES[1]*60
        # Stage 1
        INDEX_REFRESH = [idx_R,idx_TX,idx_TL,idx_dG,idx_pT]
        CONC_REFRESH = [R0,TX0,TL0,dG0,pT0]
        for j in 1:NSPECIES
            integrator.u[j] = integrator.u[j]*(1-DIL_FRAC)
        end
        for j in 1:size(INDEX_REFRESH)[1]
            integrator.u[INDEX_REFRESH[j]] = integrator.u[INDEX_REFRESH[j]] + DIL_FRAC*CONC_REFRESH[j]
        end

    elseif SWITCHTIMES[1]*60<=integrator.t<SWITCHTIMES[2]*60
        # Stage 2
        INDEX_REFRESH = [idx_R,idx_TX,idx_TL,idx_dG]
        CONC_REFRESH = [R0,TX0,TL0,dG0]
        for j in 1:NSPECIES
            integrator.u[j] = integrator.u[j]*(1-DIL_FRAC)
        end
        for j in 1:size(INDEX_REFRESH)[1]
            integrator.u[INDEX_REFRESH[j]] = integrator.u[INDEX_REFRESH[j]] + DIL_FRAC*CONC_REFRESH[j]
        end

    elseif SWITCHTIMES[2]*60<=integrator.t
        # Stage 3
        INDEX_REFRESH = [idx_R,idx_TX,idx_TL,idx_dG]
        CONC_REFRESH = [R0,TX0,TL0,dG0]
        for j in 1:NSPECIES
            integrator.u[j] = integrator.u[j]*(1-DIL_FRAC)
        end
        for j in 1:size(INDEX_REFRESH)[1]
            integrator.u[INDEX_REFRESH[j]] = integrator.u[INDEX_REFRESH[j]] + DIL_FRAC*CONC_REFRESH[j]
        end

    end
end

###############################################
####### 5. Coupled model TXTL saturation ######
###############################################

function diluteSR_1RIf!(integrator)
    if integrator.t<SWITCHTIMES[1]*60
        # Stage 1
        INDEX_REFRESH = [idx_R,idx_TX,idx_TL,idx_dG,idx_pT,idx_dT]
        CONC_REFRESH = [R0,TX0,TL0,dG0,pT0,dT0]
        for j in 1:NSPECIES
            integrator.u[j] = integrator.u[j]*(1-DIL_FRAC)
        end
        for j in 1:size(INDEX_REFRESH)[1]
            integrator.u[INDEX_REFRESH[j]] = integrator.u[INDEX_REFRESH[j]] + DIL_FRAC*CONC_REFRESH[j]
        end

    elseif SWITCHTIMES[1]*60<=integrator.t<SWITCHTIMES[2]*60
        # Stage 2
        INDEX_REFRESH = [idx_R,idx_TX,idx_TL,idx_dG,idx_dT]
        CONC_REFRESH = [R0,TX0,TL0,dG0,dT0]
        for j in 1:NSPECIES
            integrator.u[j] = integrator.u[j]*(1-DIL_FRAC)
        end
        for j in 1:size(INDEX_REFRESH)[1]
            integrator.u[INDEX_REFRESH[j]] = integrator.u[INDEX_REFRESH[j]] + DIL_FRAC*CONC_REFRESH[j]
        end

    elseif SWITCHTIMES[2]*60<=integrator.t
        # Stage 3
        INDEX_REFRESH = [idx_R,idx_TX,idx_TL,idx_dG]
        CONC_REFRESH = [R0,TX0,TL0,dG0]
        for j in 1:NSPECIES
            integrator.u[j] = integrator.u[j]*(1-DIL_FRAC)
        end
        for j in 1:size(INDEX_REFRESH)[1]
            integrator.u[INDEX_REFRESH[j]] = integrator.u[INDEX_REFRESH[j]] + DIL_FRAC*CONC_REFRESH[j]
        end

    end
end

function dilutePC_1RIf!(integrator)
    if integrator.t<SWITCHTIMES[1]*60
        # Stage 1
        INDEX_REFRESH = [idx_R,idx_TX,idx_TL,idx_dG,idx_pT]
        CONC_REFRESH = [R0,TX0,TL0,dG0,pT0]
        for j in 1:NSPECIES
            integrator.u[j] = integrator.u[j]*(1-DIL_FRAC)
        end
        for j in 1:size(INDEX_REFRESH)[1]
            integrator.u[INDEX_REFRESH[j]] = integrator.u[INDEX_REFRESH[j]] + DIL_FRAC*CONC_REFRESH[j]
        end

    elseif SWITCHTIMES[1]*60<=integrator.t<SWITCHTIMES[2]*60
        # Stage 2
        INDEX_REFRESH = [idx_R,idx_TX,idx_TL,idx_dG,idx_pT]
        CONC_REFRESH = [R0,TX0,TL0,dG0,pT0]
        for j in 1:NSPECIES
            integrator.u[j] = integrator.u[j]*(1-DIL_FRAC)
        end
        for j in 1:size(INDEX_REFRESH)[1]
            integrator.u[INDEX_REFRESH[j]] = integrator.u[INDEX_REFRESH[j]] + DIL_FRAC*CONC_REFRESH[j]
        end

    elseif SWITCHTIMES[2]*60<=integrator.t
        # Stage 3
        INDEX_REFRESH = [idx_R,idx_TX,idx_TL,idx_dG,idx_pT]
        CONC_REFRESH = [R0,TX0,TL0,dG0,pT0]
        for j in 1:NSPECIES
            integrator.u[j] = integrator.u[j]*(1-DIL_FRAC)
        end
        for j in 1:size(INDEX_REFRESH)[1]
            integrator.u[INDEX_REFRESH[j]] = integrator.u[INDEX_REFRESH[j]] + DIL_FRAC*CONC_REFRESH[j]
        end

    end
end

function diluteNC_1RIf!(integrator)
    if integrator.t<SWITCHTIMES[1]*60
        # Stage 1
        INDEX_REFRESH = [idx_R,idx_TX,idx_TL,idx_dG,idx_pT]
        CONC_REFRESH = [R0,TX0,TL0,dG0,pT0]
        for j in 1:NSPECIES
            integrator.u[j] = integrator.u[j]*(1-DIL_FRAC)
        end
        for j in 1:size(INDEX_REFRESH)[1]
            integrator.u[INDEX_REFRESH[j]] = integrator.u[INDEX_REFRESH[j]] + DIL_FRAC*CONC_REFRESH[j]
        end

    elseif SWITCHTIMES[1]*60<=integrator.t<SWITCHTIMES[2]*60
        # Stage 2
        INDEX_REFRESH = [idx_R,idx_TX,idx_TL,idx_dG]
        CONC_REFRESH = [R0,TX0,TL0,dG0]
        for j in 1:NSPECIES
            integrator.u[j] = integrator.u[j]*(1-DIL_FRAC)
        end
        for j in 1:size(INDEX_REFRESH)[1]
            integrator.u[INDEX_REFRESH[j]] = integrator.u[INDEX_REFRESH[j]] + DIL_FRAC*CONC_REFRESH[j]
        end

    elseif SWITCHTIMES[2]*60<=integrator.t
        # Stage 3
        INDEX_REFRESH = [idx_R,idx_TX,idx_TL,idx_dG]
        CONC_REFRESH = [R0,TX0,TL0,dG0]
        for j in 1:NSPECIES
            integrator.u[j] = integrator.u[j]*(1-DIL_FRAC)
        end
        for j in 1:size(INDEX_REFRESH)[1]
            integrator.u[INDEX_REFRESH[j]] = integrator.u[INDEX_REFRESH[j]] + DIL_FRAC*CONC_REFRESH[j]
        end

    end
end


text:
using DSPHelp
using Base.Test

x = [0.931056  -0.983596  -0.656435  0.400086  -0.778839  -0.0867077  0.121221  -2.58745  0.236455]
evx = [0.583755  -1.78553  -0.267607  0.156689  -0.778839  0.156689  -0.267607  -1.78553  0.583755]
odx = [0.347301  0.80193  -0.388828  0.243397  0.0  -0.243397  0.388828  -0.80193  -0.347301]

@test_approx_eq_eps ev(x) evx 1e-5
@test_approx_eq_eps od(x) odx 1e-5
@test_approx_eq_eps (ev(x) + od(x)) x 1e-5


text:
using Dojo

###############################
### Example benchmark for Atlas
###############################
mech = get_mechanism(:atlas,
    timestep=0.1,
    gravity=-9.81,
    friction_coefficient=0.5,
    damper=25.0,
    spring=1.0,
    contact_feet=true,
    contact_body=true,
    model_type=:simple)

Dojo.initialize!(mech, :atlas_stance,
    body_position=[0.0, 0.0, 1.0],
    body_orientation=[0.0, 0.2, 0.1])

SUITE["atlas"] = @benchmarkable simulate!($mech, 2.5, opts=SolverOptions(rtol=1.0e-6, btol=1.0e-6)) samples=2

###############################
### Example benchmark for Pendulum
###############################
mechanism = get_mechanism(:pendulum,
    timestep=0.01,
    gravity=-9.81,
    damper=5.0,
    spring=0.0)

function controller!(mechanism, t)
    ## Target state
    x_goal = [1.0 * π; 0.0]

    ## Current state
    x = get_minimal_state(mechanism)

    ## Gains
    K = [5.0 0.5] * 0.1

    # Control inputs
    u = -K * (x - x_goal)
    set_input!(mechanism, u)
end

initialize!(mechanism, :pendulum,
    angle=0.0 * π,
    angular_velocity=0.0);

SUITE["pendulum"] = @benchmarkable simulate!($mechanism, 2.0, $controller!) samples=2


text:
abstract type AbstractWindows end
Base.length(w::AbstractWindows) = length(w.ys);

Base.collect(W::AbstractWindows) = [copy.(w) for w in W]

# Window2 ==========================================================
struct Windows2 <: AbstractWindows
    y::AbstractVector
    t::AbstractVector
    ys
    ts
    W
    """
        Windows2(y, t, n=length(y)÷8, noverlap=-1, window_func=rect)

    noverlap = -1 sets the noverlap to n/2

    Iterating `w::Windows2` produces `(y,t)` of specified length. The window is *not* applied to the signal, instead the window array is available as `w.W`.

    # Arguments:
    - `y`: Signal
    - `t`: Time vector
    - `n`: Number of datapoints per window
    - `noverlap`: Overlap between windows
    - `window_func`: Function to create a vector of weights, e.g., `DSP.hanning, DSP.rect` etc.
    """
    function Windows2(y::AbstractArray{T},t,n::Int=length(y)>>3, noverlap::Int=n>>1,window_func=rect) where T

        noverlap < 0 && (noverlap = n>>1)
        N   = length(y)
        @assert N == length(t) "y and t has to be the same length"
        W = T.(window_func(n))
        ys = arraysplit(y,n,noverlap)
        ts = arraysplit(t,n,noverlap)
        new(y,t,ys,ts,W)
    end
end

function Base.iterate(w::Windows2, state=1)
    state > length(w) && return nothing
    ((w.ys[state],w.ts[state]), state+1)
end

"""
    mapwindows(f, W::AbstractWindows)

Apply a Function `f` over all windows represented by `W::Windows2`.
`f` must take `(y,t)->ŷ` where `y` and `ŷ` have the same length.
"""
function mapwindows(f::Function, W::AbstractWindows)
    res = map(f,W)
    merge(res,W)
end

mapwindows(f::Function, args...) = mapwindows(f, Windows2(args...))

function Base.merge(yf::AbstractVector{<:AbstractVector},w::Windows2)
    ym = zeros(eltype(yf[1]), length(w.y))
    counts = zeros(Int, length(w.y))
    dpw = length(w.ys[1])
    inds =  1:dpw

    for i in eachindex(w.ys)
        ym[inds] .+= yf[i]
        counts[inds] .+= 1
        inds = inds .+ (dpw-w.ys.noverlap)
        inds =  inds[1]:min(inds[end], length(ym))
    end
    ym ./= max.(counts,1)
end

# Window3 ==========================================================
"""
    Windows3(y, t, v, n=length(y)÷8, noverlap=-1, window_func=rect)

noverlap = -1 sets the noverlap to n/2

#Arguments:
- `y`: Signal
- `t`: Time vector
- `v`: Auxiliary vector
- `n`: Number of datapoints per window
- `noverlap`: Overlap between windows
- `window_func`: Function to apply over window
"""
struct Windows3 <: AbstractWindows
    y::AbstractVector
    t::AbstractVector
    v::AbstractVector
    ys
    ts
    vs
    W
    function Windows3(y::AbstractVector,t::AbstractVector,v::AbstractVector,n::Int=length(y)>>3, noverlap::Int=n>>1,window_func::Function=rect)
        N       = length(y)
        @assert N == length(t) == length(v) "y, t and v has to be the same length"
        noverlap < 0 && (noverlap = n>>1)
        W = window_func(n)
        ys = arraysplit(y,n,noverlap)
        ts = arraysplit(t,n,noverlap)
        vs = arraysplit(v,n,noverlap)
        new(y,t,v,ys,ts,vs,W)
    end
end

function Base.iterate(w::Windows3, state=1)
    state > length(w.ys) && return nothing
    ((w.ys[state],w.ts[state],w.vs[state]), state+1)
end


text:
using JuMP
using MathOptInterface
const MOI = MathOptInterface

using Gurobi
model = Model(Gurobi.Optimizer)
# set_optimizer_attribute(model, "Cuts",0)
set_optimizer_attribute(model, "TimeLimit", 600)

# using Cbc
# model = Model(Cbc.Optimizer)
# set_optimizer_attribute(model, "cuts","off")
# set_optimizer_attribute(model, "seconds", 60)
# set_optimizer_attribute(model, "thread", 16)

n_i = 12                # 产品种类个数
n_week = 15                # 时间周期数
n_day = 7
n_machine = 3           # 机器个数
machine = 1:n_machine   # 机器编号
product = 1:n_i         # 产品编号
weeks = 1:n_week           # 时间编号
days = 1:n_day          # 星期编号
paras = [   2 83 10 1 30;
            2 31 10 1 20;
            2 11 10 1 30;
            2 93 10 1 40;
            2 82 10 1 40;
            2 72 10 1 10;
            3 23 20 1 30;
            3 91 20 1 20;
            3 83 20 1 10;
            3 34 20 1 50;
            3 61 20 1 30;
            3 82 20 1 20]

f = paras[:,1] # 产品包装使用的机器编号、产品类别

f2 = product[f .== 2] # 谷物类产品
f3 = product[f .== 3] # 水果类产品

s_init = paras[:,2] # 初始库存
ss = paras[:, 3] # 安全库存
alpha = repeat(paras[:, 4], 1, n_machine) # 生产速度
beta = paras[:, 5] # 清洗时间
# 需求
d  = [   0 95 110 96 86 124 83 108 114 121 110 124 104 86 87;
        98 96 96 98 103 104 122 101 89 108 101 109 106 108 76;
        106 0 89 123 96 105 83 82 112 109 119 85 99 80 123;
        98 121 0 105 98 96 101 81 117 76 103 81 95 105 102;
        0 124 113 123 123 79 111 98 97 80 98 124 78 108 109;
        103 102 0 95 107 105 107 105 75 93 115 113 111 105 85;
        110 93 0 112 84 124 98 101 83 87 105 118 115 106 78;
        85 92 101 110 93 96 120 109 121 87 92 85 91 93 109;
        122 116 109 0 105 108 88 98 77 90 110 102 107 99 96;
        120 124 94 105 92 86 101 106 75 109 83 95 79 108 100;
        117 96 78 0 108 87 114 107 110 94 104 101 108 110 80;
        125 112 75 0 116 103 122 88 85 84 76 102 84 88 82]
L = [1400, 700, 700] # 机器上限

@variable(model, x[i in product, t in weeks, day in days] >= 0) # 产品生产计划
@variable(model, s[i in product, t in weeks] >= 0) # 库存
@variable(model, y[i in product, t in weeks,day in days], Bin) # 是否开工生产

# 库存平衡
@constraint(model ,dem_sat1[i in product], s_init[i] + sum(x[i, 1, day] for day in days) == d[i, 1] + s[i, 1])
@constraint(model, dem_sat[i in product, t in weeks[2:end]], s[i, t - 1] + sum(x[i, t, day] for day in days) == d[i, t] + s[i, t] )
# 安全库存
@constraint(model, ss_limit[i in product, t in weeks], s[i, t] >= ss[i])

# 最大产量约束
@constraint(model, vub[i in product, t in weeks, day in days], x[i, t, day] <= (sum(d[i, k]  for k in t:n_week) + ss[i]) * y[i, t, day])

# 混合机器上限
@constraint(model, mix_cap[i in product, t in weeks, day in days], alpha[i, 1] * x[i, t, day]  + beta[i] * y[i, t, day]  <= L[1] / 7)

# 包装机器上限
@constraint(model, pack_cap2[i in f2, t in weeks, day in days], alpha[i, 2] * x[i, t, day]  <= L[2] / 7)
@constraint(model, pack_cap3[i in f3, t in weeks, day in days], alpha[i, 3] * x[i, t, day]  <= L[3] / 7)

# 包装机器只能开启一台
@constraint(model, pcak_work2[t in weeks, day in days], sum(y[i, t, day] for i in f2) <= 1)
@constraint(model, pcak_work3[t in weeks, day in days], sum(y[i, t, day] for i in f3) <= 1)
@expression(model, inventory, sum(s[i, t] for i in product, t in weeks))
@objective(model, Min, inventory)

optimize!(model)

max_stars_repo_name: Behemoth-s/ProductPlanningByMIP

max_issues_repo_name: Behemoth-s/ProductPlanningByMIP

max_forks_repo_name: Behemoth-s/ProductPlanningByMIP


text:
struct Poisoned end

mutable struct ExprUnit
    items::Vector
    deps::OrderedSet{Symbol}
    state::Symbol
end

function ExprUnit(a::Array, deps=[]; state::Symbol=:new)
    ExprUnit(a, OrderedSet{Symbol}([target_type(dep) for dep in deps]), state)
end
function ExprUnit(e::Union{Expr,Symbol,String,Poisoned}, deps=[]; state::Symbol=:new)
    ExprUnit([e], OrderedSet{Symbol}([target_type(dep) for dep in deps]), state)
end
ExprUnit() = ExprUnit([], OrderedSet{Symbol}(), :new)

"""
Pretty-print a buffer of expressions (and comments) to an output stream
Adds blank lines at appropriate places for readability
"""
function print_buffer(out_stream, out_buffer)
    state = :comment
    for e in out_buffer
        isa(e, Poisoned) && continue
        prev_state = state
        state = (isa(e, AbstractString) ? :string :
                 isa(e, Expr) ? e.head :
                 error("output: can't handle type $(typeof(e))"))

        if state == :string && startswith(e, "# Skipping")
            state = :skipping
        end

        if ((state != prev_state && prev_state != :string) ||
            (state == prev_state && (state == :function || state == :struct)))
            println(out_stream)
        end

        if isa(e, Expr) && e.head == :macrocall && first(e.args) == Symbol("@cenum")
            println(out_stream, "@cenum $(e.args[3]) begin")
            for elem in e.args[4].args
                println(out_stream, "    $elem")
            end
            println(out_stream, "end")
            println(out_stream)
            continue
        end
        println(out_stream, e)
    end
end

function dump_to_buffer!(buffer::Vector, expr_dict::OrderedDict{Symbol,ExprUnit}, item::ExprUnit)
    (item.state == :done || item.state == :processing) && return buffer
    item.state = :processing
    for dep in item.deps
        haskey(expr_dict, dep) || continue
        dump_to_buffer!(buffer, expr_dict, expr_dict[dep])
    end
    append!(buffer, item.items)
    item.state = :done
    return buffer
end

function dump_to_buffer(expr_dict::OrderedDict{Symbol,ExprUnit})
    buffer = []
    for item in values(expr_dict)
        item.state == :done && continue
        dump_to_buffer!(buffer, expr_dict, item)
    end
    return buffer
end


text:
# Figure 3
table31 = readtable("csv/fig_3.1.csv");
table31[:line] = linspace(0,300,size(table31)[1])
p_fig3_l = plot(table31,
    layer(x=:λ_E, y=:line, Geom.line, Theme(default_color=colorant"green")),
    layer(x=:λ_H, y=:tr_prioH , Geom.line, Theme(default_color=colorant"blue")),
    layer(x=:λ_H, y=:tr_prioE, Geom.line, Theme(default_color=colorant"red")),
    Guide.manual_color_key("", ["BH", "BE",
                                "λE = 5"],
                               [colorant"blue", colorant"red", colorant"green"]),
    Guide.YLabel("wH", orientation=:horizontal),
    Guide.XLabel("λH"),
    #Guide.Title("C(d) vs B_H, vs B_E, p_E = $p_E, d = $d, p_H = $p_H, λ_H = $λ_H")
    my_theme
)
draw(PDF("figs/fig_3.1.pdf", 6inch, 4.5inch), p_fig3_l)

table32 = readtable("csv/fig_3.2.csv");
table32[:line] = linspace(0,0.8,size(table32)[1])
p_fig3_l = plot(table32,
    layer(x=:λ_E, y=:line, Geom.line, Theme(default_color=colorant"green")),
    layer(x=:λ_H, y=:tr_prioH , Geom.line, Theme(default_color=colorant"blue")),
    layer(x=:λ_H, y=:tr_prioE, Geom.line, Theme(default_color=colorant"red")),
    Guide.manual_color_key("", ["BH", "BE",
                                "λE = 5"],
                               [colorant"blue", colorant"red", colorant"green"]),
    Guide.YLabel("wH", orientation=:horizontal),
    Guide.XLabel("λH"),
    #Guide.Title("C(d) vs B_H, vs B_E, p_E = $p_E, d = $d, p_H = $p_H, λ_H = $λ_H")
    my_theme
)
draw(PDF("figs/fig_3.2.pdf", 6inch, 4.5inch), p_fig3_l)

# table31 = readtable("csv/fig_3.1.csv");
# table31[:line] = linspace(0,300,size(table31)[1])
# p_fig3_l = plot(table31,
#     layer(x=:λ_E, y=:line, Geom.line, Theme(default_color=colorant"green")),
#     layer(x=:λ_H, y=:tr_prioH , Geom.line, Theme(default_color=colorant"blue")),
#     layer(x=:λ_H, y=:tr_prioE, Geom.line, Theme(default_color=colorant"red")),
#     layer(x=:λ_H, y=:tr_random_prio, Geom.line, Theme(default_color=colorant"purple")),
#     Guide.manual_color_key("", ["BH", "BE",
#                                 "Brandom", "λE = 5"],
#                                [colorant"blue", colorant"red", colorant"purple", colorant"green"]),
#     Guide.YLabel("wH", orientation=:horizontal),
#     Guide.XLabel("λH"),
#     #Guide.Title("C(d) vs B_H, vs B_E, p_E = $p_E, d = $d, p_H = $p_H, λ_H = $λ_H")
#     my_theme
# )
# draw(PDF("figs/fig3_l_random_prio2.pdf", 6inch, 4.5inch), p_fig3_l)

# p_fig3_bis = plot(table31,
#     layer(x=:λ_H, y=:tr_chains_d1, Geom.line, Theme(default_color=colorant"purple")),
#     layer(x=:λ_H, y=:tr_chains_d50, Geom.line, Theme(default_color=colorant"green")),
#     layer(x=:λ_H, y=:tr_prioH , Geom.line, Theme(default_color=colorant"blue")),
#     Guide.manual_color_key("", ["C(1)", "C(20)", "BH"],
#                                [colorant"purple",colorant"green", colorant"blue"]),
#     Guide.YLabel("wH", orientation=:horizontal),
#     Guide.XLabel("λH"),
#     #Guide.Title("C(d) vs B_H, vs B_E, p_E = $p_E, d = $d, p_H = $p_H, λ_H = $λ_H")
#     my_theme
# )
# draw(PDF("figs/fig3_bis.pdf", 6inch, 4.5inch), p_fig3_bis)

# table32 = readtable("csv/fig_3.2.csv")
# p_fig3b_r = plot(table32,
#     layer(x=:λ_H, y=:tr_chains_d1, Geom.line, Theme(default_color=colorant"purple")),
#     layer(x=:λ_H, y=:tr_chains_d10, Geom.line, Theme(default_color=colorant"green")),
#     layer(x=:λ_H, y=:tr_prioH , Geom.line, Theme(default_color=colorant"blue")),
#     layer(x=:λ_H, y=:tr_prioE, Geom.line, Theme(default_color=colorant"red")),
#     Guide.manual_color_key("", ["C(1)", "C(10)", "BH", "BE"],
#         [colorant"purple",colorant"green", colorant"blue", colorant"red"]),
#     Guide.YLabel("wH", orientation=:horizontal),
#     Guide.XLabel("λE"),
#     #Guide.Title("C(d) vs B_H, vs B_E, p_E = $p_E, d = $d, p_H = $p_H, λ_H = $λ_H")
#     my_theme
# )
# draw(PDF("figs/fig3b_r.pdf", 6inch, 4.5inch), p_fig3b_r)


text:
using POMDPModels
using POMDPs
using POMDPSimulators
using POMDPModelTools
using BeliefUpdaters
using POMDPPolicies
using POMDPTesting
using Random
using Test
using Compose

@testset "crying" begin
    include("crying.jl")
end
@testset "gridworld" begin
    include("simple_gridworld.jl")
end
@testset "gridworld_vis" begin
    include("test_gw_vis.jl")
end
@testset "legacy_gridworld" begin
    include("legacy_gridworld.jl")
end
@testset "tiger" begin
    include("tiger.jl")
end
@testset "random" begin
    include("random.jl")
end
@testset "tmaze" begin
    include("tmaze.jl")
end
@testset "inverted" begin
    include("inverted.jl")
end
@testset "car" begin
    include("car.jl")
end
@testset "lightdark" begin
    include("lightdark.jl")
end
@testset "minihallway" begin
    include("minihallway.jl")
end


text:
using SimpleDelimitedFiles
using Documenter

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

makedocs(;
    modules = [SimpleDelimitedFiles],
    authors = "Jishnu Bhattacharya",
    repo = "https://github.com/jishnub/SimpleDelimitedFiles.jl/blob/{commit}{path}#L{line}",
    sitename = "SimpleDelimitedFiles",
    format = Documenter.HTML(;
        prettyurls = get(ENV, "CI", "false") == "true",
        canonical = "https://jishnub.github.io/SimpleDelimitedFiles.jl",
        assets = String[]
    ),
    pages = [
        "Reference" => "index.md",
    ]
)

deploydocs(;
    repo = "github.com/jishnub/SimpleDelimitedFiles.jl",
    devbranch = "main",
)


text:
# This file contains functions which help creating a Project.

"""
This function returns a dataframe containing project data if the data file exists,
otherwise it returns nothing.
"""
function extract_projectdata(inv_dir::String,
                        file_name::String)

    file_path = joinpath(inv_dir, file_name)
    if isfile(file_path)
        projectdata = read_data(file_path)
        return projectdata
    else
        return nothing
    end

end

"""
This function returns the queue cost data for a project
based on its size.
"""
function create_project_queuecost(size::Float64,
                                queue_cost_data::DataFrames.DataFrame,
                                lag_bool::Bool)
    if lag_bool
        D1_row = findfirst(x -> x == "D1", queue_cost_data[:, "Ref"])
        D2_row = findfirst(x -> x == "D2", queue_cost_data[:, "Ref"])
        M2_row = findfirst(x -> x == "M2 (per MW)", queue_cost_data[:, "Ref"])

        if size < 6
            D1 = queue_cost_data[D1_row, "less than 6 MW"]
            D2 = queue_cost_data[D2_row, "less than 6 MW"]
            M2 = queue_cost_data[M2_row, "less than 6 MW"] * size
        elseif 6 <= size <= 20
            D1 = queue_cost_data[D1_row, "6 to 20 MW"]
            D2 = queue_cost_data[D2_row, "6 to 20 MW"]
            M2 = queue_cost_data[M2_row, "6 to 20 MW"] * size
        elseif 20 < size <= 50
            D1 = queue_cost_data[D1_row, "20 to 50 MW"]
            D2 = queue_cost_data[D2_row, "20 to 50 MW"]
            M2 = queue_cost_data[M2_row, "20 to 50 MW"] * size
        elseif 50 < size <= 100
            D1 = queue_cost_data[D1_row, "50 to 100 MW"]
            D2 = queue_cost_data[D2_row, "50 to 100 MW"]
            M2 = queue_cost_data[M2_row, "50 to 100 MW"] * size
        elseif 100 < size <= 200
            D1 = queue_cost_data[D1_row, "100 to 200 MW"]
            D2 = queue_cost_data[D2_row, "100 to 200 MW"]
            M2 = queue_cost_data[M2_row, "100 to 200 MW"] * size
        elseif 200 < size <= 500
            D1 = queue_cost_data[D1_row, "200 to 500 MW"]
            D2 = queue_cost_data[D2_row, "200 to 500 MW"]
            M2 = queue_cost_data[M2_row, "200 to 500 MW"] * size
        elseif 500 < size < 1000
            D1 = queue_cost_data[D1_row, "500 to 1000 MW"]
            D2 = queue_cost_data[D2_row, "500 to 1000 MW"]
            M2 = queue_cost_data[M2_row, "500 to 1000 MW"] * size
        elseif size >= 1000
            D1 = queue_cost_data[D1_row, "greater than 1000 MW"]
            D2 = queue_cost_data[D2_row, "greater than 1000 MW"]
            M2 = queue_cost_data[M2_row, "greater than 1000 MW"] * size
        end
    else
        D1 = 0; D2 = 0; M2 = 0
    end
    queue_cost = [D1 + D2, M2]
    return queue_cost
end

"""
This function writes the project's hourly availability data
to the "DAY_AHEAD_availability.csv" file stored in
simulation_dir/timeseries_data_files/Availability
"""
function add_investor_project_availability!(simulation_dir::String,
                                              projects::Vector{Project})
    system_availability_data = DataFrames.DataFrame(CSV.File(joinpath(simulation_dir, "timeseries_data_files", "Availability", "DAY_AHEAD_availability.csv")))

    for project in projects
        project_name = get_name(project)
        tech = get_tech(project)
        type_zone_id = "$(get_type(tech))_$(get_zone(tech))"

        if (!in(project_name, names(system_availability_data)) && !in(type_zone_id, names(system_availability_data)))
            if typeof(project) == ThermalGenEMIS{Existing} || typeof(project) == BatteryEMIS{Existing}
                system_availability_data[:, project_name] = ones(DataFrames.nrow(system_availability_data))
            elseif typeof(project) == RenewableGenEMIS{Existing}
                error("Timeseries for existing renewable generator not provided")
            elseif typeof(project) == ThermalGenEMIS{Option} || typeof(project) == BatteryEMIS{Option}
                system_availability_data[:, type_zone_id] = ones(DataFrames.nrow(system_availability_data))
            elseif typeof(project) == RenewableGenEMIS{Option}
                availability = zeros(DataFrames.nrow(system_availability_data))
                bus = get_bus(tech)
                if get_type(tech) == "WT"
                    gens_in_zone = filter(g -> (first("$(bus)", 1) == first(g, 1)) && (occursin("wind", lowercase(g)) || occursin("wt", lowercase(g))),
                                                    names(system_availability_data))
                elseif get_type(tech) == "PVe"
                    gens_in_zone = filter(g -> (first("$(bus)", 1)== first(g, 1)) && occursin("pv", lowercase(g)),
                                                    names(system_availability_data))
                elseif get_type(tech) == "HY"
                    gens_in_zone = filter(g -> (first("$(bus)", 1) == first(g, 1)) && occursin("hy", lowercase(g)),
                                                    names(system_availability_data))
                end

                if length(gens_in_zone) < 1
                    error("Similar renewable generator in zone not found")
                else
                    for g in gens_in_zone
                        availability += system_availability_data[:, g]
                    end
                    system_availability_data[:, type_zone_id] = availability / length(gens_in_zone)
                end
            end
        end
    end

    write_data(joinpath(simulation_dir, "timeseries_data_files", "Availability"), "DAY_AHEAD_availability.csv", system_availability_data)

    return
end

"""
This function creates a piece-wise linear operation cost for the project according to the structure used in PSY.
This allows smooth conversion of local Project structs to PSY Device structs.
Returns PSY.TwoPartCost or PSY.ThreePartCost based on the type of the Project.
"""
function create_operation_cost(projectdata::DataFrames.DataFrameRow, size::Float64, base_power::Float64)

    output_point_fields = String[]
    heat_rate_fields = String[]

    fields = names(projectdata)

    for field in fields
        if occursin("Output_pct_", field)
            push!(output_point_fields, field)
        elseif occursin("HR_", field)
            push!(heat_rate_fields, field)
        end
    end

    @assert length(output_point_fields) > 0
    cost_colnames = zip(heat_rate_fields, output_point_fields)

    fuel_cost =  projectdata["Fuel Price \$/MMBTU"] / 1000.0
    var_cost = [(projectdata[hr], projectdata[mw]) for (hr, mw) in cost_colnames]

    var_cost = unique([
        (tryparse(Float64, string(c[1])), tryparse(Float64, string(c[2])))
        for c in var_cost if !in("NA", c)
    ])

    if length(var_cost) > 1
        var_cost[2:end] = [
            (
                var_cost[i][1] *
                (var_cost[i][2] - var_cost[i - 1][2]) *
                fuel_cost,
                var_cost[i][2] / base_power,
            ) .* size for i in 2:length(var_cost)
        ]
        var_cost[1] =
            (var_cost[1][1] * var_cost[1][2] * fuel_cost, var_cost[1][2] / base_power) .* size

        fixed = max(
            0.0,
            var_cost[1][1] -
            (var_cost[2][1] / (var_cost[2][2] - var_cost[1][2]) * var_cost[1][2]),
        )
        var_cost[1] = (var_cost[1][1] - fixed, var_cost[1][2])
        for i in 2:length(var_cost)
            var_cost[i] = (var_cost[i - 1][1] + var_cost[i][1], var_cost[i][2])
        end
    elseif length(var_cost) == 1
        # if there is only one point, use it to determine the constant $/MW cost
        var_cost = var_cost[1][1] * var_cost[1][2] * fuel_cost * base_power
        fixed = 0.0
    else
        var_cost = 0.0
        fixed = 0.0
    end

    start_up_cost = fuel_cost * projectdata["Start Heat Cold MBTU"] * 1000.0
    shut_down_cost = 0.0

    if occursin("CC", projectdata["Unit Type"]) || occursin("CT", projectdata["Unit Type"]) || occursin("ST", projectdata["Unit Type"]) || occursin("NUCLEAR", projectdata["Unit Type"])
        operation_cost = PSY.ThreePartCost(var_cost, fixed, start_up_cost, shut_down_cost)
    else
        operation_cost = PSY.TwoPartCost(var_cost, fixed)
    end

    return operation_cost
end

"""
This function populates project investment parameters and creates Finance struct.
"""
function create_investment_data(size::Float64,
                              projectdata::DataFrames.DataFrameRow,
                              investor_dir::String,
                              simulation_data::AgentSimulationData,
                              products::Vector{<: Product},
                              investor_name::String,
                              scenario_names::Vector{String},
                              lag_bool::Bool
                              )
    start_year = get_start_year(get_case(simulation_data))
    max_horizon = get_total_horizon(get_case(simulation_data))

    decision_year = 1
    queue_cost_data = get_queue_cost_data(simulation_data)
    queue_cost = create_project_queuecost(size, queue_cost_data, lag_bool)
    queue_time = length(queue_cost)
    construction_year = decision_year + (queue_time + projectdata["Lagtime"]) * lag_bool

    capex_years = projectdata["Capex Years"]

    capex_data = read_data(joinpath(investor_dir, "project_capex.csv"))

    category = projectdata["Category"]
    capex_row = findfirst(x-> x == category, capex_data[:, "Category"])

    if lag_bool
        online_year = start_year + decision_year + queue_time - 1
        age = 0
        remaining_lifetime = projectdata["Lifetime"]
        project_capex = [capex_data[capex_row, "$(start_year + y - 1)"] for y in 1:max_horizon]
    else
        online_year = projectdata["Online Year"]
        age = start_year - online_year
        remaining_lifetime = projectdata["Lifetime"] - age
        project_capex = zeros(max_horizon)

        for y = 1:min(max_horizon, (capex_years - age))
            project_capex[y] = capex_data[capex_row, "pre $(start_year)"]
        end

    end

    investment_cost, discount_rate = calculate_capcost_and_wacc(project_capex, category, investor_dir, online_year)

    investment_cost = investment_cost * size * 1000.0

    effective_investment_cost = investment_cost[queue_time + 1]
    fixedOM_cost = projectdata["Fixed OM Cost per MW"] * size

    end_life_year = (queue_time + projectdata["Lagtime"]) * lag_bool + remaining_lifetime

    retirement_year = (queue_time + projectdata["Lagtime"]) * lag_bool + remaining_lifetime

    num_profit_years = max(max_horizon, end_life_year)

    yearly_profit = [AxisArrays.AxisArray(zeros(length(products), num_profit_years),
                                                 AxisArrays.Axis{:prod}(get_name.(products)),
                                                 AxisArrays.Axis{:year}(1:1:num_profit_years))
                                                 for y in 1:num_profit_years]

    scenario_profit = Dict(scenario_name => [AxisArrays.AxisArray(zeros(length(products), num_profit_years),
                                                 AxisArrays.Axis{:prod}(get_name.(products)),
                                                 AxisArrays.Axis{:year}(1:1:num_profit_years))
                                                 for y in 1:num_profit_years]  for scenario_name in scenario_names)

    realized_profit = AxisArrays.AxisArray(zeros(length(products), num_profit_years),
                                                 AxisArrays.Axis{:prod}(get_name.(products)),
                                                 AxisArrays.Axis{:year}(1:1:num_profit_years))

    scenario_npv = Dict(scenario_name => zeros(num_profit_years) for scenario_name in scenario_names)
    expected_npv = zeros(num_profit_years)
    scenario_utility = Dict(scenario_name => zeros(num_profit_years) for scenario_name in scenario_names)
    expected_utility = zeros(num_profit_years)

    annual_cashflow = zeros(num_profit_years)

    finance_data = Finance(investment_cost,
                           effective_investment_cost,
                           projectdata["Lagtime"],
                           remaining_lifetime,
                           capex_years,
                           fixedOM_cost,
                           queue_cost,
                           scenario_profit,
                           realized_profit,
                           discount_rate,
                           scenario_npv,
                           expected_npv,
                           scenario_utility,
                           expected_utility,
                           annual_cashflow,
                           investor_name)

    return decision_year, construction_year, retirement_year, end_life_year, finance_data
end

"""
This function creates the Tech struct for projects
based only on the CSV data.
"""
function create_tech_type(name::String,
                          projectdata::DataFrames.DataFrameRow,
                          size::Float64,
                          base_power::Float64,
                          decision_year::Int64,
                          construction_year::Int64,
                          retirement_year::Int64,
                          end_life_year::Int64,
                          products::Vector{Product},
                          finance_data::Finance,
                          sys_UC_bool::Bool
                        )

    type = projectdata["Unit Type"]
    min_cap = projectdata["Min Gen pu"] * size

    active_power_limits = (min = min_cap, max = size)
    ramp_limits = (up = projectdata["Ramp Rate pu/Hr"]  * size, down = projectdata["Ramp Rate pu/Hr"] * size)
    zone = projectdata["Zone"]
    bus = projectdata["Bus ID"]

    if sys_UC_bool
        operation_cost = create_operation_cost(projectdata, size, base_power)
        up_down_time = (up = projectdata["Min Up Time Hr"], down = projectdata["Min Down Time Hr"])
    else
        operation_cost = nothing
        up_down_time = nothing
    end

    if type == "ST" || type == "CT" || type == "CC" || type == "NU_ST"
        tech = ThermalTech(type,
                       projectdata["Fuel"],
                       active_power_limits,
                       ramp_limits,
                       up_down_time,
                       operation_cost,
                       bus,
                       zone)

        return  ThermalGenEMIS{BuildPhase}(name,
                                 tech,
                                 decision_year,
                                 construction_year,
                                 retirement_year,
                                 end_life_year,
                                 products,
                                 finance_data)
    elseif type == "WT" || type == "PVe"
        tech = RenewableTech(type,
                        active_power_limits,
                        ramp_limits,
                        operation_cost,
                        bus,
                        zone)

        return  RenewableGenEMIS{BuildPhase}(name,
                                 tech,
                                 decision_year,
                                 construction_year,
                                 retirement_year,
                                 end_life_year,
                                 products,
                                 finance_data)

    elseif type == "HY"
        tech = HydroTech(type,
                       active_power_limits,
                       ramp_limits,
                       up_down_time,
                       operation_cost,
                       bus,
                       zone)

        return  HydroGenEMIS{BuildPhase}(name,
                                 tech,
                                 decision_year,
                                 construction_year,
                                 retirement_year,
                                 end_life_year,
                                 products,
                                 finance_data)

    elseif type == "BA"
        storage_capacity = parse(Float64, projectdata["Duration Hr"]) * size
        tech = BatteryTech(type,
                        (min = 0.0, max = parse(Float64, projectdata["Input Power Rating pu"]) * size),
                        (min = 0.0, max = parse(Float64, projectdata["Output Power Rating pu"]) * size),
                        ramp_limits,
                        (min = parse(Float64, projectdata["Min Storgae pu"]) * storage_capacity, max = storage_capacity),
                        0.0,
                        (in = parse(Float64, projectdata["Round Trip Efficiency pu"]), out = 1.0),
                        bus,
                        zone)
        return BatteryEMIS{BuildPhase}(name,
                        tech,
                        decision_year,
                        construction_year,
                        retirement_year,
                        end_life_year,
                        products,
                        finance_data)

    end

end

"""
This function creates the Tech struct for ThermalGen projects
based on the CSV and PSY Device data.
"""
function create_tech_type(name::String,
                          device::T,
                          projectdata::DataFrames.DataFrameRow,
                          size::Float64,
                          base_power::Float64,
                          decision_year::Int64,
                          construction_year::Int64,
                          retirement_year::Int64,
                          end_life_year::Int64,
                          products::Vector{Product},
                          finance_data::Finance
                        ) where T <: PSY.ThermalGen

    min_cap = PSY.get_active_power_limits(device)[:min] * base_power

    bus = PSY.get_bus(device)
    device_ramp_limits = PSY.get_ramp_limits(device)

    if isnothing(device_ramp_limits)
        ramp_limits = nothing
    else
        ramp_limits = (up = device_ramp_limits[:up] * base_power * 60, down = device_ramp_limits[:down] * base_power * 60)
    end

    prime_mover = string(PSY.get_prime_mover(device))
    fuel = string(PSY.get_fuel(device))

    if fuel == "NUCLEAR"
        prime_mover = "NU_ST"
    end

    tech = ThermalTech(prime_mover,
                       fuel,
                       (min = min_cap, max = size),
                       ramp_limits,
                       PSY.get_time_limits(device),
                       PSY.get_operation_cost(device),
                       PSY.get_number(bus),
                       projectdata["Zone"])

    return  ThermalGenEMIS{BuildPhase}(name,
                                 tech,
                                 decision_year,
                                 construction_year,
                                 retirement_year,
                                 end_life_year,
                                 products,
                                 finance_data)

end

"""
This function creates the Tech struct for RenewableGen projects
based on the CSV and PSY Device data.
"""
function create_tech_type(name::String,
                          device::T,
                          projectdata::DataFrames.DataFrameRow,
                          size::Float64,
                          base_power::Float64,
                          decision_year::Int64,
                          construction_year::Int64,
                          retirement_year::Int64,
                          end_life_year::Int64,
                          products::Vector{Product},
                          finance_data::Finance
                        ) where T <: PSY.RenewableGen

    bus = PSY.get_bus(device)

    type = string(PSY.get_prime_mover(device))

    tech = RenewableTech(type,
                       (min = 0.0, max = size),
                       (up = size, down = size),
                       PSY.get_operation_cost(device),
                       PSY.get_number(bus),
                       projectdata["Zone"])

    return  RenewableGenEMIS{BuildPhase}(name,
                                 tech,
                                 decision_year,
                                 construction_year,
                                 retirement_year,
                                 end_life_year,
                                 products,
                                 finance_data)

end

"""
This function creates the Tech struct for HydroGen projects
based on the CSV and PSY Device data.
"""
function create_tech_type(name::String,
                          device::P,
                          projectdata::DataFrames.DataFrameRow,
                          size::Float64,
                          base_power::Float64,
                          decision_year::Int64,
                          construction_year::Int64,
                          retirement_year::Int64,
                          end_life_year::Int64,
                          products::Vector{Product},
                          finance_data::Finance
                        ) where P <: PSY.HydroGen

    min_cap = PSY.get_active_power_limits(device)[:min] * base_power
    bus = PSY.get_bus(device)

    device_ramp_limits = PSY.get_ramp_limits(device)

    if isnothing(device_ramp_limits)
        ramp_limits = nothing
    else
        ramp_limits = (up = device_ramp_limits[:up] * base_power * 60, down = device_ramp_limits[:down] * base_power * 60)
    end

    tech = HydroTech("HY",
                       (min = min_cap, max = size),
                       ramp_limits,
                       PSY.get_time_limits(device),
                       PSY.get_operation_cost(device),
                       PSY.get_number(bus),
                       projectdata["Zone"])

        return  HydroGenEMIS{BuildPhase}(name,
                                 tech,
                                 decision_year,
                                 construction_year,
                                 retirement_year,
                                 end_life_year,
                                 products,
                                 finance_data)

end

"""
This function creates the Tech struct for Battery projects
based on the CSV and PSY Device data.
"""
function create_tech_type(name::String,
                          device::PSY.GenericBattery,
                          projectdata::DataFrames.DataFrameRow,
                          size::Float64,
                          base_power::Float64,
                          decision_year::Int64,
                          construction_year::Int64,
                          retirement_year::Int64,
                          end_life_year::Int64,
                          products::Vector{Product},
                          finance_data::Finance
                        )

    bus = PSY.get_bus(device)

    type = string(PSY.get_prime_mover(device))
    input_active_power_limits = PSY.get_input_active_power_limits(device)
    output_active_power_limits = PSY.get_output_active_power_limits(device)
    storage_capacity = PSY.get_state_of_charge_limits(device)

    tech = BatteryTech(type,
                       (min = input_active_power_limits[:min] * base_power, down = input_active_power_limits[:max] * base_power),
                       (min = output_active_power_limits[:min] * base_power, down = output_active_power_limits[:max] * base_power),
                       (up = size, down = size),
                       (min = storage_capacity[:min] * size * base_power, down = storage_capacity[:max] * size * base_power),
                       PSY.get_initial_energy(device) * size * base_power,
                       PSY.get_efficiency(device),
                       PSY.get_number(bus),
                       projectdata["Zone"])

    return  BatteryEMIS{BuildPhase}(name,
                                 tech,
                                 decision_year,
                                 construction_year,
                                 retirement_year,
                                 end_life_year,
                                 products,
                                 finance_data)

end

"""
This function returns a Project struct based only on project's CSV data.
"""
function create_project(projectdata::DataFrames.DataFrameRow,
                          simulation_data::AgentSimulationData,
                          investor_name::String,
                          investor_dir::String,
                          scenario_names::Vector{String},
                          lag_bool::Bool)

    name = projectdata["GEN_UID"]
    size = Float64(size_in_MW(investor_dir,
                     projectdata["Unit Type"],
                     projectdata["Size"]))

    base_power = size

    products = create_products(simulation_data,
                               projectdata)

    decision_year,
    construction_year,
    retirement_year,
    end_life_year,
    finance_data = create_investment_data(size, projectdata, investor_dir, simulation_data, products, investor_name, scenario_names, lag_bool)

    sys_UC_bool = !isnothing(get_system_UC(simulation_data))

    project = create_tech_type(name, projectdata, size, base_power, decision_year, construction_year, retirement_year, end_life_year, products, finance_data, sys_UC_bool)

    return project
end

"""
This function returns a Project struct based on CSV data and PSY data.
"""
function create_project(projectdata::DataFrames.DataFrameRow,
                          device::T,
                          simulation_data::AgentSimulationData,
                          investor_name::String,
                          investor_dir::String,
                          scenario_names::Vector{String},
                          lag_bool::Bool) where T <: Union{PSY.Generator, PSY.Storage}
    name = get_name(device)

    base_power = PSY.get_base_power(device)
    size = get_device_size(device) * base_power

    products = create_products(simulation_data,
                               device,
                               projectdata)

    decision_year,
    construction_year,
    retirement_year,
    end_life_year,
    finance_data = create_investment_data(size, projectdata, investor_dir, simulation_data, products, investor_name, scenario_names, lag_bool)

    project = create_tech_type(name, device, projectdata, size, base_power, decision_year, construction_year, retirement_year, end_life_year, products, finance_data)

    return project
end


text:
function CR3BPRotatingFrame(arg0::jdouble, arg1::CelestialBody, arg2::CelestialBody)
    return CR3BPRotatingFrame((jdouble, CelestialBody, CelestialBody), arg0, arg1, arg2)
end

function get_ancestor(obj::Frame, arg0::jint)
    return jcall(obj, "getAncestor", Frame, (jint,), arg0)
end

function get_depth(obj::Frame)
    return jcall(obj, "getDepth", jint, ())
end

function get_frozen_frame(obj::Frame, arg0::Frame, arg1::AbsoluteDate, arg2::JString)
    return jcall(obj, "getFrozenFrame", Frame, (Frame, AbsoluteDate, JString), arg0, arg1, arg2)
end

function get_name(obj::Frame)
    return jcall(obj, "getName", JString, ())
end

function get_parent(obj::Frame)
    return jcall(obj, "getParent", Frame, ())
end

function get_root(::Type{Frame})
    return jcall(Frame, "getRoot", Frame, ())
end

function get_transform_provider(obj::Frame)
    return jcall(obj, "getTransformProvider", TransformProvider, ())
end

function get_transform_to(obj::Frame, arg0::Frame, arg1::AbsoluteDate)
    return jcall(obj, "getTransformTo", Transform, (Frame, AbsoluteDate), arg0, arg1)
end

function get_transform_to(obj::Frame, arg0::Frame, arg1::FieldAbsoluteDate)
    return jcall(obj, "getTransformTo", FieldTransform, (Frame, FieldAbsoluteDate), arg0, arg1)
end

function is_child_of(obj::Frame, arg0::Frame)
    return jcall(obj, "isChildOf", jboolean, (Frame,), arg0)
end

function is_pseudo_inertial(obj::Frame)
    return jcall(obj, "isPseudoInertial", jboolean, ())
end

function to_string(obj::Frame)
    return jcall(obj, "toString", JString, ())
end


text:
### A Pluto.jl notebook ###
# v0.11.14

using Markdown
using InteractiveUtils

# ╔═╡ 12478d9e-01f0-11eb-05d9-51c85d57ad10
begin
	using Pkg
	Pkg.activate(mktempdir())
end

# ╔═╡ 823ed6c4-01fb-11eb-0d84-3b3f9f659ec8
begin
	pkg"add StatsBase MixedModels GLM"

	using StatsBase       # Base package that has the "fit" verb
	using GLM             # defines a linear model type and how "fit" it
	using MixedModels     # defines model types and to "fit" them
end

# ╔═╡ 7cd6090e-0263-11eb-3c32-0fcb38a68d7d
md"""
# Using functions and types: philosophy

In this notebook we will look at the way that we think about using functions and types in Julia.

"""

# ╔═╡ d2a317f2-024f-11eb-19e3-739d3bbcd9e5
md"""

## The Julian point of view

Making an analogy between computer languages and natural (human) languages, we can think of **functions** as **verbs** that act on **objects**, which correspond to **nouns**.

> In Julia we do **verb-oriented programming**.

That is, we regard the verbs as the most important part of the program. We want to *do* stuff!

In contrast, in other languages like Python, C++ and Java it is usual to do object-oriented programming (OOP):

> In traditional object-oriented programming (OOP), we do **noun-oriented programming**.

In those languages we have nouns ("Classes"), which can or can't do certain functions (methods). The objects "contain" functionality within them, which turns out to often be very awkward

In Julia we have functions (verbs), which have many different implementations, i.e. **methods**, based on the numbers and types of the arguments passed to them. Stop and think about this world view for a minute!

We will learn more about types and how to create them later in this notebook.
But remembering that the **purpose of types is to extend functions with new methods** allows effective use of Julia.
"""

# ╔═╡ 6478b0ce-0250-11eb-2101-bd419c8876a3
md"""
# Functions: verbs that take nouns
"""

# ╔═╡ 010f58ee-0264-11eb-06eb-a92345fd47a4
md"""

Let's look at an example verb.

In statistics we talk about **fitting** a model to some data.

The package `StatsBase` defines the verb (function) `fit` and a few types of model.
Other packages import `StatsBase` and define more models and how the verb `fit` works on those models:
"""

# ╔═╡ 99475ac8-01fb-11eb-359c-7b61ffc1e67e
methods(StatsBase.fit)

# ╔═╡ 3897efb6-0263-11eb-1c53-93c71a2e65fe
md"""
There are **[298 registered Julia packages that depend on StatsBase](https://juliahub.com/ui/Packages/StatsBase/EZjIG/0.33.1?t=2)**!
Some of them use `fit`, or **extend** it to act on new types!
"""

# ╔═╡ 5541133a-0200-11eb-0ad0-65c155a1ae19
md"""
# Types: nouns and how to do verbs on them

One of the reasons for creating a new type in Julia is to add functionality to a verb (function) by extending it to the new type.

The basic way of creating new types is by defining `struct`s.

## Structs

Structs have three essential roles:

- They have a type -- this can be used to define how a function acts on the struct;
- They form a container which holds many values;
- They carry type information for the values contained inside.
"""

# ╔═╡ 93872f5c-026a-11eb-077e-7bcc945b21c8
md"
We have seen the basics of defining `struct`s in previous lectures. Now let's think more about how to extend verbs (functions) to act on them, and about the types of the objects they contain.
"

# ╔═╡ 12800e28-026b-11eb-174a-6b39f8024cf9
md"Let's define an object representing a point or vector in two dimensions with integer coordinates:"

# ╔═╡ 9ff41dbe-0200-11eb-2b17-f5e1be67c003
struct IntPoint
	x::Int64
	y::Int64
end

# ╔═╡ 3b2d5ffc-0254-11eb-2133-2d643e4f8380
md"""

- The **type** of the `struct` itself is `IntPoint`.
- It **contains** fields (data) `x` and `y`.
- The **type** of `x` and `y` are both `Int`.

In the same way we can define a similar object that represents a 2D vector whose coordinates are real numbers:
"""

# ╔═╡ aa2a4ea2-0200-11eb-009a-d7e3da674f03
struct FloatPoint
	x::Float64
	y::Float64
end

# ╔═╡ 79dc50b6-0254-11eb-1643-41e1ce885bd4
md"Now suppose that we want to be able to add two vectors.

We could invent a new function called `add`. But Julia already has a suitable function, namely... `+`!

The first thing we need to do is `import` the `+` function from Julia Base (where standard functions like `+` are defined), in order to be able to **extend** it by defining how it acts on our new types:
"

# ╔═╡ e7747972-0200-11eb-02f8-57282fbb5ecf
begin
	import Base: +

	+(a::IntPoint,   b::IntPoint)   = IntPoint(  a.x + b.x, a.y + b.y)

	+(a::FloatPoint, b::FloatPoint) = FloatPoint(a.x + b.x, a.y + b.y)
end

# ╔═╡ 7706a7e2-0254-11eb-1ceb-c71c4c3d99e2
md"The first thing to do is to check that it works!"

# ╔═╡ c8d3b47e-0200-11eb-1716-af05e82b8cee
begin
	a = IntPoint(1, 2)
	b = IntPoint(3, 4)
end

# ╔═╡ d9c83660-0200-11eb-0fb3-81b0d7f942a3
a + b

# ╔═╡ 93ed61c2-026b-11eb-3a9c-1d91b8726f69
md"""Let's see what's going on "under the hood":"""

# ╔═╡ ba1ed036-0254-11eb-1d06-cb20ef1eabd8
@code_lowered a + b

# ╔═╡ 3b9efcb6-0201-11eb-3fd4-2d28974d0a69
md"At the next stage in the compilation process, Julia **propagates the types it knows about**!:"

# ╔═╡ 3490369c-0201-11eb-02d7-216bac39ff41
@code_typed a + b # Julia propagates the types of the fields!

# ╔═╡ b99faabc-026b-11eb-1e4f-bb4aee759d20
md"What happens if we do the same with the floating-point version?"

# ╔═╡ 5f3b3eda-0258-11eb-21a3-2d1b2942a5ba
begin
	p = FloatPoint(.1, .2)
	q = FloatPoint(.3, .4)

	p + q
end

# ╔═╡ 922e343c-0258-11eb-2f7a-07857d4bf2a5
@code_typed p + q

# ╔═╡ 9d434312-0258-11eb-3bf0-b39d045a13a2
md"""
Note that Julia uses `add_int` in the case of `IntPoint`, but `add_float` in the `FloatPoint` case. When Julia knows the exact types, the run time is fast!

**Exercise**: Make an `UntypedPoint` as

```julia
struct UntypedPoint
	x
	y
end
```

and check what the typed generated code looks like!
"""

# ╔═╡ 1d9de34a-0262-11eb-3315-4544128264e1
md"""
Also, since we have extended `+` from Base, we can now automatically **use other functions that themselves use `+`**! This is very powerful. For example, the `sum` function is defined in terms of `+`, so the following works:
"""

# ╔═╡ 2fab9c78-0262-11eb-20ff-db4129f9ccaf
sum( FloatPoint(i / 10, (2i) / 10) for i in 1:10 )

# ╔═╡ 8fed54ae-0259-11eb-2d68-271a9098df38
md"""
## Parametric types
"""

# ╔═╡ 31bbe554-026c-11eb-12eb-e3f381c6fbf2
md"If we look at the way we defined `IntPoint` and `FloatPoint`, and then extended `+`, you should be feeling uncomfortable! We seem to be repeating ourselves a lot, since *the code looks exactly the same in each case* -- except that we change `Int64` to `Float64`. If we wanted to add, e.g., a `Rational` point, we would have to copy and paste yet again.

Computers are supposed to be good at automating things like this! So we should be able to get Julia to do this repetition for us!
"

# ╔═╡ 8b0b9780-026c-11eb-3f58-29d7be478750
md"The way to do so is by making the field type, i.e. the type of each variable contained inside an object of the type (`Int64` or `Float64` in the above examples) into a **type parameter**:
"

# ╔═╡ 8043cea0-0201-11eb-32b1-c1ce95a3232b
struct Point{T}
	x::T
	y::T
end

# ╔═╡ e7940ec4-026c-11eb-1160-4d79007f1c71
md"We use curly braces (`{` and `}`) to tell Julia that `T` is a type parameter.
When we now create objects of this type, Julia will *automatically make objects of the correct parametrised type!*:"

# ╔═╡ 24c8416a-0205-11eb-0037-bd088c0e65c5
Point(1, 2)

# ╔═╡ a1287c6a-0259-11eb-2051-bfa9aa1e3c84
Point(0.1, 0.2)

# ╔═╡ 3ff49068-025e-11eb-0f96-bd59dc20ffbc
Point(1//2, 3//4)

# ╔═╡ 11e8e816-026d-11eb-0ac5-ed4849ec66d0
md"We can extend functions that act on these types just like we did above. For example, let's extend the `angle` function to calculate the angle of the vector:"

# ╔═╡ aa9a36aa-0259-11eb-3c29-47966f1d707e
begin
	import Base: angle
	angle(p::Point) = atan(p.y, p.x)
end

# ╔═╡ c795d5a0-0259-11eb-1ef2-e5e6e76c602c
angle(Point(1, 2))

# ╔═╡ da41443e-0259-11eb-0304-6b92fb2de7d8
@code_typed angle(Point(1, 2))  # Change the arguments to Int64s and Rationals to see what happens!

# ╔═╡ 04219ea8-01ff-11eb-3019-b9bebbe14c43
md"""
## Abstract types

We now have different types `Point{T}` with different parameters `T`. What is `Point` itself?
"""

# ╔═╡ 9e8c92cc-026d-11eb-2b9e-9f76595e822f
Point

# ╔═╡ 9ac90468-026d-11eb-0a69-7194e6966932
md"""
**`Point`** is an **abstract type**. We can think of it as the combination, or **union** (in the sense of set theory), of *all* the types `Point{T}`, where `T` is replaced by all available types in the Julia environment, that either exist now or in the future when more code may be loaded from a package!

In general, abstract types in Julia are combinations, or a "merging", of other types into a single more general type. They are "abstract" in the sense that it is not possible to actually make objects of the abstract type. Objects always have a concrete type like `Int64` or `Point{Int64}`, which may be **subtypes** of an abstract type.

The **subtype** relation is written as `<:`.

"""

# ╔═╡ 88a7e668-0255-11eb-1d19-35e764da7c63
Point{Int} <: Point

# ╔═╡ 957250d6-0255-11eb-0193-d1da1ab2c640
Point <: Point{Int}

# ╔═╡ a7124698-0255-11eb-3ef5-137d00151773
md"You can create explicit unions with `Union`:"

# ╔═╡ 83465d62-0255-11eb-1a00-dfadf9c24d3a
Int <: Union{Int, Float64}

# ╔═╡ 20587260-026e-11eb-1299-39e2ca8ff416
md"""Here, `Union{Int64, Float64}` means the union (combination) of both of those types. We can read it as "either `Int64` or `Float64`".

We can use `isa` to check if a value is of a certain type:
"""

# ╔═╡ 34717d00-026e-11eb-010d-2f9ae44eb8d6
3 isa Int64

# ╔═╡ 44edf784-026e-11eb-25a9-7d852d047f2a
3 isa Union{Int, Float64}

# ╔═╡ 48de248c-026e-11eb-32ab-e73f454aaeca
3.5 isa Int64

# ╔═╡ 4b274aae-026e-11eb-0764-71f6879c8a0e
3.5 isa Union{Int64, Float64}

# ╔═╡ a1648f62-0255-11eb-1ff9-1738a037a86c
Float64 <: Union{Int, Float64}

# ╔═╡ b17fbfe8-0255-11eb-2c58-bbed6fe2cc7f
md"""
## Different kinds of `Point`s

The mathematical concept of **point** is more abstract than just $x$ and $y$ coordinates. For example, we could represent a point using polar coordinates instead.

But we could still ask the same questions, which correspond to **verbs**:

- `x_coord`: what is the $x$ coordinate?
- `y_coord`: what is the $y$ coordinate?
- `radius`: What is the distance from the origin?
- `angle`: What is the angle from the $x$ axis?
"""

# ╔═╡ 33a3695c-026f-11eb-2481-0bc92174ee8e
md"Let's define types that correspond to the two representations of a point:"

# ╔═╡ e6a13f4e-0255-11eb-0748-dba33f2bd1be
struct CartesianPoint{T}
	x::T
	y::T
end

# ╔═╡ 46ee84ee-026f-11eb-3a00-959ae43bd920
struct PolarPoint{T}
	r::T
	θ::T
end

# ╔═╡ 05b2bbec-0256-11eb-3d92-bd5956c64cf8
begin
	x_coord(p::CartesianPoint) = p.x
	y_coord(p::CartesianPoint) = p.y

	x_coord(p::PolarPoint) = p.r * cos(p.θ)
	y_coord(p::PolarPoint) = p.r * sin(p.θ)

	theta(p) = atan(y_coord(p), x_coord(p))  # abstract
	theta(p::PolarPoint) = p.θ # -- optimization

	radius(p) = hypot(xcoord(p), ycoord(p)) # abstract
	# radius(p::PolarPoint) = p.r -- optimization
end

# ╔═╡ 85904c9e-026f-11eb-312b-79d0652db9ca
md"We have defined **generic** methods that act on any object `p`. They will work as long as the functions that we call also work. For example:"

# ╔═╡ 9b2e0054-025c-11eb-29f1-07e6a2f9f6c6
theta(CartesianPoint(1,2))

# ╔═╡ bd1d3278-025c-11eb-344e-ed13b50353c8
theta(PolarPoint(2.0,0.4636476090008061))

# ╔═╡ 05849030-025d-11eb-3c7e-07723387e94f
@code_typed theta(PolarPoint(2.0,0.4636476090008061)) # But we already know theta!

# ╔═╡ 0b044d6e-0270-11eb-20e6-9165669bb69f
md"""
But in fact, for a `PolarPoint` we already know theta -- it is crazy to do that complicated calculation. So we can define a **specialised method** that works *only* for `PolarPoint`, and leave the **generic "fallback"** that works with other types!

Julia will check the type of the argument we pass in to the `theta` function and will choose the **most specific** method of the `theta` function that applies to objects of that type.
"""

# ╔═╡ a1dac51c-025d-11eb-10ff-f7e813f9832a
md"""
Note that currently the functions `theta` and `radius` accept any type at all! But this is probably *not* what we actually intended, for example:
"""

# ╔═╡ 8f54be20-025d-11eb-2bcf-dd95f61eb477
theta("hello")

# ╔═╡ d34342a6-0270-11eb-0ea8-b141ccf5ae50
md"The function `theta` itself is happy to accept a `String`, since we didn't restrict its input arguments at all. However, the first thing that it does is call `y_coord`, which fails since it has no method which can accept a `String`.
"

# ╔═╡ 142c6db0-025d-11eb-16ff-dd0bed953f8b
md"""
A solution to this would be to create an abstract type which encompasses the different `Point` types, and define `len` and `theta` only on that.

We could use a `Union{CartesianPoint, PolarPoint}`, but we actually want to **allow future users to create new `Point` types** and have this code still work, provided they define `x_coord` and `y_coord` on the new type!

To do this we *make a new abstract type*:
"""

# ╔═╡ 83445c12-025d-11eb-06fa-91b2906b95c9
abstract type AbstractPoint end

# ╔═╡ cae32b5c-025d-11eb-365e-f1c87a3da5fc
md"""Now we can add `<: AbstractPoint` to the struct definitions to make them subtypes of this new `AbstractPoint` abstract type, and we can annotate the function definitions to only accept objects which are `AbstractPoint`s!
"""

# ╔═╡ a309eacc-0278-11eb-3cca-17a714e0d512
struct NewPoint <: AbstractPoint
	x
	y
end

# ╔═╡ b6add3ec-0278-11eb-35be-f58dc34f4d87
new_theta(p::AbstractPoint) = theta(p)

# ╔═╡ 4083d902-0271-11eb-1182-a9fc930470df
md"As a complicated example of a **type hierarchy** in Julia, here is the hierarchy of number types that are defined in base Julia:
"

# ╔═╡ 14caccfc-025e-11eb-1244-67d123320e10
md"![](https://upload.wikimedia.org/wikipedia/commons/thumb/4/40/Type-hierarchy-for-julia-numbers.png/1024px-Type-hierarchy-for-julia-numbers.png)"

# ╔═╡ 5a0b6818-0271-11eb-1475-b3e814606f77
md"The nodes with no children are concrete types. All the other nodes are abstract types. We may define new number types by making them subtypes of the relevant abstract type. For example, the [DoubleFloats.jl](https://github.com/JuliaMath/DoubleFloats.jl) package defines a new type `Double64`, having a precision roughly twice that of a `Float64`, as a new subtype of `AbstractFloat`.
"

# ╔═╡ Cell order:
# ╠═12478d9e-01f0-11eb-05d9-51c85d57ad10
# ╟─7cd6090e-0263-11eb-3c32-0fcb38a68d7d
# ╟─d2a317f2-024f-11eb-19e3-739d3bbcd9e5
# ╟─6478b0ce-0250-11eb-2101-bd419c8876a3
# ╟─010f58ee-0264-11eb-06eb-a92345fd47a4
# ╠═823ed6c4-01fb-11eb-0d84-3b3f9f659ec8
# ╠═99475ac8-01fb-11eb-359c-7b61ffc1e67e
# ╟─3897efb6-0263-11eb-1c53-93c71a2e65fe
# ╟─5541133a-0200-11eb-0ad0-65c155a1ae19
# ╟─93872f5c-026a-11eb-077e-7bcc945b21c8
# ╟─12800e28-026b-11eb-174a-6b39f8024cf9
# ╠═9ff41dbe-0200-11eb-2b17-f5e1be67c003
# ╟─3b2d5ffc-0254-11eb-2133-2d643e4f8380
# ╠═aa2a4ea2-0200-11eb-009a-d7e3da674f03
# ╟─79dc50b6-0254-11eb-1643-41e1ce885bd4
# ╠═e7747972-0200-11eb-02f8-57282fbb5ecf
# ╟─7706a7e2-0254-11eb-1ceb-c71c4c3d99e2
# ╠═c8d3b47e-0200-11eb-1716-af05e82b8cee
# ╠═d9c83660-0200-11eb-0fb3-81b0d7f942a3
# ╟─93ed61c2-026b-11eb-3a9c-1d91b8726f69
# ╠═ba1ed036-0254-11eb-1d06-cb20ef1eabd8
# ╟─3b9efcb6-0201-11eb-3fd4-2d28974d0a69
# ╠═3490369c-0201-11eb-02d7-216bac39ff41
# ╟─b99faabc-026b-11eb-1e4f-bb4aee759d20
# ╠═5f3b3eda-0258-11eb-21a3-2d1b2942a5ba
# ╠═922e343c-0258-11eb-2f7a-07857d4bf2a5
# ╟─9d434312-0258-11eb-3bf0-b39d045a13a2
# ╟─1d9de34a-0262-11eb-3315-4544128264e1
# ╠═2fab9c78-0262-11eb-20ff-db4129f9ccaf
# ╟─8fed54ae-0259-11eb-2d68-271a9098df38
# ╟─31bbe554-026c-11eb-12eb-e3f381c6fbf2
# ╟─8b0b9780-026c-11eb-3f58-29d7be478750
# ╠═8043cea0-0201-11eb-32b1-c1ce95a3232b
# ╟─e7940ec4-026c-11eb-1160-4d79007f1c71
# ╠═24c8416a-0205-11eb-0037-bd088c0e65c5
# ╠═a1287c6a-0259-11eb-2051-bfa9aa1e3c84
# ╠═3ff49068-025e-11eb-0f96-bd59dc20ffbc
# ╟─11e8e816-026d-11eb-0ac5-ed4849ec66d0
# ╠═aa9a36aa-0259-11eb-3c29-47966f1d707e
# ╠═c795d5a0-0259-11eb-1ef2-e5e6e76c602c
# ╠═da41443e-0259-11eb-0304-6b92fb2de7d8
# ╟─04219ea8-01ff-11eb-3019-b9bebbe14c43
# ╠═9e8c92cc-026d-11eb-2b9e-9f76595e822f
# ╟─9ac90468-026d-11eb-0a69-7194e6966932
# ╠═88a7e668-0255-11eb-1d19-35e764da7c63
# ╠═957250d6-0255-11eb-0193-d1da1ab2c640
# ╟─a7124698-0255-11eb-3ef5-137d00151773
# ╠═83465d62-0255-11eb-1a00-dfadf9c24d3a
# ╟─20587260-026e-11eb-1299-39e2ca8ff416
# ╠═34717d00-026e-11eb-010d-2f9ae44eb8d6
# ╠═44edf784-026e-11eb-25a9-7d852d047f2a
# ╠═48de248c-026e-11eb-32ab-e73f454aaeca
# ╠═4b274aae-026e-11eb-0764-71f6879c8a0e
# ╠═a1648f62-0255-11eb-1ff9-1738a037a86c
# ╟─b17fbfe8-0255-11eb-2c58-bbed6fe2cc7f
# ╟─33a3695c-026f-11eb-2481-0bc92174ee8e
# ╠═e6a13f4e-0255-11eb-0748-dba33f2bd1be
# ╠═46ee84ee-026f-11eb-3a00-959ae43bd920
# ╠═05b2bbec-0256-11eb-3d92-bd5956c64cf8
# ╟─85904c9e-026f-11eb-312b-79d0652db9ca
# ╠═9b2e0054-025c-11eb-29f1-07e6a2f9f6c6
# ╠═bd1d3278-025c-11eb-344e-ed13b50353c8
# ╠═05849030-025d-11eb-3c7e-07723387e94f
# ╟─0b044d6e-0270-11eb-20e6-9165669bb69f
# ╟─a1dac51c-025d-11eb-10ff-f7e813f9832a
# ╠═8f54be20-025d-11eb-2bcf-dd95f61eb477
# ╟─d34342a6-0270-11eb-0ea8-b141ccf5ae50
# ╟─142c6db0-025d-11eb-16ff-dd0bed953f8b
# ╠═83445c12-025d-11eb-06fa-91b2906b95c9
# ╟─cae32b5c-025d-11eb-365e-f1c87a3da5fc
# ╠═a309eacc-0278-11eb-3cca-17a714e0d512
# ╠═b6add3ec-0278-11eb-35be-f58dc34f4d87
# ╟─4083d902-0271-11eb-1182-a9fc930470df
# ╟─14caccfc-025e-11eb-1244-67d123320e10
# ╟─5a0b6818-0271-11eb-1475-b3e814606f77

max_stars_repo_name: mathijsvdv/ComputationalThinking

max_issues_repo_name: mathijsvdv/ComputationalThinking

max_forks_repo_name: mathijsvdv/ComputationalThinking


text:
"""
    WoodburyPDMat(
        A::AbstractMatrix{T}, D::Diagonal{T}, S::Diagonal{T},
    ) where {T<:Real}

Lazily represents matrices of the form
```julia
W = A * D * A' + S
```
`D` and `S` must have only non-negative entries.

Using this matrix type is a good idea if `size(A, 1) > size(A, 2)` as the structure in the
matrix can be exploited to accelerate operations involving `W`'s inverse [1], such as
`invquad`, and it's determinant [2], such as `logdet`.

You probably don't want to use this matrix type if `size(A, 1) < size(A, 2)`.

[1] - https://en.wikipedia.org/wiki/Woodbury_matrix_identity
[2] - https://en.wikipedia.org/wiki/Matrix_determinant_lemma
"""
struct WoodburyPDMat{
    T<:Real, TA<:AbstractMatrix{T}, TD<:Diagonal{T}, TS<:Diagonal{T},
} <: AbstractPDMat{T}
    A::TA
    D::TD
    S::TS
    function WoodburyPDMat(
        A::AbstractMatrix{T}, D::Diagonal{T}, S::Diagonal{T},
    ) where {T<:Real}
        validate_woodbury_arguments(A, D, S)
        return new{T, typeof(A), typeof(D), typeof(S)}(A, D, S)
    end
end

function WoodburyPDMat(
    A::AbstractMatrix{T}, D::AbstractVector{T}, S::AbstractVector{T},
) where {T<:Real}
    return WoodburyPDMat(A, Diagonal(D), Diagonal(S))
end

WoodburyPDMat{T, TA, TD, TS}(vals...) where {T, TA, TD, TS} = WoodburyPDMat(vals...)

PDMats.dim(W::WoodburyPDMat) = size(W.A, 1)

# Convesion method. Primarily useful for testing purposes.
Base.Matrix(W::WoodburyPDMat) = W.A * W.D * W.A' + W.S

Base.getindex(W::WoodburyPDMat, inds...) = getindex(Matrix(W), inds...)

function validate_woodbury_arguments(A, D, S)
    if size(A, 1) != size(S, 1)
        throw(ArgumentError("size(A, 1) != size(S, 1)"))
    end
    if size(A, 2) != size(D, 1)
        throw(ArgumentError("size(A, 2) != size(D, 1)"))
    end
    if any(x -> x < 0, diag(D))
        throw(ArgumentError("Detected negative element on diagonal of D: $(D)"))
    end
    if any(x -> x < 0, diag(S))
        throw(ArgumentError("Detected negative element on diagonal of S: $(S)"))
    end
end

@non_differentiable validate_woodbury_arguments(A, D, S)

function LinearAlgebra.logdet(W::WoodburyPDMat)
    C_S = cholesky(W.S)
    B = C_S.U' \ (W.A * cholesky(W.D).U')
    return logdet(C_S) + logdet(cholesky(Symmetric(I + B'B)))
end

# Utilises the matrix inversion lemma to produce an efficient implementation.
function PDMats.invquad(W::WoodburyPDMat{<:Real}, x::AbstractVector{<:Real})
    C_S = cholesky(W.S)
    B = C_S.U' \ (W.A * cholesky(W.D).U')
    α = C_S.U' \ x
    β = B' * α
    return α'α - sum(abs2, cholesky(Symmetric(I + B'B)).U' \ β)
end

# This doesn't get us the computational wins, but it's unclear how to construct a
# root for a Woodbury matrix. Consequently, if performance is very important when sampling,
# it's necessary to implement a method of `rand` or `_rand` that explicitly uses ancestral
# sampling to exploit the approximately low-rank structre in a WoodburyPDMat.
function PDMats.unwhiten!(r::DenseVecOrMat, W::WoodburyPDMat{<:Real}, x::DenseVecOrMat)
    return unwhiten!(r, PDMat(Symmetric(Matrix(W))), x)
end

# NOTE: the parameterisation to scale up the Woodbury matrix is not unique. Here we
# implement one way to scale it.
*(a::WoodburyPDMat, c::Real) = WoodburyPDMat(a.A, a.D * c, a.S * c)
*(c::Real, a::WoodburyPDMat) = a * c