text:
"""
    cla()

Clears the current axes, removing all drawings and resetting all options.
"""
function cla()
    f = gcf();
    isempty(f.axes) && return nothing
    reset!(gca(); parent=f.id);
    PreviewFigure(f)
end

"""
    clo()

Clears all objects (annotations, arrows, ...) from the current axes, leaves everything else.
"""
clo() = (gca().objects = Vector{Object2D}(); PreviewFigure(gcf()))

####
#### axes
####

function axis(short::String=""; axes=nothing)
    axes = check_axes(axes)
    if !isempty(short)
        s_lc = lowercase(short)
        if s_lc == "math"
            axes.math = true
        elseif s_lc == "nomath"
            axes.math = false
        elseif s_lc == "equal"
            throw(NotImplementedError("axis(equal)"))
            #= would need to
            1. get the size from gca() to get the aspect ratio
            2. modify the xlim/ylim so that they match that ratio
            the difficulty is to get xlim/ylim when they're not set
            explicitly. would need to figure out a way
            to declare local variables
            Maybe just xaxis min xgmin*... max xgmax*...
            =#
        else
            throw(OptionValueError("Unrecognised shorthand toggle for axis.", short))
        end
    end
    return PreviewFigure(gcf())
end

####
#### [x|y|...]axis(...)
####

for axs ∈ ("x", "y", "x2", "y2")
    f  = Symbol(axs * "axis")    # :xaxis  ...
    ex = quote
        function $f(short::String=""; axes=nothing, opts...)
            axes = check_axes(axes)
            if isempty(short)
                set_properties!(axes.$f; opts...)
            else
                # try to interpret the string as a shorthand and apply on gca
                s_lc = lowercase(short)
                if s_lc == "off"
                    set_properties!(gca().$f; off=true, opts...)
                elseif s_lc == "on"
                    set_properties!(gca().$f; off=false, opts...)
                elseif s_lc == "log"
                    set_properties!(gca().$f; log=true, opts...)
                elseif s_lc ∈ ["lin", "linear"]
                    set_properties!(gca().$f; log=false, opts...)
                else
                    throw(OptionValueError("Unrecognised shorthand toggle for axis.", short))
                end
            end
            return PreviewFigure(gcf())
        end
    end
    eval(ex)
end

"""
    math(s;...)

Toggle the (current) axes to math mode (where the axes go through (0,0)).
to reverse this (`b=false`) It is recommended to also adjust the axis limits via [`xlim`](@ref)
and [`ylim`](@ref) to make sure that the origin is somewhere in the drawn area (otherwise the
results will be rather ugly).

### Examples

```julia
math()
math("on")
math("off")
```
"""
function math(short::String=""; axes=nothing)
    axes = check_axes(axes)
    s_lc = lowercase(short)
    if isempty(short) || s_lc == "on"
         axes.math = true
    elseif s_lc == "off"
        axes.math = false
    else
        throw(OptionValueError("Unrecognised shorthand toggle for math.", short))
    end
    return PreviewFigure(gcf())
end

"""
    grid()

Set grid mode on. By default the grid will be associated with both `xticks` and `yticks` but you
can also specify one axis to only have horizontal or vertical lines by passing `axis=["x"]`.
Options can be passed to specify the color of the grid, the style of the lines or their width.

### Examples

```julia
grid(linestyle="-.")
grid(axis=["y"], linestyle="--", color="lightgray")
```
"""
function grid(short::String=""; axes=nothing, axis::Vector{String}=["x", "y"],
               opts...)
    axes = check_axes(axes)
    if !isempty(short)
        s_lc = lowercase(short)
        if s_lc == "off"
            axes.xaxis.ticks.grid = false
            axes.yaxis.ticks.grid = false
            return PreviewFigure(gcf())
        elseif s_lc == "on"
            axes.xaxis.ticks.grid = true
            axes.yaxis.ticks.grid = true
            return PreviewFigure(gcf())
        else
            throw(OptionValueError("Unrecognised shorthand toggle for grid.", short))
        end
    end
    for ax ∈ axis
        ax_lc = lowercase(ax)
        ax_lc ∈ ("x","x2","y","y2") || throw(OptionValueError("Unrecognized axis symbol(s).", axis))
        # all options affect the ticks (in particular: color and width)
        axsym = Symbol(ax_lc * "axis")
        ticks = eval(:($axes.$axsym.ticks))
        ticks.grid = true
        for optname ∈ opts.itr
            if optname ∈ (:col, :color)
                set_color!(ticks, col(opts[optname], optname))
            elseif optname ∈ (:ls, :lstyle, :linestyle)
                set_lstyle!(ticks, opts[optname])
            elseif optname ∈ (:lw, :lwidth, :linewidth)
                set_lwidth!(ticks, posfl(opts[optname], optname))
            else
                throw(UnknownOptionError(optname, ax))
            end
        end
    end
    return PreviewFigure(gcf())
end

####
#### [x|y]lim, [x|y]lim! (synonyms though with ! is preferred)
####

function _lim(axis_sym::Symbol, min::Option{Float64}, max::Option{Float64}; axes=nothing)
    axes = check_axes(axes)
    if min !== nothing && max !== nothing
        min < max || throw(ArgumentError("min must be smaller than max (got ($min, $max))"))
    end
    axis = getfield(axes, axis_sym)
    setfield!(axis, :min, min)
    setfield!(axis, :max, max)
    return PreviewFigure(gcf())
end

# Generate *lim each axis
for axs ∈ ("x", "y", "x2", "y2")
    f  = Symbol(axs * "lim")
    ex = quote
        $f(min, max; o...) = _lim(Symbol($axs * "axis"), fl(min), fl(max); o...)
        $f(; min=∅, max=∅) = _lim(Symbol($axs * "axis"), fl(min), fl(max))
    end
    eval(ex)
end

####
#### [x|y]scale
####

function _scale(axis_sym::Symbol, v::String; axes=nothing)
    axes = check_axes(axes)
    axis = getfield(axes, axis_sym)
    axis.log = get(AXSCALE, v) do
        throw(OptionValueError("axis scale", v))
    end
    return PreviewFigure(gcf())
end

# Generate *scale for each axis
for axs ∈ ("x", "y", "x2", "y2")
    f  = Symbol(axs * "scale")
    eval(:($f(v::String; o...) = _scale(Symbol($axs * "axis"), v, o...)))
end


text:
export Codecov

"""
Coverage.Codecov Module

This module provides functionality to push coverage information to the CodeCov.io
web service. It exports the `submit` and `submit_local` methods.
"""
module Codecov

using HTTP
using Coverage
using CoverageTools
using JSON
using LibGit2

export submit, submit_local, submit_generic

#=
JSON structure for Codecov.io
https://codecov.io/api#post-report

{
  "coverage": {
    "path/to/file.py": [null, 1, 0, null, true, 0, 0, 1, 1],
    "path/to/other.py": [null, 0, 1, 1, "1/3", null]
  },
  "messages": {
    "path/to/other.py": {
      "1": "custom message for line 1"
    }
  }
}
=#

# Turn vector of filename : coverage pairs into a dictionary
function to_json(fcs::Vector{FileCoverage})
    cov = Dict()
    for fc in fcs
        cov[fc.filename] = vcat(nothing, fc.coverage)
    end
    return Dict("coverage" => cov)
end

"""
kwargs provides default values to insert into args_array, only if they are
not already specified in args_array.
"""
function set_defaults(args_array::Dict; kwargs...)
    args_array = copy(args_array)
    for kwarg in kwargs
        if !haskey(args_array, kwarg[1])
            push!(args_array, kwarg)
        end
    end
    return args_array
end

"""
    submit(fcs::Vector{FileCoverage})

Takes a vector of file coverage results (produced by `process_folder`),
and submits them to Codecov.io. Assumes that this code is being run
on TravisCI or AppVeyor. If running locally, use `submit_local`.
"""
function submit(fcs::Vector{FileCoverage}; kwargs...)
    submit_generic(fcs, add_ci_to_kwargs(; kwargs...))
end

add_ci_to_kwargs(; kwargs...) = add_ci_to_kwargs(Dict{Symbol,Any}(kwargs))
function add_ci_to_kwargs(kwargs::Dict)
    if lowercase(get(ENV, "APPVEYOR", "false")) == "true"
        appveyor_pr = get(ENV, "APPVEYOR_PULL_REQUEST_NUMBER", "")
        appveyor_job = join(
            [
                ENV["APPVEYOR_ACCOUNT_NAME"],
                ENV["APPVEYOR_PROJECT_SLUG"],
                ENV["APPVEYOR_BUILD_VERSION"],
            ],
            "%2F",
        )
        kwargs = set_defaults(kwargs,
            service      = "appveyor",
            branch       = ENV["APPVEYOR_REPO_BRANCH"],
            commit       = ENV["APPVEYOR_REPO_COMMIT"],
            pull_request = appveyor_pr,
            job          = appveyor_job,
            slug         = ENV["APPVEYOR_REPO_NAME"],
            build        = ENV["APPVEYOR_JOB_ID"],
        )
    elseif lowercase(get(ENV, "TRAVIS", "false")) == "true"
        kwargs = set_defaults(kwargs,
            service      = "travis-org",
            branch       = ENV["TRAVIS_BRANCH"],
            commit       = ENV["TRAVIS_COMMIT"],
            pull_request = ENV["TRAVIS_PULL_REQUEST"],
            job          = ENV["TRAVIS_JOB_ID"],
            slug         = ENV["TRAVIS_REPO_SLUG"],
            build        = ENV["TRAVIS_JOB_NUMBER"],
        )
    elseif lowercase(get(ENV, "CIRCLECI", "false")) == "true"
        circle_slug = join(
            [
                ENV["CIRCLE_PROJECT_USERNAME"],
                ENV["CIRCLE_PROJECT_REPONAME"],
            ],
            "%2F",
        )
        kwargs = set_defaults(kwargs,
            service      = "circleci",
            branch       = ENV["CIRCLE_BRANCH"],
            commit       = ENV["CIRCLE_SHA1"],
            pull_request = get(ENV, "CIRCLE_PR_NUMBER", "false"),  # like Travis
            build_url    = ENV["CIRCLE_BUILD_URL"],
            slug         = circle_slug,
            build        = ENV["CIRCLE_BUILD_NUM"],
        )
    elseif lowercase(get(ENV, "JENKINS", "false")) == "true"
        kwargs = set_defaults(kwargs,
            service      = "jenkins",
            branch       = ENV["GIT_BRANCH"],
            commit       = ENV["GIT_COMMIT"],
            job          = ENV["JOB_NAME"],
            build        = ENV["BUILD_ID"],
            build_url    = ENV["BUILD_URL"],
            jenkins_url  = ENV["JENKINS_URL"],
        )
    elseif haskey(ENV, "BUILD_BUILDURI") # Azure Pipelines
        ref = get(ENV, "SYSTEM_PULLREQUEST_TARGETBRANCH", ENV["BUILD_SOURCEBRANCHNAME"])
        branch = startswith(ref, "refs/heads/") ? ref[12:end] : ref
        kwargs = set_defaults(kwargs,
            service      = "azure_pipelines",
            branch       = branch,
            commit       = ENV["BUILD_SOURCEVERSION"],
            pull_request = get(ENV, "SYSTEM_PULLREQUEST_PULLREQUESTNUMBER", ""),
            job          = ENV["BUILD_DEFINITIONNAME"],
            slug         = ENV["BUILD_REPOSITORY_NAME"],
            build        = ENV["BUILD_BUILDID"],
        )
    elseif haskey(ENV, "GITHUB_ACTION") # GitHub Actions
        event_path = open(JSON.Parser.parse, ENV["GITHUB_EVENT_PATH"])
        ref = ENV["GITHUB_REF"]
        if startswith(ref, "refs/heads/")
            branch = ref[12:end]
            ga_pr = "false"
        elseif startswith(ref, "refs/tags/")
            branch = ref[11:end]
            ga_pr = "false"
        elseif startswith(ref, "refs/pull/")
            branch = ENV["GITHUB_HEAD_REF"]
            ga_pr_info = get(event_path, "pull_request", Dict())
            ga_pr = get(ga_pr_info, "number", "false")
        end
        ga_build_url = "https://github.com/$(ENV["GITHUB_REPOSITORY"])/actions/runs/$(ENV["GITHUB_RUN_ID"])"
        kwargs = set_defaults(kwargs,
            service      = "github-actions",
            branch       = branch,
            commit       = ENV["GITHUB_SHA"],
            pull_request = ga_pr,
            slug         = ENV["GITHUB_REPOSITORY"],
            build        = ENV["GITHUB_RUN_ID"],
            build_url    = ga_build_url,
        )
    elseif lowercase(get(ENV, "BUILDKITE", "false")) == "true"
        kwargs = set_defaults(kwargs,
            service      = "buildkite",
            branch       = ENV["BUILDKITE_BRANCH"],
            commit       = ENV["BUILDKITE_COMMIT"],
            job          = ENV["BUILDKITE_JOB_ID"],
            build        = ENV["BUILDKITE_BUILD_NUMBER"],
            build_url    = ENV["BUILDKITE_BUILD_URL"]
        )
        if ENV["BUILDKITE_PULL_REQUEST"] != "false"
            kwargs = set_defaults(kwargs, pr = ENV["BUILDKITE_PULL_REQUEST"])
        end
    else
        error("No compatible CI platform detected")
    end

    return kwargs
end

"""
    submit_local(fcs::Vector{FileCoverage}, dir::AbstractString=pwd())

Take a `Vector` of file coverage results (produced by `process_folder`),
and submit them to Codecov.io. Assumes the submission is being made from
a local git installation, rooted at `dir`. A repository token should be specified by a
`token` keyword argument or the `CODECOV_TOKEN` environment variable.
"""
function submit_local(fcs::Vector{FileCoverage}, dir::AbstractString=pwd(); kwargs...)
    submit_generic(fcs, add_local_to_kwargs(dir; kwargs...))
end

add_local_to_kwargs(dir; kwargs...) = add_local_to_kwargs(dir, Dict{Symbol,Any}(kwargs))
function add_local_to_kwargs(dir, kwargs::Dict)
    LibGit2.with(LibGit2.GitRepoExt(dir)) do repo
        LibGit2.with(LibGit2.head(repo)) do headref
            branch_name = LibGit2.shortname(headref) # this function returns a String
            commit_oid  = LibGit2.GitHash(LibGit2.peel(headref))
            kwargs = set_defaults(kwargs,
                commit = string(commit_oid),
                branch = branch_name
                )
        end
    end

    return kwargs
end

"""
    submit_generic(fcs::Vector{FileCoverage})

Takes a vector of file coverage results (produced by `process_folder`),
and submits them to a Codecov.io instance. Keyword arguments are converted
into a generic Codecov.io API uri.  It is essential that the keywords and
values match the Codecov upload/v2 API specification.
The `codecov_url` keyword argument or the CODECOV_URL environment variable
can be used to specify the base path of the uri.
The `codecov_url_path` keyword argument or the CODECOV_URL_PATH environment variable
can be used to specify the final path of the uri.
The `dry_run` keyword can be used to prevent the http request from
being generated.
"""
submit_generic(fcs::Vector{FileCoverage}; kwargs...) =
    submit_generic(fcs, Dict{Symbol,Any}(kwargs))
