text:
function part1()
    numbers = [parse(Int, i) for i in readlines("inputD1.txt")]
    return sum(numbers[2:end] .> numbers[1:end-1])
end

function part2()
    numbers = [parse(Int, i) for i in readlines("inputD1.txt")]
    return sum(numbers[4:end] .> numbers[1:end-3])
end

println(part1())
println(part2())


text:
include("aesgf-code.jl")

global const Nks = [4, 6, 8]
global const Nbs = [4, 4, 4]
global const Nb = 4
global const Nrs = [10, 12, 14]
global const WORDLENGTH = 4

global const SBOX = [
0x63, 0x7c, 0x77, 0x7b, 0xf2, 0x6b, 0x6f, 0xc5, 0x30, 0x01, 0x67, 0x2b, 0xfe, 0xd7, 0xab, 0x76,
0xca, 0x82, 0xc9, 0x7d, 0xfa, 0x59, 0x47, 0xf0, 0xad, 0xd4, 0xa2, 0xaf, 0x9c, 0xa4, 0x72, 0xc0,
0xb7, 0xfd, 0x93, 0x26, 0x36, 0x3f, 0xf7, 0xcc, 0x34, 0xa5, 0xe5, 0xf1, 0x71, 0xd8, 0x31, 0x15,
0x04, 0xc7, 0x23, 0xc3, 0x18, 0x96, 0x05, 0x9a, 0x07, 0x12, 0x80, 0xe2, 0xeb, 0x27, 0xb2, 0x75,
0x09, 0x83, 0x2c, 0x1a, 0x1b, 0x6e, 0x5a, 0xa0, 0x52, 0x3b, 0xd6, 0xb3, 0x29, 0xe3, 0x2f, 0x84,
0x53, 0xd1, 0x00, 0xed, 0x20, 0xfc, 0xb1, 0x5b, 0x6a, 0xcb, 0xbe, 0x39, 0x4a, 0x4c, 0x58, 0xcf,
0xd0, 0xef, 0xaa, 0xfb, 0x43, 0x4d, 0x33, 0x85, 0x45, 0xf9, 0x02, 0x7f, 0x50, 0x3c, 0x9f, 0xa8,
0x51, 0xa3, 0x40, 0x8f, 0x92, 0x9d, 0x38, 0xf5, 0xbc, 0xb6, 0xda, 0x21, 0x10, 0xff, 0xf3, 0xd2,
0xcd, 0x0c, 0x13, 0xec, 0x5f, 0x97, 0x44, 0x17, 0xc4, 0xa7, 0x7e, 0x3d, 0x64, 0x5d, 0x19, 0x73,
0x60, 0x81, 0x4f, 0xdc, 0x22, 0x2a, 0x90, 0x88, 0x46, 0xee, 0xb8, 0x14, 0xde, 0x5e, 0x0b, 0xdb,
0xe0, 0x32, 0x3a, 0x0a, 0x49, 0x06, 0x24, 0x5c, 0xc2, 0xd3, 0xac, 0x62, 0x91, 0x95, 0xe4, 0x79,
0xe7, 0xc8, 0x37, 0x6d, 0x8d, 0xd5, 0x4e, 0xa9, 0x6c, 0x56, 0xf4, 0xea, 0x65, 0x7a, 0xae, 0x08,
0xba, 0x78, 0x25, 0x2e, 0x1c, 0xa6, 0xb4, 0xc6, 0xe8, 0xdd, 0x74, 0x1f, 0x4b, 0xbd, 0x8b, 0x8a,
0x70, 0x3e, 0xb5, 0x66, 0x48, 0x03, 0xf6, 0x0e, 0x61, 0x35, 0x57, 0xb9, 0x86, 0xc1, 0x1d, 0x9e,
0xe1, 0xf8, 0x98, 0x11, 0x69, 0xd9, 0x8e, 0x94, 0x9b, 0x1e, 0x87, 0xe9, 0xce, 0x55, 0x28, 0xdf,
0x8c, 0xa1, 0x89, 0x0d, 0xbf, 0xe6, 0x42, 0x68, 0x41, 0x99, 0x2d, 0x0f, 0xb0, 0x54, 0xbb, 0x16
]

global const INVSBOX = [
0x52, 0x09, 0x6a, 0xd5, 0x30, 0x36, 0xa5, 0x38, 0xbf, 0x40, 0xa3, 0x9e, 0x81, 0xf3, 0xd7, 0xfb,
0x7c, 0xe3, 0x39, 0x82, 0x9b, 0x2f, 0xff, 0x87, 0x34, 0x8e, 0x43, 0x44, 0xc4, 0xde, 0xe9, 0xcb,
0x54, 0x7b, 0x94, 0x32, 0xa6, 0xc2, 0x23, 0x3d, 0xee, 0x4c, 0x95, 0x0b, 0x42, 0xfa, 0xc3, 0x4e,
0x08, 0x2e, 0xa1, 0x66, 0x28, 0xd9, 0x24, 0xb2, 0x76, 0x5b, 0xa2, 0x49, 0x6d, 0x8b, 0xd1, 0x25,
0x72, 0xf8, 0xf6, 0x64, 0x86, 0x68, 0x98, 0x16, 0xd4, 0xa4, 0x5c, 0xcc, 0x5d, 0x65, 0xb6, 0x92,
0x6c, 0x70, 0x48, 0x50, 0xfd, 0xed, 0xb9, 0xda, 0x5e, 0x15, 0x46, 0x57, 0xa7, 0x8d, 0x9d, 0x84,
0x90, 0xd8, 0xab, 0x00, 0x8c, 0xbc, 0xd3, 0x0a, 0xf7, 0xe4, 0x58, 0x05, 0xb8, 0xb3, 0x45, 0x06,
0xd0, 0x2c, 0x1e, 0x8f, 0xca, 0x3f, 0x0f, 0x02, 0xc1, 0xaf, 0xbd, 0x03, 0x01, 0x13, 0x8a, 0x6b,
0x3a, 0x91, 0x11, 0x41, 0x4f, 0x67, 0xdc, 0xea, 0x97, 0xf2, 0xcf, 0xce, 0xf0, 0xb4, 0xe6, 0x73,
0x96, 0xac, 0x74, 0x22, 0xe7, 0xad, 0x35, 0x85, 0xe2, 0xf9, 0x37, 0xe8, 0x1c, 0x75, 0xdf, 0x6e,
0x47, 0xf1, 0x1a, 0x71, 0x1d, 0x29, 0xc5, 0x89, 0x6f, 0xb7, 0x62, 0x0e, 0xaa, 0x18, 0xbe, 0x1b,
0xfc, 0x56, 0x3e, 0x4b, 0xc6, 0xd2, 0x79, 0x20, 0x9a, 0xdb, 0xc0, 0xfe, 0x78, 0xcd, 0x5a, 0xf4,
0x1f, 0xdd, 0xa8, 0x33, 0x88, 0x07, 0xc7, 0x31, 0xb1, 0x12, 0x10, 0x59, 0x27, 0x80, 0xec, 0x5f,
0x60, 0x51, 0x7f, 0xa9, 0x19, 0xb5, 0x4a, 0x0d, 0x2d, 0xe5, 0x7a, 0x9f, 0x93, 0xc9, 0x9c, 0xef,
0xa0, 0xe0, 0x3b, 0x4d, 0xae, 0x2a, 0xf5, 0xb0, 0xc8, 0xeb, 0xbb, 0x3c, 0x83, 0x53, 0x99, 0x61,
0x17, 0x2b, 0x04, 0x7e, 0xba, 0x77, 0xd6, 0x26, 0xe1, 0x69, 0x14, 0x63, 0x55, 0x21, 0x0c, 0x7d
]

