text:
requiredTypes = [ "AcademicLevel", "Affiliation", "Assignment", "CareerStatus", "Job", "Rank", "SubJob", "CareerPath" ]

for reqType in requiredTypes
    if !isdefined( Symbol( uppercase( string( reqType[ 1 ] ) ) * reqType[ 2:end ] ) )
        include( joinpath( dirname( Base.source_path() ), "..", "functions", reqType * ".jl" ) )
    end  # if !isdefined( Symbol( ...
end
#println(joinpath( dirname( Base.source_path() ), "..", "..", "..", GuyCareerPathsAll.xlsx ) )

#_______________________________________________________________________________
# Input variables

MaxCareerPathLength = 40

SimulationYear = 2018

XlsxFile = joinpath( dirname( Base.source_path() ), "..", "..", "..", "GuyCareerPathsAllPop.xlsx" )

requiredData = ["CPByGuyAllpop"] # GuyCareerPaths , InitManpower , MPObjectives

for reqData in requiredData
    if !isdefined( Symbol( uppercase( string( reqData[ 1 ] ) ) * reqData[ 2:end ] ) )
        include( joinpath( dirname( Base.source_path() ), "..", "Data", reqData * ".jl" ) )
    end  # if !isdefined( Symbol( ...
end

#_______________________________________________________________________________
# Variables number

    # Compute which careers could each Initial manpower subset go with

InitMPPartCP = Array{Vector{Int}}(length(InitManpower))
InitMPDivisionNb = 0
for i in 1:length(InitManpower)
    InitIpartCP = Vector{Int}()
    for j in 1:length(GuyCareerPaths)
        index = get_status_index_after_duration(GuyCareerPaths[j], InitManpower[i].Seniority + 1)
        if index == -1
            continue
        end
        if contains(get_Name_Level(GuyCareerPaths[j].Path[index], AcademicLevel), InitManpower[i].Academiclvl)
            push!(InitIpartCP, j)
        end
    end
    InitMPPartCP[i] = InitIpartCP
    InitMPDivisionNb += length(InitIpartCP)
end
InitMPDivNbOnePop = InitMPDivisionNb
InitMPDivisionNb *= length(SubPopulations)

AnnualRecDivNbOnePop = length(GuyCareerPaths) * NBYears
AnnualRecDivNb = AnnualRecDivNbOnePop * length(SubPopulations)

DeviationVarNb = length(MPObjectives) * NBYears * 2

DeviationSubPop = length(SubPopulations)* NBYears * 2

VarNb = InitMPDivisionNb + AnnualRecDivNb + DeviationVarNb + DeviationSubPop

#_______________________________________________________________________________
# Constraints number

InitConst = length(InitManpower)

RecruitmentConst = NBYears

GoalsConst = NBYears * length(MPObjectives)

RecruitmentDeviation = (NBYears - 1) * 2

SubPopConst = length(SubPopulations)* NBYears

ConstNb = InitConst + RecruitmentConst + GoalsConst + RecruitmentDeviation + SubPopConst

#_______________________________________________________________________________
# MIP construction
    # Definition of Arrays

Cost = Array{Float64}(VarNb) # Done

A = zeros(Float64, ConstNb, VarNb) # Init Done, Rec Done, Goal Done

sense = Array{Char}(ConstNb) # Done

b = zeros(Float64, ConstNb)#Array{Float64}(ConstNb) # Done

vartypes = Array{Symbol}(VarNb) # Done

    # Fill Cost
Cost[1:(InitMPDivisionNb)] = 0
Cost[1 + InitMPDivisionNb:(InitMPDivisionNb + AnnualRecDivNb)] = 0
for i in 1:length(MPObjectives)
    for j in 1:NBYears
        Cost[(InitMPDivisionNb + AnnualRecDivNb) +                                  (j - 1)*length(MPObjectives) + i ] = MPObjectives[i].priority
        Cost[(InitMPDivisionNb + AnnualRecDivNb) + (length(MPObjectives)*NBYears) + (j - 1)*length(MPObjectives) + i ] = MPObjectives[i].priority
    end
end
for i in 1:length(SubPopulations)
    for j in 1:NBYears
        Cost[(InitMPDivisionNb + AnnualRecDivNb) + DeviationVarNb +                                    (j - 1)*length(SubPopulations) + i ] = SubPopulations[i].priority
        Cost[(InitMPDivisionNb + AnnualRecDivNb) + DeviationVarNb + (length(SubPopulations)*NBYears) + (j - 1)*length(SubPopulations) + i ] = SubPopulations[i].priority
    end
end
    # Fill A
        # Initial Manpower Constraints
AdvacmentY = 0
for i in 1:InitConst
    for j in 1:length(InitMPPartCP[i])
        for k in 1:length(SubPopulations)
            A[i, InitMPDivNbOnePop*(k-1) + AdvacmentY + j] = 1
        end
    end
    AdvacmentY += length(InitMPPartCP[i])
end

        # Annual Recruitment

for i in 1:(RecruitmentConst)
    for j in 1:length(GuyCareerPaths)
        for k in 1:length(SubPopulations)
            A[InitConst + i, InitMPDivisionNb + AnnualRecDivNbOnePop*(k-1) + (i-1)*length(GuyCareerPaths) + j] = 1
        end
    end
end

        # Goals

            # compute which careers go with each goal

SetsPartCPinGoals = Vector{Vector{Vector{Int}}}()

for i in 1:length(GuyCareerPaths)
    AllParticipations = Vector{Vector{Int}}()
    for j in 1:length(GuyCareerPaths[i].Path)
        Participation = Vector{Int}()

        for k in 1:length(MPObjectives)
            if contains(get_Name_Level(GuyCareerPaths[i].Path[j], MPObjectives[k].targetLevel[1]), MPObjectives[k].Objectives[1])#_______________________________________________________________________________
                push!(Participation, k)
            end
        end
        push!(AllParticipations, Participation)
    end
    push!(SetsPartCPinGoals, AllParticipations)
end

            # Compute Career Paths Length

for i in 1:length(GuyCareerPaths)
    compute_Length(GuyCareerPaths[i])
end

            # Initial Manpower participation

TempActualVar = 1
for imp in 1:length(InitManpower)
    for i in 1:length(InitMPPartCP[imp])
        for j in (1 + InitManpower[imp].Seniority):GuyCareerPaths[InitMPPartCP[imp][i]].Length
            indx = get_status_index_after_duration(GuyCareerPaths[InitMPPartCP[imp][i]], j)
            if indx == -1
                break
            end
            for k in 1:length(SetsPartCPinGoals[InitMPPartCP[imp][i]][indx])
                if (j - 1) <= NBYears
                    for l in 1:length(SubPopulations)
                        A[InitConst + RecruitmentConst + (j - InitManpower[imp].Seniority - 1)*length(MPObjectives) + SetsPartCPinGoals[InitMPPartCP[imp][i]][indx][k], InitMPDivNbOnePop*(l-1) + TempActualVar] = 1
                    end
                end
            end
        end
        TempActualVar += 1
    end
end

            # Annual Recruitment participation

for y in 1:NBYears
    for i in 1:length(GuyCareerPaths)
        for j in 1:GuyCareerPaths[i].Length # j => year for const
            indx = get_status_index_after_duration(GuyCareerPaths[i], j)
            if indx == -1
                break
            end
            for k in 1:length(SetsPartCPinGoals[i][indx])
                if (y + j - 1) <= NBYears
                    for l in 1:length(SubPopulations)
                        A[InitConst + RecruitmentConst + (y - 1)*length(MPObjectives) + (j - 1)*length(MPObjectives) + SetsPartCPinGoals[i][indx][k], InitMPDivisionNb + AnnualRecDivNbOnePop*(l-1) + (y - 1)*length(GuyCareerPaths) + i] = 1 - AttritionAtYear(GuyCareerPaths[i], j)
                    end
                end
            end
        end
    end
end

            # deviation variables
for y in 1:NBYears
    for eq in 1:length(MPObjectives)

        A[InitConst + RecruitmentConst + (y - 1)*length(MPObjectives) + eq, InitMPDivisionNb + AnnualRecDivNb +                                  (y - 1)*length(MPObjectives) + eq] =  1
        A[InitConst + RecruitmentConst + (y - 1)*length(MPObjectives) + eq, InitMPDivisionNb + AnnualRecDivNb + length(MPObjectives) * NBYears + (y - 1)*length(MPObjectives) + eq] = -1

    end
end

            # Recruitment deviation

for i in 1:NBYears-1
    for j in 1:length(GuyCareerPaths)
        for k in 1:length(SubPopulations)
            A[InitConst + RecruitmentConst + GoalsConst + i, InitMPDivisionNb + AnnualRecDivNbOnePop*(k-1) + (i  )*length(GuyCareerPaths) + j] = 1
            A[InitConst + RecruitmentConst + GoalsConst + i, InitMPDivisionNb + AnnualRecDivNbOnePop*(k-1) + (i-1)*length(GuyCareerPaths) + j] = -1.10
        end
    end
end

for i in 1:NBYears-1
    for j in 1:length(GuyCareerPaths)
        for k in 1:length(SubPopulations)
            A[InitConst + RecruitmentConst + GoalsConst + (NBYears - 1) + i, InitMPDivisionNb + AnnualRecDivNbOnePop*(k-1) + (i  )*length(GuyCareerPaths) + j] = 1
            A[InitConst + RecruitmentConst + GoalsConst + (NBYears - 1) + i, InitMPDivisionNb + AnnualRecDivNbOnePop*(k-1) + (i-1)*length(GuyCareerPaths) + j] = -0.90
        end
    end
end

# Subpopulations fulfillement

    # Participation of the recruitment
for i in 1:NBYears
    for j in 1:length(SubPopulations)
        for k in 1:length(GuyCareerPaths)
            A[InitConst + RecruitmentConst + GoalsConst + RecruitmentDeviation + (i-1)*length(SubPopulations) + j, InitMPDivisionNb + ] = 1 #- AttritionAtYear(GuyCareerPaths[k], )
        end
    end
end

    # Fill sense

sense[1:InitConst] = '='
sense[(InitConst + 1):(InitConst + RecruitmentConst)] = '<'
sense[(InitConst + RecruitmentConst + 1):(InitConst + RecruitmentConst + GoalsConst)] = '='
sense[(InitConst + RecruitmentConst + GoalsConst + 1):(InitConst + RecruitmentConst + GoalsConst + (NBYears - 1))] = '<'
sense[(InitConst + RecruitmentConst + GoalsConst + (NBYears - 1) + 1):end] = '>'

    # Fill b
for i in 1:InitConst
    b[i] = InitManpower[i].Nb
end
b[(InitConst + 1):(InitConst + RecruitmentConst)] = MaxRecruitment
index = InitConst + RecruitmentConst + 1
for y in 1:NBYears
    for i in 1:length(MPObjectives)
        b[index] = MPObjectives[i].Number
        index += 1
    end
end

    # Fill vartypes

vartypes[1:(InitMPDivisionNb + AnnualRecDivNb)] = :Int
vartypes[(InitMPDivisionNb + AnnualRecDivNb + 1):end] = :Cont

writedlm("matrix.txt", A)

#_______________________________________________________________________________
# MIP resolution
println("Solving LP...")
using MathProgBase, CPLEX#, Gurobi
stattoendStart = now()
sol = mixintprog(Cost, A, sense, b, vartypes, 0, Inf, CplexSolver(CPXPARAM_MIP_Tolerances_MIPGap=0.01))#CbcSolver(allowableGap=0.8)) #GurobiSolver(Presolve=0)
stattoendEnd = now()
println( "ended with: $(sol.status). Elapsed time: $(Dates.canonicalize(Dates.CompoundPeriod(Dates.Millisecond(stattoendEnd - stattoendStart))))." )
writedlm("solution.txt", sol.sol)

#_______________________________________________________________________________
# Computing results

RequirementFulfillment = zeros(Float64, NBYears, length(MPObjectives))

for i in 1:NBYears
    for j in 1:length(MPObjectives)
        for k in 1:(InitMPDivisionNb + AnnualRecDivNb)
            if A[InitConst + RecruitmentConst + (i - 1)*length(MPObjectives) + j, k] != 0
                RequirementFulfillment[i, j] += round(Int, sol.sol[k]*A[InitConst + RecruitmentConst + (i - 1)*length(MPObjectives) + j, k])
            end
        end
    end
end

YearlyRecruitment = zeros(Float64, NBYears)

for i in 1:NBYears
    for j in 1:length(GuyCareerPaths)
        for k in 1:length(SubPopulations)
            YearlyRecruitment[i] += sol.sol[InitMPDivisionNb + AnnualRecDivNbOnePop*(k-1) + (i - 1)*length(GuyCareerPaths) + j]
        end
    end
end

include("WrittingResultsAllPop.jl")

max_stars_repo_name: Omazaria/CareerPathBeDef

max_issues_repo_name: Omazaria/CareerPathBeDef

max_forks_repo_name: Omazaria/CareerPathBeDef