function submit_generic(fcs::Vector{FileCoverage}, kwargs::Dict)
    @assert length(kwargs) > 0
    dry_run = get(kwargs, :dry_run, false)

    uri_str = construct_uri_string(kwargs)

    @info "Submitting data to Codecov..."
    @debug "Codecov.io API URL:\n" * mask_token(uri_str)

    is_black_hole_server = parse(Bool, strip(get(ENV, "JULIA_COVERAGE_IS_BLACK_HOLE_SERVER", "false")))::Bool
    if !dry_run
        # Tell Codecov we have an upload for them
        response = HTTP.post(uri_str; headers=Dict("Accept" => "text/plain"))
        # Get the temporary URL to use for uploading to S3
        repr = String(response)
        s3url = get(split(String(response.body), '\n'), 2, "")
        repr = chomp(replace(repr, s3url => ""))
        @debug "Result of submission:" * repr
        !is_black_hole_server && upload_to_s3(; s3url=s3url, fcs=fcs)
    end
end

function upload_to_s3(; s3url, fcs)
    startswith(s3url, "https://") || error("Invalid codecov response: $s3url")
    # Upload to S3
    request = HTTP.put(s3url; body=json(to_json(fcs)),
                       header=Dict("Content-Type" => "application/json",
                                   "x-amz-storage-class" => "REDUCED_REDUNDANCY"))
    @debug "Result of submission:" * mask_token(String(request))
end

function construct_uri_string(kwargs::Dict)
    url = get(ENV, "CODECOV_URL", "")
    isempty(url) || (kwargs = set_defaults(kwargs, codecov_url = url))

    path = get(ENV, "CODECOV_URL_PATH", "")
    isempty(path) || (kwargs = set_defaults(kwargs, codecov_url_path = path))

    token = get(ENV, "CODECOV_TOKEN", "")
    isempty(token) || (kwargs = set_defaults(kwargs, token = token))

    flags = get(ENV, "CODECOV_FLAGS", "")
    isempty(flags) || (kwargs = set_defaults(kwargs; flags = flags))

    name = get(ENV, "CODECOV_NAME", "")
    isempty(name) || (kwargs = set_defaults(kwargs; name = name))

    codecov_url = get(kwargs, :codecov_url, "https://codecov.io")
    if isempty(codecov_url) || codecov_url[end] == '/'
        error("the codecov_url should not end with a /, given url $(repr(codecov_url))")
    end

    codecov_url_path = get(kwargs, :codecov_url_path, "/upload/v4")
    if isempty(codecov_url_path) || codecov_url_path[1] != '/' || codecov_url_path[end] == '/'
        error("the codecov_url_path should begin with, but not end with, a /, given url $(repr(codecov_url_path))")
    end

    uri_str = "$(codecov_url)$(codecov_url_path)?"
    for (k, v) in kwargs
        # add all except a few special key/value pairs to the URL
        # (:verbose is there for backwards compatibility with versions
        # of this code that treated it in a special way)
        if k != :codecov_url && k != :dry_run && k != :verbose
            uri_str = "$(uri_str)$(k)=$(v)&"
        end
    end

    return uri_str
end

function mask_token(uri_string)
    return replace(uri_string, r"token=[^&]*" => "token=<HIDDEN>")
end

end # module

ext: jl


text:
# read the output current
export psp_readcurrent

const READ_CURRENT = [0x41, 0x0d]

"""
    psp_readcurrent(io)

Read the output current.
"""
function psp_readcurrent(io_psp::IO)
  write(io_psp, READ_CURRENT)
  currentstring = bytestring(read(io_psp,UInt8, 9))
  if  currentstring[1] != 'A' ||
      currentstring[3] != '.' ||
      currentstring[7:9] != "\r\r\n"
        error("PSP Invalid Response: ",currentstring)
  end
  parse(Float64,currentstring[2:6])
end


text:
module Sensorfusion

export Statevector, KalmanFilter

include("statevec.jl")
include("kalman.jl")

end


text:
module TestCausalAlgorithms
using Base.Test
using Causality
using DataFrames
using Requests

sample_data = replace(replace(Requests.get("http://data.princeton.edu/wws509/datasets/smoking.dat").data, r" +", "\t"), r"\r", "")
io = IOBuffer(sample_data)
df = DataFrames.readtable(io, separator='\t', makefactors=true, header=true)
vars = names(df)

# each of these


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiEventsV1beta1EventSeries <: SwaggerModel
    count::Any # spec type: Union{ Nothing, Int32 } # spec name: count
    lastObservedTime::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1MicroTime } # spec name: lastObservedTime
    state::Any # spec type: Union{ Nothing, String } # spec name: state

    function IoK8sApiEventsV1beta1EventSeries(;count=nothing, lastObservedTime=nothing, state=nothing)
        o = new()
        validate_property(IoK8sApiEventsV1beta1EventSeries, Symbol("count"), count)
        setfield!(o, Symbol("count"), count)
        validate_property(IoK8sApiEventsV1beta1EventSeries, Symbol("lastObservedTime"), lastObservedTime)
        setfield!(o, Symbol("lastObservedTime"), lastObservedTime)
        validate_property(IoK8sApiEventsV1beta1EventSeries, Symbol("state"), state)
        setfield!(o, Symbol("state"), state)
        o
    end
end # type IoK8sApiEventsV1beta1EventSeries

const _property_map_IoK8sApiEventsV1beta1EventSeries = Dict{Symbol,Symbol}(Symbol("count")=>Symbol("count"), Symbol("lastObservedTime")=>Symbol("lastObservedTime"), Symbol("state")=>Symbol("state"))
const _property_types_IoK8sApiEventsV1beta1EventSeries = Dict{Symbol,String}(Symbol("count")=>"Int32", Symbol("lastObservedTime")=>"IoK8sApimachineryPkgApisMetaV1MicroTime", Symbol("state")=>"String")
Base.propertynames(::Type{ IoK8sApiEventsV1beta1EventSeries }) = collect(keys(_property_map_IoK8sApiEventsV1beta1EventSeries))
Swagger.property_type(::Type{ IoK8sApiEventsV1beta1EventSeries }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiEventsV1beta1EventSeries[name]))}
Swagger.field_name(::Type{ IoK8sApiEventsV1beta1EventSeries }, property_name::Symbol) =  _property_map_IoK8sApiEventsV1beta1EventSeries[property_name]

function check_required(o::IoK8sApiEventsV1beta1EventSeries)
    (getproperty(o, Symbol("count")) === nothing) && (return false)
    (getproperty(o, Symbol("lastObservedTime")) === nothing) && (return false)
    (getproperty(o, Symbol("state")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiEventsV1beta1EventSeries }, name::Symbol, val)
end


text:
using Test
using Miletus
using Dates: today, Day

using CSV, DataFrames

s = SingleStock()
d1 = today()

bs = CSV.read("bs.csv")

# 1. Test with floats
calls = EuropeanCall.(d1 .+ Day.(bs[:days]), [s], float.(bs[:strikeprice]))
puts = EuropeanPut.(d1 .+ Day.(bs[:days]), [s], float.(bs[:strikeprice]))
m = GeomBMModel.([d1], float.(bs[:startprice]), bs[:interestrate], bs[:carryrate], bs[:sigma])

atol=1e-6; rtol=1e-3

@test all(isapprox.(value.(m,calls), bs[:c], rtol=rtol,atol=atol))
@test all(isapprox.(value.(m,puts), bs[:p],  rtol=rtol,atol=atol))

# @test all(isapprox.(ivol.(m,calls,value.(m,calls)), bs[:sigma]))

crr = CSV.read("crr.csv")

eucalls = EuropeanCall.(d1 .+ Dates.Day.(crr[:days]), [s], float.(crr[:strikeprice]))
euputs  = EuropeanPut.(d1 .+ Dates.Day.(crr[:days]), [s], float.(crr[:strikeprice] ))
amcalls = AmericanCall.(d1 .+ Dates.Day.(crr[:days]), [s], float.(crr[:strikeprice]))
amputs  = AmericanPut.(d1 .+ Dates.Day.(crr[:days]), [s], float.(crr[:strikeprice] ))

m = CRRModel.(d1, d1 .+ Dates.Day.(crr[:days]), crr[:nsteps], float.(crr[:startprice]), crr[:interestrate], crr[:carryrate], crr[:sigma])

@test all(isapprox.(value.(m,eucalls) , crr[:ce], rtol=rtol,atol=atol))
@test all(isapprox.(value.(m,euputs)  , crr[:pe], rtol=rtol,atol=atol))
@test all(isapprox.(value.(m,amcalls) , crr[:ca], rtol=rtol,atol=atol))
@test all(isapprox.(value.(m,amputs)  , crr[:pa], rtol=rtol,atol=atol))

tian = CSV.read("tian.csv")

eucalls = EuropeanCall.(d1 .+ Dates.Day.(tian[:days]), [s], float.(tian[:strikeprice]))
euputs  = EuropeanPut.(d1 .+ Dates.Day.(tian[:days]), [s], float.(tian[:strikeprice] ))
amcalls = AmericanCall.(d1 .+ Dates.Day.(tian[:days]), [s], float.(tian[:strikeprice]))
amputs  = AmericanPut.(d1 .+ Dates.Day.(tian[:days]), [s], float.(tian[:strikeprice] ))

m = TianModel.([d1],d1 .+ Dates.Day.(tian[:days]), tian[:nsteps], float.(tian[:startprice]), tian[:interestrate], tian[:carryrate], tian[:sigma])

@test all(isapprox.(value.(m,eucalls) , tian[:ce], rtol=rtol,atol=atol))
@test all(isapprox.(value.(m,euputs)  , tian[:pe], rtol=rtol,atol=atol))
@test all(isapprox.(value.(m,amcalls) , tian[:ca], rtol=rtol,atol=atol))
@test all(isapprox.(value.(m,amputs)  , tian[:pa], rtol=rtol,atol=atol))

# 2. Test with currencies
calls = EuropeanCall.(d1 .+ Day.(bs[:days]), [s], float.(bs[:strikeprice].*[USD]))
puts = EuropeanPut.(d1 .+ Day.(bs[:days]), [s], float.(bs[:strikeprice].*[USD]))
m = GeomBMModel.([d1], float.(bs[:startprice].*[USD]), bs[:interestrate], bs[:carryrate], bs[:sigma])

atol=1e-6; rtol=1e-3

@test all(isapprox.(value.(m,calls) ./ [USD], bs[:c], rtol=rtol,atol=atol))
@test all(isapprox.(value.(m,puts) ./ [USD], bs[:p],  rtol=rtol,atol=atol))

crr = CSV.read("crr.csv")

eucalls = EuropeanCall.(d1 .+ Dates.Day.(crr[:days]), [s], float.(crr[:strikeprice] .* [USD]))
euputs  = EuropeanPut.(d1 .+ Dates.Day.(crr[:days]), [s], float.(crr[:strikeprice] .* [USD]))
amcalls = AmericanCall.(d1 .+ Dates.Day.(crr[:days]), [s], float.(crr[:strikeprice] .* [USD]))
amputs  = AmericanPut.(d1 .+ Dates.Day.(crr[:days]), [s], float.(crr[:strikeprice] .* [USD]))

m = CRRModel.([d1],d1 .+ Dates.Day.(crr[:days]), crr[:nsteps], float.(crr[:startprice].*[USD]), crr[:interestrate], crr[:carryrate], crr[:sigma])

@test all(isapprox.(value.(m,eucalls) ./ [USD], crr[:ce], rtol=rtol,atol=atol))
@test all(isapprox.(value.(m,euputs)  ./ [USD], crr[:pe], rtol=rtol,atol=atol))
@test all(isapprox.(value.(m,amcalls) ./ [USD], crr[:ca], rtol=rtol,atol=atol))
@test all(isapprox.(value.(m,amputs)  ./ [USD], crr[:pa], rtol=rtol,atol=atol))

tian = CSV.read("tian.csv")

eucalls = EuropeanCall.(d1 .+ Dates.Day.(tian[:days]), [s], float.(tian[:strikeprice] .* [USD]))
euputs  = EuropeanPut.(d1 .+ Dates.Day.(tian[:days]), [s], float.(tian[:strikeprice] .* [USD]))
amcalls = AmericanCall.(d1 .+ Dates.Day.(tian[:days]), [s], float.(tian[:strikeprice] .* [USD]))
amputs  = AmericanPut.(d1 .+ Dates.Day.(tian[:days]), [s], float.(tian[:strikeprice] .* [USD]))

m = TianModel.([d1],d1 .+ Dates.Day.(tian[:days]), tian[:nsteps], float.(tian[:startprice].*[USD]), tian[:interestrate], tian[:carryrate], tian[:sigma])

@test all(isapprox.(value.(m,eucalls) ./ [USD], tian[:ce], rtol=rtol,atol=atol))
@test all(isapprox.(value.(m,euputs)  ./ [USD], tian[:pe], rtol=rtol,atol=atol))
@test all(isapprox.(value.(m,amcalls) ./ [USD], tian[:ca], rtol=rtol,atol=atol))
@test all(isapprox.(value.(m,amputs)  ./ [USD], tian[:pa], rtol=rtol,atol=atol))

jr = CSV.read("jr.csv")

eucalls = EuropeanCall.(d1 .+ Dates.Day.(jr[:days]), [s], float.(jr[:strikeprice] .* [USD]))
euputs  = EuropeanPut.(d1 .+ Dates.Day.(jr[:days]), [s], float.(jr[:strikeprice] .* [USD]))
amcalls = AmericanCall.(d1 .+ Dates.Day.(jr[:days]), [s], float.(jr[:strikeprice] .* [USD]))
amputs  = AmericanPut.(d1 .+ Dates.Day.(jr[:days]), [s], float.(jr[:strikeprice] .* [USD]))

m = JRModel.([d1],d1 .+ Dates.Day.(jr[:days]), jr[:nsteps], float.(jr[:startprice].*[USD]), jr[:interestrate], jr[:carryrate], jr[:sigma])

@test all(isapprox.(value.(m,eucalls) ./ [USD], jr[:ce], rtol=rtol,atol=atol))
@test all(isapprox.(value.(m,euputs)  ./ [USD], jr[:pe], rtol=rtol,atol=atol))
@test all(isapprox.(value.(m,amcalls) ./ [USD], jr[:ca], rtol=rtol,atol=atol))
@test all(isapprox.(value.(m,amputs)  ./ [USD], jr[:pa], rtol=rtol,atol=atol))

# Test Black-76 model
d2 = d1 + Day(120)
cm = CoreForwardModel(d1, d2, 110.00USD, 0.02)
m = GeomBMModel(cm, 0.3)
o = EuropeanCall( d2, s, 100.00USD)
@test isapprox(value(m,o), 13.1833USD, rtol=1e-5)


text:
# setup test file
dir = mkpath(abspath(@__FILE__, "..", "data", "ball"))
cd(dir)
mkpath("sub")
mkpath("empty")
for (fn, size) in (("a", 10), ("b", 100), ("c", 0), ("d", 3)), sub in ("", "sub",)
    open(joinpath(sub, fn * ".dat"), "w") do io
        write(io, fn^size)
    end
end
cd(abspath(dir, ".."))
tarfile = "test.tar"
run(`sh -c "tar -cf $tarfile ball"`)
run(`sh -c "gzip <$tarfile >$tarfile.gz"`)

function opener(file::AbstractString, f=nothing; close_stream::Bool=false)
    source = open(file)
    if endswith(file, ".tar.gz") || endswith(file, ".tgz")
        source = GzipDecompressorStream(source)
    end
    TarIterator(source, f, close_stream=close_stream)
end

@testset "reading all elements of $tarfile" for tarfile = (tarfile, tarfile*".gz")
    ti = opener(tarfile)
    @test ti !== nothing
    s = iterate(ti)
    @test s isa Tuple
    (h, io), st = s
    @test h isa Tar.Header
    @test io isa BoundedInputStream
    @test h.type == :directory
    @test h.path == "ball/"
    @test close(ti) === nothing
    ti = opener(tarfile)
    @test open(ti) isa BoundedInputStream

    seekstart(ti)
    res1 = Any[]
    for (h, io) in ti
        push!(res1, h)
        push!(res1, read(io, String))
    end
    @test length(res1) == 22

    seekstart(ti)
    res2 = Any[]
    open(ti) do h, io
        push!(res2, h)
        push!(res2, read(io, String))
    end
    @test length(res2) == 22
    @test res1 == res2