global const MIXCOLUMNSMATRIX = [
0x02, 0x03, 0x01, 0x01,
0x01, 0x02, 0x03, 0x01,
0x01, 0x01, 0x02, 0x03,
0x03, 0x01, 0x01, 0x02
]

global const INVMIXCOLUMNSMATRIX = [
0x0e, 0x0b, 0x0d, 0x09,
0x09, 0x0e, 0x0b, 0x0d,
0x0d, 0x09, 0x0e, 0x0b,
0x0b, 0x0d, 0x09, 0x0e
]

# Encrypts the given plaintext block using the given key
# and returns the resulting ciphertext.
# Both the plaintext and key must be represented as strings consisting of
# only hexadecimal characters.
# The returned ciphertext is a string consisting of only hexadecimal
# characters.
function AESEncrypt(plain::String, key::String)
	bytes2hex(AESEncrypt(hex2bytes(plain), hex2bytes(key)))
end

# Decrypts the given ciphertext block using the given key
# and returns the resulting plaintext.
# Both the ciphertext and key must be represented as strings consisting of
# only hexadecimal characters.
# The returned plaintext is a string consisting of only hexadecimal
# characters.
function AESDecrypt(cipher::String, key::String)
	bytes2hex(AESDecrypt(hex2bytes(cipher), hex2bytes(key)))
end

# Encrypts the given plaintext block using the given key
# and returns the resulting ciphertext.
# Both the plaintext and key are arrays holding bytes of type UInt8.
# The returned ciphertext is an array holding bytes of type UInt8.
function AESEncrypt(plain::Array{UInt8, 1}, key::Array{UInt8, 1})
	(w, Nr) = AEScrypt(plain, key)
	return AESCipher(plain, w, Nr)
end

# Decrypts the given ciphertext block using the given key
# and returns the resulting plaintext.
# Both the ciphertext and key are arrays holding bytes of type UInt8.
# The returned plaintext is an array holding bytes of type UInt8.
function AESDecrypt(cipher::Array{UInt8, 1}, key::Array{UInt8, 1})
	(w, Nr) = AEScrypt(cipher, key)
	return AESInvCipher(cipher, w, Nr)
end

# General function for checking the parameters of the encryption/decryption
# functions and expanding the key.
function AEScrypt(input::Array{UInt8, 1}, key::Array{UInt8, 1})
	if length(input) != (WORDLENGTH * Nb)
		error("input must be a 16-byte block!")
	end
	(Nk, Nr) = AESParameters(key)
	w = KeyExpansion(key, Nk, Nr)
	return (w, Nr)
end

# Return a tuple (Nk, Nr) where Nk is the number of key blocks
# and Nr is the number of rounds for the key size.
function AESParameters(key::Array{UInt8, 1})
	if mod(length(key), WORDLENGTH) != 0
		error("the key length must be a multiple of four!")
	end

	Nk = div(length(key), WORDLENGTH)
	i = first(indexin([Nk], Nks))

	if i == 0
		error("key length is non-standard!")
	end

	Nr = Nrs[i]
	return (Nk, Nr)
end

function KeyExpansion(key::Array{UInt8, 1}, Nk::Int, Nr::Int)
	@assert(length(key) == (WORDLENGTH * Nk))

	w = Array{UInt8}(undef, WORDLENGTH * Nb * (Nr + 1))
	w[1:(WORDLENGTH * Nk)] = copy(key)
	i = Nk

	while i < (Nb * (Nr + 1))
		temp = w[((i - 1) * WORDLENGTH + 1):(i * WORDLENGTH)]
		if mod(i, Nk) == 0
			temp = xor.(SubWord(RotWord(temp)) , Rcon(div(i, Nk)))
		elseif (Nk > 6) && (mod(i, Nk) == Nb)
			temp = SubWord(temp)
		end
		w[(i * WORDLENGTH + 1):((i + 1) * WORDLENGTH)] = xor.(w[((i - Nk) * WORDLENGTH + 1):((i - Nk + 1) * WORDLENGTH)] , temp)
		i += 1
	end

	return w
end

function SubWord(w::Array{UInt8, 1})
	@assert(length(w) == WORDLENGTH)
	map!(x -> SBOX[Int(x) + 1], w, w)
	return w
end

function RotWord(w::Array{UInt8, 1})
	@assert(length(w) == WORDLENGTH)
	permute!(w, [2, 3, 4, 1])
end

function Rcon(i::Int)
	@assert(i > 0)
	x = 0x01
	for j=1:(i-1)
		x = gmul(x, 0x02)
	end
	return [x, 0x00, 0x00, 0x00]
end

function AESCipher(inBytes::Array{UInt8, 1}, w::Array{UInt8, 1}, Nr::Int)
	@assert(WORDLENGTH == Nb)
	@assert(length(inBytes) == (WORDLENGTH * Nb))
	@assert(length(w) == (WORDLENGTH * Nb * (Nr + 1)))

	state = copy(inBytes)
	AddRoundKey(state, w[1:(Nb * WORDLENGTH)])

	for round=1:(Nr-1)
		SubBytes(state)
		ShiftRows(state)
		MixColumns(state)
		AddRoundKey(state, w[(round * Nb * WORDLENGTH + 1):((round + 1) * Nb * WORDLENGTH)])
	end

	SubBytes(state)
	ShiftRows(state)
	AddRoundKey(state, w[(Nr * Nb * WORDLENGTH + 1):((Nr + 1) * Nb * WORDLENGTH)])

	return state
end

function AESInvCipher(inBytes::Array{UInt8, 1}, w::Array{UInt8, 1}, Nr::Int)
	@assert(WORDLENGTH == Nb)
	@assert(length(inBytes) == (WORDLENGTH * Nb))
	@assert(length(w) == (WORDLENGTH * Nb * (Nr + 1)))

	state = copy(inBytes)
	AddRoundKey(state, w[(Nr * Nb * WORDLENGTH + 1):((Nr + 1) * Nb * WORDLENGTH)])

	for round=(Nr-1):-1:1
		InvShiftRows(state)
		InvSubBytes(state)
		AddRoundKey(state, w[(round * Nb * WORDLENGTH + 1):((round + 1) * Nb * WORDLENGTH)])
		InvMixColumns(state)
	end

	InvShiftRows(state)
	InvSubBytes(state)
	AddRoundKey(state, w[1:(Nb * WORDLENGTH)])

	return state
end

# Returns one-based column indices of the given column (one-based).
function columnIndices(column::Int)
	return column:Nb:Nb*Nb
end

# Returns row indices of the given row (one-based).
function rowIndices(row::Int)
	return ((row - 1) * Nb + 1):(row * Nb)
end

function SubBytes(a::Array{UInt8, 1})
	SubBytesGen(a, false)
end

function InvSubBytes(a::Array{UInt8, 1})
	SubBytesGen(a, true)
end

function SubBytesGen(a::Array{UInt8, 1}, inv::Bool)
	box = inv ? INVSBOX : SBOX
	f = (x) -> box[Int(x) + 1]
	map!(f, a, a)
	return a
end

function ShiftRows(a::Array{UInt8, 1})
	ShiftRowsGen(a, false)
end

function InvShiftRows(a::Array{UInt8, 1})
	ShiftRowsGen(a, true)
end

function ShiftRowsGen(a::Array{UInt8, 1}, inv::Bool)
	# permute columns in _last three_ rows
	# note that a column is actually a row in memory
	for r=2:Nb
		indices = columnIndices(r)
		step = inv ? Nb + 2 - r : r
		p = map(c -> mod(c + step - 1, Nb) + 1, 0:(Nb-1))
		a[indices] = a[indices][p]
	end
	return a
end

function MixColumns(a::Array{UInt8, 1})
	MixColumnsGen(a, false)
end

function InvMixColumns(a::Array{UInt8, 1})
	MixColumnsGen(a, true)
end

function MixColumnsGen(a::Array{UInt8, 1}, inv::Bool)
	# note that columns are actually rows in memory
	matrix = inv ? INVMIXCOLUMNSMATRIX : MIXCOLUMNSMATRIX
	for c=1:Nb
		indices = rowIndices(c)
		ai = copy(a[indices])
		@assert(length(ai) == Nb)
		# Matrix multiplication with Galois field operations
		for r=1:Nb
			mi = matrix[rowIndices(r)]
			a[indices[r]] = gadd(map(gmul, ai, mi))
		end
	end
	return a
end

function AddRoundKey(s::Array{UInt8, 1}, w::Array{UInt8, 1})
	@assert(length(w) == (WORDLENGTH * Nb) && (WORDLENGTH == Nb))
	return map!(gadd, s, s, w)
end

ext: jl


text:
# Things that are useful to know across platforms
if Sys.iswindows()
    const LIBPATH_env = "PATH"
    const LIBPATH_default = ""
    const pathsep = ';'
elseif Sys.isapple()
    const LIBPATH_env = "DYLD_FALLBACK_LIBRARY_PATH"
    const LIBPATH_default = "~/lib:/usr/local/lib:/lib:/usr/lib"
    const pathsep = ':'
else
    const LIBPATH_env = "LD_LIBRARY_PATH"
    const LIBPATH_default = ""
    const pathsep = ':'
end

function adjust_ENV!(env::Dict, PATH::String, LIBPATH::String, adjust_PATH::Bool, adjust_LIBPATH::Bool)
    if adjust_PATH
        if !isempty(get(env, "PATH", ""))
            env["PATH"] = string(PATH, pathsep, env["PATH"])
        else
            env["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(env, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env[LIBPATH_env] = string(LIBPATH, pathsep, LIBPATH_base)
        else
            env[LIBPATH_env] = LIBPATH
        end
    end
    return env
end

function withenv_executable_wrapper(f::Function,
                                    executable_path::String,
                                    PATH::String,
                                    LIBPATH::String,
                                    adjust_PATH::Bool,
                                    adjust_LIBPATH::Bool)
    env = Dict{String,String}(
        "PATH" => get(ENV, "PATH", ""),
        LIBPATH_env => get(ENV, LIBPATH_env, ""),
    )
    env = adjust_ENV!(env, PATH, LIBPATH, adjust_PATH, adjust_LIBPATH)
    withenv(env...) do
        f(executable_path)
    end
end

function dev_jll(src_name)
    # Grab `Pkg` module
    Pkg = first(filter(p-> p[1].name == "Pkg", Base.loaded_modules))[2]

    # First, `dev` out the package, but don't affect the current project
    mktempdir() do temp_env
        Pkg.activate(temp_env) do
            Pkg.develop("$(src_name)_jll")
            Pkg.instantiate()
        end
    end

    # Create the override directory by populating it with the artifact contents
    override_dir = joinpath(Pkg.devdir(), "$(src_name)_jll", "override")
    if !isdir(override_dir)
        artifacts_toml = joinpath(Pkg.devdir(), "$(src_name)_jll", "Artifacts.toml")
        art_hash = Pkg.Artifacts.artifact_hash(src_name, artifacts_toml)
        art_location = Pkg.Artifacts.artifact_path(art_hash)
        cp(art_location, override_dir)
    end
    # Force recompilation of that package, just in case it wasn't dev'ed before
    touch(joinpath(Pkg.devdir(), "$(src_name)_jll", "src", "$(src_name)_jll.jl"))
    @info("$(src_name)_jll dev'ed out to $(joinpath(Pkg.devdir(), "$(src_name)_jll")) with pre-populated override directory")
end

const JULIA_LIBDIRS = String[]
"""
    get_julia_libpaths()

Return the library paths that e.g. libjulia and such are stored in.
"""
function get_julia_libpaths()
    if isempty(JULIA_LIBDIRS)
        append!(JULIA_LIBDIRS, [joinpath(Sys.BINDIR, Base.LIBDIR, "julia"), joinpath(Sys.BINDIR, Base.LIBDIR)])
    end
    return JULIA_LIBDIRS
end


text:
include("stages.jl")
include("stop_conditions.jl")
include("hooks.jl")
include("run.jl")


text:
using Gadfly

set_default_plot_size(6inch, 6inch)

plot(x=rand(100), y=rand(100), size=rand(1:8, 100)./100, Geom.point)


text:
"""
    Fermi.CoupledCluster.RCCSD{T}(Alg::CTF)

Compute a RCCSD wave function using the Compiled time factorization algorithm (CTF)
"""
function ecRCCSD{T}(Alg::CTF) where T <: AbstractFloat
    @output "Calling CASCI module...\n"
    # Call CASCI
    cas = Fermi.ConfigurationInteraction.CASCI{T}()
    ecRCCSD{T}(cas, Alg)
end

function ecRCCSD{T}(cas::Fermi.ConfigurationInteraction.CASCI, Alg::CTF) where T <: AbstractFloat

    # Save reference wavefunction and process CAS data. Modify Ref (if not HF)
    @output "Processing CAS data...\n"
    refdet, Casdata = process_cas(cas)
    refwfn = cas.ref
    ecRCCSD{T}(refdet, refwfn, Casdata, Alg)
end

function ecRCCSD{T}(refdet::Determinant, refwfn::Fermi.HartreeFock.RHF, Casdata::Dict{String,Array}, Alg::CTF) where T <: AbstractFloat
    # Print intro
    Fermi.CoupledCluster.print_header()
    @output "\n    • Computing Externally Corrected CCSD with the ecRCCSD module.\n\n"

    # Recover integrals object
    ints = refwfn.ints

    # Get MO Integrals
    drop_occ = Fermi.CurrentOptions["drop_occ"]
    drop_vir = Fermi.CurrentOptions["drop_vir"]

    # Process CAS data to get ecT1 and ecT2 (Cluster Decomposition step)
    frozen = Fermi.CurrentOptions["cas_frozen"]
    active = Fermi.CurrentOptions["cas_active"] ≢ -1 ? Fermi.CurrentOptions["cas_active"] : refwfn.nvir+refwfn.ndocc-frozen

    if drop_occ > frozen
        error("\nFrozen orbitals in the CC step ($drop_occ) cannot be greater than the number of frozen orbitals in the CAS ($frozen).")
    end

    if drop_vir > (refwfn.ndocc+refwfn.nvir) - active - frozen
        error("\nToo many virtual orbitals dropped ($drop_vir) for the active space.")
    end

    @output "\n   • CAS Decomposition started:\n"
    t = @elapsed T1, T2, ecT1, ecT2 = cas_decomposition(refdet, Casdata, refwfn.ndocc, drop_occ, ints["FOV"], ints["OOVV"], ints["OVVV"], ints["OOOV"])
    @output "Finished in {:5.5} seconds.\n" t
    Casdata = nothing

    RCCSD{T}(refwfn, ints, T1, T2, Alg, ecT1=ecT1, ecT2=ecT2)
end

function process_cas(cas::Fermi.ConfigurationInteraction.CASCI)

    # This function process the CAS wave function to return a dictionary with C coefficients
    # organized by excitation level, along with the corresponding determinant
    # CI coefficients are intermediate normalized.
    # The reference determinant is taken as the HF one.

    dets = cas.dets
    Ccas = cas.coef

    ref = Determinant(repeat("1",cas.ref.ndocc),repeat("1",cas.ref.ndocc))
    z = 0
    acf = 0

    @output "   • CAS Composition\n"
    @output "Excitation      N of dets\n"
    while acf < length(dets)
        x = count(d->excitation_level(ref,d)==z, dets)
        @output "{}           {}\n" z x
        z += 1
        acf += x
    end

    ref = dets[1]
    zeroth = repeat('1', cas.ref.ndocc)*repeat('0', cas.ref.nvir)
    hf = Fermi.ConfigurationInteraction.Determinant(zeroth, zeroth)
    if ref == hf
        @output "Dominant configuration is the RHF determinant.\n"
    else
        error("Dominant determinant is not the RHF.")
    end
    C0 = Ccas[1]

    # Intermediate Normalization
    Ccas = Ccas ./ C0

    # Split the Cas data into excitation level
    Casdata = Dict{String,Array}()

    Casdata["C1"] = Float64[]
    Casdata["C2"] = Float64[]
    Casdata["C3"] = Float64[]
    Casdata["C4"] = Float64[]

    Casdata["C1dets"] = Determinant[]
    Casdata["C2dets"] = Determinant[]
    Casdata["C3dets"] = Determinant[]
    Casdata["C4dets"] = Determinant[]

    for i in eachindex(dets)

        exc = excitation_level(ref, dets[i])

        if exc == 1
            push!(Casdata["C1"], Ccas[i])
            push!(Casdata["C1dets"], dets[i])

        elseif exc == 2
            push!(Casdata["C2"], Ccas[i])
            push!(Casdata["C2dets"], dets[i])

        elseif  exc == 3
            push!(Casdata["C3"], Ccas[i])
            push!(Casdata["C3dets"], dets[i])

        elseif  exc == 4
            push!(Casdata["C4"], Ccas[i])
            push!(Casdata["C4dets"], dets[i])
        end
    end
    return ref, Casdata
end


text:
module TestPkg

greet() = print("Hello World!")
greeting() = print("Hello")

end # module


text:
module CLI

include("mcmc.jl")
include("sample.jl")
include("select.jl")

end # module

max_stars_repo_name: biobakery/GaPLAC

max_issues_repo_name: biobakery/GaPLAC

max_forks_repo_name: biobakery/GaPLAC


text:
module StochasticAiry

using LinearAlgebra
using Random, Distributions
using SpecialFunctions

import Distributions: ContinuousMatrixDistribution
import Random: rand

export SAi, rand

struct SAi <: ContinuousMatrixDistribution
    beta::Real
end

function rand(d::SAi, argtype::Symbol, basepoint::Number, argrange::AbstractRange{T}; kwargs...) where T <: Real

    if argtype == :time
        t_grid = collect(argrange)
        dt = step(argrange)

        U0 = rand(AiryU0(d.beta), argrange)

        return solve_sai(U0, t_grid, dt, basepoint)
    elseif argtype == :space

        return rand(d, :space, basepoint, collect(argrange); kwargs...)
    else
        throw("argtype must be :time or :space")
    end
end

function rand(d::SAi, argtype::Symbol, basepoint::Real, argrange::Array{T,dim};
              verbose=false, brownian_gridlen=3*10^3, brownian_gridmax=10) where {T <: Number, dim}

    if argtype != :space
        throw("Array argranges only defined for argtype == :space")
    end

    t = basepoint
    s = brownian_gridmax
    t_grid = range(s, t, length=brownian_gridlen);

    dt = (t - s)/brownian_gridlen

    U0 = rand(AiryU0(d.beta), t_grid)

    solve_phi(U0, collect(t_grid), dt, argrange; verbose=verbose)

end

include("AiryU0.jl")
include("Kernels.jl")
include("AirySolvers.jl")

end


text:
mutable struct MSetLPSolver <: MSetSolver
  solver
  m # Int
  n_arms # Int

  # Optimisation model (reused by solve_linear).
  model # ::Model
  x # ::Vector{Variable}

  function MSetLPSolver(solver)
    return new(solver, nothing, nothing, nothing, nothing)
  end
end

has_lp_formulation(::MSetLPSolver) = true
supports_solve_budgeted_linear(::MSetLPSolver) = true
supports_solve_all_budgeted_linear(::MSetLPSolver) = false
approximation_ratio(::MSetLPSolver) = 1.0
approximation_term(::MSetLPSolver) = 0.0
approximation_ratio_budgeted(::MSetLPSolver) = 1.0
approximation_term_budgeted(::MSetLPSolver) = 0.0

copy(solver::MSetLPSolver) = MSetLPSolver(solver.solver)

function build!(solver::MSetLPSolver, m::Int, n_arms::Int)
  solver.m = m
  solver.n_arms = n_arms
  solver.model = Model(solver.solver)
  solver.x = @variable(solver.model, [1:n_arms], binary=true, lower_bound=0.0)

  for i in 1:n_arms
    set_name(solver.x[i], "x_$i")
  end

  @constraint(solver.model, sum(solver.x) <= m)
end

function get_lp_formulation(solver::MSetLPSolver, rewards::Dict{T, Float64}) where T
  return solver.model,
    sum(rewards[i] * solver.x[i] for i in keys(rewards)),
    Dict{Int, JuMP.VariableRef}(i => solver.x[i] for i in keys(rewards))
end


text:
using Darwin
using RoboGrid
import Darwin: uniform_mutation
import RoboGrid: water_meta, cross_meta
import Base.isless
import JSON.json

include("cgp.jl")
include("darwin.jl")

function NeurodevoInd(genes::Array{Array{Float64}}, fitness::Array{Float64})
    NeurodevoInd(genes, fitness, RoboGrid.coverage_fitness, 1)
end

function generation(e::Evolution)
    fits = [RoboGrid.coverage_fitness, RoboGrid.water_memorize_fitness,
            RoboGrid.cross_memorize_fitness, RoboGrid.cross_strategy_fitness]
    for i in e.population
        i.seed = floor(Int64, (e.gen - 1) / 10)
        func = mod(floor(Int64, (e.gen - 1) / e.cfg["goal_gen"]), length(fits)) + 1
        i.func = fits[func]
    end
end

function robo_controller(m::Model, inputs::Array{Float64})
    Neurodevo.reward!(m, [inputs[1]])
    Neurodevo.step!(m, inputs[2:end])
end

function robo_eval(ind::Individual)
    cfg = cgp_cfg(Config("cfg/evo.yaml"), ind.genes[1])
    c = cgp_controller(cfg, ind.genes[2:end])
    m = Model(cfg, c)
    nin = 9
    nout = length(RoboGrid.ACTIONS)
    if cfg["init_method"] == 0
        random_init!(m, nin, nout; nreward=1, nhidden=cfg["nhidden"])
    else
        layered_init!(m, nin, nout; nreward=1, nhidden=cfg["nhidden"])
    end
    develop!(m)
    f(x::Array{Float64}) = robo_controller(m, x)
    [ind.func(f; seed=ind.seed)]
end


text:
"""
	traversal(potree::String, params::ParametersOrthophoto)

Trie traversal.
If entire point cloud falls in volume process all files of Potree project
else travers trie, depth search first, and process nodes falling in region of interest.

Input:
 - params: initial parameters
 - potree: potree hierarchy
"""
function traversal(potree::String, params::ParametersOrthophoto)
    flushprintln("= ")
    flushprintln("= PROJECT: $potree")
    flushprintln("= ")

    metadata = CloudMetadata(potree) # metadata of current potree project
    trie = potree2trie(potree)
    params.numNodes = length(keys(trie))

    # if model contains the whole point cloud ( == 2)
    #	process all files
    # else
    # 	navigate potree

    intersection =
        Common.modelsdetection(params.model, metadata.tightBoundingBox)

    if intersection == 2
        flushprintln("FULL model")
        for k in keys(trie)
            params.numFilesProcessed = params.numFilesProcessed + 1
            if params.numFilesProcessed % 100 == 0
                flushprintln(
                    params.numFilesProcessed,
                    " files processed of ",
                    params.numNodes,
                )
            end

            file = trie[k]
            updateWithoutControl!(params, file)

        end
    elseif intersection == 1
        flushprintln("DFS")
        dfs(trie, params)

        if params.numNodes - params.numFilesProcessed > 0
            flushprintln("$(params.numNodes-params.numFilesProcessed) file of $(params.numNodes) not processed - out of region of interest")
        end
    elseif intersection == 0
        flushprintln("OUT OF REGION OF INTEREST")
    end

end

"""
	dfs(t::DataStructures.Trie{String},
	params::Union{ParametersOrthophoto,ParametersExtraction},
	s::Union{Nothing,IOStream},n::Int64,nfiles::Int64,l::Int64)

Depth search first.
"""
function dfs(trie::DataStructures.Trie{String}, params::ParametersOrthophoto)# due callback: 1 con controllo e 1 senza controllo

    file = trie.value # path to node file
    nodebb = FileManager.las2aabb(file) # aabb of current octree
    inter = Common.modelsdetection(params.model, nodebb)

    if inter == 1
        # intersecato ma non contenuto
        # alcuni punti ricadono nel modello altri no
        params.numFilesProcessed = params.numFilesProcessed + 1
        if params.numFilesProcessed % 100 == 0
            flushprintln(
                params.numFilesProcessed,
                " files processed of ",
                params.numNodes,
            )
        end

        updateWithControl!(params, file) # update with check
        for key in collect(keys(trie.children)) # for all children
            dfs(trie.children[key], params)
        end
    elseif inter == 2
        # contenuto: tutti i punti del albero sono nel modello
        for k in keys(trie)
            params.numFilesProcessed = params.numFilesProcessed + 1
            if params.numFilesProcessed % 100 == 0
                flushprintln(
                    params.numFilesProcessed,
                    " files processed of ",
                    params.numNodes,
                )
            end
            file = trie[k]
            updateWithoutControl!(params, file) # update without check
        end
    end

end


text:
using Pkg
@info "Setting up ACRLab. This can take several minutes the first time."
Pkg.activate("ACRL")
# Pkg.add(url="https://github.com/janbruedigam/ACRLab.jl") # Needed when updating the manifest file
Pkg.instantiate()

@info "Loading ACRLab. This can take several minutes the first time."
using ACRLab

@info "Running main()"
main()


text:
# Defines the core Node abstraction
# Defines the core Node and operations on the graph

using Base
importall Base

typealias Float Float64
typealias TensorValue Union{Real, Array}

###############
# Basic types #
###############
abstract OpType

type Node
    op::OpType
    inputs::Vector{Node}
    outputs::Vector{Node}
    name::AbstractString
end

pred(node::Node) = node.inputs
succ(node::Node) = node.outputs

# TODO: Actual scoping on naming
function Node(data::OpType, name::AbstractString="")
    newName = length(name) == 0 ? "$(gensym())" : name
    Node(data, Vector{Node}(), Vector{Node}(), newName)
end

function Node(op::OpType, inputs::Vector{Node}, name::AbstractString="")
    newName = length(name) == 0 ? "$(gensym())" : name
    node = Node(op, inputs, Vector{Node}(), newName)
    for i in inputs
        push!(i.outputs, node)
    end
    node
end

immutable Graph
    nodes::Set{Node}
end

# May use to have a notion of "default" graph
# const Dict{Symbol, Graph} context
# context[:default] = Graph(Set{Node}())
########################
# Gradient computation #
########################

# out w.r.t. each element of wrt
function grad(out::Node, wrt::Vector{Node})
    # input gradients = derivative w.r.t. input - may be different for the nth and n+1th argument
    # output gradient = derivative w.r.t. output (i.e. before we apply this node's grad)

    # When we process a node, associate its input gradient with the appropriate input
    node_to_grad_vec = Dict{Node, Vector{Node}}(out => [fill_like(cg.constant(1.0), out)])

    # Map a node to its single output (usually sum of all output gradients)
    node_to_grad = Dict{Node, Node}()

    # Set of nodes on all paths between the set `wrt` and `out`
    for node = reverse(toposort_on_path([out], wrt))
        # Compute the gradient w.r.t. the node's output
        @assert haskey(node_to_grad_vec, node)
        grads = node_to_grad_vec[node]
        if (length(grads) == 1)
            output_grad = node_to_grad_vec[node][1]
        else
            # Sum up the gradients of all the outputs
            output_grad = foldr((sum, next) -> sum + next, grads)
        end
        if length(node.name) > 0
            output_grad.name = "Gin:$(node.name)"
        end
        node_to_grad[node] = output_grad

        if (length(node.inputs) > 0)
            # If the node has inputs, calculate the gradient w.r.t their outputs along this path
            input_grads = grad(node.op, node, output_grad, node.inputs...)
            for (input_node, grad) = zip(node.inputs, input_grads)
                if (!haskey(node_to_grad_vec, input_node))
                    node_to_grad_vec[input_node] = Vector{Node}()
                end
                if length(grad.name) > 0
                    output_grad.name = "Gout:$(input_node.name)"
                end
                push!(node_to_grad_vec[input_node], grad)
            end
        end
    end
    [node_to_grad[x] for x in wrt]
end

####################
# Graph operations #
####################

# Node hashing for use in deduping computation graphs
# Hash node itself, then recursively add in hashes of parents
function hash_node(node::Node)
    res = UInt64(node.op)
    for input = node.inputs
        res = hash(hash_node(input), res)
    end
    res
end

# Return set of nodes that are influenced in the DAG from any in set `nodes`
# i.e. that would be influenced by in a computation
# child=true means go to children in dag, false means go to parents
function influenced_by(nodes::Vector{Node}, child::Bool)
    queue = copy(nodes)
    influenced = Set{Node}(queue)
    next_method = child ? succ : pred
    while !isempty(queue)
        next = pop!(queue)
        for node = next_method(next)
            if node in influenced
                continue
            end
            push!(queue, node)
            push!(influenced, node)
        end
    end
    influenced
end

function nodes_on_path(sink::Vector{Node}, source::Vector{Node}=Vector{Node}())
    downstream = influenced_by(source, true)
    upstream = influenced_by(sink, false)
    if length(source) > 0
        intersect(upstream, downstream)
    else
        upstream
    end

end

function toposort_on_path(sink::Vector{Node}, source::Vector{Node}=Vector{Node}())
    nodes = nodes_on_path(sink, source)
    toposorted = toposort(get_graph(union(sink, source)))
    intersect(toposorted, nodes)
end

# Returns nodes in topological order
function toposort(graph::Graph)
    result = Vector{Node}()
    marks = Dict{Node, Symbol}() # :marked, :temp
    function visit(cur::Node)
        mark = get(marks, cur, :unmarked)
        if isequal(mark, :temp)
            @assert false && "Graph is not a DAG!"
            return
        elseif isequal(mark, :unmarked)
            marks[cur] = :temp
            for child = succ(cur)
                visit(child)
            end
            marks[cur] = :marked
            push!(result, cur)
        end
    end

    for node = graph.nodes
        if isequal(get(marks, node, :unmarked), :unmarked)
            visit(node)
        end
    end
    reverse!(result)
    result
end

# Return graph consisting of all nodes connected to given Variables
function get_connected(nodes::Set{Node})
    stack = collect(Node, nodes)
    seen = Set{Node}(nodes)
    while !isempty(stack)
        cur = pop!(stack)
        prev = pred(cur)
        next = succ(cur)
        # Would like the chain(.) function
        for n in prev
            if !(n in seen)
                push!(seen, n)
                push!(stack, n)
            end
        end
        for n in next
            if !(n in seen)
                push!(seen, n)
                push!(stack, n)
            end
        end
    end
    seen
end

function get_graph(nodes::Vector{Node})
    Graph(get_connected(Set(nodes)))
end

function group_nodes(nodes::Vector{Node}, group::AbstractString)
    for node in nodes
        node.name = string(group, "/", node.name)
    end
end

function group_between(inputs::Vector{Node}, outputs::Vector{Node}, group::AbstractString ; include_in=true, include_out=true)
    set = nodes_on_path(outputs, inputs)
    if !include_in
        set = setdiff(set, Set(inputs))
    end
    if !include_out
        set = setdiff(set, Set(outputs))
    end
    group_nodes(collect(set), group)
end


text:
module Rfam

using Downloads: download
using Scratch: @get_scratch!, clear_scratchspaces!
using FASTX: FASTA
import Gzip_jll

include("files.jl")
include("fetch.jl")
include("util.jl")

const rfam_version = "14.7"
const rfam_baseurl = "http://ftp.ebi.ac.uk/pub/databases/Rfam/$rfam_version"
cache = "" # directory to store downloaded files (set during __init__)

function __init__()
    # Directory to store downloaded files.
    # This space will be unique for each Rfam database version.
    global cache = @get_scratch!("Rfam-$rfam_version")
end

"""
    clear()

Removes all downloaded data.
"""
function clear()
    clear_scratchspaces!(Rfam)
end

end


text:
# Image resizing utilities

export anti_aliasing_filter

function anti_aliasing_filter(n::Integer; fact::Integer=1, flat::Bool=false, coeff::Real=1, T::DataType=Float32)
    kmax = T(pi)
    k = kmax*coord_norm(n)
    filt = zeros(T, n)
    if ~flat
        i1 = abs.(k) .<= kmax/2^fact-kmax/2^(fact+1)
        i2 = (abs.(k) .> kmax/2^fact-kmax/2^(fact+1)) .&& (abs.(k) .<= kmax/2^fact+kmax/2^(fact+1))
        filt[i2] .= ((cos.(range(-T(pi), T(pi); length=length(findall(i2)))).+1)/2).^T(coeff)
    else
        i1 = abs.(k) .<= kmax/2^fact
    end
    filt[i1] .= 1
    return filt
end

function anti_aliasing_filter(K::KSpaceCartesianSampling{T}; fact::Integer=1, flat::Bool=false, coeff::Real=1) where {T<:Real}
    nt, nk = size(K)
    perm = (K.phase_encoding..., K.readout)
    filt = reshape(permutedims(anti_aliasing_filter(K.X.n; fact=fact, flat=flat, coeff=coeff), perm), nt, nk)[K.subsampling,:]
    i_pe, i_r = downscale_phase_encode_index(K; fact=fact, readout=true)
    return filt[i_pe, i_r]
end

function anti_aliasing_filter(n::NTuple{3,Integer}; fact::Integer=1, flat::Bool=false, coeff::Real=1, T::DataType=Float32)
    filt_x = reshape(anti_aliasing_filter(n[1]; fact=fact, flat=flat, coeff=coeff, T=T), :, 1, 1)
    filt_y = reshape(anti_aliasing_filter(n[2]; fact=fact, flat=flat, coeff=coeff, T=T), 1, :, 1)
    filt_z = reshape(anti_aliasing_filter(n[3]; fact=fact, flat=flat, coeff=coeff, T=T), 1, 1, :)
    return filt_x.*filt_y.*filt_z
end

anti_aliasing_filter(X::RegularCartesianSpatialSampling{T}; fact::Integer=1, flat::Bool=false, coeff::Real=1) where {T<:Real} = anti_aliasing_filter(X.n; fact=fact, flat=flat, coeff=coeff, T=T)

function downscale(u::AbstractArray{CT,3}, X::RegularCartesianSpatialSampling{T}; fact::Integer=1, flat::Bool=false, coeff::Real=1) where {T<:Real,CT<:RealOrComplex{T}}
    (fact == 0) && (return u)
    x, y, z = coord(X)
    itp = interpolate((vec(x), vec(y), vec(z)), u, Gridded(Linear()))
    u_ = ifft(ifftshift(anti_aliasing_filter(X; fact=fact, flat=flat, coeff=coeff)).*fft(u))
    X_h = downscale(X; fact=fact)
    u_h = similar(u_, size(X_h))
    x_h, y_h, z_h = coord(X_h)
    @inbounds for i = 1:X_h.n[1], j = 1:X_h.n[2], k = 1:X_h.n[3]
        u_h[i, j, k] = itp(x_h[i], y_h[j], z_h[k])
    end
    return u_h
end

function upscale(u::AbstractArray{CT,3}, X::RegularCartesianSpatialSampling{T}; fact::Integer=1) where {T<:Real,CT<:RealOrComplex{T}}
    (fact == 0) && (return u)
    x, y, z = coord(X)
    itp = extrapolate(interpolate((vec(x), vec(y), vec(z)), u, Gridded(Linear())), T(0))
    X_h = upscale(X; fact=fact)
    u_h = similar(u, size(X_h))
    x_h, y_h, z_h = coord(X_h)
    @inbounds for i = 1:X_h.n[1], j = 1:X_h.n[2], k = 1:X_h.n[3]
        u_h[i, j, k] = itp(x_h[i], y_h[j], z_h[k])
    end
    return u_h
end

max_stars_repo_name: grizzuti/UtilitiesForMRI

max_issues_repo_name: grizzuti/UtilitiesForMRI

max_forks_repo_name: grizzuti/UtilitiesForMRI


text:
using Distributions
using DifferentialEquations
using DataFrames
using CSV
using PyPlot
using LinearAlgebra

include("LinearModels.jl")

"""
	function mcmc(data,cols,p,N,total::Int,sicrfn!)

Metropolis-Hastings MCMC of the SIR ODE model assuming a Poisson likelihood
for the cumulative number of cases and deaths, and a noninformative prior
data = DataFrame of
cols = [active cases, deaths, recovered cases]
p[1:end-1] = array of parameters used in the model
p[end] = number infected on day 0 of data (Jan 22, 2020)
pfit = array of parameters to be fit
N = population of data region
total = number of MCMC steps
output = (loglikelihoods, parameters, maximum loglikelihood, max loglikelihood parameters)
sicrfn! = ODE file to use

"""
function mcmc(data::DataFrame,cols::Array,p::Array,pfit::Array,N::Float64,total::Int,sicrfn!,temp::Float64 = 1.,sigmaguess::Float64 = .002)

	chi = ll(data,cols,p,N,sicrfn!)
	accept = 0
	pt = copy(p)
	chiml = chi
	pml = copy(p)

	chiout = Array{Float64,1}(undef,total)
	pout = Array{Float64,2}(undef,total,length(p))

	for step in 1:total
		sampleparams!(pt,p,pfit,sigmaguess)
		chit = ll(data,cols,pt,N,sicrfn!)
		if rand() < exp((chi - chit)/temp)
			chi = chit
			p = copy(pt)
			accept += 1
			if chi < chiml
				chiml = chi
				pml = copy(p)
			end
		end
		chiout[step] = chi
		pout[step,:] = p
	end
	println("maxll: ",chiml)
	return chiout, pout, accept/total, chiml, pml
end

# generate MCMC step guess from current parameter value r using logNormal distirbution for all
# parameters except start day with scale parameter manually adjusted so approximately a third of guesses are accepted
function sampleparams!(pt,p,pfit,sigmaguess)
	for i in pfit
		d = Distributions.LogNormal(log(p[i]),sigmaguess)
		pt[i] = rand(d)
	end
end

# loglikelihood of data given parameters at r using a Poisson likelihood function
function ll(data,cols,p,N,sicrfn!)
	firstday = findfirst(data[:,cols[1]] .> 0)
	lastday = length(data[:,cols[1]])
	tspan = (firstday,lastday)
	uin = uinitial(data,cols,firstday,N)
	prediction,days = model(p,cols,uin,float.(tspan),N,sicrfn!)
	chi = 0
	for set in eachindex(cols)
		for i in eachindex(days)
			d = Distributions.Poisson(max(prediction[i,set],eps(Float64)))
			chi -= loglikelihood(d,[max(data[days[i],cols[set]],0)])
		end
	end
	return chi
end

function null_ll(data,cols,p,N,sicrfn!)
	firstday = findfirst(data[:,cols[1]] .> 0)
	lastday = length(data[:,cols[1]])
	tspan = (firstday,lastday)
	uin = uinitial(data,cols,firstday,N)
	prediction,days = model(p,cols,uin,float.(tspan),N,sicrfn!)
	chi = 0
	ent = 0
	for set in eachindex(cols)
		for i in eachindex(days)
			d = Distributions.Poisson(max(prediction[i,set],eps(Float64)))
			psample = rand(d)
			chi -= loglikelihood(d,[psample])
			ent += entropy(d)
		end
	end
	return chi,ent
end

function uinitial(data,cols,firstday,N)
	uin = Array{Float64,1}(undef,length(cols))
	for i in eachindex(cols)
		uin[i] = data[firstday,cols[i]]/N
	end
	return uin
end

# produce arrays of cases and deaths
function model(pin,cols,uin,tspan,N,sicrfn!)
	sol = modelsol(pin,uin,tspan,N,sicrfn!)
	prediction = Array{Float64,2}(undef,length(sol[1,:]),length(cols))
	prediction[:,1] = pin[2]*sol[1,:]*N
	for i = 2:length(cols)
		prediction[:,i] = pin[i+1]*sol[2,:]*N
	end
	return prediction, round.(Int,sol.t)
end

function modelsol(pin,uin,tspan,N,sicrfn!)
	p = pin[1:end-1]
	# I0 = (p[1]- p[3])/p[2] * uin[1] + pin[end]/N
	u0 = [pin[end]/N,uin[1]]
	prob = ODEProblem(sicrfn!,u0,tspan,p)
	solve(prob,saveat=1.)
end

function makehistograms(out)
	for i in 1:size(out)[2]
		subplot(2,3,i)
		hist(out[:,i],bins=100,density=true)
		title("param$(i)",fontsize=15)
	end
	tight_layout()
end

# returns predictions for C,R,D, matching data, and array of days
# solution array is in same row col order as data
# data columns have been ordered to match solutions.
function modelprediction(data,cols,p,N,sicrfn!)
	firstday = findfirst(data[:,cols[1]] .>0)
	lastday = length(data[:,cols[1]])
	sol = modelprediction(data,cols,p,N,sicrfn!,lastday)
	prediction = Array{Float64,2}(undef,length(sol[1,:]),length(cols))
	prediction[:,1] = p[2]*sol[1,:]*N
	for i = 2:length(cols)
		prediction[:,i] = p[i+1]*sol[2,:]*N
	end
	return prediction, data[firstday:end,cols], Int.(sol.t)
end

# returns solutions to model from day of first case to requested last day
function modelprediction(data,cols,p,N,sicrfn!,lastday)
	firstday = findfirst(data[:,cols[1]] .>0)
	tspan = (firstday,lastday)
	uin = uinitial(data,cols,firstday,N)
	modelsol(p,float.(uin),float.(tspan),N,sicrfn!)
end

# returns statistics for posteriors and other measures
function measures(psamples,pml)
	nparams = size(psamples,2)
	pout = Array{Tuple,1}(undef,nparams)
	for i in 1:nparams
		p = psamples[:,i]
		pout[i] = pml[i],mean(p),std(p),quantile(p,[.025,.5,.975])
	end

	R0 = r0(pml), mean(r0(psamples)),std(r0(psamples)),quantile(r0(psamples),[.025,.5,.975])
	return pout, R0
end

function r0(p::Array{Float64,1})
	beta,sigmac,sigmar,sigmad,q,f = p[1:6]
	sigma = sigmar + sigmad
	return (beta * (sigma + q * sigmac)) / (sigma * (sigmac + f * sigma))
end

function r0(p::Array{Float64,2})
	beta = p[:,1]
	sigmac = p[:,2]
	sigmar = p[:,3]
	sigmad = p[:,4]
	q = p[:,5]
	f = p[:,6]
	sigma = sigmar .+ sigmad
	return (beta .* (sigma .+ q .* sigmac)) ./ (sigma .* (sigmac .+ f .* sigma))
end

max_stars_repo_name: cannin/covid-sicr-test

max_issues_repo_name: cannin/covid-sicr-test

max_forks_repo_name: cannin/covid-sicr-test


text:
module FlexLinearAlgebra
using LinearAlgebra

import Base: (+), (-), (*), (==), sum, #(.*),
    getindex, setindex!, hash, show, keys, values, size,
    keytype, valtype, length, haskey, hash, Vector, Matrix

#    dot, adjoint

#import LinearAlgebra: dot, adjoint

include("FlexVector.jl")
include("FlexMatrix.jl")
end # module