text:
module TontiTest
  using Decapods.TontiDiagrams
  using CombinatorialSpaces
  using GeometryBasics
  using Random
  using Test
  td = TontiDiagram();
  add_variables!(td, (:C,0,true),  (:dC,0,true),  (:ΔC,1,true),
                 (:∑ϕ,2,false), (:ϕ,1,false))

  s = EmbeddedDeltaSet2D{Bool, Point{3,Float64}}()
  points = [(0,0,0),(0,0,1),(0,1,0),(0,1,1)]

  add_vertices!(s, 4, point=points)
  glue_sorted_triangle!(s, 1,2,3)
  glue_sorted_triangle!(s, 2,3,4)

  sp = Space(s)

  @test gen_form(s, x->1) == [1,1,1,1]

  add_derivatives!(td, sp, :C=>:ΔC, :ϕ =>:∑ϕ)
  add_time_dep!(td, :dC, :C)

  const k = 0.003

  add_transition!(td, [:ΔC], (ϕ,ΔC)->(ϕ .= sp.hodge[1,2] * (k .* ΔC)), [:ϕ])
  add_transition!(td, [:∑ϕ], (dC, ∑ϕ)->(dC .= sp.hodge[2,3] * ∑ϕ), [:dC]);

  data_range, sim= vectorfield(td, sp)

  @test data_range[:C] == (1,4)

  Random.seed!(42)
  u  = gen_form(s,x->rand())
  du = zeros(Float64,nv(s))

  tempu = copy(u)
  dt = 0.001
  for i in 10000
    sim(du, tempu, [],0)
    tempu .+= du
    # Check that mass is conserved
  end
  @test all(tempu .!= u)
  @test all(sum(sp.hodge[1,1]*tempu)-sum(sp.hodge[1,1]*u) < 1e-6)
end


text:
#
# Datatypes
#
#TODO: make sure the writing is not changed
const DT_FIXED_POINT = UInt8(0) #| (UInt8(3) << 4)
const DT_FLOATING_POINT = UInt8(1) #| (UInt8(3) << 4)
const DT_TIME = UInt8(2) #| (UInt8(3) << 4)
const DT_STRING = UInt8(3) #| (UInt8(3) << 4)
const DT_BITFIELD = UInt8(4) #| (UInt8(3) << 4)
const DT_OPAQUE = UInt8(5) #| (UInt8(3) << 4)
const DT_COMPOUND = UInt8(6) #| (UInt8(3) << 4)
const DT_REFERENCE = UInt8(7) #| (UInt8(3) << 4)
const DT_ENUMERATED = UInt8(8) #| (UInt8(3) << 4)
const DT_VARIABLE_LENGTH = UInt8(9) #| (UInt8(3) << 4)
const DT_ARRAY = UInt8(10) #| (UInt8(3) << 4)

const DATATYPES = Dict{UInt8, String}(
    0 => "DT_FIXED_POINT",
    1 => "DT_FLOATING_POINT",
    2 => "DT_TIME",
    3 => "DT_STRING",
    4 => "DT_BITFIELD",
    5 => "DT_OPAQUE",
    6 => "DT_COMPOUND",
    7 => "DT_REFERENCE",
    8 => "DT_ENUMERATED",
    9 => "DT_VARIABLE_LENGTH",
    10=> "DT_ARRAY")

# This is the description for:
#    Strings
#    Opaque datatypes
#    References
struct BasicDatatype <: H5Datatype
    class::UInt8
    bitfield1::UInt8
    bitfield2::UInt8
    bitfield3::UInt8
    size::UInt32
end
define_packed(BasicDatatype)
StringDatatype(::Type{String}, size::Integer) =
    BasicDatatype(DT_STRING, 0x11, 0x00, 0x00, size)
OpaqueDatatype(size::Integer) =
    BasicDatatype(DT_OPAQUE, 0x00, 0x00, 0x00, size) # XXX make sure ignoring the tag is OK
ReferenceDatatype() =
    BasicDatatype(DT_REFERENCE, 0x00, 0x00, 0x00, jlsizeof(RelOffset))

function Base.:(==)(dt1::BasicDatatype, dt2::BasicDatatype)
    ret = true
    ret &= (dt1.class << 4) == (dt2.class << 4)
    ret &= dt1.bitfield1 == dt2.bitfield1
    ret &= dt1.bitfield2 == dt2.bitfield2
    ret &= dt1.bitfield3 == dt2.bitfield3
    ret &= dt1.size == dt2.size
    ret
end
# Reads a datatype message and returns a (offset::RelOffset, class::UInt8)
# tuple. If the datatype is committed, the offset is the offset of the
# committed datatype and the class is typemax(UInt8). Otherwise, the
# offset is the offset of the datatype in the file, and the class is
# the corresponding datatype class.
function read_datatype_message(io::IO, f::JLDFile, committed)
    if committed
        # Shared datatype
        jlread(io, UInt8) == 3 || throw(UnsupportedVersionException())
        jlread(io, UInt8) == 2 || throw(UnsupportedFeatureException())
        (typemax(UInt8), Int64(fileoffset(f, jlread(io, RelOffset))))
    else
        # Datatype stored here # bottom four bits
        class = jlread(io, UInt8) << 4 >> 4
        #println("datatype class $class")
        (class, Int64(position(io)-1))
    end
end

# Replace a symbol or expression in an AST with a new one
function replace_expr(x, from, to)
    ex = copy(x)
    for i = 1:length(ex.args)
        x = ex.args[i]
        if x == from
            ex.args[i] = to
        elseif isa(x, Expr)
            ex.args[i] = replace_expr(x, from, to)
        end
    end
    ex
end

# Macro to read an entire datatype from a file, avoiding dynamic
# dispatch for all but variable length types
macro read_datatype(io, datatype_class, datatype, then)
    esc(quote
        if $datatype_class == DT_FIXED_POINT
            $(replace_expr(then, datatype, :(jlread($io, FixedPointDatatype))))
        elseif $datatype_class == DT_FLOATING_POINT
            $(replace_expr(then, datatype, :(jlread($io, FloatingPointDatatype))))
        elseif $datatype_class == DT_STRING || $datatype_class == DT_OPAQUE || $datatype_class == DT_REFERENCE
            $(replace_expr(then, datatype, :(jlread($io, BasicDatatype))))
        elseif $datatype_class == DT_COMPOUND
            $(replace_expr(then, datatype, :(jlread($io, CompoundDatatype))))
        elseif $datatype_class == DT_VARIABLE_LENGTH
            $(replace_expr(then, datatype, :(jlread($io, VariableLengthDatatype))))
        elseif $datatype_class == DT_BITFIELD
            $(replace_expr(then, datatype, :(jlread($io, BitFieldDatatype))))
        else
            throw(UnsupportedFeatureException("invalid datatype class $datatype_class"))
        end
    end)
end

struct FixedPointDatatype <: H5Datatype
    class::UInt8
    bitfield1::UInt8
    bitfield2::UInt8
    bitfield3::UInt8
    size::UInt32
    bitoffset::UInt16
    bitprecision::UInt16
end
define_packed(FixedPointDatatype)
FixedPointDatatype(size::Integer, signed::Bool) =
    FixedPointDatatype(DT_FIXED_POINT, ifelse(signed, 0x08, 0x00), 0x00, 0x00, size, 0, 8*size)

struct BitFieldDatatype <: H5Datatype
    class::UInt8
    bitfield1::UInt8
    bitfield2::UInt8
    bitfield3::UInt8
    size::UInt32
    bitoffset::UInt16
    bitprecision::UInt16
end
define_packed(BitFieldDatatype)
BitFieldDatatype(size) =
    BitFieldDatatype(DT_BITFIELD, 0x00, 0x00, 0x00, size, 0, 8*size)

struct FloatingPointDatatype <: H5Datatype
    class::UInt8
    bitfield1::UInt8
    bitfield2::UInt8
    bitfield3::UInt8
    size::UInt32
    bitoffset::UInt16
    bitprecision::UInt16
    exponentlocation::UInt8
    exponentsize::UInt8
    mantissalocation::UInt8
    mantissasize::UInt8
    exponentbias::UInt32
end
define_packed(FloatingPointDatatype)

class(dt::Union{BasicDatatype,FixedPointDatatype,FloatingPointDatatype}) = dt.class

struct CompoundDatatype <: H5Datatype
    size::UInt32
    names::Vector{Symbol}
    offsets::Vector{Int}
    members::Vector{H5Datatype}

    function CompoundDatatype(size, names, offsets, members)
        length(names) == length(offsets) == length(members) ||
            throw(ArgumentError("names, offsets, and members must have same length"))
        new(size, names, offsets, members)
    end
end

Base.:(==)(x::CompoundDatatype, y::CompoundDatatype) =
    x.size == y.size && x.names == y.names && x.offsets == y.offsets &&
    x.members == y.members
Base.hash(::CompoundDatatype) = throw(ArgumentError("hash not defined for CompoundDatatype"))

class(dt::CompoundDatatype) = DT_COMPOUND
function jlsizeof(dt::CompoundDatatype)
    sz = jlsizeof(BasicDatatype) + size_size(dt.size)*length(dt.names)
    for i = 1:length(dt.names)
        # Extra byte for null padding of name
        sz += symbol_length(dt.names[i]) + 1 + jlsizeof(dt.members[i])
    end
    sz::Int
end

function jlwrite(io::IO, dt::CompoundDatatype)
    n = length(dt.names)
    jlwrite(io, BasicDatatype(DT_COMPOUND, n % UInt8, (n >> 8) % UInt8, 0x00, dt.size))
    for i = 1:length(dt.names)
        # Name
        name = dt.names[i]
        name_ptr = Base.unsafe_convert(Ptr{UInt8}, name)
        unsafe_write(io, name_ptr, symbol_length(name)+1)

        # Byte offset of member
        if dt.size <= typemax(UInt8)
            jlwrite(io, UInt8(dt.offsets[i]))
        elseif dt.size <= typemax(UInt16)
            jlwrite(io, UInt16(dt.offsets[i]))
        else
            jlwrite(io, UInt32(dt.offsets[i]))
        end

        # Member type message
        jlwrite(io, dt.members[i])
    end
end

function jlread(io::IO, ::Type{CompoundDatatype})
    dt = jlread(io, BasicDatatype)
    nfields = UInt16(dt.bitfield1) | UInt16(dt.bitfield2 << 8)
    dt.bitfield3 == 0 || throw(UnsupportedFeatureException())

    names = Vector{Symbol}(undef, nfields)
    offsets = Vector{Int}(undef, nfields)
    members = Vector{H5Datatype}(undef, nfields)
    for i = 1:nfields
        # Name
        names[i] = Symbol(read_bytestring(io))

        # Byte offset of member
        if dt.size <= typemax(UInt8)
            offsets[i] = jlread(io, UInt8)
        elseif dt.size <= typemax(UInt16)
            offsets[i] = jlread(io, UInt16)
        else
            offsets[i] = jlread(io, UInt32)
        end

        # Member type message
        datatype_class = jlread(io, UInt8)
        skip(io, -1)
        @read_datatype io datatype_class member begin
            members[i] = member
        end
    end

    CompoundDatatype(dt.size, names, offsets, members)
end

struct VariableLengthDatatype{T<:H5Datatype} <: H5Datatype
    class::UInt8
    bitfield1::UInt8
    bitfield2::UInt8
    bitfield3::UInt8
    size::UInt32
    basetype::T
end
VariableLengthDatatype(basetype::H5Datatype) =
    VariableLengthDatatype{typeof(basetype)}(DT_VARIABLE_LENGTH, 0x00, 0x00, 0x00, 8+jlsizeof(RelOffset), basetype)
VariableLengthDatatype(class, bitfield1, bitfield2, bitfield3, size, basetype::H5Datatype) =
    VariableLengthDatatype{typeof(basetype)}(class, bitfield1, bitfield2, bitfield3, size, basetype)

Base.:(==)(x::VariableLengthDatatype, y::VariableLengthDatatype) =
    x.class == y.class && x.bitfield1 == y.bitfield1 &&
    x.bitfield2 == y.bitfield2 && x.size == y.size &&
    x.basetype == y.basetype
Base.hash(::VariableLengthDatatype) = throw(ArgumentError("hash not defined for CompoundDatatype"))

class(dt::VariableLengthDatatype) = dt.class
jlsizeof(dt::VariableLengthDatatype) =
    jlsizeof(BasicDatatype) + jlsizeof(dt.basetype)

function jlwrite(io::IO, dt::VariableLengthDatatype)
    jlwrite(io, BasicDatatype(DT_VARIABLE_LENGTH, dt.bitfield1, dt.bitfield2, dt.bitfield3, dt.size))
    jlwrite(io, dt.basetype)
end

function jlread(io::IO, ::Type{VariableLengthDatatype})
    dtype = jlread(io, BasicDatatype)
    datatype_class = jlread(io, UInt8)
    skip(io, -1)
    @read_datatype io datatype_class dt begin
        VariableLengthDatatype(dtype.class, dtype.bitfield1, dtype.bitfield2, dtype.bitfield3, dtype.size, dt)
    end
end

jlsizeof(dt::CommittedDatatype) = 2 + jlsizeof(RelOffset)

function jlwrite(io::IO, dt::CommittedDatatype)
    jlwrite(io, UInt8(3))
    jlwrite(io, UInt8(2))
    jlwrite(io, dt.header_offset)
end