end

@testset "iterations predicate $p" for (p,m) in [(nothing, 11),
                                                 ("ball/b.dat", 1),
                                                 (r"^ball/(|.*/)a\.dat", 2),
                                                 (:file, 8),
                                                 (h->h.type == :file, 8),
                                                 ((:file, h->h.size>10), 2),
                                                 ([:directory, h->h.size>10], 5)]
    n = 0
    open(TarIterator(tarfile, p)) do h, io
        n += 1
    end
    @test n == m
end


text:
module ThirdWay

import Base: +, isempty, isless, length, merge!, push!, show
import Base: Callable, call

include("activation_point.jl")
include("event.jl")
include("schedule.jl")
include("schedule_elaborations.jl")
include("util.jl")

end


text:
module BinaryDecisionDiagrams

import Base: print, println, show, string, haskey, in, ==, !=, |, &, ~

export Ordering, ListOrdering, BDD, OBDD, restrict, descendents, ancestors,
       applyoperator, changeordering, set_dynamic_ordering,
       set_static_ordering, SAT_assignments

include("Ordering.jl")
include("BDD.jl")
include("OBDD.jl")
include("parsing.jl")

end # module


text:
"""
Unit tests for ControlUnit type.

------------------------------------------------------------------------------
COPYRIGHT/LICENSE. This file is part of the XYZ package. It is subject to
the license terms in the LICENSE file found in the top-level directory of
this distribution. No part of the XYZ package, including this file, may be
copied, modified, propagated, or distributed except according to the terms
contained in the LICENSE file.
------------------------------------------------------------------------------
"""
# --- Imports

# Standard library
using Test

# External packages
using ZMQ

# Boid.jl package
using Boid

# --- Test fixtures

# Concrete subtype of AbstractControlLogicCore
include("fixtures/TestControlLogicCore.jl")

# Constants
control_url = construct_ipc_url(".", "control-unit-test.zmq")
control_ipc_path = control_url[7:end]

# --- Set up/tear down methods

function _tearDown()
    rm(control_ipc_path, force=true)
end

# --- Tests

@testset "ControlUnit: constructor" begin
    # --- Preparations

    logic_core = TestControlLogicCore()

    # --- Tests

    # copy_logic_core == true
    control_unit = ControlUnit(logic_core, control_url)

    @test control_unit.logic_core isa TestControlLogicCore
    @test control_unit.logic_core !== logic_core
    @test control_unit.url == control_url
    @test control_unit.socket isa Socket

    # copy_logic_core == false
    control_unit = ControlUnit(logic_core, control_url, copy_logic_core=false)

    @test control_unit.logic_core isa TestControlLogicCore
    @test control_unit.logic_core === logic_core
    @test control_unit.url == control_url
    @test control_unit.socket isa Socket

    # --- Clean up

    _tearDown()
end

@testset "ControlUnit: run!()" begin
    # --- Preparations

    logic_core = TestControlLogicCore()
    control_unit = ControlUnit(logic_core, control_url)

    # --- Tests

    # Start control signal processing
    run!(control_unit)

    # Send Boid.START signal
    socket = Socket(REQ)
    connect(socket, control_url)
    send(socket, Boid.START)
    response = recv(socket, ControlUnitResponse)
    @test is_running(control_unit)
    @test !is_terminated(control_unit)
    @test control_unit.logic_core.count == 0
    @test response == Boid.SUCCESS

    # Send INCREMENT signal
    socket = Socket(REQ)
    connect(socket, control_url)
    send(socket, INCREMENT)
    response = recv(socket, ControlUnitResponse)
    @test is_running(control_unit)
    @test !is_terminated(control_unit)
    @test control_unit.logic_core.count == 1
    @test response == Boid.SUCCESS

    # Send unknown signal
    socket = Socket(REQ)
    connect(socket, control_url)
    send(socket, 100)
    response = recv(socket, ControlUnitResponse)
    @test response == Boid.UNKNOWN_SIGNAL

    # Send Boid.STOP signal
    socket = Socket(REQ)
    connect(socket, control_url)
    send(socket, Boid.STOP)
    response = recv(socket, ControlUnitResponse)
    @test !is_running(control_unit)
    @test !is_terminated(control_unit)
    @test control_unit.logic_core.count == 1
    @test response == Boid.SUCCESS

    # Send Boid.TERMINATED signal
    socket = Socket(REQ)
    connect(socket, control_url)
    send(socket, Boid.TERMINATE)
    response = recv(socket, ControlUnitResponse)
    @test !is_running(control_unit)
    @test is_terminated(control_unit)
    @test response == Boid.SUCCESS

    # --- Clean up

    _tearDown()
end

@testset "ControlUnit: process_control_signal!()" begin
    @test_skip 1
end

@testset "ControlUnit: is_running()" begin
    @test_skip 1
end

@testset "ControlUnit: is_terminated()" begin
    @test_skip 1
end

@testset "ControlUnit: wait_for_input_ready()" begin
    @test_skip 1
end


text:
export CuPtr

struct CuPtr{T}
    ptr::Ptr{T}
end

Base.cconvert(::Type{Ptr{T}}, x::CuPtr{T}) where T = x.ptr

#=
mutable struct CuPtr{T}
    ptr::Ptr{T}
    size::Int
    dev::Int
end

CuPtr(ptr::Ptr{T}, size::Int) where T = CuPtr(ptr, size, getdevice())
CuPtr{T}() where T = CuPtr(Ptr{T}(C_NULL), 0)

Base.cconvert(::Type{Ptr{T}}, x::CuPtr) where T = Ptr{T}(x.ptr)
Base.pointer(x::CuPtr{T}, index::Int=1) where T = x.ptr + sizeof(T)*(index-1)
=#


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

const line_width = 60

const alu = string(
   "GGCCGGGCGCGGTGGCTCACGCCTGTAATCCCAGCACTTTGG",
   "GAGGCCGAGGCGGGCGGATCACCTGAGGTCAGGAGTTCGAGA",
   "CCAGCCTGGCCAACATGGTGAAACCCCGTCTCTACTAAAAAT",
   "ACAAAAATTAGCCGGGCGTGGTGGCGCGCGCCTGTAATCCCA",
   "GCTACTCGGGAGGCTGAGGCAGGAGAATCGCTTGAACCCGGG",
   "AGGCGGAGGTTGCAGTGAGCCGAGATCGCGCCACTGCACTCC",
   "AGCCTGGGCGACAGAGCGAGACTCCGTCTCAAAAA")

const iub1 = b"acgtBDHKMNRSVWY"
const iub2 = [0.27, 0.12, 0.12, 0.27, 0.02,0.02, 0.02, 0.02, 0.02, 0.02,0.02, 0.02, 0.02, 0.02, 0.02]

const homosapiens1 = b"acgt"
const homosapiens2 = [0.3029549426680, 0.1979883004921,0.1975473066391, 0.3015094502008]

const IM  = 139968.0
const IA  =   3877.0
const IC  =  29573.0

function gen_random()
    global rng_state = ((rng_state::Float64 * IA + IC) % IM) / IM
end
function repeat_fasta(src, n)
    k = length(src)
    return string(src, src, src[1:n % k])
end
function choose_char(cs)
    k = length(cs)
    r = gen_random()
    r < cs[1] && return 1
    a = 1
    b = k
    while b > a + 1
        c = fld(a + b, 2)
        if r < cs[c]; b = c; else a = c; end
    end
    b
end
function random_fasta(symb, pr, n)
    cs = cumsum(pr)
    line = Array{UInt8}(line_width)
    k = n
    while k > 0
        m = min(k, line_width)
        for i = 1:m
            line[i] = symb[choose_char(cs)]
        end
        k -= line_width
    end
end

rng_state = 42.0
function fasta(n=25000000)
    repeat_fasta(alu, 2n)
    random_fasta(iub1, iub2, 3n)
    random_fasta(homosapiens1, homosapiens2, 5n)
end


text:
using ..Continuations: as,
    ContinuationSweep, ContinuationSolution,
    ContinuationCache, ContinuationSolver

abstract type AbstractSpecialPoint{tkind <: TimeKind, ptType} end
TimeKind(::Type{<: AbstractSpecialPoint{tkind}}) where tkind = tkind()

struct SpecialPoint{tkind, ptType, uType, JType,
                    } <: AbstractSpecialPoint{tkind, ptType}
    timekind::tkind
    point_type::ptType
    point_index::Int
    u::uType
    J::JType
    sweep::WeakRef
end
# TODO: put BifurcationProblem in SpecialPoint/Interval struct so that
# init(point) can initiate codim+1 bifurcation solver.

struct SpecialPointInterval{tkind, ptType, uType, JType,
                            } <: AbstractSpecialPoint{tkind, ptType}
    timekind::tkind
    point_type::ptType
    point_index::Int
    u0::uType
    u1::uType
    # J0::JType
    J1::JType
    sweep::WeakRef
end

struct BifurcationSweep{
        tkind <: TimeKind,
        ptType,
        S <: ContinuationSweep,
        JType <: AbstractArray,
        eType <: AbstractArray,
        siType <: SpecialPointInterval{tkind, ptType},
        spType <: SpecialPoint{tkind, ptType},
        }
    timekind::tkind
    super::S
    jacobians::Vector{JType}
    eigvals::Vector{eType}
    special_intervals::Vector{siType}
    special_points::Vector{spType}
end
# See: [[../continuations/solution.jl::ContinuationSweep]]

BifurcationSweep{tkind, ptType, S, JType, eType, siType, spType}(super::S) where{
    tkind <: TimeKind,
    ptType,
    S <: ContinuationSweep,
    JType <: AbstractArray,
    eType <: AbstractArray,
    siType <: SpecialPointInterval{tkind, ptType},
    spType <: SpecialPoint{tkind, ptType},
} = BifurcationSweep{tkind, ptType, S, JType, eType, siType, spType}(
    tkind(),
    super,
    JType[],
    eType[],
    siType[],
    spType[],
)

TimeKind(::Type{<: BifurcationSweep{tkind}}) where tkind = tkind()

Base.length(sweep::BifurcationSweep) = length(as(sweep, ContinuationSweep))

PointTypeType(::Type{<: BifurcationSweep{_tk, ptType}}) where {_tk, ptType} =
    ptType

point_type_type(::T) where {T <: BifurcationSweep} = PointTypeType(T)

contkind(sweep::BifurcationSweep) = contkind(problem_of(sweep)) # TODO: don't
# this is terrible!

function check_sweep_length(sweep)
    @assert length(sweep) == length(sweep.jacobians) == length(sweep.eigvals)
end

function sweeptype(prob::BifurcationProblem,
                   solver::ContinuationSolver,
                   JType::Type = typeof(solver.cache.J),
                   eType::Type = typeof(eigvals_prototype(prob, solver.cache)),
                   )
    tkind = typeof(timekind(solver.cache))
    S = eltype(solver.sol.sweeps)
    uType = eltype(S)
    ptType = point_type_type(prob)
    siType = SpecialPointInterval{tkind, ptType, uType, JType}
    spType = SpecialPoint{tkind, ptType, uType, JType}
    return BifurcationSweep{tkind, ptType, S, JType, eType, siType, spType}
end

function BifurcationSweep(super::ContinuationSweep, solver::ContinuationSolver)
    return sweeptype(solver.prob, solver)(super)
end

function push_special_point!(sweep::BifurcationSweep, point::SpecialPoint)
    @assert point.point_index < length(sweep)

    push!(sweep.special_points, point)

    # TODO: Add two new points rather than replacing the old one.
    # At the moment, just "move" the recorded points at the resolved
    # `point` so that there would be no discrepancy between the
    # plotted special points and stability of the lines.
    cont_sweep = as(sweep, ContinuationSweep)
    cont_sweep.u[point.point_index] = point.u
    cont_sweep.u[point.point_index + 1] = point.u

    return sweep
end

struct BifurcationSolution{S <: ContinuationSolution,
                           W <: BifurcationSweep}
    super::S
    sweeps::Vector{W}
end

BifurcationSolution(super::ContinuationSolution,
                    SweesiType::Type{<: BifurcationSweep}) =
    BifurcationSolution(super, SweesiType[])

TimeKind(::Type{<: BifurcationSolution{_, W}}) where {_, W} = TimeKind(W)


text:
using MultiResolutionFilters
using Test
using DataFrames, CSV
using Distributions
using StaticArrays

@testset "MultiResolutionFilters.jl" begin
    data = CSV.read("test_data.csv", DataFrame)
    locations = [SVector(row...) for row in eachrow(data[:, [:x, :y]])]
    observations = data.z
    nparticles = 1000
    state_particles = randn(nparticles)

    loglik(state, obs) = loglikelihood(Normal(state, 0.1), obs)
    predict_downscale(state, child_area) = state + 0.5randn()
    observe(state) = state + randn()
    needsrefinement(cell) = area(cell.boundary) > 250#dA[6]
    jitter(state) = state

    refinery = ParticleRefinery(predict_downscale, observe, loglik, needsrefinement, jitter)
    mrf = MultiResolutionPF(locations, observations, refinery, state_particles)
    adaptive_filter!(mrf)

    @test 1 == 1
end


text:
using Documenter, nGraph

makedocs(
    modules = [nGraph],
    format = :html,
    checkdocs = :exports,
    html_prettyurls = get(ENV, "CI", nothing) == "true",
    sitename = "nGraph.jl",
    pages = Any[
        "index.md",
        "ngraph" => [
            "ngraph/types.md",
        ],
    ]
)

deploydocs(
    repo = "github.com/hildebrandmw/nGraph.jl.git",
)


text:
using Documenter
using ReadCoverageDistributions

makedocs(
    sitename = "ReadCoverageDistributions",
    format = Documenter.HTML(prettyurls = true),
    modules = [ReadCoverageDistributions]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "github.com:Dinesh-Adhithya-H/ReadCoverageDistributions.jl.git"
)


text:
using ModernGL
using Images
using Serialization

function main()
  # write ascii to file
  fp = open("ascii.txt", "w")
  char_per_line = 32
  text = ""
  for j in 0:255
    c = Char(j)
    text = string(text, iscntrl(c) ? '?' : c)
    if rem(j + 1, char_per_line) == 0 && j != 255
      text = string(text, '\n')
    end
  end
  write(fp, text)
  close(fp)

  # convert ascii to a bitmap with ImageMagick
  run(`convert -background transparent -font Courier -pointsize 64 label:@ascii.txt font.bmp`)

  # write the bmp to row major binary data for OpenGL
  img = load("font.bmp")
  (height, width) = size(img)

  pixels = fill(GLfloat(0.0), (4 * width, height))
  for i in 1:height
    for j in 1:width
      r = GLfloat(img[i, j].r)
      g = GLfloat(img[i, j].g)
      b = GLfloat(img[i, j].b)
      a = GLfloat(img[i, j].alpha)
      pixels[4 * (j - 1) + 1, i] = r
      pixels[4 * (j - 1) + 2, i] = g
      pixels[4 * (j - 1) + 3, i] = b
      pixels[4 * (j - 1) + 4, i] = a
    end
  end

  fp = open("font.bin", "w")
  serialize(fp, pixels)
  close(fp)
  return
end

main()


text:
#export @mkdeepconvert3

isstringint(ex::Expr) = ex.head == :macrocall &&
    (ex.args[1] == Symbol("@int128_str") || ex.args[1] == Symbol("@bigint_str"))