text:
function mle_Sigma_L!(Sigma_L::AbstractMatrix, S::AbstractMatrix, Sigma_X::SparseMatrixCSC, A::SparseMatrixCSC)

    # invert Sigma_X efficiently
    Theta_X = deepcopy(Sigma_X)
    Theta_X.nzval[:] = 1 ./ Theta_X.nzval

    B = A'*Theta_X*S*Theta_X*A
    C = A'Theta_X*A
    Cinv = deepcopy(C)
    Cinv.nzval[:] = 1 ./ Cinv.nzval # assumes that C is diagonal!!
    Sigma_L[:,:] = Cinv*(B - C)*Cinv
end

# not correct! you can't just take the diagonal to get the optimal Sigma_X
# function mle_Sigma_X!(Sigma_X::SparseMatrixCSC, S::AbstractMatrix, A::SparseMatrixCSC, Sigma_L::AbstractMatrix)
#     Sigma_X.nzval[:] = diag(S - A*Sigma_L*A')
# end
# function mle_Sigma_XandL!(Sigma_X::SparseMatrixCSC, Sigma_L::AbstractMatrix, S::AbstractMatrix, A::SparseMatrixCSC; iterations=1000)
#     lastll = -Inf
#     for i in 1:iterations
#         #println(ll_score(S, Sigma_Xest, Sigma_Lest, Aest))
#         Sigma_Lest[:,:] = mle_Sigma_L(S, Sigma_Xest, Aest)
#         #println("md ", minimum(diag(Sigma_Lest)))
#         Sigma_Xest[:,:] = mle_Sigma_X(S, Sigma_Lest, Aest)
#         #println("md ", minimum(diag(Sigma_Xest)))
#         ll = loglikelihood(S, Sigma_X, A, Sigma_L)
#         @assert ll > lastll
#         if ll - lastll < 1e-8 break end
#         println(ll)
#     end
# end

ext: jl


text:
mutable struct Scatter{N}
    positions::Vector{Point{N, Float32}}
    color::TorVector{RGBAf0}
    marker::TorVector{<: Union{Symbol, Char, Type{Circle}}}
    markersize::TorVector{Union{Float32, Vec{N, Float32}}}
    markeroffset::TorVector{Vec{N, Float32}}
    strokecolor::TorVector{RGBAf0}
    strokewidth::TorVector{Float32}
    markerspace::Space
    transform_marker::Bool
    on_update::Observable{Dict{Symbol, Any}}
    on_event::Observable{Tuple{Symbol, Any}}
    camera::Camera
    transformation::Transformation
end

optional_fields(T::Type{Scatter}) = keys(defaults(T))

function defaults(::Type{Scatter})
    return (
        color = :black,
        marker = :circle,
        markersize = 10f0,
        markeroffset = Vec2f0(0),
        strokecolor = :white,
        strokewidth = 1f0,
        markerspace = Pixel,
        transform_marker = false
    )
end

function Scatter(positions; kw...)
    sdefaults = defaults(Scatter)
    pos = convert_arguments(Scatter, positions)
    fields = map(optional_fields(Scatter)) do name
        value = get(kw, name, getfield(sdefaults, name))
        return convert_attribute(Scatter, Key(name), value)
    end
    on_update = Observable{Dict{Symbol, Any}}()
    on_event = Observable{Tuple{Symbol, Any}}()
    camera = Camera()
    transformation = Transformation()
    return Scatter(pos, fields..., on_update, on_event, camera, transformation)
end

function update!(scatter::Scatter, values::Dict{Symbol})
    scatter.on_update[] = values
end

function Base.setproperty!(scatter::Scatter, field::Symbol, value::Any)
    converted = if field === :positions
        convert_arguments(Scatter, value)
    else
        convert_attribute(Scatter, Key(field), value)
    end
    setfield!(scatter, field, converted)
    update!(scatter, Dict(field => converted))
end


text:
using RealAlg, UnitSystems

CellData = Cell(Constants(),Geometry(),Negative(),Positive(),Seperator(),RealisationAlgorthim())
loc = (Number[0, 128e-6, 204e-6, 394e-6],Number[0,1], Number[0,1], Number[128e-6, 204e-6, 394e-6],Number[1],Number[1],Number[0,1],Number[0,1],Number[0,1],Number[0,1])

Arr_Factor = (1/CellData.Const.T_ref-1/CellData.Const.T)/universal(SI2019)
FCall= FCalls(kappa(CellData.Const.ce0))
κ_dra=FCall.Kap.κ*exp(CellData.Const.Ea_κ*Arr_Factor)
TransferFuns = TransferFun()

const A_DRA = Array{Float64}(undef,0,CellData.RA.M)
const B_DRA = Array{Float64}(undef,0,CellData.RA.M)
const C_DRA = Array{Float64}(undef,0,CellData.RA.M)
const D_DRA = Array{Float64}(undef,0,CellData.RA.M)
for Temperature in 5.0:5.0:5.0
    CellData.Const.T = 273.15+Temperature
    for SOC in 0:0.05:0.05
        CellData.Const.SOC = SOC
        A, B, C, D = DRA(CellData,FCall,loc,TransferFuns)
        println("DRA")
    end
end

max_stars_repo_name: BradyPlanden/RealisationAlgorithms

max_issues_repo_name: BradyPlanden/RealisationAlgorithms

max_forks_repo_name: BradyPlanden/RealisationAlgorithms


text:
# HttpServer module
#
# Serve HTTP requests in Julia.
#
__precompile__()
module HttpServer

using HttpCommon
using MbedTLS

include("RequestParser.jl")

export HttpHandler,
       Server,
       run,
       write,
       # from HttpCommon
       Request,
       Response,
       FileResponse,
       escapeHTML,
       parsequerystring,
       setcookie!

import Base: run, listen

defaultevents = Dict{ASCIIString, Function}()
defaultevents["error"]  = ( client, err )->println( err )
defaultevents["listen"] = ( saddr )      ->println("Listening on $saddr...")

"""
`HttpHandler` types are used to instantiate a `Server`.

An `HttpHandler` is responsible for the entirety of the request / response
cycle. Instantiate an `HttpHandler` with a single `Function` argument, which
becomes `HttpHandler.handle`. This handler is called on every incoming
request and passed `req::Request, res::Response`. The return value of
`HttpHandler.handle` is the response sent to the client for the given `req`.

`HttpHandler.handle` can return a `String`:

```
handler = HttpHandler() do req, res
  "Hello world!"
end
```

an HTTP status code as `Int`:

```
handler = HttpHandler() do req, res
  404
end
```

or a full `Response` instance:

```
handler = HttpHandler() do req, res
  Response(200, "Success", "Hello World!"
end
```

`HttpHandler.sock` is used internally to store the `TcpSocket` for incoming
connections.

`HttpHandler.events` is a dictionary of functions to call when certain
server events occur.
Set these functions by direct assignment:

```
handler.events["listen"] = (port) "Listening on \$port..."
```

All default events and their arguments:

* `"listen"  => (port::Int)`
* `"connect" => (client::Client)`
* `"write"   => (client::Client, res::Response)`
* `"close"   => (client::Client)`

If you want to trigger custom events on your server, use the `event`
function:

```
# listen for "foo"
handler.events["foo"] = (bar) "Hello \$bar"
# trigger "foo"
HttpServer.event(server, "foo", "Julia")
```
"""
immutable HttpHandler
    handle::Function
    sock::Base.UVServer
    events::Dict

    HttpHandler(handle::Function, sock::Base.UVServer) = new(handle, sock, defaultevents)
    HttpHandler(handle::Function) = new(handle, Base.TCPServer(), defaultevents)
end
handle(handler::HttpHandler, req::Request, res::Response) = handler.handle(req, res)

""" Client encapsulates a single connection

 When new connections are initialized a `Client` is created with a new id and
 the connection socket. `Client.parser` will store a `ClientParser` to handle
 all HTTP parsing for the connection lifecycle.
"""
type Client{T<:IO}
    id::Int
    sock::T
    parser::ClientParser

    Client(id::Int,sock::T) = new(id,sock)
end
Client{T<:IO}(id::Int,sock::T) = Client{T}(id,sock)

""" `WebSocketInterface` defines the abstract protocol for a WebSocketHandler.

The methods `is_websocket_handshake` and `handle` will be called if `Server.
sockets` is populated. Concrete types of `WebSocketInterface` are required
to define these methods.
"""
abstract WebSocketInterface

""" `is_websocket_handshake` should determine if `req` is a valid websocket
upgrade request.
"""
function is_websocket_handshake(handler::WebSocketInterface, req::Request)
    throw("`$(typeof(handler))` does not implement `is_websocket_handshake`.")
end

""" `handle` is called when `is_websocket_handshake` returns true, and takes
full control of the connection.
"""
function handle(handler::WebSocketInterface, req::Request, client::Client)
    throw("`$(typeof(handler))` does not implement `handle`.")
end

""" `Server` types encapsulate an `HttpHandler` and optional
`WebSocketInterface` to serve requests.

* Instantiate with both an `HttpHandler` and `WebSocketInterface` to serve
  both protocols.
* Instantiate with just an `HttpHandler` to serve only standard `Http`
* Instantiate with just a `Function` to create an `HttpHandler` automatically
* Instantiate with just a `WebSocketInterface` to only serve websockets
  requests and `404` all others.
"""
immutable Server
    http::HttpHandler
    websock::Union(Nothing, WebSocketInterface)
end
Server(http::HttpHandler)           = Server(http, nothing)
Server(handler::Function)           = Server(HttpHandler(handler))
Server(websock::WebSocketInterface) = Server(HttpHandler((req, res)->Response(404)), websock)

"""Triggers `event` on `server`.

If there is a function bound to `event` in `server.events` it will be called
with `args...`
"""
function event(event::String, server::Server, args...)
    haskey(server.http.events, event) && server.http.events[event](args...)
end

"""
Sets a cookie with the given name, value, and attributes on the given response object.
"""
setcookie!(r::Response, name, value=utf8(""), attrs=Dict{UTF8String, UTF8String}()) =
  setcookie!(r, Cookie(name, value, attrs))

setcookie!(r::Response, cookie::Cookie) = (r.cookies[cookie.name] = cookie; r)

import Base.write
"Converts a `Response` to an HTTP response string"
function write(io::IO, response::Response)
    write(io, join(["HTTP/1.1", response.status, HttpCommon.STATUS_CODES[response.status], "\r\n"], " "))

    response.headers["Content-Length"] = string(sizeof(response.data))
    for (header,value) in response.headers
        write(io, string(join([ header, ": ", value ]), "\r\n"))
    end
    for (cookie_name, cookie) in response.cookies
        write(io, "Set-Cookie: ", cookie_name, "=", cookie.value)
        for (attr_name, attr_val) in cookie.attrs
            write(io, "; ", attr_name)
            if !isempty(attr_val)
                write(io, "=", attr_val)
            end
        end
        write(io, "\r\n")
    end

    write(io, "\r\n")
    write(io, response.data)
end

""" Start `server` to listen on specified socket address.

    listen(server::Server, host::Base.IpAddr, port::Integer) Server

    Setup "server" so it listens on "port" on the address specified by "host".
    To listen on all interfaces pass, "IPv4(0)" or "IPv6(0)" as appropriate.
"""
function listen(server::Server, host::Base.IpAddr, port::Integer)
    Base.uv_error("listen", !bind(server.http.sock, host, UInt16(port)))
    listen(server.http.sock)
    inet = "$host:$port"
    event("listen", server, inet)
    return server
end
listen(server::Server, port::Integer) = listen(server, IPv4(0), port)

""" Start `server` to listen on named pipe/domain socket.

    listen(server::Server, path::String) Server

    Setup "server" to listen on named pipe/domain socket specified by "path".
"""
@unix_only function listen(server::Server, path::String)
    bind(server.http.sock, path) || throw(ArgumentError("could not listen on path $path"))
    Base.uv_error("listen", Base._listen(server.http.sock))
    event("listen", server, path)
    return server
end

""" Handle HTTP request from client """
function handle_http_request(server::Server)
    id_pool = 0 # Increments for each connection
    websockets_enabled = server.websock != nothing
    while true # handle requests, Base.wait_accept blocks until a connection is made
        client = Client(id_pool += 1, accept(server.http.sock))
        client.parser = ClientParser(message_handler(server, client, websockets_enabled))
        @async process_client(server, client, websockets_enabled)
    end
end

""" Handle HTTPS request from client """
function handle_https_request(server::Server, ssl_config::MbedTLS.SSLConfig)
    id_pool = 0 # Increments for each connection
    websockets_enabled = server.websock != nothing
    while true
        sess = MbedTLS.SSLContext()
        MbedTLS.setup!(sess, ssl_config)
        # set_priority_string!(sess)
        # set_credentials!(sess, cert_store)
        client = accept(server.http.sock)
        try
            MbedTLS.set_bio!(sess, client)
            MbedTLS.handshake(sess)
            # associate_stream(sess, client)
            # handshake!(sess)
        catch e
            println("Error establishing SSL connection: ", e)
            close(client)
            continue
        end
        client = Client(id_pool += 1, sess)
        client.parser = ClientParser(message_handler(server, client, websockets_enabled))
        @async process_client(server, client, websockets_enabled)
    end
