text:
# This file defines the MPcondition type, used to define the extra conditions on
#   transitions.

export MPcondition
"""
The `MPcondition` type defines a condition on a single attribute of a personnel member, to be used for checking if a transition will take place.

The type contains the following fields:
* `attribute::String`: the attribute upon which the prerequisite is defined.
* `value::Union{Real, String, Vector{String}}`: the value with which the attribute is compared.
* `operator::Function`: the relation between the attribute and the value. Valid relations are `==`, `!=`, `<`, `>`, `<=`, `>=`, `∈`, and `∉`. Naturally, the operator must make sense for the type of the value. This means that each type has the following restrictions:
1. `Real`: the operator must be `==`, `!=`, `<`, `>`, `<=`, or `>=`
2. `String`: the operator must be `==` or `!=`
3. `Vector{String}`: the operator must be `∈` or `∉`.

A condition is always of the form `attribute operator value`. For example, if a personnel member requires a physical fitness score of 14.0, the condition is written as "fitness" >= 14.0.

Constructor:
```
MPcondition(
    attribute::String,
    operator::Function,
    value::Union{Real, String, Vector{String}} )
```
This constructor creates a `MPcondition` object of the form `attribute` `operator` `value`.
"""
struct MPcondition

    attribute::String
    value::Union{Real, String, Vector{String}}
    operator::Function

    # Constructor
    function MPcondition( attribute::String, operator::Function,
        value::Union{Real, String, Vector{String}} )::MPcondition

        if ( value isa Real ) && ( operator ∉ [==, !=, <, <=, >, >=] )
            error( "A numeric value permits only ==, !=, <, <=, >, and >= as operators." )
        elseif ( value isa String ) && ( operator ∉ [==, !=] )
            error( "A string value permits only == and != as operators." )
        elseif ( value isa Vector{String} ) && ( operator ∉ [∈, ∉] )
            error( "A vector of string values permits only ∈ and ∉ as operators." )
        end  # if ( value isa Real ) && ...

        return new( attribute, value, operator )

    end  # MPcondition( attribute, operator, value )

end  # struct MPcondition


text:
using Test, Astro

# equinox tests from meeus chapter

print("27.a Approximate solstice")
eqa = equinox_approx(1962, "summer")
y, m , d = jd_to_cal(eqa)

@test y == 1962
@test m == 6
@test isapprox(d, 21, atol=1)

println(" passed")

print("27.a Finding nearest Equinox for 12 January 2007")
jd = cal_to_jd(2007,1,12)
eq = equinox(jd, "Spring", 1) # return jd of equinox to the nearest day

y, m , d = jd_to_cal(eq)
@test y == 2007
@test m == 3
@test isapprox(d, 21, atol=1)

println(" passed")

ext: jl


text:
module Diffusion

import MPI
using ImplicitGlobalGrid
using CUDA
using Plots

@views d_xa(A) = A[2:end  , :     ] .- A[1:end-1, :     ];
@views d_xi(A) = A[2:end  ,2:end-1] .- A[1:end-1,2:end-1];
@views d_ya(A) = A[ :     ,2:end  ] .- A[ :     ,1:end-1];
@views d_yi(A) = A[2:end-1,2:end  ] .- A[2:end-1,1:end-1];
@views  inn(A) = A[2:end-1,2:end-1]

@views function diffusion2D(nx, ny, nt, lx, ly, lam, cp_min, amp_T, Cp, comm)
    # Numerics
    me, dims = init_global_grid(nx, ny, 1; init_MPI=false, comm=comm, reorder=0)  # Initialize the implicit global grid
    dx       = lx/(nx_g()-1)                                                      # Space step in dimension x
    dy       = ly/(ny_g()-1)                                                      # ...        in dimension y

    # Array initializations
    T     = CUDA.zeros(Float64, nx,   ny, )
    dTedt = CUDA.zeros(Float64, nx-2, ny-2)
    qx    = CUDA.zeros(Float64, nx-1, ny-2)
    qy    = CUDA.zeros(Float64, nx-2, ny-1)

    # Initialize temperature with a Gaussian anomaly of amplitude amp_T
    T    .= CuArray([amp_T*exp(-((x_g(ix,dx,T)-lx/2)/2)^2-((y_g(iy,dy,T)-ly/2)/2)^2) for ix=1:size(T,1), iy=1:size(T,2)])

    # Preparation of visualisation
    gr()
    ENV["GKSwstype"]="nul"
    anim     = Animation();
    nx_v     = (nx-2)*dims[1]
    ny_v     = (ny-2)*dims[2]
    T_v      = zeros(nx_v, ny_v)
    T_nohalo = zeros(nx-2, ny-2)

    # Time loop
    nframes = 50
    dt = min(dx,dy)^2*cp_min/lam/4.1                            # Time step for the 2D Heat diffusion
    for it = 1:nt
        qx    .= -lam.*d_xi(T)./dx                              # Fourier's law of heat conduction: q_x   = -λ δT/δx
        qy    .= -lam.*d_yi(T)./dy                              # ...                               q_y   = -λ δT/δy
        dTedt .= 1.0./inn(Cp).*(-d_xa(qx)./dx .- d_ya(qy)./dy)  # Conservation of energy:           δT/δt = 1/cₚ (-δq_x/δx - δq_y/dy)
        T[2:end-1,2:end-1] .= inn(T) .+ dt.*dTedt               # Update of temperature             T_new = T_old + δT/δt
        update_halo!(T)                                         # Update the halo of T
        if mod(it, nt/nframes) == 1
            T_nohalo .= T[2:end-1,2:end-1]                      # Copy data to CPU removing the halo.
            gather!(T_nohalo, T_v)                              # Gather the global temperature field on process 0.
            if (me==0)                                          # Visualize it on process 0.
                println("Time step $it: maximum temperature is $(maximum(T_v))")
                heatmap(T_v, title="Temperature at time step $it", aspect_ratio=1, xlims=(1,nx_v), ylims=(1,ny_v), clims=(0,amp_T)); frame(anim)
            end
        end
    end

    # Postprocessing
    if (me==0) gif(anim, "diffusion2D.gif", fps = 15) end       # Create an animated gif on process 0.
    finalize_global_grid(;finalize_MPI=false)
end

Base.@ccallable function julia_diffusion(nx::Csize_t, ny::Csize_t, nt::Csize_t, lx::Cdouble, ly::Cdouble, lam::Cdouble, cp_min::Cdouble, amp_T::Cdouble, Cp_c::CuPtr{Cdouble}, comm_c::MPI.MPI_Comm)::Cint
    try
        comm = MPI.Comm_dup(MPI.Comm(comm_c))
        Cp   = unsafe_wrap(CuArray, Cp_c, (Int64(nx), Int64(ny)))
        diffusion2D(nx, ny, nt, lx, ly, lam, cp_min, amp_T, Cp, comm)
    catch
        Base.invokelatest(Base.display_error, Base.catch_stack())
        return 1
    end
    return 0
end

end # module

max_stars_repo_name: omlins/libdiffusion

max_issues_repo_name: omlins/libdiffusion

max_forks_repo_name: omlins/libdiffusion


text:
@testset "Get line" begin
    table = init_log_table(
        (id=:open_nodes, name="#Open"),
        (id=:closed_nodes, name="#Closed", alignment=:right),
        (id=:time, name="#Time", alignment=:left),
    )
    set_value!(table, :open_nodes, 10)
    set_value!(table, :closed_nodes, 20)
    set_value!(table, :time, 1.1)
    line = TableLogger.get_line(table)
    @test occursin("10", line)
    @test occursin("20", line)
    @test occursin("1.1", line)
    @test occursin("20  1.1", line)
    @test length(line) == sum(table.setup.widths)
end

@testset "Shall print line" begin
    table = init_log_table(
        (id=:open_nodes, name="#Open"),
        (id=:closed_nodes, name="#Closed"),
    )
    set_value!(table, :open_nodes, 10)
    set_value!(table, :closed_nodes, 20)
    @test TableLogger.shall_print_line(table; force=true)
    @test TableLogger.shall_print_line(table)

    TableLogger.update_for_new_row(table)
    TableLogger.fill_from_prev!(table)

    @test !TableLogger.shall_print_line(table)
    @test TableLogger.shall_print_line(table; force=true)
end

@testset "Diff5 shall print line" begin
    table = init_log_table(
        (id=:open_nodes, name="#Open"),
    )
    set_value!(table, :open_nodes, Diff5(10))
    @test TableLogger.shall_print_line(table; force=true)
    @test TableLogger.shall_print_line(table)

    TableLogger.update_for_new_row(table)
    set_value!(table, :open_nodes, Diff5(12))
    @test !TableLogger.shall_print_line(table)
    set_value!(table, :open_nodes, Diff5(15))
    @test TableLogger.shall_print_line(table)
end

@testset "DiffX shall print line" begin
    table = init_log_table(
        (id=:open_nodes, name="#Open"),
    )
    Diff7 = DiffX(7)
    set_value!(table, :open_nodes, Diff7(10))
    @test TableLogger.shall_print_line(table; force=true)
    @test TableLogger.shall_print_line(table)

    TableLogger.update_for_new_row(table)
    set_value!(table, :open_nodes, Diff7(12))
    @test !TableLogger.shall_print_line(table)
    set_value!(table, :open_nodes, Diff7(15))
    @test !TableLogger.shall_print_line(table)
    set_value!(table, :open_nodes, Diff7(17))
    @test TableLogger.shall_print_line(table)
end

@testset "Diffd1 shall print line" begin
    table = init_log_table(
        (id=:open_nodes, name="#Open"),
    )
    Diffd1 = DiffX(0.1)
    set_value!(table, :open_nodes, Diffd1(10.0))
    @test TableLogger.shall_print_line(table; force=true)
    @test TableLogger.shall_print_line(table)

    TableLogger.update_for_new_row(table)
    set_value!(table, :open_nodes, Diffd1(10.02))
    @test !TableLogger.shall_print_line(table)
    set_value!(table, :open_nodes, Diffd1(10.09))
    @test !TableLogger.shall_print_line(table)
    set_value!(table, :open_nodes, Diffd1(10.101))
    @test TableLogger.shall_print_line(table)
end


text:
#=
The code below has been ported to Julia from the original C source:
https://github.com/nexusformat/HDF5-External-Filter-Plugins/blob/master/BZIP2/src/H5Zbzip2.c
The filter function  H5Z_filter_bzip2 was adopted from:
PyTables http://www.pytables.org.
The plugin can be used with the HDF5 library version 1.8.11+ to read HDF5 datasets compressed with bzip2 created by PyTables.
License: licenses/H5Zbzip2_LICENSE.txt

The following license applies to the Julia port.
Copyright (c) 2021 Mark Kittisopikul and Howard Hughes Medical Institute. License MIT, see LICENSE.txt
=#
module H5Zbzip2

using CodecBzip2
import CodecBzip2: libbzip2
using HDF5.API
import HDF5.Filters: Filter, filterid, register_filter, filtername, filter_func, filter_cfunc

export H5Z_FILTER_BZIP2, H5Z_filter_bzip2, Bzip2Filter

const H5Z_FILTER_BZIP2 = API.H5Z_filter_t(307)
const bzip2_name = "HDF5 bzip2 filter; see http://www.hdfgroup.org/services/contributions.html"