macro mkdeepconvert(ff, ccfunc)
    f = esc(ff)
    cfunc = esc(ccfunc)
    quote
        function ($f)(ex::Expr)
            if ex.args[1] == ://
                evex = eval(ex)
                return convert(promote_type($cfunc, typeof(evex)), evex)
            else
              eval(Expr(ex.head, map(
              (x) ->
                begin
                    tx = typeof(x)
                    if tx <: Real
                        return ($cfunc)(x)
                    elseif  tx == Expr
                        return ($f)(x)
                    else
                        return x
                    end
                end,
               ex.args)...))
            end
        end
        ($f)(x::AbstractString) = ($f)(Meta.parse(x))
        ($f)(x) = ($cfunc)(x)
    end
end

macro mkdeepconvert1(ff, ccfunc)
    f = esc(ff)
    cfunc = esc(ccfunc)
    quote
        function $(f)(ex::Expr)
             eval(Expr(ex.head, map(
              (x) ->
             begin
               tx = typeof(x)
               if tx <: Real
                   return ($cfunc)(x)
               elseif  tx == Expr
                       return ($f)(x)
               else
                   return x
               end
             end,
             ex.args)...))
        end
        ($f)(x::AbstractString) = ($f)(Meta.parse(x))
        ($f)(x) = ($cfunc)(x)
    end
end

macro mkdeepconvert2(ff, ccfunc, targtype)
    f = esc(ff)
    cfunc = esc(ccfunc)
    quote
        function $(f)(ex::Expr)
             Expr(ex.head, map(
              (x) ->
             begin
               tx = typeof(x)
               if tx <: $targtype
                   return ($cfunc)(x)
               elseif  tx == Expr
                       return ($f)(x)
               else
                   return x
               end
             end,
             ex.args)...)
        end
        ($f)(x::AbstractString) = ($f)(Meta.parse(x))
        ($f)(x) = ($cfunc)(x)
    end
end

# ff is the name of the function to be generated.
# ccfunc is the existing function that does the conversion.
macro mkdeepconvert3(ff, ccfunc, targtype)
   f = esc(ff); cfunc = esc(ccfunc)
   quote
     function $(f)(ex::Expr)
       isstringint(ex) && return Expr(:call, $cfunc, ex)
         Expr(ex.head, map((x) -> (typeof(x) <: $targtype ?
           ($cfunc)(x) : typeof(x) == Expr ? ($f)(x) : x), ex.args)...)
     end
     ($f)(x::AbstractString) = ($f)(Meta.parse(x))
     ($f)(x) = ($cfunc)(x)
   end
end


text:
if Base.HOME_PROJECT[] !== nothing
    Base.HOME_PROJECT[] = abspath(Base.HOME_PROJECT[])
end

push!(LOAD_PATH, joinpath("..", "src"))

using Documenter, RunLengthArrays

makedocs(;
    modules = [RunLengthArrays],
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    pages = [
        "Home" => "index.md",
    ],
    repo = "https://github.com/ziotom78/RunLengthArrays.jl/blob/{commit}{path}#L{line}",
    sitename = "RunLengthArrays.jl",
    authors = "Maurizio Tomasi",)

deploydocs(;
    repo = "github.com/ziotom78/RunLengthArrays.jl",)


text:
import Blowfish

using Test

mutable struct CryptTest
	key::Array{UInt8, 1}
	in::Array{UInt8, 1}
	out::Array{UInt8, 1}
end

const encryptTests = (
	CryptTest(
		UInt8[0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00],
		UInt8[0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00],
		UInt8[0x4E, 0xF9, 0x97, 0x45, 0x61, 0x98, 0xDD, 0x78]),
	CryptTest(
		UInt8[0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF],
		UInt8[0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF],
		UInt8[0x51, 0x86, 0x6F, 0xD5, 0xB8, 0x5E, 0xCB, 0x8A]),
	CryptTest(
		UInt8[0x30, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00],
		UInt8[0x10, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x01],
		UInt8[0x7D, 0x85, 0x6F, 0x9A, 0x61, 0x30, 0x63, 0xF2]),
	CryptTest(
		UInt8[0x11, 0x11, 0x11, 0x11, 0x11, 0x11, 0x11, 0x11],
		UInt8[0x11, 0x11, 0x11, 0x11, 0x11, 0x11, 0x11, 0x11],
		UInt8[0x24, 0x66, 0xDD, 0x87, 0x8B, 0x96, 0x3C, 0x9D]),
	CryptTest(
		UInt8[0x01, 0x23, 0x45, 0x67, 0x89, 0xAB, 0xCD, 0xEF],
		UInt8[0x11, 0x11, 0x11, 0x11, 0x11, 0x11, 0x11, 0x11],
		UInt8[0x61, 0xF9, 0xC3, 0x80, 0x22, 0x81, 0xB0, 0x96]),
	CryptTest(
		UInt8[0x11, 0x11, 0x11, 0x11, 0x11, 0x11, 0x11, 0x11],
		UInt8[0x01, 0x23, 0x45, 0x67, 0x89, 0xAB, 0xCD, 0xEF],
		UInt8[0x7D, 0x0C, 0xC6, 0x30, 0xAF, 0xDA, 0x1E, 0xC7]),
	CryptTest(
		UInt8[0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00],
		UInt8[0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00],
		UInt8[0x4E, 0xF9, 0x97, 0x45, 0x61, 0x98, 0xDD, 0x78]),
	CryptTest(
		UInt8[0xFE, 0xDC, 0xBA, 0x98, 0x76, 0x54, 0x32, 0x10],
		UInt8[0x01, 0x23, 0x45, 0x67, 0x89, 0xAB, 0xCD, 0xEF],
		UInt8[0x0A, 0xCE, 0xAB, 0x0F, 0xC6, 0xA0, 0xA2, 0x8D]),
	CryptTest(
		UInt8[0x7C, 0xA1, 0x10, 0x45, 0x4A, 0x1A, 0x6E, 0x57],
		UInt8[0x01, 0xA1, 0xD6, 0xD0, 0x39, 0x77, 0x67, 0x42],
		UInt8[0x59, 0xC6, 0x82, 0x45, 0xEB, 0x05, 0x28, 0x2B]),
	CryptTest(
		UInt8[0x01, 0x31, 0xD9, 0x61, 0x9D, 0xC1, 0x37, 0x6E],
		UInt8[0x5C, 0xD5, 0x4C, 0xA8, 0x3D, 0xEF, 0x57, 0xDA],
		UInt8[0xB1, 0xB8, 0xCC, 0x0B, 0x25, 0x0F, 0x09, 0xA0]),
	CryptTest(
		UInt8[0x07, 0xA1, 0x13, 0x3E, 0x4A, 0x0B, 0x26, 0x86],
		UInt8[0x02, 0x48, 0xD4, 0x38, 0x06, 0xF6, 0x71, 0x72],
		UInt8[0x17, 0x30, 0xE5, 0x77, 0x8B, 0xEA, 0x1D, 0xA4]),
	CryptTest(
		UInt8[0x38, 0x49, 0x67, 0x4C, 0x26, 0x02, 0x31, 0x9E],
		UInt8[0x51, 0x45, 0x4B, 0x58, 0x2D, 0xDF, 0x44, 0x0A],
		UInt8[0xA2, 0x5E, 0x78, 0x56, 0xCF, 0x26, 0x51, 0xEB]),
	CryptTest(
		UInt8[0x04, 0xB9, 0x15, 0xBA, 0x43, 0xFE, 0xB5, 0xB6],
		UInt8[0x42, 0xFD, 0x44, 0x30, 0x59, 0x57, 0x7F, 0xA2],
		UInt8[0x35, 0x38, 0x82, 0xB1, 0x09, 0xCE, 0x8F, 0x1A]),
	CryptTest(
		UInt8[0x01, 0x13, 0xB9, 0x70, 0xFD, 0x34, 0xF2, 0xCE],
		UInt8[0x05, 0x9B, 0x5E, 0x08, 0x51, 0xCF, 0x14, 0x3A],
		UInt8[0x48, 0xF4, 0xD0, 0x88, 0x4C, 0x37, 0x99, 0x18]),
	CryptTest(
		UInt8[0x01, 0x70, 0xF1, 0x75, 0x46, 0x8F, 0xB5, 0xE6],
		UInt8[0x07, 0x56, 0xD8, 0xE0, 0x77, 0x47, 0x61, 0xD2],
		UInt8[0x43, 0x21, 0x93, 0xB7, 0x89, 0x51, 0xFC, 0x98]),
	CryptTest(
		UInt8[0x43, 0x29, 0x7F, 0xAD, 0x38, 0xE3, 0x73, 0xFE],
		UInt8[0x76, 0x25, 0x14, 0xB8, 0x29, 0xBF, 0x48, 0x6A],
		UInt8[0x13, 0xF0, 0x41, 0x54, 0xD6, 0x9D, 0x1A, 0xE5]),
	CryptTest(
		UInt8[0x07, 0xA7, 0x13, 0x70, 0x45, 0xDA, 0x2A, 0x16],
		UInt8[0x3B, 0xDD, 0x11, 0x90, 0x49, 0x37, 0x28, 0x02],
		UInt8[0x2E, 0xED, 0xDA, 0x93, 0xFF, 0xD3, 0x9C, 0x79]),
	CryptTest(
		UInt8[0x04, 0x68, 0x91, 0x04, 0xC2, 0xFD, 0x3B, 0x2F],
		UInt8[0x26, 0x95, 0x5F, 0x68, 0x35, 0xAF, 0x60, 0x9A],
		UInt8[0xD8, 0x87, 0xE0, 0x39, 0x3C, 0x2D, 0xA6, 0xE3]),
	CryptTest(
		UInt8[0x37, 0xD0, 0x6B, 0xB5, 0x16, 0xCB, 0x75, 0x46],
		UInt8[0x16, 0x4D, 0x5E, 0x40, 0x4F, 0x27, 0x52, 0x32],
		UInt8[0x5F, 0x99, 0xD0, 0x4F, 0x5B, 0x16, 0x39, 0x69]),
	CryptTest(
		UInt8[0x1F, 0x08, 0x26, 0x0D, 0x1A, 0xC2, 0x46, 0x5E],
		UInt8[0x6B, 0x05, 0x6E, 0x18, 0x75, 0x9F, 0x5C, 0xCA],
		UInt8[0x4A, 0x05, 0x7A, 0x3B, 0x24, 0xD3, 0x97, 0x7B]),
	CryptTest(
		UInt8[0x58, 0x40, 0x23, 0x64, 0x1A, 0xBA, 0x61, 0x76],
		UInt8[0x00, 0x4B, 0xD6, 0xEF, 0x09, 0x17, 0x60, 0x62],
		UInt8[0x45, 0x20, 0x31, 0xC1, 0xE4, 0xFA, 0xDA, 0x8E]),
	CryptTest(
		UInt8[0x02, 0x58, 0x16, 0x16, 0x46, 0x29, 0xB0, 0x07],
		UInt8[0x48, 0x0D, 0x39, 0x00, 0x6E, 0xE7, 0x62, 0xF2],
		UInt8[0x75, 0x55, 0xAE, 0x39, 0xF5, 0x9B, 0x87, 0xBD]),
	CryptTest(
		UInt8[0x49, 0x79, 0x3E, 0xBC, 0x79, 0xB3, 0x25, 0x8F],
		UInt8[0x43, 0x75, 0x40, 0xC8, 0x69, 0x8F, 0x3C, 0xFA],
		UInt8[0x53, 0xC5, 0x5F, 0x9C, 0xB4, 0x9F, 0xC0, 0x19]),
	CryptTest(
		UInt8[0x4F, 0xB0, 0x5E, 0x15, 0x15, 0xAB, 0x73, 0xA7],
		UInt8[0x07, 0x2D, 0x43, 0xA0, 0x77, 0x07, 0x52, 0x92],
		UInt8[0x7A, 0x8E, 0x7B, 0xFA, 0x93, 0x7E, 0x89, 0xA3]),
	CryptTest(
		UInt8[0x49, 0xE9, 0x5D, 0x6D, 0x4C, 0xA2, 0x29, 0xBF],
		UInt8[0x02, 0xFE, 0x55, 0x77, 0x81, 0x17, 0xF1, 0x2A],
		UInt8[0xCF, 0x9C, 0x5D, 0x7A, 0x49, 0x86, 0xAD, 0xB5]),
	CryptTest(
		UInt8[0x01, 0x83, 0x10, 0xDC, 0x40, 0x9B, 0x26, 0xD6],
		UInt8[0x1D, 0x9D, 0x5C, 0x50, 0x18, 0xF7, 0x28, 0xC2],
		UInt8[0xD1, 0xAB, 0xB2, 0x90, 0x65, 0x8B, 0xC7, 0x78]),
	CryptTest(
		UInt8[0x1C, 0x58, 0x7F, 0x1C, 0x13, 0x92, 0x4F, 0xEF],
		UInt8[0x30, 0x55, 0x32, 0x28, 0x6D, 0x6F, 0x29, 0x5A],
		UInt8[0x55, 0xCB, 0x37, 0x74, 0xD1, 0x3E, 0xF2, 0x01]),
	CryptTest(
		UInt8[0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01],
		UInt8[0x01, 0x23, 0x45, 0x67, 0x89, 0xAB, 0xCD, 0xEF],
		UInt8[0xFA, 0x34, 0xEC, 0x48, 0x47, 0xB2, 0x68, 0xB2]),
	CryptTest(
		UInt8[0x1F, 0x1F, 0x1F, 0x1F, 0x0E, 0x0E, 0x0E, 0x0E],
		UInt8[0x01, 0x23, 0x45, 0x67, 0x89, 0xAB, 0xCD, 0xEF],
		UInt8[0xA7, 0x90, 0x79, 0x51, 0x08, 0xEA, 0x3C, 0xAE]),
	CryptTest(
		UInt8[0xE0, 0xFE, 0xE0, 0xFE, 0xF1, 0xFE, 0xF1, 0xFE],
		UInt8[0x01, 0x23, 0x45, 0x67, 0x89, 0xAB, 0xCD, 0xEF],
		UInt8[0xC3, 0x9E, 0x07, 0x2D, 0x9F, 0xAC, 0x63, 0x1D]),
	CryptTest(
		UInt8[0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00],
		UInt8[0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF],
		UInt8[0x01, 0x49, 0x33, 0xE0, 0xCD, 0xAF, 0xF6, 0xE4]),
	CryptTest(
		UInt8[0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF],
		UInt8[0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00],
		UInt8[0xF2, 0x1E, 0x9A, 0x77, 0xB7, 0x1C, 0x49, 0xBC]),
	CryptTest(
		UInt8[0x01, 0x23, 0x45, 0x67, 0x89, 0xAB, 0xCD, 0xEF],
		UInt8[0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00],
		UInt8[0x24, 0x59, 0x46, 0x88, 0x57, 0x54, 0x36, 0x9A]),
	CryptTest(
		UInt8[0xFE, 0xDC, 0xBA, 0x98, 0x76, 0x54, 0x32, 0x10],
		UInt8[0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF],
		UInt8[0x6B, 0x5C, 0x5A, 0x9C, 0x5D, 0x9E, 0x0A, 0x5A]
	),
)::NTuple{34, CryptTest}

@testset "test cipher encrypt" begin
	for cryptTest in encryptTests
			c = Blowfish.NewCipher(cryptTest.key)
			ct = zeros(UInt8, length(cryptTest.out))
			Blowfish.Encrypt!(c, ct, cryptTest.in)
			for i in 1:length(ct)
				@test ct[i] == cryptTest.out[i]
			end
	end
end

@testset "test cipher decrypt" begin
	for cryptTest in encryptTests
			c = Blowfish.NewCipher(cryptTest.key)
			pt = zeros(UInt8, length(cryptTest.in))
			Blowfish.Decrypt!(c, pt, cryptTest.out)
			for i in 1:length(pt)
				@test pt[i] == cryptTest.in[i]
			end
	end