end

function default_ssl_config(ssl_cert, key)
    conf = MbedTLS.SSLConfig()
    entropy = MbedTLS.Entropy()
    rng = MbedTLS.CtrDrbg()
    MbedTLS.config_defaults!(conf, endpoint=MbedTLS.MBEDTLS_SSL_IS_SERVER)
    MbedTLS.seed!(rng, entropy)
    MbedTLS.rng!(conf, rng)
    MbedTLS.own_cert!(conf, ssl_cert, key)
    MbedTLS.dbg!(conf, (level, filename, number, msg)->begin
        warn("MbedTLS emitted debug info: $msg in $filename:$number")
    end)
    conf
end

function handle_https_request(server::Server, ssl_cert::Tuple{MbedTLS.CRT, MbedTLS.PKContext})
    handle_https_request(server, default_ssl_config(ssl_cert...))
end

""" `run` starts `server`

Functionality:

* Accepts incoming connections and instantiates each `Client`.
* Manages the `client.id` pool.
* Spawns a new `Task` for each connection.
* Blocks forever.

Method accepts following keyword arguments:

* host - binding address
* port - binding port
* ssl  - SSL configuration. Use this argument to enable HTTPS support.
Can be either an MbedTLS.SSLConfig object that already
has associated certificates, or a tuple of an MbedTLS.CRT (certificate) and
MbedTLS.PKContext (private key)). In the latter case, a configuration with reasonable
defaults will be used.

* socket - named pipe/domain socket path. Use this argument to enable Unix socket support.
It's available only on Unix. Network options are ignored.

Compatibility:

* for backward compatibility use `run(server::Server, port::Integer)`

Example:
```
server = Server() do req, res
  "Hello world"
end

# start server on localhost
run(server, host=IPv4(127,0,0,1), port=8000)
# or
run(server, 8000)
```
"""
function run(server::Server; args...)
    params = Dict(args)

    # parse parameters
    port = get(params, :port, 0)
    host = get(params, :host, IPv4(0))
    use_https = haskey(params, :ssl)
    use_sockets = @unix? (haskey(params, :socket) && !haskey(params, :port)) : false

    server = if use_sockets
        listen(server, params[:socket])  # start server on Unix socket
    else
        host, port
        listen(server, host, port)       # start server on network socket
    end

    if use_https
        handle_https_request(server, params[:ssl])
    else
        handle_http_request(server)
    end
end
# backward compatibility method
run(server::Server, port::Integer) = run(server, port=port)

"""Handles live connections, runs inside a `Task`.

Blocks ( yields ) until a line can be read, then passes it into `client.
parser`.
"""
function process_client(server::Server, client::Client, websockets_enabled::Bool)
    event("connect", server, client)

    while !eof(client.sock)#isopen(client.sock)
        # try
            if !upgrade(client.parser.parser)
                # IMPORTANT NOTE: This is technically incorrect as there may be data
                # in the buffer that we have not yet read. The way to deal with this
                # is to manually adjust the position of the buffer, but for that to
                # happen, we need use the return value of http_parser_execute, which
                # we don't have, since we launch websocket handlers in the callback.
                # Anyway, since there always needs to be a handshake this is probably
                # fine for now.
                data = readavailable(client.sock)
                if length(data) > 0
                    add_data(client.parser, data)
                end
            else
                wait(client.sock.closenotify)
            end
        # catch e
        #     if isa(e,GnuTLS.GnuTLSException) && e.code == -110
        #         # Because Chrome is buggy on OS X, ignore E_PREMATURE_TERMINATION
        #     else
        #         rethrow()
        #     end
        # end
    end
    event("close", server, client)
end

"Callback factory for providing `on_message_complete` for `client.parser`"
function message_handler(server::Server, client::Client, websockets_enabled::Bool)

    # Called when the `ClientParser` has finished parsing a `Request`.
    #
    # If websockets are enabled and claim this request, defer to `server.
    # websock` and return. Otherwise, call `server.http.handle` to get a
    # `Response` for the `client`. Catches any errors with `server.http.
    # handle`, returns a `500` and writes a stacktrace to `stdout`.
    #
    function on_message_complete(req::Request)
        if websockets_enabled && is_websocket_handshake(server.websock, req)
            handle(server.websock, req, client)
            return true
        end

        local response

        try
            response = handle(server.http, req, Response()) # Run the server handler
            if !isa(response, Response)                     # Promote return to Response
                response = Response(response)
            end
        catch err
            response = Response(500)
            event("error", server, client, err)
            Base.display_error(err, catch_backtrace())      # Prints backtrace without throwing
        end

        # We have a Keep-Alive header set to 0 or 1 in RequestParser
        # to denote if the client wants keep-alive or not
        # We also have to check if the server wants to close the connection
        # If this is keep-alive, we need to add Connection header with keep-alive
        local keep_alive = get(req.headers, "Keep-Alive", nothing) == "1" &&
                           get(response.headers, "Connection", nothing) != "close"
        if keep_alive
            response.headers["Connection"] = "keep-alive"
        end

        write(client.sock, response)
        event("write", server, client, response)

        if !keep_alive
            close(client.sock)
        end
    end
end

include("mimetypes.jl")

function FileResponse(filename)
    if isfile(filename)
        s = open(readbytes,filename)
        (_, ext) = splitext(filename)
        mime = length(ext)>1 && haskey(mimetypes,ext[2:end]) ? mimetypes[ext[2:end]] : "application/octet-stream"
        Response(200, Dict{String,String}([("Content-Type",mime)]), s)
    else
        Response(404, "Not Found - file $filename could not be found")
    end
end

function __init__()
    # Turn all the callbacks into C callable functions.
    global const on_message_begin_cb = cfunction(on_message_begin, HTTP_CB...)
    global const on_url_cb = cfunction(on_url, HTTP_DATA_CB...)
    global const on_status_complete_cb = cfunction(on_status_complete, HTTP_CB...)
    global const on_header_field_cb = cfunction(on_header_field, HTTP_DATA_CB...)
    global const on_header_value_cb = cfunction(on_header_value, HTTP_DATA_CB...)
    global const on_headers_complete_cb = cfunction(on_headers_complete, HTTP_CB...)
    global const on_body_cb = cfunction(on_body, HTTP_DATA_CB...)
    global const on_message_complete_cb = cfunction(on_message_complete, HTTP_CB...)
end

end # module HttpServer


text:
# Voice Conversion

### basic setup

# Load source speaker's (`clb`) speech signal.
path = joinpath(Pkg.dir("VoiceConversion", "test", "data", "clb_a0028.txt"))
x = vec(readdlm(path))
fs = 16000
period = 5.0
order = 40
alpha = mcepalpha(fs)

f0, timeaxis = dio(x, fs, DioOption(period=period))
f0 = stonemask(x, fs, timeaxis, f0)
spectrogram = cheaptrick(x, fs, timeaxis, f0)
src_clb28 = sp2mc(spectrogram, order, alpha)
@test all(isfinite.(src_clb28))
ap = d4c(x, fs, timeaxis, f0)
@test all(isfinite.(ap))

x_clb28 = copy(x)

# peform conversion and return synthesized waveform
function vc_base(src, mapper)
    converted = vc(mapper, src)
    @test all(isfinite.(converted))
    fftlen = size(spectrogram,1)*2-1
    converted_spectrogram = mc2sp(converted, alpha, fftlen)
    synthesis(f0, converted_spectrogram, ap, period, fs, length(x_clb28))
end

println("testing: voice conversion using the WORLD vocoder.")

let
    println("Female (`clb`) to female (`slt`) voice conversion")
    println("GMM-based frame-by-frame mapping")
    x = copy(src_clb28)

    # Load GMM to convert speech signal of `clb` to that of `slt` and vise versa,
    # mixture: 32, order of mel-cepstrum: 40
    modelpath = joinpath(Pkg.dir("VoiceConversion"), "test", "models",
                         "clb_and_slt_gmm32_order40.jld")
    gmm = load(modelpath)
    @assert !gmm["diff"]

    # Construct GMM-based frame-by-frame mapping
    mapper = GMMMap(gmm["weights"], gmm["means"], gmm["covars"])

    y = vc_base(x, mapper)
    @test all(isfinite.(y))
end

let
    println("Female (`clb`) to female (`slt`) voice conversion")
    println("GMM-based trajectory parameter mapping")
    x = copy(src_clb28)

    # add dynamic feature
    x = [reshape(x[1,:], 1, size(x,2)); push_delta(x[2:end,:])]

    # Load GMM to convert speech signal of `clb` to that of `slt`,
    # mixture: 32, order of mel-cepstrum: 40+40 (with delta feature)
    modelpath = joinpath(Pkg.dir("VoiceConversion"), "test", "models",
                         "clb_and_slt_gmm32_order40_with_delta.jld")
    gmm = load(modelpath)
    @assert !gmm["diff"]

    # Construct trajectory-based GMM parameter mapping
    mapper = GMMMap(gmm["weights"], gmm["means"], gmm["covars"])
    mapper = TrajectoryGMMMap(mapper, 70)
    y = vc_base(x, mapper)
    @test all(isfinite.(y))
end


text:
using ITensors
using ITensorInfiniteMPS
using Random

Random.seed!(1234)

function _siteind(site_tag, n::Int; space)
  return addtags(Index(space, "Site,n=$n"), site_tag)
end
function _siteinds(site_tag, N::Int; space)
  return [_siteind(site_tag, n; space=space) for n in 1:N]
end

# Number of sites in the unit cell
N = 2
space = [QN("SzParity", 1, 2) => 1, QN("SzParity", 0, 2) => 1]
s = _siteinds("S=1/2", N; space=space)
model_kwargs = (J=1.0, h=1.5)
model = Model(:ising)
initstate(n) = isodd(n) ? "↑" : "↓"

# Finite MPS
Nfinite = 100
sfinite = siteinds("S=1/2", Nfinite; conserve_szparity=true)
Hfinite = MPO(model, sfinite; model_kwargs...)
ψfinite = randomMPS(sfinite, initstate)
sweeps = Sweeps(20)
setmaxdim!(sweeps, 1)
setcutoff!(sweeps, 1E-10)
energy_finite, ψfinite = dmrg(Hfinite, ψfinite, sweeps)
@show energy_finite / Nfinite

# Form the Hamiltonian
Σ∞h = InfiniteSum{MPO}(model, s; model_kwargs...)

d = 1
χ = [QN("SzParity", 1, 2) => d, QN("SzParity", 0, 2) => d]
ψ = InfiniteMPS(s; space=χ)
for n in 1:N
  for b in nzblocks(QN(), inds(ψ[n]))
    # Need to overwrite the tensor, since modifying
    # it in-place doesn' work when not modifying the storage directly.
    ψ[n] = insertblock!(ψ[n], b)
  end
end
randn!.(ψ)
ψ = orthogonalize(ψ, :)

# Check translational invariance
@show norm(contract(ψ.AL[1:N]..., ψ.C[N]) - contract(ψ.C[0], ψ.AR[1:N]...))

ψ = vumps(Σ∞h, ψ; niter=30, environment_iterations=20)

# Check translational invariance
@show norm(contract(ψ.AL[1:N]..., ψ.C[N]) - contract(ψ.C[0], ψ.AR[1:N]...))

function energy(ψ1, ψ2, h)
  ϕ = ψ1 * ψ2
  return (noprime(ϕ * h) * dag(ϕ))[]
end

function expect(ψ, o)
  return (noprime(ψ * op(o, filterinds(ψ, "Site")...)) * dag(ψ))[]
end

nfinite = Nfinite ÷ 2
orthogonalize!(ψfinite, nfinite)
hnfinite = ITensor(model, sfinite[nfinite], sfinite[nfinite + 1]; model_kwargs...)
energy_finite = energy(ψfinite[nfinite], ψfinite[nfinite + 1], hnfinite)
energy_infinite = energy(ψ.AL[1], ψ.AL[2] * ψ.C[2], Σ∞h[(1, 2)])
@show energy_finite, energy_infinite
@show abs(energy_finite - energy_infinite)

@show expect(ψfinite[nfinite], "Sz")
@show expect(ψfinite[nfinite + 1], "Sz")
@show expect(ψ[1], "Sz")
@show expect(ψ[2], "Sz")


text:
# calc(rows, cols, colfunc)
#
# calc((:a,:b), r -> (r.x+2, r.y-3))
# calc(:width, r -> r[3] - r[2] + 1)

export mutate, calc

struct Mutate{I,F,C,O}
    rows::I
    func::F
end

"""
    mutate(rows, cols::Tuple, func)
    mutate(rows, col::Symbol, func)
    rows |> mutate(cols::Tuple, func)
    rows |> mutate(col::Symbol, func)

Add or replace columns `cols` computed as `cols = func(row)` to genome ordered stream `rows`.
"""
function mutate(rows, cols::Tuple, func)
    # figure out eltype of merge(r, NamedTuple{cols}(func(r)))
    sch = Tables.schema(rows)
    names = sch.names
    types = sch.types
    name2type = Dict{Symbol, Any}(zip(names, types))
    fnames = cols
    ftypes = Base.return_types(func, (eltype(rows),))[1].parameters
    for (n,t) in zip(fnames, ftypes)
        name2type[n] = t
    end
    onames = (names..., setdiff(fnames, names)...)
    otypes = Tuple((name2type[n] for n in onames))
    O = NamedTuple{onames, Tuple{otypes...}}

    Mutate{typeof(rows), typeof(func), fnames, O}(rows, func)