function H5Z_filter_bzip2(flags::Cuint, cd_nelmts::Csize_t,
                        cd_values::Ptr{Cuint}, nbytes::Csize_t,
                        buf_size::Ptr{Csize_t}, buf::Ptr{Ptr{Cvoid}})::Csize_t
    outbuf = C_NULL
    outdatalen = Cuint(0)

    # Prepare the output buffer

    try

    if flags & API.H5Z_FLAG_REVERSE != 0
        # Decompress

        outbuflen = nbytes * 3 + 1
        outbuf = Libc.malloc(outbuflen)
        if outbuf == C_NULL
            error("H5Zbzip2: memory allocation failed for bzip2 decompression.")
        end

        stream = CodecBzip2.BZStream()
        # Just use default malloc and free
        stream.bzalloc = C_NULL
        stream.bzfree = C_NULL
        # BZ2_bzDecompressInit
        ret = CodecBzip2.decompress_init!(stream, 0, false)
        if ret != CodecBzip2.BZ_OK
            errror("H5Zbzip2: bzip2 decompress start failed with error $ret.")
        end

        stream.next_out = outbuf
        stream.avail_out = outbuflen
        stream.next_in = unsafe_load(buf)
        stream.avail_in = nbytes

        cont = true

        while cont
            # BZ2_bzDecompress
            ret = CodecBzip2.decompress!(stream)
            if ret < 0
                error("H5Zbzip2: bzip2 decompression failed with error $ret.")
            end
            cont = ret != CodecBzip2.BZ_STREAM_END
            if cont && stream.avail_out == 0
                # Grow the output buffer
                newbuflen = outbuflen * 2
                newbuf = Libc.realloc(outbuf, newbuflen)
                if newbuf == C_NULL
                    error("H5Zbzip2: memory allocation failed for bzip2 decompression.")
                end
                stream.next_out = newbuf + outbuflen
                stream.avail_out = outbuflen
                outbuf = newbuf
                outbuflen = newbuflen
            end
        end

        outdatalen = stream.total_out_lo32
        # BZ2_bzDecompressEnd
        ret = CodecBzip2.decompress_end!(stream)
        if ret != CodecBzip2.BZ_OK
            error("H5Zbzip2: bzip2 compression end failed with error $ret.")
        end
    else
        # Compress data

        # Maybe not the same size as outdatalen
        odatalen = Cuint(0)
        blockSize100k = 9

        # Get compression blocksize if present
        if cd_nelmts > 0
            blockSize100k = unsafe_load(cd_values)
            if blockSize100k < 1 || blockSize100k > 9
                error("H5Zbzip2: Invalid compression blocksize: $blockSize100k")
            end
        end

        # Prepare the output buffer
        outbuflen = nbytes + nbytes ÷ 100 + 600 # worse case (bzip2 docs)
        outbuf = Libc.malloc(outbuflen)
        @debug "Allocated" outbuflen outbuf
        if outbuf == C_NULL
            error("H5Zbzip2: Memory allocation failed for bzip2 compression")
        end

        # Compress data
        odatalen = outbuflen
        r_odatalen = Ref{Cuint}(odatalen)
        ret = BZ2_bzBuffToBuffCompress(outbuf, r_odatalen, unsafe_load(buf), nbytes,
                                       blockSize100k, 0, 0)
        outdatalen = r_odatalen[]
        if ret != CodecBzip2.BZ_OK
            error("H5Zbzip2: bzip2 compression failed with error $ret.")
        end
    end # if flags & API.H5Z_FLAG_REVERSE != 0
    Libc.free(unsafe_load(buf))
    unsafe_store!(buf, outbuf)
    unsafe_store!(buf_size, outbuflen)

    catch err
        #  "In the case of failure, the return value is 0 (zero) and all pointer arguments are left unchanged."
        outdatalen = Csize_t(0)
        if outbuf != C_NULL
            Libc.free(outbuf)
        end
        @error "H5Zbzip2.jl Non-Fatal ERROR: " err
        display(stacktrace(catch_backtrace()))
    end # try - catch

    return Csize_t(outdatalen)
end # function H5Z_filter_bzip2

# Need stdcall for 32-bit Windows?
function BZ2_bzBuffToBuffCompress(dest, destLen, source, sourceLen, blockSize100k, verbosity, workFactor)
    @static if CodecBzip2.WIN32
        return ccall(
            ("BZ2_bzBuffToBuffCompress@28", libbzip2),
            stdcall,
            Cint,
            (Ptr{Cchar}, Ptr{Cuint}, Ptr{Cchar}, Cuint, Cint, Cint, Cint),
            dest,
            destLen,
            source,
            sourceLen,
            blockSize100k,
            verbosity,
            workFactor
        )
    else
        return ccall(
            (:BZ2_bzBuffToBuffCompress, libbzip2),
            Cint,
            (Ptr{Cchar}, Ptr{Cuint}, Ptr{Cchar}, Cuint, Cint, Cint, Cint),
            dest,
            destLen,
            source,
            sourceLen,
            blockSize100k,
            verbosity,
            workFactor
        )
    end
end

function BZ2_bzBuffToBuffDecompress(dest, destLen, source, sourceLen, small, verbosity)
    @static if CodecBzip2.WIN32
        return ccall(
            ("BZ2_bzBuffToBuffDecompress@24", libbzip2),
            stdcall,
            Cint,
            (Ptr{Cchar}, Ptr{Cuint}, Ptr{Cchar}, Cuint, Cint, Cint),
            dest,
            destLen,
            source,
            sourceLen,
            small,
            verbosity
        )
    else
        return ccall(
            (:BZ2_bzBuffToBuffDecompress, libbzip2),
            Cint,
            (Ptr{Cchar}, Ptr{Cuint}, Ptr{Cchar}, Cuint, Cint, Cint),
            dest,
            destLen,
            source,
            sourceLen,
            small,
            verbosity
        )
    end
end

# Filters Module

struct Bzip2Filter <: Filter
    blockSize100k::Cuint
end
Bzip2Filter() = Bzip2Filter(9)

filterid(::Type{Bzip2Filter}) = H5Z_FILTER_BZIP2
filtername(::Type{Bzip2Filter}) = bzip2_name
filter_func(::Type{Bzip2Filter}) = H5Z_filter_bzip2
filter_cfunc(::Type{Bzip2Filter}) = @cfunction(H5Z_filter_bzip2, Csize_t,
                                               (Cuint, Csize_t, Ptr{Cuint}, Csize_t,
                                               Ptr{Csize_t}, Ptr{Ptr{Cvoid}}))

function __init__()
    register_filter(Bzip2Filter)
end

end # module H5Zbzip2


text:
module WooldridgeDatasets
    using CSV: lazystrings

    include("data.jl")