end

@testset "test salted key length" begin
	try
		Blowfish.NewSaltedCipher(UInt8[], UInt8['a'])
	catch ex
		@test isa(ex, Blowfish.KeySizeError)
	end
	key = Array{UInt8, 1}("012345678901234567890123456789012345678901234567890123456")
	Blowfish.NewSaltedCipher(key, UInt8['a'])
end

const saltedVectors = (
	UInt8[0x0c, 0x82, 0x3b, 0x7b, 0x8d, 0x01, 0x4b, 0x7e],
	UInt8[0xd1, 0xe1, 0x93, 0xf0, 0x70, 0xa6, 0xdb, 0x12],
	UInt8[0xfc, 0x5e, 0xba, 0xde, 0xcb, 0xf8, 0x59, 0xad],
	UInt8[0x8a, 0x0c, 0x76, 0xe7, 0xdd, 0x2c, 0xd3, 0xa8],
	UInt8[0x2c, 0xcb, 0x7b, 0xee, 0xac, 0x7b, 0x7f, 0xf8],
	UInt8[0xbb, 0xf6, 0x30, 0x6f, 0xe1, 0x5d, 0x62, 0xbf],
	UInt8[0x97, 0x1e, 0xc1, 0x3d, 0x3d, 0xe0, 0x11, 0xe9],
	UInt8[0x06, 0xd7, 0x4d, 0xb1, 0x80, 0xa3, 0xb1, 0x38],
	UInt8[0x67, 0xa1, 0xa9, 0x75, 0x0e, 0x5b, 0xc6, 0xb4],
	UInt8[0x51, 0x0f, 0x33, 0x0e, 0x4f, 0x67, 0xd2, 0x0c],
	UInt8[0xf1, 0x73, 0x7e, 0xd8, 0x44, 0xea, 0xdb, 0xe5],
	UInt8[0x14, 0x0e, 0x16, 0xce, 0x7f, 0x4a, 0x9c, 0x7b],
	UInt8[0x4b, 0xfe, 0x43, 0xfd, 0xbf, 0x36, 0x04, 0x47],
	UInt8[0xb1, 0xeb, 0x3e, 0x15, 0x36, 0xa7, 0xbb, 0xe2],
	UInt8[0x6d, 0x0b, 0x41, 0xdd, 0x00, 0x98, 0x0b, 0x19],
	UInt8[0xd3, 0xce, 0x45, 0xce, 0x1d, 0x56, 0xb7, 0xfc],
	UInt8[0xd9, 0xf0, 0xfd, 0xda, 0xc0, 0x23, 0xb7, 0x93],
	UInt8[0x4c, 0x6f, 0xa1, 0xe4, 0x0c, 0xa8, 0xca, 0x57],
	UInt8[0xe6, 0x2f, 0x28, 0xa7, 0x0c, 0x94, 0x0d, 0x08],
	UInt8[0x8f, 0xe3, 0xf0, 0xb6, 0x29, 0xe3, 0x44, 0x03],
	UInt8[0xff, 0x98, 0xdd, 0x04, 0x45, 0xb4, 0x6d, 0x1f],
	UInt8[0x9e, 0x45, 0x4d, 0x18, 0x40, 0x53, 0xdb, 0xef],
	UInt8[0xb7, 0x3b, 0xef, 0x29, 0xbe, 0xa8, 0x13, 0x71],
	UInt8[0x02, 0x54, 0x55, 0x41, 0x8e, 0x04, 0xfc, 0xad],
	UInt8[0x6a, 0x0a, 0xee, 0x7c, 0x10, 0xd9, 0x19, 0xfe],
	UInt8[0x0a, 0x22, 0xd9, 0x41, 0xcc, 0x23, 0x87, 0x13],
	UInt8[0x6e, 0xff, 0x1f, 0xff, 0x36, 0x17, 0x9c, 0xbe],
	UInt8[0x79, 0xad, 0xb7, 0x40, 0xf4, 0x9f, 0x51, 0xa6],
	UInt8[0x97, 0x81, 0x99, 0xa4, 0xde, 0x9e, 0x9f, 0xb6],
	UInt8[0x12, 0x19, 0x7a, 0x28, 0xd0, 0xdc, 0xcc, 0x92],
	UInt8[0x81, 0xda, 0x60, 0x1e, 0x0e, 0xdd, 0x65, 0x56],
	UInt8[0x7d, 0x76, 0x20, 0xb2, 0x73, 0xc9, 0x9e, 0xee],
)::NTuple{32, Array{UInt8, 1}}

@testset "test salted cipher" begin
	key, salt = zeros(UInt8, 32), zeros(UInt8, 32)
	for i in 1:32
		key[i] = UInt8(i-1 & 0xff)
		salt[i] = UInt8((i-1 + 32) & 0xff)
	end
	for (i, saltedVector) in enumerate(saltedVectors)
		c = Blowfish.NewSaltedCipher(key, salt[1:i-1])
		buf = zeros(UInt8, 8)
		Blowfish.Encrypt!(c, buf, buf)
		@test buf == saltedVector
	end
end


text:
#=
smoother.jl

    Kalman smoother routines for a linear Gaussian State Space model

@author: Quint Wiersma <q.wiersma@vu.nl>

@date: 2020/07/29
=#

# Struct
struct Smoother{Tα, TV, TL, Tr, TN} <: KalmanSmoother
	α::Tα	# smoothed state
	V::TV	# smoothed state covariances
	L::TL	# innovation error smoother transition matrix
	r::Tr	# backward smoothing recursion (state)
	N::TN	# backward smoothing recursion (variance)
end
# Constructor
function Smoother(p::Integer, h::Integer, T_len::Integer, T::Type)
    # filter output
    α= Matrix{T}(undef, p, T_len)
    V= Array{T,4}(undef, p, p, h+1, T_len)
    L= Matrix{T}(undef, p, p)
    r= Vector{T}(undef, p)
    N= Matrix{T}(undef, p, p)

    return Smoother(α, V, L, r, N)
end

"""
	computeL!(L_t, K_t, Z, T)

Compute the ``L`` matrix at time ``t``, storing the result in `L_t`.

#### Arguments
  - `K_t::AbstractMatrix`	: Kalman gain (p x n)
  - `Z::AbstractMatrix`		: system matrix ``Z`` (n x p)
  - `T::AbstractMatrix`		: system matrix ``T`` (p x p)

#### Returns
  - `L_t::AbstractMatrix`	: ``Lₜ`` (p x p)
"""
function computeL!(L_t::AbstractMatrix, K_t::AbstractMatrix, Z::AbstractMatrix, T::AbstractMatrix)
	# -KₜｘZ
	mul!(L_t, K_t, Z, -1., .0)
	# Lₜ = T - KₜｘZ
	L_t.+= T

	return nothing
end

function computeL!(L_t::AbstractMatrix, K_t::AbstractMatrix, Z::AbstractMatrix, T::Diagonal)
	# -KₜｘZ
	mul!(L_t, K_t, Z, -1., .0)
	# Lₜ = T - KₜｘZ
	@inbounds @fastmath for i in axes(T,1)
		L_t[i,i]+= T.diag[i]
	end

	return nothing
end

"""
	computeL_eq!(L_it, K_it, Z_i)

Compute the ``L`` matrix at time ``t`` for series ``i``, storing the result in
`L_it`.

#### Arguments
  - `K_it::AbstractVector`	: Kalman gain for series ``i``(p x 1)
  - `Z_i::AbstractVector`	: system vector ``Z`` for series ``i`` (p x 1)

#### Returns
  - `L_it::AbstractMatrix`	: ``Lₜᵢ`` (p x p)
"""
function computeL_eq!(L_it::AbstractMatrix, K_it::AbstractVector, Z_i::AbstractVector)
	# -KₜᵢｘZᵢ'
	L_it.= -K_it .* transpose(Z_i)
	# Lₜᵢ = I - KₜᵢｘZᵢ'
	@inbounds @fastmath for i in axes(L_it,1)
		L_it[i,i]+= 1.
	end

	return nothing
end

"""
	back_state!(r, L_t, fac, v_t, Z, tmp_n, tmp_p)

Compute the backward recursion for the smoothed state ``r`` at time ``t`` using
the inverse of ``F``, storing the result in `r`.

#### Arguments
  - `L_t::AbstractMatrix`	: ``Lₜ`` (p x p)
  - `Fi_t::AbstractMatrix`	: inverse of ``Fₜ`` (n x n)
  - `v_t::AbstractVector`	: forecast error (n x 1)
  - `Z::AbstractMatrix`		: system matrix ``Z`` (n x p)
  - `tmp_n::AbstractVector`	: tmp storage array (n x 1)
  - `tmp_p::AbstractVector`	: tmp storage array (p x 1)

#### Returns
  - `r::AbstractVector`	: Backward recursion for the smoothed state (p x 1)
"""
function back_state!(r::AbstractVector, L_t::AbstractMatrix, Fi_t::AbstractMatrix,
					v_t::AbstractVector, Z::AbstractMatrix, tmp_n::AbstractVector,
					tmp_p::AbstractVector)
	# Lₜ'ｘrₜ
	mul!(tmp_p, transpose(L_t), r)
	# Fₜ⁻¹ｘvₜ
	mul!(tmp_n, Fi_t, v_t)
	# Z'ｘFₜ⁻¹ｘvₜ
	mul!(r, transpose(Z), tmp_n)
	# rₜ₋₁ = Z'ｘFₜ⁻¹ｘvₜ + Lₜ'ｘrₜ
	r.+= tmp_p

	return nothing
end

"""
	back_state!(r, L_t, fac, v_t, Z, tmp_n, tmp_p)

Compute the backward recursion for the smoothed state ``r`` at time ``t`` using
the factorization of ``F``, storing the result in `r`.
"""
function back_state!(r::AbstractVector, L_t::AbstractMatrix, fac::Factorization,
					v_t::AbstractVector, Z::AbstractMatrix, tmp_n::AbstractVector,
					tmp_p::AbstractVector)
	# Lₜ'ｘrₜ
	mul!(tmp_p, transpose(L_t), r)
	# Fₜ⁻¹ｘvₜ
	ldiv!(tmp_n, fac, v_t)
	# Z'ｘFₜ⁻¹ｘvₜ
	mul!(r, transpose(Z), tmp_n)
	# rₜ₋₁ = Z'ｘFₜ⁻¹ｘvₜ + Lₜ'ｘrₜ
	r.+= tmp_p

	return nothing
end

"""
	back_state_eq!(r, L_i, F, v, Z_i, tmp_p)

Compute the backward recursion for the smoothed state ``r`` at time ``t`` for
series ``i-1`` using univariate treatment, storing the result in `r`.

#### Arguments
  - `L_i::AbstractMatrix`	: ``Lₜᵢ`` (p x p)
  - `F::Real`				: forecst error variance for series ``i``
  - `v::Real`				: forecast error for series ``i``
  - `Z_i::AbstractVector`	: system vector ``Z`` for series ``i`` (p x 1)
  - `tmp_p::AbstractVector`	: tmp storage array (p x 1)

#### Returns
  - `r::AbstractVector`	: Backward recursion for the smoothed state (p x 1)
"""
function back_state_eq!(r::AbstractVector, L_i::AbstractMatrix, F::Real, v::Real,
						Z_i::AbstractVector, tmp_p::AbstractVector)
	# Lᵢ'ｘrᵢ
	mul!(tmp_p, transpose(L_i), r)
	# rᵢ₋₁ = Zᵢ'ｘF⁻¹ｘv + Lₜ'ｘrᵢ
	r.= inv(F) .* v .* Z_i .+ tmp_p

	return nothing
end

"""
	back_state_var!(N, L_t, Fi_t, Z, tmp_pn, tmp_p)

Compute the backward recursion for the smoothed state variance ``N`` at time
``t`` using the inverse of ``F``, storing the result in `N`.

#### Arguments
  - `L_t::AbstractMatrix`	: ``Lₜ`` (p x p)
  - `Fi_t::AbstractMatrix`	: inverse of ``Fₜ`` (n x n)
  - `Z::AbstractMatrix`		: system matrix ``Z`` (n x p)
  - `tmp_np::AbstractMatrix`: tmp storage array (n x p)
  - `tmp_pp::AbstractMatrix`: tmp storage array (p x p)

#### Returns
  - `N::AbstractMatrix`	: Backward recursion for the smoothed state variance (p x p)
"""
function back_state_var!(N::AbstractMatrix, L_t::AbstractMatrix, Fi_t::AbstractMatrix,
						Z::AbstractMatrix, tmp_np::AbstractMatrix, tmp_pp::AbstractMatrix)
	# NₜｘLₜ
	mul!(tmp_pp, N, L_t)
	# LₜｘNₜｘLₜ
	mul!(N, transpose(L_t), tmp_pp)
	# Fₜ⁻¹ｘZ
	mul!(tmp_np, Fi_t, Z)
	# Nₜ₋₁ = Z'ｘFₜ⁻¹ｘZ + LₜｘNₜｘLₜ
	mul!(N, transpose(Z), tmp_np, 1., 1.)

	return nothing
end

"""
	back_state_var!(N, L_t, fac, Z, tmp_np, tmp_pp)

Compute the backward recursion for the smoothed state variance ``N`` at time
``t`` using the factorization of ``F``, storing the result in `N`.
"""
function back_state_var!(N::AbstractMatrix, L_t::AbstractMatrix, fac::Factorization,
						Z::AbstractMatrix, tmp_np::AbstractMatrix, tmp_pp::AbstractMatrix)
	# NₜｘLₜ
	mul!(tmp_pp, N, L_t)
	# LₜｘNₜｘLₜ
	mul!(N, transpose(L_t), tmp_pp)
	# Fₜ⁻¹ｘZ
	ldiv!(tmp_np, fac, Z)
	# Nₜ₋₁ = Z'ｘFₜ⁻¹ｘZ + LₜｘNₜｘLₜ
	mul!(N, transpose(Z), tmp_np, 1., 1.)

	return nothing
end

"""
	back_state_var_eq!(N, L_i, F, Z_i, tmp_p)

Compute the backward recursion for the smoothed state variance ``N`` at time
``t`` for series ``i-1`` using the univariate treatment, storing the result in
`N`.

#### Arguments
  - `L_i::AbstractMatrix`	: ``Lₜᵢ`` (p x p)
  - `F::Real`				: forecst error variance for series ``i``
  - `Z_i::AbstractVector`	: system vector ``Z`` for series ``i`` (p x 1)
  - `tmp_pp::AbstractMatrix`: tmp storage array (p x p)

#### Returns
  - `N::AbstractMatrix`	: Backward recursion for the smoothed state variance (p x p)
"""
function back_state_var_eq!(N::AbstractMatrix, L_i::AbstractMatrix, F::Real,
						Z_i::AbstractVector, tmp_pp::AbstractMatrix)
	# NᵢｘLᵢ
	mul!(tmp_pp, N, L_i)
	# LᵢｘNᵢｘLᵢ
	mul!(N, transpose(L_i), tmp_pp)
	# Nᵢ₋₁ = Zᵢ'ｘF⁻¹ｘZᵢ + LᵢｘNᵢｘLᵢ
	N.+= inv(F) .* Z_i .* transpose(Z_i)

	return nothing
end