end

mutate(rows, col::Symbol, func) = mutate(rows, (col,), row -> (func(row),))

mutate(cols::Tuple, func) = rows -> mutate(rows, cols, func)
mutate(col::Symbol, func) = rows -> mutate(rows, col, func)

calc = mutate

Base.IteratorEltype(::Type{Mutate{I,F,C,M}}) where {I,F,C,M} = Base.HasEltype()
Base.eltype(::Type{Mutate{I,F,C,M}}) where {I,F,C,M} = M
Base.IteratorSize(::Type{Mutate{I,F,C,M}}) where {I,F,C,M} = Base.IteratorSize(I)
Base.length(m::Mutate{I,F,C,M}) where {I,F,C,M} = Base.length(m.rows)
#Base.size(m::Mutate{I,F,C,M}) where {I,F,C,M} = Base.size(m.rows)

Tables.istable(m::Mutate{I,F,C,M}) where {I,F,C,M} = true
Tables.rowaccess(m::Mutate{I,F,C,M}) where {I,F,C,M} = true
Tables.schema(m::Mutate{I,F,C,M}) where {I,F,C,M} = Tables.Schema(M)

function Base.iterate(m::Mutate{I,F,C,M}) where {I,F,C,M}
    elt_s = iterate(m.rows)
    state = elt_s
    iterate(m, state)
end

function Base.iterate(m::Mutate{I,F,C,M}, state) where {I,F,C,M}
    elt_s = state
    elt_s == nothing && return nothing
    M(Base.merge(elt_s[1], NamedTuple{C}(m.func(elt_s[1])))), iterate(m.rows, elt_s[2])
end


text:
# This needs renaming to avoid collision with the DecionModelResults/EmulationModelResults
mutable struct ProblemResults <: IS.Results
    base_power::Float64
    timestamps::StepRange{Dates.DateTime, Dates.Millisecond}
    system::Union{Nothing, PSY.System}
    system_uuid::Base.UUID
    aux_variable_values::Dict{AuxVarKey, DataFrames.DataFrame}
    variable_values::Dict{VariableKey, DataFrames.DataFrame}
    dual_values::Dict{ConstraintKey, DataFrames.DataFrame}
    parameter_values::Dict{ParameterKey, DataFrames.DataFrame}
    expression_values::Dict{ExpressionKey, DataFrames.DataFrame}
    optimizer_stats::DataFrames.DataFrame
    optimization_container_metadata::OptimizationContainerMetadata
    model_type::String
    output_dir::String
end

list_aux_variable_keys(res::ProblemResults) = collect(keys(res.aux_variable_values))
list_aux_variable_names(res::ProblemResults) =
    encode_keys_as_strings(keys(res.aux_variable_values))
list_variable_keys(res::ProblemResults) = collect(keys(res.variable_values))
list_variable_names(res::ProblemResults) = encode_keys_as_strings(keys(res.variable_values))
list_parameter_keys(res::ProblemResults) = collect(keys(res.parameter_values))
list_parameter_names(res::ProblemResults) =
    encode_keys_as_strings(keys(res.parameter_values))
list_dual_keys(res::ProblemResults) = collect(keys(res.dual_values))
list_dual_names(res::ProblemResults) = encode_keys_as_strings(keys(res.dual_values))
list_expression_keys(res::ProblemResults) = collect(keys(res.expression_values))
list_expression_names(res::ProblemResults) =
    encode_keys_as_strings(keys(res.expression_values))
get_timestamps(res::ProblemResults) = res.timestamps
get_model_base_power(res::ProblemResults) = res.base_power
get_dual_values(res::ProblemResults) = res.dual_values
get_expression_values(res::ProblemResults) = res.expression_values
get_variable_values(res::ProblemResults) = res.variable_values
get_aux_variable_values(res::ProblemResults) = res.aux_variable_values
get_total_cost(res::ProblemResults) = get_objective_value(res)
get_optimizer_stats(res::ProblemResults) = res.optimizer_stats
get_parameter_values(res::ProblemResults) = res.parameter_values
get_resolution(res::ProblemResults) = res.timestamps.step
get_system(res::ProblemResults) = res.system
get_forecast_horizon(res::ProblemResults) = length(get_timestamps(res))

function get_objective_value(res::ProblemResults, execution=1)
    return res.optimizer_stats[execution, :objective_value]
end

"""
Construct ProblemResults from a solved DecisionModel.
"""
function ProblemResults(model::DecisionModel)
    status = get_run_status(model)
    status != RunStatus.SUCCESSFUL && error("problem was not solved successfully: $status")

    model_store = get_store(model)

    if isempty(model_store)
        error("Model Solved as part of a Simulation.")
    end

    timestamps = get_timestamps(model)
    optimizer_stats = to_dataframe(get_optimizer_stats(model))

    aux_variable_values =
        Dict(x => read_aux_variable(model, x) for x in list_aux_variable_keys(model))
    variable_values = Dict(x => read_variable(model, x) for x in list_variable_keys(model))
    dual_values = Dict(x => read_dual(model, x) for x in list_dual_keys(model))
    parameter_values =
        Dict(x => read_parameter(model, x) for x in list_parameter_keys(model))
    expression_values =
        Dict(x => read_expression(model, x) for x in list_expression_keys(model))

    sys = get_system(model)

    return ProblemResults(
        get_problem_base_power(model),
        timestamps,
        sys,
        IS.get_uuid(sys),
        aux_variable_values,
        variable_values,
        dual_values,
        parameter_values,
        expression_values,
        optimizer_stats,
        get_metadata(get_optimization_container(model)),
        IS.strip_module_name(typeof(model)),
        mkpath(joinpath(get_output_dir(model), "results")),
    )
end

"""
Construct ProblemResults from a solved EmulationModel.
"""
function ProblemResults(model::EmulationModel)
    status = get_run_status(model)
    status != RunStatus.SUCCESSFUL && error("problem was not solved successfully: $status")

    model_store = get_store(model)

    if isempty(model_store)
        error("Model Solved as part of a Simulation.")
    end

    aux_variables =
        Dict(x => read_aux_variable(model, x) for x in list_aux_variable_keys(model))
    variables = Dict(x => read_variable(model, x) for x in list_variable_keys(model))
    duals = Dict(x => read_dual(model, x) for x in list_dual_keys(model))
    parameters = Dict(x => read_parameter(model, x) for x in list_parameter_keys(model))
    expression = Dict(x => read_expression(model, x) for x in list_expression_keys(model))
    optimizer_stats = read_optimizer_stats(model)
    initial_time = get_initial_time(model)
    container = get_optimization_container(model)
    sys = get_system(model)

    return ProblemResults(
        get_problem_base_power(model),
        StepRange(initial_time, get_resolution(model), initial_time),
        sys,
        IS.get_uuid(sys),
        aux_variables,
        variables,
        duals,
        parameters,
        expression,
        optimizer_stats,
        get_metadata(container),
        IS.strip_module_name(typeof(model)),
        mkpath(joinpath(get_output_dir(model), "results")),
    )
end

"""
Exports all results from the operations problem.
"""
function export_results(results::ProblemResults; kwargs...)
    exports = ProblemResultsExport(
        "Problem",
        store_all_duals=true,
        store_all_parameters=true,
        store_all_variables=true,
        store_all_aux_variables=true,
    )
    return export_results(results, exports; kwargs...)
end

function export_results(
    results::ProblemResults,
    exports::ProblemResultsExport;
    file_type=CSV.File,
)
    file_type != CSV.File && error("only CSV.File is currently supported")
    export_path = mkpath(joinpath(results.output_dir, "variables"))
    for (key, df) in results.variable_values
        if should_export_variable(exports, key)
            export_result(file_type, export_path, key, df)
        end
    end

    export_path = mkpath(joinpath(results.output_dir, "aux_variables"))
    for (key, df) in results.aux_variable_values
        if should_export_aux_variable(exports, key)
            export_result(file_type, export_path, key, df)
        end
    end

    export_path = mkpath(joinpath(results.output_dir, "duals"))
    for (key, df) in results.dual_values
        if should_export_dual(exports, key)
            export_result(file_type, export_path, key, df)
        end
    end

    export_path = mkpath(joinpath(results.output_dir, "parameters"))
    for (key, df) in results.parameter_values
        if should_export_parameter(exports, key)
            export_result(file_type, export_path, key, df)
        end
    end

    export_path = mkpath(joinpath(results.output_dir, "expressions"))
    for (key, df) in results.expression_values
        if should_export_expression(exports, key)
            export_result(file_type, export_path, key, df)
        end
    end

    if exports.optimizer_stats
        export_result(
            file_type,
            joinpath(results.output_dir, "optimizer_stats.csv"),
            results.optimizer_stats,
        )
    end

    @info "Exported ProblemResults to $(results.output_dir)"
end

function _deserialize_key(
    ::Type{<:OptimizationContainerKey},
    results::ProblemResults,
    name::AbstractString,
)
    return deserialize_key(results.optimization_container_metadata, name)
end

function _deserialize_key(
    ::Type{T},
    ::ProblemResults,
    args...,
) where {T <: OptimizationContainerKey}
    return make_key(T, args...)
end

read_optimizer_stats(res::ProblemResults) = res.optimizer_stats

"""
Set the system in the results instance.

Throws InvalidValue if the system UUID is incorrect.
"""
function set_system!(res::ProblemResults, system::PSY.System)
    sys_uuid = IS.get_uuid(system)
    if sys_uuid != res.system_uuid
        throw(
            IS.InvalidValue(
                "System mismatch. $sys_uuid does not match the stored value of $(res.system_uuid)",
            ),
        )
    end

    res.system = system
    return
end

const _PROBLEM_RESULTS_FILENAME = "problem_results.bin"

"""
Serialize the results to a binary file.

It is recommended that `directory` be the directory that contains a serialized
OperationModel. That will allow automatic deserialization of the PowerSystems.System.
The `ProblemResults` instance can be deserialized with `ProblemResults(directory)`.
"""
function serialize_results(res::ProblemResults, directory::AbstractString)
    mkpath(directory)
    filename = joinpath(directory, _PROBLEM_RESULTS_FILENAME)
    isfile(filename) && rm(filename)
    Serialization.serialize(filename, _copy_for_serialization(res))
    @info "Serialize ProblemResults to $filename"
end

"""
Construct a ProblemResults instance from a serialized directory.

If the directory contains a serialized PowerSystems.System then it will deserialize that
system and add it to the results. Otherwise, it is up to the caller to call
[`set_system!`](@ref) on the returned instance to restore it.
"""
function ProblemResults(directory::AbstractString)
    filename = joinpath(directory, _PROBLEM_RESULTS_FILENAME)
    if !isfile(filename)
        error("No results file exists in $directory")
    end

    results = Serialization.deserialize(filename)
    possible_sys_file = joinpath(directory, make_system_filename(results.system_uuid))
    if isfile(possible_sys_file)
        set_system!(results, PSY.System(possible_sys_file))
    else
        @info "$directory does not contain a serialized System, skipping deserialization."
    end

    return results
end

function _copy_for_serialization(res::ProblemResults)
    return ProblemResults(
        res.base_power,
        res.timestamps,
        nothing,
        res.system_uuid,
        res.aux_variable_values,
        res.variable_values,
        res.dual_values,
        res.parameter_values,
        res.expression_values,
        res.optimizer_stats,
        res.optimization_container_metadata,
        res.model_type,
        res.output_dir,
    )
end

function _read_results(
    result_values::Dict{<:OptimizationContainerKey, DataFrames.DataFrame},
    container_keys,
    timestamps,
    time_ids,
    base_power,
)
    existing_keys = keys(result_values)
    container_keys = container_keys === nothing ? existing_keys : container_keys
    _validate_keys(existing_keys, container_keys)
    results = Dict{OptimizationContainerKey, DataFrames.DataFrame}()
    for (k, v) in result_values
        if k in container_keys
            results[k] =
                convert_result_to_natural_units(k) ? v[time_ids, :] .* base_power :
                v[time_ids, :]
            DataFrames.insertcols!(results[k], 1, :DateTime => timestamps)
        end
    end
    return results
end

function _process_timestamps(
    res::ProblemResults,
    start_time::Union{Nothing, Dates.DateTime},
    len::Union{Int, Nothing},
)
    if start_time === nothing
        start_time = first(get_timestamps(res))
    elseif start_time ∉ get_timestamps(res)
        throw(IS.InvalidValue("start_time not in result timestamps"))
    end

    if startswith(res.model_type, "EmulationModel{")
        def_len = DataFrames.nrow(get_optimizer_stats(res))
        requested_range =
            collect(findfirst(x -> x >= start_time, get_timestamps(res)):def_len)
        timestamps = repeat(get_timestamps(res), def_len)
    else
        timestamps = get_timestamps(res)
        requested_range = findall(x -> x >= start_time, timestamps)
        def_len = length(requested_range)
    end
    len = len === nothing ? def_len : len
    if len > def_len
        throw(IS.InvalidValue("requested results have less than $len values"))
    end
    timestamp_ids = requested_range[1:len]
    return timestamp_ids, timestamps[timestamp_ids]