export wooldridge
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:2] <= 1)

                     @variable(m, -3.357520876571253 <= q <= 3.7351340364874854)

                     add_NL_constraint(m, :(sigmoid(0.618782057809276 + -0.4828154761503023*sigmoid(0.8842828756523886 + -0.6089179668773923*sigmoid(0.019018219753813348 + 0.8244023518746615*sigmoid(-0.7862452642307147 + -0.6499989200319587*$(x[1]) + -0.39995929749666725*$(x[2])) + 0.3296833580573324*sigmoid(0.2947581350949262 + -0.3447038526107806*$(x[1]) + -0.5996590697259889*$(x[2]))) + 0.44138625638024864*sigmoid(-0.9126485111761089 + 0.9176394232569307*sigmoid(-0.7862452642307147 + -0.6499989200319587*$(x[1]) + -0.39995929749666725*$(x[2])) + 0.371990627187599*sigmoid(0.2947581350949262 + -0.3447038526107806*$(x[1]) + -0.5996590697259889*$(x[2])))) + 0.3197591788003944*sigmoid(0.7251254505916935 + -0.923704860259535*sigmoid(0.019018219753813348 + 0.8244023518746615*sigmoid(-0.7862452642307147 + -0.6499989200319587*$(x[1]) + -0.39995929749666725*$(x[2])) + 0.3296833580573324*sigmoid(0.2947581350949262 + -0.3447038526107806*$(x[1]) + -0.5996590697259889*$(x[2]))) + 0.28767280168552967*sigmoid(-0.9126485111761089 + 0.9176394232569307*sigmoid(-0.7862452642307147 + -0.6499989200319587*$(x[1]) + -0.39995929749666725*$(x[2])) + 0.371990627187599*sigmoid(0.2947581350949262 + -0.3447038526107806*$(x[1]) + -0.5996590697259889*$(x[2]))))) + sigmoid(-0.0992210488643468 + 0.772244110024173*sigmoid(0.8842828756523886 + -0.6089179668773923*sigmoid(0.019018219753813348 + 0.8244023518746615*sigmoid(-0.7862452642307147 + -0.6499989200319587*$(x[1]) + -0.39995929749666725*$(x[2])) + 0.3296833580573324*sigmoid(0.2947581350949262 + -0.3447038526107806*$(x[1]) + -0.5996590697259889*$(x[2]))) + 0.44138625638024864*sigmoid(-0.9126485111761089 + 0.9176394232569307*sigmoid(-0.7862452642307147 + -0.6499989200319587*$(x[1]) + -0.39995929749666725*$(x[2])) + 0.371990627187599*sigmoid(0.2947581350949262 + -0.3447038526107806*$(x[1]) + -0.5996590697259889*$(x[2])))) + 0.07492471407403789*sigmoid(0.7251254505916935 + -0.923704860259535*sigmoid(0.019018219753813348 + 0.8244023518746615*sigmoid(-0.7862452642307147 + -0.6499989200319587*$(x[1]) + -0.39995929749666725*$(x[2])) + 0.3296833580573324*sigmoid(0.2947581350949262 + -0.3447038526107806*$(x[1]) + -0.5996590697259889*$(x[2]))) + 0.28767280168552967*sigmoid(-0.9126485111761089 + 0.9176394232569307*sigmoid(-0.7862452642307147 + -0.6499989200319587*$(x[1]) + -0.39995929749666725*$(x[2])) + 0.371990627187599*sigmoid(0.2947581350949262 + -0.3447038526107806*$(x[1]) + -0.5996590697259889*$(x[2]))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
## resolve the ground truth trajectory from ordinary differential equations

# ddy = g
# dx/dt = Fx + Gu
# [ddy;   = [0  0][dy  + [1][g=-10]
#   dy ]    [1  0]  y]   [0]
function falling(du,u,p,t)
  du[1] = 0    # drag
  du[2] = -10  # gravity
  du[3] = u[1]
  du[4] = u[2]
  nothing
end

u0 = [5,15,-20.0,0.0]
tspan = (0.0,5.0)
prob = ODEProblem(falling,u0,tspan)
sol = solve(prob)
# sol(t)

Plots.plot(sol, vars=(3,4))

## extract the desired measurement information from the physical simulation

tstop = 5.0
nposes = 10
ts = range(0, stop=tstop, length=nposes)
XY = zeros(2,nposes)

for i in 1:nposes
  t = ts[i]
  XY[:,i] = sol(t)[3:4]
end

gx = (t)->sol(t)[3]
gy = (t)->sol(t)[4]

# Calculate the true radar range and angle

ranges = sqrt.(vec(sum(XY.^2, dims=1)))
angles = atan.(XY[2,:], XY[1,:])

XYdense = [gx.(0:0.01:tstop) gy.(0:0.01:tstop)]';
randense = sqrt.(vec(sum(XYdense.^2, dims=1)))
angdense = atan.(XYdense[2,:], XYdense[1,:])


text:
function largest_palindrome()
    large = 0
    for i in 990:-11:100

        if i*999 < large
            return large
        end

        for j in 999:-1:100
            num = i*j

            if num < large
                break
            end

            if string(num) == reverse(string(num))
                large = num
                break
            end

        end
    end
end

println(largest_palindrome())

max_stars_repo_name: HarrisonGreen/Project-Euler-Solutions

max_issues_repo_name: HarrisonGreen/Project-Euler-Solutions

max_forks_repo_name: HarrisonGreen/Project-Euler-Solutions


text:
response = Dict(
        :question => "Sure.",
        :yelling => "Whoa, chill out!",
        :silence => "Fine. Be that way!",
        :misc => "Whatever.",
        :forceful_question => "Calm down, I know what I'm doing!",
)

struct Stimulus
    last_char
    all_caps::Bool
    empty::Bool
    function Stimulus(text)
        striped_text = strip(text)
        if isempty(striped_text)
            new(striped_text, false, true)
        else
            last_char = striped_text[end]
            uppercase_test = [isuppercase(char) for char in filter(isletter, text)]
            all_caps = all(uppercase_test) & any(uppercase_test)
            new(last_char, all_caps, false)
        end
    end
end

function bob(stimulus)
    stimulus = Stimulus(stimulus)
    if stimulus.last_char == '?'
        if stimulus.all_caps == true
            key = :forceful_question
        else
            key = :question
        end
    elseif stimulus.all_caps == true
        key = :yelling
    elseif stimulus.empty == true
        key = :silence
    else
        key = :misc
    end
    return response[key]
end

max_stars_repo_name: tfoliva/exercism

max_issues_repo_name: tfoliva/exercism

max_forks_repo_name: tfoliva/exercism


text:
using YaoArrayRegister, YaoBlocks, YaoBlocks.ConstGate, Random
using BenchmarkTools, PkgBenchmark

function bench(block)
    N = nqubits(block)
    r = rand_state(N)
    return @benchmarkable apply!($r, $block)
end

const SUITE = BenchmarkGroup()

SUITE["primitive"] = BenchmarkGroup()
SUITE["composite"] = BenchmarkGroup()

for G in [:X, :Y, :Z, :S, :T, :Sdag, :Tdag]
    GT = Expr(:(.), :ConstGate, QuoteNode(Symbol(G, :Gate)))
    GN = Expr(:(.), :ConstGate, QuoteNode(G))

    SUITE["primitive"][string(G)] = BenchmarkGroup()
    for n in 1:4:25
        SUITE["primitive"][string(G)][n] = @eval bench(repeat($n, $GN))
    end
end

for n in 5:5:25
    SUITE["composite"]["kron(rand_const)"] = bench(kron(rand([X, Y, Z, H]) for _ in 1:n))
    SUITE["composite"]["kron(sparse_const)"] =
        bench(kron(n, k => rand([X, Y, Z, H]) for k in randperm(n)[1:n÷5]))
    SUITE["composite"]["swap"] = bench(swap(n, 2, 4))
    SUITE["composite"]["swap"] = bench(pswap(n, 2, 4, 0.5))
end

function heisenberg(n::Int; periodic::Bool = true)
    Sx(i) = put(n, i => X)
    Sy(i) = put(n, i => Y)
    Sz(i) = put(n, i => Z)

    return sum(1:(periodic ? n : n - 1)) do i
        j = mod1(i, n)
        Sx(i) * Sx(j) + Sy(i) * Sy(j) + Sz(i) * Sz(j)
    end
end

for n in 5:5:25
    SUITE["primitive"]["TimeEvolution(heisenberg(n), 0.2)"] = bench(TimeEvolution(heisenberg(n), 0.2))
end


text:
module TestCasing
using SymSpellChecker: get_opcodes, transfer_casing_for_similar_text, transfer_casing_for_matching_text

include("preamble.jl")

@testset "opcode" begin
    @test get_opcodes("a", "") == [("delete", 1, 1, 1, 0)]
    @test get_opcodes("", "a") == [("insert", 1, 0, 1, 1)]
    @test get_opcodes("a", "b") == [("replace", 1, 1, 1, 1)]
    @test get_opcodes("a", "bc") == [("replace", 1, 1, 1, 2)]
    @test get_opcodes("a", "ab") == [("equal", 1, 1, 1, 1), ("insert", 2, 1, 2, 2)]
    @test get_opcodes("a", "bac") == [("insert", 1, 0, 1, 1), ("equal", 1, 1, 2, 2), ("insert", 2, 1, 3, 3)]
    @test get_opcodes("ac", "abc") == [("equal", 1, 1, 1, 1), ("insert", 2, 1, 2, 2), ("equal", 2, 2, 3, 3)]
    @test get_opcodes("qabxcd", "abycdf") == [("delete", 1, 1, 1, 0), ("equal", 2, 3, 1, 2),
                                                        ("replace", 4, 4, 3, 3), ("equal", 5, 6, 4, 5),
                                                        ("insert", 7, 6, 6, 6)]
end

@testset "transfer_casing" begin
    text_w_casing  = "Haw is the eeather in New York?"
    text_wo_casing = "how is the weather in new york?"
    @test transfer_casing_for_matching_text(text_w_casing, text_wo_casing) == "How is the weather in New York?"

    text_w_casing = "Haaw is the weeather in New York?"
    text_wo_casing = "how is the weather in new york?"
    @test transfer_casing_for_similar_text(text_w_casing, text_wo_casing) == "How is the weather in New York?"

    text_w_casing  = "BcDefg J kLn"
    text_wo_casing = "abcdefg hj kmn"
    @test transfer_casing_for_similar_text(text_w_casing, text_wo_casing) == "ABcDefg HJ kMn"
end
end # module


text:
import StatsBase as SB
using BenchmarkTools
using Normalization
using Statistics
using Test

@testset "1D normalization" begin
    # * 1D array
    _X = rand(100)
    X = copy(_X)
    T = fit(ZScore, X)
    Y = normalize(X, T)
    @test !isnothing(T.p)
    @test length(T.p) == 2
    @test length(T.p[1]) == 1 == length(T.p[2])
    @test Y ≈ (X.-mean(X))./std(X)
    @test denormalize(Y, T) ≈ X
    @test_nowarn normalize!(X, T)
    @test X == Y
    @test_nowarn denormalize!(Y, T)
    @test Y ≈ _X
end

@testset "2D normalization" begin # Adapted from https://github.com/JuliaStats/StatsBase.jl/blob/e8ab26500d9a34ef358b2d3cf619ae41c71785fc/test/transformations.jl

    #* 2D array
    _X = rand(10, 5)
    X = copy(_X)

    # * ZScore a 2D array over the first dim.
    T = fit(ZScore, X, dims=1)
    Y = normalize(X, T)
    @test !isnothing(T.p)
    @test length(T.p) == 2
    @test length(T.p[1]) == length(T.p[2]) == size(X, 2)
    @test Y ≈ (X.-mean(X, dims=1))./std(X, dims=1)
    @test denormalize(Y, T) ≈ X
    @test_nowarn normalize!(X, T)
    @test X == Y
    @test_nowarn denormalize!(Y, T)
    @test Y ≈ _X

    # X = copy(_X)
    # t = fit(ZScoreTransform, X, dims=1, scale=false)
    # Y = transform(t, X)
    # @test length(t.mean) == 8
    # @test isempty(t.scale)
    # @test Y ≈ X .- mean(X, dims=1)
    # @test reconstruct(t, Y) ≈ X
    # @test transform!(t, X) === X
    # @test isequal(X, Y)
    # @test reconstruct!(t, Y) === Y
    # @test Y ≈ _X

    # X = copy(_X)
    # t = fit(ZScoreTransform, X, dims=1)
    # Y = transform(t, X)
    # @test length(t.mean) == 8
    # @test length(t.scale) == 8
    # @test Y ≈ (X .- mean(X, dims=1)) ./ std(X, dims=1)
    # @test reconstruct(t, Y) ≈ X
    # @test Y ≈ standardize(ZScoreTransform, X, dims=1)
    # @test transform!(t, X) === X
    # @test isequal(X, Y)
    # @test reconstruct!(t, Y) === Y
    # @test Y ≈ _X

    # X = copy(_X)
    # t = fit(ZScoreTransform, X, dims=2)
    # Y = transform(t, X)
    # @test length(t.mean) == 5
    # @test length(t.scale) == 5
    # @test Y ≈ (X .- mean(X, dims=2)) ./ std(X, dims=2)
    # @test reconstruct(t, Y) ≈ X
    # @test Y ≈ standardize(ZScoreTransform, X, dims=2)
    # @test transform!(t, X) === X
    # @test isequal(X, Y)
    # @test reconstruct!(t, Y) === Y
    # @test Y ≈ _X

    # X = copy(_X)
    # t = fit(UnitRangeTransform, X, dims=1, unit=false)
    # Y = transform(t, X)
    # @test length(t.min) == 8
    # @test length(t.scale) == 8
    # @test Y ≈ X ./ (maximum(X, dims=1) .- minimum(X, dims=1))
    # @test reconstruct(t, Y) ≈ X
    # @test transform!(t, X) === X
    # @test isequal(X, Y)
    # @test reconstruct!(t, Y) === Y
    # @test Y ≈ _X

    # X = copy(_X)
    # t = fit(UnitRangeTransform, X, dims=1)
    # Y = transform(t, X)
    # @test isa(t, AbstractDataTransform)
    # @test length(t.min) == 8
    # @test length(t.scale) == 8
    # @test Y ≈ (X .- minimum(X, dims=1)) ./ (maximum(X, dims=1) .- minimum(X, dims=1))
    # @test reconstruct(t, Y) ≈ X
    # @test Y ≈ standardize(UnitRangeTransform, X, dims=1)
    # @test transform!(t, X) === X
    # @test isequal(X, Y)
    # @test reconstruct!(t, Y) === Y
    # @test Y ≈ _X

    # X = copy(_X)
    # t = fit(UnitRangeTransform, X, dims=2)
    # Y = transform(t, X)
    # @test isa(t, AbstractDataTransform)
    # @test length(t.min) == 5
    # @test length(t.scale) == 5
    # @test Y ≈ (X .- minimum(X, dims=2)) ./ (maximum(X, dims=2) .- minimum(X, dims=2))
    # @test reconstruct(t, Y) ≈ X
    # @test transform!(t, X) === X
    # @test isequal(X, Y)
    # @test reconstruct!(t, Y) === Y
    # @test Y ≈ _X

    # # vector
    # X = rand(10)
    # _X = copy(X)

    # t = fit(ZScoreTransform, X, dims=1, center=false, scale=false)
    # Y = transform(t, X)
    # @test transform(t, X) ≈ Y
    # @test reconstruct(t, Y) ≈ X
    # @test transform!(t, X) === X
    # @test isequal(X, Y)
    # @test reconstruct!(t, Y) === Y
    # @test Y ≈ _X

    # X = copy(_X)
    # t = fit(ZScoreTransform, X, dims=1, center=false)
    # Y = transform(t, X)
    # @test Y ≈ X ./ std(X, dims=1)
    # @test transform(t, X) ≈ Y
    # @test reconstruct(t, Y) ≈ X
    # @test transform!(t, X) === X
    # @test isequal(X, Y)
    # @test reconstruct!(t, Y) === Y
    # @test Y ≈ _X

    # X = copy(_X)
    # t = fit(ZScoreTransform, X, dims=1, scale=false)
    # Y = transform(t, X)
    # @test Y ≈ X .- mean(X, dims=1)
    # @test transform(t, X) ≈ Y
    # @test reconstruct(t, Y) ≈ X
    # @test transform!(t, X) === X
    # @test isequal(X, Y)
    # @test reconstruct!(t, Y) === Y
    # @test Y ≈ _X

    # X = copy(_X)
    # t = fit(ZScoreTransform, X, dims=1)
    # Y = transform(t, X)
    # @test Y ≈ (X .- mean(X, dims=1)) ./ std(X, dims=1)
    # @test transform(t, X) ≈ Y
    # @test reconstruct(t, Y) ≈ X
    # @test Y ≈ standardize(ZScoreTransform, X, dims=1)
    # @test transform!(t, X) === X
    # @test isequal(X, Y)
    # @test reconstruct!(t, Y) === Y
    # @test Y ≈ _X

    # X = copy(_X)
    # t = fit(UnitRangeTransform, X, dims=1)
    # Y = transform(t, X)
    # @test Y ≈ (X .- minimum(X, dims=1)) ./ (maximum(X, dims=1) .- minimum(X, dims=1))
    # @test transform(t, X) ≈ Y
    # @test reconstruct(t, Y) ≈ X
    # @test transform!(t, X) === X
    # @test isequal(X, Y)
    # @test reconstruct!(t, Y) === Y
    # @test Y ≈ _X

    # X = copy(_X)
    # t = fit(UnitRangeTransform, X, dims=1, unit=false)
    # Y = transform(t, X)
    # @test Y ≈ X ./ (maximum(X, dims=1) .- minimum(X, dims=1))
    # @test transform(t, X) ≈ Y
    # @test reconstruct(t, Y) ≈ X
    # @test Y ≈ standardize(UnitRangeTransform, X, dims=1, unit=false)
    # @test transform!(t, X) === X
    # @test isequal(X, Y)
    # @test reconstruct!(t, Y) === Y
    # @test Y ≈ _X
end

@testset "3D normalization" begin
    _X = randn(10, 10, 100)
    X = copy(_X)
    T = fit(ZScore, X, dims=3)
    Y = normalize(X, T)
    Z = copy(_X)
    for i ∈ CartesianIndices((axes(Z, 1), axes(Z, 2)))
        x = @view Z[i, :]
        x .= (x.-mean(x))./std(x)
    end
    @test !isnothing(T.p)
    @test length(T.p) == 2
    @test size(T.p[1])[1:2] == size(T.p[2])[1:2] == size(X)[1:2]
    @test Y ≈ Z
    @test denormalize(Y, T) ≈ X
    @test_nowarn normalize!(X, T)
    @test X == Y
    @test_nowarn denormalize!(Y, T)
    @test Y ≈ _X
end

@testset "Sigmoid" begin
    _X = randn(10, 10, 100)
    X = copy(_X)
    T = fit(Sigmoid, X, dims=3)
    Y = normalize(X, T)
    Z = copy(_X)
    for i ∈ CartesianIndices((axes(Z, 1), axes(Z, 2)))
        x = @view Z[i, :]
        x .= 1.0./(1 .+ exp.(.-(x.-mean(x))./std(x)))
    end
    @test !isnothing(T.p)
    @test length(T.p) == 2
    @test size(T.p[1])[1:2] == size(T.p[2])[1:2] == size(X)[1:2]
    @test Y ≈ Z
    @test denormalize(Y, T) ≈ X
    @test_nowarn normalize!(X, T)
    @test X == Y
    @test_nowarn denormalize!(Y, T)
    @test Y ≈ _X
end

@testset "RobustSigmoid" begin
    _X = randn(10, 10, 100)
    X = copy(_X)
    T = fit(RobustSigmoid, X, dims=3)
    Y = normalize(X, T)
    Z = copy(_X)
    for i ∈ CartesianIndices((axes(Z, 1), axes(Z, 2)))
        x = @view Z[i, :]
        x .= 1.0./(1 .+ exp.(.-(x.-median(x))./(SB.iqr(x)./1.35)))
    end
    @test !isnothing(T.p)
    @test length(T.p) == 2
    @test size(T.p[1])[1:2] == size(T.p[2])[1:2] == size(X)[1:2]
    @test Y ≈ Z
    @test denormalize(Y, T) ≈ X
    @test_nowarn normalize!(X, T)
    @test X == Y
    @test_nowarn denormalize!(Y, T)
    @test Y ≈ _X
end

@testset "NaNZScore" begin
    _X = randn(10, 10, 100)
    idxs = rand(1:prod(size(_X)), 100)
    _X[idxs] .= NaN
    X = copy(_X)
    T = fit(nansafe(ZScore), X, dims=3)
    Y = normalize(X, T)
    Z = copy(_X)
    for i ∈ CartesianIndices((axes(Z, 1), axes(Z, 2)))
        x = @view Z[i, :]
        x .= (x.-mean(filter(!isnan, x)))./std(filter(!isnan, x))
    end
    @test !isnothing(T.p)
    @test length(T.p) == 2
    @test size(T.p[1])[1:2] == size(T.p[2])[1:2] == size(X)[1:2]
    @test filter!(!isnan, Y[:]) ≈ filter!(!isnan, Z[:])
    @test filter!(!isnan, denormalize(Y, T)[:]) ≈ filter!(!isnan, X[:])
    @test_nowarn normalize!(X, T)
    @test filter!(!isnan, X[:]) == filter!(!isnan, Y[:])
    @test_nowarn denormalize!(Y, T)
    @test filter!(!isnan, Y[:]) ≈ filter!(!isnan, _X[:])
end

@testset "NaNSigmoid" begin
    _X = randn(10, 10, 100)
    idxs = rand(1:prod(size(_X)), 100)
    _X[idxs] .= NaN
    X = copy(_X)
    T = fit(nansafe(Sigmoid), X, dims=3)
    Y = normalize(X, T)
    Z = copy(_X)
    for i ∈ CartesianIndices((axes(Z, 1), axes(Z, 2)))
        x = @view Z[i, :]
        x .= 1.0./(1 .+ exp.(.-(x.-mean(filter(!isnan, x)))./(std(filter(!isnan, x)))))
    end
    @test !isnothing(T.p)
    @test length(T.p) == 2
    @test size(T.p[1])[1:2] == size(T.p[2])[1:2] == size(X)[1:2]
    @test filter!(!isnan, Y[:]) ≈ filter!(!isnan, Z[:])
    @test filter!(!isnan, denormalize(Y, T)[:]) ≈ filter!(!isnan, X[:])
    @test_nowarn normalize!(X, T)
    @test filter!(!isnan, X[:]) == filter!(!isnan, Y[:])
    @test_nowarn denormalize!(Y, T)
    @test filter!(!isnan, Y[:]) ≈ filter!(!isnan, _X[:])
end

@testset "ND normalization" begin
    Nmax = 5
    for _ = 1:10
        ds = rand(2:Nmax)
        sz = rand(2:ds, ds)
        _X = rand(sz...)
        Nnorm = rand(1:Nmax÷2)
        normdims = unique(rand(1:ds, Nnorm))
        notnormdims = setdiff(1:ndims(_X), normdims)
        X = copy(_X)
        T = fit(ZScore, X, dims=normdims)
        Y = normalize(X, T)
        @test !isnothing(T.p)
        @test length(T.p) == 2
        @test size(T.p[1])[notnormdims] == size(T.p[2])[notnormdims] == size(X)[notnormdims]
        @test denormalize(Y, T) ≈ X
        @test_nowarn normalize!(X, T)
        @test X == Y
        @test_nowarn denormalize!(Y, T)
        @test Y ≈ _X
    end
end

@testset "StatsBase comparison" begin
    X = rand(1000, 50)
    Y = normalize(X, ZScore; dims=1)
    Z = SB.standardize(SB.ZScoreTransform, X; dims=1)
    @test Y ≈ Z
end

println("\nNormalization.jl")
display(@benchmark normalize(rand(1000, 10000), ZScore; dims=1))
println("\nStatsBase.jl")
display(@benchmark SB.standardize(SB.ZScoreTransform, rand(1000, 10000); dims=1))


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

## Basic functions ##

"""
    AbstractQuasiArray{T,N}

Supertype for `N`-dimensional arrays (or array-like types) with elements of type `T`.
[`Array`](@ref) and other types are subtypes of this. See the manual section on the
[`AbstractQuasiArray` interface](@ref man-interface-array).
"""
AbstractQuasiArray

convert(::Type{T}, a::T) where {T<:AbstractQuasiArray} = a
convert(::Type{AbstractQuasiArray{T}}, a::AbstractQuasiArray) where {T} = AbstractQuasiArray{T}(a)
convert(::Type{AbstractQuasiArray{T,N}}, a::AbstractQuasiArray{<:Any,N}) where {T,N} = AbstractQuasiArray{T,N}(a)

Array{T}(a::AbstractQuasiArray) where T = T[a[k] for k in eachindex(a)]
Array{T,N}(a::AbstractQuasiArray{<:Any,N}) where {T,N} = Array{T}(a)
Array(a::AbstractQuasiArray{T}) where T = Array{T}(a)
Matrix(a::AbstractQuasiMatrix{T}) where T = Array{T}(a)
Vector(a::AbstractQuasiVector{T}) where T = Array{T}(a)
AbstractArray(a::AbstractQuasiArray) = Array(a)
AbstractMatrix(a::AbstractQuasiMatrix) = Matrix(a)
AbstractVector(a::AbstractQuasiVector) = Vector(a)
AbstractArray{T}(a::AbstractQuasiArray) where T = Array{T}(a)
AbstractMatrix{T}(a::AbstractQuasiMatrix) where T = Matrix{T}(a)
AbstractVector{T}(a::AbstractQuasiVector) where T = Vector{T}(a)

convert(::Type{Array{T}}, a::AbstractQuasiArray) where T = Array{T}(a)
convert(::Type{Array{T,N}}, a::AbstractQuasiArray{<:Any,N}) where {T,N} = Array{T}(a)
convert(::Type{Array}, a::AbstractQuasiArray{T}) where T = Array{T}(a)
convert(::Type{Matrix}, a::AbstractQuasiMatrix{T}) where T = Array{T}(a)
convert(::Type{Vector}, a::AbstractQuasiVector{T}) where T = Array{T}(a)

convert(::Type{AbstractArray{T}}, a::AbstractQuasiArray) where T = Array{T}(a)
convert(::Type{AbstractArray}, a::AbstractQuasiArray{T}) where T = convert(AbstractArray{T}, a)
convert(::Type{AbstractArray{T,N}}, a::AbstractQuasiArray{<:Any,N}) where {T,N} = convert(AbstractArray{T}, a)
convert(::Type{AbstractMatrix}, a::AbstractQuasiMatrix) = convert(AbstractArray, a)
convert(::Type{AbstractVector}, a::AbstractQuasiVector) = convert(AbstractArray, a)

"""
    size(A::AbstractQuasiArray, [dim])

Return a tuple containing the dimensions of `A`. Optionally you can specify a
dimension to just get the length of that dimension.

Note that `size` may not be defined for arrays with non-standard indices, in which case [`axes`](@ref)
may be useful. See the manual chapter on [arrays with custom indices](@ref man-custom-indices).

# Examples
```jldoctest
julia> A = fill(1, (2,3,4));

julia> size(A)
(2, 3, 4)

julia> size(A, 2)
3
```
"""
size(t::AbstractQuasiArray{T,N}, d) where {T,N} = d <= N ? size(t)[d] : 1

"""
    axes(A, d)

Return the valid range of indices for array `A` along dimension `d`.

See also [`size`](@ref), and the manual chapter on [arrays with custom indices](@ref man-custom-indices).

# Examples
```jldoctest
julia> A = fill(1, (5,6,7));

julia> axes(A, 2)
Base.OneTo(6)
```
"""
function axes(A::AbstractQuasiArray{T,N}, d) where {T,N}
    @_inline_meta
    d <= N ? axes(A)[d] : OneTo(1)
end

# Performance optimization: get rid of a branch on `d` in `axes(A, d)`
# for d=1. 1d arrays are heavily used, and the first dimension comes up
# in other applications.
axes1(A::AbstractQuasiArray{<:Any,0}) = OneTo(1)
axes1(A::AbstractQuasiArray) = (@_inline_meta; axes(A)[1])

keys(a::AbstractQuasiArray) = CartesianIndices(axes(a))
keys(a::AbstractQuasiVector) = LinearIndices(a)

eltype(::Type{<:AbstractQuasiArray{E}}) where {E} = @isdefined(E) ? E : Any
elsize(A::AbstractQuasiArray) = elsize(typeof(A))

"""
    ndims(A::AbstractQuasiArray) -> Integer

Return the number of dimensions of `A`.

# Examples
```jldoctest
julia> A = fill(1, (3,4,5));

julia> ndims(A)
3
```
"""
ndims(::AbstractQuasiArray{T,N}) where {T,N} = N
ndims(::Type{<:AbstractQuasiArray{T,N}}) where {T,N} = N

"""
    length(collection) -> Integer

Return the number of elements in the collection.

Use [`lastindex`](@ref) to get the last valid index of an indexable collection.

# Examples
```jldoctest
julia> length(1:5)
5

julia> length([1, 2, 3, 4])
4

julia> length([1 2; 3 4])
4
```
"""
length

"""
    length(A::AbstractQuasiArray)

Return the number of elements in the array, defaults to `prod(size(A))`.

# Examples
```jldoctest
julia> length([1, 2, 3, 4])
4

julia> length([1 2; 3 4])
4
```
"""
length(t::AbstractQuasiArray) = (@_inline_meta; prod(size(t)))

# eachindex iterates over all indices. IndexCartesian definitions are later.
eachindex(A::AbstractQuasiVector) = (@_inline_meta(); axes1(A))

"""
    eachindex(A...)

Create an iterable object for visiting each index of an `AbstractQuasiArray` `A` in an efficient
manner. For array types that have opted into fast linear indexing (like `Array`), this is
simply the range `1:length(A)`. For other array types, return a specialized Cartesian
range to efficiently index into the array with indices specified for every dimension. For
other iterables, including strings and dictionaries, return an iterator object
supporting arbitrary index types (e.g. unevenly spaced or non-integer indices).

If you supply more than one `AbstractQuasiArray` argument, `eachindex` will create an
iterable object that is fast for all arguments (a [`UnitRange`](@ref)
if all inputs have fast linear indexing, a [`CartesianIndices`](@ref)
otherwise).
If the arrays have different sizes and/or dimensionalities, `eachindex` will return an
iterable that spans the largest range along each dimension.

# Examples
```jldoctest
julia> A = [1 2; 3 4];

julia> for i in eachindex(A) # linear indexing
           println(i)
       end
1
2
3
4

julia> for i in eachindex(view(A, 1:2, 1:1)) # Cartesian indexing
           println(i)
       end
CartesianIndex(1, 1)
CartesianIndex(2, 1)
```
"""
eachindex(A::AbstractQuasiArray) = (@_inline_meta(); eachindex(IndexStyle(A), A))

function eachindex(A::AbstractQuasiArray, B::AbstractQuasiArray)
    @_inline_meta
    eachindex(IndexStyle(A,B), A, B)
end
function eachindex(A::AbstractQuasiArray, B::AbstractQuasiArray...)
    @_inline_meta
    eachindex(IndexStyle(A,B...), A, B...)
end
eachindex(::IndexLinear, A::AbstractQuasiArray) = (@_inline_meta; OneTo(length(A)))
eachindex(::IndexLinear, A::AbstractQuasiVector) = (@_inline_meta; axes1(A))
function eachindex(::IndexLinear, A::AbstractQuasiArray, B::AbstractQuasiArray...)
    @_inline_meta
    indsA = eachindex(IndexLinear(), A)
    _all_match_first(X->eachindex(IndexLinear(), X), indsA, B...) ||
        throw_eachindex_mismatch(IndexLinear(), A, B...)
    indsA
end

# keys with an IndexStyle
keys(s::IndexStyle, A::AbstractQuasiArray, B::AbstractQuasiArray...) = eachindex(s, A, B...)

"""
    lastindex(collection) -> index
    lastindex(collection, d) -> index

Return the last index of `collection`. If `d` is given, return the last index of `collection` along dimension `d`.

The syntaxes `A[end]` and `A[end, end]` lower to `A[lastindex(A)]` and
`A[lastindex(A, 1), lastindex(A, 2)]`, respectively.

# Examples
```jldoctest
julia> lastindex([1,2,4])
3

julia> lastindex(rand(3,4,5), 2)
4
```
"""
lastindex(a::AbstractQuasiArray) = (@_inline_meta; last(eachindex(IndexLinear(), a)))
lastindex(a::AbstractQuasiArray, d) = (@_inline_meta; last(axes(a, d)))

"""
    firstindex(collection) -> index
    firstindex(collection, d) -> index

Return the first index of `collection`. If `d` is given, return the first index of `collection` along dimension `d`.

# Examples
```jldoctest
julia> firstindex([1,2,4])
1

julia> firstindex(rand(3,4,5), 2)
1
```
"""
firstindex(a::AbstractQuasiArray) = (@_inline_meta; first(eachindex(IndexLinear(), a)))
firstindex(a::AbstractQuasiArray, d) = (@_inline_meta; first(axes(a, d)))

first(a::AbstractQuasiArray) = a[first(eachindex(a))]
stride(A::AbstractQuasiArray, k::Integer) = strides(A)[k]

function isassigned(a::AbstractQuasiArray, i...)
    try
        a[i...]
        true
    catch e
        if isa(e, BoundsError) || isa(e, UndefRefError)
            return false
        else
            rethrow(e)
        end
    end
end

function checkbounds(::Type{Bool}, A::AbstractQuasiArray, I...)
    @_inline_meta
    checkbounds_indices(Bool, axes(A), I)
end

# As a special extension, allow using logical arrays that match the source array exactly
function checkbounds(::Type{Bool}, A::AbstractQuasiArray{<:Any,N}, I::AbstractQuasiArray{Bool,N}) where N
    @_inline_meta
    axes(A) == axes(I)
end

function checkbounds(A::AbstractQuasiArray, I...)
    @_inline_meta
    checkbounds(Bool, A, to_indices(A,I)...) || throw_boundserror(A, I)
    nothing
end

const QuasiDimOrInd = Union{Number, AbstractQuasiOrVector{<:Number}}

similar(a::AbstractQuasiArray{T}) where {T}                             = similar(a, T)
similar(a::AbstractQuasiArray, ::Type{T}) where {T}                     = similar(a, T, axes(a))
similar(a::AbstractQuasiArray{T}, dims::Tuple) where {T}                = similar(a, T, dims)
similar(a::AbstractQuasiArray{T}, dims::QuasiDimOrInd...) where {T}          = similar(a, T, dims)
similar(a::AbstractQuasiArray, ::Type{T}, dims::QuasiDimOrInd...) where {T}  = similar(a, T, dims)
similar(::Type{<:AbstractQuasiArray{T}}, shape::NTuple{N,AbstractQuasiOrVector}) where {N,T} =
    QuasiArray{T,N}(undef, convert.(AbstractVector, shape))
similar(a::AbstractQuasiArray, ::Type{T}, dims::NTuple{N,AbstractQuasiOrVector}) where {T,N} =
    QuasiArray{T,N}(undef, convert.(AbstractVector, dims))
similar(a::AbstractQuasiArray, ::Type{T}, dims::Vararg{AbstractQuasiOrVector,N}) where {T,N} =
    QuasiArray{T,N}(undef, convert.(AbstractVector, dims))

similar(a::AbstractQuasiArray{T}, m::Int) where {T}              = Vector{T}(undef, m)
similar(a::AbstractQuasiArray, T::Type, dims::Dims{N}) where {N} = Array{T,N}(undef, dims)
similar(a::AbstractQuasiArray{T}, dims::Dims{N}) where {T,N}     = Array{T,N}(undef, dims)

similar(::Type{T}, dims::QuasiDimOrInd...) where {T<:AbstractQuasiArray} = similar(T, dims)
similar(::Type{T}, shape::Tuple{AbstractQuasiVector,Vararg{QuasiDimOrInd}}) where {N,T<:AbstractArray{<:Any,N}} =
    QuasiArray{T,N}(undef, convert.(AbstractVector, shape))

empty(a::AbstractQuasiVector{T}, ::Type{U}=T) where {T,U} = Vector{U}()

# like empty, but should return a mutable collection, a Vector by default
emptymutable(a::AbstractQuasiVector{T}, ::Type{U}=T) where {T,U} = Vector{U}()

## copy between abstract arrays - generally more efficient
## since a single index variable can be used.

copyto!(dest::AbstractQuasiArray, src::AbstractQuasiArray) =
    copyto!(IndexStyle(dest), dest, IndexStyle(src), src)

function copyto!(::IndexStyle, dest::AbstractQuasiArray, ::IndexStyle, src::AbstractQuasiArray)
    destinds, srcinds = LinearIndices(dest), LinearIndices(src)
    isempty(srcinds) || (checkbounds(Bool, destinds, first(srcinds)) && checkbounds(Bool, destinds, last(srcinds))) ||
        throw(BoundsError(dest, srcinds))
    @inbounds for i in srcinds
        dest[i] = src[i]
    end
    return dest
end

function copyto!(::IndexStyle, dest::AbstractQuasiArray, ::IndexCartesian, src::AbstractQuasiArray)
    axes(dest) == axes(src) || throw(BoundsError(dest, axes(src)))
    i = 0
    @inbounds for i in eachindex(src)
        dest[i] = src[i]
    end
    return dest
end

function copy(a::AbstractQuasiArray)
    @_propagate_inbounds_meta
    copymutable(a)
end

function copymutable(a::AbstractQuasiArray)
    @_propagate_inbounds_meta
    copyto!(similar(a), a)
end

## iteration support for arrays by iterating over `eachindex` in the array ##
# Allows fast iteration by default for both IndexLinear and IndexCartesian arrays

# While the definitions for IndexLinear are all simple enough to inline on their
# own, IndexCartesian's CartesianIndices is more complicated and requires explicit
# inlining.
function iterate(A::AbstractQuasiArray, state=(eachindex(A),))
    y = iterate(state...)
    y === nothing && return nothing
    A[y[1]], (state[1], tail(y)...)
end

isempty(a::AbstractQuasiArray) = (length(a) == 0)

getindex(A::AbstractQuasiArray, I...) = _getindex(indextype(A), A, I)

function _getindex(::Type{IND}, A::AbstractQuasiArray, I) where IND
    @_propagate_inbounds_meta
    error_if_canonical_getindex(IndexStyle(A), A, I)
    _getindex(IND, IndexStyle(A), A, to_indices(A, I))
end
function unsafe_getindex(A::AbstractQuasiArray, I...)
    @_inline_meta
    @inbounds r = getindex(A, I...)
    r
end

## Internal definitions
_getindex(_, ::IndexStyle, A::AbstractQuasiArray, I) = layout_getindex(A, I...)

## IndexCartesian Scalar indexing: Canonical method is full dimensionality of indices
function _getindex(::Type{IND}, ::IndexCartesian, A::AbstractQuasiArray, I::IND) where {M,IND}
    @_inline_meta
    @boundscheck checkbounds(A, I...) # generally _to_subscript_indices requires bounds checking
    @inbounds r = getindex(A, _to_subscript_indices(A, I...)...)
    r
end

error_if_canonical_getindex(::IndexCartesian, A::AbstractQuasiArray{T,N}, I::Tuple) where {T,N} =
    _error_if_canonical_getindex(indextype(A), A, I)

error_if_canonical_getindex(::IndexLinear, A::AbstractQuasiArray{T,N}, I::Tuple) where {T,N} =
    _error_if_canonical_getindex(indextype(A), A, I)

_error_if_canonical_getindex(::Type{IND}, A::AbstractQuasiArray{T,N}, I::IND) where {T,N,IND} =
    error("getindex not defined for ", typeof(A))

_error_if_canonical_getindex(::Type, ::AbstractQuasiArray, ::Any...) = nothing

_to_subscript_indices(A::AbstractQuasiArray, i) = (@_inline_meta; _unsafe_ind2sub(A, i))
_to_subscript_indices(A::AbstractQuasiArray{T,N}) where {T,N} = (@_inline_meta; fill_to_length((), 1, Val(N)))
_to_subscript_indices(A::AbstractQuasiArray{T,0}) where {T} = ()
_to_subscript_indices(A::AbstractQuasiArray{T,0}, i) where {T} = ()
_to_subscript_indices(A::AbstractQuasiArray{T,0}, I...) where {T} = ()
function _to_subscript_indices(A::AbstractQuasiArray{T,N}, I...) where {T,N}
    @_inline_meta
    J, Jrem = IteratorsMD.split(I, Val(N))
    _to_subscript_indices(A, J, Jrem)
end
_to_subscript_indices(A::AbstractQuasiArray, J::Tuple, Jrem::Tuple{}) =
    __to_subscript_indices(A, axes(A), J, Jrem)
function __to_subscript_indices(A::AbstractQuasiArray,
        ::Tuple{AbstractUnitRange,Vararg{AbstractUnitRange}}, J::Tuple, Jrem::Tuple{})
    @_inline_meta
    (J..., map(first, tail(_remaining_size(J, axes(A))))...)
end
_to_subscript_indices(A::AbstractQuasiArray{T,N}, I::Vararg{Any,N}) where {T,N} = I

## Setindex! is defined similarly. We first dispatch to an internal _setindex!
# function that allows dispatch on array storage

setindex!(A::AbstractQuasiArray, v, I...) = _setindex!(indextype(A), A, v, I)

function _setindex!(::Type{IND}, A::AbstractQuasiArray, v, I) where IND
    @_propagate_inbounds_meta
    error_if_canonical_setindex(IndexStyle(A), A, I)
    _setindex!(IND, IndexStyle(A), A, v, to_indices(A, I))
end

error_if_canonical_setindex(::IndexCartesian, A::AbstractQuasiArray{T,N}, I::Tuple) where {T,N} =
    _error_if_canonical_setindex(indextype(A), A, I)

error_if_canonical_setindex(::IndexLinear, A::AbstractQuasiArray{T,N}, I::Tuple) where {T,N} =
    _error_if_canonical_setindex(indextype(A), A, I)

_error_if_canonical_setindex(::Type{IND}, A::AbstractQuasiArray{T,N}, I::IND) where {T,N,IND} =
    error("setindex! not defined for ", typeof(A))

_error_if_canonical_setindex(::Type, ::AbstractQuasiArray, ::Any...) = nothing
## Internal definitions
_setindex!(::Type, ::IndexStyle, A::AbstractQuasiArray, v, I) =
    error("setindex! for $(typeof(A)) with types $(typeof(I)) is not supported")

# IndexCartesian Scalar indexing
function _setindex!(::Type{IND}, ::IndexCartesian, A::AbstractQuasiArray{T,N}, v, I::NTuple{N,Any}) where {T,N,IND}
    @_propagate_inbounds_meta
    setindex!(A, v, I...)
end
function _setindex!(::Type{IND}, ::IndexCartesian, A::AbstractQuasiArray, v, I::NTuple{M,Any}) where {M,IND}
    @_inline_meta
    @boundscheck checkbounds(A, I...)
    @inbounds r = setindex!(A, v, _to_subscript_indices(IND, A, I...)...)
    r
end

parent(a::AbstractQuasiArray) = a

unalias(dest, A::AbstractQuasiArray) = mightalias(dest, A) ? unaliascopy(A) : A
unaliascopy(A::AbstractQuasiArray)::typeof(A) = (@_noinline_meta; _unaliascopy(A, copy(A)))

"""
    Base.mightalias(A::AbstractQuasiArray, B::AbstractQuasiArray)

Perform a conservative test to check if quasi arrays `A` and `B` might share the same memory.

By default, this simply checks if either of the arrays reference the same memory
regions, as identified by their [`Base.dataids`](@ref).
"""
mightalias(A::AbstractQuasiArray, B::AbstractQuasiArray) = !_isdisjoint(dataids(A), dataids(B))

"""
    Base.dataids(A::AbstractQuasiArray)

Return a tuple of `UInt`s that represent the mutable data segments of an array.

Custom arrays that would like to opt-in to aliasing detection of their component
parts can specialize this method to return the concatenation of the `dataids` of
their component parts.  A typical definition for an array that wraps a parent is
`Base.dataids(C::CustomArray) = dataids(C.parent)`.
"""
dataids(A::AbstractQuasiArray) = (UInt(objectid(A)),)

## structured matrix methods ##
replace_in_print_matrix(A::AbstractQuasiMatrix,i,j,s::AbstractString) = s
replace_in_print_matrix(A::AbstractQuasiVector,i,j,s::AbstractString) = s

## Concatenation ##
eltypeof(x::AbstractQuasiArray) = eltype(x)

## Reductions and accumulates ##

function isequal(A::AbstractQuasiArray, B::AbstractQuasiArray)
    if A === B return true end
    if axes(A) != axes(B)
        return false
    end
    for (a, b) in zip(A, B)
        if !isequal(a, b)
            return false
        end
    end
    return true
end

function cmp(A::AbstractQuasiVector, B::AbstractQuasiVector)
    for (a, b) in zip(A, B)
        if !isequal(a, b)
            return isless(a, b) ? -1 : 1
        end
    end
    return cmp(length(A), length(B))
end

isless(A::AbstractQuasiVector, B::AbstractQuasiVector) = cmp(A, B) < 0

function (==)(A::AbstractQuasiArray, B::AbstractQuasiArray)
    if axes(A) != axes(B)
        return false
    end
    anymissing = false
    for (a, b) in zip(A, B)
        eq = (a == b)
        if ismissing(eq)
            anymissing = true
        elseif !eq
            return false
        end
    end
    return anymissing ? missing : true
end

_lookup(ind, r::Inclusion) = ind

_ind2sub(dims::NTuple{N,Any}, ind) where N = (@_inline_meta; _ind2sub_recurse(dims, ind-1))
_ind2sub(inds::QuasiIndices, ind)     = (@_inline_meta; _ind2sub_recurse(inds, ind-1))
_ind2sub(inds::Tuple{Inclusion{<:Any},AbstractUnitRange{<:Integer}}, ind)     = (@_inline_meta; _ind2sub_recurse(inds, ind-1))
_ind2sub(inds::Tuple{AbstractUnitRange{<:Integer},Inclusion{<:Any}}, ind)     = (@_inline_meta; _ind2sub_recurse(inds, ind-1))

function _ind2sub(inds::Union{NTuple{N,Any},QuasiIndices{N}}, ind::AbstractQuasiVector) where N
    M = length(ind)
    t = ntuple(n->similar(ind),Val(N))
    for (i,idx) in pairs(IndexLinear(), ind)
        sub = _ind2sub(inds, idx)
        for j = 1:N
            t[j][i] = sub[j]
        end
    end
    t
end

## map over arrays ##

## transform any set of dimensions
## dims specifies which dimensions will be transformed. for example
## dims==1:2 will call f on all slices A[:,:,...]

function mapslices(f, A::AbstractQuasiArray; dims)
    if isempty(dims)
        return map(f,A)
    end
    if !isa(dims, AbstractQuasiVector)
        dims = [dims...]
    end

    dimsA = [axes(A)...]
    ndimsA = ndims(A)
    alldims = [1:ndimsA;]

    otherdims = setdiff(alldims, dims)

    idx = Any[first(ind) for ind in axes(A)]
    itershape   = tuple(dimsA[otherdims]...)
    for d in dims
        idx[d] = Slice(axes(A, d))
    end

    # Apply the function to the first slice in order to determine the next steps
    Aslice = A[idx...]
    r1 = f(Aslice)
    # In some cases, we can re-use the first slice for a dramatic performance
    # increase. The slice itself must be mutable and the result cannot contain
    # any mutable containers. The following errs on the side of being overly
    # strict (#18570 & #21123).
    safe_for_reuse = isa(Aslice, StridedArray) &&
                     (isa(r1, Number) || (isa(r1, AbstractQuasiArray) && eltype(r1) <: Number))

    # determine result size and allocate
    Rsize = copy(dimsA)
    # TODO: maybe support removing dimensions
    if !isa(r1, AbstractQuasiArray) || ndims(r1) == 0
        # If the result of f on a single slice is a scalar then we add singleton
        # dimensions. When adding the dimensions, we have to respect the
        # index type of the input array (e.g. in the case of OffsetArrays)
        tmp = similar(Aslice, typeof(r1), reduced_indices(Aslice, 1:ndims(Aslice)))
        tmp[firstindex(tmp)] = r1
        r1 = tmp
    end
    nextra = max(0, length(dims)-ndims(r1))
    if eltype(Rsize) == Int
        Rsize[dims] = [size(r1)..., ntuple(d->1, nextra)...]
    else
        Rsize[dims] = [axes(r1)..., ntuple(d->OneTo(1), nextra)...]
    end
    R = similar(r1, tuple(Rsize...,))

    ridx = Any[map(first, axes(R))...]
    for d in dims
        ridx[d] = axes(R,d)
    end

    concatenate_setindex!(R, r1, ridx...)

    nidx = length(otherdims)
    indices = Iterators.drop(CartesianIndices(itershape), 1) # skip the first element, we already handled it
    inner_mapslices!(safe_for_reuse, indices, nidx, idx, otherdims, ridx, Aslice, A, f, R)
end

concatenate_setindex!(R, X::AbstractQuasiArray, I...) = (R[I...] = X)

## 1 argument

function map!(f::F, dest::AbstractQuasiArray, A::AbstractQuasiArray) where F
    for (i,j) in zip(eachindex(dest),eachindex(A))
        dest[i] = f(A[j])
    end
    return dest
end

# map on collections
map(f, A::AbstractQuasiArray) = collect_similar(A, Generator(f,A))

## 2 argument
function map!(f::F, dest::AbstractQuasiArray, A::AbstractQuasiArray, B::AbstractQuasiArray) where F
    for (i, j, k) in zip(eachindex(dest), eachindex(A), eachindex(B))
        dest[i] = f(A[j], B[k])
    end
    return dest
end

function map_n!(f::F, dest::AbstractQuasiArray, As) where F
    for i = LinearIndices(As[1])
        dest[i] = f(ith_all(i, As)...)
    end
    return dest
end

map!(f::F, dest::AbstractQuasiArray, As::AbstractQuasiArray...) where {F} = map_n!(f, dest, As)

## hashing AbstractQuasiArray ##

function hash(A::AbstractQuasiArray, h::UInt)
    h = hash(AbstractQuasiArray, h)
    # Axes are themselves AbstractQuasiArrays, so hashing them directly would stack overflow
    # Instead hash the tuple of firsts and lasts along each dimension
    h = hash(map(first, axes(A)), h)
    h = hash(map(last, axes(A)), h)
    isempty(A) && return h

    # Goal: Hash approximately log(N) entries with a higher density of hashed elements
    # weighted towards the end and special consideration for repeated values. Colliding
    # hashes will often subsequently be compared by equality -- and equality between arrays
    # works elementwise forwards and is short-circuiting. This means that a collision
    # between arrays that differ by elements at the beginning is cheaper than one where the
    # difference is towards the end. Furthermore, blindly choosing log(N) entries from a
    # sparse array will likely only choose the same element repeatedly (zero in this case).

    # To achieve this, we work backwards, starting by hashing the last element of the
    # array. After hashing each element, we skip `fibskip` elements, where `fibskip`
    # is pulled from the Fibonacci sequence -- Fibonacci was chosen as a simple
    # ~O(log(N)) algorithm that ensures we don't hit a common divisor of a dimension
    # and only end up hashing one slice of the array (as might happen with powers of
    # two). Finally, we find the next distinct value from the one we just hashed.

    # This is a little tricky since skipping an integer number of values inherently works
    # with linear indices, but `findprev` uses `keys`. Hoist out the conversion "maps":
    ks = keys(A)
    key_to_linear = LinearIndices(ks) # Index into this map to compute the linear index
    linear_to_key = vec(ks)           # And vice-versa

    # Start at the last index
    keyidx = last(ks)
    linidx = key_to_linear[keyidx]
    fibskip = prevfibskip = oneunit(linidx)
    n = 0
    while true
        n += 1
        # Hash the current key-index and its element
        elt = A[keyidx]
        h = hash(keyidx=>elt, h)

        # Skip backwards a Fibonacci number of indices -- this is a linear index operation
        linidx = key_to_linear[keyidx]
        linidx <= fibskip && break
        linidx -= fibskip
        keyidx = linear_to_key[linidx]

        # Only increase the Fibonacci skip once every N iterations. This was chosen
        # to be big enough that all elements of small arrays get hashed while
        # obscenely large arrays are still tractable. With a choice of N=4096, an
        # entirely-distinct 8000-element array will have ~75% of its elements hashed,
        # with every other element hashed in the first half of the array. At the same
        # time, hashing a `typemax(Int64)`-length Float64 range takes about a second.
        if rem(n, 4096) == 0
            fibskip, prevfibskip = fibskip + prevfibskip, fibskip
        end

        # Find a key index with a value distinct from `elt` -- might be `keyidx` itself
        keyidx = findprev(!isequal(elt), A, keyidx)
        keyidx === nothing && break
    end

    return h
end

##
# show
##

show(io::IO, A::AbstractQuasiArray) = summary(io, A)


text:
module Mrate

export modelm, modelm_j, modelm_h
export sir_Γ!, sird_Γ!

using SpecialFunctions, MyFunctions

function modelm_def(t;p=[0.25,0],ρ=[2,1])::Float64
    #compute m and verify that length of p is even
    m,r = fldmod( length(p), 2) .-(1,0)
    @assert r==0 "Length of p must be even"
    @assert iseven( length( ρ ) ) "Length of ρ must be even"

    a0 = p[1]*t -p[2] # evaluate common part
    if m == 0
        return σ( a0 )
    elseif m == 1   # hand typed for m=1, but the general case would works
        am = (p[3]-p[1])*IΓ(t-p[4];ρ=ρ)
        return σ( a0 +am )
    else
        am = 0.0
        for i=1:m
            j=2*i
            am += (p[j+1]-p[j-1])*IΓ(t-p[j+2];ρ=ρ)
        end
        return σ( a0 +am )
    end
end
"""
    modelm(t,p,ρ)

Compute the solution for the Losgistic differential equation with m+1
    constant rates ``r = {r_0, r_1, \\ldots, r_m }`` and off-set
    ``t_0``. The switching between the constant rates occures at
    times ``t = {r_1, \\ldots, t_m }`` and is modelled according to a
    Gamma distribution

Requires SpecialFunctions package and makes use of gamma_inc
    function.

modelm invokes the source function modelm_def to allow for trivial
    broadcast over the variable x, model_def is not exported by
    default.

Note that LsqFit.curve_fit requires a model function that supports
    two arguments. The first, an abstract array for the indepenendent
    variables. The second, another abstract array for the parameters
    that will be optimized.

# Sources:
1. private notes contact me via https://github.com/f-hipolito

# Examples
```julia-repl
julia> p0 = [0.25, 0. ];

julia> ρ0=[6.2,.97];

julia> modelm([-10, 0.0, 10.0], p0, ρ0)
3-element Array{Float64,1}:
 0.07585818002124355
 0.5
 0.9241418199787566

julia> p2 = [0.25, 18.8, 0.075, 45. ];

julia> modelm([0.0, 50.0, 200.0], p2, ρ0)
3-element Array{Float64,1}:
 6.8432709753876305e-9
 0.0017103219749579988
 0.9944784329234642
```
"""
@. function modelm(t,p,ρ)
    modelm_def(t;p=p,ρ=ρ)
end

"""
    Rt!(R,m,x,p,ρ)

Evaluates the argument of the m-rate logistic function, where R is an
array (mutable by default).

# Sources:
1. private notes contact me via https://github.com/f-hipolito

# Examples
```julia-repl
julia> Rt!(R0,0,[-10, 0.0, 10.0],[0.25, 0. ],[6.2,.97])
3-element Array{Float64,1}:
 -2.5
  0.0
  2.5
```
"""
function Rt!(R,m,x::Array,p,ρ)
    for i=0:m
        if i==0
            R[:] += p[1] .*x .-p[2]
        elseif i>0
            j=2*i
            R[:] += ( p[j+1] -p[j-1] ) .*IΓ.( x .-p[j+2]; ρ=ρ )
        end
    end
    return R
end
"""
    Rt(m,x,p,ρ)

Evaluates the argument of the m-rate logistic function, where R is a
scalar Float64.

# Sources:
1. private notes contact me via https://github.com/f-hipolito

# Examples
```julia-repl
julia> Rt(0,10,[.25, 0],[6.2,0.97])
2.5
```
"""
function Rt(m,x::Number,p,ρ)
    R::Float64 = 0
    for i=0:m
        if i==0
            R += p[1] .*x .-p[2]
        elseif i>0
            j=2*i
            R += ( p[j+1] -p[j-1] ) .*IΓ.( x .-p[j+2]; ρ=ρ )
        end
    end
    return R
end

"""
    dr( i,x,p,ρ )

Evaluates the ith partial derivative of the argument of the m-rate
logistic function, with respect to the 'rate' parameters

# Sources:
1. private notes contact me via https://github.com/f-hipolito

# Examples
```julia-repl
julia> dr( 0,0,[0, 25, 100],p0,ρ0 )
3-element Array{Float64,1}:
   0.0
  25.0
 100.0

julia> dr( 0,0,25,p0,ρ0 )
25.0
```
"""
function dr( i,m,x::Array,p,ρ )
    out = zeros(Float64, length(x) )
    if i == 0 && m == 0
        @. out = x
    elseif i == 0 && m >= 1
        out = x -IΓ.( x .-p[4]; ρ=ρ )
    elseif i >= 1 && i < m
        j=2*i
        out = IΓ.( x .-p[j+2]; ρ=ρ ) -IΓ.( x .-p[j+4]; ρ=ρ )
    elseif i != 0 && i == m
        out = IΓ.( x .-p[2*i+2]; ρ=ρ )
    end
    return out
end
function dr( i,m,x::Number,p,ρ )
    out::Float64 = 0
    if i == 0 && m == 0
        out = x
    elseif i == 0 && m >= 1
        out = x -IΓ( x -p[4]; ρ=ρ )
    elseif i >= 1 && i < m
        j=2*i
        out = IΓ( x -p[j+2]; ρ=ρ ) -IΓ( x -p[j+4]; ρ=ρ )
    elseif i != 0 && i == m
        out = IΓ( x -p[2*i+2]; ρ=ρ )
    end
    return out
end

"""
    dt( i,x,p,ρ )

Evaluates the ith partial derivative of the argument of the m-rate
logistic function, with respect to the 'time' parameters

# Sources:
1. private notes contact me via https://github.com/f-hipolito

# Examples
```julia-repl
julia> dt( 0,[0.0, 50.0, 200.0],[0.25, 0. ],[6.2,.97] )
3-element Array{Float64,1}:
 -1.0
 -1.0
 -1.0

julia> dt( 0,25,p0,ρ0 )
-1.0
```
"""
function dt( i,x::Array,p,ρ )
    out = zeros(Float64, length(x) )
    if i == 0
        @. out = -1
    else
        j=2*i
        @. out = (p[j-1]-p[j+1]) *FΓ( x-p[j+2]; p=ρ, IND=0 )
    end
    return out
end
function dt( i,x::Number,p,ρ )
    out::Float64 = 0
    if i == 0
        out = -1
    else
        j=2*i
        out = (p[j-1]-p[j+1]) *FΓ( x-p[j+2]; p=ρ, IND=0 )
    end
    return out
end

"""
    modelm_j(t,p,ρ)

Compute the Jacobian for modelm(t,p,ρ) with respect to parameters p,
    in the form:
    ``p = \\{ r_0, t_0, r_1, t_1, \\ldots , r_m, t_m \\}``, thus
    generating a  ``n_t \\, \\times \\, 2m+1`` matrix, where ``n_t``
    is the number of elements in input argument ``t``.

# Sources:
1. private notes contact me via https://github.com/f-hipolito

# Examples
```julia-repl
julia> x1 = [0.0, 50.0, 200.0]

julia> p0 = [0.25, 0. ];

julia> ρ0=[6.2,.97];

julia> modelm_j(x1, p0, ρ0)
3×2 Array{Float64,2}:
 0.0          -0.25
 0.000186331  -3.72663e-6
 3.8575e-20   -1.92875e-22

julia> modelm_j(x1, p1, ρ0)
3×4 Array{Float64,2}:
 0.0        -6.84327e-9  0.0          0.0
 0.0846931  -0.0017074   0.000676733  9.74565e-5
 0.282196   -0.00549108  0.81602      0.000960939

julia> modelm_j(x1, p2, ρ0)
3×6 Array{Float64,2}:
 0.0        -6.84327e-9   0.0          0.0         0.0        -0.0
 0.0846931  -0.0017074    0.000676733  9.74565e-5  0.0        -0.0
 0.0166569  -0.000324115  0.011344     5.67202e-5  0.0368222  -8.10288e-6
```
"""
function modelm_j(x,p,ρ)
    #compute m and verify that length of p is even
    lp::Int = length(p)
    lx::Int = length(x)
    m::Int,r::Int = fldmod( lp, 2) .- (1,0)
    # @assert r==0 "Length of p must be even"
    # @assert iseven( length( ρ ) ) "Length of ρ must be even"

    J = zeros(Float64, lx, lp )
    R  = zeros(Float64, lx )
    dn = zeros(Float64, lx )

    j::Int = 0

    Rt!(R,m,x,p,ρ)

    dn = σ.(R) .* σ.(.-R)

    for i = 0:m
        j=2*i
        J[:,j+1] = dn .* dr(i,m,x,p,ρ)
        J[:,j+2] = dn .* dt(i,x,p,ρ)
    end

    return J
end

"""
    d2rr( k,i,m,x,p,ρ )

Evaluates the k,i-th 2nd order partial derivative of the argument of
the m-rate logistic function, with respect to the 'rate' parameters.

# Sources:
1. private notes contact me via https://github.com/f-hipolito

# Examples
```julia-repl
julia> d2rr(0,0,0,-10,[0.25, 0. ],[6.2,.97])
0.0

julia> d2rr(0,0,0,[-10, 0.0, 10.0],[0.25, 0. ],[6.2,.97])
3-element Array{Float64,1}:
 0.0
 0.0
 0.0
```
"""
function d2rr( k,i,m,x::Array,p,ρ )
    out = zeros(Float64, length(x) )
    # all are identically zero in the current model
    if k == 0 && i == 0
        #@. out = 0
    elseif k != 0 && i == 0
        #@. out = 0
    elseif k == 0 && i != 0
        #@. out = 0
    elseif k != 0 && i != 0
        #@. out = 0
    end
    return out
end
function d2rr( k,i,m,x::Number,p,ρ )
    out::Float64 = 0
    # all are identically zero in the current model
    if k == 0 && i == 0
        #@. out = 0
    elseif k != 0 && i == 0
        #@. out = 0
    elseif k == 0 && i != 0
        #@. out = 0
    elseif k != 0 && i != 0
        #@. out = 0
    end
    return out
end

"""
    d2tr( k,i,m,x,p,ρ )

Evaluates the k,i-th 2nd order partial derivative of the argument of
the m-rate logistic function, with respect to 'time' and 'rate' parameters.
Note, the partial derivative is defined as follows:
``\\partial^2 R / \\partial t_k \\partial r_i \\equiv
\\partial / \\partial t_k  \\big( \\partial R / \\partial r_i \\big)``

# Sources:
1. private notes contact me via https://github.com/f-hipolito

# Examples
```julia-repl
julia> d2tr( 1,1,1,200.0,[0.25, 18.8, 0.075, 45. ],[6.2,.97] )
-1.0

julia> d2tr( 1,1,1,[0.0, 50.0, 200.0],[0.25, 18.8, 0.075, 45. ],[6.2,.97] )
3-element Array{Float64,1}:
 -0.0
 -0.3261656196046895
 -1.0
```
"""
function d2tr( k,i,m,x::Array,p,ρ )
    out = zeros(Float64, length(x) )
    if k == 0 && i == 0
        #@. out = 0     #no instruction required for the present model
    elseif k != 0 && i == 0
        #@. out = 0     #no instruction required for the present model
    elseif k == 0 && i != 0
        #@. out = 0     #no instruction required for the present model
    elseif k != 0 && i != 0  # ∂^2 R / ( ∂t_k ∂r_i )
        if k == i       # (-1) δ_{k,i}
            @. out = -FΓ( x -p[2*k+2];p=ρ, IND=0 )
        elseif k == i+1 # δ_{k,i+1}
            @. out = FΓ( x -p[2*k+2];p=ρ, IND=0 )
        end
    end
    return out
end
function d2tr( k,i,m,x::Number,p,ρ )
    out::Float64 = 0
    if k == 0 && i == 0
        #@. out = 0     #no instruction required for the present model
    elseif k != 0 && i == 0
        #@. out = 0     #no instruction required for the present model
    elseif k == 0 && i != 0
        #@. out = 0     #no instruction required for the present model
    elseif k != 0 && i != 0  # ∂^2 R / ( ∂t_k ∂r_i )
        if k == i       # (-1) δ_{k,i}
            out = -FΓ( x -p[2*k+2];p=ρ, IND=0 )
        elseif k == i+1 # δ_{k,i+1}
            out = FΓ( x -p[2*k+2];p=ρ, IND=0 )
        end
    end
    return out
end

"""
    d2rt( k,i,m,x,p,ρ )

Evaluates the k,i-th 2nd order partial derivative of the argument of
the m-rate logistic function, with respect to 'rate' and 'time' parameters.
Note, the partial derivative is defined as follows:
``\\partial^2 R / \\partial r_k \\partial t_i \\equiv
\\partial / \\partial r_k  \\big( \\partial R / \\partial t_i \\big)``

# Sources:
1. private notes contact me via https://github.com/f-hipolito

# Examples
```julia-repl
julia> d2rt( 1,1,1,200,[0.25, 18.8, 0.075, 45. ],[6.2,.97] )
-1.0

julia> d2rt( 1,1,1,[0.0, 50.0, 200.0],[0.25, 18.8, 0.075, 45. ],[6.2,.97] )
3-element Array{Float64,1}:
 -0.0
 -0.3261656196046895
 -1.0
```
"""
function d2rt( k,i,m,x::Array,p,ρ )
    out = zeros(Float64, length(x) )
    if k == 0 && i == 0
        #@. out = 0     #no instruction required for the present model
    elseif k != 0 && i == 0
        #@. out = 0     #no instruction required for the present model
    elseif k == 0 && i != 0     # ∂^2 R / ( ∂r_0 ∂t_i )
        if i==1
            @. out = FΓ( x -p[2*i+2] ; p=ρ, IND=0)
        end
    elseif k != 0 && i != 0     # ∂^2 R / ( ∂r_k ∂t_i )
        if k == i       # (-1) δ_{k,i}
            @. out = -FΓ( x -p[2*i+2] ; p=ρ, IND=0)
        elseif k+1 == i # δ_{k+1,i}
            @. out = FΓ( x -p[2*i+2] ; p=ρ, IND=0)
        end
    end
    return out
end
function d2rt( k,i,m,x::Number,p,ρ )
    out::Float64 = 0
    if k == 0 && i == 0
        #@. out = 0     #no instruction required for the present model
    elseif k != 0 && i == 0
        #@. out = 0     #no instruction required for the present model
    elseif k == 0 && i != 0     # ∂^2 R / ( ∂r_0 ∂t_i )
        if i==1
            out = FΓ( x -p[2*i+2] ; p=ρ, IND=0)
        end
    elseif k != 0 && i != 0     # ∂^2 R / ( ∂r_k ∂t_i )
        if k == i       # (-1) δ_{k,i}
            out = -FΓ( x -p[2*i+2] ; p=ρ, IND=0)
        elseif k+1 == i # δ_{k+1,i}
            out = FΓ( x -p[2*i+2] ; p=ρ, IND=0)
        end
    end
    return out
end

"""
    d2tt( k,i,m,x,p,ρ )

Evaluates the k,i-th 2nd order partial derivative of the argument of
the m-rate logistic function, with respect to the 'time' parameters.

# Sources:
1. private notes contact me via https://github.com/f-hipolito

# Examples
```julia-repl
julia> d2tt( 1,1,1,50,[0.25, 18.8, 0.075, 45. ],[6.2,.97] )
-0.028867524522792173

julia> d2tt( 1,1,1,[0.0, 50.0, 200.0],[0.25, 18.8, 0.075, 45. ],[6.2,.97] )
3-element Array{Float64,1}:
 -0.0
 -0.028867524522792173
 -1.0608237504005124e-57
```
"""
function d2tt( k,i,m,x::Number,p,ρ )
    out::Float64 = 0
    if k == 0 && i == 0
        #@. out = 0     #no instruction required for the present model
    elseif k != 0 && i == 0
        #@. out = 0     #no instruction required for the present model
    elseif k == 0 && i != 0
        #@. out =       #no instruction required for the present model
    elseif k != 0 && i == k     # ∂^2 R / ( ∂t_k ∂t_i )
        j = 2*i
        out = ( p[j+1] -p[j-1] ) * fΓ( x -p[j+2]; p=ρ )
    end
    return out
end
function d2tt( k,i,m,x::Array,p,ρ )
    out = zeros(Float64, length(x) )
    if k == 0 && i == 0
        #@. out = 0     #no instruction required for the present model
    elseif k != 0 && i == 0
        #@. out = 0     #no instruction required for the present model
    elseif k == 0 && i != 0
        #@. out =       #no instruction required for the present model
    elseif k != 0 && i == k     # ∂^2 R / ( ∂t_k ∂t_i )
        j = 2*i
        @. out = ( p[j+1] -p[j-1] ) * fΓ( x -p[j+2]; p=ρ )
    end
    return out
end

"""
    modelm_h(t,p,ρ)

Compute the Hessian for modelm(t,p,ρ).
Compute the Hessian for modelm(t,p,ρ) with respect to parameters p,
    in the form:
    ``p = \\{ r_0, t_0, r_1, t_1, \\ldots , r_m, t_m \\}``, thus
    generating a ``n_t \\, \\times \\, 2m+1 \\, \\times \\, 2m+1``
    matrix, where ``n_t`` is the number of elements in input
    argument ``t``.

# Sources:
1. private notes contact me via https://github.com/f-hipolito

# Examples
```julia-repl
julia> modelm_h( 10, [0.25, 0. ], [6.2,.97] )
2×2 Array{Float64,2}:
-5.94678    0.594678
 0.594678  -0.0594678

julia> modelm_h( [-10, 0.0, 10.0], [0.25, 0. ], [6.2,.97] )
3×2×2 Array{Float64,3}:
[:, :, 1] =
  5.94678  0.594678
  0.0      0.0
 -5.94678  0.594678

[:, :, 2] =
 0.594678   0.0594678
 0.0        0.0
 0.594678  -0.0594678
```
"""
function modelm_h(x::Array,p,ρ)
    #compute m and verify that length of p is even
    lp::Int = length(p)
    lx::Int = length(x)
    m::Int,r::Int = fldmod( lp, 2) .-(1,0)
    # @assert r==0 "Length of p must be even"
    # @assert iseven( length( ρ ) ) "Length of ρ must be even"

    H   = zeros(Float64, lx, lp, lp )
    d2R = zeros(Float64, lx, lp, lp )
    dR = zeros(Float64, lx, lp )
    R  = zeros(Float64, lx )
    dn = zeros(Float64, lx )
    δn = zeros(Float64, lx )
    j::Int=0

    # evaluate R and all terms involving the Logistic function
    Rt!(R,m,x,p,ρ)
    dn = σ.(  R) .* σ.(.-R)
    δn = σ.(.-R) .- σ.(  R)

    # evaluate the first and second order derivatives
    for i=0:m
        j=2*i
        dR[:,j+1] = dr(i,m,x,p,ρ)
        dR[:,j+2] = dt(i,x,p,ρ)
        for k=0:m
            l=2*k
            # ∂^2 R / ∂r^2   # all are identically zero in the current model
            # d2R[:,l+1,j+1] = d2rr( k,i,m,x,p,ρ )

            # ∂^2 R / ( ∂t ∂r )
            d2R[:,l+2,j+1] = d2tr( k,i,m,x,p,ρ )

            # ∂^2 R / ( ∂t ∂r )
            d2R[:,l+1,j+2] = d2rt( k,i,m,x,p,ρ )

            # ∂^2 R / ∂r^2
            d2R[:,l+2,j+2] = d2tt( k,i,m,x,p,ρ )
        end
    end

    for j=1:lp
        for l=1:lp
            @. H[:,l,j] = dn *( d2R[:,l,j] +δn*dR[:,l]*dR[:,j] )
        end
    end

    return H
end

function modelm_h(x::Number,p,ρ)
    #compute m and verify that length of p is even
    lp::Int = length(p)
    m::Int,r::Int = fldmod( lp, 2) .-(1,0)
    # @assert r==0 "Length of p must be even"
    # @assert iseven( length( ρ ) ) "Length of ρ must be even"

    H   = zeros(Float64, lp, lp )
    d2R = zeros(Float64, lp, lp )
    dR = zeros(Float64, lp )
    R::Float64  = 0
    dn::Float64 = 0
    δn::Float64 = 0
    j::Int=0

    # evaluate R and all terms involving the Logistic function
    R = Rt(m,x,p,ρ)
    dn = σ( R) * σ(-R)
    δn = σ(-R) - σ( R)

    # evaluate the first and second order derivatives
    for i=0:m
        j=2*i
        dR[j+1] = dr(i,m,x,p,ρ)
        dR[j+2] = dt(i,x,p,ρ)
        for k=0:m
            l=2*k
            # ∂^2 R / ∂r^2   # all are identically zero in the current model
            # d2R[:,l+1,j+1] = d2rr( k,i,m,x,p,ρ )

            # ∂^2 R / ( ∂t ∂r )
            d2R[l+2,j+1] = d2tr( k,i,m,x,p,ρ )

            # ∂^2 R / ( ∂t ∂r )
            d2R[l+1,j+2] = d2rt( k,i,m,x,p,ρ )

            # ∂^2 R / ∂r^2
            d2R[l+2,j+2] = d2tt( k,i,m,x,p,ρ )
        end
    end

    for j=1:lp
        for l=1:lp
            H[l,j] = dn *( d2R[l,j] +δn*dR[l]*dR[j] )
        end
    end

    return H
end

x0 = [-10, 0.0, 10.0]
x1 = [0.0, 50.0, 200.0]

p0 = [0.25, 0. ]
p1 = [0.25, 18.8, 0.075, 45. ]
p2 = [0.25, 18.8, 0.075, 45., 0.1, 80 ]

ρ0 = [6.2,.97]

R0 = zeros(Float64, length(x0))
R1 = zeros(Float64, length(x0))
R2 = zeros(Float64, length(x0))

@assert Rt!(R0,0,x0,p0,ρ0) == [ -2.5, 0.0, 2.5 ]
@assert Rt!(R1,1,x0,p1,ρ0) == [ -21.3, -18.8, -16.3 ]
@assert Rt!(R2,2,x0,p2,ρ0) == [ -21.3, -18.8, -16.3 ]

@assert dr( 1,1,x1,p1,ρ0 ) == [ 0.0, 0.3963534314129777, 148.60824742268042 ]
@assert dt( 1,x1,p2,ρ0 ) == [ 0.0, 0.05707898343082066, 0.175 ]

@assert d2rr( 0,0,0,x0,p0,ρ0 ) == zeros(Float64,length(x0))
@assert d2tr( 1,1,1,x1,p1,ρ0 ) == [ -0.0, -0.3261656196046895, -1.0 ]
@assert d2rt( 1,1,1,x1,p1,ρ0 ) == [ -0.0, -0.3261656196046895, -1.0 ]
@assert d2tt( 1,1,1,x1,p1,ρ0 ) == [ -0.0, -0.028867524522792173, -1.0608237504005124e-57 ]

modelm_j(x1, p2, ρ0);
modelm_h( x1, p2, ρ0 );

"""
    sir_Γ!(du,u,t; p, ν=0.01, ρ=[2.0,1.0])

Computes the differential equations for susceptible, infected and removed
    fractions of the population, using a Gamma distribution for the incubation
    times.
Computes the differentials ``du = [ ds, di, dr ]`` for ``u = [ 1-s, i, r ]``
    evaluated at time ``t`` and parametrized via
    ``p_0 = \\{ \\omega_0, \\omega_1, t_1, \\ldots \\omega_m, t_m \\}`` where
    ``\\omega_i`` and ``t_i`` are the ``i^\\mathrm{th}`` growth rate and
    transition times, respectively. ``\\nu`` is the combined recovery/death
    decay-like time scale.
    The distribution of incubation times is parametrized with default values
    ``\\rho_0 = [ \\alpha_0, \\beta_0 ] = [2,1]``.

Note: we normalized variables with respect to the total population ``N_t``, i.e.
    population densities, namely
1. ``s+i+r = 1``
2. ``s = S/N_t`` susceptible individuals
3. ``i = I/N_t`` infected    individuals
4. ``r = R/N_t`` removed     individuals

Furthermore, the equations are written considering for ``1-s`` rather than
    ``s``, to facilitate the optimization process, as at the initial instant
    ``s`` is, typically, ``N_t \\sim 10^{6-8}`` times larger than ``i`` and
    ``r``.

# Sources:
1. see https://github.com/f-hipolito/PopulationDynamics/blob/master/Population_dynamics.ipynb

# Examples
```julia-repl
julia> du0 = zeros(3); u0 = [ 1e-7, 1e-7, 0. ]; t0 = 1.0;

julia> p0 = [ 0.26, 0.01, 34.9 ]; νμ0 = 1/21.0; ρ0 = [ 6.25, 0.96 ];

julia> sir_Γ!(du0,u0,t0; p=p0, ν=νμ0, ρ=ρ0)

julia> @assert isapprox( du0, [2.59999974e-8,2.1238092638e-8,4.761904761e-9] )
```
"""
function sir_Γ!(du,u,t; p, ν=0.01, ρ=[2.0,1.0])

    # basic testing of parameters
    #
    len = length(p) +1

    if len < 4
        error( "insufficient parameters for sir_Γ!" )
    elseif isodd( len )
        println( len )
        error( "requires even number of parameters!" )
    end

    # loading parameters
    #
    #   the growth rates ωᵢ for the m-steady state model
    m = fld( len, 2) -1     # number m-steady states
    ω = zeros( len )
    ω[1] = p[1]
    ω[3:end] = p[2:end]

    #   the incubation time paramters
    ρ0 = ρ # ρ₀ = ρ = [ α₀, β₀ ]

    #   compute the "m-steady state" rates
    #
    Ω = ω[1]
    for i=1:m
        j=2*i
        Ω += (ω[j+1]-ω[j-1])*FΓ(t-ω[j+2];p=ρ0)
    end

    #   the differential equations
    #
    dr = ν *u[2]
    ds = Ω *u[2] *(1.0 -u[1])
    du[1] =  ds                 #   u[1] ≣ 1 -susceptible
    du[2] =  di = ds -dr        #   u[2] ≣ infected
    du[3] =  dr                 #   u[3] ≣ removed
end

du0 = zeros(3); u0 = [ 1e-7, 1e-7, 0. ]; t0 = 1.0;

p0 = [ 0.26, 0.01, 34.9 ]; νμ0 = 1/21.0; ρ0 = [ 6.25, 0.96 ];

sir_Γ!(du0,u0,t0; p=p0, ν=νμ0, ρ=ρ0)

@assert isapprox( du0, [ 2.59999974e-8, 2.1238092638e-8, 4.761904761e-9 ] )

"""
    sird_Γ!(du,u,t; p, νμ=[1e-2,1e-4], ρ=[2.0,1.0] )

Computes the differential equations for susceptible, infected, recovered and
    deceased fractions of the population, using a Gamma distribution for the
    incubation times.
Computes the differentials ``du = [ ds, di, dr, dd ]`` for
    ``u = [ 1-s, i, r, d ]`` evaluated at time ``t`` and parametrized via
    ``p = \\{ \\omega_0, \\omega_1, t_1, \\ldots \\omega_m, t_m \\}``
    where ``\\omega_i`` and ``t_i`` are the ``i^\\mathrm{th}`` growth rate and
    transition times, respectively. ``\\nu,\\, \\mu`` are the recovery and death
    decay-like time scales.
    The distribution of incubation times is parametrized with default values
    ``\\rho_0 = [ \\alpha_0, \\beta_0 ] = [2,1]``.

Note: we normalized variables with respect to the total population ``N_t``, i.e.
    population densities, namely
1. ``s+i+r = 1``
2. ``s = S/N_t`` susceptible individuals
3. ``i = I/N_t`` infected    individuals
4. ``r = R/N_t`` recovered   individuals
5. ``d = D/N_t`` deceased    individuals

Furthermore, the equations are written considering for ``1-s`` rather than
    ``s``, to facilitate the optimization process, as at the initial instant
    ``s`` is, typically, ``N_t \\sim 10^{6-8}`` times larger than ``i`` and
    ``r``.

# Sources:
1. see https://github.com/f-hipolito/PopulationDynamics/blob/master/Population_dynamics.ipynb

# Examples
```julia-repl
julia> du0 = zeros(4); u0 = [ 1e-7,1e-7,0.,0. ]; t0 = 1.0;

julia> p0 = [ 0.26, 0.01, 34.9 ]; νμ0 = [ 1/21.0, 0.5/21 ]; ρ0 = [ 6.25, 0.96 ];

julia> sird_Γ!(du0,u0,t0; p=p0, νμ=νμ0, ρ=ρ0)

julia> @assert isapprox(du0,[2.59999974e-8,1.885714025e-8,4.76190476e-9,2.38095238e-9])
```
"""
function sird_Γ!(du,u,t; p, νμ = [ 1e-2, 1e-4 ], ρ=[2.0,1.0])
        # basic testing of parameters
    #
    len = length(p) +1

    if len < 4
        error( "insufficient parameters for sird_Γ!" )
    elseif isodd( len )
        println( len )
        error( "requires even number of p parameters!" )
    end

    # loading parameters
    #
    #   the growth rates ωᵢ for the m-steady state model
    m = fld( len, 2) -1     # number m-steady states
    ω = zeros( len )
    ω[1] = p[1]
    ω[3:end] = p[2:end]

    #   the incubation time paramters
    ρ0 = ρ # ρ₀ = ρ = [ α₀, β₀ ]

    #   the infection decay time
    ν = νμ[1]
    μ = νμ[2]

    #   compute the "m-steady state" rates
    #
    Ω = ω[1]
    for i=1:m
        j=2*i
        Ω += (ω[j+1]-ω[j-1])*FΓ(t-ω[j+2];p=ρ0)
    end

    #   the differential equations
    #
    dr = ν *u[2]
    dd = μ *u[2]
    ds = Ω *u[2] *( 1.0 -u[1] )
    du[1] =  ds                 #   u[1] ≣ 1 -susceptible
    du[2] =  di = ds -dr -dd    #   u[2] ≣ infected
    du[3] =  dr                 #   u[3] ≣ recovered
    du[4] =  dd                 #   u[3] ≣ deceased
end

du0 = zeros(4); u0 = [ 1e-7, 1e-7, 0., 0. ]; t0 = 1.0;

p0 = [ 0.26, 0.01, 34.9 ]; νμ0 = [ 1/21.0, 0.5/21 ]; ρ0 = [ 6.25, 0.96 ];

sird_Γ!(du0,u0,t0; p=p0, νμ=νμ0, ρ=ρ0)

@assert isapprox(du0,[2.59999974e-8,1.885714025e-8,4.76190476e-9,2.38095238e-9])

end

max_stars_repo_name: f-hipolito/PopulationDynamics

max_issues_repo_name: f-hipolito/PopulationDynamics

max_forks_repo_name: f-hipolito/PopulationDynamics


text:
#=
This is a comment block
=#

# This is an inline comment

# Variable declaration
iVar1 = 1
iVar2 = 2

# Variable arithmatic
iVar3 = iVar1+iVar2
iVar4 = iVar1*iVar2

# Print tests
println(iVar3)
println(iVar4)
println()

# Loop format and test
for int = 1:10
    println(int)
end

println()

iVar5 = 1
while iVar5 <= 10
    println(iVar5)
    # Must include 'global' for variables in loops
    global iVar5 += 1
end