"""
	transition!(X, T, tT, tmp)

Compute the linear time transition for smoother recursion objects from time
``t`` to ``t-1`` using the equation-by-equation or univariate version of the
Kalman smoother. Storing the results in `X`.

#### Arguments
  - `T::AbstractMatrix`		: system matrix ``T`` (p x p)
  - `tT::AbstractChar`		: transpose indicator
  - `tmp::AbstractVecOrMat`	: tmp storage array (p x 1) or (p x p)

#### Returns
  - `X::AbstractVecOrMat`	: smoother object (p x 1) or (p x p)
"""
function transition!(X::AbstractVecOrMat, T::AbstractMatrix, tT::AbstractChar, tmp::AbstractVecOrMat)
	# Store
	copyto!(tmp, X)

	if tT == 'T'
		# Xₜ₋₁ = T'ｘXₜ
		mul!(X, transpose(T), tmp)
	elseif tT == 'N'
		# Xₜ₋₁ = TｘXₜ
		mul!(X, T, tmp)
	end

	return nothing
end

function transition!(X::AbstractVecOrMat, T::Diagonal, tT::AbstractChar, tmp::AbstractVecOrMat)
	# Xₜ₋₁ = TｘXₜ
	X.= T.diag .* X

	return nothing
end

"""
	transition2!(X, T, tmp)

Compute the quadratic time transition for smoother recursion objects from time
``t`` to ``t-1`` using the equation-by-equation or univariate version of the
Kalman smoother. Storing the results in `X`.

#### Arguments
  - `T::AbstractMatrix`		: system matrix ``T`` (p x p)
  - `tmp::AbstractMatrix`	: tmp storage array (p x p)

#### Returns
  - `X::AbstractMatrix`		: smoother object (p x p)
"""
function transition2!(X::AbstractMatrix, T::AbstractMatrix, tmp::AbstractMatrix)
	# T'ｘXₜ
	mul!(tmp, transpose(T), X)
	# Xₜ₋₁ = T'ｘXₜｘT
	mul!(X, tmp, T)

	return nothing
end

function transition2!(X::AbstractMatrix, T::Diagonal, tmp::AbstractMatrix)
	# Xₜ₋₁ = T'ｘXₜｘT
	X.= T.diag .* X .* transpose(T.diag)

	return nothing
end

"""
	smooth_state!(α_t, a_t, P_t, r)

Compute the smoothed state ``α`` at time ``t``, storing the result in `α_t`.

#### Arguments
  - `a_t::AbstractVector`	: filtered state at time ``t`` (p x 1)
  - `P_t::AbstractMatrix`	: filtered state variance at time ``t`` (p x p)
  - `r::AbstractVector`		: Backward recursion for states (p x 1)

#### Returns
  - `α_t::AbstractVector`	: smoothed state at time ``t`` (p x 1)
"""
function smooth_state!(α_t::AbstractVector, a_t::AbstractVector, P_t::AbstractMatrix, r::AbstractVector)
	# Pₜｘrₜ₋₁
	mul!(α_t, P_t, r)
	# αₜ = aₜ + Pₜｘrₜ₋₁
	α_t.+= a_t

	return nothing
end

"""
	smooth_state_var!(V_t, P_t, N, tmp)

Compute the smoothed state variance ``V`` at time ``t``, storing the result in
`V_t`.

#### Arguments
  - `P_t::AbstractMatrix`	: filtered state variance at time ``t`` (p x p)
  - `N::AbstractMatrix`		: Backward recursion for states variance (p x p)
  - `tmp::AbstractMatrix`	: tmp storage array (p x p)

#### Returns
  - `V_t::AbstractMatrix`	: smoothed state variance at time ``t`` (p x p)
"""
function smooth_state_var!(V_t::AbstractMatrix, P_t::AbstractMatrix, N::AbstractMatrix, tmp::AbstractMatrix)
	# Nₜ₋₁ｘPₜ
	mul!(tmp, N, P_t)
	# -PₜｘNₜ₋₁ｘPₜ
	mul!(V_t, P_t, tmp, -1., .0)
	# Vₜ = Pₜ - PₜｘNₜ₋₁ｘPₜ
	V_t.+= P_t

	return nothing
end

"""
	smooth_state_cov!(V, N, L, P_t, t, h, tmp_0, tmp_1)

Compute the smoothed state variance and autocovariances ``V(j)`` at lags 0 to
``h`` at time ``t``, storing the result in `V`.

#### Arguments
  - `N::AbstractMatrix`		: Backward recursion for states variance (p x p)
  - `L::AbstractMatrix`		: ``L`` (p x p)
  - `P_t::AbstractMatrix`	: filtered state variance at time ``t`` (p x p)
  - `t::Integer`			: time index
  - `h::Integer`			: lag length
  - `tmp_0::AbstractMatrix`	: tmp storage array (p x p)
  - `tmp_1::AbstractMatrix`	: tmp storage array (p x p)

#### Returns
  - `V::AbstractArray`	: smoothed state autocovariances at time ``t`` (p x p x h+1 x T)
"""
function smooth_state_cov!(V::AbstractArray, N::AbstractMatrix, L::AbstractMatrix,
							P_t::AbstractMatrix, t::Integer, h::Integer,
							tmp_0::AbstractMatrix, tmp_1::AbstractMatrix)
	# Get dims
	T_len= size(V, 4)

	# -PₜｘNₜ₋₁
	mul!(tmp_0, P_t, N, -1., .0)
	# I - PₜｘNₜ₋₁
	@inbounds @fastmath for i in axes(P_t,1)
		tmp_0[i,i]+= 1.
	end

	# Loop through lags
	@inbounds @fastmath for j in 0:h
		# Store
		V[:,:,j+1,t].= tmp_0

		for s in t+1:min(t+j,T_len)
			# Store view
			V_js= view(V,:,:,j+1,s)

			# Propagate backwards with Lₜ₋₁
			mul!(tmp_1, V_js, L)
			V_js.= tmp_1
		end

		# Complete autocovariance
		if t+j <= T_len
			# Store view
			V_jt= view(V,:,:,j+1,t+j)

			# Vₜ(j) = (I - Pₜ₊ⱼｘNₜ₊ⱼ₋₁)ｘLₜ₊ⱼ₋₁ｘ...ｘLₜ₊₁ｘLₜｘPₜ
			mul!(tmp_1, V_jt, P_t)
			V_jt.= tmp_1
		end
	end

	return nothing
end

"""
	kalman_smoother!(smoother, filter, sys)

Compute smoothed states ``α`` and corresponding variances ``V`` for a linear
Gaussian State Space model with system matrices `sys` and Kalman filter output
`filter`, storing the result in `smoother`.

#### Arguments
  - `filter::KalmanFilter`	: Kalman filter output
  - `sys::StateSpaceSystem`	: state space system matrices

#### Returns
  - `smoother::Smoother`	: Kalman smoother output
"""
function kalman_smoother!(smoother::Smoother, filter::MultivariateFilter, sys::LinearTimeInvariant)
	# Get dims
	(p,n,T_len)= size(filter.K)

    # Initialize temp. containers
    tmp_n= Vector{Float64}(undef, n)
    tmp_p= Vector{Float64}(undef, p)
	tmp_nn= Matrix{Float64}(undef, (n,n))
	tmp_np= Matrix{Float64}(undef, (n,p))
	tmp_pp= Matrix{Float64}(undef, (p,p))

	# Initialize smoother
	smoother.r.= zero(Float64)
	smoother.N.= zero(Float64)

	# Smoother
	@inbounds @fastmath for t in T_len:-1:1
        # Store views
		α_t= view(smoother.α,:,t)
		V_t= view(smoother.V,:,:,t)
        a_t= view(filter.a,:,t)
		P_t= view(filter.P,:,:,t)
        v_t= view(filter.v,:,t)
        F_t= view(filter.F,:,:,t)
        K_t= view(filter.K,:,:,t)

		# Lₜ
		computeL!(smoother.L, K_t, sys.Z, sys.T)

		# Cholesky factorization of Fₜ
		copyto!(tmp_nn, F_t)
		fac= cholesky!(Hermitian(tmp_nn))

		# Backward recursion state
		back_state!(smoother.r, smoother.L, fac, v_t, sys.Z, tmp_n, tmp_p)

		# Backward recursion state variance
		back_state_var!(smoother.N, smoother.L, fac, sys.Z, tmp_np, tmp_pp)

		# Smoothed states
		smooth_state!(α_t, a_t, P_t, smoother.r)

		# Smoothed states variance
		smooth_state_var!(V_t, P_t, smoother.N, tmp_pp)
	end

	return nothing
end

function kalman_smoother!(smoother::Smoother, filter::MultivariateFilter, sys::LinearTimeVariant)
	# Get dims
	(p,n,T_len)= size(filter.K)

    # Initialize temp. containers
    tmp_n= Vector{Float64}(undef, n)
    tmp_p= Vector{Float64}(undef, p)
	tmp_nn= Matrix{Float64}(undef, (n,n))
	tmp_np= Matrix{Float64}(undef, (n,p))
	tmp_pp= Matrix{Float64}(undef, (p,p))

	# Initialize smoother
	smoother.r.= zero(Float64)
	smoother.N.= zero(Float64)

	# Smoother
	@inbounds @fastmath for t in T_len:-1:1
        # Store views
		α_t= view(smoother.α,:,t)
		V_t= view(smoother.V,:,:,t)
        a_t= view(filter.a,:,t)
		P_t= view(filter.P,:,:,t)
        v_t= view(filter.v,:,t)
        F_t= view(filter.F,:,:,t)
        K_t= view(filter.K,:,:,t)

		# Lₜ
		computeL!(smoother.L, K_t, sys.Z[t], sys.T[t])

		# Cholesky factorization of Fₜ
		copyto!(tmp_nn, F_t)
		fac= cholesky!(Hermitian(tmp_nn))

		# Backward recursion state
		back_state!(smoother.r, smoother.L, fac, v_t, sys.Z[t], tmp_n, tmp_p)

		# Backward recursion state variance
		back_state_var!(smoother.N, smoother.L, fac, sys.Z[t], tmp_np, tmp_pp)

		# Smoothed states
		smooth_state!(α_t, a_t, P_t, smoother.r)

		# Smoothed states variance
		smooth_state_var!(V_t, P_t, smoother.N, tmp_pp)
	end

	return nothing
end

function kalman_smoother!(smoother::Smoother, filter::WoodburyFilter, sys::LinearTimeInvariant)
	# Get dims
	(p,n,T_len)= size(filter.K)

    # Initialize temp. containers
    tmp_n= Vector{Float64}(undef, n)
    tmp_p= Vector{Float64}(undef, p)
	tmp_np= Matrix{Float64}(undef, (n,p))
	tmp_pp= Matrix{Float64}(undef, (p,p))

	# Initialize smoother
	smoother.r.= zero(Float64)
	smoother.N.= zero(Float64)

	# Smoother
	@inbounds @fastmath for t in T_len:-1:1
        # Store views
		α_t= view(smoother.α,:,t)
		V_t= view(smoother.V,:,:,t)
        a_t= view(filter.a,:,t)
		P_t= view(filter.P,:,:,t)
        v_t= view(filter.v,:,t)
        Fi_t= view(filter.Fi,:,:,t)
        K_t= view(filter.K,:,:,t)

		# Lₜ
		computeL!(smoother.L, K_t, sys.Z, sys.T)

		# Backward recursion state
		back_state!(smoother.r, smoother.L, Fi_t, v_t, sys.Z, tmp_n, tmp_p)

		# Backward recursion state variance
		back_state_var!(smoother.N, smoother.L, Fi_t, sys.Z, tmp_np, tmp_pp)

		# Smoothed states
		smooth_state!(α_t, a_t, P_t, smoother.r)

		# Smoothed states variance
		smooth_state_var!(V_t, P_t, smoother.N, tmp_pp)
	end

	return nothing
end

function kalman_smoother!(smoother::Smoother, filter::WoodburyFilter, sys::LinearTimeVariant)
	# Get dims
	(p,n,T_len)= size(filter.K)

    # Initialize temp. containers
    tmp_n= Vector{Float64}(undef, n)
    tmp_p= Vector{Float64}(undef, p)
	tmp_np= Matrix{Float64}(undef, (n,p))
	tmp_pp= Matrix{Float64}(undef, (p,p))

	# Initialize smoother
	smoother.r.= zero(Float64)
	smoother.N.= zero(Float64)

	# Smoother
	@inbounds @fastmath for t in T_len:-1:1
        # Store views
		α_t= view(smoother.α,:,t)
		V_t= view(smoother.V,:,:,t)
        a_t= view(filter.a,:,t)
		P_t= view(filter.P,:,:,t)
        v_t= view(filter.v,:,t)
        Fi_t= view(filter.Fi,:,:,t)
        K_t= view(filter.K,:,:,t)

		# Lₜ
		computeL!(smoother.L, K_t, sys.Z[t], sys.T[t])

		# Backward recursion state
		back_state!(smoother.r, smoother.L, Fi_t, v_t, sys.Z[t], tmp_n, tmp_p)

		# Backward recursion state variance
		back_state_var!(smoother.N, smoother.L, Fi_t, sys.Z[t], tmp_np, tmp_pp)

		# Smoothed states
		smooth_state!(α_t, a_t, P_t, smoother.r)

		# Smoothed states variance
		smooth_state_var!(V_t, P_t, smoother.N, tmp_pp)
	end

	return nothing
end

function kalman_smoother!(smoother::Smoother, filter::UnivariateFilter, sys::LinearTimeInvariant)
	# Get dims
	(p,n,T_len)= size(filter.K)

    # Initialize temp. containers
    tmp_p= Vector{Float64}(undef, p)
	tmp_pp= Matrix{Float64}(undef, (p,p))

	# Initialize smoother
	smoother.r.= zero(Float64)
	smoother.N.= zero(Float64)

	# Tranpose
	Zt= transpose(sys.Z)

	# Smoother
	@inbounds @fastmath for	t in T_len:-1:1
		# Store views
		α_t= view(smoother.α,:,t)
		V_t= view(smoother.V,:,:,t)
        a_t= view(filter.a,:,1,t)
		P_t= view(filter.P,:,:,1,t)

		# Time transition
		transition!(smoother.r, sys.T, 'T', tmp_p)
		transition2!(smoother.N, sys.T, tmp_pp)

		for i in n:-1:1
			# Store views
	        K_it= view(filter.K,:,i,t)
			Z_i= view(Zt,:,i)

			# Lₜᵢ
			computeL_eq!(smoother.L, K_it, Z_i)

			# Backward recursion state
			back_state_eq!(smoother.r, smoother.L, filter.F[i,t], filter.v[i,t], Z_i, tmp_p)

			# Backward recursion state variance
			back_state_var_eq!(smoother.N, smoother.L, filter.F[i,t], Z_i, tmp_pp)
		end

		# Smoothed states
		smooth_state!(α_t, a_t, P_t, smoother.r)

		# Smoothed states variance
		smooth_state_var!(V_t, P_t, smoother.N, tmp_pp)
	end

	return nothing
end

function kalman_smoother!(smoother::Smoother, filter::UnivariateFilter, sys::LinearTimeVariant)
	# Get dims
	(p,n,T_len)= size(filter.K)

    # Initialize temp. containers
    tmp_p= Vector{Float64}(undef, p)
	tmp_pp= Matrix{Float64}(undef, (p,p))

	# Initialize smoother
	smoother.r.= zero(Float64)
	smoother.N.= zero(Float64)

	# Smoother
	@inbounds @fastmath for	t in T_len:-1:1
		# store
		Z_t= sys.Z[t]

		# Store views
		α_t= view(smoother.α,:,t)
		V_t= view(smoother.V,:,:,t)
        a_t= view(filter.a,:,1,t)
		P_t= view(filter.P,:,:,1,t)

		# Time transition
		transition!(smoother.r, sys.T[t], 'T', tmp_p)
		transition2!(smoother.N, sys.T[t], tmp_pp)

		for i in n:-1:1
			# Store views
	        K_it= view(filter.K,:,i,t)
			Z_it= view(Z_t,i,:)

			# Lₜᵢ
			computeL_eq!(smoother.L, K_it, Z_it)

			# Backward recursion state
			back_state_eq!(smoother.r, smoother.L, filter.F[i,t], filter.v[i,t], Z_it, tmp_p)

			# Backward recursion state variance
			back_state_var_eq!(smoother.N, smoother.L, filter.F[i,t], Z_it, tmp_pp)
		end

		# Smoothed states
		smooth_state!(α_t, a_t, P_t, smoother.r)

		# Smoothed states variance
		smooth_state_var!(V_t, P_t, smoother.N, tmp_pp)
	end

	return nothing