end

"""
Return the values for the requested variable key for a problem.
Accepts a vector of keys for the return of the values. If the time stamps and keys are
loaded using the [load_results!](@ref) function it will read from memory.

# Arguments

  - `variable::Tuple{Type{<:VariableType}, Type{<:PSY.Component}` : Tuple with variable type and device type for the desired results
  - `start_time::Dates.DateTime` : start time of the requested results
  - `len::Int`: length of results
"""
function read_variable(res::ProblemResults, args...; kwargs...)
    key = VariableKey(args...)
    return read_variable(res, key; kwargs...)
end

function read_variable(res::ProblemResults, key::AbstractString; kwargs...)
    return read_variable(res, _deserialize_key(VariableKey, res, key); kwargs...)
end

function read_variable(
    res::ProblemResults,
    key::VariableKey;
    start_time::Union{Nothing, Dates.DateTime}=nothing,
    len::Union{Int, Nothing}=nothing,
)
    return read_variables_with_keys(res, [key]; start_time=start_time, len=len)[key]
end

"""
Return the values for the requested variable keys for a problem.
Accepts a vector of keys for the return of the values. If the time stamps and keys are
loaded using the [load_results!](@ref) function it will read from memory.

# Arguments

  - `variables::Vector{Tuple{Type{<:VariableType}, Type{<:PSY.Component}}` : Tuple with variable type and device type for the desired results
  - `start_time::Dates.DateTime` : initial time of the requested results
  - `len::Int`: length of results
"""
function read_variables(res::ProblemResults, variables; kwargs...)
    return read_variables(res, [VariableKey(x...) for x in variables]; kwargs...)
end

function read_variables(res::ProblemResults, variables::Vector{<:AbstractString}; kwargs...)
    return read_variables(
        res,
        [_deserialize_key(VariableKey, res, x) for x in variables];
        kwargs...,
    )
end

function read_variables(
    res::ProblemResults,
    variables::Vector{<:OptimizationContainerKey};
    start_time::Union{Nothing, Dates.DateTime}=nothing,
    len::Union{Int, Nothing}=nothing,
)
    result_values = read_variables_with_keys(res, variables; start_time=start_time, len=len)
    return Dict(encode_key_as_string(k) => v for (k, v) in result_values)
end

"""
Return the values for all variables.
"""
function read_variables(res::IS.Results)
    variables = Dict(x => read_variable(res, x) for x in list_variable_names(res))
end

function read_variables_with_keys(
    res::ProblemResults,
    variables::Vector{<:OptimizationContainerKey};
    start_time::Union{Nothing, Dates.DateTime}=nothing,
    len::Union{Int, Nothing}=nothing,
)
    (timestamp_ids, timestamps) = _process_timestamps(res, start_time, len)
    return _read_results(
        res.variable_values,
        variables,
        timestamps,
        timestamp_ids,
        get_model_base_power(res),
    )
end

"""
Return the values for the requested dual key for a problem.
Accepts a vector of keys for the return of the values. If the time stamps and keys are
loaded using the [load_results!](@ref) function it will read from memory.

# Arguments

  - `dual::Tuple{Type{<:ConstraintType}, Type{<:PSY.Component}` : Tuple with dual type and device type for the desired results
  - `start_time::Dates.DateTime` : initial time of the requested results
  - `len::Int`: length of results
"""
function read_dual(res::ProblemResults, args...; kwargs...)
    key = ConstraintKey(args...)
    return read_dual(res, key; kwargs...)
end

function read_dual(res::ProblemResults, key::AbstractString; kwargs...)
    return read_dual(res, _deserialize_key(ConstraintKey, res, key); kwargs...)
end

function read_dual(
    res::ProblemResults,
    key::ConstraintKey;
    start_time::Union{Nothing, Dates.DateTime}=nothing,
    len::Union{Int, Nothing}=nothing,
)
    return read_duals_with_keys(res, [key]; start_time=start_time, len=len)[key]
end

"""
Return the values for the requested dual keys for a problem.
Accepts a vector of keys for the return of the values. If the time stamps and keys are
loaded using the [load_results!](@ref) function it will read from memory.

# Arguments

  - `duals::Vector{Tuple{Type{<:ConstraintType}, Type{<:PSY.Component}}` : Tuple with dual type and device type for the desired results
  - `start_time::Dates.DateTime` : initial time of the requested results
  - `len::Int`: length of results
"""
function read_duals(res::ProblemResults, duals; kwargs...)
    return read_duals(res, [ConstraintKey(x...) for x in duals]; kwargs...)
end

function read_duals(res::ProblemResults, duals::Vector{<:AbstractString}; kwargs...)
    return read_duals(
        res,
        [_deserialize_key(ConstraintKey, res, x) for x in duals];
        kwargs...,
    )
end

function read_duals(
    res::ProblemResults,
    duals::Vector{<:OptimizationContainerKey};
    start_time::Union{Nothing, Dates.DateTime}=nothing,
    len::Union{Int, Nothing}=nothing,
)
    result_values = read_duals_with_keys(res, duals; start_time=start_time, len=len)
    return Dict(encode_key_as_string(k) => v for (k, v) in result_values)
end

"""
Return the values for all duals.
"""
function read_duals(res::IS.Results)
    duals = Dict(x => read_dual(res, x) for x in list_dual_names(res))
end

function read_duals_with_keys(
    res::ProblemResults,
    duals::Vector{<:OptimizationContainerKey};
    start_time::Union{Nothing, Dates.DateTime}=nothing,
    len::Union{Int, Nothing}=nothing,
)
    (timestamp_ids, timestamps) = _process_timestamps(res, start_time, len)
    return _read_results(
        res.dual_values,
        duals,
        timestamps,
        timestamp_ids,
        get_model_base_power(res),
    )
end

"""
Return the values for the requested parameter key for a problem.
Accepts a vector of keys for the return of the values. If the time stamps and keys are
loaded using the [load_results!](@ref) function it will read from memory.

# Arguments

  - `parameter::Tuple{Type{<:ParameterType}, Type{<:PSY.Component}` : Tuple with parameter type and device type for the desired results
  - `start_time::Dates.DateTime` : initial time of the requested results
  - `len::Int`: length of results
"""
function read_parameter(res::ProblemResults, args...; kwargs...)
    key = ParameterKey(args...)
    return read_parameter(res, key; kwargs...)
end

function read_parameter(res::ProblemResults, key::AbstractString; kwargs...)
    return read_parameter(res, _deserialize_key(ParameterKey, res, key); kwargs...)
end

function read_parameter(
    res::ProblemResults,
    key::ParameterKey;
    start_time::Union{Nothing, Dates.DateTime}=nothing,
    len::Union{Int, Nothing}=nothing,
)
    return read_parameters_with_keys(res, [key]; start_time=start_time, len=len)[key]
end

"""
Return the values for the requested parameter keys for a problem.
Accepts a vector of keys for the return of the values. If the time stamps and keys are
loaded using the [load_results!](@ref) function it will read from memory.

# Arguments

  - `parameters::Vector{Tuple{Type{<:ParameterType}, Type{<:PSY.Component}}` : Tuple with parameter type and device type for the desired results
  - `start_time::Dates.DateTime` : initial time of the requested results
  - `len::Int`: length of results
"""
function read_parameters(res::ProblemResults, parameters; kwargs...)
    return read_parameters(res, [ParameterKey(x...) for x in parameters]; kwargs...)
end

function read_parameters(
    res::ProblemResults,
    parameters::Vector{<:AbstractString};
    kwargs...,
)
    return read_parameters(
        res,
        [_deserialize_key(ParameterKey, res, x) for x in parameters];
        kwargs...,
    )
end

function read_parameters(
    res::ProblemResults,
    parameters::Vector{<:OptimizationContainerKey};
    start_time::Union{Nothing, Dates.DateTime}=nothing,
    len::Union{Int, Nothing}=nothing,
)
    result_values =
        read_parameters_with_keys(res, parameters; start_time=start_time, len=len)
    return Dict(encode_key_as_string(k) => v for (k, v) in result_values)
end

"""
Return the values for all parameters.
"""
function read_parameters(res::IS.Results)
    parameters = Dict(x => read_parameter(res, x) for x in list_parameter_names(res))
end

function read_parameters_with_keys(
    res::ProblemResults,
    parameters::Vector{<:OptimizationContainerKey};
    start_time::Union{Nothing, Dates.DateTime}=nothing,
    len::Union{Int, Nothing}=nothing,
)
    (timestamp_ids, timestamps) = _process_timestamps(res, start_time, len)
    return _read_results(
        res.parameter_values,
        parameters,
        timestamps,
        timestamp_ids,
        get_model_base_power(res),
    )
end

"""
Return the values for the requested aux_variable key for a problem.
Accepts a vector of keys for the return of the values. If the time stamps and keys are
loaded using the [load_results!](@ref) function it will read from memory.

# Arguments

  - `aux_variable::Tuple{Type{<:AuxVariableType}, Type{<:PSY.Component}` : Tuple with aux_variable type and device type for the desired results
  - `start_time::Dates.DateTime` : initial time of the requested results
  - `len::Int`: length of results
"""
function read_aux_variable(res::ProblemResults, args...; kwargs...)
    key = AuxVarKey(args...)
    return read_aux_variable(res, key; kwargs...)
end

function read_aux_variable(res::ProblemResults, key::AbstractString; kwargs...)
    return read_aux_variable(res, _deserialize_key(AuxVarKey, res, key); kwargs...)
end

function read_aux_variable(
    res::ProblemResults,
    key::AuxVarKey;
    start_time::Union{Nothing, Dates.DateTime}=nothing,
    len::Union{Int, Nothing}=nothing,
)
    return read_aux_variables_with_keys(res, [key]; start_time=start_time, len=len)[key]
end

"""
Return the values for the requested aux_variable keys for a problem.
Accepts a vector of keys for the return of the values. If the time stamps and keys are
loaded using the [load_results!](@ref) function it will read from memory.

# Arguments

  - `aux_variables::Vector{Tuple{Type{<:AuxVariableType}, Type{<:PSY.Component}}` : Tuple with aux_variable type and device type for the desired results
  - `start_time::Dates.DateTime` : initial time of the requested results
  - `len::Int`: length of results
"""
function read_aux_variables(res::ProblemResults, aux_variables; kwargs...)
    return read_aux_variables(res, [AuxVarKey(x...) for x in aux_variables]; kwargs...)
end

function read_aux_variables(
    res::ProblemResults,
    aux_variables::Vector{<:AbstractString};
    kwargs...,
)
    return read_aux_variables(
        res,
        [_deserialize_key(AuxVarKey, res, x) for x in aux_variables];
        kwargs...,
    )
end

function read_aux_variables(
    res::ProblemResults,
    aux_variables::Vector{<:OptimizationContainerKey};
    start_time::Union{Nothing, Dates.DateTime}=nothing,
    len::Union{Int, Nothing}=nothing,
)
    result_values =
        read_aux_variables_with_keys(res, aux_variables; start_time=start_time, len=len)
    return Dict(encode_key_as_string(k) => v for (k, v) in result_values)
end

"""
Return the values for all auxiliary variables.
"""
function read_aux_variables(res::IS.Results)
    variables = Dict(x => read_aux_variable(res, x) for x in list_aux_variable_names(res))
end

function read_aux_variables_with_keys(
    res::ProblemResults,
    aux_variables::Vector{<:OptimizationContainerKey};
    start_time::Union{Nothing, Dates.DateTime}=nothing,
    len::Union{Int, Nothing}=nothing,
)
    (timestamp_ids, timestamps) = _process_timestamps(res, start_time, len)
    return _read_results(
        res.aux_variable_values,
        aux_variables,
        timestamps,
        timestamp_ids,
        get_model_base_power(res),
    )
end

"""
Return the values for the requested expression key for a problem.
Accepts a vector of keys for the return of the values. If the time stamps and keys are
loaded using the [load_results!](@ref) function it will read from memory.

# Arguments

  - `expression::Tuple{Type{<:ExpressionType}, Type{<:PSY.Component}` : Tuple with expression type and device type for the desired results
  - `start_time::Dates.DateTime` : initial time of the requested results
  - `len::Int`: length of results
"""
function read_expression(res::ProblemResults, args...; kwargs...)
    key = ExpressionKey(args...)
    return read_expression(res, key; kwargs...)
end

function read_expression(res::ProblemResults, key::AbstractString; kwargs...)
    return read_expression(res, _deserialize_key(ExpressionKey, res, key); kwargs...)
end

function read_expression(
    res::ProblemResults,
    key::ExpressionKey;
    start_time::Union{Nothing, Dates.DateTime}=nothing,
    len::Union{Int, Nothing}=nothing,
)
    return read_expressions_with_keys(res, [key]; start_time=start_time, len=len)[key]
end

"""
Return the values for the requested expression keys for a problem.
Accepts a vector of keys for the return of the values. If the time stamps and keys are
loaded using the [load_results!](@ref) function it will read from memory.

# Arguments

  - `expressions::Vector{Tuple{Type{<:ExpressionType}, Type{<:PSY.Component}}` : Tuple with expression type and device type for the desired results
  - `start_time::Dates.DateTime` : initial time of the requested results
  - `len::Int`: length of results
"""
function read_expressions(res::ProblemResults; kwargs...)
    return read_expressions(res, collect(keys(res.expression_values)); kwargs...)
end

function read_expressions(res::ProblemResults, expressions; kwargs...)
    return read_expressions(res, [ExpressionKey(x...) for x in expressions]; kwargs...)
end

function read_expressions(
    res::ProblemResults,
    expressions::Vector{<:AbstractString};
    kwargs...,
)
    return read_expressions(
        res,
        [_deserialize_key(ExpressionKey, res, x) for x in expressions];
        kwargs...,
    )
end

function read_expressions(
    res::ProblemResults,
    expressions::Vector{<:OptimizationContainerKey};
    start_time::Union{Nothing, Dates.DateTime}=nothing,
    len::Union{Int, Nothing}=nothing,
)
    result_values =
        read_expressions_with_keys(res, expressions; start_time=start_time, len=len)
    return Dict(encode_key_as_string(k) => v for (k, v) in result_values)
end

"""
Return the values for all expressions.
"""
function read_expressions(res::IS.Results)
    expressions = Dict(x => read_expression(res, x) for x in list_expression_names(res))
end

function read_expressions_with_keys(
    res::ProblemResults,
    expressions::Vector{<:OptimizationContainerKey};
    start_time::Union{Nothing, Dates.DateTime}=nothing,
    len::Union{Int, Nothing}=nothing,
)
    (timestamp_ids, timestamps) = _process_timestamps(res, start_time, len)
    return _read_results(
        res.expression_values,
        expressions,
        timestamps,
        timestamp_ids,
        get_model_base_power(res),
    )
end

function export_realized_results(res::ProblemResults)
    save_path = mkpath(joinpath(res.output_dir, "export"))
    return export_realized_results(res, save_path)
end


text:
import NSG: make_ref, make_bed
function tst_make_ref()
    lmp = ["dat/2021-05/data/maps/HD.map",
           "dat/2021-05/data/maps/18k_v1.map",
           "dat/2021-05/data/maps/18k_v2.map",
           "dat/2021-05/data/maps/18k_v3.map",
           "dat/2021-05/data/maps/18k_v4.map",
           "dat/2021-05/data/maps/8k_rudi.map"]

    frd = ["dat/2021-05/data/genotypes/600k", # final report directory
           "dat/2021-05/data/genotypes/b17k",
           "dat/2021-05/data/genotypes/a17k",
           "dat/2021-05/data/genotypes/g17k",
           "dat/2021-05/data/genotypes/S18k_v4",
           "dat/2021-05/data/genotypes/7327"]

    tgt = ["dat/2021-05/result/hd",
           "dat/2021-05/result/v1",
           "dat/2021-05/result/v2",
           "dat/2021-05/result/v3",
           "dat/2021-05/result/v4",
           "dat/2021-05/result/8k"]

    for i in 1:6
        # Note, here I used the simplified version of `make_bed`
        # You can later change it to:
        # make_bed(dir, lmap, sid, nf, target)
        make_bed(frd[i], lmp[i], tgt[i])
    end

    make_ref("dat/2021-05/result/v2ref", lmp[3], tgt[1:5]...)
    make_ref("dat/2021-05/result/v4ref", lmp[5], tgt[1:5]...)
end


text:
# Test residual block
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: January 2020

using LinearAlgebra, InvertibleNetworks, Test

# Input
nx = 28
ny = 28
n_in = 4
n_hidden = 8
batchsize = 2
k1 = 3
k2 = 3

# Input
X = glorot_uniform(nx, ny, n_in, batchsize)
X0 = glorot_uniform(nx, ny, n_in, batchsize)
dX = X - X0

# Weights
W1 = glorot_uniform(k1, k1, n_in, n_hidden)
W2 = glorot_uniform(k2, k2, n_hidden, n_hidden)
W3 = glorot_uniform(k1, k1, 2*n_in, n_hidden)
b1 = glorot_uniform(n_hidden)
b2 = glorot_uniform(n_hidden)

W01 = glorot_uniform(k1, k1, n_in, n_hidden)
W02 = glorot_uniform(k2, k2, n_hidden, n_hidden)
W03 = glorot_uniform(k1, k1, 2*n_in, n_hidden)
b01 = glorot_uniform(n_hidden)
b02 = glorot_uniform(n_hidden)

dW1 = W1 - W01
dW2 = W2 - W02
dW3 = W3 - W03
db1 = b1 - b01
db2 = b2 - b02

# Residual blocks
RB = ResidualBlock(W1, W2, W3, b1, b2, nx, ny, batchsize)   # true weights

# Observed data
Y = RB.forward(X)

function loss(RB, X, Y)
    Y_ = RB.forward(X)
    ΔY = Y_ - Y
    f = .5f0*norm(ΔY)^2
    ΔX = RB.backward(ΔY, X)
    return f, ΔX, RB.W1.grad, RB.W2.grad, RB.W3.grad, RB.b1.grad, RB.b2.grad
end

###################################################################################################
# Gradient tests

# Gradient test w.r.t. input
f0, ΔX = loss(RB, X0, Y)[1:2]
h = 0.1f0
maxiter = 5
err1 = zeros(Float32, maxiter)
err2 = zeros(Float32, maxiter)

print("\nGradient test convolutions\n")
for j=1:maxiter
    f = loss(RB, X0 + h*dX, Y)[1]
    err1[j] = abs(f - f0)
    err2[j] = abs(f - f0 - h*dot(dX, ΔX))
    print(err1[j], "; ", err2[j], "\n")
    global h = h/2f0
end

@test isapprox(err1[end] / (err1[1]/2^(maxiter-1)), 1f0; atol=1f1)
@test isapprox(err2[end] / (err2[1]/4^(maxiter-1)), 1f0; atol=1f1)

# Gradient test for weights
RB0 = ResidualBlock(W01, W02, W03, b1, b2, nx, ny, batchsize)   # initial weights
f0, ΔX, ΔW1, ΔW2, ΔW3 = loss(RB0, X, Y)[1:5]
h = 0.1f0
maxiter = 5
err3 = zeros(Float32, maxiter)
err4 = zeros(Float32, maxiter)

print("\nGradient test convolutions\n")
for j=1:maxiter
    RB0.W1.data = W01 + h*dW1
    RB0.W2.data = W02 + h*dW2
    RB0.W3.data = W03 + h*dW3
    f = loss(RB0, X, Y)[1]
    err3[j] = abs(f - f0)
    err4[j] = abs(f - f0 - h*dot(dW1, ΔW1) - h*dot(dW2, ΔW2) - h*dot(dW3, ΔW3))
    print(err3[j], "; ", err4[j], "\n")
    global h = h/2f0
end

@test isapprox(err3[end] / (err3[1]/2^(maxiter-1)), 1f0; atol=1f1)
@test isapprox(err4[end] / (err4[1]/4^(maxiter-1)), 1f0; atol=1f1)

# Gradient test for bias
RB0 = ResidualBlock(W1, W2, W3, b01, b02, nx, ny, batchsize)
f0, ΔX, ΔW1, ΔW2, ΔW3, Δb1, Δb2 = loss(RB0, X, Y)
h = 0.1f0
maxiter = 5
err5 = zeros(Float32, maxiter)
err6 = zeros(Float32, maxiter)
print("\nGradient test convolutions\n")
for j=1:maxiter
    RB0.b1.data = b01 + h*db1
    RB0.b2.data = b02 + h*db2
    f = loss(RB0, X, Y)[1]
    err5[j] = abs(f - f0)
    err6[j] = abs(f - f0 - h*dot(db1, Δb1) - h*dot(db2, Δb2))
    print(err5[j], "; ", err6[j], "\n")
    global h = h/2f0
end

@test isapprox(err5[end] / (err5[1]/2^(maxiter-1)), 1f0; atol=1f1)
@test isapprox(err6[end] / (err6[1]/4^(maxiter-1)), 1f0; atol=1f1)


text:
function applycnsics!(system::CVSystem,old=Dict("a"=>0),restart="no")
    if restart == "no"
        system.cns.H[1] = 1/system.heart.activation.th[1]*minTos;
        system.cns.Emaxlv[1] = system.heart.lv.Emax[1];
        system.cns.Emaxrv[1] = system.heart.rv.Emax[1];

        system.cns.R2L[1] = system.branches.term[13].R[2];
        system.cns.R3L[1] = system.branches.term[13].R[3];
        system.cns.R2U[1] = system.branches.term[54].R[2];
        system.cns.R3U[1] = system.branches.term[54].R[3];
        system.cns.C4L[1] = system.branches.term[13].C[4];
        system.cns.C5L[1] = system.branches.term[13].C[5];
        system.cns.C4U[1] = system.branches.term[54].C[4];
        system.cns.C5U[1] = system.branches.term[54].C[5];
        system.cns.V4L[1] = system.branches.term[13].V0[4];
        system.cns.V5L[1] = system.branches.term[13].V0[5];
        system.cns.V4U[1] = system.branches.term[54].V0[4];
        system.cns.V5U[1] = system.branches.term[54].V0[5];

        system.cns.gammalv = 0.8*system.cns.Emaxlv[1];
        system.cns.alphalv = 1.6*system.cns.Emaxlv[1] - system.cns.gammalv;
        system.cns.gammarv = 0.8*system.cns.Emaxrv[1];
        system.cns.alpharv = 1.6*system.cns.Emaxrv[1] - system.cns.gammarv;
        system.cns.gammaR2L = 0.6*system.cns.R2L[1];
        system.cns.alphaR2L = 2.2*system.cns.R2L[1] - system.cns.gammaR2L;
        system.cns.gammaR3L = 0.6*system.cns.R3L[1];
        system.cns.alphaR3L = 2.2*system.cns.R3L[1] - system.cns.gammaR3L;
        system.cns.gammaR2U = 0.6*system.cns.R2U[1];
        system.cns.alphaR2U = 2.2*system.cns.R2U[1] - system.cns.gammaR2U;
        system.cns.gammaR3U = 0.6*system.cns.R3U[1];
        system.cns.alphaR3U = 2.2*system.cns.R3U[1] - system.cns.gammaR3U;
        system.cns.gammaC4L = 1.1*system.cns.C4L[1];
        system.cns.alphaC4L = -0.7*system.cns.C4L[1] + system.cns.gammaC4L;
        system.cns.gammaC5L = 1.1*system.cns.C5L[1];
        system.cns.alphaC5L = -0.7*system.cns.C5L[1] + system.cns.gammaC5L;
        system.cns.gammaC4U = 1.1*system.cns.C4U[1];
        system.cns.alphaC4U = -0.7*system.cns.C4U[1] + system.cns.gammaC4U;
        system.cns.gammaC5U = 1.1*system.cns.C5U[1];
        system.cns.alphaC5U = -0.7*system.cns.C5U[1] + system.cns.gammaC5U;
        system.cns.gammaV4L = 1.05*system.cns.V4L[1];
        system.cns.alphaV4L = -0.85*system.cns.V4L[1] + system.cns.gammaV4L;
        system.cns.gammaV5L = 1.05*system.cns.V5L[1];
        system.cns.alphaV5L = -0.85*system.cns.V5L[1] + system.cns.gammaV5L;
        system.cns.gammaV4U = 1.05*system.cns.V4U[1];
        system.cns.alphaV4U = -0.85*system.cns.V4U[1] + system.cns.gammaV4U;
        system.cns.gammaV5U = 1.05*system.cns.V5U[1];
        system.cns.alphaV5U = -0.85*system.cns.V5U[1] + system.cns.gammaV5U;

        system.cns.ns = [0.25];
        system.cns.np = [0.25];
    elseif restart == "yes"
        system.cns.H[1] = old["H"][end];
        system.cns.Emaxlv[1] = old["Emaxlv"][end];
        system.cns.Emaxrv[1] = old["Emaxrv"][end];

        system.cns.R2L[1] = old["R2L"][end];
        system.cns.R3L[1] = old["R3L"][end];
        system.cns.R2U[1] = old["R2U"][end];
        system.cns.R3U[1] = old["R3U"][end];
        system.cns.C4L[1] = old["C4L"][end];
        system.cns.C5L[1] = old["C5L"][end];
        system.cns.C4U[1] = old["C4U"][end];
        system.cns.C5U[1] = old["C5U"][end];
        system.cns.V4L[1] = old["V4L"][end];
        system.cns.V5L[1] = old["V5L"][end];
        system.cns.V4U[1] = old["V4U"][end];
        system.cns.V5U[1] = old["V5U"][end];

        system.cns.gammalv = old["gammalv"];
        system.cns.alphalv = old["alphalv"];
        system.cns.gammarv = old["gammarv"];
        system.cns.alpharv = old["alpharv"];
        system.cns.gammaR2L = old["gammaR2L"];
        system.cns.alphaR2L = old["alphaR2L"];
        system.cns.gammaR3L = old["gammaR3L"];
        system.cns.alphaR3L = old["alphaR3L"];
        system.cns.gammaR2U = old["gammaR2U"];
        system.cns.alphaR2U = old["alphaR2U"];
        system.cns.gammaR3U = old["gammaR3U"];
        system.cns.alphaR3U = old["alphaR3U"];
        system.cns.gammaC4L = old["gammaC4L"];
        system.cns.alphaC4L = old["alphaC4L"];
        system.cns.gammaC5L = old["gammaC5L"];
        system.cns.alphaC5L = old["alphaC5L"];
        system.cns.gammaC4U = old["gammaC4U"];
        system.cns.alphaC4U = old["alphaC4U"];
        system.cns.gammaC5U = old["gammaC5U"];
        system.cns.alphaC5U = old["alphaC5U"];
        system.cns.gammaV4L = old["gammaV4L"];
        system.cns.alphaV4L = old["alphaV4L"];
        system.cns.gammaV5L = old["gammaV5L"];
        system.cns.alphaV5L = old["alphaV5L"];
        system.cns.gammaV4U = old["gammaV4U"];
        system.cns.alphaV4U = old["alphaV4U"];
        system.cns.gammaV5U = old["gammaV5U"];
        system.cns.alphaV5U = old["alphaV5U"];

        system.cns.ns = [old["ns"][end]];
        system.cns.np = [old["np"][end]];
    end
    system.cns.alphaH = 2.1563;
    system.cns.betaH = 0.5313;
    system.cns.gammaH = 0.8438;
end

max_stars_repo_name: dcanuto/CardioModeling

max_issues_repo_name: dcanuto/CardioModeling

max_forks_repo_name: dcanuto/CardioModeling


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Functions to obtain a view of a text.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export textview

"""
    textview([io::IO, ][text::AbstractString | lines::Vector{AbstractString}], view::NTuple{4, Int}; kwargs...)

Create a view of `text` or `lines` considering a `view` configuration. The
latter is a tuple with four integers that has the following meaning:

- Top line;
- Number of lines;
- Left column; and
- Number of columns.

If a value equal or lower than 0 is passed to any of those options, then its
extreme value is used.

# Keywords

- `active_match::Int`: The match number that is considered active. This match is
    highlighted using `active_highlight` instead of `highlight`.
    (**Default** = 0)
- `highlight::String`: ANSI escape sequence that contains the decoration of the
    highlight. (**Default** = `\\e[7m`)
- `active_highlight::String`: ANSI escape sequence that contains the decoration
    of the active highlight. (**Default** = `\\e[30;43m`.)
- `frozen_lines_at_beginning::Int`: Number of frozen lines that are drawn in the
    beginning. (**Default** = 0)
- `frozen_columns_at_beginning::Int`: Number of frozen columns that are drawn in
    the beginning. (**Default** = 0)

If `text::AbstractString` is passed, then the following keyword is available:

- `search_regex::Uniont{Nothing, Regex}`: A regex used to highlight matches in
    the text view. (**Default** = `nothing`).

If `lines::Vector{AbstractString}` is passed, then the following keyword is
available:

- `search_matches::Union{Nothing, Dict{Int, Vector{Tuple{Int, Int}}}}`: The
    search matches that are highlighted in the text view. This dictionary must
    be created using the function [`string_search_per_line`](@ref).
    (**Default** = `nothing`)

# Returns

If `io` is passed, the view is written to it and the function returns:

- The number of cropped lines at the end;
- The maximum number cropped characters in a row.

However, if `io` is not passed, the function returns:

- The string with the required view;
- The number of cropped lines at the end;
- The maximum number cropped characters in a row.

!!! note
    If only frozen lines are printed, the second returned value is set to 0.

    If only frozen columns are printed, the third returned value is set to 0.
"""
function textview(
    text::AbstractString,
    view::NTuple{4, Int};
    kwargs...
)
    buf = IOBuffer()
    num_cropped_lines, max_cropped_chars = textview(buf, text, view; kwargs...)
    return String(take!(buf)), num_cropped_lines, max_cropped_chars
end

function textview(
    buf::IO,
    text::AbstractString,
    view::NTuple{4, Int};
    search_regex::Union{Nothing, Regex} = nothing,
    kwargs...
)
    lines = split(text, '\n')

    search_matches = !isnothing(search_regex) ?
        string_search_per_line(lines, search_regex) :
        nothing

    return textview(
        buf,
        lines,
        view;
        search_matches,
        kwargs...
    )
end

function textview(
    lines::Vector{T},
    view::NTuple{4, Int};
    kwargs...
) where T<:AbstractString
    buf = IOBuffer()
    num_cropped_lines, max_cropped_chars = textview(buf, lines, view; kwargs...)
    return String(take!(buf)), num_cropped_lines, max_cropped_chars
end

function textview(
    buf::IO,
    lines::Vector{T},
    view::NTuple{4, Int};
    active_highlight::String = _CSI * "30;43m",
    active_match::Int = 0,
    frozen_columns_at_beginning::Int = 0,
    frozen_lines_at_beginning::Int = 0,
    highlight::String = _CSI * "7m",
    maximum_number_of_columns::Int = -1,
    maximum_number_of_lines::Int = -1,
    parse_decorations_before_view::Bool = false,
    ruler_decoration::String = _CSI * "90m",
    search_matches::Union{Nothing, Dict{Int, Vector{Tuple{Int, Int}}}} = nothing,
    show_ruler::Bool = false,
    title_lines::Int = 0
) where T<:AbstractString

    # Verification of the input parameters
    # ==========================================================================

    start_line  = view[1]
    num_lines   = view[2]
    total_lines = length(lines)

    start_line = clamp(start_line, 1, total_lines)

    if num_lines ≥ 0
        num_lines = clamp(num_lines, 0, total_lines - start_line + 1)
    else
        num_lines = total_lines
    end

    if start_line + num_lines - 1 > total_lines
        num_lines = total_lines - start_line + 1
    end

    if frozen_lines_at_beginning > 0
        start_line = clamp(start_line, frozen_lines_at_beginning + 1, total_lines)
    end

    start_column = view[3]

    if start_column < 1
        start_column = 1
    end

    num_columns = view[4]

    if frozen_columns_at_beginning > 0
        if start_column ≤ frozen_columns_at_beginning
            start_column = frozen_columns_at_beginning + 1
        end
    end

    # If the user wants a ruler, compute the required size here.
    ruler_spacing = 0

    if show_ruler
        ruler_spacing = floor(Int, total_lines |> abs |> log10) + 1

        # If the user selected a maximum number of columns, we need to decrease
        # it to take into account the ruler.
        if maximum_number_of_columns ≥ 0
            maximum_number_of_columns = maximum_number_of_columns - ruler_spacing - 3

            if maximum_number_of_columns < 0
                maximum_number_of_columns = 0
            end
        end
    end

    # Internal variables
    # ==========================================================================

    # Buffers used to store the entire the line view rendering, created here to
    # reduce the number of allocations.
    line_buf = IOBuffer()

    # Count how many matches we passed in the current line that is being
    # processed.
    num_matches = 0

    # Variable to store the maximum number of cropped lines at the end.
    num_cropped_lines_at_end = 0

    # Variable to compute the maximum number of characters cropped in the right.
    max_cropped_chars = 0

    # Variable to store the decorations before the view, including those in the
    # frozen lines. It is used if the option `parse_decorations_before_view`
    # is `true`.
    pre_decorations = ""

    # Check if we have a maximum number of lines.
    if maximum_number_of_lines ≥ 0
        if frozen_lines_at_beginning ≥ maximum_number_of_lines
            frozen_lines_at_beginning = maximum_number_of_lines
            start_lines = frozen_lines_at_beginning + 1
            num_lines = 0

        else
            if num_lines > maximum_number_of_lines - frozen_lines_at_beginning
                num_lines = maximum_number_of_lines - frozen_lines_at_beginning
            end

            num_cropped_lines_at_end = total_lines - start_line - (num_lines - 1)
        end
    else
       num_cropped_lines_at_end = total_lines - start_line - (num_lines - 1)
    end

    # Check if we have a maximum number of columns.
    if maximum_number_of_columns ≥ 0
        # Check if we can only draw the frozen columns.
        if frozen_columns_at_beginning ≥ maximum_number_of_columns
            frozen_columns_at_beginning = maximum_number_of_columns
            start_column = frozen_columns_at_beginning + 1
            num_columns = 0

        else
            if num_columns < 0
                num_columns =
                    maximum_number_of_columns - frozen_columns_at_beginning
            else
                num_columns = clamp(
                    num_columns,
                    -1,
                    maximum_number_of_columns - frozen_columns_at_beginning
                )
            end
        end
    end

    # Frozen lines
    # ==========================================================================

    for l in 1:frozen_lines_at_beginning
        line_active_match = active_match - num_matches

        if !isnothing(search_matches) && haskey(search_matches, l)
            line_search_matches = search_matches[l]
        else
            line_search_matches = nothing
        end

        if show_ruler
            line_number_str = lpad(l, ruler_spacing)
            write(buf, ruler_decoration, " ")
            write(buf, line_number_str, " │")
            write(buf, _CSI, "0m")
        end

        cropped_chars_in_line = 0

        if l ≤ title_lines
            _draw_line_view!(
                buf,
                line_buf,
                lines[l],
                line_search_matches,
                line_active_match,
                highlight,
                active_highlight,
                1,
                0,
                num_columns + frozen_columns_at_beginning
            )

        else
            cropped_chars_in_line = _draw_line_view!(
                buf,
                line_buf,
                lines[l],
                line_search_matches,
                line_active_match,
                highlight,
                active_highlight,
                start_column,
                num_columns,
                frozen_columns_at_beginning
            )
        end

        # We should not compute the number of cropped chars if we are only
        # printing frozen columns.
        if frozen_columns_at_beginning != maximum_number_of_columns
            max_cropped_chars = max(max_cropped_chars, cropped_chars_in_line)
        end

        # At the last frozen line, we must reset all the decorations.
        if l != frozen_lines_at_beginning
            write(buf, '\n')
        else
            write(buf, _CSI, "0m")
            num_lines > 0 && write(buf, '\n')
        end

        if !isnothing(line_search_matches)
            num_matches += length(line_search_matches)
        end

        # If the user wants, we need to accumulate all the decorations from the
        # beginning of the text up to the first line in the view. Here, we
        # accumulate those related to the frozen lines.
        if parse_decorations_before_view
            pre_decorations *= get_decorations(lines[l])
        end
    end

    for l = (frozen_lines_at_beginning + 1):(start_line - 1)
        # Sum the number of matches between the frozen line and the displayed
        # line. This computation is important to find which match is active.
        if !isnothing(search_matches) && haskey(search_matches, l)
            num_matches += length(search_matches[l])
        end

        # If we need to parse the decorations before the view, then obtain the
        # decorations of the current hidden line, and merge with the decorations
        # of the other lines.
        if parse_decorations_before_view
            pre_decorations *= get_decorations(lines[l])
        end
    end

    if parse_decorations_before_view
        d = parse_decoration(pre_decorations)

        # If the pre_decoration is a reset, then we just need to reinitialize
        # it since the reset escape sequence was already written to the buffer.
        !d.reset && write(buf, d |> String)
    end

    # Line views
    # ==========================================================================

    for k in 1:num_lines
        # Get the current line number.
        l = start_line + (k - 1)

        line_active_match = active_match - num_matches

        if !isnothing(search_matches) && haskey(search_matches, l)
            line_search_matches = search_matches[l]
        else
            line_search_matches = nothing
        end

        if show_ruler
            line_number_str = lpad(l, ruler_spacing)
            write(buf, ruler_decoration, " ")
            write(buf, line_number_str, " │")
            write(buf, _CSI, "0m")
        end

        cropped_chars_in_line = _draw_line_view!(
            buf,
            line_buf,
            lines[l],
            line_search_matches,
            line_active_match,
            highlight,
            active_highlight,
            start_column,
            num_columns,
            frozen_columns_at_beginning
        )

        # We should not compute the number of cropped chars if we are only
        # printing frozen columns.
        if frozen_columns_at_beginning != maximum_number_of_columns
            max_cropped_chars = max(max_cropped_chars, cropped_chars_in_line)
        end

        k != num_lines && write(buf, '\n')

        if !isnothing(line_search_matches)
            num_matches += length(line_search_matches)
        end
    end

    return num_cropped_lines_at_end, max_cropped_chars
end

################################################################################
#                              Private functions
################################################################################