function commit(f::JLDFile,
        @nospecialize(dt::H5Datatype),
        attrs::Tuple{Vararg{WrittenAttribute}}=())
    psz = jlsizeof(HeaderMessage) * (length(attrs) + 1) + jlsizeof(dt)
    for attr in attrs
        psz += jlsizeof(attr)
    end
    io = f.io

    sz = jlsizeof(ObjectStart) + size_size(psz) + psz
    offset = f.end_of_data
    seek(io, offset)
    f.end_of_data = offset + sz + 4

    cio = begin_checksum_write(io, sz)
    jlwrite(cio, ObjectStart(size_flag(psz)))
    write_size(cio, psz)
    jlwrite(cio, HeaderMessage(HM_DATATYPE, jlsizeof(dt), 64))
    jlwrite(cio, dt)
    for attr in attrs
        jlwrite(cio, HeaderMessage(HM_ATTRIBUTE, jlsizeof(attr), 0))
        write_attribute(cio, f, attr, f.datatype_wsession)
    end
    jlwrite(io, end_checksum(cio))
end

# Read the actual datatype for a committed datatype
function read_committed_datatype(f::JLDFile, cdt::CommittedDatatype)
    io = f.io
    seek(io, fileoffset(f, cdt.header_offset))
    cio = begin_checksum_read(io)
    sz = read_obj_start(cio)
    pmax = position(cio) + sz

    # Messages
    datatype_class::UInt8 = 0
    datatype_offset::Int = 0
    attrs = ReadAttribute[]
    while position(cio) < pmax
        msg = jlread(cio, HeaderMessage)
        endpos = position(cio) + msg.size
        if msg.msg_type == HM_DATATYPE
            # Datatype stored here
            datatype_offset = position(cio)
            datatype_class = jlread(cio, UInt8)
        elseif msg.msg_type == HM_ATTRIBUTE
            push!(attrs, read_attribute(cio, f))
        end
        seek(cio, endpos)
    end
    seek(cio, pmax)

    # Checksum
    end_checksum(cio) == jlread(io, UInt32) || throw(InvalidDataException())

    seek(io, datatype_offset)
    @read_datatype io datatype_class dt begin
        return (dt, attrs)
    end
end

struct FixedLengthString{T<:AbstractString}
    length::Int
end

ext: jl


text:
const _clblast_status_codes = Dict{Int, Symbol}(
    # Status codes in common with the OpenCL standard
        0=>:CLBlastSuccess                   , # CL_SUCCESS
       -3=>:CLBlastOpenCLCompilerNotAvailable, # CL_COMPILER_NOT_AVAILABLE
       -4=>:CLBlastTempBufferAllocFailure    , # CL_MEM_OBJECT_ALLOCATION_FAILURE
       -5=>:CLBlastOpenCLOutOfResources      , # CL_OUT_OF_RESOURCES
       -6=>:CLBlastOpenCLOutOfHostMemory     , # CL_OUT_OF_HOST_MEMORY
      -11=>:CLBlastOpenCLBuildProgramFailure , # CL_BUILD_PROGRAM_FAILURE: OpenCL compilation error
      -30=>:CLBlastInvalidValue              , # CL_INVALID_VALUE
      -36=>:CLBlastInvalidCommandQueue       , # CL_INVALID_COMMAND_QUEUE
      -38=>:CLBlastInvalidMemObject          , # CL_INVALID_MEM_OBJECT
      -42=>:CLBlastInvalidBinary             , # CL_INVALID_BINARY
      -43=>:CLBlastInvalidBuildOptions       , # CL_INVALID_BUILD_OPTIONS
      -44=>:CLBlastInvalidProgram            , # CL_INVALID_PROGRAM
      -45=>:CLBlastInvalidProgramExecutable  , # CL_INVALID_PROGRAM_EXECUTABLE
      -46=>:CLBlastInvalidKernelName         , # CL_INVALID_KERNEL_NAME
      -47=>:CLBlastInvalidKernelDefinition   , # CL_INVALID_KERNEL_DEFINITION
      -48=>:CLBlastInvalidKernel             , # CL_INVALID_KERNEL
      -49=>:CLBlastInvalidArgIndex           , # CL_INVALID_ARG_INDEX
      -50=>:CLBlastInvalidArgValue           , # CL_INVALID_ARG_VALUE
      -51=>:CLBlastInvalidArgSize            , # CL_INVALID_ARG_SIZE
      -52=>:CLBlastInvalidKernelArgs         , # CL_INVALID_KERNEL_ARGS
      -53=>:CLBlastInvalidLocalNumDimensions , # CL_INVALID_WORK_DIMENSION: Too many thread dimensions
      -54=>:CLBlastInvalidLocalThreadsTotal  , # CL_INVALID_WORK_GROUP_SIZE: Too many threads in total
      -55=>:CLBlastInvalidLocalThreadsDim    , # CL_INVALID_WORK_ITEM_SIZE: ... or for a specific dimension
      -56=>:CLBlastInvalidGlobalOffset       , # CL_INVALID_GLOBAL_OFFSET
      -57=>:CLBlastInvalidEventWaitList      , # CL_INVALID_EVENT_WAIT_LIST
      -58=>:CLBlastInvalidEvent              , # CL_INVALID_EVENT
      -59=>:CLBlastInvalidOperation          , # CL_INVALID_OPERATION
      -61=>:CLBlastInvalidBufferSize         , # CL_INVALID_BUFFER_SIZE
      -63=>:CLBlastInvalidGlobalWorkSize     , # CL_INVALID_GLOBAL_WORK_SIZE
    #
    # Status codes in common with the clBLAS library
    -1024=>:CLBlastNotImplemented            , # Routine or functionality not implemented yet
    -1022=>:CLBlastInvalidMatrixA            , # Matrix A is not a valid OpenCL buffer
    -1021=>:CLBlastInvalidMatrixB            , # Matrix B is not a valid OpenCL buffer
    -1020=>:CLBlastInvalidMatrixC            , # Matrix C is not a valid OpenCL buffer
    -1019=>:CLBlastInvalidVectorX            , # Vector X is not a valid OpenCL buffer
    -1018=>:CLBlastInvalidVectorY            , # Vector Y is not a valid OpenCL buffer
    -1017=>:CLBlastInvalidDimension          , # Dimensions M, N, and K have to be larger than zero
    -1016=>:CLBlastInvalidLeadDimA           , # LD of A is smaller than the matrix's first dimension
    -1015=>:CLBlastInvalidLeadDimB           , # LD of B is smaller than the matrix's first dimension
    -1014=>:CLBlastInvalidLeadDimC           , # LD of C is smaller than the matrix's first dimension
    -1013=>:CLBlastInvalidIncrementX         , # Increment of vector X cannot be zero
    -1012=>:CLBlastInvalidIncrementY         , # Increment of vector Y cannot be zero
    -1011=>:CLBlastInsufficientMemoryA       , # Matrix A's OpenCL buffer is too small
    -1010=>:CLBlastInsufficientMemoryB       , # Matrix B's OpenCL buffer is too small
    -1009=>:CLBlastInsufficientMemoryC       , # Matrix C's OpenCL buffer is too small
    -1008=>:CLBlastInsufficientMemoryX       , # Vector X's OpenCL buffer is too small
    -1007=>:CLBlastInsufficientMemoryY       , # Vector Y's OpenCL buffer is too small
    #
    # Custom additional status codes for CLBlast
    -2050=>:CLBlastInsufficientMemoryTemp    , # Temporary buffer provided to GEMM routine is too small
    -2049=>:CLBlastInvalidBatchCount         , # The batch count needs to be positive
    -2048=>:CLBlastInvalidOverrideKernel     , # Trying to override parameters for an invalid kernel
    -2047=>:CLBlastMissingOverrideParameter  , # Missing override parameter(s) for the target kernel
    -2046=>:CLBlastInvalidLocalMemUsage      , # Not enough local memory available on this device
    -2045=>:CLBlastNoHalfPrecision           , # Half precision (16-bits) not supported by the device
    -2044=>:CLBlastNoDoublePrecision         , # Double precision (64-bits) not supported by the device
    -2043=>:CLBlastInvalidVectorScalar       , # The unit-sized vector is not a valid OpenCL buffer
    -2042=>:CLBlastInsufficientMemoryScalar  , # The unit-sized vector's OpenCL buffer is too small
    -2041=>:CLBlastDatabaseError             , # Entry for the device was not found in the database
    -2040=>:CLBlastUnknownError              , # A catch-all error code representing an unspecified error
    -2039=>:CLBlastUnexpectedError              # A catch-all error code representing an unexpected exception
)

struct CLBlastError <: Exception
    code::Int
    desc::Symbol

    function CLBlastError(c::Integer)
        new(c, get(_clblast_status_codes, Int(c), :CL_UNKNOWN_ERROR_CODE))
    end
end

Base.show(io::IO, err::CLBlastError) = print(io, "CLBlastError(code=$(err.code), $(err.desc))")


text:
# multi-bytes numeric types:
const mbntypes = Union{Int16, UInt16, Int32, UInt32, Int64, UInt64, Int128, UInt128, Float16, Float32, Float64}

"""
    LittleEndian{T, TSubCon<:Construct{T}} <: Adapter{T, T}

Little endian data adapter for serializing and deserializing.
"""
struct LittleEndian{T<:mbntypes, TSubCon<:Construct{T}} <: Adapter{T, T}
    subcon::TSubCon
end

LittleEndian(::Type{T}) where {T<:mbntypes} = LittleEndian(Construct(T))

subcon(wrapper::LittleEndian) = wrapper.subcon
encode(::LittleEndian{T, TSubCon}, obj::T; contextkw...) where {T, TSubCon} = htol(obj)
decode(::LittleEndian{T, TSubCon}, obj::T; contextkw...) where {T, TSubCon} = ltoh(obj)

"""
    BigEndian{T, TSubCon<:Construct{T}} <: Adapter{T, T}

Big endian data adapter for serializing and deserializing.
"""
struct BigEndian{T<:mbntypes, TSubCon<:Construct{T}} <: Adapter{T, T}
    subcon::TSubCon
end

BigEndian(::Type{T}) where {T<:mbntypes} = BigEndian(Construct(T))

subcon(wrapper::BigEndian) = wrapper.subcon
encode(::BigEndian{T, TSubCon}, obj::T; contextkw...) where {T, TSubCon} = hton(obj)
decode(::BigEndian{T, TSubCon}, obj::T; contextkw...) where {T, TSubCon} = ntoh(obj)


text:
push!(LOAD_PATH, joinpath(homedir(),"GitHub/KenyaCoV/src"))
using Plots,Parameters,Distributions,DifferentialEquations,JLD2,DataFrames,StatsPlots,FileIO,DelimitedFiles
using Statistics: median, quantile
using LinearAlgebra: eigen

#Load data
@load joinpath(homedir(),"Github/KenyaCoVOutputs/sims_consensus_baseline.jld2") sims_baseline
@load joinpath(homedir(),"Github/KenyaCoVOutputs/sims_consensus_baseline_scaled.jld2") sims_baseline_scaled
@load joinpath(homedir(),"Github/KenyaCoVOutputs/sims_consensus_control.jld2") sims_controls
@load joinpath(homedir(),"Github/KenyaCoVOutputs/sims_consensus_control_scaled.jld2") sims_controls_scaled

nametag = "controls"

## Functions --- to be used

include("hospitalisations.jl");

## Cumulative incidence
cum_incidence_A,cum_incidence_M,cum_incidence_V,cum_incidence_H = cum_incidence_for_each_sim_by_type(sims_controls)
cum_incidence_total = cum_incidence_A.+cum_incidence_M.+cum_incidence_V
total_incidence_peaks = find_peak_time_by_sim(sims_controls)
peak_days = [x[2] for x in total_incidence_peaks]

#Cumulative values for spreadsheet days 30,45,60,90,180,365(end)
#Total infections
median_total_infections = [ median(cum_incidence_total[:,t+1]) for t in [30,45,60,90,180,360,365] ]
lb_total_infections = [ quantile(cum_incidence_total[:,t+1],0.025) for t in [30,45,60,90,180,360,365] ]
ub_total_infections = [ quantile(cum_incidence_total[:,t+1],0.975) for t in [30,45,60,90,180,360,365] ]
total_inf_array = hcat(median_total_infections,lb_total_infections,ub_total_infections)
writedlm("output/save_total_inf_"*nametag*".csv",total_inf_array,',')
median_total_infections_at_peak = median([ cum_incidence_total[k,t+1] for (k,t) in enumerate(peak_days) ])
lb_total_infections_at_peak = quantile([ cum_incidence_total[k,t+1] for (k,t) in enumerate(peak_days) ],0.025)
ub_total_infections_at_peak = quantile([ cum_incidence_total[k,t+1] for (k,t) in enumerate(peak_days) ],0.975)
writedlm("output/save_total_inf_atpeak"*nametag*".csv",[median_total_infections_at_peak,lb_total_infections_at_peak,ub_total_infections_at_peak],',')

#Mild infection
median_mild_infections = [ median(cum_incidence_M[:,t+1]) for t in [30,45,60,90,180,360,365] ]
lb_mild_infections = [ quantile(cum_incidence_M[:,t+1],0.025) for t in [30,45,60,90,180,360,365] ]
ub_mild_infections = [ quantile(cum_incidence_M[:,t+1],0.975) for t in [30,45,60,90,180,360,365] ]
total_mild_array = hcat(median_mild_infections,lb_mild_infections,ub_mild_infections)
writedlm("output/save_total_mild_inf_"*nametag*".csv",total_mild_array,',')

median_mild_infections_at_peak = median([ cum_incidence_M[k,t+1] for (k,t) in enumerate(peak_days) ])
lb_mild_infections_at_peak = quantile([ cum_incidence_M[k,t+1] for (k,t) in enumerate(peak_days) ],0.025)
ub_mild_infections_at_peak = quantile([ cum_incidence_M[k,t+1] for (k,t) in enumerate(peak_days) ],0.975)
writedlm("output/save_total_mild_inf_atpeak"*nametag*".csv",[median_mild_infections_at_peak,lb_mild_infections_at_peak,ub_mild_infections_at_peak],',')

#Incidence
incidence_A = diff(cum_incidence_A,dims = 2)
incidence_M = diff(cum_incidence_M,dims = 2)
incidence_V = diff(cum_incidence_V,dims = 2)
incidence_H = diff(cum_incidence_H,dims = 2)
total_incidence = incidence_A.+incidence_M.+incidence_V

median_total_incidence = [ median(total_incidence[:,t]) for t in [30,45,60,90,180,360,365] ]
lb_total_incidence = [ quantile(total_incidence[:,t],0.025) for t in [30,45,60,90,180,360,365] ]
ub_total_incidence = [ quantile(total_incidence[:,t],0.975) for t in [30,45,60,90,180,360,365] ]
incidence_array = hcat(median_total_incidence,lb_total_incidence,ub_total_incidence)
writedlm("output/save_incidence_inf_"*nametag*".csv",incidence_array,',')
median_total_incidence_at_peak = median([ total_incidence[k,t] for (k,t) in enumerate(peak_days) ])
lb_total_incidence_at_peak = quantile([ total_incidence[k,t] for (k,t) in enumerate(peak_days) ],0.025)
ub_total_incidence_at_peak = quantile([ total_incidence[k,t] for (k,t) in enumerate(peak_days) ],0.975)
writedlm("output/save_incidence_inf_atpeak"*nametag*".csv",[median_total_incidence_at_peak,lb_total_incidence_at_peak,ub_total_incidence_at_peak],',')

#Mild infection
median_mild_incidence = [ median(incidence_M[:,t]) for t in [30,45,60,90,180,360,365] ]
lb_mild_incidence = [ quantile(incidence_M[:,t],0.025) for t in [30,45,60,90,180,360,365] ]
ub_mild_incidence = [ quantile(incidence_M[:,t],0.975) for t in [30,45,60,90,180,360,365] ]
mild_incidence_array = hcat(median_mild_incidence,lb_mild_incidence,ub_mild_incidence)
writedlm("output/save_mild_incidence_inf_"*nametag*".csv",mild_incidence_array,',')
median_mild_incidence_at_peak = median([ incidence_M[k,t] for (k,t) in enumerate(peak_days) ])
lb_mild_incidence_at_peak = quantile([ incidence_M[k,t] for (k,t) in enumerate(peak_days) ],0.025)
ub_mild_incidence_at_peak = quantile([ incidence_M[k,t] for (k,t) in enumerate(peak_days) ],0.975)
writedlm("output/save_mild_incidence_inf_atpeak"*nametag*".csv",[median_mild_incidence_at_peak,lb_mild_incidence_at_peak,ub_mild_incidence_at_peak],',')

#Severe infection
median_severe_incidence = [ median(incidence_H[:,t]) for t in [30,45,60,90,180,360,365] ]
lb_severe_incidence = [ quantile(incidence_H[:,t],0.025) for t in [30,45,60,90,180,360,365] ]
ub_severe_incidence = [ quantile(incidence_H[:,t],0.975) for t in [30,45,60,90,180,360,365] ]
median_severe_incidence_at_peak = median([ incidence_H[k,t] for (k,t) in enumerate(peak_days) ])
lb_severe_incidence_at_peak = quantile([ incidence_H[k,t] for (k,t) in enumerate(peak_days) ],0.025)
ub_severe_incidence_at_peak = quantile([ incidence_H[k,t] for (k,t) in enumerate(peak_days) ],0.975)

## Hospitalisation outcomes
hosp_occup_per_sim,ICU_occup_per_sim,new_ICU_per_sim,death_incidence_per_sim = total_hospital_outcomes_per_sim(sims_controls)
cum_death_per_sim = similar(death_incidence_per_sim)
for k = 1:1000
    cum_death_per_sim[k,:] .= cumsum(death_incidence_per_sim[k,:])
end
#Peak quantities

death_incidence_peaks = find_peaks(death_incidence_per_sim)
D_ts_at_peak = [median(death_incidence_peaks),quantile(death_incidence_peaks,0.025),quantile(death_incidence_peaks,0.975)]
writedlm("output/save_death_incidence_atpeak"*nametag*".csv",D_ts_at_peak,',')
ICU_incidence_peaks = find_peaks(new_ICU_per_sim)
ICU_incidence_at_peak = [median(ICU_incidence_peaks),quantile(ICU_incidence_peaks,0.025),quantile(ICU_incidence_peaks,0.975)]
writedlm("output/save_ICU_incidence_atpeak"*nametag*".csv",ICU_incidence_at_peak,',')

H_incidence_peaks = find_peaks(incidence_H)
H_incidence_at_peak = [median(H_incidence_peaks),quantile(H_incidence_peaks,0.025),quantile(H_incidence_peaks,0.975)]
writedlm("output/save_hosp_incidence_atpeak"*nametag*".csv",H_incidence_at_peak,',')

H_occup_peaks = find_peaks(hosp_occup_per_sim)
H_occup_at_peak = [median(H_occup_peaks),quantile(H_occup_peaks,0.025),quantile(H_occup_peaks,0.975)]
writedlm("output/save_hosp_occup_atpeak"*nametag*".csv",H_occup_at_peak,',')

ICU_occup_peaks = find_peaks(ICU_occup_per_sim)
ICU_occup_at_peak = [median(ICU_occup_peaks),quantile(ICU_occup_peaks,0.025),quantile(ICU_occup_peaks,0.975)]
writedlm("output/save_ICU_occup_atpeak"*nametag*".csv",ICU_occup_at_peak,',')

median_hosp_occup = [ median(hosp_occup_per_sim[:,t]) for t in [30,45,60,90,180,360,365] ]
lb_hosp_occup  = [ quantile(hosp_occup_per_sim[:,t],0.025) for t in [30,45,60,90,180,360,365] ]
ub_hosp_occup  = [ quantile(hosp_occup_per_sim[:,t],0.975) for t in [30,45,60,90,180,360,365] ]
hosp_occup_array = hcat(median_hosp_occup,lb_hosp_occup,ub_hosp_occup)
writedlm("output/save_hosp_occup_"*nametag*".csv",hosp_occup_array,',')

median_H = [ median(incidence_H[:,t]) for t in [30,45,60,90,180,360,365]]
lb_H = [ quantile(incidence_H[:,t],0.025) for t in [30,45,60,90,180,360,365] ]
ub_H = [ quantile(incidence_H[:,t],0.975) for t in [30,45,60,90,180,360,365] ]
hosp_incidence_array = hcat(median_H,lb_H,ub_H)
writedlm("output/save_hosp_incidence_"*nametag*".csv",hosp_incidence_array',',')

median_ICU_incidence = [ median(new_ICU_per_sim[:,t]) for t in [30,45,60,90,180,360,365] ]
lb_ICU_incidence  = [ quantile(new_ICU_per_sim[:,t],0.025) for t in [30,45,60,90,180,360,365] ]
ub_ICU_incidence  = [ quantile(new_ICU_per_sim[:,t],0.975) for t in [30,45,60,90,180,360,365] ]
ICU_incidence_array = vcat(median_ICU_incidence',lb_ICU_incidence',ub_ICU_incidence')
writedlm("output/save_ICU_incidence_"*nametag*".csv",ICU_incidence_array,',')
median_ICU_incidence_ts = [ median(new_ICU_per_sim[:,t]) for t in 1:365 ]

median_ICU_occup = [ median(ICU_occup_per_sim[:,t]) for t in [30,45,60,90,180,360,365] ]
lb_ICU_occup  = [ quantile(ICU_occup_per_sim[:,t],0.025) for t in [30,45,60,90,180,360,365] ]
ub_ICU_occup  = [ quantile(ICU_occup_per_sim[:,t],0.975) for t in [30,45,60,90,180,360,365] ]
ICU_occup_array = vcat(median_ICU_occup',lb_ICU_occup',ub_ICU_occup')
writedlm("output/save_ICU_occup_"*nametag*".csv",ICU_occup_array,',')

median_death_incidence = [ median(death_incidence_per_sim[:,t]) for t in [30,45,60,90,180,360,365] ]
lb_death_incidence  = [ quantile(death_incidence_per_sim[:,t],0.025) for t in [30,45,60,90,180,360,365] ]
ub_death_incidence  = [ quantile(death_incidence_per_sim[:,t],0.975) for t in [30,45,60,90,180,360,365] ]
death_incidence_array = vcat(median_death_incidence',lb_death_incidence',ub_death_incidence')
writedlm("output/save_death_incidence_"*nametag*".csv",death_incidence_array,',')

median_death_cum = [ median(cum_death_per_sim[:,t]) for t in [30,45,60,90,180,360,365] ]
lb_cum_death_incidence  = [ quantile(cum_death_per_sim[:,t],0.025) for t in [30,45,60,90,180,360,365] ]
ub_cum_death_incidence  = [ quantile(cum_death_per_sim[:,t],0.975) for t in [30,45,60,90,180,360,365] ]

cum_death_incidence_array = vcat(median_death_cum',lb_cum_death_incidence',ub_cum_death_incidence')
writedlm("output/save_death_cum_"*nametag*".csv",cum_death_incidence_array,',')

## Plotting
using Dates
gr()
first_of_months = [dayofyear(2020,m,1) for m = 4:12] .- dayofyear(2020,3,13)
median_deaths_ts = [ median(death_incidence_per_sim[:,t]) for t in 1:365 ]
lb_deaths_ts = [ quantile(death_incidence_per_sim[:,t],0.025) for t in 1:365 ]
ub_deaths_ts = [ quantile(death_incidence_per_sim[:,t],0.975) for t in 1:365 ]

median_H_ts = [ median(incidence_H[:,t]) for t in 1:365 ]
lb_H_ts = [ quantile(incidence_H[:,t],0.025) for t in 1:365 ]
ub_H_ts = [ quantile(incidence_H[:,t],0.975) for t in 1:365 ]

median_Hocc_ts = [ median(hosp_occup_per_sim[:,t]) for t in 1:365 ]
lb_HOcc_ts = [ quantile(hosp_occup_per_sim[:,t],0.025) for t in 1:365 ]
ub_HOcc_ts = [ quantile(hosp_occup_per_sim[:,t],0.975) for t in 1:365 ]

median_ICUocc_ts = [ median(ICU_occup_per_sim[:,t]) for t in 1:365 ]
lb_ICUOcc_ts = [ quantile(ICU_occup_per_sim[:,t],0.025) for t in 1:365 ]
ub_ICUOcc_ts = [ quantile(ICU_occup_per_sim[:,t],0.975) for t in 1:365 ]

# plt_incidence = plot(0:364,median_H_ts,
#                         lab = "Hospitalisations",
#                         lw = 3,color = :red,
#                         ribbon = (median_H_ts.-lb_H_ts,ub_H_ts .- median_H_ts),
#                         fillalpha = 0.25,
#                         xticks = (first_of_months,["Apr","May","June","July","Aug","Sept","Oct","Nov","Dec"]),
#                         xlims = (0.,275.),
#                         title = "Incidence of disease")
#
#  plot!(plt_incidence,0:364,median_deaths_ts,
#                         lab = "Deaths",
#                         lw = 3,color = :black,
#                         ribbon = (median_deaths_ts.-lb_deaths_ts,ub_deaths_ts .- median_deaths_ts),
#                         fillalpha = 0.5,
#                         xticks = (first_of_months,["Apr","May","June","July","Aug","Sept","Oct","Nov","Dec"]),
#                         xlims = (0.,275.),
#                         ylabel = "Daily Incidence" )

plot!(plt_incidence_both,0:364,median_H_ts,
                        lab = "",
                        lw = 3,color = :red,ls=:dot,
                        ribbon = (median_H_ts.-lb_H_ts,ub_H_ts .- median_H_ts),
                        fillalpha = 0.15,
                        xticks = (first_of_months,["Apr","May","June","July","Aug","Sept","Oct","Nov","Dec"]),
                        xlims = (0.,275.),
                        title = "Incidence of disease")
plot!(plt_incidence_both,0:364,median_deaths_ts,
                        lab = "",
                        lw = 3,color = :black,ls=:dot,
                        ribbon = (median_deaths_ts.-lb_deaths_ts,ub_deaths_ts .- median_deaths_ts),
                        fillalpha = 0.15,
                        xticks = (first_of_months,["Apr","May","June","July","Aug","Sept","Oct","Nov","Dec"]),
                        xlims = (0.,275.),
                        ylabel = "Daily Incidence" )

# plt_health_usage = plot(0:364,median_Hocc_ts,
#                         lab = "Hospital beds",
#                         lw = 3,color = :blue,
#                         ribbon = (median_Hocc_ts.-lb_HOcc_ts,ub_HOcc_ts .- median_Hocc_ts),
#                         fillalpha = 0.25,
#                         xticks = (first_of_months,["Apr","May","June","July","Aug","Sept","Oct","Nov","Dec"]),
#                         xlims = (0.,275.),
#                         ylabel = "Daily Occupancy",
#                         title = "Health system usage")
# plot!(plt_health_usage,0:364,median_ICUocc_ts,
#                         lab = "ICU beds",
#                         lw = 3,color = :green,
#                         ribbon = (median_ICUocc_ts.-lb_ICUOcc_ts,ub_ICUOcc_ts .- median_ICUocc_ts),
#                         fillalpha = 0.25,
#                         xticks = (first_of_months,["Apr","May","June","July","Aug","Sept","Oct","Nov","Dec"]),
#                         xlims = (0.,275.),
#                         ylabel = "Daily Occupancy",
#                         title = "Health system usage")

plot!(plt_health_usage_both,0:364,median_Hocc_ts,
                        lab = "",
                        lw = 3,color = :blue,ls=:dot,
                        ribbon = (median_Hocc_ts.-lb_HOcc_ts,ub_HOcc_ts .- median_Hocc_ts),
                        fillalpha = 0.25,
                        xticks = (first_of_months,["Apr","May","June","July","Aug","Sept","Oct","Nov","Dec"]),
                        xlims = (0.,275.),
                        ylabel = "Daily Occupancy",
                        title = "Health system usage")
plot!(plt_health_usage_both,0:364,median_ICUocc_ts,
                        lab = "",ls = :dot,
                        lw = 3,color = :green,
                        ribbon = (median_ICUocc_ts.-lb_ICUOcc_ts,ub_ICUOcc_ts .- median_ICUocc_ts),
                        fillalpha = 0.25,
                        xticks = (first_of_months,["Apr","May","June","July","Aug","Sept","Oct","Nov","Dec"]),
                        xlims = (0.,275.),
                        ylabel = "Daily Occupancy",
                        title = "Health system usage")

plot(0:364,[ median(cum_incidence_total[:,t+1]) for t in 1:365],lab = "Cum. incidence",
                        lw = 3,color = :blue,
                        fillalpha = 0.25,
                        xticks = (first_of_months,["Apr","May","June","July","Aug","Sept","Oct","Nov","Dec"]),
                        xlims = (0.,275.),
                        ylabel = "Total infections",title = "Controlled epidemic")
savefig("output/controlled_epidemic.png")

median([total_incidence[:,t] for t = 1:365])
savefig(plt_incidence,"output/plt_incidence.png")
savefig(plt_incidence_both,"output/plt_incidence_with_controls.png")

savefig(plt_health_usage,"output/plt_health_usage.png")
savefig(plt_health_usage_both,"output/plt_health_usage_with_controls.png")

# @save "output/plt_incidence.jld2" plt_incidence
# @save "output/plt_health_usage.jld2" plt_health_usage

@load "output/plt_incidence.jld2" plt_incidence
@load "output/plt_health_usage.jld2" plt_health_usage

plt_incidence_both = deepcopy(plt_incidence)
plt_health_usage_both = deepcopy(plt_health_usage)

function cum_incidence_for_each_sim_by_type(sims)
    n = length(sims.u)
    T = length(sims.u[1])
    cum_incidence_A = zeros(n,T)
    cum_incidence_M = zeros(n,T)
    cum_incidence_V = zeros(n,T)
    cum_incidence_H = zeros(n,T)

    for k = 1:n,t = 1:T
        cum_incidence_A[k,t] = sum(sims.u[k][t][:,:,1])
        cum_incidence_M[k,t] = sum(sims.u[k][t][:,:,2])
        cum_incidence_V[k,t] = sum(sims.u[k][t][:,:,3])
        cum_incidence_H[k,t] = sum(sims.u[k][t][:,:,4])
    end
    return cum_incidence_A,cum_incidence_M,cum_incidence_V,cum_incidence_H
end


text:
# parallel loops:
function buffon(n)
	hit = 0
	for i = 1:n
		mp = rand()
		phi = (rand() * pi) - pi / 2 # angle at which needle falls
		xright = mp + cos(phi)/2  # x-location of needle
		xleft = mp - cos(phi)/2
		# if xright >= 1 || xleft <= 0
		# 	hit += 1
		# end
		# Does needle cross either x == 0 or x == 1?
		p = (xright >= 1 || xleft <= 0) ? 1 : 0
		hit += p
	end
	miss = n - hit
	piapprox = n / hit * 2
end

@time buffon(100000)
#   0.208500 seconds (504.79 k allocations: 25.730 MiB, 7.10% gc time)
# 3.1441597233139444
@time buffon(100000000)
#   4.112683 seconds (5 allocations: 176 bytes)
# 3.141258861373451

function buffon_par(n)
	hit = @distributed (+) for i = 1:n
			mp = rand()
			phi = (rand() * pi) - pi / 2
			xright = mp + cos(phi)/2
			xleft = mp - cos(phi)/2
			(xright >= 1 || xleft <= 0) ? 1 : 0
		end
	miss = n - hit
	piapprox = n / hit * 2
end
@time buffon_par(100000)
#   1.058487 seconds (951.35 k allocations: 48.192 MiB, 2.04% gc time)
# 3.15059861373661

@time buffon_par(100000000)
#   0.735853 seconds (1.84 k allocations: 133.156 KiB)
# 3.1418106012575633

# arr is not correctly initialized:
arr = zeros(100000)
@distributed for i=1:100000
    arr[i] = i
end #> Task (queued) @0x00000000147ad430
#> arr only contains zeros!

# parallel maps:
using LinearAlgebra
function rank_marray()
	marr = [rand(1000,1000) for i=1:10]
	for arr in marr
	   println(LinearAlgebra.rank(arr))
	end
end

@time rank_marray()
# 1000
# 1000
# 1000
# 1000
# 1000
# 1000
# 1000
# 1000
# 1000
# 1000
# 7.310404 seconds (91.33 k allocations: 162.878 MiB, 1.15% gc time)

function prank_marray()
	marr = [rand(1000,1000) for i=1:10]
	println(pmap(LinearAlgebra.rank, marr))
end
@time prank_marray()
# {1000,1000,1000,1000,1000,1000,1000,1000,1000,1000}
# 5.966216 seconds (4.15 M allocations: 285.610 MiB, 2.15% gc time)


text:
using Plots, Statistics, MLDatasets, LinearAlgebra, Clustering, Printf, Dates, Random
path1=realpath(dirname(@__FILE__)*"/..")
include(string(path1,"/src/ANNBN.jl"))
using SpecialFunctions, Distributions
using PyCall, ScikitLearn
@sk_import ensemble: AdaBoostRegressor
using XGBoost, DecisionTree, LaTeXStrings
plot()
#


text:
# ==========================================================================
# This script compares the facilities offered by CovarianceEstimation.jl
# to existing packages and compares the performances for a set of reasonably
# large matrices.
#
# Libraries considered:
#
# - [python] sklearn
# - [R] corpcor
# - [Matlab/Octave] Ledoit & Wolf's implementations
#
# ==========================================================================

using CovarianceEstimation
using Random: seed!
using DelimitedFiles
using Statistics
using LinearAlgebra

seed!(65489) # random seed for reproducibility

p = [20, 40, 200, 400]
n = [40, 20, 400, 200]
L = [randn(pᵢ, pᵢ) for pᵢ ∈ p]
C = [Lᵢ' * Lᵢ for Lᵢ ∈ L]

X = [randn(n[i], p[i]) * L[i] for i ∈ eachindex(p)]

for Cᵢ ∈ C
    pᵢ = size(Cᵢ, 1)
    writedlm("benchmark/tmpmat/C_$(pᵢ).csv", Cᵢ)
end

for Xᵢ ∈ X
    nᵢ, pᵢ = size(Xᵢ)
    writedlm("benchmark/tmpmat/X_$(nᵢ)x$(pᵢ).csv", Xᵢ)
    # octave is unstable for eigenvalue decomposition so we do everything here
    F = eigen(Xᵢ' * Xᵢ / size(Xᵢ, 1))
    writedlm("benchmark/tmpmat/S_lambda_$(nᵢ)x$(pᵢ).csv", F.values)
    writedlm("benchmark/tmpmat/S_u_$(nᵢ)x$(pᵢ).csv", F.vectors)
end

# =========
# SKLEARN
# =========
oas = LinearShrinkage(
        target=DiagonalCommonVariance(),
        shrinkage=:oas)
lw = LinearShrinkage(
        target=DiagonalCommonVariance(),
        shrinkage=:lw)
times = zeros(length(p))
res = zeros(length(p))
for i ∈ eachindex(p)
    Xᵢ = X[i]
    Cᵢ = C[i]
    start = time()
    C_oas = cov(oas, Xᵢ)
    times[i] = time() - start
    res[i] = norm(C_oas - Cᵢ)
end
res_ours = res
times_ours = times
res_sklearn = [  51.05167435  214.21678806 1602.04944525 6533.94732887][:]
times_sklearn = [0.00026703 0.00023913 0.00225997 0.00262403][:]

δ = mean(res_ours - res_sklearn)
speedup = mean(times_sklearn ./ times_ours)

for i ∈ eachindex(p)
    Xᵢ = X[i]
    Cᵢ = C[i]
    start = time()
    C_lw = cov(lw, Xᵢ)
    times[i] = time() - start
    res[i] = norm(C_lw - Cᵢ)
end
res_ours = res
times_ours = times
res_sklearn = [  50.95012365  215.37414448 1602.01890586 6534.65450682][:]
times_sklearn = [0.00076485 0.00042701 0.00249386 0.00562191][:]

δ = mean(res_ours - res_sklearn)
speedup = mean(times_sklearn ./ times_ours)

# =========
# CORPCOR
# =========
ss = LinearShrinkage(
        target=DiagonalUnequalVariance(),
        shrinkage=:ss)
times = zeros(length(p))
res = zeros(length(p))
for i ∈ eachindex(p)
    Xᵢ = X[i]
    Cᵢ = C[i]
    start = time()
    C_ss = cov(ss, Xᵢ)
    times[i] = time() - start
    res[i] = norm(C_ss - Cᵢ)
end
res_ours = res
times_ours = times
res_corpcor = [  51.39624  219.00562 1606.20378 6562.85887][:]
times_corpcor = [0.001765013, 0.001574039, 0.1169169, 0.122371]

δ = mean(res_ours .- res_corpcor)

# roughly 22x speedup
speedup = mean(times_corpcor ./ times_ours)

# ===============
# LEDOIT WOLF 1
# ===============
lw = LinearShrinkage(
        target=ConstantCorrelation(),
        shrinkage=:lw)
times = zeros(length(p))
res = zeros(length(p))
for i ∈ eachindex(p)
    Xᵢ = X[i]
    Cᵢ = C[i]
    start = time()
    C_lw = cov(lw, Xᵢ)
    times[i] = time() - start
    res[i] = norm(C_lw - Cᵢ)
end
res_ours = res
times_ours = times
res_lw1 = [  51.503, 221.721, 1606.199, 6563.666]
times_lw1 = [0.0017309, 0.0015059, 0.0107820, 0.0208941]

δ = mean(res_ours .- res_lw1)

# roughly 10x speedup
speedup = mean(times_lw1 ./ times_ours)

# ===============
# LEDOIT WOLF 2
# ===============
lw = AnalyticalNonlinearShrinkage()
times = zeros(length(p))
res = zeros(length(p))
for i ∈ [1, 3] # fat matrix don't work reliably atm, see #38
    Xᵢ = X[i]
    Cᵢ = C[i]
    start = time()
    C_lw = cov(lw, Xᵢ)
    times[i] = time() - start
    res[i] = norm(C_lw - Cᵢ)
end
res_ours = res
times_ours = times
res_lw2 = [ 49.10317, 6985.10751]
times_lw2 = [0.015798, 0.007848]

# to our advantage
δ = mean(res_ours[[1, 3]] .- res_lw2)

# roughly 25x speedup
speedup = mean(times_lw2 ./ times_ours[[1, 3]])

# =========
# ASTROPY
# =========
bw = BiweightMidcovariance()
times = zeros(length(p))
res = zeros(length(p))
for i ∈ eachindex(p)
    Xᵢ = X[i]
    Cᵢ = C[i]
    start = time()
    C_bw = cov(bw, Xᵢ)
    times[i] = time() - start
    res[i] = norm(C_bw - Cᵢ)
end
res_ours = res
times_ours = times
res_astropy = [   70.41980995   331.84897257  1975.77842564 11556.0246496][:]
times_astropy = [0.00025606 0.00023079 0.00453758 0.00508738][:]

δ = mean(res_ours .- res_astropy)

# speedup varies from x5 for small matrices to ~similar speed for large matrices.
# The average speedup is ~x3 by this definition.
speedups = times_astropy ./ times_ours
speedup = mean(speedups)


text:
module Jot

# IMPORTS
using Base.Filesystem
using Dates
using JSON3
using OrderedCollections
using Parameters
using Pkg
using PrettyTables
using IsURL
using Random
using Setfield
using StructTypes
using TOML
import Base.delete!

# EXPORTS
export AWSConfig, LambdaException
export get_responder, AbstractResponder, LocalPackageResponder
export LocalImage, Container, RemoteImage, ECRRepo, AWSRole, LambdaFunction
export get_aws_role, get_user_labels
export LambdaFunctionState, pending, active
export get_dockerfile, build_definition
export run_image_locally, create_local_image, get_local_image
export send_local_request
export run_test
export stop_container, is_container_running
export create_ecr_repo, get_ecr_repo, push_to_ecr!
export get_remote_image
export create_aws_role
export create_lambda_function, get_lambda_function, invoke_function
export create_lambda_components, with_remote_image!, with_lambda_function!
export delete!
export show_lambdas
export get_labels, get_lambda_name
export get_all_local_images, get_all_remote_images, get_all_ecr_repos, get_all_lambda_functions
export get_all_containers, get_all_aws_roles
export LambdaComponents

# CONSTANTS
const docker_hash_limit = 12

abstract type LambdaComponent end

include("Responder.jl")
include("LocalImage.jl")
include("ECRRepo.jl")
include("RemoteImage.jl")
include("Container.jl")
include("LambdaFunction.jl")
include("AWS.jl")
include("LambdaComponents.jl")
include("Labels.jl")

function get_commit(path::String)::String
  try
    readchomp(`set -euo pipefail bash -c "PWD=pwd; cd $path; git rev-parse HEAD; cd \$PWD"`)
  catch e
    "None"
  end
end

function get_commit(mod::Module)::Union{Missing, String}
  mod_path = get_package_path(mod)
  get_commit(mod_path)
end

function get_package_path(mod::Module)::String
  module_path = Base.moduleroot(mod) |> pathof
  joinpath(splitpath(module_path)[begin:end-2]...)
end

function get_package_name(mod::Module)::String
  splitpath(get_package_path(mod))[end]
end

# -- get_tree_hash --

function get_tree_hash(res::LocalPackageResponder)::String
  get_tree_hash(joinpath(res.build_dir, res.package_name))
end

function get_tree_hash(path::String)::String
  Pkg.GitTools.tree_hash(path) |> bytes2hex
end

function get_tree_hash(i::Union{LocalImage, RemoteImage})::String
  get_labels(i).RESPONDER_TREE_HASH
end

function get_tree_hash(l::LambdaFunction)::String
  get_labels(l).RESPONDER_TREE_HASH
end

function get_tree_hash(lc::LambdaComponents)::String
  if !isnothing(lc.local_image)
    get_tree_hash(lc.local_image)
  elseif !isnothing(lc.remote_image)
    get_tree_hash(lc.remote_image)
  else
    get_tree_hash(lc.lambda_function)
  end
end

function get_image_full_name(
    aws_config::AWSConfig,
    image_suffix::String,
  )::String
  registry = "$(aws_config.account_id).dkr.ecr.$(aws_config.region).amazonaws.com"
  "$registry/$image_suffix"
end

function get_image_full_name_plus_tag(
    aws_config::AWSConfig,
    image_suffix::String,
    image_tag::String,
  )::String
  "$(get_image_full_name(aws_config, image_suffix)):$image_tag"
end

function get_image_full_name_plus_tag(image::LocalImage)::String
  "$(image.Repository):$(image.Tag)"
end

include("BuildDockerfile.jl")
include("Runtime.jl")
include("Scripts.jl")

function get_response_function_name(res::LocalPackageResponder)::String
  "$(res.package_name).$(res.response_function)"
end

function move_local_to_build_directory(build_dir::String, path::String, name::String)
  Base.Filesystem.cp(path, joinpath(build_dir, name))
end

function create_build_directory()::String
  build_dir = mktempdir()
  build_dir
end

function write_to_build_dir(
    content::String,
    build_dir::String,
    filename::String,
  )
  full_fpath = joinpath(build_dir, filename)
  open(full_fpath, "w") do f
    write(f, content)
  end
end

function add_scripts_to_build_dir(
    package_compile::Bool,
    julia_cpu_target::String,
    responder::AbstractResponder,
  )
  add_to_build(content, fname) = write_to_build_dir(content, responder.build_dir, fname)
  bootstrap_script |> x -> add_to_build(x, "bootstrap")
  get_init_script(package_compile, julia_cpu_target) |> x -> add_to_build(x, "init.jl")
  if package_compile
    get_precompile_jl(responder.package_name) |> x -> add_to_build(x, "precompile.jl")
  end
end

"""
    get_dockerfile(
        responder::AbstractResponder,
        julia_base_version::String,
        package_compile::Bool;
        user_defined_labels::AbstractDict{String, String} = AbstractDict{String, String}(),
        dockerfile_update::Function = x -> x,
      )::String

Returns contents for a Dockerfile. This function is called in `create_local_image` in order to
create a local docker image.
"""
function get_dockerfile(
    responder::AbstractResponder,
    julia_base_version::String,
    package_compile::Bool;
    user_defined_labels::AbstractDict{String, String} = AbstractDict{String, String}(),
    dockerfile_update::Function = x -> x,
  )::String
  overlapped_keys = [key for key in user_defined_labels if key in map(String, fieldnames(Labels))]
  if length(overlapped_keys) > 0
    error("User-defined labels $(join(overlapped_keys, ", ")) overlap with Jot-defined labels")
  end
  combined_labels = add_user_defined_labels(get_labels(responder), user_defined_labels)
  generated_dockerfile = foldl(
    *, [
    dockerfile_add_julia_image(julia_base_version),
    dockerfile_add_utilities(),
    dockerfile_add_runtime_directories(),
    dockerfile_add_additional_registries(responder.registry_urls),
    dockerfile_copy_build_dir(),
    dockerfile_add_responder(responder),
    dockerfile_add_labels(combined_labels),
    dockerfile_add_jot(),
    dockerfile_add_aws_rie(),
    dockerfile_add_bootstrap(responder.package_name,
                             String(responder.response_function),
                             responder.response_function_param_type),
    dockerfile_add_precompile(package_compile),
  ]; init = "")
  dockerfile_update(generated_dockerfile)
end

"""
    create_local_image(
        responder::AbstractResponder;
        image_suffix::Union{Nothing, String} = nothing,
        aws_config::Union{Nothing, AWSConfig} = nothing,
        image_tag::String = "latest",
        no_cache::Bool = false,
        julia_base_version::String = "1.6.1",
        julia_cpu_target::String = "x86-64",
        package_compile::Bool = false,
        user_defined_labels::AbstractDict{String, String} = OrderedDict{String, String}(),
        dockerfile_update::Function = x -> x,
        build_args::AbstractDict{String, String} = OrderedDict{String, String}(),
      )::LocalImage

Creates a locally-stored docker image containing the specified responder. This can be tested
tested locally, or directly uploaded to an AWS ECR Repo for use as an AWS Lambda function.

`package_compile` indicates whether `PackageCompiler.jl` should be used to compile the image -
this is not necessarily for testing/exploration but is highly recommended for production use.

Use `no_cache` to construct the local image without using a cache; this is sometimes necessary
if nothing locally has changed, but the image is querying a remote object (say, a github repo)
which has changed.

If `user_defined_labels` are defined, these will be added to the generated `LocalImage`, as well
as all subsequent types based on the `LocalImage`, like the remote image or the ultimate Lambda
function. They can be retrieved using the `get_labels` function, alongside the Jot-generated
labels.

`dockerfile_update` is a function that accepts the pre-generated Dockerfile as its only argument,
and returns a new Dockerfile. The most likely use case for this is to add customized extensions
to the generated Dockerfile. For example, passing
`(dockerfile) -> dockerfile * "RUN ssh-keygen -t rsa -f .ssh/id_rsa -N"`
will cause an SSH key to be created within the docker image.

`build_args` are arguments that are appended to the `docker build` command using the `--build-arg`
command-line argument.
"""
function create_local_image(
    responder::AbstractResponder;
    image_suffix::Union{Nothing, String} = nothing,
    aws_config::Union{Nothing, AWSConfig} = nothing,
    image_tag::String = "latest",
    no_cache::Bool = false,
    julia_base_version::String = "1.6.1",
    julia_cpu_target::String = "x86-64",
    package_compile::Bool = false,
    user_defined_labels::AbstractDict{String, String} = OrderedDict{String, String}(),
    dockerfile_update::Function = x -> x,
    build_args::AbstractDict{String, String} = OrderedDict{String, String}(),
  )::LocalImage
  # TODO check if the image suffix already exists
  image_suffix = isnothing(image_suffix) ? get_lambda_name(responder) : image_suffix
  aws_config = isnothing(aws_config) ? get_aws_config() : aws_config
  add_scripts_to_build_dir(package_compile, julia_cpu_target, responder)
  dockerfile = get_dockerfile(responder,
                              julia_base_version,
                              package_compile;
                              user_defined_labels,
                              dockerfile_update
                             )
  @debug dockerfile
  open(joinpath(responder.build_dir, "Dockerfile"), "w") do f
    write(f, dockerfile)
  end
  image_name_plus_tag = get_image_full_name_plus_tag(aws_config,
                                                     image_suffix,
                                                     image_tag)
  # Build the actual image
  build_cmd = get_dockerfile_build_cmd(dockerfile,
                                       image_name_plus_tag,
                                       no_cache,
                                       build_args)
  run(Cmd(build_cmd, dir=responder.build_dir))
  # Locate it and return it
  id_fname = joinpath(responder.build_dir, "id")
  @debug id_fname
  image_id = open(id_fname, "r") do f
    full_str = String(read(f))
    @debug full_str
    split(full_str, ':')[2]
  end
  @debug image_id
  this_image = get_local_image(image_id)
  isnothing(this_image) ? error("Unable to locate created local image with image ID $image_id") : this_image
end

function parse_docker_ls_output(::Type{T}, raw_output::AbstractString)::Vector{T} where {T}
  if raw_output == ""
    Vector{T}()
  else
    # as_lower = lowercase(raw_output)
    by_line = split(raw_output, "\n")
    as_dicts = [JSON3.read(js_str) for js_str in by_line]
    # type_args = [[dict[Symbol(fn)] for fn in fieldnames(T)] for dict in as_dicts]
    [T(; filter(p -> p.first in fieldnames(T), dict)..., exists=true) for dict in as_dicts]
  end
end

# -- get_labels --

function get_labels(res::LocalPackageResponder)::Labels
  Labels( RESPONDER_PACKAGE_NAME=res.package_name,
          RESPONDER_FUNCTION_NAME=get_responder_function_name(res),
          RESPONDER_COMMIT=get_commit(res),
          RESPONDER_TREE_HASH=get_tree_hash(res),
          RESPONDER_PKG_SOURCE=get_responder_path(res),
          IS_JOT_GENERATED="true",
         )
end

"""
    get_user_labels(l::Union{LocalImage, ECRRepo, RemoteImage, LambdaFunction})::Dict

Retrieves any user_defined labels for the given resource as a Dict of key=>value pairs.
"""
function get_user_labels(l::Union{LocalImage, ECRRepo, RemoteImage, LambdaFunction})::Dict
  get_labels(l) |> filter_user_defined_only
end

function get_labels(image::LocalImage)::Labels
  image_inspect_json = readchomp(`docker inspect $(image.ID)`)
  iis = JSON3.read(image_inspect_json, Vector{Dict{String, Any}})
  if length(iis) == 0
    error("Unable to find labels for image $(image.Repository)")
  else
    get_labels(iis[1])
  end
end

function get_labels(ecr_repo::ECRRepo)::Labels
  tags_json = readchomp(
    `aws ecr list-tags-for-resource
       --resource-arn $(ecr_repo.repositoryArn)`
  )
  tags_raw = JSON3.read(tags_json)
  haskey(tags_raw, "tags") || error("ECR Repo $(ecr_repo.repositoryName) has no tags")
  tags_dict = Dict(d["Key"] => d["Value"] for d in tags_raw["tags"])
  Labels(tags_dict)
end

function get_labels(image::RemoteImage)::Labels
  batch_image_json = readchomp(
    `aws ecr batch-get-image
      --repository-name $(image.ecr_repo.repositoryName)
      --image-id imageTag=$(image.imageTag)
      --accepted-media-types "application/vnd.docker.distribution.manifest.v1+json" --output json`
     )
  batch_image = JSON3.read(batch_image_json)
  try
    manifest = JSON3.read(batch_image["images"][1]["imageManifest"])
    v1_compat = JSON3.read(manifest["history"][1]["v1Compatibility"])
    labels = v1_compat["config"]["Labels"]
    labels = OrderedDict{Symbol, String}(labels)
    labels = OrderedDict(String(k) => v for (k, v) in labels)
    Labels(labels)
  catch e
    if isa(e, BoundsError) || isa(e, KeyError)
      error("Unable to find labels for image $(image.ecr_repo.repositoryName)")
    else
      throw(e)
    end
  end
end

function get_labels(image_inspect::AbstractDict{String, Any})::Labels
  ii = image_inspect["Config"]["Labels"]
  isnothing(ii) && error("Unable to get labels")
  ii_sym = Dict(k => v for (k, v) in ii)
  Labels(ii_sym)
end

function get_labels(lambda_function::LambdaFunction)::Labels
  get_lf_tags_script = get_lambda_function_tags_script(lambda_function)
  tags_json = readchomp(`bash -c $get_lf_tags_script`)
  try
    JSON3.read(tags_json, Dict{String, Labels})["Tags"]
  catch e
    error("Unable to find labels for lambda function $(lambda_function.FunctionName)")
  end
end

function get_labels(lc::LambdaComponents)::Labels
  if !isnothing(lc.local_image)
    get_labels(lc.local_image)
  elseif !isnothing(lc.remote_image)
    get_labels(lc.remote_image)
  else
    get_labels(lc.lambda_function)
  end
end

function is_jot_generated(c::Union{ECRRepo, LambdaComponent})::Bool
  try
    labels = get_labels(c)
    Meta.parse(labels.IS_JOT_GENERATED)
  catch e
    false
  end
end

"""
    run_test(
      image::LocalImage,
      function_argument::Any = "",
      expected_response::Any = nothing;
      then_stop::Bool = false,
    )::Tuple{Bool, Float64}

Runs a test of the given local docker image, passing `function_argument` (if given), and expecting
`expected_response`(if given). If a function_argument is not given, then it will merely test
that any kind of response is received - this response may be an error JSON and the test will still
pass, establishing only that the image can be contacted. Returns a tuple of (test result, time)
where time is the time taken for a response to be received, in seconds.

The test will use an already-running docker container, if one exists. If this is the case then the
`then_stop` parameter tells the function whether to stop the docker container after running the
test. If the run_test function finds no docker container already running, it will start one, and
then shut it down afterwards, regardless of the value of `then_stop`.
"""
function run_test(
    image::LocalImage,
    function_argument::Any = "",
    expected_response::Any = nothing;
    then_stop::Bool = false,
  )::Tuple{Bool, Float64}
  running = get_all_containers(image)
  con = length(running) == 0 ? run_image_locally(image) : nothing
  time_taken = @elapsed actual = send_local_request(function_argument)
  !isnothing(con) && (stop_container(con); delete!(con))
  then_stop && foreach(stop_container, running)
  if !isnothing(expected_response)
    passed = actual == expected_response
    passed && @info "Local test passed in $time_taken seconds; result received matched expected $actual"
    !passed && @info "Local test failed; actual: $actual was not equal to expected: $expected_response"
    (passed, time_taken)
  else
    @info "Response received from container; test passed"
    (true, time_taken)
  end
end

"""
    run_test(
      func::LambdaFunction,
      function_argument::Any = "",
      expected_response::Any = nothing;
    )::Tuple{Bool, Float64}

Runs a test of the given Lambda Function, passing `function_argument` (if given), and expecting
`expected_response`(if given). If a function_argument is not given, then it will merely test
that any kind of response is received - this response may be an error JSON and the test will still
pass, establishing only that the function can be contacted. Returns a tuple of (test result, time)
where time is the time taken for a response to be received, in seconds.
"""
function run_test(
    func::LambdaFunction,
    function_argument::Any = "",
    expected_response::Any = nothing;
    check_function_state::Bool = false,
  )::Tuple{Bool, Union{Missing, Float64}}
  try
    time_taken = @elapsed actual = invoke_function(function_argument, func; check_state = check_function_state)
    passed = actual == expected_response
    passed && @info "Remote test passed in $time_taken seconds; result received matched expected $actual"
    !passed && @info "Remote test failed; actual: $actual was not equal to expected: $expected_response"
    (passed, time_taken)
  catch e
    if isa(e, LambdaException) && isnothing(expected_response)
      @info "Error response received from lambda function; test passed"
      (true, missing)
    else
      rethrow()
    end
  end
end

function ecr_login_for_image(aws_config::AWSConfig, image_suffix::String)
  script = get_ecr_login_script(aws_config, image_suffix)
  run(`bash -c $script`)
end

function ecr_login_for_image(image::LocalImage)
  ecr_login_for_image(get_aws_config(image), get_lambda_name(image))
end

"""
    push_to_ecr!(image::LocalImage)::RemoteImage
Pushes the given local docker image to an AWS ECR Repo, a prerequisite of creating an AWS Lambda
Function. If an ECR Repo for the given local image does not exist, it will be created
automatically. Returns a RemoteImage object that represents the docker image that is hosted on the
ECR Repo. The ECR Repo itself is an attribute of the RemoteImage.
"""
function push_to_ecr!(image::LocalImage)::RemoteImage
  image.exists || error("Image does not exist")
  ecr_login_for_image(image)
  existing_repo = get_ecr_repo(image)
  repo = if isnothing(existing_repo)
    create_ecr_repo(image)
  else
    existing_repo
  end
  push_script = get_image_full_name_plus_tag(image) |> get_docker_push_script
  readchomp(`bash -c $push_script`)
  all_images = get_all_local_images()
  img_idx = findfirst(img -> img.ID[1:docker_hash_limit] == image.ID[1:docker_hash_limit], all_images)
  image.Digest = all_images[img_idx].Digest
  out = get_remote_image(image)
  @debug out
  out
end

"""
    create_lambda_function(
        remote_image::RemoteImage;
        role::Union{AWSRole, Nothing} = nothing,
        function_name::Union{Nothing, String} = nothing,
        timeout::Int64 = 60,
        memory_size::Int64 = 2000,
      )::LambdaFunction
Creates a function that exists on the AWS Lambda service. The function will use the given
`RemoteImage`, and runs using the given AWS Role.
"""
function create_lambda_function(
    remote_image::RemoteImage;
    role::Union{AWSRole, Nothing} = nothing,
    function_name::Union{Nothing, String} = nothing,
    timeout::Int64 = 60,
    memory_size::Int64 = 2000,
  )::LambdaFunction
  function_name = isnothing(function_name) ? remote_image.ecr_repo.repositoryName : function_name
  role = isnothing(role) ? create_aws_role(create_role_name(function_name)) : role
  image_uri = "$(remote_image.ecr_repo.repositoryUri):$(remote_image.imageTag)"
  @debug image_uri
  labels = get_labels(remote_image)
  out = create_lambda_function(image_uri, role, function_name, timeout, memory_size, labels)
  @debug out
  out
end

function create_role_name(function_name::String)::String
  function_name * "-generated-lambda-role"
end

"""
    create_lambda_function(
        repo::ECRRepo;
        role::AWSRole = nothing,
        function_name::Union{Nothing, String} = nothing,
        image_tag::String = "latest",
        timeout::Int64 = 60,
        memory_size::Int64 = 2000,
      )::LambdaFunction
Creates a function that exists on the AWS Lambda service. The function will use the given ECR Repo,
and runs using the given AWS Role. If given, the image_tag will decide which of the images in the
ECR Repo is used.
"""
function create_lambda_function(
    repo::ECRRepo;
    role::AWSRole = nothing,
    function_name::Union{Nothing, String} = nothing,
    image_tag::String = "latest",
    timeout::Int64 = 60,
    memory_size::Int64 = 2000,
  )::LambdaFunction
  role = isnothing(role) ? create_aws_role(get_lambda_name(repo) * "-lambda-role") : role
  function_name = isnothing(function_name) ? repo.repositoryName : function_name
  image_uri = "$(repo.repositoryUri):$image_tag"
  labels = get_labels(repo)
  create_lambda_function(image_uri, role, function_name, timeout, memory_size, labels)
end

function create_lambda_function(
    image_uri::String,
    role::AWSRole,
    function_name::String,
    timeout::Int64,
    memory_size::Int64,
    labels::Labels,
  )::LambdaFunction
  existing_lf = get_lambda_function(function_name)
  if !isnothing(existing_lf)
    @info "Lambda function $function_name already exists; overwriting"
    delete!(existing_lf)
  end
  aws_role_has_lambda_execution_permissions(role) || error("Role $role does not have permission to execute Lambda functions")
  create_script = get_create_lambda_function_script(function_name,
                                                    image_uri,
                                                    role.Arn,
                                                    timeout,
                                                    memory_size;
                                                    labels=labels,
                                                   )

  out = Pipe(); err = Pipe()
  proc = run(pipeline(ignorestatus(`bash -c $create_script`), stdout=out, stderr=err), wait=true)
  close(out.in); close(err.in)
  if proc.exitcode != 0
    @info read(err, String)
    @info read(out, String)
    error("proc exited with $(proc.exitcode)")
  end

  while true
    sleep(1)
    if get_function_state(function_name) == Active
      break
    end
  end

  func_json = read(out, String)
  return JSON3.read(func_json, LambdaFunction)
end

"""
    run_image_locally(local_image::LocalImage; detached::Bool=true)::Container

Runs the given local image, starting a docker container. If `detached`, the container will run in
the background. The container can be stopped/deleted by eg `stop_container`, `delete!`.
"""
function run_image_locally(local_image::LocalImage; detached::Bool=true)::Container
  args = ["-p", "9000:8080"]
  detached && push!(args, "-d")
  container_id = readchomp(`docker run $args $(local_image.ID)`)
  Container(ID=container_id, Image=local_image.ID)
end

"""
    send_local_request(request::Any)

Make a function call to a locally-running docker container and returns the value. A container can
be initiated by eg `run_image_locally`.
"""
function send_local_request(request::Any)
  endpoint = "http://localhost:9000/2015-03-31/functions/function/invocations"
  http = HTTP.post(
                   endpoint,
                   [],
                   JSON3.write(request),
                  )
  @debug http
  JSON3.read(http.body)
end

function get_environment_variables(image_inspect::AbstractDict{String, Any})::Dict{String, String}
  envs = image_inspect["ContainerConfig"]["Env"]
  Dict(
       env_split[1] => env_split[2] for env_split in map(x -> split(x, "="), envs)
      )
end

# -- get_lambda_name --

function get_lambda_name(res::AbstractResponder)::String
  lowercase(res.package_name)
end

function get_lambda_name(local_image::LocalImage)::Union{Nothing, String}
  @debug local_image.Repository
  is_jot_generated(local_image) ? split(local_image.Repository, '/')[2] : nothing
end

function get_lambda_name(remote_image::RemoteImage)::Union{Nothing, String}
  is_jot_generated(remote_image) ?  get_lambda_name(remote_image.ecr_repo) : nothing
end

function get_lambda_name(repo::ECRRepo)::Union{Nothing, String}
  is_jot_generated(repo) ? repo.repositoryName : nothing
end

function get_lambda_name(lf::LambdaFunction)::Union{Nothing, String}
  is_jot_generated(lf) ? lf.FunctionName : nothing
end

function get_lambda_name(l::LambdaComponents)::String
  get_from_any_component(get_lambda_name, l)
end

# -- get_response_function_name --

function get_response_function_name(
    env_vars::AbstractDict{String, String}
  )::Union{Nothing, String}
  try
    env_vars["FUNC_FULL_NAME"]
  catch e
    isa(e, KeyError) ? nothing : throw(e)
  end
end

function get_response_function_name(c::LambdaComponent)::String
  labels = get_labels(c)
  "$(labels.RESPONDER_PACKAGE_NAME).$(labels.RESPONDER_FUNCTION_NAME)"
end

function get_response_function_name(lambda::LambdaComponents)::String
  get_from_any_component(get_response_function_name, lambda)
end

end


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export
    connectedcomponents, connectedmembership,
    cutvertices, bridges, edgebiconnectedcomponents

"""
    connectedcomponents(graph::UndirectedGraph) -> Vector{Set{Int}}

Compute connectivity and return sets of the connected components.
"""
@cachefirst function connectedcomponents(graph::UndirectedGraph)
    nodes = nodeset(graph)
    components = Set{Int}[]
    while !isempty(nodes)
        root = pop!(nodes)
        tree = Set(keys(dfstree(adjacencies, graph, root)))
        push!(components, tree)
        setdiff!(nodes, tree)
    end
    return components
end

connectedcomponents(view::SubgraphView) = connectedcomponents(view.graph)

"""
    connectedmembership(graph::OrderedGraph) -> Vector{Int}

Return connected component membership array.
"""
@cachefirst function connectedmembership(graph::OrderedGraph)
    mem = zeros(Int, nodecount(graph))
    for (i, conn) in enumerate(connectedcomponents(graph))
        for n in conn
            mem[n] = i
        end
    end
    return mem
end

connectedmembership(view::SubgraphView) = connectedmembership(view.graph)

struct BiconnectedState{T<:UndirectedGraph}
    graph::T

    pred::Dict{Int,Int}
    level::Dict{Int,Int}
    low::Dict{Int,Int}

    cutvertices::Set{Int}
    bridges::Vector{Int}
    biconnected::Vector{Set{Int}} # biconnected edges

    function BiconnectedState{T}(graph) where {T<:UndirectedGraph}
        new(graph, Dict(), Dict(), Dict(), Set(), [], [])
    end
end

function dfs!(state::BiconnectedState, n::Int)
    state.pred[n] = n
    dfs!(state, 1, n)
end

function dfs!(state::BiconnectedState, depth::Int, n::Int)
    state.level[n] = depth
    state.low[n] = depth
    childcnt = 0
    compbuf = Set{Int}()
    for (ninc, nadj) in neighbors(state.graph, n)
        if !haskey(state.level, nadj)
            # New node
            childcnt += 1
            state.pred[nadj] = n
            comp = dfs!(state, depth + 1, nadj)
            push!(comp, ninc)
            if state.low[nadj] >= state.level[n]
                # Articulation point
                if state.low[nadj] > state.level[n]
                    push!(state.bridges, ninc) # except for bridgehead
                end
                push!(state.biconnected, comp)
                push!(state.cutvertices, n)
            else
                union!(compbuf, comp)
            end
            state.low[n] = min(state.low[n], state.low[nadj])
        elseif state.pred[n] != nadj
            # Cycle found
            state.low[n] = min(state.low[n], state.level[nadj])
            state.pred[nadj] = n
            push!(compbuf, ninc)
        end
    end
    if depth == 1 && childcnt < 2
        delete!(state.cutvertices, n)
    end
    return compbuf
end

function findbiconnected(graph::T, sym::Symbol) where {T<:UndirectedGraph}
    state = BiconnectedState{T}(graph)
    nodes = nodeset(graph)
    while !isempty(nodes)
        dfs!(state, pop!(nodes))
        setdiff!(nodes, keys(state.level))
    end
    return getproperty(state, sym)
end

"""
    cutvertices(graph::UndirectedGraph) -> Set{Int}

Compute biconnectivity and return cut vertices (articulation points).
"""
@cachefirst function cutvertices(graph::UndirectedGraph)
    return findbiconnected(graph, :cutvertices)
end

cutvertices(view::SubgraphView) = cutvertices(view.graph)

"""
    bridges(graph::UndirectedGraph) -> Set{Int}

Compute biconnectivity and return bridges.
"""
@cachefirst function bridges(graph::UndirectedGraph)
    return findbiconnected(graph, :bridges)
end

bridges(view::SubgraphView) = bridges(view.graph)

"""
    edgebiconnectedcomponents(graph::UndirectedGraph) -> Vector{Vector{Int}}

Compute sets of biconnected component edges.
"""
@cachefirst function edgebiconnectedcomponents(graph::UndirectedGraph)
    return findbiconnected(graph, :biconnected)
end

edgebiconnectedcomponents(view::SubgraphView
    ) = edgebiconnectedcomponents(view.graph)


text:
function gradient(f::Function,uh::FEFunction)
  fuh = f(uh)
  _gradient(f,uh,fuh)
end

function _gradient(f,uh,fuh::AbstractArray)
  @unreachable """\n
  In order to perform AD on a Function taking a FEFunction as argument, such Function
  has to return a DomainContribution.

  Make sure that you are using a Measure instead of a CellQuadrature to perform integration.
  """
end

function _gradient(f,uh,fuh::DomainContribution)
  terms = DomainContribution()
  for trian in get_domains(fuh)
    g = _change_argument(gradient,f,trian,uh)
    cell_u = get_cell_dof_values(uh)
    glue = get_glue(trian,Val(num_cell_dims(get_triangulation(uh))))
    cell_id = _compute_cell_ids(uh,trian)
    cell_grad = autodiff_array_gradient(g,cell_u,cell_id)
    add_contribution!(terms,trian,cell_grad)
  end
  terms
end

function _compute_cell_ids(uh,ttrian)
  strian = get_triangulation(uh)
  if strian === ttrian
    return collect(IdentityVector(Int32(num_cells(strian))))
  end
  @check is_change_possible(strian,ttrian)
  D = num_cell_dims(strian)
  sglue = get_glue(strian,Val(D))
  tglue = get_glue(ttrian,Val(D))
  @notimplementedif !isa(sglue,FaceToFaceGlue)
  @notimplementedif !isa(tglue,FaceToFaceGlue)
  scells = IdentityVector(Int32(num_cells(strian)))
  mcells = extend(scells,sglue.mface_to_tface)
  tcells = lazy_map(Reindex(mcells),tglue.tface_to_mface)
  collect(tcells)
end

function jacobian(f::Function,uh::FEFunction)
  fuh = f(uh)
  _jacobian(f,uh,fuh)
end

function _jacobian(f,uh,fuh::AbstractArray)
  @unreachable """\n
  In order to perform AD on a Function taking a FEFunction as argument, such Function
  has to return a DomainContribution.

  Make sure that you are using a Measure instead of a CellQuadrature to perform integration.
  """
end

function _jacobian(f,uh,fuh::DomainContribution)
  terms = DomainContribution()
  for trian in get_domains(fuh)
    g = _change_argument(jacobian,f,trian,uh)
    cell_u = get_cell_dof_values(uh)
    cell_id = _compute_cell_ids(uh,trian)
    cell_grad = autodiff_array_jacobian(g,cell_u,cell_id)
    add_contribution!(terms,trian,cell_grad)
  end
  terms
end

function hessian(f::Function,uh::FEFunction)
  fuh = f(uh)
  _hessian(f,uh,fuh)
end

function _hessian(f,uh,fuh::AbstractArray)
  @unreachable """\n
  In order to perform AD on a Function taking a FEFunction as argument, such Function
  has to return a DomainContribution.

  Make sure that you are using a Measure instead of a CellQuadrature to perform integration.
  """
end

function _hessian(f,uh,fuh::DomainContribution)
  terms = DomainContribution()
  for trian in get_domains(fuh)
    g = _change_argument(hessian,f,trian,uh)
    cell_u = get_cell_dof_values(uh)
    cell_id = _compute_cell_ids(uh,trian)
    cell_grad = autodiff_array_hessian(g,cell_u,cell_id)
    add_contribution!(terms,trian,cell_grad)
  end
  terms
end

function _change_argument(op,f,trian,uh::SingleFieldFEFunction)
  U = get_fe_space(uh)
  function g(cell_u)
    cf = CellField(U,cell_u)
    cell_grad = f(cf)
    get_contribution(cell_grad,trian)
  end
  g
end


text:
include("./imports.jl")

include("metrics.jl")


text:
instructions = String[]
global line = chomp(readline());
while line != ""
    global line
    push!(instructions, line)
    global line = chomp(readline());
end

instructions2 = copy(instructions)
for change in 1:length(instructions)

    global instructions = copy(instructions2)
    changeElems = split(instructions[change])
    changeOp = changeElems[1]
    changeAmount = changeElems[2]
    if changeOp == "nop"
        changeOp = "jmp"
    elseif changeOp == "jmp"
        changeOp = "nop"
    end
    global instructions[change] = changeOp * " " * changeAmount

    called = Int[]
    global inst = 0
    global accumulator = 0
    while !(inst in called) && inst < length(instructions)
        push!(called, inst)
        elems = split(instructions[inst + 1])
        op = elems[1]
        amount = parse(Int, elems[2])
        if op == "acc"
            global accumulator += amount
        end

        if op == "jmp"
            global inst += amount
        else
            global inst += 1
        end
    end

    if inst >= length(instructions)
        println(accumulator)
    end

end

max_stars_repo_name: noeppi-noeppi/aoc

max_issues_repo_name: noeppi-noeppi/aoc

max_forks_repo_name: noeppi-noeppi/aoc


text:
using Compat

# create type MatrixNetwork
type MatrixNetwork{T}
    n::Int64 # number of columns/rows
    rp::Vector{Int64} # row pointers
    ci::Vector{Int64} # column indices
    vals::Vector{T} # corresponding values
end

function MatrixNetwork{T}(A::SparseMatrixCSC{T,Int64})
    At = A'
    return MatrixNetwork(size(At,2),At.colptr,At.rowval,At.nzval)
end

MatrixNetwork(ei::Vector{Int64},ej::Vector{Int64}) =
    MatrixNetwork(ei,ej,max(maximum(ei),maximum(ej)))

function MatrixNetwork(ei::Vector{Int64},ej::Vector{Int64},n::Int64)
    At = sparse(ej,ei,true,n,n);
    return MatrixNetwork(size(At,2),At.colptr,At.rowval,At.nzval)
end

function _matrix_network_direct{T}(A::SparseMatrixCSC{T,Int64})
    return MatrixNetwork(size(A,2),A.colptr,A.rowval,A.nzval)
end

function _matrix_network_direct{T}(A::SparseMatrixCSC{T,Int64},v)
    nzval = ones(typeof(v),length(A.nzval))
    return MatrixNetwork(size(A,2),A.colptr,A.rowval,nzval)
end

import Base.sparse, Base.size

"""
Return back an adjacency matrix representation
of the transpose. This requires no work.
"""
function sparse_transpose{T}(A::MatrixNetwork{T})
    return SparseMatrixCSC(A.n,A.n,A.rp,A.ci,A.vals)
end

"""
Return back an adjacency matrix representation
of the current MatrixNetwork
"""
function sparse{T}(A::MatrixNetwork{T})
    return sparse_transpose(A)'
end

function size(A::MatrixNetwork)
    return (A.n,A.n)
end

function size(A::MatrixNetwork, dim::Integer)
    if dim == 1 || dim == 2
        return A.n
    elseif dim > 2
        return 1
    else
        throw(DomainError())
    end
end

"""
`is_undirected`
===============

Check the matrix associated with a matrix network
for symmetry.

Input
-----
- `A`: a matrix network

Returns
-------
- `bool` with true indicating the network is undirected
    and the matrix is symmetric
"""
function is_undirected end

function is_undirected(A::MatrixNetwork)
   M = sparse_transpose(A)
   return issymmetric(M)
end

function is_undirected(A::SparseMatrixCSC)
   return issymmetric(A)
end

"""
`is_connected`
==============

Check the matrix associated with a matrix network
for (strong) connectivity

Usage
-----
- `is_connected(A)`

Input
-----
- `A`: a `MatrixNetwork` or `SparseMatrixCSC` class

Returns
-------
- `bool` with true indicating the matrix is strongly connected
and false indicating
"""
function is_connected end

function is_connected(A::MatrixNetwork)
    return maximum(scomponents(A).map) == 1
end

function is_connected(A::SparseMatrixCSC)
    return maximum(scomponents(A).map) == 1
end


text:
println("Testing...")

using Physt

h = h1([1,2,3,4,5,6,7], 2)

h = h1(rand(100), 12)
println(h)

h_ = h << .4
println(h_)
h__ = h_ << rand(20)

println(h)
update!(h, .40)
println(h)


text:
import SparseArrays

"""
    abstract type AbstractIndexing end

Indexing to be used for storing the row and column indices of
`MutableSparseMatrixCSC`. See [`ZeroBasedIndexing`](@ref) and
[`OneBasedIndexing`](@ref).
"""
abstract type AbstractIndexing end

"""
    struct ZeroBasedIndexing <: AbstractIndexing end

Zero-based indexing: the `i`th row or column has index `i - 1`. This is useful
when the vectors of row and column indices need to be communicated to a
library using zero-based indexing such as C libraries.
"""
struct ZeroBasedIndexing <: AbstractIndexing end

"""
    struct ZeroBasedIndexing <: AbstractIndexing end

One-based indexing: the `i`th row or column has index `i`. This enables an
allocation-free conversion of [`MutableSparseMatrixCSC`](@ref) to
`SparseArrays.SparseMatrixCSC`.
"""
struct OneBasedIndexing <: AbstractIndexing end

_first_index(::ZeroBasedIndexing) = 0
_first_index(::OneBasedIndexing) = 1

_shift(x, ::T, ::T) where {T<:AbstractIndexing} = x
_shift(x::Integer, ::ZeroBasedIndexing, ::OneBasedIndexing) = x + 1
_shift(x::Integer, ::OneBasedIndexing, ::ZeroBasedIndexing) = x - 1
_shift(x::Array{<:Integer}, ::ZeroBasedIndexing, ::OneBasedIndexing) = x .+ 1

"""
    mutable struct MutableSparseMatrixCSC{Tv,Ti<:Integer,I<:AbstractIndexing}
        indexing::I
        m::Int
        n::Int
        colptr::Vector{Ti}
        rowval::Vector{Ti}
        nzval::Vector{Tv}
    end

Matrix type loading sparse matrices in the Compressed Sparse Column format.
The indexing used is `indexing`, see [`AbstractIndexing`](@ref). The other
fields have the same meaning than for `SparseArrays.SparseMatrixCSC` except
that the indexing is different unless `indexing` is `OneBasedIndexing`.

The matrix is loaded in 5 steps:
1) `MOI.empty!` is called.
2) `MOI.Utilities.add_column` and `MOI.Utilities.allocate_terms` are called in
   any order.
3) `MOI.Utilities.set_number_of_rows` is called.
4) `MOI.Utilities.load_terms` is called for each affine function.
5) `MOI.Utilities.final_touch` is called.
"""
mutable struct MutableSparseMatrixCSC{Tv,Ti<:Integer,I<:AbstractIndexing}
    indexing::I
    m::Int # Number of rows
    n::Int # Number of columns
    colptr::Vector{Ti}
    rowval::Vector{Ti}
    nzval::Vector{Tv}
    function MutableSparseMatrixCSC{Tv,Ti,I}() where {Tv,Ti<:Integer,I}
        return new{Tv,Ti,I}(I(), 0, 0, [zero(Ti)], Ti[], Tv[])
    end
end

function SparseArrays.nzrange(A::MutableSparseMatrixCSC, col)
    start = _shift(A.colptr[col], A.indexing, OneBasedIndexing())
    stop = _shift(A.colptr[col+1], A.indexing, ZeroBasedIndexing())
    return start:stop
end

function MOI.empty!(A::MutableSparseMatrixCSC{Tv,Ti}) where {Tv,Ti}
    A.m = 0
    A.n = 0
    resize!(A.colptr, 1)
    A.colptr[1] = zero(Ti)
    empty!(A.rowval)
    empty!(A.nzval)
    return
end

function add_column(A::MutableSparseMatrixCSC{Tv,Ti}) where {Tv,Ti}
    A.n += 1
    push!(A.colptr, zero(Ti))
    return
end

function add_columns(A::MutableSparseMatrixCSC{Tv,Ti}, n) where {Tv,Ti}
    A.n += n
    for _ in 1:n
        push!(A.colptr, zero(Ti))
    end
    return
end

function set_number_of_rows(A::MutableSparseMatrixCSC, num_rows)
    A.m = num_rows
    for i in 3:length(A.colptr)
        A.colptr[i] += A.colptr[i-1]
    end
    resize!(A.rowval, A.colptr[end])
    resize!(A.nzval, A.colptr[end])
    return
end

function final_touch(A::MutableSparseMatrixCSC)
    for i in length(A.colptr):-1:2
        A.colptr[i] = _shift(A.colptr[i-1], ZeroBasedIndexing(), A.indexing)
    end
    A.colptr[1] = _first_index(A.indexing)
    return
end

_variable(term::MOI.ScalarAffineTerm) = term.variable
_variable(term::MOI.VectorAffineTerm) = _variable(term.scalar_term)

function allocate_terms(
    A::MutableSparseMatrixCSC,
    index_map,
    func::Union{MOI.ScalarAffineFunction,MOI.VectorAffineFunction},
)
    for term in func.terms
        A.colptr[index_map[_variable(term)].value+1] += 1
    end
    return
end

_output_index(::MOI.ScalarAffineTerm) = 1
_output_index(term::MOI.VectorAffineTerm) = term.output_index

function load_terms(
    A::MutableSparseMatrixCSC,
    index_map,
    func::Union{MOI.ScalarAffineFunction,MOI.VectorAffineFunction},
    offset::Int,
)
    new_offset = _shift(offset, OneBasedIndexing(), A.indexing)
    for term in func.terms
        ptr = A.colptr[index_map[_variable(term)].value] += 1
        A.rowval[ptr] = new_offset + _output_index(term)
        A.nzval[ptr] = MOI.coefficient(term)
    end
    return
end

"""
    Base.convert(
        ::Type{SparseMatrixCSC{Tv,Ti}},
        A::MutableSparseMatrixCSC{Tv,Ti,I},
    ) where {Tv,Ti,I}

Converts `A` to a `SparseMatrixCSC`. Note that the field `A.nzval` is **not
copied** so if `A` is modified after the call of this function, it can still
affect the value returned. Moreover, if `I` is `OneBasedIndexing`, `colptr`
and `rowval` are not copied either, i.e., the conversion is allocation-free.
"""
function Base.convert(
    ::Type{SparseArrays.SparseMatrixCSC{Tv,Ti}},
    A::MutableSparseMatrixCSC{Tv,Ti},
) where {Tv,Ti}
    return SparseArrays.SparseMatrixCSC{Tv,Ti}(
        A.m,
        A.n,
        _shift(A.colptr, A.indexing, OneBasedIndexing()),
        _shift(A.rowval, A.indexing, OneBasedIndexing()),
        A.nzval,
    )
end

function _first_in_column(
    A::MutableSparseMatrixCSC{Tv,Ti},
    row::Integer,
    col::Integer,
) where {Tv,Ti}
    range = SparseArrays.nzrange(A, col)
    row = _shift(row, OneBasedIndexing(), A.indexing)
    idx = searchsortedfirst(view(A.rowval, range), row)
    return get(range, idx, last(range) + 1)
end

function extract_function(
    A::MutableSparseMatrixCSC{T},
    row::Integer,
    constant::T,
) where {T}
    func = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm{T}[], constant)
    for col in 1:A.n
        idx = _first_in_column(A, row, col)
        if idx > last(SparseArrays.nzrange(A, col))
            continue
        end
        r = _shift(A.rowval[idx], A.indexing, OneBasedIndexing())
        if r == row
            push!(
                func.terms,
                MOI.ScalarAffineTerm(A.nzval[idx], MOI.VariableIndex(col)),
            )
        end
    end
    return func
end

function extract_function(
    A::MutableSparseMatrixCSC{T},
    rows::UnitRange,
    constants::Vector{T},
) where {T}
    func = MOI.VectorAffineFunction(MOI.VectorAffineTerm{T}[], constants)
    if isempty(rows)
        return func
    end
    idx = [_first_in_column(A, first(rows), col) for col in 1:A.n]
    for output_index in eachindex(rows)
        for col in 1:A.n
            if idx[col] > last(SparseArrays.nzrange(A, col))
                continue
            end
            row = _shift(A.rowval[idx[col]], A.indexing, OneBasedIndexing())
            if row != rows[output_index]
                continue
            end
            push!(
                func.terms,
                MOI.VectorAffineTerm(
                    output_index,
                    MOI.ScalarAffineTerm(
                        A.nzval[idx[col]],
                        MOI.VariableIndex(col),
                    ),
                ),
            )
            idx[col] += 1
        end
    end
    return func
end


text:
using HttpServer

http = HttpHandler() do req::Request, res::Response
    Response(strig("Hello, World!"))
end

http.events["error"]  = ( client, err ) -> println( err )
http.events["listen"] = ( port )        -> println("Listening on $port...")

server = Server( http )
run( server, 2000 )

max_stars_repo_name: est/snippets

max_issues_repo_name: est/snippets

max_forks_repo_name: est/snippets


text:
module TestMutator
  using FactCheck
  import Creature
  import Mutator
  import Config
  import Helper
  import Code

  # #
  # # Just a container for testing data
  # #
  # type Data
  #   d::Array{Any, 1}
  # end
  # #
  # # Global data container instance
  # #
  # global _d = Data(Array{Any, 1}(10))
  # #
  # # Changes mutation parameters to supress any mutation type
  # # except adding
  # #
  # function setAddMutation()
  #   _d.d[1] = Config.val(:orgMutationProbs)
  #   Config.val(:orgMutationProbs, [1,0,0,0,0,0,0])
  # end
  # function restoreAddMutation()	Config.val(:orgMutationProbs, _d.d[1]) end
  # #
  # # Sets maximum function parameters to 1. Related to add mutation.
  # #
  # function setMaxFuncParams()
  #   _d.d[2] = Config.val(:codeFuncParamAmount)
  #   Config.val(:codeFuncParamAmount, 1)
  # end
  # function restoreMaxFuncParams() Config.val(:codeFuncParamAmount, _d.d[2]) end
  # #
  # # pos - "1,2,3" fn - 1, block - 2, line - 3
  # #
  # function setMutationPos(pos::String) Mutator.Code.eval(parse("function getRandPos(org::Creature.Organism) Helper.CodePos(" * pos * ") end")) end
  # function setCodePart(part::String) Mutator.Code.eval(parse("const CODE_PARTS = [CodePart(" * part * ")]")) end
  # function restoreMutator() include("../src/organism/Mutator.jl") end
  # #
  # # Checks if one of supported code part types was added into AST
  # #
  # function isVarAdded(org)
  #   local hasType = false
  #   for i in [Float16] if length(org.funcs[1].blocks[1].vars[i]) > 0 hasType = true end end
  #   hasType
  # end

  # facts("Adding var to empty code") do
  #   setAddMutation()
  #   setMutationPos("1,1,1")
  #   setCodePart("Code.plus, false")

  #   org = Creature.create()
  #   @fact Mutator._onAdd(org)                         --> true
  #   @fact org.code.args[2].args[1].head               --> :local
  #   @fact Helper.getLines(org.code, [2,1,1,1])[1]     --> :var_1

  #   restoreAddMutation()
  #   restoreMutator()
  # end
  # facts("Adding two vars to empty code") do
  #   setAddMutation()
  #   setMutationPos("1,1,1")
  #   setCodePart("Code.plus, false")

  #   org = Creature.create()
  #   @fact Mutator._onAdd(org)                         --> true
  #   @fact Mutator._onAdd(org)                         --> true
  #   @fact org.code.args[2].args[1].head               --> :local
  #   @fact org.code.args[2].args[2].head               --> :local
  #   @fact Helper.getLines(org.code, [2,1,1,1])[1]     --> :var_2
  #   @fact Helper.getLines(org.code, [2,2,1,1])[1]     --> :var_1

  #   restoreAddMutation()
  #   restoreMutator()
  # end
  # facts("Adding var into the function") do
  #   setAddMutation()
  #   setMutationPos("1,1,1")
  #   setCodePart("Code.fn, true")
  #   setMaxFuncParams()

  #   org = Creature.create()
  #   @fact Mutator._onAdd(org)                         --> true
  #   @fact Helper.getLines(org.code, [2,1,1])[1]       --> :func_2
  #   setMutationPos("2,1,1")
  #   println(Code.getRandPos(org))
  #   setCodePart("Code.plus, false")
  #   @fact Mutator._onAdd(org)                         --> true
  #   #println(Helper.getLines(org.code, [2,2,1,1])[1]     --> :var_1)

  #   restoreMaxFuncParams()
  #   restoreAddMutation()
  #   restoreMutator()
  # end
  # facts("Adding function declaration") do
  #   setAddMutation()
  #   setMutationPos("1,1,1")
  #   setCodePart("Code.fn, true")
  #   setMaxFuncParams()

  #   org = Creature.create()
  #   @fact Mutator._onAdd(org)                         --> true
  #   @fact Helper.getLines(org.code, [2,1,1])[1]       --> :func_2

  #   restoreMaxFuncParams()
  #   restoreAddMutation()
  #   restoreMutator()
  # end
end

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


text:
# This file was generated, do not modify it. # hide
using Gnuplot
x=[0,1,2]
y=[0,1,2]
Z=[10 10 10; 10 3 10; 10 2 10]
@gsp x y Z "w l lc 'red'"
save(term="png", output="plt3d_ex1.png")


text:
module hex_ip_constr_examples
using FinEtools
using FinEtools.MeshExportModule.CSV: savecsv
using FinEtoolsDeforLinear
using FinEtoolsDeforLinear.AlgoDeforLinearModule
using FinEtools.AlgoBaseModule: evalconvergencestudy
using FinEtoolsDeforLinear.AlgoDeforLinearModule: linearstatics, exportstresselementwise, exportstress
using Statistics: mean
using LinearAlgebra: Symmetric, cholesky, eigen, svd
using PGFPlotsX

mesh() = (FinEtools.FENodeSetModule.FENodeSet([
	0.0 0.0 0.0;
	1.0 0.0 0.0;
	1.0 1.0 0.0;
	0.0 1.0 0.0;
	0.0 0.0 1.0;
	1.0 0.0 1.0;
	1.0 1.0 1.0
	0.0 1.0 1.0;
	]), FinEtools.FESetModule.FESetH8(reshape([1, 2, 3, 4, 5, 6, 7, 8], 1, 8)))

function hex_ip_constr_par()
	fens,fes = mesh()

	# integration_rule = TrapezoidalRule(3)
	integration_rule = GaussRule(3, 2)
	pc = integration_rule.param_coords;
	w  =  integration_rule.weights
	npts = integration_rule.npts;

	Q = fill(0.0, npts, 24)
	for j in 1:npts
		gradNpar = bfundpar(fes, vec(pc[j,:]))
		k = 1
		for i in 1:size(gradNpar, 1)
			Q[j, k:(k+2)] .= gradNpar[i, :]
			k = k + 3
		end
	end

	decomp = svd(Q)
	@show decomp.S

		true
end # hex_ip_constr_par

function hex_ip_constr_xyz()
	xyzperturbation = [
	0.0767656  -0.983206    -0.14343
	0.45767     0.981479     0.450997
	-0.295854    0.542922     0.321333
	-0.85204    -0.97824     -0.772874
	-0.0539756   0.994907     0.822798
	0.447173    0.528742     0.0445352
	-0.468417    0.00673427   0.953151
	-0.898513   -0.915871     0.933237  ] ./ 5.0
	fens,fes = mesh()
	fens.xyz .+= xyzperturbation

	# integration_rule = TrapezoidalRule(3)
	integration_rule = GaussRule(3, 2)
	pc = integration_rule.param_coords;
	w  =  integration_rule.weights
	npts = integration_rule.npts;

	Q = fill(0.0, npts, 24)
	for j in 1:npts
		gradNpar = bfundpar(fes, vec(pc[j,:]))
		J = (transpose(fens.xyz) * gradNpar)
		gradN = gradNpar / J
		k = 1
		for i in 1:size(gradN, 1)
			Q[j, k:(k+2)] .= gradN[i, :]
			k = k + 3
		end
	end

	decomp = svd(Q)
	@show decomp.S

	true
end # hex_ip_constr_xyz

function allrun()
	println("#####################################################")
	println("# hex_ip_constr_par ")
	hex_ip_constr_par()
	println("#####################################################")
	println("# hex_ip_constr_xyz ")
	hex_ip_constr_xyz()
	return true
end # function allrun

end # module hex_ip_constr_examples


text:
using DrWatson
quickactivate(@__DIR__, "Watson")
using DynamicalBilliards, PyPlot, LinearAlgebra

include(srcdir("plot_perturbationgrowth.jl"))
include(srcdir("unitcells.jl"))

@tagsave(datadir("mushrooms, "Λ_N=$N.bson"), (@dict Λ Λσ ws hs description))


text:
### BasicMCRange

# BasicMCRange is used for sampling a single parameter via serial Monte Carlo
# It is the most elementary and typical Markov chain Monte Carlo (MCMC) range
# It contains nagivational info (burnin, thinning, number of steps)

struct BasicMCRange{T<:Integer} <: MCRange{T}
  burnin::T
  thinning::T
  nsteps::T
  postrange::StepRange{T, T}
  npoststeps::T

  function BasicMCRange{T}(postrange::StepRange{T, T}) where T<:Integer
    burnin = first(postrange)-1
    thinning = postrange.step
    nsteps = last(postrange)

    @assert burnin >= 0 "Number of burn-in iterations should be non-negative"
    @assert thinning >= 1 "Thinning should be >= 1"
    @assert nsteps > burnin "Total number of MCMC iterations should be greater than number of burn-in iterations"

    npoststeps = length(postrange)

    new(burnin, thinning, nsteps, postrange, npoststeps)
  end
end

BasicMCRange(postrange::StepRange{T, T}) where T<:Integer = BasicMCRange{T}(postrange)

BasicMCRange(postrange::UnitRange{T}) where {T<:Integer} = BasicMCRange{T}(first(postrange):1:last(postrange))

BasicMCRange(; burnin::T=0, thinning::T=1, nsteps::T=100) where {T<:Integer} = BasicMCRange{T}((burnin+1):thinning:nsteps)

show(io::IO, r::BasicMCRange) =
  print(io, "BasicMCRange: number of steps = $(r.nsteps), burnin = $(r.burnin), thinning = $(r.thinning)")