end

"""
	kalman_smoother_cov!(smoother, filter, sys; h=1)

Compute smoothed states ``α`` and corresponding variances ``V`` and
autocovariances ``V(h)`` up until lag `h` for a linear Gaussian State Space
model with system matrices `sys` and Kalman filter output `filter`, storing the
results in `smoother`.

#### Arguments
  - `filter::KalmanFilter`	: Kalman filter output
  - `sys::StateSpaceSystem`	: state space system matrices
  - `h::Integer`			: lag length

#### Returns
  - `smoother::Smoother`	: Kalman smoother output
"""
function kalman_smoother_cov!(smoother::Smoother, filter::MultivariateFilter,
								sys::LinearTimeInvariant; h::Integer=1)
	# Get dims
	(p,n,T_len)= size(filter.K)

    # Initialize temp. containers
    tmp_n= Vector{Float64}(undef, n)
    tmp_p= Vector{Float64}(undef, p)
	tmp_nn= Matrix{Float64}(undef, (n,n))
	tmp_np= Matrix{Float64}(undef, (n,p))
	tmp_pp= Matrix{Float64}(undef, (p,p))
	tmp1_pp= Matrix{Float64}(undef, (p,p))

	# Initialize smoother
	smoother.r.= zero(Float64)
	smoother.N.= zero(Float64)

	# Smoother
	@inbounds @fastmath for t in T_len:-1:1
        # Store views
		α_t= view(smoother.α,:,t)
        a_t= view(filter.a,:,t)
		P_t= view(filter.P,:,:,t)
        v_t= view(filter.v,:,t)
        F_t= view(filter.F,:,:,t)
        K_t= view(filter.K,:,:,t)

		# Lₜ
		computeL!(smoother.L, K_t, sys.Z, sys.T)

		# Cholesky factorization of Fₜ
		copyto!(tmp_nn, F_t)
		fac= cholesky!(Hermitian(tmp_nn))

		# Backward recursion state
		back_state!(smoother.r, smoother.L, fac, v_t, sys.Z, tmp_n, tmp_p)

		# Backward recursion state variance
		back_state_var!(smoother.N, smoother.L, fac, sys.Z, tmp_np, tmp_pp)

		# Smoothed states
		smooth_state!(α_t, a_t, P_t, smoother.r)

		# Smoothed states autocovariance
		smooth_state_cov!(smoother.V, smoother.N, smoother.L, P_t, t, h, tmp_pp, tmp1_pp)
	end

	return nothing
end

function kalman_smoother_cov!(smoother::Smoother, filter::MultivariateFilter,
								sys::LinearTimeVariant; h::Integer=1)
	# Get dims
	(p,n,T_len)= size(filter.K)

    # Initialize temp. containers
    tmp_n= Vector{Float64}(undef, n)
    tmp_p= Vector{Float64}(undef, p)
	tmp_nn= Matrix{Float64}(undef, (n,n))
	tmp_np= Matrix{Float64}(undef, (n,p))
	tmp_pp= Matrix{Float64}(undef, (p,p))
	tmp1_pp= Matrix{Float64}(undef, (p,p))

	# Initialize smoother
	smoother.r.= zero(Float64)
	smoother.N.= zero(Float64)

	# Smoother
	@inbounds @fastmath for t in T_len:-1:1
        # Store views
		α_t= view(smoother.α,:,t)
        a_t= view(filter.a,:,t)
		P_t= view(filter.P,:,:,t)
        v_t= view(filter.v,:,t)
        F_t= view(filter.F,:,:,t)
        K_t= view(filter.K,:,:,t)

		# Lₜ
		computeL!(smoother.L, K_t, sys.Z[t], sys.T[t])

		# Cholesky factorization of Fₜ
		copyto!(tmp_nn, F_t)
		fac= cholesky!(Hermitian(tmp_nn))

		# Backward recursion state
		back_state!(smoother.r, smoother.L, fac, v_t, sys.Z[t], tmp_n, tmp_p)

		# Backward recursion state variance
		back_state_var!(smoother.N, smoother.L, fac, sys.Z[t], tmp_np, tmp_pp)

		# Smoothed states
		smooth_state!(α_t, a_t, P_t, smoother.r)

		# Smoothed states autocovariance
		smooth_state_cov!(smoother.V, smoother.N, smoother.L, P_t, t, h, tmp_pp, tmp1_pp)
	end

	return nothing
end

function kalman_smoother_cov!(smoother::Smoother, filter::WoodburyFilter,
								sys::LinearTimeInvariant; h::Integer=1)
	# Get dims
	(p,n,T_len)= size(filter.K)

    # Initialize temp. containers
    tmp_n= Vector{Float64}(undef, n)
    tmp_p= Vector{Float64}(undef, p)
	tmp_np= Matrix{Float64}(undef, (n,p))
	tmp_pp= Matrix{Float64}(undef, (p,p))
	tmp1_pp= Matrix{Float64}(undef, (p,p))

	# Initialize smoother
	smoother.r.= zero(Float64)
	smoother.N.= zero(Float64)

	# Smoother
	@inbounds @fastmath for t in T_len:-1:1
        # Store views
		α_t= view(smoother.α,:,t)
        a_t= view(filter.a,:,t)
		P_t= view(filter.P,:,:,t)
        v_t= view(filter.v,:,t)
        Fi_t= view(filter.Fi,:,:,t)
        K_t= view(filter.K,:,:,t)

		# Lₜ
		computeL!(smoother.L, K_t, sys.Z, sys.T)

		# Backward recursion state
		back_state!(smoother.r, smoother.L, Fi_t, v_t, sys.Z, tmp_n, tmp_p)

		# Backward recursion state variance
		back_state_var!(smoother.N, smoother.L, Fi_t, sys.Z, tmp_np, tmp_pp)

		# Smoothed states
		smooth_state!(α_t, a_t, P_t, smoother.r)

		# Smoothed states autocovariance
		smooth_state_cov!(smoother.V, smoother.N, smoother.L, P_t, t, h, tmp_pp, tmp1_pp)
	end

	return nothing
end

function kalman_smoother_cov!(smoother::Smoother, filter::WoodburyFilter,
								sys::LinearTimeVariant; h::Integer=1)
	# Get dims
	(p,n,T_len)= size(filter.K)

    # Initialize temp. containers
    tmp_n= Vector{Float64}(undef, n)
    tmp_p= Vector{Float64}(undef, p)
	tmp_np= Matrix{Float64}(undef, (n,p))
	tmp_pp= Matrix{Float64}(undef, (p,p))
	tmp1_pp= Matrix{Float64}(undef, (p,p))

	# Initialize smoother
	smoother.r.= zero(Float64)
	smoother.N.= zero(Float64)

	# Smoother
	@inbounds @fastmath for t in T_len:-1:1
        # Store views
		α_t= view(smoother.α,:,t)
        a_t= view(filter.a,:,t)
		P_t= view(filter.P,:,:,t)
        v_t= view(filter.v,:,t)
        Fi_t= view(filter.Fi,:,:,t)
        K_t= view(filter.K,:,:,t)

		# Lₜ
		computeL!(smoother.L, K_t, sys.Z[t], sys.T[t])

		# Backward recursion state
		back_state!(smoother.r, smoother.L, Fi_t, v_t, sys.Z[t], tmp_n, tmp_p)

		# Backward recursion state variance
		back_state_var!(smoother.N, smoother.L, Fi_t, sys.Z[t], tmp_np, tmp_pp)

		# Smoothed states
		smooth_state!(α_t, a_t, P_t, smoother.r)

		# Smoothed states autocovariance
		smooth_state_cov!(smoother.V, smoother.N, smoother.L, P_t, t, h, tmp_pp, tmp1_pp)
	end

	return nothing
end

function kalman_smoother_cov!(smoother::Smoother, filter::UnivariateFilter,
								sys::LinearTimeInvariant; h::Integer=1)
	# Get dims
	(p,n,T_len)= size(filter.K)

    # Initialize temp. containers
	L_i= Matrix{Float64}(undef, (p,p))
    tmp_p= Vector{Float64}(undef, p)
	tmp_pp= Matrix{Float64}(undef, (p,p))
	tmp1_pp= Matrix{Float64}(undef, (p,p))

	# Initialize smoother
	smoother.r.= zero(Float64)
	smoother.N.= zero(Float64)

	# Smoother
	@inbounds @fastmath for	t in T_len:-1:1
		# Store views
		α_t= view(smoother.α,:,t)
        a_t= view(filter.a,:,1,t)
		P_t= view(filter.P,:,:,1,t)

		# reset
		smoother.L.= zero(Float64)
		for i in 1:p
			smoother.L[i,i]= one(Float64)
		end

		# Time transitions
		transition!(smoother.r, sys.T, 'T', tmp_p)
		transition2!(smoother.N, sys.T, tmp_pp)

		for i in n:-1:1
			# Store views
	        K_it= view(filter.K,:,i,t)
			Z_i= view(sys.Z,i,:)

			# Lₜᵢ
			computeL_eq!(L_i, K_it, Z_i)

			# Lₜᵢｘ...ｘLₜ₁
			copyto!(tmp_pp, smoother.L)
			mul!(smoother.L, tmp_pp, L_i)

			# Backward recursion state
			back_state_eq!(smoother.r, L_i, filter.F[i,t], filter.v[i,t], Z_i, tmp_p)

			# Backward recursion state variance
			back_state_var_eq!(smoother.N, L_i, filter.F[i,t], Z_i, tmp_pp)
		end

		# Lₜ = TｘLₜₙｘ...ｘLₜ₁
		transition!(smoother.L, sys.T, 'N', tmp_pp)

		# Smoothed states
		smooth_state!(α_t, a_t, P_t, smoother.r)

		# Smoothed states autocovariance
		smooth_state_cov!(smoother.V, smoother.N, smoother.L, P_t, t, h, tmp_pp, tmp1_pp)
	end

	return nothing
end

function kalman_smoother_cov!(smoother::Smoother, filter::UnivariateFilter,
								sys::LinearTimeVariant; h::Integer=1)
	# Get dims
	(p,n,T_len)= size(filter.K)

    # Initialize temp. containers
	L_i= Matrix{Float64}(undef, (p,p))
    tmp_p= Vector{Float64}(undef, p)
	tmp_pp= Matrix{Float64}(undef, (p,p))
	tmp1_pp= Matrix{Float64}(undef, (p,p))

	# Initialize smoother
	smoother.r.= zero(Float64)
	smoother.N.= zero(Float64)

	# Smoother
	@inbounds @fastmath for	t in T_len:-1:1
		# store
		Z_t= sys.Z[t]

		# Store views
		α_t= view(smoother.α,:,t)
        a_t= view(filter.a,:,1,t)
		P_t= view(filter.P,:,:,1,t)

		# reset
		smoother.L.= zero(Float64)
		for i in 1:p
			smoother.L[i,i]= one(Float64)
		end

		# Time transitions
		transition!(smoother.r, sys.T[t], 'T', tmp_p)
		transition2!(smoother.N, sys.T[t], tmp_pp)

		for i in n:-1:1
			# Store views
	        K_it= view(filter.K,:,i,t)
			Z_it= view(Z_t,i,:)

			# Lₜᵢ
			computeL_eq!(L_i, K_it, Z_it)

			# Lₜᵢｘ...ｘLₜ₁
			copyto!(tmp_pp, smoother.L)
			mul!(smoother.L, tmp_pp, L_i)

			# Backward recursion state
			back_state_eq!(smoother.r, L_i, filter.F[i,t], filter.v[i,t], Z_it, tmp_p)

			# Backward recursion state variance
			back_state_var_eq!(smoother.N, L_i, filter.F[i,t], Z_it, tmp_pp)
		end

		# Lₜ = TｘLₜₙｘ...ｘLₜ₁
		transition!(smoother.L, sys.T[t], 'N', tmp_pp)

		# Smoothed states
		smooth_state!(α_t, a_t, P_t, smoother.r)

		# Smoothed states autocovariance
		smooth_state_cov!(smoother.V, smoother.N, smoother.L, P_t, t, h, tmp_pp, tmp1_pp)
	end

	return nothing
end


text:
#Julia Unit Testing

using FactCheck
f(x) = x^3
facts("cubes") do
    @fact f(2) --> 8
    @fact f(2) --> 7
end

using Base.Test
f(x) = x^3
@test f(2) == 8
@test f(2) == 7

max_stars_repo_name: PacktPublishing/Learning-Jupyter

max_issues_repo_name: PacktPublishing/Learning-Jupyter

max_forks_repo_name: PacktPublishing/Learning-Jupyter


text:
# some views are in colorchannels.jl
using Colors, ImageCore, Base.Test

@testset "rawview" begin
    a = map(U8, rand(3,5))
    v = rawview(a)
    @test v[2,2] === a[2,2].i
    v[1,3] = 0xff
    @test a[1,3] === U8(1)
    v[1,3] = 0x01
    @test a[1,3] === U8(1/255)
end

@testset "permuteddimsview" begin
    a = [1 3; 2 4]
    v = permuteddimsview(a, (1,2))
    @test v == a
    v = permuteddimsview(a, (2,1))
    @test v == a'
    a = rand(3,7,5)
    v = permuteddimsview(a, (2,3,1))
    @test v == permutedims(a, (2,3,1))
end

nothing


text:
using MultiwayNumberPartitioning, HiGHS, JuMP
using DataFrames, UUIDs, StableRNGs, Chain

# In this example, we have a collection of individual animals, identified by UUID, each with a species.
# We wish to obtain partition our animals into 4 groups,
n_partitions = 4
# such that:
# * the groups are of roughly equal size
# * all animals of the same species are in the same group

# First, let's generate our data:
function make_table(n_animals)
    rng = StableRNG(324)
    possible_species = ["Aardvark", "Albatross", "Alligator", "Alpaca", "Anole", "Ant",
                        "Anteater"]
    return DataFrame(:species => [rand(rng, possible_species) for i in 1:n_animals],
                     :id => [uuid4(rng) for i in 1:n_animals])
end

n_animals = 100
df = make_table(n_animals)

# Now, we will group them by species and count how many individuals we have:
count_by_species = combine(groupby(df, :species), nrow => :n_individuals)

# Now, we will partition them.
# First, we choose our optimizer. We will use the new MIT licensed HiGHS solver (https://github.com/ERGO-Code/HiGHS)
optimizer = optimizer_with_attributes(HiGHS.Optimizer, MOI.Silent() => true) # quiet, please!

# Now, perform the partitioning:
transform!(count_by_species,
           :n_individuals => (v -> partition(v, n_partitions; optimizer)) => :partition)

# Already here we can look at `count_by_species` to see the results. However, we can also
# `leftjoin!` these partitions back onto our original DataFrame so that we can see which individual is in which partition:
partitioned_df = leftjoin(df, count_by_species; on=:species)

#####
##### Comparing objective functions
#####

# We can also explore trying different objective functions. Let's start again with our `count_by_species` DataFrame:
count_by_species = combine(groupby(df, :species), nrow => :n_individuals)