# Draw a line view and return the number of right characters that was cropped.
function _draw_line_view!(
    buf::IO,
    line_buf::IOBuffer,
    line::AbstractString,
    line_search_matches::Union{Nothing, Vector{Tuple{Int, Int}}},
    line_active_match::Int,
    highlight::String,
    active_highlight::String,
    start_column::Int,
    num_columns::Int,
    frozen_columns_at_beginning::Int
)

    # Check if we need to highlight the current line.
    if !isnothing(line_search_matches)
        write(line_buf, highlight_search(
            line,
            line_search_matches;
            active_highlight,
            active_match = line_active_match,
            highlight
        ))
    else
        write(line_buf, line)
    end

    line_str = String(take!(line_buf))

    # Frozen columns.
    if frozen_columns_at_beginning > 0
        left, frozen_str = split_string(line_str, frozen_columns_at_beginning)
        write(buf, left, _RESET_DECORATIONS)
    end

    left_ansi  = ""
    right_ansi = ""

    if start_column > 0
        # TODO: Can we improve this?
        # Maybe we can improve the performance by creating a function that
        # remove a certain number of characters from the left/right and also
        # returns the ANSI escape sequence.
        left, line_str = split_string(line_str, start_column - 1)

        # Here we simplify the decorations to avoid too many escape sequences.
        left_ansi = get_decorations(left) |> parse_decoration |> String
    end

    if num_columns ≥ 0
        line_str, right = split_string(line_str, num_columns)

        # Here we simplify the decorations to avoid too many escape sequences.
        right_ansi = get_decorations(right) |> parse_decoration |> String

        # Compute the amount of printable cropped characters in the right
        # string.
        cropped_chars = printable_textwidth(right)
    else
        cropped_chars = 0
    end

    # Write to the buffer.
    write(buf, left_ansi)
    write(buf, line_str)
    write(buf, right_ansi)

    return cropped_chars
end


text:
export
    PopulationManager,
    release!,
    createCandidate,
    randomize!,
    chooseParentIndex,
    keepElites!,
    set!

type PopulationManager
    population::Population
    populationSize::Int
    dnaSize::Int
    best::BestEntity
end

PopulationManager(ctx, populationSize::Int, dnaSize::Int, comparer::AbstractComparer, decode::Function) =
    PopulationManager(Population(ctx, comparer, decode), populationSize, dnaSize, BestEntity(comparer))

release!(popMan::PopulationManager) = (release!(popMan.population); release!(popMan.best))

createCandidate(popMan::PopulationManager) =
    Population(popMan.population.ctx, popMan.population.comparer, popMan.population.decode)

function randomize!(popMan::PopulationManager)
    set!(popMan.best, randomize!(popMan.population, popMan.populationSize, popMan.dnaSize))
end

function chooseParentIndices(popMan::PopulationManager, stdDev::Float32)
    idx1 = chooseParentIndex(popMan, stdDev)
    idx2 = chooseParentIndex(popMan, stdDev)
    while idx1 == idx2
        idx2 = chooseParentIndex(popMan, stdDev)
    end
    idx1, idx2
end

chooseParentIndex(popMan::PopulationManager, stdDev::Float32) =
    chooseParentIndex(stdDev, length(popMan.population))

function chooseParentIndex(stdDev::Float32, size::Int)
    v = abs(randn() * stdDev)
    v = v - floor(v)
    Int(floor(v * size)) + 1
end

function keepElites!(popMan::PopulationManager, candidatePop::Population, rate::Float32)
    keep = Int(rate < 1.0 ? round(length(popMan.populationSize) * rate) : round(rate))
    for i in 1:keep
        push!(candidatePop, copy(popMan.population[i]))
    end
end

function set!(popMan::PopulationManager, candidatePop::Population, noSort = false)
    if noSort == false # sort
        update!(popMan.best, sort!(candidatePop))
    else
        bestOf = candidatePop[1]
        for entity in candidatePop.entities # TODO: support iteration in Population
            if entity != bestOf && popMan.population.comparer(entity.body, bestOf.body)
                bestOf = entity
            end
        end
        update!(popMan.best, bestOf)
    end
    release!(popMan.population)
    popMan.population = candidatePop
end

max_stars_repo_name: unbornchikken/julia-ml-proto

max_issues_repo_name: unbornchikken/julia-ml-proto

max_forks_repo_name: unbornchikken/julia-ml-proto


text:
module test_julia_event

using Test

tc = Ref(0)
t = Timer(0; interval=1) do t
    tc[] += 1
#    @info :tc tc[]
    tc[] >= 5 && close(t)
#    @test isopen(t)
end
@test tc[] == 0
@test isopen(t)

end # module test_julia_event


text:
using CompScienceMeshes, BEAST

trc = X->ntrace(X,γ)
dvg = divergence

h1 = h2 = h3 = 1/15; h = max(h1, h2, h3)
width, height = 1.0, 0.5
Γ1 = meshrectangle(width, height, h1)
Γ2 = meshrectangle(width, height, h2); CompScienceMeshes.rotate!(Γ2, 0.5π*x̂)
Γ3 = meshrectangle(width, height, h3); CompScienceMeshes.rotate!(Γ3, 1.0π*x̂)
γ = meshsegment(width, width, 3)

κ = 1.0
S, T, I = SingleLayerTrace(κ*im), MWSingleLayer3D(κ), Identity()
St = transpose(S)
E = Maxwell3D.planewave(direction=ẑ, polarization=x̂, wavenumber=κ)
e = (n×E)×n

in_interior1 = CompScienceMeshes.interior_tpredicate(Γ1)
in_interior2 = CompScienceMeshes.interior_tpredicate(Γ2)
in_interior3 = CompScienceMeshes.interior_tpredicate(Γ3)

on_junction = CompScienceMeshes.overlap_gpredicate(γ)
edges1 = skeleton(e -> in_interior1(e) || on_junction(chart(Γ1,e)), Γ1,1)
edges2 = skeleton(e -> in_interior2(e) || on_junction(chart(Γ2,e)), Γ2,1)
edges3 = skeleton(e -> in_interior3(e) || on_junction(chart(Γ3,e)), Γ3,1)

X1 = raviartthomas(Γ1, edges1)
X2 = raviartthomas(Γ2, edges2)
X3 = raviartthomas(Γ3, edges3)
# X1, X2, X3 = raviartthomas(Γ1, γ), raviartthomas(Γ2, γ), raviartthomas(Γ3, γ)
X = X1 × X2 × X3

@hilbertspace j
@hilbertspace k
α, β = 1/(im*κ), log(abs(h))/(im*κ)
Eq = @varform T[k,j] + α*S[trc(k), dvg(j)] + α*St[dvg(k), trc(j)] + β*I[trc(k), trc(j)] == e[k]
eq = @discretise Eq  j∈X k∈X
u = solve(eq)

Θ, Φ = range(0.0,stop=π,length=100), 0.0
ffpoints = [point(cos(ϕ)*sin(θ), sin(ϕ)*sin(θ), cos(θ)) for θ in Θ for ϕ in Φ]
farfield = potential(MWFarField3D(κ*im), ffpoints, u, X)

using Plots
using LinearAlgebra
plot(Θ,norm.(farfield))


text:
module day17

data = Array{Bool}(hcat(map(line -> map(c -> c == '#', collect(line)), readlines("day17_input.txt"))...)')
steps = 6

# Part 1
function simulate(current)
    next = copy(current)
    region = CartesianIndices(current)
    ifirst, ilast = first(region), last(region)
    iunit = oneunit(ifirst)
    for i in region
        neighbours = [current[j] for j in max(ifirst, i - iunit):min(ilast, i + iunit) if j != i]
        if (current[i] && count(neighbours) ∉ [2, 3]) next[i] = false end
        if (!current[i] && count(neighbours) == 3) next[i] = true end
    end
    next
end

function runall(dim)
    grid = zeros(Bool, map(+, (size(data)..., fill(1, dim-2)...), tuple(fill(steps*2, dim)...)))
    s = steps + 1
    e = s + first(size(data)) - 1
    grid[s:e, s:e, fill(s:s, dim-2)...] = data
    for i in 1:steps
        grid = simulate(grid)
    end
    grid
end
println("Part 1: active count = $(sum(runall(3)))")

# Part 2
println("Part 2: active count = $(sum(runall(4)))")

end

max_stars_repo_name: ericgreveson/adventofcode

max_issues_repo_name: ericgreveson/adventofcode

max_forks_repo_name: ericgreveson/adventofcode


text:
adj = [0 1 0 1;
       1 0 1 0;
       0 1 0 1;
       1 0 1 0]
deg = [2 0 0 0;
       0 2 0 0;
       0 0 2 0;
       0 0 0 2]
lap = [2 -1 0 -1;
       -1 2 -1 0;
       0 -1 2 -1;
       -1 0 -1 2]
norm_lap = [1. -.5 0. -.5;
           -.5 1. -.5 0.;
           0. -.5 1. -.5;
           -.5 0. -.5 1.]
scaled_lap =   [0 -0.5 0 -0.5;
                -0.5 0 -0.5 -0;
                0 -0.5 0 -0.5;
                -0.5 0 -0.5 0]

@testset "linalg" begin
    fg = FeaturedGraph(adj)
    @test fg.matrix_type == :adjm
    for T in [Int8, Int16, Int32, Int64, Int128]
        @test GraphSignals.adjacency_matrix(adj, T) == T.(adj)
        @test GraphSignals.adjacency_matrix(fg, T) == T.(adj)
        @test GraphSignals.degrees(fg; dir=:both) == [2, 2, 2, 2]
        dm = GraphSignals.degree_matrix(fg, T; dir=:out)
        @test dm == T.(deg)
        @test GraphSignals.degree_matrix(adj, T; dir=:in) == dm
        @test GraphSignals.degree_matrix(adj, T; dir=:both) == dm
        @test GraphSignals.laplacian_matrix(fg, T) == T.(lap)
    end

    fg = FeaturedGraph(Float64.(adj))
    @test fg.matrix_type == :adjm
    for T in [Float16, Float32, Float64]
        dm = GraphSignals.degree_matrix(fg, T; dir=:out)
        @test dm == T.(deg)
        @test GraphSignals.degree_matrix(adj, T; dir=:in) == dm
        @test GraphSignals.degree_matrix(adj, T; dir=:both) == dm
        @test GraphSignals.laplacian_matrix(fg, T) == T.(lap)
        @test GraphSignals.normalized_laplacian(fg, T) ≈ T.(norm_lap)
        @test GraphSignals.scaled_laplacian(fg, T) ≈ T.(scaled_lap)

        fg_ = GraphSignals.laplacian_matrix!(deepcopy(fg), T)
        @test fg_.graph == T.(lap)
        @test fg_.matrix_type == :laplacian
        fg_ = GraphSignals.normalized_laplacian!(deepcopy(fg), T)
        @test fg_.graph ≈ T.(norm_lap)
        @test fg_.matrix_type == :normalized
        fg_ = GraphSignals.scaled_laplacian!(deepcopy(fg), T)
        @test fg_.graph ≈ T.(scaled_lap)
        @test fg_.matrix_type == :scaled
    end
end


text:
export resample

"""
	resample(A::AbstractGeoArray, resolution::Number;
			 crs::GeoFormat=crs(A), method::String="near")
	resample(A::AbstractGeoArray, snap::AbstractGeoArray; method::String="near")

`resample` uses `ArchGDAL.gdalwarp` to resample an `AbstractGeoArray`.

## Arguments
- `A`: The `AbstractGeoArray` to resample.
- `resolution`: A `Number` specifying the resolution for the output.
  If the keyword argument `crs` (described below) is specified, `resolution` must be in units of the `crs`.
- `snap`: an `AbstractGeoArray` whos resolution, crs and bounds will be snapped to.
  For best results it should roughly cover the same extent, or a subset of `A`.

## Keyword Arguments
- `crs`: A `GeoFormatTypes.GeoFormat` specifying an output crs
  (`A` with be reprojected to `crs` in addition to being resampled). Defaults to `crs(A)`
- `method`: A `String` specifying the method to use for resampling. Defaults to `"near"`
  (nearest neighbor resampling). See [resampling method](https://gdal.org/programs/gdalwarp.html#cmdoption-gdalwarp-r)
  in the gdalwarp docs for a complete list of possible values.

"""
function resample end
function resample(A::AbstractGeoArray, resolution::Number;
				  crs::GeoFormat=crs(A),
                  method::String="near")
    wkt = convert(String, convert(WellKnownText, crs))
    flags = ["-t_srs", "$(wkt)",
             "-tr", "$(resolution)", "$(resolution)",
             "-r", "$(method)"]
    AG.Dataset(A) do dataset
        AG.gdalwarp([dataset], flags) do warped
            GeoArray(warped)
        end
    end
end
function resample(A::AbstractGeoArray, snap::AbstractGeoArray; method::String="near")
    wkt = convert(String, convert(WellKnownText, crs(snap)))
    latres, lonres = map(abs ∘ step, span(snap, (Lat(), Lon())))
    (latmin, latmax), (lonmin, lonmax) = bounds(snap, (Lat(), Lon()))
    flags = ["-t_srs", "$(wkt)",
             "-tr", "$(latres)", "$(lonres)",
             "-te", "$(lonmin)", "$(latmin)", "$(lonmax)", "$(latmax)",
             "-r", "$(method)"]
    AG.Dataset(A) do dataset
        AG.gdalwarp([dataset], flags) do warped
            GeoArray(warped)
        end
    end
end