text:
@testset "LTI systems tests" begin

    s = Romeo.LTI.Diff(1)
    G = (s + 1) / (s + 2)
    zerosys = Romeo.LTI.ZeroSys()
    idsys = Romeo.LTI.UnitSys()

    @testset "Diff tests" begin
        @test s.α == 1
        @test Romeo.LTI.Diff(0) == idsys
    end

    @testset "Unary operators tests" begin
        @test +s == s
        @test -s == Romeo.LTI.ScaledSystem(-1, s)
        @test -(-s) == s
    end

    @testset "Series connection tests" begin
        for sys in [s, G]
            @test 1 * sys == sys * 1 == sys
            @test 1.0 * sys == sys * 1.0 == sys
            @test 2 * sys == 2.0 * sys == sys * 2.0 == sys * 2
            @test 0 * sys == sys * 0 == 0.0 * sys == sys * 0.0 == zerosys
            @test -0 * sys == -sys * 0 == -0.0 * sys == -sys * 0.0 == zerosys
            @test idsys * sys == sys * idsys == sys
            @test zerosys * sys == sys * zerosys == zerosys
            @test 2 * sys == sys * 2 == Romeo.LTI.ScaledSystem(2, sys)
            @test 4 * (sys * 2) == Romeo.LTI.ScaledSystem(8, sys)
            @test -(8 * sys) == Romeo.LTI.ScaledSystem(-8, sys)
            @test (2 * sys) * 0.5 == sys
            @test (2 * sys) * (4 * sys) == 8 * (sys * sys)
        end
        @test (2 * s) * G == s * (2 * G) == 2 * (s * G)
    end

    @testset "Parallel connection tests" begin
        for sys in [s, G]
            @test sys + sys == 2 * sys
            @test 2 * sys + 4 * sys == 6 * sys
            @test sys + 2 * sys == 2 * sys + sys == 3 * sys
            @test sys - sys == zerosys
            @test sys + zerosys == zerosys + sys == sys
            @test sys + 0 == 0 + sys == sys + 0.0 == 0.0 + sys == sys
        end
    end

    @testset "Power system tests" begin
        orders = [0.5, 2, 2.5, 3]
        @test idsys^0 == idsys^0.0 == idsys
        @test zerosys^0 == zerosys^0.0 == idsys
        @test sqrt(idsys) == idsys
        @test sqrt(zerosys) == zerosys
        @test sqrt(s^0.5) == s^0.25
        for α in orders
            @test idsys^α == idsys
            @test zerosys^α == zerosys
        end
        for sys in [s, G]
            @test sys^0 == sys^0.0 == idsys
            @test sys^1 == sys^1.0 == sys
            @test sys^0.5 == sqrt(sys)
            @test sqrt(4 * sys) == 2 * sqrt(sys)
            for α in orders
                @test (2 * sys)^α == 2^α * sys^α
            end
        end
    end

    @testset "Rational systems tests" begin
        for sys in [s, G]
            @test sys / 1 == sys / 1.0 == sys == 1 \ sys == 1.0 \ sys
            @test sys / idsys == sys == idsys \ sys
            # @test_throws ErrorException sys / zerosys
            # @test_throws ErrorException zerosys \ sys
            @test sys / sys == idsys
            @test sys^(15 // 100) / sys^(1 // 10) == sys^(5 // 100)
            @test sys / sys^0.4 == sys^0.6
            @test sys^(12 // 10) / sys == sys^(2 // 10)
        end
        @test s / (4 * G) == 1 / 4 * s / G
        @test (6 * s) / G == 6 * (s / G)
    end

end # testset


text:
export getData
import jInv.ForwardShare.getData
function getData(m,pFor::ConvDiffFEMParam)
	pFor.Ainv.doClear = 0 # PDE does not depend on m, so we shouldn't clear i

	# solve PDE for non-Dirichlet nodes
	Ut,pFor.Ainv      = solveLinearSystem(pFor.A, pFor.MassConst*m - pFor.rbc, pFor.Ainv)

	U                 = pFor.ubc
	U[pFor.idint]     = Ut
	dobs 		      = pFor.P*U
	pFor.Fields       = U

	return dobs,pFor
end

import jInv.ForwardShare.getSensMatVec
function getSensMatVec(v::Vector{Float64},m::Vector{Float64},pFor::ConvDiffFEMParam)
	dUt = solveLinearSystem(pFor.A, pFor.MassConst*v, pFor.Ainv)[1]
	dU             = zeros(prod(pFor.M.n.+1),1)
	dU[pFor.idint] = dUt
	return  pFor.P*dU
end

import jInv.ForwardShare.getSensTMatVec
function getSensTMatVec(v::Vector{Float64},m::Vector{Float64},pFor::ConvDiffFEMParam)
	return  pFor.MassConst'*solveLinearSystem(pFor.A,(pFor.P'*v)[pFor.idint], pFor.Ainv,1)[1]
end

ext: jl

max_stars_repo_name: JuliaInv/ConvDiffMIPDECO

max_issues_repo_name: JuliaInv/ConvDiffMIPDECO

max_forks_repo_name: JuliaInv/ConvDiffMIPDECO


text:
module AbstractDAQs

export AbstractDAQ, AbstractPressureScanner
export DAQTask, isreading, samplesread, issamplesavailable
export isdaqfinished
export stoptask, stoptask!, cleartask!
export samplingrate, settiming!
export setdaqthread!, daqthread
export setdaqtask!, daqtask
export daqaddinput, daqacquire, daqacquire!, daqstart
export daqread, daqread!, daqstop, devname, devtype

export daqreference, daqzero, daqconfig, daqconfigdev
export numchannels, daqchannels
export CircMatBuffer, bufwidth
export nextbuffer, isfull, isempty, flatten, flatten!, capacity
export DAQConfig, iparameters, fparameters, sparameters
export daqdevip, daqdevmodel, daqdevserialnum, daqdevtag
export savedaqdata, savedaqconfig, readdaqdata
export TestDev

using Dates

import DataStructures: OrderedDict
"Abstract type to handle any kind of device"
abstract type AbstractDevice end

"Abstract type to handle data acquisition devices"
abstract type AbstractDAQ <: AbstractDevice end

"Abstract type to handle pressure scanners"
abstract type AbstractPressureScanner <: AbstractDAQ end

"Returns the type of device"
devtype(dev) = string(typeof(dev))

"""
`devname(dev::AbstractDevice)`

The device name is a string that is used to refer to a specific device.

This string is used when saving data and post processing.
"""
devname(dev::AbstractDevice) = dev.devname

include("utils.jl")
include("daqconfig.jl")
include("daqtask.jl")
include("circbuffer.jl")
include("measdata.jl")
include("hdf5io.jl")
include("interface.jl")
include("deviceset.jl")
include("testdevice.jl")
end


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

create_org(name::String, address::String, zip::String, city::String, country::String) =
    Organization(
        name::String,
        address::String,
        zip::String,
        city::String,
        country::String)

create_training(id::String, date::DateTime, days::Int64, name::String, price::Float64) =
    Training(id, date, days, name, price)::Training

create_order(org::Organization, training::Training, order_ref::String, contact_name::String, contact_email::String, students::Array{String,1}) =
    Order(org, training, order_ref, contact_name, contact_email, students)::Order


text:
using Dates
using Printf

using .Networks

###############################################
# Sequence Generation
###############################################

struct SequenceGenerator <: AbstractDataGenerator
    network::Networks.Network #Not used for now (but potentially we can have specific operations at DMO horizons)
    target_timepoints::Vector{Dates.DateTime}
    horizon_timepoints::Vector{Dates.DateTime}
    management_mode::ManagementMode
end

function launch(seq_generator::SequenceGenerator)
    if seq_generator.management_mode == PSCOPF_MODE_1
        return gen_seq_mode1(seq_generator)

    elseif seq_generator.management_mode == PSCOPF_MODE_2
        return gen_seq_mode2(seq_generator)

    elseif seq_generator.management_mode == PSCOPF_MODE_3
        return gen_seq_mode3(seq_generator)
    end

    error("unsuppported mode : ", seq_generator.management_mode)
end

"""
    generate_sequence
"""
function generate_sequence(network::Networks.Network, target_timepoints::Vector{Dates.DateTime},
                            horizon_timepoints::Vector{Dates.DateTime}, management_mode::ManagementMode)
    generator = SequenceGenerator(network, target_timepoints, horizon_timepoints, management_mode)
    return launch(generator)
end

function gen_seq_mode1(seq_generator::SequenceGenerator)
    sequence = Sequence()
    first_ts = seq_generator.target_timepoints[1]
    fo_startpoint = first_ts - get_fo_length(seq_generator.management_mode)

    for ech in seq_generator.horizon_timepoints
        if ech <  first_ts
            if ech < fo_startpoint
                add_step!(sequence, EnergyMarket, ech)
                add_step!(sequence, TSOOutFO, ech)

            elseif ech == fo_startpoint
                add_step!(sequence, EnergyMarketAtFO, ech)
                add_step!(sequence, EnterFO, ech)
                add_step!(sequence, TSOInFO, ech)

            else
                add_step!(sequence, TSOInFO, ech)
            end
        elseif first_ts < ech
            msg = @sprintf(("Error when generating sequence: ech (%s) is after target timepoint (%s)."),
                            ech, first_ts)
            throw( error(msg) )
        end
    end

    add_step!(sequence, Assessment, seq_generator.horizon_timepoints[end])

    return sequence
end

function gen_seq_mode2(seq_generator::SequenceGenerator)
    sequence = Sequence()
    first_ts = seq_generator.target_timepoints[1]
    fo_startpoint = first_ts - get_fo_length(seq_generator.management_mode)

    preceding_ech = nothing
    for ech in seq_generator.horizon_timepoints
        if ech <  first_ts
            if ech < fo_startpoint
                add_step!(sequence, EnergyMarket, ech)
                add_step!(sequence, TSOOutFO, ech)

            elseif ech == fo_startpoint
                add_step!(sequence, EnergyMarketAtFO, ech)
                add_step!(sequence, EnterFO, ech)
                add_step!(sequence, TSOOutFO, ech)

            elseif preceding_ech == fo_startpoint
                add_step!(sequence, EnergyMarket(EnergyMarketConfigs(REF_SCHEDULE_TYPE=TSO())), ech)
                add_step!(sequence, TSOBilevel(TSOBilevelConfigs(REF_SCHEDULE_TYPE_IN_TSO=TSO())), ech) #Ref can be market cause it respects the preceding TSO
                add_step!(sequence, BalanceMarket(EnergyMarketConfigs(REF_SCHEDULE_TYPE=TSO())), ech)
            else
                add_step!(sequence, EnergyMarket(EnergyMarketConfigs(REF_SCHEDULE_TYPE=Market())), ech)
                # NOTE: The TSOBilevel considers the EnergyMarket both as a reference for deltas and for decided schedules
                #we can keep this cause the EnergyMarket respects the decisions of the previous BalanceMarket
                add_step!(sequence, TSOBilevel(TSOBilevelConfigs(REF_SCHEDULE_TYPE_IN_TSO=Market())), ech)
                add_step!(sequence, BalanceMarket(EnergyMarketConfigs(REF_SCHEDULE_TYPE=Market())), ech)
            end
        elseif first_ts < ech
            msg = @sprintf(("Error when generating sequence: ech (%s) is after target timepoint (%s)."),
                            ech, first_ts)
            throw( error(msg) )
        end
        preceding_ech = ech
    end

    add_step!(sequence, Assessment, seq_generator.horizon_timepoints[end])

    return sequence
end

function gen_seq_mode3(seq_generator::SequenceGenerator)
    sequence = Sequence()
    first_ts = seq_generator.target_timepoints[1]
    fo_startpoint = first_ts - get_fo_length(seq_generator.management_mode)

    for ech in seq_generator.horizon_timepoints
        if ech <  first_ts
            if ech < fo_startpoint
                add_step!(sequence, EnergyMarket, ech)
                add_step!(sequence, TSOOutFO, ech)

            elseif ech == fo_startpoint
                add_step!(sequence, EnergyMarket, ech)
                add_step!(sequence, EnterFO, ech)
                add_step!(sequence, TSOBilevel(TSOBilevelConfigs(REF_SCHEDULE_TYPE_IN_TSO=TSO())), ech)

            else
                #TODO : check if this is an EnergyMarket or BalanceMarket or another implem
                add_step!(sequence, BalanceMarket, ech)
            end
        elseif first_ts < ech
            msg = @sprintf(("Error when generating sequence: ech (%s) is after target timepoint (%s)."),
                            ech, first_ts)
            throw( error(msg) )
        end
    end

    add_step!(sequence, Assessment, seq_generator.horizon_timepoints[end])

    return sequence
end

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
"""
    genmonomials(dim::Integer, maxdegree:Integer)

Generates an array of monomials.

# Arguments
- `dim::Int`: the amount of variables, aka dimension of the input data,
- `maxdegree::Int`: maximal degree of each term in the polynomial.

Returns the monomial variables and the monomials themselves.
"""
function genmonomials(dim::Integer, maxdegree::Integer)
  vars = @polyvar x[1:dim]
  if dim < 0 || maxdegree < 0
    error("Cannot generate a polynomial with negative dimensions")
  end
  monomials::Vector{Monomial{true}} = [constantterm(1, sum(vars))]

  for v in vars
    tmpMonomials::Vector{Monomial{true}} = []
    for done in monomials
      @inbounds for d in 1:maxdegree
        tmp::Monomial{true} = done * v^d
        if degree(tmp) > maxdegree
          break
        else
          push!(tmpMonomials, tmp)
        end
      end
    end
    for m in tmpMonomials
      push!(monomials, m)
    end
  end

  return vars, monomials
end


text:
#=
  These routines return +/-Inf when given +/-Inf.
  That differs from nextfloat(-Inf) == -realmax(), prevfloat(Inf) == realmax()
  (prevfloat(Inf)==Inf makes more sense to me, and likely is more helpful).
  And steps twice when given values of very small magnitude (see paper).
  The alternative implementation, converting to [U]Int and adding/subtracting 1,
  returns NaN when given +/-Inf -- and checking for Inf adds branching.
  ref:
  'Predecessor and Successor In Rounding To Nearest' by Siegried M. Rump
  "The routines deliver the exact answer except for a small range near underflow,
   in which case the true result is overestimated by eta [the value added/subtracted below]."
=#
# exact for |x| > 8.900295434028806e-308
nextNearerToZero(x::Float64)   = (0.9999999999999999*x)-5.0e-324 # (x-1.1102230246251568e-16*x)-5.0e-324
nextAwayFromZero(x::Float64)   = (1.0000000000000002*x)+5.0e-324 # (x+1.1102230246251568e-16*x)+5.0e-324
# exact for |x| > 4.814825f-35
nextNearerToZero(x::Float32)   = (0.99999994f0*x)-1.435f-42      # (x-5.960465f-8*x)-1.435f-42
nextAwayFromZero(x::Float32)   = (1.00000010f0*x)+1.435f-42      # (x+5.960465f-8*x)+1.435f-42
# the multiplicative formulation for Float16 is exact for |x| > Float16(0.25)
# which is quite coarse, we do not use that here
nextNearerToZero(x::Union{BigFloat,Float16}) = signbit(x) ? nextfloat(x) : prevfloat(x)
nextAwayFromZero(x::Union{BigFloat,Float16}) = signbit(x) ? prevfloat(x) : nextfloat(x)

@inline nextFloat{T<:AbstractFloat}(x::T) = signbit(x) ? nextNearerToZero(x) : nextAwayFromZero(x)
@inline prevFloat{T<:AbstractFloat}(x::T) = signbit(x) ? nextAwayFromZero(x) : nextNearerToZero(x)

@inline nextFloat(x::Union{BigFloat,Float16}) = nextfloat(x)
@inline prevFloat(x::Union{BigFloat,Float16}) = prevfloat(x)

function nextFloat{T<:AbstractFloat}(x::T, n::Int)
    if !signbit(n)
        for i in 1:n
            x = nextFloat(x)
        end
    else
        x = prevFloat(x,-n)
    end
    x
end

function prevFloat{T<:AbstractFloat}(x::T, n::Int)
    if !signbit(n)
        for i in 1:n
            x = prevFloat(x)
        end
    else
        x = nextFloat(x,-n)
    end
  x
end

for (F,U,I) in ((:Float16,:UInt16,:Int16), (:Float32,:UInt32,:Int32), (:Float64,:UInt64,:Int64))
  @eval begin
    function nFloatsSeparate{T<:$F}(a::T, b::T)
      isneg = (a > b)
      a,b = minmax(a,b)

      if signbit(a) == signbit(b)
        if signbit(a)
           z = reinterpret($U,abs(a))-reinterpret($U,abs(b))
        else
           z = reinterpret($U,b)-reinterpret($U,a)
        end
      else
        z = reinterpret($U,b)+reinterpret($U,-a)+one($U)
      end

      if z > reinterpret($U, typemax($I))
        throw(ArgumentError(string("separation exceeds typemax(",$I,")")))
      else
        isneg ? -reinterpret($I,z) : reinterpret($I,z)
      end
    end
  end
end

#= presumes xpa == frexp(a)[2] == frexp(b)[2]
function nFloatsSeparate{T<:Float64}(a::T, b::T,xpa::Int)
    isneg = a > b
    x,y = minmax(a,b)
    z = floor(Int64, ((y-x)/2^xpa)*2^53)
    isneg ? -z : z
end
=#


text:
# push!(LOAD_PATH, joinpath(@__DIR__, "..", ".."))
using JLD2

# jlfile = jldopen("data/local.jld2", "a+")
# jlfile = jldopen("data/nonlocal.jld2", "a+")
jlfile = jldopen("data/nonlocal_symmetric.jld2", "a+")
# jlfile = jldopen("data/nonlocal_more_velocity.jld2")
EF¹¹ = jlfile["diffusivity"]["K11"]
EF¹² = jlfile["diffusivity"]["K12"]
EF²¹ = jlfile["diffusivity"]["K21"]
EF²² = jlfile["diffusivity"]["K22"]
x = jlfile["grid"]["x"]
z = jlfile["grid"]["z"]
N = length(x)
M = length(z)
tmpE = copy(EF¹¹)
EF¹¹ = reshape(permutedims(reshape(EF¹¹, (N, M, N, M)), (2, 1, 4, 3)), (N * M, N * M))
EF¹² = reshape(permutedims(reshape(EF¹², (N, M, N, M)), (2, 1, 4, 3)), (N * M, N * M))
EF²¹ = reshape(permutedims(reshape(EF²¹, (N, M, N, M)), (2, 1, 4, 3)), (N * M, N * M))
EF²² = reshape(permutedims(reshape(EF²², (N, M, N, M)), (2, 1, 4, 3)), (N * M, N * M))
# 2 goes to 1, 1 goes to 2, 4 goes to 3, 3 goes to 4
##

using GLMakie
fig = Figure(resolution = (1800, 1300), title = "Nonlocal Operators")
titlestring = "Kˣˣ"
ax1 = Axis(fig[1, 1], title = titlestring, titlesize = 30)
titlestring = "Kˣᶻ"
ax2 = Axis(fig[1, 3], title = titlestring, titlesize = 30)
titlestring = "Kᶻˣ"
ax3 = Axis(fig[2, 1], title = titlestring, titlesize = 30)
titlestring = "Kᶻᶻ"
ax4 = Axis(fig[2, 3], title = titlestring, titlesize = 30)

colormap = :thermal
colormap2 = :balance
hm1 = heatmap!(ax1, EF¹¹, colormap = colormap)
ax1.yreversed = true

hm2 = heatmap!(ax2, EF¹², colormap = colormap2)
ax2.yreversed = true

hm3 = heatmap!(ax3, EF²¹, colormap = colormap2)
ax3.yreversed = true

hm4 = heatmap!(ax4, EF²², colormap = colormap)
ax4.yreversed = true

Colorbar(fig[1, 2], hm1, height = Relative(3 / 4), width = 25, ticklabelsize = 30,
    labelsize = 30, ticksize = 25, tickalign = 1,)
Colorbar(fig[2, 2], hm2, height = Relative(3 / 4), width = 25, ticklabelsize = 30,
    labelsize = 30, ticksize = 25, tickalign = 1,)
Colorbar(fig[1, 4], hm3, height = Relative(3 / 4), width = 25, ticklabelsize = 30,
    labelsize = 30, ticksize = 25, tickalign = 1,)
Colorbar(fig[2, 4], hm4, height = Relative(3 / 4), width = 25, ticklabelsize = 30,
    labelsize = 30, ticksize = 25, tickalign = 1,)
display(fig)

max_stars_repo_name: sandreza/StatisticalNonlocality

max_issues_repo_name: sandreza/StatisticalNonlocality

max_forks_repo_name: sandreza/StatisticalNonlocality


text:
#=
Tools for matching algorithms.

Author: Daisuke Oyama

=#
using Random

# random_prefs

"""
    random_prefs([rng, ]m, n[, ReturnCaps]; allow_unmatched=true)

Generate random preference order lists for two groups, say, m males and n
females.

Each male has a preference order over femals [1, ..., n] and "unmatched" which
is represented by 0, while each female has a preference order over males
[1, ..., m] and "unmatched" which is again represented by 0.

The argument `ReturnCaps` should be supplied in the context of college
admissions, in which case "males" and "females" should be read as "students"
and "colleges", respectively, where each college has its capacity.

The optional `rng` argument specifies a random number generator.

# Arguments

* `m::Integer` : Number of males.
* `n::Integer` : Number of females.
* `::Type{ReturnCaps}` : If supplied, `caps` is also returned.
* `;allow_unmatched::Bool(true)` : If false, return preference order lists of
  males and females where 0 is always placed in the last place, (i.e.,
  "unmatched" is least preferred by every individual).

# Returns

* `m_prefs::Matrix{Int}` :  Array of shape (n+1, m), where each column contains
  a random permutation of 0, 1, ..., n.
* `f_prefs::Matrix{Int}` :  Array of shape (m+1, n), where each column contains
  a random permutation of 0, 1, ..., m.
* `caps::Vector{Int}` : Vector of length n containing each female's (or
  college's) capacity. Returned only when `ReturnCaps` is supplied.

# Examples

```julia
julia> m_prefs, f_prefs = random_prefs(4, 3);

julia> m_prefs
4x4 Array{Int64,2}:
 3  3  1  3
 0  2  3  1
 2  1  2  0
 1  0  0  2

julia> f_prefs
5x3 Array{Int64,2}:
 1  2  4
 4  3  1
 3  4  2
 0  0  0
 2  1  3

julia> m_prefs, f_prefs = random_prefs(4, 3, allow_unmatched=false);

julia> m_prefs
4x4 Array{Int64,2}:
 1  3  1  2
 2  1  3  3
 3  2  2  1
 0  0  0  0

julia> f_prefs
5x3 Array{Int64,2}:
 1  2  3
 2  3  4
 4  1  1
 3  4  2
 0  0  0

julia> s_prefs, c_prefs, caps = random_prefs(4, 3, ReturnCaps);

julia> s_prefs
4x4 Array{Int64,2}:
 2  1  2  1
 1  3  1  0
 3  2  3  3
 0  0  0  2

julia> c_prefs
5x3 Array{Int64,2}:
 3  4  1
 0  1  4
 4  3  0
 1  2  3
 2  0  2

julia> caps
3-element Array{Int64,1}:
 1
 2
 2
```
"""
function random_prefs(rng::AbstractRNG,
                      m::Integer, n::Integer;
                      allow_unmatched::Bool=true)
    m_prefs = _random_prefs(rng, m, n, allow_unmatched)
    f_prefs = _random_prefs(rng, n, m, allow_unmatched)
    return m_prefs, f_prefs
end

random_prefs(m::Integer, n::Integer; allow_unmatched::Bool=true) =
    random_prefs(Random.GLOBAL_RNG, m, n, allow_unmatched=allow_unmatched)

struct ReturnCaps end

function random_prefs(rng::AbstractRNG,
                      m::Integer, n::Integer, ::Type{ReturnCaps};
                      allow_unmatched::Bool=true)
    s_prefs = _random_prefs(rng, m, n, allow_unmatched)
    c_prefs = _random_prefs(rng, n, m)

    if allow_unmatched
        unmatched_rankings = Array{Int}(undef, n) #rand(r, 2:n+1, m)
        _random_unmatched!(rng, c_prefs, unmatched_rankings)
        caps = _random_caps(rng, unmatched_rankings)
    else
        caps = _random_caps(rng, m, n)
    end
    return s_prefs, c_prefs, caps
end

random_prefs(m::Integer, n::Integer, T::Type{ReturnCaps};
             allow_unmatched::Bool=true) =
    random_prefs(Random.GLOBAL_RNG, m, n, T, allow_unmatched=allow_unmatched)

function _random_prefs(rng::AbstractRNG, m::Integer, n::Integer)
    prefs = Array{Int}(undef, n+1, m)
    for j in 1:m
        prefs[end, j] = 0
    end

    _randperm2d!(rng, view(prefs, 1:n, :))

    return prefs
end

function _random_unmatched!(rng::AbstractRNG, prefs::Matrix{Int},
                            unmatched_rankings::Vector{Int})
    n = size(prefs, 1) - 1
    m = size(prefs, 2)
    rand!(rng, unmatched_rankings, 2:n+1)
    for j in 1:m
        prefs[end, j] = prefs[unmatched_rankings[j], j]
        prefs[unmatched_rankings[j], j] = 0
    end
end

function _random_prefs(rng::AbstractRNG, m::Integer, n::Integer,
                       allow_unmatched::Bool)
    prefs = _random_prefs(rng, m, n)

    if allow_unmatched
        unmatched_rankings = Array{Int}(undef, m)
        _random_unmatched!(rng, prefs, unmatched_rankings)
    end

    return prefs
end

function _random_caps(rng::AbstractRNG, unmatched_rankings::Vector{Int})
    n = length(unmatched_rankings)
    u = rand(rng, n)
    for i in 1:n
        u[i] *= unmatched_rankings[i]
        u[i] += 1
    end
    return floor.(Int, u)
end

function _random_caps(rng::AbstractRNG, m::Int, n::Int)
    u = rand(rng, n)
    for i in 1:n
        u[i] *= m
        u[i] += 1
    end
    return floor.(Int, u)
end

# From stdlib/Random/src/misc.jl
rand_lt(r::AbstractRNG, n::Int, mask::Int=nextpow(2, n)-1) =
    rand(r, Random.ltm52(n, mask))

# In-place version of randperm in base/random.jl
"""
Given a vector `a` of length n, generate a random permutation of length n and
store it in `a`.
"""
function _randperm!(r::AbstractRNG, a::AbstractVector{T}) where T<:Integer
    # a = Array{typeof(n)}(n)
    n = length(a)
    @assert n <= Int64(2)^52
    if n == 0
       return a
    end
    a[1] = 1
    mask = 3
    @inbounds for i = 2:Int(n)
        j = 1 + rand_lt(r, i, mask)
        if i != j # a[i] is uninitialized (and could be #undef)
            a[i] = a[j]
        end
        a[j] = i
        i == 1+mask && (mask = 2mask + 1)
    end
    return a
end

"""
Given an m x n matrix `a`, generate n random permutations of length m and store
them in columns of `a`.
"""
function _randperm2d!(r::AbstractRNG, a::AbstractMatrix{T}) where T<:Integer
    m, n = size(a)
    for j in :1:n
        _randperm!(r, view(a, :, j))
    end
    a
end


text:
using OrdinaryDiffEq, Test, LinearAlgebra

mutable struct SimType{T} <: DEDataVector{T}
  x::Array{T,1}
  f1::T
end

mutable struct SimType2{T} <: DEDataVector{T}
  x::Vector{T}
  y::Vector{T}
  u::Vector{T}
end

@testset "DEDataVector" begin
  f = function (du,u,p,t)
    du[1] = -0.5*u[1] + u.f1
    du[2] = -0.5*u[2]
  end

  tstop1 = [5.]
  tstop2 = [8.]
  tstop = [5.;8.]

  condition = function (u,t,integrator)
    t in tstop1
  end

  affect! = function (integrator)
    for c in OrdinaryDiffEq.full_cache(integrator)
      c.f1 = 1.5
    end
  end

  save_positions = (true,true)

  cb = DiscreteCallback(condition, affect!; save_positions=save_positions)

  condition2 = function (u,t,integrator)
    t in tstop2
  end

  affect2! = function (integrator)
    for c in OrdinaryDiffEq.full_cache(integrator)
      c.f1 = -1.5
    end
  end

  save_positions = (true,true)

  cb2 = DiscreteCallback(condition2, affect2!, save_positions=save_positions)

  cbs = CallbackSet(cb,cb2)

  u0 = SimType{Float64}([10;10], 0.0)
  prob = ODEProblem(f,u0,(0.0,10.0))
  sol = solve(prob,Tsit5(),callback = cbs, tstops=tstop)

  sol(1.5:0.5:2.5)

  @test [sol[i].f1 for i in eachindex(sol)] == [fill(0., 9);1.5*ones(5);-1.5*ones(4)]

  A = Matrix(Diagonal([0.3,0.6,0.9]))
  B = transpose([1 2 3])
  C = [1/3 1/3 1/3]

  function mysystem(t,x,dx,p,u)
      ucalc = u(x,p,t)
      x.u = ucalc
      x.y = C*x.x
      dx .= A*x.x + B*x.u
  end

  input = (x,p,t)->(1*one(t)≤t≤2*one(t) ? [one(t)] : [zero(t)])
  prob = DiscreteProblem((dx,x,p,t)->mysystem(t,x,dx,p,input), SimType2(fill(0., 3), fill(0., 1), fill(0., 1)), (0//1,4//1))
  sln = solve(prob, FunctionMap(scale_by_time=false), dt = 1//10)

  u1 = [sln[idx].u for idx in 1:length(sln)]
  u2 = [sln(t).u for t in range(0,stop=4,length=41)]
  @test any(x->x[1]>0, u1)
  @test any(x->x[1]>0, u2)

  sln = solve(prob, FunctionMap(scale_by_time=true), dt = 1//10)

  u1 = [sln[idx].u for idx in 1:length(sln)]
  u2 = [sln(t).u for t in range(0,stop=4,length=41)]
  @test any(x->x[1]>0, u1)
  @test any(x->x[1]>0, u2)

  sln = solve(prob, Euler(), dt = 1//10)

  @test u1 == [sln[idx].u for idx in 1:length(sln)] # Show that discrete is the same
  u1 = [sln[idx].u for idx in 1:length(sln)]
  u2 = [sln(t).u for t in range(0,stop=4,length=41)]
  @test any(x->x[1]>0, u1)
  @test any(x->x[1]>0, u2)

  sln = solve(prob, DP5(), dt = 1//10, adaptive=false)

  u1 = [sln[idx].u for idx in 1:length(sln)]
  u2 = [sln(t).u for t in range(0,stop=4,length=41)]
  @test any(x->x[1]>0, u1)
  @test any(x->x[1]>0, u2)
end

######################
# DEDataMatrix
mutable struct SimTypeg{T,T2,N} <: DEDataArray{T,N}
  x::Array{T,N} # two dimensional
  f1::T2
end

@testset "DEDataMatrix" begin
  tstop1 = [10.0]
  tstop2 = [300.]

  function mat_condition(u,t,integrator)
    t in tstop1
  end

  function mat_condition2(u,t,integrator)
    t in tstop2
  end

  function mat_affect!(integrator)
    for c in OrdinaryDiffEq.full_cache(integrator)
      c.f1 = +1.0
    end
  #  integrator.u[1,1] = 0.001
  end

  function mat_affect2!(integrator)
    for c in OrdinaryDiffEq.full_cache(integrator)
      c.f1 = 0.0
    end
  end

  save_positions = (true,true)
  cb = DiscreteCallback(mat_condition, mat_affect!, save_positions=save_positions)
  save_positions = (false,true)
  cb2 = DiscreteCallback(mat_condition2, mat_affect2!, save_positions=save_positions)
  cbs = CallbackSet(cb,cb2)

  function sigmoid(du,u,p,t)
    du[1,1] = 0.01*u[1,1]*(1-u[1,1]/20)
    du[1,2] = 0.01*u[1,2]*(1-u[1,2]/20)
    du[2,1] = 0.01*u[2,1]*(1-u[2,1]/20)
    du[2,2] = u.f1*du[1,1]
  end

  u0 = SimTypeg(fill(0.00001,2,2),0.0)
  tspan = (0.0,3000.0)
  prob = ODEProblem(sigmoid,u0,tspan)

  tstop =[tstop1;tstop2]
  sol = solve(prob,Tsit5(),callback = cbs, tstops=tstop)
  @test_broken solve(prob,Rodas4(),callback = cbs, tstops=tstop) isa ODESolution
  @test_broken sol = solve(prob,Kvaerno3(),callback = cbs, tstops=tstop) isa ODESolution
  sol = solve(prob,Rodas4(autodiff=false),callback = cbs, tstops=tstop)
  sol = solve(prob,Kvaerno3(autodiff=false),callback = cbs, tstops=tstop)
end

mutable struct CtrlSimTypeScalar{T,T2} <: DEDataVector{T}
    x::Vector{T}
    ctrl_x::T2 # controller state
    ctrl_u::T2 # controller output
end

@testset "Interpolation Sides" begin
  function f(x,p,t)
      x.ctrl_u .- x
  end
  function f(dx,x,p,t)
      dx[1] = x.ctrl_u - x[1]
  end

  ctrl_f(e, x) = 0.85(e + 1.2x)

  x0 = CtrlSimTypeScalar([0.0], 0.0, 0.0)
  prob = ODEProblem{false}(f, x0, (0.0, 8.0))

  function ctrl_fun(int)
      e = 1 - int.u[1] # control error

      # pre-calculate values to avoid overhead in iteration over cache
      x_new = int.u.ctrl_x + e
      u_new = ctrl_f(e, x_new)

      # iterate over cache...
      if DiffEqBase.isinplace(int.sol.prob)
        for c in full_cache(int)
            c.ctrl_x = x_new
            c.ctrl_u = u_new
        end
      end
  end

  integrator = init(prob, Tsit5())

  # take 8 steps with a sampling time of 1s
  Ts = 1.0
  for i in 1:8
      ctrl_fun(integrator)
      DiffEqBase.step!(integrator,Ts,true)
  end

  sol = integrator.sol
  @test [u.ctrl_u for u in sol.u[2:end]] == [sol(t).ctrl_u for t in sol.t[2:end]]

  prob = ODEProblem{true}(f, x0, (0.0, 8.0))

  integrator = init(prob, Rosenbrock23())

  # take 8 steps with a sampling time of 1s
  Ts = 1.0
  for i in 1:8
      ctrl_fun(integrator)
      DiffEqBase.step!(integrator,Ts,true)
  end

  sol = integrator.sol
  @test [u.ctrl_u for u in sol.u[2:end]] == [sol(t).ctrl_u for t in sol.t[2:end]]
end


text:
## check_input.jl : common functions to perform input argument checking
#
# Common functions for input argument checking.
#
# This file is part of MultilevelEstimators.jl - A Julia toolbox for
# Multilevel Monte Carlo Methods (c) Pieterjan Robbe, 2019

## checks ##
for (descr, sym) in zip(("equal to", "larger than", "larger than or equal to", "smaller than", "smaller than or equal to"), (:(==),:>, :≥, :<, :≤))
    @eval begin
        $(Symbol("check_", join(split(descr), "_")))(type_name, x, parameter_name, y) = all($(sym).(x, y)) || throw(ArgumentError(string("in ", type_name, ", ", parameter_name, " must be ", $(descr), " ", y, ".")))
    end
end
function check_ndims(x, type_name, parameter_name, m, msg)
    ndims(x) == m || throw(ArgumentError(string("in ", type_name, ", ", msg, ", got ", ndims(x), ", expected ", m)))
end

check_finite(type_name, x, parameter_name) = isfinite(x) || throw(ArgumentError(string("in ", type_name, ", ", parameter_name, " must be finite.")))

check_ordered(type_name, a, b, parameter_name_1, parameter_name_2) = a < b || throw(ArgumentError(string("in ", type_name, ", ", parameter_name_1, " must be smaller than ", parameter_name_2)))

check_ordered_or_equal(type_name, a, b, parameter_name_1, parameter_name_2) = a ≤ b || throw(ArgumentError(string("in ", type_name, ", ", parameter_name_1, " must be smaller than or equal to ", parameter_name_2)))

check_type(type_name, x, parameter_name, requested_type) = x isa requested_type || throw(ArgumentError(string("in ", type_name, ", ", parameter_name, " must be of type ", requested_type, ", got ", typeof(x))))


text:
using Vega
using LDA
using RDatasets

# Try with Iris
iris = data("datasets", "iris")

# Create data matrix
data = matrix(iris[:, 2:5])'

# Create labels vector
labels = Array(Int, size(data, 2))
for i = 1:size(data, 2)
    x = iris[:, 6][i]
    if  x == "setosa"
        labels[i] = 1
    elseif x == "virginica"
        labels[i] = 2
    else
        labels[i] = 3
    end
end

# lineaire
rbf_res = lda(data, labels, 3, linear)
println("Printing results")
print_2Ddecision(rbf_res, data, labels)

# Poly 2
rbf_res = lda(data, labels, 3, polynomial, 2)
println("Printing results")
print_2Ddecision(rbf_res, data, labels)

# Poly 3
rbf_res = lda(data, labels, 3, polynomial, 3)
println("Printing results")
print_2Ddecision(rbf_res, data, labels)

# Poly 4
rbf_res = lda(data, labels, 3, polynomial, 4)
println("Printing results")
print_2Ddecision(rbf_res, data, labels)

# RBF 1.0
rbf_res = lda(data, labels, 3, rbf, 1.0)
println("Printing results")
print_2Ddecision(rbf_res, data, labels)

# RBF 0.4
rbf_res = lda(data, labels, 3, rbf, 0.4)
println("Printing results")
print_2Ddecision(rbf_res, data, labels)

# RBF 0.1
rbf_res = lda(data, labels, 3, rbf, 0.1)
println("Printing results")
print_2Ddecision(rbf_res, data, labels)


text:
"""
    count_nucleotides(strand)

The count of each nucleotide within `strand` as a dictionary.

Invalid strands raise a `DomainError`.

"""
function count_nucleotides(strand)
    res = Dict('A' => 0, 'C' => 0, 'G' => 0, 'T' => 0)
    for c in strand
        if c in ['A','C','G','T']
            res[c] = get(res,c,0)+1
        else
            throw(DomainError(c, "expects A,C,G, or T"))
        end
    end
    return res
end

max_stars_repo_name: stepinski/julia-in

max_issues_repo_name: stepinski/julia-init

max_forks_repo_name: stepinski/julia-init


text:
using JuliaVariables
using MLStyle
import MacroTools

unwrap_scoped(ex) = @match ex begin
    Expr(:scoped, _, a) => unwrap_scoped(a)
    Expr(head, args...) => Expr(head, map(unwrap_scoped, args)...)
    a => a
end

globals(s::Symbol) = [s]
globals(x) = []

function globals(ex::Expr)
    branch(head, newargs) = union(newargs...)

    function leaf(v::JuliaVariables.Var)
        v.is_global ? [v.name] : Symbol[]
    end

    leaf(x) = []

    solved_ex = unwrap_scoped(solve_from_local!(simplify_ex(ex)))

    return foldall(leaf, branch)(solved_ex)
end

function graph(T, m::Model)
    g = T()

    mvars = variables(m)
    for v in mvars
        SimplePosets.add!(g, v)
    end

    for (v, expr) in pairs(m.vals)
        for p in globals(expr) ∩ variables(m)
            SimplePosets.add!(g, p, v)
        end
    end

    for (v, expr) in pairs(m.dists)
        for p in globals(expr) ∩ variables(m)
            SimplePosets.add!(g, p, v)
        end
    end

    g
end

export digraph
digraph(m::Model) = graph(SimpleDigraph{Symbol}, m)

export poset
poset(m::Model) = graph(SimplePoset{Symbol}, m)


text:
module TimeToLive

export TTL

using Dates: Period

struct Node{T}
    v::T
    id::Symbol

    Node{T}(v::T) where T = new(v, gensym())
end

Base.:(==)(a::Node, b::Node) = a.v == b.v

"""
    TTL(ttl::Period; refresh_on_access::Bool=true) -> TTL{Any, Any}
    TTL{K, V}(ttl::Period; refresh_on_access::Bool=true) -> TTL{K, V}

A [TTL](https://en.wikipedia.org/wiki/Time_to_live) cache.
If `refresh_on_access` is set, expiries are reset whenever they are accessed.
"""
struct TTL{K, V} <: AbstractDict{K, V}
    d::Dict{K, Node{V}}
    ttl::Period
    refresh::Bool

    function TTL{K, V}(ttl::Period; refresh_on_access::Bool=true) where {K, V}
        return new(Dict{K, Node{V}}(), ttl, refresh_on_access)
    end
    function TTL(ttl::Period; refresh_on_access::Bool=true)
        return TTL{Any, Any}(ttl; refresh_on_access=refresh_on_access)
    end
end

# TODO: These functions have race conditions.
function delete_later(t::TTL, k, v::Node)
    id = v.id
    sleep(t.ttl)
    haskey(t, k) && t.d[k].id === id && delete!(t, k)
end
Base.get(t::TTL, key, default) = haskey(t.d, key) ? t.d[key].v : default
Base.delete!(t::TTL, key) = (delete!(t.d, key); t)
Base.empty!(t::TTL) = (empty!(t.d); t)
Base.getindex(t::TTL, k) = (t.refresh && touch(t, k); t.d[k].v)
Base.getkey(t::TTL, key, default) = getkey(t.d, k)
Base.length(t::TTL) = length(t.d)
Base.pop!(t::TTL) = (p = pop!(t.d); p.first => p.second.v)
Base.pop!(t::TTL, key) = pop!(t.d, key).v
Base.push!(t::TTL, p::Pair) = (t[p.first] =  p.second; t)
Base.sizehint!(t::TTL, newsz) = (sizehint!(t.d, newsz); t)
Base.sort(t::TTL) = sort(Dict(collect(t)))

function Base.iterate(t::TTL, ks=keys(t.d))
    isempty(ks) && return nothing
    k = first(ks)
    return k => t.d[k].v, Iterators.drop(ks, 1)
end

function Base.setindex!(t::TTL{K, V}, v, k) where {K, V}
    node = Node{V}(v)
    @async delete_later(t, k, node)
    return t.d[k] = node
end

# Extras + restrictions.

"""
    touch(t::TTL, k)

Reset the expiry time for the value at `t[k]`.
"""
function Base.touch(t::TTL{K, V}, k) where {K, V}
    t.d[k] = Node{V}(t.d[k].v)  # Change the ID.
    @async delete_later(t, k, t.d[k])
    return nothing
end

# There's no way to properly copy the deletion tasks.
Base.copy(::TTL) = error("TTL cannot be copied")
Base.deepcopy(::TTL) = error("TTL cannot be copied")

end


text:
@testset "Implicit broadcasting" begin
    # test that tensor broadcasting behaves exactly like array broadcasting
    a1 = reshape(Vector(1:9), 3, 3)
    @test Tensor(a1) .^ 2 == Tensor(a1 .^ 2)

    a2 = rand(3, 3, 3)
    @test (Tensor(a2) .< 0.5) == Tensor(a2 .< 0.5)

    a3 = rand(3, 3, 3)
    @test (Tensor(a2) .< Tensor(a3)) == Tensor(a2 .< a3)

    a4 = randn(4, 2, 5, 3)
    @test abs.(Tensor(a4)) == Tensor(abs.(a4))

    # more complex broadcasting example taken from a blog post at
    # https://julialang.org/blog/2018/05/extensible-broadcast-fusion/
    a5 = [1, 2, 3]
    a6 = [10 20 30 40]
    a7 = 10

    # test broadcasting behaviour of the expression (a5 .+ a6) ./ a7
    @test (Tensor(a5) .+ a6) ./ a7 == Tensor((a5 .+ a6) ./ a7)
    @test (a5 .+ Tensor(a6)) ./ a7 == Tensor((a5 .+ a6) ./ a7)
    @test (a5 .+ a6) ./ Tensor(a7) == Tensor((a5 .+ a6) ./ a7)
    @test (Tensor(a5 .+ a6)) ./ a7 == Tensor((a5 .+ a6) ./ a7)

    @test (Tensor(a5 .+ a6)) ./ Tensor(a7) == Tensor((a5 .+ a6) ./ a7)
    @test (Tensor(a5) .+ a6) ./ Tensor(a7) == Tensor((a5 .+ a6) ./ a7)
    @test (a5 .+ Tensor(a6)) ./ Tensor(a7) == Tensor((a5 .+ a6) ./ a7)
    @test (Tensor(a5) .+ Tensor(a6)) ./ a7 == Tensor((a5 .+ a6) ./ a7)

    @test (Tensor(a5) .+ Tensor(a6)) ./ Tensor(a7) == Tensor((a5 .+ a6) ./ a7)
end

@testset "Explicit broadcasting" begin
    a1 = [1 2 3 4 5]
    a2 = [
        1 3 5 7 9
        2 4 6 8 10
    ]
    t1 = Tensor(a1)
    t2 = Tensor(a2)

    # explicit broadcasting works as expected
    @test broadcast(+, t1, a2) == Tensor([
        2 5 8 11 14
        3 6 9 12 15
    ])
    @test broadcast(+, a1, t2) == Tensor([
        2 5 8 11 14
        3 6 9 12 15
    ])
    @test broadcast(+, t1, t2) == Tensor([
        2 5 8 11 14
        3 6 9 12 15
    ])

    # explicit in-place broadcasting is also possible
    @test t2 != Tensor([
        2 5 8 11 14
        3 6 9 12 15
    ])
    broadcast!(+, t2, t1, t2)
    @test t2 == Tensor([
        2 5 8 11 14
        3 6 9 12 15
    ])
end


text:
using JuliaFormatter

include(joinpath(@__DIR__, "..", "formatter_options.jl"))

format(@__FILE__; formatter_options...)


text:
#==============================================================================#
# SQS/test/runtests.jl
#
# Copyright OC Technology Pty Ltd 2014 - All rights reserved
#==============================================================================#

using AWSSQS
using AWSCore
using Test
using Dates

AWSCore.set_debug_level(1)
aws = AWSCore.aws_config(region="ap-southeast-2")

@testset "AWSSQS.jl" begin
    include("queue.jl")
end


text:
mutable struct History
    default::Vector{Organism}
    detailed::Union{Vector{DetailedOrganism}, Nothing}

    ordered_by::Dict{PropertyType, Dict{Any, Vector{DetailedOrganism}}}

    History(organisms) = new(organisms, nothing, Dict())
end

get_default_history(history::History) = history.default

function get_detailed_history(history::History)
    if history.detailed === nothing
        history.detailed = [DetailedOrganism(org, 0, 0, Dict()) for org in history.default]

        id_to_organism = Dict([get_id(o) => o for o in history.detailed])

        # calculate the numer of offsprings
        for organism in history.detailed
            if organism.parent != -1
                parent = id_to_organism[organism.parent]
                parent.num_offsprings += 1

                if parent.genotype == organism.genotype
                    parent.num_identical_offsprings += 1
                end
            end
        end

        # calculate populations
        currently_alive = []

        current_time = 0
        born_at_current_time = []

        for organism in history.detailed
            if length(born_at_current_time) == 0
                current_time = organism.time_birth
                born_at_current_time = [organism]

            elseif organism.time_birth == current_time
                push!(born_at_current_time, organism)
                continue

            else
                # remove dead from currently_alive and add newborn
                currently_alive = [o for o in currently_alive if o.time_death < current_time]
                append!(currently_alive, born_at_current_time)

                for born in born_at_current_time
                    population = born.population

                    for alive in currently_alive
                        genotype = alive.genotype

                        if haskey(population, genotype)
                            population[genotype] += 1
                        else
                            population[genotype] = 1
                        end
                    end
                end

                current_time = organism.time_birth
                born_at_current_time = []
            end
        end

    end

    return history.detailed
end

function get_grouped_by(history::History, properties::PropertyType...)
    if !haskey(history.ordered_by, properties)
        history.ordered_by[properties] = _group_by(get_detailed_history(history), get_group_function(properties))
    end

    history.ordered_by[properties]
end

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
# set to true to support intel fortran compiler

try
	println("************************Trying to build ParSpMatVec ******************************")
	useIntelFortran = false

	# construct absolute path
	depsdir  = splitdir(Base.source_path())[1]
	builddir = joinpath(depsdir,"builds")
	srcdir   = joinpath(depsdir,"src")

	println("=== Building ParSpMatVec ===")
	println("depsdir  = $depsdir")
	println("builddir = $builddir")
	println("srcdir   = $srcdir")
	println("useIntel = $useIntelFortran")

	if !isdir(builddir)
		println("creating build directory")
		mkdir(builddir)
		if !isdir(builddir)
			error("Could not create build directory")
		end
	end

	@static if Sys.isunix()

		src1 = joinpath(srcdir,"A_mul_B.f90")
		src2 = joinpath(srcdir,"Ac_mul_B.f90")
		outfile = joinpath(builddir,"ParSpMatVec.so")

		if useIntelFortran
			run(`ifort -O3 -xHost -fPIC -fpp -openmp -integer-size 64 -diag-disable=7841 -shared  $src1 $src2 -o $outfile`)
		else
			println("fortran version")
			# run(`gfortran --version`)
			run(`gfortran -v -O3 -fPIC -cpp -fopenmp -fdefault-integer-8 -shared $src1 $src2 -o $outfile`)
			println("Done compiling.")
		end
	end

	@static if Sys.iswindows()
		src1 = joinpath(srcdir,"A_mul_B.f90")
		src2 = joinpath(srcdir,"Ac_mul_B.f90")
		outfile = joinpath(builddir,"ParSpMatVec.dll")
		run(`gfortran --version`)
		run(`gfortran -v -O3 -cpp -fopenmp -fdefault-integer-8 -shared -DBUILD_DLL  $src1 $src2 -o $outfile`)
	end
catch
	@warn "Warning: Unable to build ParSpMatVec"
end


text:
import Pkg
Pkg.activate("emergency_services")

# load all necessary packages and functions
include("functions/start_load_packages.jl")

# state the parts of the framework that should be executed
# stage 1: only the optimisation will be executed
# stage 2: only the simulation will be executed (requires at least one run of stage 1)
# both:    both stages will be executed
    framework   = "both"::String

# state the name of the problem instance that should be solved, possible examples
# are "0510","1008","1508","2040" where the name corresponds to the number of BAs
    problem     = "2040"::String

# state the name of the subproblem that should be solved, for example "one_location",
# "basic", "two_locations", eg.
    subproblem  = "nine_locations"::String

# state the strength of the compactness and contiguity constraints
# C0 = no contiguity constraints (no compactness)
# C1 = solely contiguity constrains (no compactness)
# C2 = contiguity and normal ompactness constraints
# C3 = contiguity and strong compactness constraints
# For more details take a look at the article this program is based on
    compactness = "C2"::String

# state the number of simulations that should be executed due to the variability of the driving time
    sim_number = 100::Int64

# state the main input parameters for the optimisation (framework stage 1)
    number_districts = 9::Int64       # number of districts that should be opened
    max_drive        = 30.0::Float64  # maximum driving distance (minutes) to district border
    nearby_districts = 2::Int64       # minimal number of districts within nearby radius
    nearby_radius    = 30.0::Float64  # maximal driving time to nearby district center
    fixed_locations  = 0::Int64       # number of current locations that should not be moved
    plot_district    = true::Bool     # state whether the resulting district should be plotted

# state the optimisation options
    solver = "gurobi"
    const GRB_ENV = Gurobi.Env()
    const optcr   = 0.000::Float64         # allowed gap
    const reslim  = 21600::Int64           # maximal duration of optimisation in seconds
    const cores   = 8::Int64               # number of CPU cores
    const nodlim  = 1000000::Int64         # maximal number of nodes
    const iterlim = 100000000::Int64       # maximal number of iterations
    const silent  = false::Bool            # state whether to surpress the optimisation log

# state the weight of each priority for the driving time
# important: the weight has to equal the number of priorities
# in the incident data stage
    prio_weight = [1, 1, 1, 1, 1]

# state the main parameters for the simulation (framework stage 2)
    const min_capacity   = 2::Int64        # minimal capacity for each district during each weekhour
    const exchange_prio  = 5::Int64        # till which priority can cars be exchanged to foreign districts
    const patrol_prio    = 5::Int64        # till which priority can patrol cars be dispatched as backup
    const patrol_ratio   = 0.00::Float64   # which proportion of the ressources is assigned to be on patrol
    const patrol_time    = 5.00::Float64   # maximal driving time for a patrol car dispatch as help
    const patrol_area    = 5::Int64        # time spent patrolling per basic area
    const backlog_max    = 60::Int64       # maximal average backlog (minutes) per car in district for exchange
    const max_queue      = 80::Int64       # maximal length of the queue of incidents per district and priority
    const real_capacity  = false::Bool     # state whether a predefined capacity plan should be loaded
    const drop_incident  = 360::Int64      # total number of minutes after which an incident
                                           # will leave the queue even if it's not fully fulfilled

# state the main parameters for the capacity estimation if no capacity plan is given
    const total_capacity   = 50.0::Float64   # average capacity per hour in the area over the incident timeframe
    const capacity_service = 0.90::Float64    # alpha service level for weekhour workload estimation

# state how many cars should be reserved for the own district per incident priority
# during exchanges to other districts this threshold will not be crossed
# important: the exchange_reserve has to equal the number of priorities
    exchange_reserve = [0,0,0,0,0]

# state whether to save the plotted results
    const save_plots = true::Bool # save the plots

# lock to avoid racing conditions during the simulation
    ren_lock = ReentrantLock()

# load the input data
    include("2_define_data.jl")
    print("\n Input data sucessfully loaded.")

# prepare the input data for stage 1 (also neccessary for stage 2!)
# append the weekhour to the incident dataset
    incidents.weekhour = epoch_weekhour.(incidents.epoch)

# prepare the workload for each BA if allocated to a potential center
    if isfile("data/$problem/workload_$problem.csv")
        workload = readdlm("data/$problem/workload_$problem.csv", Float64)
        we = readdlm("data/$problem/weight_$problem.csv", Float64)[:,1]
        print("\n Workload and weight recovered from saved file.")
    else
        dur = @elapsed workload, we = workload_calculation(incidents::DataFrame,
                                                        prio_weight::Vector{Int64},
                                                        drivingtime::Array{Float64,2},
                                                        traffic::Array{Float64,2},
                                                        size(airdist,1)::Int64)
        workload = round.(workload, digits = 3)
        writedlm("data/$problem/workload_$problem.csv", workload)
        writedlm("data/$problem/weight_$problem.csv", we)
        print("\n Finished workload and weight calculation after ",dur," seconds.")
    end

# prepare the sets for the contiguity and compactness constraints
    dur = @elapsed N, M, card_n, card_m = sets_m_n(airdist::Array{Float64,2}, size(airdist,1)::Int64)
    print("\n Finished set calculation after ", dur," seconds.")
    print("\n Input sucessfully prepared for both stages.")

# state which problem and constraints will be the aim of the framework
    print("\n Framework will be applied to the problem ",subproblem," with the constraint set ", compactness,".")

# district optimisation (framework stage 1)
if framework != "stage 2"
    print("\n Starting optimisation.")

# Create a Dataframe to save the main results
    opt_out = DataFrame(prbl = String[], cmpct = String[], objv = Float64[], gp = Float64[], time = Float64[])

# Start the optimisation model
    X, gap, objval, scnds = model_and_optimise(solver, number_districts, drivingtime, we,
                                                workload, adjacent, N, M, card_n, card_m, max_drive,
                                                compactness, potential_locations, nearby_radius,
                                                nearby_districts, current_locations, fixed_locations)
    districts =  create_frame(X, size(airdist,1))
    print("\n Duration of the optimisation: ", scnds, " seconds")

    # Save the optimisation results and district layouts
    CSV.write("results_stage1/$problem/district_layout_$(problem)_$(subproblem)_$(compactness).csv", districts)
    push!(opt_out, (prbl = problem, cmpct = compactness, objv = objval, gp = gap, time = scnds))
    CSV.write("results_stage1/$problem/optimisation_$(problem)_$(subproblem)_$(compactness).csv", opt_out)
    print("\n Results from stage 1 written to file.")

# Plot the resulting district layout
    if plot_district && hexshape !== nothing
        district_plot = plot_generation(districts, hexshape)
        display(district_plot)
        if save_plots == 1
            savefig("graphs/$problem/district_layout_$(problem)_$(subproblem)_$(compactness).pdf")
        end
    end

# End stage 1 of the framework
end

# emergency service simulation (framework stage 2)
if framework != "stage 1"

# prepare the input data for the simulation
# load the district layout in case solely stage 2 is executed
    if isfile("results_stage1/$problem/district_layout_$(problem)_$(subproblem)_$(compactness).csv")
        districts = CSV.read("results_stage1/$problem/district_layout_$(problem)_$(subproblem)_$(compactness).csv", DataFrame)
        districts[!,:index] = convert.(Int64,districts[!,:index])
        districts[!,:location] = convert.(Int64,districts[!,:location])
        print("\n Districts loaded from file.")
    else
        error("Perform stage 1 first, to determine the district layouts.")
    end

# prepare the input data for stage 2
    sim_data = prepare_simulation_data(incidents, districts)

## simulation_capacity: apply a heuristic to derive the capacity of each location in each weekhour
    simulation_capacity, capacity_plot = capacity_heuristic(incidents,sim_data,total_capacity,min_capacity)

## ressource_flow: fill the intial ressource flow matrix of the simulation
    ressource_flow = ressource_flow_matrix(sim_data, incidents, simulation_capacity)

    print("\n Input sucessfully prepared for simulation.")
    print("\n Ressource matrix build.","\n")

# pepare the execution of multiple simulation runs
    weekly_location = DataFrame()
    weekly_priority = DataFrame()
    capacity_status = DataFrame()
    main_results    = DataFrame()

# start the simulation
    print("\n Starting ", sim_number," emergency service simulations.")
    Threads.@threads for sim_run = 1:sim_number
        rfw = copy(ressource_flow)
        smd = copy(sim_data)
        scnds = @elapsed smd,rfw = part2_simulation!(districts::DataFrame,
                                                     incidents::DataFrame,
                                                     smd::DataFrame,
                                                     rfw::Array{Int64,3},
                                                     drivingtime::Array{Float64,2},
                                                     traffic::Array{Float64,2},
                                                     max_drive::Float64,
                                                     drop_incident::Int64,
                                                     exchange_reserve::Vector{Int64})
        main_results_local, weekly_location_local, weekly_priority_local,
        capacity_status_local = evaluate_results(incidents, smd, rfw, simulation_capacity)

    # save the results for later evaluation
        weekly_location_local[:,:simulation_run] .= sim_run
        weekly_priority_local[:,:simulation_run] .= sim_run
        capacity_status_local[:,:simulation_run] .= sim_run
        main_results_local[:,:simulation_run]    .= sim_run

        Threads.lock(ren_lock) do
            append!(weekly_location, weekly_location_local)
            append!(weekly_priority, weekly_priority_local)
            append!(capacity_status, capacity_status_local)
            append!(main_results,    main_results_local)
        end
        print("\n Simulation ",nrow(main_results)," of ", sim_number," completed after ", scnds, " seconds on thread ",Threads.threadid())
    end

    # Aggregate the results of all simulations
    weekly_location = groupby(weekly_location,[:location_responsible, :weekhour])
    weekly_location = combine(weekly_location, [n => mean => n for n in names(weekly_location) if n != "simulation_run"])
    weekly_location = round.(weekly_location, digits = 5)
    weekly_priority = groupby(weekly_priority,[:weekhour,:priority])
    weekly_priority = combine(weekly_priority, [n => mean => n for n in names(weekly_priority) if n != "simulation_run"])
    weekly_priority = round.(weekly_priority, digits = 5)
    capacity_status = groupby(capacity_status,[:weekhour,:variable])
    capacity_status = combine(capacity_status, :value => mean => :value)
    main_results    = combine(main_results, [n => mean => n for n in names(main_results) if n != "simulation_run"])
    main_results    = round.(main_results, digits = 5)
    main_priority   = groupby(weekly_priority, :priority)
    main_priority   = combine(main_priority, [n => mean => n for n in names(main_priority) if n != "weekhour"])
    main_priority   = round.(main_priority, digits = 5)

    print("\n\n Main results of simulation:")
    print("\n Average dispatch time first car:    ", main_results[1,:dispatch_time_first])
    for prio = 1:nrow(main_priority)
        print("\n   Average dispatch time priority ", prio,":   ", main_priority[prio,:dispatch_time_first])
    end
    print("\n Average dispatch time all cars:     ", main_results[1,:dispatch_time_all])
    print("\n Average response time first car:    ", main_results[1,:response_time_first])
    for prio = 1:nrow(main_priority)
        print("\n   Average response time priority ", prio,":   ", main_priority[prio,:response_time_first])
    end
    print("\n Average response time all cars:     ", main_results[1,:response_time_all])
    print("\n Average exchange ratio:             ", main_results[1,:exchange_ratio])
    print("\n Average ratio undispatched cars:    ", main_results[1,:ratio_cars_undispatched])
    print("\n Unfullfilled calls for service:     ", main_results[1,:incidents_unfulfilled])

    # Plot the aggregated results
    plot_simulation_results(weekly_location, weekly_priority, capacity_status, capacity_plot, subproblem, compactness)

    # Save the main results of the simulation
    CSV.write("results_stage2/$problem/simulation_$(problem)_$(subproblem)_$(compactness).csv", main_results)
    print("\n Results from stage 2 written to file. Framework ended.")
    print("\n")
end
print("\n")


text:
""""
    `dare(A, B, Q, R)`

Compute `X`, the solution to the discrete-time algebraic Riccati equation,
defined as A'XA - X - (A'XB)(B'XB + R)^-1(B'XA) + Q = 0, where Q>=0 and R>0

Algorithm taken from:
Laub, "A Schur Method for Solving Algebraic Riccati Equations."
[http://dspace.mit.edu/bitstream/handle/1721.1/1301/R-0859-05666488.pdf](http://dspace.mit.edu/bitstream/handle/1721.1/1301/R-0859-05666488.pdf)

Function taken from the official `ControlSystems.jl` repository
[https://github.com/JuliaControl/ControlSystems.jl](https://github.com/JuliaControl/ControlSystems.jl)
"""
function dare(A, B, Q, R)
    if (!ishermitian(Q) || minimum(eigvals(real(Q))) < 0)
        error("Q must be positive-semidefinite.");
    end
    if (!isposdef(R))
        error("R must be positive definite.");
    end

    n = size(A, 1);

    E = [
        Matrix{Float64}(I, n, n) B/R*B';
        zeros(size(A)) A'
    ];
    F = [
        A zeros(size(A));
        -Q Matrix{Float64}(I, n, n)
    ];

    QZ = schur(F, E);
    QZ = ordschur(QZ, abs.(QZ.alpha./QZ.beta) .< 1);

    return QZ.Z[(n+1):end, 1:n]/QZ.Z[1:n, 1:n];
end

"""
    `dlqr(A, B, Q, R)`, `dlqr(sys, Q, R)`

Calculate the optimal gain matrix `K` for the state-feedback law `u[k] = K*x[k]` that
minimizes the cost function:

J = sum(x'Qx + u'Ru, 0, inf).

For the discrte time model `x[k+1] = Ax[k] + Bu[k]`. See also `lqg`.

Usage example:

```julia
using LinearAlgebra # For identity matrix I
h = 0.1
A = [1 h; 0 1]
B = [0;1]
C = [1 0]
sys = ss(A,B,C,0, h)
Q = I
R = I
L = dlqr(A,B,Q,R) # lqr(sys,Q,R) can also be used
u(x,t) = -L*x # Form control law,
t=0:h:5
x0 = [1,0]
y, t, x, uout = lsim(sys,u,t,x0=x0)
plot(t,x, lab=["Position"  "Velocity"], xlabel="Time [s]")
```

Function taken from the official `ControlSystems.jl` repository
[https://github.com/JuliaControl/ControlSystems.jl](https://github.com/JuliaControl/ControlSystems.jl)
"""
function dlqr(A, B, Q, R)
    S = dare(A, B, Q, R)
    K = (B'*S*B + R)\(B'S*A)
    return K
end

"""
    `dkalman(A, C, R1, R2)` kalman(sys, R1, R2)`

Calculate the optimal Kalman gain for discrete time systems.

Function taken from the official `ControlSystems.jl` repository
[https://github.com/JuliaControl/ControlSystems.jl](https://github.com/JuliaControl/ControlSystems.jl)
"""
dkalman(A, C, R1,R2) = Matrix(dlqr(A',C',R1,R2)')

"""
`place(A, B, p)`, `place(sys::StateSpace, p)`

Calculate gain matrix `K` such that the poles of `(A-BK)` in are in `p`. Uses
Ackermann's formula.

Function taken from the official `ControlSystems.jl` repository
[https://github.com/JuliaControl/ControlSystems.jl](https://github.com/JuliaControl/ControlSystems.jl)
"""
function place(A, B, p)
    n = length(p)
    n != size(A,1) && error("Must define as many poles as states")
    n != size(B,1) && error("A and B must have same number of rows")
    if size(B,2) == 1
        acker(A,B,p)
    else
        error("place only implemented for SISO systems")
    end
end


text:
struct BlockBasisCoo <: Field end

function evaluate_field!(cache,b::BlockBasisCoo,x)
  msg =
  """
  Not implemented! We not support iteration of vectors of block bases before
  evaluation. This is unlikely to be be needed in the future.
  """
  @notimplemented msg
end

struct VectorOfBlockBasisCoo <: AbstractVector{BlockBasisCoo}
  blocks::Tuple
  blockids
  axes
end

Base.size(v::VectorOfBlockBasisCoo) = (length(first(v.blocks)) ,)
Base.getindex(v::VectorOfBlockBasisCoo,i::Integer) = BlockBasisCoo()

function evaluate_field_array(v::VectorOfBlockBasisCoo,x::AbstractArray)
  blocks = evaluate_field_arrays(v.blocks,x)
  axs = apply(_new_axes,x,v.axes)
  blockids = map(i->(1,i...),v.blockids)
  VectorOfBlockArrayCoo(blocks,blockids,axs)
end

function field_array_gradient(v::VectorOfBlockBasisCoo)
  blocks = map(field_array_gradient,v.blocks)
  VectorOfBlockBasisCoo(blocks,v.blockids,v.axes)
end

function reindex(v::VectorOfBlockBasisCoo,j_to_i::AbstractArray)
  blocks = map(b->reindex(b,j_to_i),v.blocks)
  axs = reindex(v.axes,j_to_i)
  VectorOfBlockBasisCoo(blocks,v.blockids,axs)
end

function compose_field_arrays(v::VectorOfBlockBasisCoo,f)
  blocks = map(b->compose_field_arrays(b,f),v.blocks)
  VectorOfBlockBasisCoo(blocks,v.blockids,v.axes)
end

function insert_array_of_bases_in_block(i::Integer,a,ax1,ax2...)
  blocks = (a,)
  blockids = _compute_blockids(eltype(ax1),i)
  axs = create_array_of_blocked_axes(ax1,ax2...)
  VectorOfBlockBasisCoo(blocks,blockids,axs)
end

#TODO perhaps this should go into Arrays/BlockArraysCoo.jl
function create_array_of_blocked_axes(axs...)
  axs = apply(_cat_axes,axs...)
end

function _compute_blockids(::Type{<:NTuple{1}},i)
  [(i,)]
end

function _compute_blockids(::Type{<:NTuple{2}},i)
  [(1,i)]
end

function _cat_axes(a::NTuple{1},b::NTuple{1})
  (_blockedrange([a[1],b[1]]),)
end

function _cat_axes(a::NTuple{1}...)
  (_blockedrange([map(i->i[1],a)...]),)
end

function _cat_axes(a::NTuple{2},b::NTuple{2})
  @assert length(a[1]) == 1
  @assert length(b[1]) == 1
  ran = (_blockedrange([a[2],b[2]]),)
  _add_singleton_block(ran)
end

function _cat_axes(a::NTuple{2}...)
  ran = (_blockedrange([map(i->i[2],a)...]),)
  _add_singleton_block(ran)
end

function _new_axes(x,ran::NTuple{N,<:BlockedUnitRange} where N)
  np = length(x)
  (blockedrange([np]),ran...)
end

function _new_axes(x,ran::NTuple{N,<:TwoLevelBlockedUnitRange} where N)
  np = length(x)
  r = blockedrange([np])
  (blockedrange([r]),ran...)
end

function _add_singleton_block(ran::NTuple{N,<:Base.OneTo} where N)
  (Base.OneTo(1),ran...)
end

function _add_singleton_block(ran::NTuple{N,<:BlockedUnitRange} where N)
  (blockedrange([1]),ran...)
end

function _add_singleton_block(ran::NTuple{N,<:TwoLevelBlockedUnitRange} where N)
  r = blockedrange([1])
  (blockedrange([r]),ran...)
end

function _blockedrange(a::Vector{<:Base.OneTo})
  blockedrange(map(length,a))
end

function _blockedrange(a::Vector{<:BlockedUnitRange})
  blockedrange(a)
end


text:
module ExaPF

# Standard library
using Printf
using LinearAlgebra
using SparseArrays

import CUDA
import CUDA.CUBLAS
import CUDA.CUSPARSE
import CUDA.CUSOLVER

import ForwardDiff
using KernelAbstractions
const KA = KernelAbstractions
import MathOptInterface
const MOI = MathOptInterface
using TimerOutputs: @timeit, TimerOutput

import Base: show, get

const VERBOSE_LEVEL_HIGH = 3
const VERBOSE_LEVEL_MEDIUM = 2
const VERBOSE_LEVEL_LOW = 1
const VERBOSE_LEVEL_NONE = 0
const TIMER = TimerOutput()

include("utils.jl")
include("architectures.jl")
# Templates
include("models.jl")

# Import submodules
include("autodiff.jl")
using .AutoDiff
include("LinearSolvers/LinearSolvers.jl")
using .LinearSolvers
include("PowerSystem/PowerSystem.jl")
using .PowerSystem

const PS = PowerSystem

# Polar formulation
include("Polar/polar.jl")

# Evaluators
include("Evaluators/Evaluators.jl")

end


text:
"""
Abstract basis

  T: data type
  N: number of nodes
"""
abstract type Basis{T<:AbstractFloat, N} end

nbasis(b::Basis{T,N}) where {T,N} = error("Not implemented!")
nnodes(b::Basis{T,N}) where {T,N} = error("Not implemented!")
degree(b::Basis{T,N}) where {T,N} = error("Not implemented!")
nodes(b::Basis{T,N})  where {T,N} = error("Not implemented!")

eval_basis(b::Basis, x, y)  = error("Not implemented!")
deriv_basis(b::Basis, x, y) = error("Not implemented!")
int_basis(b::Basis, x, y)   = error("Not implemented!")

"""
Abstract polynomial basis.
"""
abstract type PolynomialBasis{T,N} <: Basis{T,N} end

nbasis(b::PolynomialBasis{T,N}) where {T,N} = N
nnodes(b::PolynomialBasis{T,N}) where {T,N} = N
degree(b::PolynomialBasis{T,N}) where {T,N} = N-1


text:
struct Region
    slug::String
    name::String
    sizes::Array{String, 1}
    available::Union{Nothing, Bool}
    features::Array{String, 1}

    function Region(data::Dict{String})
        new(
            data["slug"],
            data["name"],
            data["sizes"],
            data["available"],
            data["features"]
        )
    end
end

function show(io::IO, r::Region)
    print(io, "Region ($(r.name))")
end

function getallregions!(client::AbstractClient)
    response = getdata!(client, joinpath(ENDPOINT, "regions"))

    if floor(response.status/100) == 2 # OK
        body = parse(String(response.body))
        meta = body["meta"]
        links = body["links"]
        data = body["regions"]

        regions = Array{Region, 1}(UndefInitializer(), meta["total"])

        for (i, region) in enumerate(data)
            regions[i] = Region(region)
        end
    else
        error("Received error $(response.status)")
    end

    regions
end


text:
const MATRIX_TYPES = [:nonmatrix, :adjm, :laplacian, :normalized, :scaled]
const DIRECTEDS = [:auto, :directed, :undirected]

abstract type AbstractFeaturedGraph end

"""
    NullGraph()

Null object for `FeaturedGraph`.
"""
struct NullGraph <: AbstractFeaturedGraph end

"""
    FeaturedGraph(graph, node_feature, edge_feature, global_feature, mt, directed)

A feature-equipped graph structure for passing graph to layer in order to provide graph dynamically.
References to graph or features are hold in this type.

# Arguments
- `graph`: should be a adjacency matrix, `SimpleGraph`, `SimpleDiGraph` (from LightGraphs) or `SimpleWeightedGraph`,
`SimpleWeightedDiGraph` (from SimpleWeightedGraphs).
- `node_feature`: node features attached to graph.
- `edge_feature`: edge features attached to graph.
- `gloabl_feature`: gloabl graph features attached to graph.
- `mask`: mask for `graph`.
- `mt`: matrix type for `graph` in matrix form. if `graph` is in matrix form, `mt` is recorded as one of `:adjm`,
`:laplacian`, `:normalized` or `:scaled`. Otherwise, `:nonmatrix` is recorded.
- `directed`: the direction of `graph`. it is `true` for directed graph; it is `false` for undirected graph.
"""
mutable struct FeaturedGraph{T,S<:AbstractMatrix,R<:AbstractMatrix,Q<:AbstractVector} <: AbstractFeaturedGraph
    graph::T
    nf::S
    ef::R
    gf::Q
    mask
    matrix_type::Symbol
    directed::Bool

    function FeaturedGraph(graph::T, nf::S, ef::R, gf::Q,
            mask, mt::Symbol, directed::Bool) where {T,S<:AbstractMatrix,R<:AbstractMatrix,Q<:AbstractVector}
        _check_precondition(graph, nf, ef, mt)
        new{T,S,R,Q}(graph, nf, ef, gf, mask, mt, directed)
    end
    function FeaturedGraph{T,S,R,Q}(graph, nf, ef, gf,
            mask, mt, directed) where {T,S<:AbstractMatrix,R<:AbstractMatrix,Q<:AbstractVector}
        _check_precondition(graph, nf, ef, mt)
        new{T,S,R,Q}(T(graph), S(nf), R(ef), Q(gf), mask, mt, directed)
    end
end

FeaturedGraph() = NullGraph()

## Graph from JuliaGraphs

function FeaturedGraph(graph; directed::Symbol=:auto, T=eltype(graph), N=nv(graph), E=ne(graph),
                       nf=Fill(zero(T), (0, N)), ef=Fill(zero(T), (0, E)), gf=Fill(zero(T), 0))
    @assert directed ∈ DIRECTEDS "directed must be one of :auto, :directed and :undirected"
    dir = (directed == :auto) ? is_directed(graph) : directed == :directed
    mask = Fill(zero(T), (N, N))
    FeaturedGraph(graph, nf, ef, gf, mask, :nonmatrix, dir)
end

## Graph in adjacency list

function FeaturedGraph(graph::AbstractVector{T}; directed::Symbol=:auto, ET=eltype(graph[1]), N=nv(graph), E=ne(graph),
                       nf=Fill(zero(ET), (0, N)), ef=Fill(zero(ET), (0, E)), gf=Fill(zero(ET), 0)) where {T<:AbstractVector}
    @assert directed ∈ DIRECTEDS "directed must be one of :auto, :directed and :undirected"
    dir = (directed == :auto) ? is_directed(graph) : directed == :directed
    mask = Fill(zero(ET), (N, N))
    FeaturedGraph(graph, nf, ef, gf, mask, :nonmatrix, dir)
end

## Graph in adjacency matrix

function FeaturedGraph(graph::AbstractMatrix{T}; directed::Symbol=:auto, N=nv(graph), E=ne(graph),
                       nf=Fill(zero(T), (0, N)), ef=Fill(zero(T), (0, E)), gf=Fill(zero(T), 0)) where {T<:Real}
    @assert directed ∈ DIRECTEDS "directed must be one of :auto, :directed and :undirected"
    dir = (directed == :auto) ? !issymmetric(graph) : directed == :directed
    graph = promote_graph(graph, nf)
    mask = Fill(zero(T), (N, N))
    FeaturedGraph(graph, nf, ef, gf, mask, :adjm, dir)
end

function check_num_node(graph_nv::Real, N::Real)
    if graph_nv != N
        throw(DimensionMismatch("number of nodes must match between graph ($graph_nv) and node features ($N)"))
    end
end

function check_num_edge(graph_ne::Real, E::Real)
    # allow for the number of edge in directed and undirected graph
    if graph_ne != E && 2*graph_ne != E
        throw(DimensionMismatch("number of edges must match between graph ($graph_ne) and edge features ($E)"))
    end
end

check_num_node(graph_nv::Real, nf) = check_num_node(graph_nv, size(nf, 2))
check_num_edge(graph_ne::Real, ef) = check_num_edge(graph_ne, size(ef, 2))

check_num_node(g, nf) = check_num_node(nv(g), nf)
check_num_edge(g, ef) = check_num_edge(ne(g), ef)

function _check_precondition(graph, nf, ef, mt)
    @assert mt ∈ MATRIX_TYPES "matrix_type must be one of :nonmatrix, :adjm, :laplacian, :normalized or :scaled"
    check_num_edge(ne(graph), ef)
    check_num_node(nv(graph), nf)
    return
end

function Base.setproperty!(fg::FeaturedGraph, prop::Symbol, x)
    if prop == :graph
        check_num_node(x, fg.nf)
        check_num_edge(x, fg.ef)
    elseif prop == :nf
        check_num_node(fg.graph, x)
    elseif prop == :ef
        check_num_edge(fg.graph, x)
    end
    setfield!(fg, prop, x)
end

"""
    graph(::AbstractFeaturedGraph)

Get referenced graph.
"""
graph(::NullGraph) = nothing
graph(fg::FeaturedGraph) = fg.graph

"""
    node_feature(::AbstractFeaturedGraph)

Get node feature attached to graph.
"""
node_feature(::NullGraph) = nothing
node_feature(fg::FeaturedGraph) = fg.nf

"""
    edge_feature(::AbstractFeaturedGraph)

Get edge feature attached to graph.
"""
edge_feature(::NullGraph) = nothing
edge_feature(fg::FeaturedGraph) = fg.ef

"""
    global_feature(::AbstractFeaturedGraph)

Get global feature attached to graph.
"""
global_feature(::NullGraph) = nothing
global_feature(fg::FeaturedGraph) = fg.gf

mask(::NullGraph) = nothing
mask(fg::FeaturedGraph) = fg.mask

"""
    has_graph(::AbstractFeaturedGraph)

Check if graph is available or not.
"""
has_graph(::NullGraph) = false
has_graph(fg::FeaturedGraph) = fg.graph != Fill(0., (0,0))
Zygote.@nograd has_graph

"""
    has_node_feature(::AbstractFeaturedGraph)

Check if node feature is available or not.
"""
has_node_feature(::NullGraph) = false
has_node_feature(fg::FeaturedGraph) = !isempty(fg.nf)
Zygote.@nograd has_node_feature

"""
    has_edge_feature(::AbstractFeaturedGraph)

Check if edge feature is available or not.
"""
has_edge_feature(::NullGraph) = false
has_edge_feature(fg::FeaturedGraph) = !isempty(fg.ef)
Zygote.@nograd has_edge_feature

"""
    has_global_feature(::AbstractFeaturedGraph)

Check if global feature is available or not.
"""
has_global_feature(::NullGraph) = false
has_global_feature(fg::FeaturedGraph) = !isempty(fg.gf)
Zygote.@nograd has_global_feature


text:
"""
Module UniformTriMesh

Includes mesh related utility functions
"""

module UniformTriMesh

using SparseArrays
using CommonUtils # for meshgrid

export uniform_tri_mesh
export tri_face_vertices

"""
uniform_tri_mesh(Kx::Int,Ky::Int)

Matlab uniform triangular mesh.

# Examples
```jldoctest
```
"""

function uniform_tri_mesh(Kx,Ky)
        (VY, VX) = meshgrid(LinRange(-1,1,Ky+1),LinRange(-1,1,Kx+1))
        sk = 1
        EToV = zeros(Int,2*Kx*Ky,3)
        for ey = 1:Ky
                for ex = 1:Kx
                        id(ex,ey) = ex + (ey-1)*(Kx+1) # index function
                        id1 = id(ex,ey);
                        id2 = id(ex+1,ey);
                        id3 = id(ex+1,ey+1);
                        id4 = id(ex,ey+1);
                        EToV[2*sk-1,:] = [id1 id2 id3];
                        EToV[2*sk,:] = [id3 id4 id1];
                        sk += 1
                end
        end
        return (VX[:],VY[:],EToV)
end

function uniform_tri_mesh(Kx)
        return uniform_tri_mesh(Kx,Kx)
end

function tri_face_vertices()
        return [1,2],[2,3],[3,1]
end

end

max_stars_repo_name: jlchan/ESDG-CNS

max_issues_repo_name: jlchan/ESDG-CNS

max_forks_repo_name: jlchan/ESDG-CNS


text:
using Documenter, Persa

makedocs(
    modules = [Persa],
    analytics = "UA-128580038-1",
    sitename = "Persa.jl"
)

# Only build plots in travis if we are deploying
# And dont install the dependencies unless we are deploying
function myDeps()
    if get(ENV, "TRAVIS", "") != ""
        println("Installing deploy dependencies")
        run(`pip install --user pygments mkdocs`)
    end
end

deploydocs(
    repo = "github.com/JuliaRecsys/Persa.jl.git",
    latest = "master",
    julia = "1.0",
    deps = myDeps
)


text:
# # Padua transform
# This demonstrates the Padua transform and inverse transform,
# explaining precisely the normalization and points

using FastTransforms

# We define the Padua points and extract Cartesian components:
N = 15
pts = paduapoints(N)
x = pts[:,1]
y = pts[:,2];

# We take the Padua transform of the function:
f = (x,y) -> exp(x + cos(y))
f̌ = paduatransform(f.(x , y));

# and use the coefficients to create an approximation to the function $f$:
f̃ = (x,y) -> begin
    j = 1
    ret = 0.0
    for n in 0:N, k in 0:n
        ret += f̌[j]*cos((n-k)*acos(x)) * cos(k*acos(y))
        j += 1
    end
    ret
end

# At a particular point, is the function well-approximated?
f̃(0.1,0.2) ≈ f(0.1,0.2)

# Does the inverse transform bring us back to the grid?
ipaduatransform(f̌) ≈ f̃.(x,y)


text:
# Reference/Adapted from:
# https://www.mortality.org/Public/HMD_4th_Symposium/Pascariu_poster.pdf
# https://github.com/mpascariu/MortalityLaws/blob/master/R/MortalityLaw_models.R

abstract type ParametricMortality end
"""
    Makeham(;a,b,c)

Construct a mortality model following Makeham's law.

``
\\mathrm{hazard} \\left( {\\rm age} \\right) =  ae^{bx} + c
``

Default args:

    a = 0.0002
    b = 0.13
    c = 0.001

"""
Base.@kwdef struct Makeham <: ParametricMortality
    a = 0.0002
    b = 0.13
    c = 0.001
end

"""
    hazard(model,age)

The force of mortality at `age`. More precisely: the ratio of the probability of failure/death to the survival function.
"""
function hazard(m::Makeham,age)
    @unpack a,b,c = m
    return a*exp(b*age) + c
end

"""
    cumhazard(model,age)

The cumulative force of mortality at `age`. More precisely: the ratio of the cumulative probability of failure/death to the survival function.
"""
function cumhazard(m::Makeham,age)
    @unpack a,b,c = m
    return a / b * (exp(b*age) - 1) + age * c
end

survival(m::Makeham,age) = exp(-cumhazard(m,age))

"""
    Gompertz(;a,b)

Construct a mortality model following Gompertz' law of mortality.

``
\\mathrm{hazard} \\left( {\\rm age} \\right) =  ae^{bx}
``

This is a special case of Makeham's law and will `Makeham` model where `c=0`.

Default args:

    a = 0.0002
    b = 0.13

"""
function Gompertz(;a=0.0002, b=0.13)
    return Makeham(a=a, b=b, c=0)
end

"""
    InverseGompertz(;a,b,c)

Construct a mortality model following InverseGompertz's law.

```math
\\begin{aligned}
\\mathrm{hazard} \\left( {\\rm age} \\right) &= \\frac{1}{\\sigma}e^\\frac{age-m}{\\sigma}/e^{e^\\frac{-(age-m)}{\\sigma}-1}``
\\\\
\\mathrm{survival} \\left( {\\rm age} \\right) &= \\frac{1 - e^{ - e^{\\frac{ - \\left( {\\rm age} - m \\right)}{\\sigma}}}}{1 - e^{ - e^{\\frac{m}{\\sigma}}}}``
\\end{aligned}
```

Default args:

    m = 49
    σ = 7.7

"""
Base.@kwdef struct InverseGompertz <: ParametricMortality
    m = 49
    σ = 7.7
end

function hazard(m::InverseGompertz,age)
    @unpack m,σ = m
    return 1 / σ * exp(-(age - m)/σ) / (exp(exp(-(age - m)/σ)) - 1)
end

cumhazard(m::InverseGompertz,age) = -log(survival(m,age))

function survival(m::InverseGompertz,age)
    @unpack m,σ = m
    return (1 - exp(-exp(-(age - m)/σ))) / (1 - exp(-exp(m/σ)))
end

"""
    Opperman(;a,b,c)

Construct a mortality model following Opperman's law of mortality.

``
\\mathrm{hazard} \\left( {\\rm age} \\right) = \\frac{a}{\\sqrt{age}} + b +c\\sqrt[3]{age}
``

Default args:

    a = 0.04
    b = 0.0004
    c = 0.001
"""
Base.@kwdef struct Opperman <: ParametricMortality
    a = 0.04
    b = 0.0004
    c = 0.001
end

function hazard(m::Opperman,age)
    @unpack a,b,c = m
    return max(a / √(age+1) - b + c * √(age+1),0.0)
end

"""
    Thiele(;a,b,c,d,e,f,g)

Construct a mortality model following Opperman's law of mortality.

```math
\\begin{aligned}
\\mu_1 &= a \\cdot e^{\\left(  - b \\right) \\cdot {\\rm age}}
\\\\
\\mu_2 &= c \\cdot e^{-0.5 \\cdot d \\cdot \\left( {\\rm age} - e \\right)^{2}}
\\\\
\\mu_3 &= f \\cdot e^{g \\cdot {\\rm age}}
\\\\
\\mathrm{hazard} \\left( {\\rm age} \\right) &= \\begin{cases}
\\mu_1 + \\mu_3 & \\text{if } \\left( {\\rm age} = 0 \\right)\\\\
\\mu_1 + \\mu_2 + \\mu_3 & \\text{otherwise}
\\end{cases}
\\end{aligned}
```
Default args:

    a = 0.02474
    b = 0.3
    c = 0.004
    d = 0.5
    e = 25
    f = 0.0001
    g = 0.13
"""
Base.@kwdef struct Thiele <: ParametricMortality
    a = 0.02474
    b = 0.3
    c = 0.004
    d = 0.5
    e = 25
    f = 0.0001
    g = 0.13
end

function hazard(m::Thiele,age)
    @unpack a,b,c,d,e,f,g = m
    μ₁ = a * exp(-b * age)
    μ₂ = c * exp(-0.5 * d * (age - e)^2)
    μ₃ = f * exp(g * age)

    if age == 0
        return μ₁ + μ₃
    else
        return  μ₁ + μ₂ + μ₃
    end
end

"""
    Wittstein(;a,b,m,n)

Construct a mortality model following Wittstein's law of mortality.

``\\mathrm{hazard} \\left( {\\rm age} \\right) = \\frac{1}{b} \\cdot a^{ - \\left( b \\cdot {\\rm age} \\right)^{n}} + a^{ - \\left( m - {\\rm age} \\right)^{n}}``

Default args:

    a = 1.5
    b = 1.
    n = 0.5
    m = 100

"""
Base.@kwdef struct Wittstein <: ParametricMortality
    a = 1.5
    b = 1.
    n = 0.5
    m = 100
end

function hazard(m::Wittstein,age)
    @unpack a,b,m,n = m
    return (1/b) * a ^ -((b * age) ^ n) + a^ -((m -  age) ^ n)
end

"""
    Weibull(;m,σ)

Construct a mortality model following Weibull's law of mortality.

Note that if σ > m, then the mode of the density is 0 and hx is a non-increasing function of x, while if σ < m, then the mode is greater than 0 and hx is an increasing function.
 - `m >0` is a measure of location
 - `σ >0` is measure of dispersion

```math
\\begin{aligned}
\\mathrm{hazard} \\left( {\\rm age} \\right) = \\frac{1}{\\sigma} \\cdot \\left( \\frac{{\\rm age}}{m} \\right)^{\\frac{m}{\\sigma} - 1}
\\\\
\\mathrm{cumhazard} \\left( {\\rm age} \\right) = \\left( \\frac{{\\rm age}}{m} \\right)^{\\frac{m}{\\sigma}}
\\\\
\\mathrm{survival} \\left( {\\rm age} \\right) =  e^{ - \\mathrm{cumhazard} \\left( m, {\\rm age} \\right)}
\\end{aligned}
```

 Default args:

    m = 1
    σ = 2
"""
Base.@kwdef struct Weibull <: ParametricMortality
    m = 1
    σ = 2
end

function hazard(m::Weibull,age)
    @unpack m,σ =m
    if age == 0
        return 1.0
    else
        return 1 / σ * (age / m)^(m / σ - 1)
    end
end

function cumhazard(m::Weibull,age)
    @unpack m,σ =m
    return (age / m) ^ (m / σ)
end

function survival(m::Weibull,age)
    return exp(-cumhazard(m,age))
end

"""
    InverseWeibull(;m,σ)

Construct a mortality model following Weibull's law of mortality.

The Inverse-Weibull proves useful for modelling the childhood and teenage years, because the logarithm of h(x) is a concave function.
 - `m >0` is a measure of location
 - `σ >0` is measure of dispersion

```math
\\begin{aligned}
\\mathrm{hazard} \\left( {\\rm age} \\right) &= \\frac{\\frac{1}{\\sigma} \\cdot \\left( \\frac{{\\rm age}}{m} \\right)^{\\frac{ - m}{\\sigma} - 1}}{e^{\\left( \\frac{{\\rm age}}{m} \\right)^{\\frac{ - m}{\\sigma}}} - 1}
\\\\
\\mathrm{cumhazard}\\left( {\\rm age} \\right) &=  - \\log\\left( 1 - e^{ - \\left( \\frac{{\\rm age}}{m} \\right)^{\\frac{ - m}{\\sigma}}} \\right)
\\\\
\\mathrm{survival}\\left( {\\rm age} \\right) &=  e^{ - \\mathrm{cumhazard}\\left( m, {\\rm age} \\right)}
\\end{aligned}
```

 Default args:

    m = 5
    σ = 10

"""
Base.@kwdef struct InverseWeibull <: ParametricMortality
    m = 5
    σ = 10
end

function hazard(m::InverseWeibull,age)
    @unpack m,σ = m
    return (1/σ) * (age/m)^(-m/σ - 1) / (exp((age/m)^(-m/σ)) - 1)
end

function cumhazard(m::InverseWeibull,age)
    @unpack m,σ = m
    return -log(1 - exp(-(age/m)^(-m/σ)))
end

function survival(m::InverseWeibull,age)
    return exp(-cumhazard(m,age))
end

"""
    Perks(;a,b,c,d)

Construct a mortality model following Perks' law of mortality.

``
\\mathrm{hazard} \\left( {\\rm age} \\right) = \\frac{a + b \\cdot c^{{\\rm age}}}{b \\cdot c^{ - {\\rm age}} + 1 + d \\cdot c^{{\\rm age}}}
``

Default args:

    a = 0.002
    b = 0.13
    c = 0.01
    d = 0.01
"""
Base.@kwdef struct Perks <: ParametricMortality
    a = 0.002
    b = 0.13
    c = 0.01
    d = 0.01
end

function hazard(m::Perks,age)
    @unpack a,b,c,d = m
    return (a + b*c^age) / (b*(c^-age) + 1 + d*c^age)
end

"""
    VanderMaen(;a,b,c,i,n)

Construct a mortality model following VanderMaen's law of mortality.

``
\\mathrm{hazard} \\left( {\\rm age} \\right) = a + b \\cdot {\\rm age} + c \\cdot {\\rm age}^{2} + \\frac{i}{n - {\\rm age}}
``

Default args:

    a = 0.01
    b = 1
    c = 0.01
    i = 100
    n = 200
"""
Base.@kwdef struct VanderMaen <: ParametricMortality
    a = 0.01
    b = 1.
    c = 0.01
    i = 100.
    n = 200.
end

function hazard(m::VanderMaen,age)
    @unpack a,b,c,i,n = m
    return a + b*age + c*(age^2) + i/(n - age)
end

"""
    VanderMaen2(;a,b,i,n)

Construct a mortality model following VanderMaen2's law of mortality.

``
\\mathrm{hazard} \\left( {\\rm age} \\right) = a + b \\cdot {\\rm age} + \\frac{i}{n - {\\rm age}}
``

Default args:

    a = 0.01
    b = 1
    i = 100
    n = 200

"""
Base.@kwdef struct VanderMaen2 <: ParametricMortality
    a = 0.01
    b = 1.
    i = 100.
    n = 200.
end

function hazard(m::VanderMaen2,age)
    @unpack a,b,i,n = m
    return a + b * age + i/(n - age)
end

"""
    StrehlerMildvan(;k,v₀,b,d)

Construct a mortality model following StrehlerMildvan's law of mortality.

``
\\mathrm{hazard} \\left( {\\rm age} \\right) = k \\cdot e^{\\frac{\\left(  - v_0 \\right) \\cdot \\left( 1 - b \\cdot {\\rm age} \\right)}{d}}
``

Default args:

    k   = 0.01
    v₀  = 2.5
    b   = 0.2
    d   = 6.0

"""
Base.@kwdef struct StrehlerMildvan <: ParametricMortality
    k   = 0.01
    v₀  = 2.5
    b   = 0.2
    d   = 6.0
end

function hazard(m::StrehlerMildvan,age)
    @unpack k,v₀,b,d = m
    return  k * exp(-v₀ * (1 - b * age) / d)
end

"""
    Beard(;a,b,k)

Construct a mortality model following Beard's law of mortality.

``
\\mathrm{hazard} \\left( {\\rm age} \\right) = \\frac{a \\cdot e^{b \\cdot {\\rm age}}}{1 + k \\cdot a \\cdot e^{b \\cdot {\\rm age}}}
``

Default args:

    a = 0.002
    b = 0.13
    k = 1.
"""
Base.@kwdef struct Beard <: ParametricMortality
    a = 0.002
    b = 0.13
    k = 1.
end

function hazard(m::Beard,age)
    @unpack a,b,k = m
    return  a * exp(b*age) / (1 + k * a * exp(b*age))
end

"""
    MakehamBeard(;a,b,c,k)

Construct a mortality model following MakehamBeard's law of mortality.

``
\\mathrm{hazard} \\left( {\\rm age} \\right) =\\left( {\\rm age} \\right) = \\frac{a \\cdot e^{b \\cdot {\\rm age}}}{1 + k \\cdot a \\cdot e^{b \\cdot {\\rm age}}} + c
``

Default args:

    a = 0.002
    b = 0.13
    c = 0.01
    k = 1.
"""
Base.@kwdef struct MakehamBeard <: ParametricMortality
    a = 0.002
    b = 0.13
    c = 0.01
    k = 1.
end

function hazard(m::MakehamBeard,age)
    @unpack a,b,c,k = m
    return  a * exp(b*age) / (1 + k * a * exp(b*age)) + c
end

"""
    Quadratic(;a,b,c)

Construct a mortality model following Quadratic law of mortality.

``
\\mathrm{hazard} \\left( {\\rm age} \\right) = a + b \\cdot {\\rm age} + c \\cdot {\\rm age}^{2}
``

Default args:

    a = 0.01
    b = 1.
    c = 0.01
"""
Base.@kwdef struct Quadratic <: ParametricMortality
    a = 0.01
    b = 1.
    c = 0.01
end

function hazard(m::Quadratic,age)
    @unpack a,b,c = m
    return  a + b * age + c * age^2
end

"""
    GammaGompertz(;a,b,γ)

Construct a mortality model following GammaGompertz law of mortality.

``
\\mathrm{hazard} \\left( {\\rm age} \\right) = \\frac{a \\cdot e^{b \\cdot {\\rm age}}}{1 + \\frac{a \\cdot \\gamma}{b} \\cdot \\left( e^{b \\cdot {\\rm age}} - 1 \\right)}
``

Default args:

    a = 0.002
    b = 0.13
    γ = 1
"""
Base.@kwdef struct GammaGompertz <: ParametricMortality
    a = 0.002
    b = 0.13
    γ = 1
end

function hazard(m::GammaGompertz,age)
    @unpack a,b,γ = m
    return  (a * exp(b * age)) / (1 + ( a * γ / b) * (exp(b * age) - 1))
end

"""
    Siler(;a,b,c,d,e)

Construct a mortality model following Siler law of mortality.

``
\\mathrm{hazard} \\left( {\\rm age} \\right) = a \\cdot e^{\\left(  - b \\right) \\cdot {\\rm age}} + c + d \\cdot e^{e \\cdot {\\rm age}}
``

Default args:

    a = 0.0002
    b = 0.13
    c = 0.001
    d = 0.001
    e = 0.013
"""
Base.@kwdef struct Siler <: ParametricMortality
    a = 0.0002
    b = 0.13
    c = 0.001
    d = 0.001
    e = 0.013
end

function hazard(m::Siler,age)
    @unpack a,b,c,d,e = m
    return  a * exp(-b* age) + c + d * exp(e * age)
end

"""
    HeligmanPollard(;a,b,c,d,e,f,g,h)

Construct a mortality model following HeligmanPollard law of mortality with 8 parameters.

``
\\mathrm{hazard} \\left( {\\rm age} \\right) = a \\cdot e^{\\left(  - b \\right) \\cdot {\\rm age}} + c + d \\cdot e^{e \\cdot {\\rm age}}
``

Default args:

    a = 0.0002
    b = 0.13
    c = 0.001
    d = 0.001
    e = 0.013
"""
Base.@kwdef struct HeligmanPollard <: ParametricMortality
    a = .0005
    b = .004
    c = .08
    d = .001
    e = 10
    f = 17
    g = .00005
    h = 1.1
end

function hazard(m::HeligmanPollard,age)
    @unpack a,b,c,d,e,f,g,h = m
    μ₁ = a^((age + b)^c) + g * h^age
    μ₂ = d * exp(-e * (log(age/f))^2)
    η = age == 0 ?  μ₁ :  μ₁ + μ₂
    return  η / (1 + η)
end

"""
    HeligmanPollard2(;a,b,c,d,e,f,g,h)

Construct a mortality model following HeligmanPollard (alternate) law of mortality with 8 parameters.

```math
\\begin{aligned}
\\mu_1 &= a^{\\left( {\\rm age} + b \\right)^{c}} + \\frac{g \\cdot h^{{\\rm age}}}{1 + g \\cdot h^{{\\rm age}}}
\\\\
\\mu_2 &= d \\cdot e^{\\left(  - e \\right) \\cdot \\left( \\log\\left( \\frac{{\\rm age}}{f} \\right) \\right)^{2}}
\\\\
\\mathrm{hazard}\\left( {\\rm age} \\right) &= \\begin{cases}
\\mu_1 & \\text{if } \\left( {\\rm age} = 0 \\right)\\\\
\\mu_1 + \\mu_2 & \\text{otherwise}
\\end{cases}
\\end{aligned}
```

Default args:

    a = .0005
    b = .004
    c = .08
    d = .001
    e = 10
    f = 17
    g = .00005
    h = 1.1
"""
Base.@kwdef struct HeligmanPollard2 <: ParametricMortality
    a = 0.0005
    b = 0.004
    c = 0.08
    d = 0.001
    e = 10.
    f = 17.
    g = 0.00005
    h = 1.1
end

function hazard(m::HeligmanPollard2,age)
    @unpack a,b,c,d,e,f,g,h = m
    μ₁ = a^((age + b)^c) + (g * h^age) / (1 + g * h ^ age)
    μ₂ = d * exp(-e * (log(age/f))^2)
    return age == 0 ?  μ₁ :  μ₁ + μ₂
end

"""
    HeligmanPollard3(;a,b,c,d,e,f,g,h,k)

Construct a mortality model following HeligmanPollard (alternate) law of mortality with 9 parameters.

```math
\\begin{aligned}
\\mu_1 &= a^{\\left( {\\rm age} + b \\right)^{c}} + \\frac{g \\cdot h^{{\\rm age}}}{1 + k \\cdot g \\cdot h^{{\\rm age}}}
\\\\
\\mu_2 &= d \\cdot e^{\\left(  - e \\right) \\cdot \\left( \\log\\left( \\frac{{\\rm age}}{f} \\right) \\right)^{2}}
\\\\
\\mathrm{hazard}\\left( {\\rm age} \\right) &= \\begin{cases}
\\mu_1 & \\text{if } \\left( {\\rm age} = 0 \\right)\\\\
\\mu_1 + \\mu_2 & \\text{otherwise}
\\end{cases}
\\end{aligned}
```

Default args:

    a = .0005
    b = .004
    c = .08
    d = .001
    e = 10
    f = 17
    g = .00005
    h = 1.1
    k= 1.
"""
Base.@kwdef struct HeligmanPollard3 <: ParametricMortality
    a = .0005
    b = .004
    c = .08
    d = .001
    e = 10
    f = 17
    g = .00005
    h = 1.1
    k = 1.
end

function hazard(m::HeligmanPollard3,age)
    @unpack a,b,c,d,e,f,g,h,k = m
    μ₁ = a^((age + b)^c) + (g * h^age) / (1 + k * g * h ^ age)
    μ₂ = d * exp(-e * (log(age/f))^2)
    return age == 0 ?  μ₁ :  μ₁ + μ₂
end

"""
    HeligmanPollard4(;a,b,c,d,e,f,g,h,k)

Construct a mortality model following HeligmanPollard (alternate) law of mortality with 9 parameters.

```math
\\begin{aligned}
\\mu_1 &= a^{\\left( {\\rm age} + b \\right)^{c}} + \\frac{g \\cdot h^{{\\rm age}^{k}}}{1 + g \\cdot h^{{\\rm age}^{k}}}
\\\\
\\mu_2 &= d \\cdot e^{\\left(  - e \\right) \\cdot \\left( \\log\\left( \\frac{{\\rm age}}{f} \\right) \\right)^{2}}
\\\\
\\mathrm{hazard}\\left( {\\rm age} \\right) &= \\begin{cases}
\\mu_1 & \\text{if } \\left( {\\rm age} = 0 \\right)\\\\
\\mu_1 + \\mu_2 & \\text{otherwise}
\\end{cases}
\\end{aligned}
```

Default args:

    a = .0005
    b = .004
    c = .08
    d = .001
    e = 10
    f = 17
    g = .00005
    h = 1.1
    k= 1.
"""
Base.@kwdef struct HeligmanPollard4 <: ParametricMortality
    a = .0005
    b = .004
    c = .08
    d = .001
    e = 10
    f = 17
    g = .00005
    h = 1.1
    k = 1.
end

function hazard(m::HeligmanPollard4,age)
    @unpack a,b,c,d,e,f,g,h,k = m
    μ₁ = a^((age + b)^c) + (g * h^(age ^ k)) / (1 + g * h ^ (age ^ k))
    μ₂ = d * exp(-e * (log(age/f))^2)
    return age == 0 ?  μ₁ :  μ₁ + μ₂
end

"""
    RogersPlanck(;a₀, a₁, a₂, a₃, a, b, c, d, u)

Construct a mortality model following RogersPlanck law of mortality.

``
\\mathrm{hazard}\\left( {\\rm age} \\right) = a_0 + a_1 \\cdot e^{\\left(  - a \\right) \\cdot {\\rm age}} + a_2 \\cdot e^{b \\cdot \\left( {\\rm age} - u \\right) - e^{\\left(  - c \\right) \\cdot \\left( {\\rm age} - u \\right)}} + a_3 \\cdot e^{d \\cdot {\\rm age}}
``

Default args:

    a₀ = 0.0001
    a₁ = 0.02
    a₂ = 0.001
    a₃ = 0.0001
    a  = 2.
    b  = 0.001
    c  = 100.
    d  = 0.1
    u  = 0.33

"""
Base.@kwdef struct RogersPlanck <: ParametricMortality
    a₀ = 0.0001
    a₁ = 0.02
    a₂ = 0.001
    a₃ = 0.0001
    a  = 2.
    b  = 0.001
    c  = 100.
    d  = 0.1
    u  = 0.33
end

function hazard(m::RogersPlanck,age)
    @unpack a₀,a₁,a₂,a₃,a,b,c,d,u = m
    return  a₀ + a₁ * exp(-a * age) + a₂ * exp(b*(age - u) - exp(-c*(age - u))) + a₃*exp(d*age)
end

"""
    Martinelle(;a,b,c,d,k)

Construct a mortality model following Martinelle's law of mortality.

``
\\mathrm{hazard}\\left( {\\rm age} \\right) = \\frac{a \\cdot e^{b \\cdot {\\rm age}} + c}{1 + d \\cdot e^{b \\cdot {\\rm age}}} + k \\cdot e^{b \\cdot {\\rm age}}
``

Default args:

    a = 0.001
    b = 0.13
    c = 0.001
    d = 0.1
    k = 0.001
"""
Base.@kwdef struct Martinelle <: ParametricMortality
    a = 0.001
    b = 0.13
    c = 0.001
    d = 0.1
    k = 0.001
end

function hazard(m::Martinelle,age)
    @unpack a,b,c,d,k = m
    return  (a*exp(b*age) + c) / (1 + d*exp(b * age)) + k*exp(b * age)
end

"""
    Kostaki(;a,b,c,d,e1,e2,f,g,h)

Construct a mortality model following Kostaki's law of mortality. A nine-parameter adaptation of `HeligmanPollard`.

```math
\\begin{aligned}
\\mu_1 &= a^{\\left( {\\rm age} + b \\right)^{c}} + g \\cdot h^{{\\rm age}}
\\\\
\\mu_2 &= \\begin{cases}
d \\cdot e^{ - \\left( e1 \\cdot \\log\\left( \\frac{{\\rm age}}{f} \\right) \\right)^{2}} & \\text{if } \\left( {\\rm age} \\leq f \\right)\\\\
d \\cdot e^{ - \\left( e2 \\cdot \\log\\left( \\frac{{\\rm age}}{f} \\right) \\right)^{2}} & \\text{otherwise}
\\end{cases}
\\\\
\\eta &= \\begin{cases}
\\mu_1 & \\text{if } \\left( {\\rm age} = 0 \\right)\\\\
\\mu_1 + \\mu_2 & \\text{otherwise}
\\end{cases}
\\\\
\\mathrm{hazard}\\left( {\\rm age} \\right) &= \\frac{\\eta}{1 + \\eta}

\\end{aligned}
```

Default args:

    a = 0.0005
    b = 0.01
    c = 0.10
    d = 0.001
    e1 = 3.
    e2 = 0.1
    f = 25.
    g = .00005
    h = 1.1

> Kostaki, A. (1992). A nine‐parameter version of the Heligman‐Pollard formula. Mathematical Population Studies, 3(4), 277–288. doi:10.1080/08898489209525346
"""
Base.@kwdef struct Kostaki <: ParametricMortality
    a = 0.0005
    b = 0.01
    c = 0.10
    d = 0.001
    e1 = 3.
    e2 = 0.1
    f = 25.
    g = .00005
    h = 1.1
end

function hazard(m::Kostaki,age)
    @unpack a,b,c,d,e1,e2,f,g,h = m
    μ₁ = a^((age + b)^c) + g*h^age
    if age <= f
        μ₂ =d *exp(-(e1*log(age/f))^2)
    else
        μ₂ =d *exp(-(e2*log(age/f))^2)
    end

    η = age == 0 ? μ₁ : μ₁ + μ₂

    return η / (1+η)
end

"""
    Kannisto(;a,b)

Construct a mortality model following Kannisto's law of mortality.

```math
\\begin{aligned}
\\mathrm{hazard}\\left( {\\rm age} \\right) &= \\frac{a \\cdot e^{b \\cdot {\\rm age}}}{1 + a \\cdot e^{b \\cdot {\\rm age}}}
\\\\
\\mathrm{cumhazard}\\left( {\\rm age} \\right) &= 1/a * log((1 + b*exp(b*age)) / (1 + a))
\\\\
\\mathrm{survival}\\left( {\\rm age} \\right) &= e^{ - \\mathrm{cumhazard}\\left( m, {\\rm age} \\right)}
\\end{aligned}
```

Default args:

    a = 0.5
    b = 0.13
"""
Base.@kwdef struct Kannisto <: ParametricMortality
    a = 0.5
    b = 0.13
end

function hazard(m::Kannisto,age)
    @unpack a,b = m
    return  a * exp(b * age) / (1 + a * exp(b*age))
end

function cumhazard(m::Kannisto,age)
    @unpack a,b = m
    return  1/a * log((1 + b*exp(b*age)) / (1 + a))
end

function  survival(m::Kannisto,age)
    return exp(-cumhazard(m,age))
end

"""
    KannistoMakeham(;a,b,c)

Construct a mortality model following KannistoMakeham's law of mortality.

``
\\mathrm{hazard}\\left( {\\rm age} \\right) = \\frac{a \\cdot e^{b \\cdot {\\rm age}}}{1 + a \\cdot e^{b \\cdot {\\rm age}}} + c
``

Default args:

    a = 0.5
    b = 0.13
    c = 0.001
"""
Base.@kwdef struct KannistoMakeham <: ParametricMortality
    a = 0.5
    b = 0.13
    c = 0.001
end

function hazard(m::KannistoMakeham,age)
    @unpack a,b,c = m
    return  a * exp(b * age) / (1 + a * exp(b*age)) + c
end

### Generic Functions

"""
    μ(;m::ParametricMortality,age)

``\\mu_x``: Return the force of mortality at the given age.
"""
function μ(m::ParametricMortality, age)
    return hazard(m,age)
end

# # use the integral to calculate the one-year survival
# function survival(m::ParametricMortality, from_age, to_age)
#     if from_age == to_age
#         return 1.0
#     else
#         return exp(-quadgk(age->μ(m, age), from_age, to_age)[1])
#     end
# end
# survival(m::ParametricMortality,to_age) = survival(m, 0, to_age)

survival(m::ParametricMortality,to_age) = exp(-quadgk(age->μ(m, age), 0, to_age)[1])
survival(m::ParametricMortality,from,to) = survival(m,to) / survival(m,from)
decrement(m::ParametricMortality,from_age,to_age) = 1 - survival(m, from_age, to_age)
decrement(m::ParametricMortality,to_age) = 1 - survival(m, to_age)

(m::ParametricMortality)(x) = μ(m, x)
Base.getindex(m::ParametricMortality,x) = m(x)
Base.broadcastable(pm::ParametricMortality) = Ref(pm)