# Now, let us try all three algorithms supported by the package:
for alg in (partition_min_largest!, partition_max_smallest!, partition_min_range!)
    alg_name = repr(alg; context=:compact => true)
    transform!(count_by_species,
               :n_individuals => (v -> partition(alg, v, n_partitions; optimizer)) => alg_name)
end
# Let's also generate two random partitions to compare to:
count_by_species.partition_random_1 = rand(StableRNG(431), 1:n_partitions,
                                           nrow(count_by_species))
count_by_species.partition_random_2 = rand(StableRNG(891), 1:n_partitions,
                                           nrow(count_by_species))

# We can summarize the results with a little DataFrames manipulation.
# See <https://bkamins.github.io/julialang/2021/05/28/pivot.html> for more on `stack` and `unstack`.
summary_df = @chain count_by_species begin
    stack(r"partition", :n_individuals; variable_name=:algorithm, value_name=:group)
    groupby([:algorithm, :group])
    combine(:n_individuals => sum => :group_size)
    unstack(:group, :algorithm, :group_size)
end

# We can see the random partitions looks much worse in terms of having groups of equal size; e.g. the first
# has a group of size 12 and one of size 47, while the other partitions have groups of closer sizes.

# We can do some simple checks:
@testset "Partitions add up to the whole" begin
    for col in names(summary_df, r"partition")
        @test sum(summary_df[:, col]) == n_animals
    end
end

@testset "Partitions are optimal (rough check)" begin
    for col in names(summary_df, r"partition")
        # "max smallest" really holds, compared to the other algorithms and to the random partition:
        @test minimum(summary_df[:, :partition_max_smallest!]) >=
              minimum(summary_df[:, col])

        # "min largest" really holds, compared to the other algorithms and to the random partition:
        @test maximum(summary_df[:, :partition_min_largest!]) <= maximum(summary_df[:, col])

        # "min range" really holds:
        r = v -> maximum(v) - minimum(v)
        @test r(summary_df[:, :partition_min_range!]) <= r(summary_df[:, col])
    end
end


text:
module GmshSDK

const PROJECT_ROOT =  pkgdir(GmshSDK)

using Artifacts
using StaticArrays
using RecipesBase
using Printf
import WavePropBase as WPB

const GMSH_VERSION   = "4.8.4"

const GMSH_PATH = artifact"gmsh4.8.4"

const dirs      = readdir(GMSH_PATH,join=true)

@assert length(dirs)==1 "there should be only one directory for the untared gmsh file. Got $dirs"

include(joinpath(first(dirs),"lib","gmsh.jl"))

include("gmshIO.jl")

export
    gmsh,
    @gmsh

end # module

max_stars_repo_name: maltezfaria/GmshSDK

max_issues_repo_name: maltezfaria/GmshSDK

max_forks_repo_name: maltezfaria/GmshSDK


text:
module SynchronicBallot

using PeaceCypher: CypherSuite, Signer, Secret, Layer, secure, id

#abstract type Layer end
#function secure end

using Multiplexers
using Sockets

struct Regulation
    N::UInt8
    M::UInt8
    metadata::Vector{UInt8}
end

struct Ballot
    metadata::Vector{UInt8}
    votes::Array{UInt8,2}
end

abstract type Officer end

function regulation end
function audit! end

##### For Debugging #####

import Base.sync_varname
import Base.@async

macro async(expr)

    tryexpr = quote
        try
            $expr
        catch err
            @warn "error within async" exception=err # line $(__source__.line):
            @show stacktrace(catch_backtrace())
        end
    end

    thunk = esc(:(()->($tryexpr)))

    var = esc(sync_varname)
    quote
        local task = Task($thunk)
        if $(Expr(:isdefined, var))
            push!($var, task)
        end
        schedule(task)
        task
    end
end

const OPEN = UInt8(0)

function stack(io::IO,msg::Vector{UInt8})
    frontbytes = reinterpret(UInt8,Int16[length(msg)])
    item = UInt8[frontbytes...,msg...]
    write(io,item)
end

function unstack(io::IO)
    sizebytes = [read(io,UInt8),read(io,UInt8)]
    size = reinterpret(Int16,sizebytes)[1]

    msg = UInt8[]
    for i in 1:size
        push!(msg,read(io,UInt8))
    end
    return msg
end

####### BallotBox #######

struct Mix
    port
    id
    crypto::CypherSuite
end

function Sockets.connect(mix::Mix)
    layer = Layer(mix.crypto, mix.id)
    mixsocket = connect(mix.port)
    mixsl = secure(mixsocket, layer)
    return mixsl
end

function serve(secureserversocket::IO, ballotmember::Layer)
    N = read(secureserversocket,UInt8)
    M = read(secureserversocket,UInt8)

    mux = Multiplexer(secureserversocket,N)

    susersockets = []
    for i in 1:N
        securesocket = secure(mux.lines[i],ballotmember)
        push!(susersockets,securesocket)
    end

    messages = UInt8[]
    for i in 1:N
        write(susersockets[i],OPEN)

        msg = unstack(susersockets[i])  ### One needs to give a msg type here. ID

        @assert length(msg) == M

        push!(messages,msg...)
    end
    close(mux)

    shapedmessages = reshape(messages,(M,N))

    ### I need to sort whoole collumns
    #sort!(shapedmessages, dims=2) # We may need to look into it
    sortedmessages = sortslices(shapedmessages, dims=2)
    stack(secureserversocket, reshape(sortedmessages,:))
end

function serve(mix::Mix, signer::Union{Signer, Secret})
    @assert mix.id == id(signer)
    layer = Layer(mix.crypto, signer)

    server = listen(mix.port)

    while true
        client = secure(accept(server), layer)

        @async while isopen(client)
            serve(client, layer)
        end
    end
end

struct GateKeeper
    port
    id
    crypto::CypherSuite
    mix::Mix
end

function run(b::Regulation, usersockets::Vector{IO}, secureroutersocket::IO)
    N, M, metadata = b.N, b.M, b.metadata

    write(secureroutersocket,UInt8[N,M])

    mux = Multiplexer(secureroutersocket,usersockets)
    wait(mux)

    messages = unstack(secureroutersocket) ### A single list of bytes.

    ballot = reshape(messages,(M,N))

    signatures = Vector{UInt8}[]

    ### Grouping of data.
    ### Perhaps that could be done with in a Sockets layer which
    ### would only send data if flushed!
    ### The flush may be done automatically when one reads a byte
    io = IOBuffer()
    stack(io,metadata)
    stack(io,messages)
    rawballot = take!(io)

    for i in 1:N
        write(usersockets[i], rawballot)
        s = unstack(usersockets[i]) ### One should verify also the signatures
        push!(signatures, s)
    end

    ### User then can use metadata to form a Braid one needs
    return Ballot(metadata, ballot), signatures
end

function serve!(gk::GateKeeper, ar::Officer, signer::Union{Signer, Secret})

    @assert id(signer) == gk.id

    mix = connect(gk.mix)
    server = listen(gk.port)
    gatemember = Layer(gk.crypto, signer, ar)

    while true

        config = regulation(ar)

        usersockets = IO[]

        while length(usersockets) < config.N
            connection = accept(server)
            secureusersocket = secure(connection, gatemember)
            push!(usersockets, secureusersocket)
            #lock(ar, id(secureusersocket))
        end

        ### Also a partial set of signatures could be returned
        ### In that way it would be possible to figure out who did not sign the ballot
        ### and whether that happens to be a consistent behaviour which can be qualified as DDOS attack
        ballot, signatures = run(config, usersockets, mix)

        ### shall it actually error?
        audit!(ar, ballot, signatures)
    end
end

function vote(port, membergate::Layer, memberballot::Layer, msg::Vector{UInt8}, sign::Function)
    securesocket = secure(connect(port),membergate)
    sroutersocket = secure(securesocket,memberballot)

    @assert read(sroutersocket,UInt8)==OPEN
    stack(sroutersocket,msg)

    metadata = unstack(securesocket)
    rawballot = unstack(securesocket)

    M = length(msg)

    @assert mod(length(rawballot),M) == 0
    N = div(length(rawballot),M)

    ballot = reshape(rawballot,(M,N))

    s = sign(Ballot(metadata,ballot))
    stack(securesocket,s)
end

function vote(msg::Vector{UInt8}, gk::GateKeeper, s::Union{Signer, Secret}, sign::Function)
    membergate = Layer(gk.crypto, s, gk.id)
    memberballot = Layer(gk.mix.crypto, gk.mix.id)
    return vote(gk.port, membergate, memberballot, msg, sign)
end

export Mix, GateKeeper, Officer, Ballot, Regulation, vote

end


text:
module regression

using LinearAlgebra
using Statistics
using ..utils

include("metrics.jl"); export max_error, mean_absolute_error, mean_squared_error, mean_squared_log_error, median_absolute_error, mean_absolute_percentage_error

end


text:
function pressure(model::EoSModel, v, T,  z=SA[1.])
    return -∂f∂v(model,v,T,z)
end

function vt_entropy(model::EoSModel, v, T,   z=SA[1.])
    return -∂f∂t(model,v,T,z)
end

function vt_chemical_potential(model::EoSModel, v, T, z= SA[1.])
    fun(x) = eos(model, v, T,x)
    return ForwardDiff.gradient(fun,z)
end

function vt_internal_energy(model::EoSModel, v, T,  z=SA[1.])
    _df,f =  ∂f(model,v,T,z)
    dv,dt = _df
    return f  - dt*T
end

function vt_enthalpy(model::EoSModel, v, T,  z=SA[1.])
    _df,f =  ∂f(model,v,T,z)
    dv,dt = _df
    return f  - dv*v - dt*T
end

function vt_gibbs_free_energy(model::EoSModel, v, T,  z=SA[1.])
    _df,f =  ∂f(model,v,T,z)
    dv,dt = _df
    return f  - dv*v
end

function vt_helmholtz_free_energy(model::EoSModel, v, T,  z=SA[1.])
    return eos(model, v, T,z)
end

function vt_isochoric_heat_capacity(model::EoSModel, v, T,  z=SA[1.])
    fun(x)  = eos(model, v, x,z)
    df(x)   = ForwardDiff.derivative(fun,x)
    d2f(x)  = ForwardDiff.derivative(df,x)
    return -T*d2f(T)
end

function vt_isobaric_heat_capacity(model::EoSModel, v, T,  z=SA[1.])
    d2f = f_hess(model,v,T,z)
    return T*(d2f[1,2]^2/d2f[1]-d2f[2,2])
end

function vt_isothermal_compressibility(model::EoSModel, v, T,  z=SA[1.])
    p0,dpdv = p∂p∂v(model,v,T,z)
    return -1/v*dpdv^-1
end

function vt_isentropic_compressibility(model::EoSModel, v, T,  z=SA[1.])
    d2f = f_hess(model,v,T,z)
    return 1/v*(d2f[1]-d2f[1,2]^2/d2f[2,2])^-1
end

function vt_speed_of_sound(model::EoSModel, v, T,  z=SA[1.])
    Mr      = molecular_weight(model,z)
    d2f = f_hess(model,v,T,z)
    return v*sqrt((d2f[1]-d2f[1,2]^2/d2f[2,2])/Mr)
end

function vt_isobaric_expansivity(model::EoSModel, v, T,  z=SA[1.])
    d2f = f_hess(model,v,T,z)
    return d2f[1,2]/(v*d2f[1])
end

function vt_joule_thomson_coefficient(model::EoSModel, v, T,  z=SA[1.])

    d2f = f_hess(model,v,T,z)
    return -(d2f[1,2]-d2f[1]*((T*d2f[2,2]+v*d2f[1,2])/(T*d2f[1,2]+v*d2f[1])))^-1
end

function second_virial_coefficient(model::EoSModel, T,  z=SA[1.])
    TT = promote_type(eltype(z),typeof(T))
    V = 1/sqrt(eps(TT))
    _∂2f = ∂2f(model,V,T,z)
    hessf,gradf,f = _∂2f
    _p,dpdv = p∂p∂v(model,V,T,z)
    return -V^2/(R̄*T)*(_p+V*dpdv)
end

function vt_compressibility_factor(model::EoSModel, v, T,  z=SA[1.])
    p = pressure(model,v,T,z)
    return p*v/(R̄*T)
end

"""
    pip(model::EoSModel,v,T,z=[1.0])

Phase identification parameter `Π`. as described in _1_. If `Π > 1`, then the phase is clasified as a liquid or a liquid-like vapor, being a vapor or vapor-like liquid otherwise.

This identification parameter fails at temperatures and pressures well above the critical point.

Calculated as:
```
Π = v*((∂²p/∂v∂T)/(∂p/∂T) - (∂²p/∂v²)/(∂p/∂v))
```

1.  G. Venkatarathnama, L.R. Oellrich, Identification of the phase of a fluid using partial derivatives of pressure, volume,and temperature without reference to saturation properties: Applications in phase equilibria calculations, Fluid Phase Equilibria 301 (2011) 225–233

"""
function pip(model::EoSModel,v,T,z=SA[1.0])
    _∂2p = ∂2p(model,v,T,z)
    hess_p, grad_p, _ = _∂2p
    Π = v*(hess_p[1,2]/grad_p[2]  - hess_p[1,1]/grad_p[1])
end

export second_virial_coefficient


text:
"""
GenX: An Configurable Capacity Expansion Model
Copyright (C) 2021,  Massachusetts Institute of Technology
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
A complete copy of the GNU General Public License v2 (GPLv2) is available
in LICENSE.txt.  Users uncompressing this from an archive may not have
received this license file.  If not, see <http://www.gnu.org/licenses/>.
"""

@doc raw"""
	write_nse(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)

Function for reporting non-served energy for every model zone, time step and cost-segment.
"""
function write_nse(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)
	dfGen = inputs["dfGen"]
	T = inputs["T"]     # Number of time steps (hours)
	Z = inputs["Z"]     # Number of zones
    	SEG = inputs["SEG"] # Number of load curtailment segments
	# Non-served energy/demand curtailment by segment in each time step
	dfNse = DataFrame()
	dfTemp = Dict()
	for z in 1:Z
		dfTemp = DataFrame(Segment=zeros(SEG), Zone=zeros(SEG), AnnualSum = Array{Union{Missing,Float32}}(undef, SEG))
		dfTemp[!,:Segment] = (1:SEG)
		dfTemp[!,:Zone] = fill(z,(SEG))
		if setup["ParameterScale"]==1
			for i in 1:SEG
				dfTemp[!,:AnnualSum][i] = sum(inputs["omega"].* (value.(EP[:vNSE])[i,:,z]))* ModelScalingFactor
			end
			dfTemp = hcat(dfTemp, DataFrame((value.(EP[:vNSE])[:,:,z])* ModelScalingFactor, :auto))
		else
			for i in 1:SEG
				dfTemp[!,:AnnualSum][i] = sum(inputs["omega"].* (value.(EP[:vNSE])[i,:,z]))
			end
			dfTemp = hcat(dfTemp, DataFrame(value.(EP[:vNSE])[:,:,z], :auto))
		end
		if z == 1
			dfNse = dfTemp
		else
			dfNse = vcat(dfNse,dfTemp)
		end
	end

	auxNew_Names=[Symbol("Segment");Symbol("Zone");Symbol("AnnualSum");[Symbol("t$t") for t in 1:T]]
	rename!(dfNse,auxNew_Names)
	total = DataFrame(["Total" 0 sum(dfNse[!,:AnnualSum]) fill(0.0, (1,T))], :auto)
	for t in 1:T
		total[:,t+3] .= sum(dfNse[:,Symbol("t$t")][1:Z])
	end
	rename!(total,auxNew_Names)
	dfNse = vcat(dfNse, total)

	CSV.write(joinpath(path, "nse.csv"),  dftranspose(dfNse, false), writeheader=false)
	return dfTemp
end

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM