text:
@testset "ConvVAE" begin
	using DrWatson
	@quickactivate
	using Test
	using GenerativeAD
	import StatsBase: fit!, predict
	using ConditionalDists, GenerativeModels, Distributions
	using EvalMetrics
	using IPMeasures
	using Flux

	# toy example
	M = (8,8,1)
	M2 = (4,8,1)
	N1 = 100
	N2 = 50
	X = cat(randn(Float32, M2..., N1), ones(Float32, M2..., N1), dims=1)
	Y = cat(randn(Float32, M2..., N2), ones(Float32, M2..., N2), dims=1)
	Z = cat(ones(Float32, M2..., N2), randn(Float32, M2..., N2), dims=1)
	data = ((X,),(Y,zeros(size(Y,4))))

	# first test conv_encoder and decoder
	idim = size(X)[1:3]
	zdim = 4
	kernelsizes = (3,5)
	channels = (2,4)
	scalings = (2,1)
	encoder = GenerativeAD.Models.conv_encoder(idim, zdim, kernelsizes, channels, scalings)
	@test size(encoder(X)) == (zdim, N1)
	decoder = GenerativeAD.Models.conv_decoder(idim, zdim, reverse(kernelsizes),
		reverse(channels), reverse(scalings))
	z = randn(zdim, N1)
	@test size(decoder(z)) == (idim..., N1)
	hdim = 32
	decoder = GenerativeAD.Models.conv_decoder(idim, zdim, reverse(kernelsizes), reverse(channels),
		reverse(scalings), activation="tanh", vec_output=true, vec_output_dim=hdim)
	@test size(decoder(z)) == (hdim, N1)
	function basic_convergence_test(model, loss, scoref; kwargs...)
		history, iterations, model = fit!(model, data, loss; patience = 100, max_train_time=600,
			kwargs...)
		scores = map(x->scoref(model, x), (X,Y,Z))
		@test mean(scores[1]) < mean(scores[3])
		@test mean(scores[2]) < mean(scores[3])
		model
	end
	# test convolutional vae
	parameters = (idim=idim, zdim=zdim, activation="swish", kernelsizes=(5,3),
		channels=(4,8), scalings=(2,2))
	model = GenerativeAD.Models.conv_vae_constructor(;parameters...)
	loss(m,x) = -elbo(m,x)
	loss(m, x, batchsize::Int) =
		mean(map(y->loss(m,y), Flux.Data.DataLoader(x, batchsize=batchsize)))
	model = basic_convergence_test(model, loss, GenerativeAD.Models.reconstruction_score)
	@test size(GenerativeAD.Models.reconstruct(model, X)) == size(X)
	@test size(GenerativeAD.Models.generate(model, 10, idim)) == (idim..., 10)
	# gpu
	gmodel = GenerativeAD.Models.conv_vae_constructor(;parameters...) |> gpu
	gX = X |> gpu
	# the Array is there since gpu computations dont work on views created by MLDataPattern iterators
	gloss(m,x) = -elbo(m,gpu(Array(x)))
	gloss(m, x, batchsize::Int) =
		mean(map(y->gloss(m,y), Flux.Data.DataLoader(x, batchsize=batchsize)))
	gscoref(model, x) = GenerativeAD.Models.reconstruction_score(model, gpu(x))
	gmodel = basic_convergence_test(gmodel, gloss, gscoref; parameters...)
end

ext: jl


text:
module UnmatchedBrackets

"!!var(unmatched_brackets:This will throw an error (Unmatched Brackets)"
unmatched_brackets = ()

end


text:
module LifeContingencies

using ActuaryUtilities
using MortalityTables
using Transducers
using Dates
using Yields

const mt = MortalityTables

export LifeContingency,
    Insurance, AnnuityDue, AnnuityImmediate,
    APV,
    SingleLife, Frasier, JointLife,
    LastSurvivor,
    survival,
    reserve_premium_net,
    insurance,
    annuity_due,
    annuity_immediate,
    premium_net,
    omega,
    survival,
    discount,
    benefit,
    probability,
    cashflows,
    cashflows,
    timepoints,
    present_value

# 'actuarial objects' that combine multiple forms of decrements (lapse, interest, death, etc)
abstract type Life end

"""
    struct SingleLife
        mortality
        issue_age::Int
        alive::Bool
        fractional_assump::MortalityTables.DeathDistribution
    end

A `Life` object containing the necessary assumptions for contingent maths related to a single life. Use with a `LifeContingency` to do many actuarial present value calculations.

Keyword arguments:
- `mortality` pass a mortality vector, which is an array of applicable mortality rates indexed by attained age
- `issue_age` is the assumed issue age for the `SingleLife` and is the basis of many contingency calculations.
- `alive` Default value is `true`. Useful for joint insurances with different status on the lives insured.
- `fractional_assump`. Default value is `Uniform()`. This is a `DeathDistribution` from the `MortalityTables.jl` package and is the assumption to use for non-integer ages/times.

# Examples
    using MortalityTables
    mortality = MortalityTables.table("2001 VBT Residual Standard Select and Ultimate - Male Nonsmoker, ANB")

    SingleLife(
        mort       = mort.select[30],
        issue_age  = 30
    )
"""
struct SingleLife{M,D} <: Life
    mortality::M
    issue_age::Int
    alive::Bool
    fractional_assump::D
end

function SingleLife(; mortality, issue_age = nothing, alive = true, fractional_assump = mt.Uniform())
    return SingleLife(mortality; issue_age, alive, fractional_assump)
end

function SingleLife(mortality; issue_age = nothing, alive = true, fractional_assump = mt.Uniform())
    if isnothing(issue_age)
        issue_age = firstindex(mortality)
    end

    if !(eltype(mortality) <: Real)
        # most likely case is that mortality is an array of vectors
        # use issue age to select the right one (assuming indexed with issue age
        return SingleLife(mortality[issue_age], issue_age, alive, fractional_assump)
    else
        return SingleLife(mortality, issue_age, alive, fractional_assump)
    end

end

"""
    JointAssumption()

An abstract type representing the different assumed relationship between the survival of the lives on a JointLife. Available options to use include:
- `Frasier()`
"""
abstract type JointAssumption end

"""
    Frasier()

The assumption of independent lives in a joint life calculation.
Is a subtype of `JointAssumption`.
"""
struct Frasier <: JointAssumption end

"""
    Contingency()

An abstract type representing the different triggers for contingent benefits. Available options to use include:
- `LastSurvivor()`
"""
abstract type Contingency end

"""
    LastSurvivor()
The contingency whereupon benefits are payable upon both lives passing.
Is a subtype of `Contingency`
"""
struct LastSurvivor <: Contingency end

# TODO: Not Implemented
# """
#     FirstToDie()
# The contingency whereupon benefits are payable upon the first life passing.

# Is a subtype of `Contingency`
# """
# struct FirstToDie <: Contingency end

"""
    struct JointLife
        lives
        contingency
        joint_assumption
    end

    A `Life` object containing the necessary assumptions for contingent maths related to a joint life insurance. Use with a `LifeContingency` to do many actuarial present value calculations.

Keyword arguments:
- `lives` is a tuple of two `SingleLife`s
- `contingency` default is `LastSurvivor()`. It is the trigger for contingent benefits. See `?Contingency`.
- `joint_assumption` Default value is `Frasier()`. It is the assumed relationship between the mortality of the two lives. See `?JointAssumption`.

# Examples
    using MortalityTables
    mortality = MortalityTables.table("2001 VBT Residual Standard Select and Ultimate - Male Nonsmoker, ANB")

    l1 = SingleLife(
        mortality       = mortality.select[30],
        issue_age  = 30
    )
    l2 = SingleLife(
        mortality       = mortality.select[30],
        issue_age  = 30
    )

    jl = JointLife(
        lives = (l1,l2),
        contingency = LastSurvivor(),
        joint_assumption = Frasier()
    )
"""
Base.@kwdef struct JointLife{C<:Contingency,J<:JointAssumption} <: Life
    lives::Tuple{SingleLife,SingleLife}
    contingency::C = LastSurvivor()
    joint_assumption::J = Frasier()
end

"""
    struct LifeContingency
        life::Life
"""
struct LifeContingency{L,Y}
    life::L
    int::Y
end

Base.broadcastable(lc::LifeContingency) = Ref(lc)

"""
    omega(lc::LifeContingency)
    omega(l::Life)
    omega(i::InterestRate)

# `Life`s and `LifeContingency`s

Returns the last defined time_period for both the interest rate and mortality table.
Note that this is *different* than calling `omega` on a `MortalityTable`, which will give you the last `attained_age`.

Example: if the `LifeContingency` has issue age 60, and the last defined attained age for the `MortalityTable` is 100, then `omega` of the `MortalityTable` will be `100` and `omega` of the
`LifeContingency` will be `40`.

# `InterestRate`s

The last period that the interest rate is defined for. Assumed to be infinite (`Inf`) for
    functional and constant interest rate types. Returns the `lastindex` of the vector if
    a vector type.
"""
function mt.omega(lc::LifeContingency)
    # if one of the omegas is infinity, that's a Float so we need
    # to narrow the type with Int
    return Int(omega(lc.life))
end

function mt.omega(l::SingleLife)
    return mt.omega(l.mortality) - l.issue_age + 1
end

function mt.omega(l::JointLife)
    return minimum(omega.(l.lives))
end

###################
## COMMUTATIONS ###
###################

"""
    D(lc::LifeContingency, to_time)

``D_x`` is a retrospective actuarial commutation function which is the product of the survival and discount factor.
"""
function D(lc::LifeContingency, to_time)
    return discount(lc.int, to_time) * survival(lc, to_time)
end

"""
    l(lc::LifeContingency, to_time)

``l_x`` is a retrospective actuarial commutation function which is the survival up to a certain point in time. By default, will have a unitary basis (ie `1.0`), but you can specify `basis` keyword argument to use something different (e.g. `1000` is common in the literature.)
"""
function l(lc::LifeContingency, to_time; basis = 1.0)
    return survival(lc.life, to_time) * basis
end

"""
    C(lc::LifeContingency, to_time)

``C_x`` is a retrospective actuarial commutation function which is the product of the discount factor and the difference in `l` (``l_x``).
"""
function C(lc::LifeContingency, to_time)
    discount(lc.int, to_time + 1) * (l(lc, to_time) - l(lc, to_time + 1))

end

"""
    N(lc::LifeContingency, from_time)

``N_x`` is a prospective actuarial commutation function which is the sum of the `D` (``D_x``) values from the given time to the end of the mortality table.
"""
function N(lc::LifeContingency, from_time)
    range = from_time:(omega(lc)-1)
    vals = Iterators.map(from_time -> D(lc, from_time), range)
    return sum(vals)
end

"""
    M(lc::LifeContingency, from_time)

The ``M_x`` actuarial commutation function where the `from_time` argument is `x`.
Issue age is based on the issue_age in the LifeContingency `lc`.
"""
function M(lc::LifeContingency, from_time)
    range = from_time:omega(lc)-1
    vals = Iterators.map(from_time -> C(lc, from_time), range)
    return sum(vals)
end

E(lc::LifeContingency, t, x) = D(lc, x + t) / D(lc, x)

##################
### Insurances ###
##################

abstract type Insurance end

function LifeContingency(ins::I) where {I<:Insurance}
    LifeContingency(ins.life, ins.int)
end

struct WholeLife{L,Y} <: Insurance
    life::L
    int::Y
end

struct Term{L,Y} <: Insurance
    life::L
    int::Y
    term::Int
end

"""
    Insurance(lc::LifeContingency, term)
    Insurance(life,interest, term)
    Insurance(lc::LifeContingency)
    Insurance(life,interest)

Life insurance with a term period of `term`. If `term` is `nothing`, then whole life insurance.

Issue age is based on the `issue_age` in the LifeContingency `lc`.

# Examples

```
ins = Insurance(
    SingleLife(mortality = UltimateMortality([0.5,0.5]),issue_age = 0),
    Yields.Constant(0.05),
    1           # 1 year term
)
```
"""
Insurance(lc::LifeContingency, term) = Insurance(lc.life, lc.int, term)
Insurance(lc::LifeContingency) = Insurance(lc.life, lc.int)

function Insurance(life, int, term::Int)
    term < 1 && return ZeroBenefit(life, int)
    return Term(life, int, term)
end
function Insurance(life, int)
    return WholeLife(life, int)
end

abstract type AnnuityKind end
struct Due <: AnnuityKind end
struct Immediate <: AnnuityKind end

abstract type AnnuityPayable end
abstract type AnnuityCertain <: AnnuityPayable end

struct TermCertain <: AnnuityCertain
    term::Int
    certain::Int
end
struct LifeCertain <: AnnuityCertain
    certain::Int
end
struct TermAnnuity <: AnnuityPayable
    term::Int
end
struct LifeAnnuity <: AnnuityPayable end

struct Annuity{L,Y,K<:AnnuityKind,P<:AnnuityPayable} <: Insurance
    life::L
    int::Y
    kind::K
    payable::P
    start_time::Int
    frequency::Int
end
struct ZeroBenefit{L,Y} <: Insurance
    life::L
    int::Y
end

function ZeroBenefit(lc::LifeContingency)
    return ZeroBenefit(lc.life, lc.int)
end

"""
    AnnuityDue(lc::LifeContingency; n=nothing, start_time=0; certain=nothing,frequency=1)
    AnnuityDue(life, interest; n=nothing, start_time=0; certain=nothing,frequency=1)

Annuity due with the benefit period starting at `start_time` and ending after `n` periods with `frequency` payments per year of `1/frequency` amount and a `certain` period with non-contingent payments.

# Examples

```
ins = AnnuityDue(
    SingleLife(mortality = UltimateMortality([0.5,0.5]),issue_age = 0),
    Yields.Constant(0.05),
    1, # term of policy
)
```
"""
function AnnuityDue(life, int, term; certain = nothing, start_time = 0, frequency = 1)
    term < 1 && return ZeroBenefit(life, int)
    if isnothing(certain)
        Annuity(life, int, Due(), TermAnnuity(term), start_time, frequency)
    else
        Annuity(life, int, Due(), TermCertain(term, certain), start_time, frequency)
    end
end

function AnnuityDue(life, int; certain = nothing, start_time = 0, frequency = 1)
    if isnothing(certain)
        Annuity(life, int, Due(), LifeAnnuity(), start_time, frequency)
    else
        Annuity(life, int, Due(), LifeCertain(certain), start_time, frequency)
    end
end

function AnnuityDue(lc::L, term; certain = nothing, start_time = 0, frequency = 1) where {L<:LifeContingency}
    return AnnuityDue(lc.life, lc.int, term; certain, start_time, frequency)
end

function AnnuityDue(lc::L; certain = nothing, start_time = 0, frequency = 1) where {L<:LifeContingency}
    return AnnuityDue(lc.life, lc.int; certain, start_time, frequency)
end

"""
    AnnuityImmediate(lc::LifeContingency; term=nothing, start_time=0; certain=nothing,frequency=1)
    AnnuityImmediate(life, interest; term=nothing, start_time=0; certain=nothing,frequency=1)

Annuity immediate with the benefit period starting at `start_time` and ending after `term` periods with `frequency` payments per year of `1/frequency` amount and a `certain` period with non-contingent payments.

# Examples

```
ins = AnnuityImmediate(
    SingleLife(mortality = UltimateMortality([0.5,0.5]),issue_age = 0),
    Yields.Constant(0.05),
    1 # term of policy
)
```

"""
function AnnuityImmediate(life, int, term; certain = nothing, start_time = 0, frequency = 1)
    term < 1 && return ZeroBenefit(life, int)
    if isnothing(certain)
        Annuity(life, int, Immediate(), TermAnnuity(term), start_time, frequency)
    else
        Annuity(life, int, Immediate(), TermCertain(term, certain), start_time, frequency)
    end
end

function AnnuityImmediate(life, int; certain = nothing, start_time = 0, frequency = 1)
    if isnothing(certain)
        Annuity(life, int, Immediate(), LifeAnnuity(), start_time, frequency)
    else
        Annuity(life, int, Immediate(), LifeCertain(certain), start_time, frequency)
    end
end

function AnnuityImmediate(lc::L, term; certain = nothing, start_time = 0, frequency = 1) where {L<:LifeContingency}
    return AnnuityImmediate(lc.life, lc.int, term; certain, start_time, frequency)
end

function AnnuityImmediate(lc::L; certain = nothing, start_time = 0, frequency = 1) where {L<:LifeContingency}
    return AnnuityImmediate(lc.life, lc.int; certain, start_time, frequency)
end

"""
    survival(Insurance)

The survivorship vector for the given insurance.

To get the fully computed and allocated vector, call `collect(survival(...))`.
"""
function MortalityTables.survival(ins::I) where {I<:Insurance}
    return Iterators.map(t -> survival(ins.life, t - 1), timepoints(ins))
end

function MortalityTables.survival(ins::Annuity)
    return Iterators.map(t -> survival(ins.life, t), timepoints(ins))
end

"""
    discount(Insurance)

The discount vector for the given insurance.

To get the fully computed and allocated vector, call `collect(discount(...))`.
"""
function Yields.discount(ins::I) where {I<:Insurance}
    return Iterators.map(t -> Yields.discount.(ins.int, t), timepoints(ins))
end

"""
    benefit(Insurance)

The unit benefit for the given insurance.
"""
function benefit(ins::I) where {I<:Insurance}
    return 1.0
end

function benefit(ins::ZeroBenefit)
    return 0.0
end

function benefit(ins::Annuity)
    return 1.0 / ins.frequency
end

"""
    probability(Insurance)

The vector of contingent benefit probabilities for the given insurance.

To get the fully computed and allocated vector, call `collect(probability(...))`.
"""
function probability(ins::I) where {I<:Insurance}
    return Iterators.map(timepoints(ins)) do t
        survival(ins.life, t - 1) * decrement(ins.life, t - 1, t)
    end
end

function probability(ins::ZeroBenefit)
    return Iterators.repeated(1.0, length(timepoints(ins)))
end

function probability(ins::Annuity)
    return probability(ins.payable, ins)
end

function probability(ap::AP, ins::Annuity) where {AP<:AnnuityPayable}
    return Iterators.map(t -> survival(ins.life, t), timepoints(ins))
end

function probability(ap::AP, ins::Annuity) where {AP<:AnnuityCertain}
    return Iterators.map(timepoints(ins)) do t
        t <= ap.certain + ins.start_time ? 1.0 : survival(ins.life, t)
    end
end

"""
    cashflows(Insurance)

The vector of decremented benefit cashflows for the given insurance.

To get the fully computed and allocated vector, call `collect(cashflows(...))`.
"""
function cashflows(ins::I) where {I<:Insurance}
    b = benefit(ins)
    return Iterators.map(p -> p * b, probability(ins))
end

"""
    timepoints(Insurance)

The vector of times corresponding to the cashflow vector for the given insurance.

To get the fully computed and allocated vector, call `collect(timepoints(...))`.
"""
function timepoints(ins::Insurance)::UnitRange{Int64}
    return 1:omega(ins.life)
end

function timepoints(ins::Term)::UnitRange{Int64}
    return 1:min(omega(ins.life), ins.term)
end

function timepoints(ins::ZeroBenefit)
    return Iterators.repeated(0.0, 1)
end

function timepoints(ins::Annuity)
    return timepoints(ins, ins.kind)
end

function timepoints(ins::Annuity, kind::K) where {K<:AnnuityKind}
    return timepoints(ins.payable, ins, kind)
end

function timepoints(ap::LifeCertain, ins::Annuity, ::Due)
    end_time = omega(ins.life)
    timestep = 1 / ins.frequency
    return ins.start_time:timestep:end_time
end

function timepoints(ap::LifeAnnuity, ins::Annuity, ::Due)
    # same timepoints as LifeCertain
    end_time = omega(ins.life)
    timestep = 1 / ins.frequency
    return ins.start_time:timestep:end_time
end

function timepoints(ap::TermCertain, ins::Annuity, ::Due)
    end_time = ap.term + ins.start_time - 1 / ins.frequency
    timestep = 1 / ins.frequency
    return ins.start_time:timestep:end_time
end

function timepoints(ap::TermAnnuity, ins::Annuity, ::Due)
    # same timepoints as
    end_time = ap.term + ins.start_time - 1 / ins.frequency
    timestep = 1 / ins.frequency
    return ins.start_time:timestep:end_time
end

function timepoints(ap::LifeCertain, ins::Annuity, ::Immediate)
    end_time = omega(ins.life)
    timestep = 1 / ins.frequency
    end_time = max(ins.start_time + timestep, end_time) # return at least one timepoint to avoid returning empty array
    return (ins.start_time+timestep):timestep:end_time
end

function timepoints(ap::LifeAnnuity, ins::Annuity, ::Immediate)
    # same timepoints as LifeCertain
    end_time = omega(ins.life)
    timestep = 1 / ins.frequency
    end_time = max(ins.start_time + timestep, end_time) # return at least one timepoint to avoid returning empty array
    return (ins.start_time+timestep):timestep:end_time
end

function timepoints(ap::TermCertain, ins::Annuity, ::Immediate)
    end_time = ap.term + ins.start_time
    timestep = 1 / ins.frequency
    end_time = max(ins.start_time + timestep, end_time) # return at least one timepoint to avoid returning empty array
    return (ins.start_time+timestep):timestep:end_time
end

function timepoints(ap::TermAnnuity, ins::Annuity, ::Immediate)
    # same timepoints as
    end_time = ap.term + ins.start_time
    timestep = 1 / ins.frequency
    end_time = max(ins.start_time + timestep, end_time) # return at least one timepoint to avoid returning empty array
    return (ins.start_time+timestep):timestep:end_time
end

"""
    present_value(Insurance)

The actuarial present value of the given insurance.
"""
function ActuaryUtilities.present_value(ins::T) where {T<:Insurance}
    cfs = cashflows(ins)
    times = timepoints(ins)
    yield = ins.int
    pv = present_value(yield, cfs, times)
    return pv
end

"""
    premium_net(lc::LifeContingency)
    premium_net(lc::LifeContingency,to_time)

The net premium for a whole life insurance (without second argument) or a term life insurance through `to_time`.

The net premium is based on 1 unit of insurance with the death benfit payable at the end of the year and assuming annual net premiums.
"""
function premium_net(lc::LifeContingency)
    return A(lc) / ä(lc)
end

premium_net(lc::LifeContingency, to_time) = A(lc, to_time) / ä(lc, to_time)

"""
     reserve_premium_net(lc::LifeContingency,time)

The net premium reserve at the end of year `time`.
"""
function reserve_premium_net(lc::LifeContingency, time)
    PVFB = A(lc) - A(lc, time)
    PVFP = premium_net(lc) * (ä(lc) - ä(lc, time))
    return (PVFB - PVFP) / APV(lc, time)
end

"""
    APV(lc::LifeContingency,to_time)

The **actuarial present value** which is the survival times the discount factor for the life contingency.
"""
function APV(lc::LifeContingency, to_time)
    return survival(lc, to_time) * discount(lc.int, to_time)
end

"""
    decrement(lc::LifeContingency,to_time)
    decrement(lc::LifeContingency,from_time,to_time)

Return the probablity of death for the given LifeContingency.
"""
mt.decrement(lc::LifeContingency, from_time, to_time) = 1 - survival(lc.life, from_time, to_time)

"""
    survival(lc::LifeContingency,from_time,to_time)
    survival(lc::LifeContingency,to_time)

Return the probability of survival for the given LifeContingency.
"""
mt.survival(lc::LifeContingency, to_time) = survival(lc.life, 0, to_time)
mt.survival(lc::LifeContingency, from_time, to_time) = survival(lc.life, from_time, to_time)

mt.survival(l::SingleLife, to_time) = survival(l, 0, to_time)
mt.survival(l::SingleLife, from_time, to_time) = survival(l.mortality, l.issue_age + from_time, l.issue_age + to_time, l.fractional_assump)

"""
    survival(life)

Return a survival vector for the given life.
"""
function mt.survival(l::L) where {L<:Life}
    ω = omega(l)
    return Iterators.map(t -> survival(l, t), 0:ω)
end

mt.survival(l::JointLife, to_time) = survival(l::JointLife, 0, to_time)
function mt.survival(l::JointLife, from_time, to_time)
    return survival(l.contingency, l.joint_assumption, l::JointLife, from_time, to_time)
end

function mt.survival(ins::LastSurvivor, assump::JointAssumption, l::JointLife, from_time, to_time)
    to_time == 0 && return 1.0

    l1, l2 = l.lives
    ₜpₓ = survival(l1.mortality, l1.issue_age + from_time, l1.issue_age + to_time, l1.fractional_assump)
    ₜpᵧ = survival(l2.mortality, l2.issue_age + from_time, l2.issue_age + to_time, l2.fractional_assump)
    return ₜpₓ + ₜpᵧ - ₜpₓ * ₜpᵧ
end

Yields.discount(lc::LifeContingency, t) = discount(lc.int, t)
Yields.discount(lc::LifeContingency, t1, t2) = discount(lc.int, t1, t2)

# unexported aliases
const V = reserve_premium_net
const v = Yields.discount
# A(args) = present_value(Insurance(args))
# a(args, kwargs) = present_value(AnnuityImmediate(args...; kwargs...))
# ä(args) = present_value(AnnuityDue(args))
const A = present_value ∘ Insurance
const a = present_value ∘ AnnuityImmediate
const ä = present_value ∘ AnnuityDue
const P = premium_net
const ω = omega

end # module


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*--------------------------------------------------------------------
#* Author: Gerson Balbueno Bicca
#* $Id$
#*--------------------------------------------------------------------
type DoublePipe_PressureDrop
	DoublePipe_PressureDrop()=begin
		new(
			length(Dict{Symbol,Any}(
				:Brief=>"Hydraulic Diameter of Pipe for Pressure Drop",
				:Lower=>1e-6
			)),
			press_delta(Dict{Symbol,Any}(
				:Brief=>"Total Pressure Drop",
				:Default=>0.01,
				:Lower=>0,
				:DisplayUnit=>"kPa",
				:Symbol=>"\\Delta P"
			)),
			press_delta(Dict{Symbol,Any}(
				:Brief=>"Pressure Drop for friction",
				:Default=>0.01,
				:Lower=>0,
				:DisplayUnit=>"kPa",
				:Symbol=>"\\Delta P_{fric}"
			)),
			press_delta(Dict{Symbol,Any}(
				:Brief=>"Pressure Drop due to return",
				:Default=>0.01,
				:Lower=>0,
				:DisplayUnit=>"kPa",
				:Symbol=>"\\Delta P_{return}"
			)),
			fricfactor(Dict{Symbol,Any}(
				:Brief=>"Friction Factor",
				:Default=>0.05,
				:Lower=>1e-10,
				:Upper=>2000
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Reynolds Number",
				:Default=>100,
				:Lower=>1
			)),
			[:Dh,],
			[:Pdrop,:Pd_fric,:Pd_ret,:fi,:Re,]
		)
	end
	Dh::length
	Pdrop::press_delta
	Pd_fric::press_delta
	Pd_ret::press_delta
	fi::fricfactor
	Re::positive
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export DoublePipe_PressureDrop
function atributes(in::DoublePipe_PressureDrop,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=false
	fields[:Brief]="to be documented"
	fields[:Info]="to be documented"
	drive!(fields,_)
	return fields
end
DoublePipe_PressureDrop(_::Dict{Symbol,Any})=begin
	newModel=DoublePipe_PressureDrop()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
# Autogenerated wrapper script for SDPA_jll for x86_64-w64-mingw32-libgfortran4-cxx11
export libsdpa, libsdpawrap, sdpa

using libcxxwrap_julia_jll
using OpenBLAS32_jll
using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("SDPA")
JLLWrappers.@declare_library_product(libsdpa, "libsdpa-0.dll")
JLLWrappers.@declare_library_product(libsdpawrap, "libsdpawrap.dll")
JLLWrappers.@declare_executable_product(sdpa)
function __init__()
    JLLWrappers.@generate_init_header(libcxxwrap_julia_jll, OpenBLAS32_jll, CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libsdpa,
        "bin\\libsdpa-0.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsdpawrap,
        "bin\\libsdpawrap.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        sdpa,
        "bin\\sdpa.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
rootdir = Pkg.dir("Dyno")
filename = joinpath(rootdir,"examples","rbc.mod")

import Dolang

import Dyno
model = Dyno.import_modfile(filename)

sol = Dyno.solve(model)
Dyno.simulate(sol)


text:
@testset "832.flipping-an-image.jl" begin
    @test flip_and_invert_image([[1, 1, 0], [1, 0, 1], [0, 0, 0]]) ==
          [[1, 0, 0], [0, 1, 0], [1, 1, 1]]
    @test flip_and_invert_image([[1, 1, 0, 0], [1, 0, 0, 1], [0, 1, 1, 1], [1, 0, 1, 0]]) ==
          [[1, 1, 0, 0], [0, 1, 1, 0], [0, 0, 0, 1], [1, 0, 1, 0]]
    @test flip_and_invert_image([[1]]) == [[0]]
    @test flip_and_invert_image([[1, 0, 0], [1, 1, 1], [0, 0, 0]]) ==
          [[1, 1, 0], [0, 0, 0], [1, 1, 1]]
end


text:
#*************************************************************************#
#
# Initializer Abstract Type
#
#*************************************************************************#
abstract type Initializer end

"""
    initial_guess(f::StaticModel[, init::Initializer]) -> params

For a given static model `f` provide a suitable initial guess
for the parameters vectors `params`. This could be used
as start point for an optimization procedure. The provided
`Initializer` will determine what strategy will be used.
"""
function initial_guess(f::ParametricModel)
    initial_guess(f, DefaultInitializer())
end

"""
    DefaultInitializer()

Construct `Initializer` for static model using the default settings
for the given class.

Valid for any ParametricModel.
"""
struct DefaultInitializer <: Initializer
end

function initial_guess(f::ParametricModel, init::DefaultInitializer)
    initial_guess(f, ZeroInitializer())
end

"""
    ZeroInitializer()

Construct `Initializer` for static model. The `ZeroInitializer`
will initialize any parameter vector with zeros.

Valid for any ParametricModel.
"""
struct ZeroInitializer <: Initializer
end

function initial_guess(f::ParametricModel, init::ZeroInitializer)
    return zeros(f.nparams)
end

"""
    CustomRandomInitializer([s::Sampleable])

Construct `Initializer` for static model. The `CustomRandomInitializer`
will initialize any parameter vector with of random values picked
from the given distribution.

Valid for any ParametricModel.
"""
struct CustomRandomInitializer{F<:Distributions.VariateForm,
                               S<:Distributions.ValueSupport} <: Initializer
    s::Distributions.Sampleable{F, S}
end

function initial_guess(f::ParametricModel, init::CustomRandomInitializer)
    return rand(init.s, f.nparams)
end


text:
using SparseArrays, LinearAlgebra, Statistics

function kkt_col_norms!(P::AbstractMatrix{T}, A::AbstractMatrix{T}, norm_LHS::AbstractVector{T}, norm_RHS::AbstractVector{T}) where {T <: AbstractFloat}
	col_norms!(norm_LHS, P, reset = true);   #start from zero
	col_norms!(norm_LHS, A, reset = false);  #incrementally from P norms
	row_norms!(norm_RHS, A)                 #same as column norms of A'
	return nothing
end

function limit_scaling!(s::Vector{T}, set::COSMO.Settings{T}) where {T <: AbstractFloat}
	@.s = clip(s, set.MIN_SCALING, set.MAX_SCALING, one(T))
	return nothing
end

function limit_scaling(s::T, set::COSMO.Settings{T}) where {T <: AbstractFloat}
	s = clip(s, set.MIN_SCALING, set.MAX_SCALING, one(T))
	return s
end

function scale_ruiz!(ws::COSMO.Workspace{T}) where {T <: AbstractFloat}

	set = ws.settings
	#references to scaling matrices from workspace
	D    = ws.sm.D
	E    = ws.sm.E
	c    = ws.sm.c[]

	#unit scaling to start
	D.diag .= one(T)
	E.diag .= one(T)
	c       = one(T)

	#use the inverse scalings as intermediate
	#work vectors as well, since we don't
	#compute the inverse scaling until the
	#final step
	Dwork = ws.sm.Dinv
	Ework = ws.sm.Einv

	#references to QP data matrices
	P = ws.p.P
	A = ws.p.A
	q = ws.p.q
	b = ws.p.b

	#perform scaling operations for a fixed
	#number of steps, i.e. no tolerance or
	#convergence check
	for i = 1:set.scaling

		kkt_col_norms!(P, A, Dwork.diag, Ework.diag)
		limit_scaling!(Dwork.diag, set)
		limit_scaling!(Ework.diag, set)

		inv_sqrt!(Dwork.diag)
		inv_sqrt!(Ework.diag)

		# Scale the problem data and update the
		# equilibration matrices
		scale_data!(P, A, q, b, Dwork, Ework , one(T))
		LinearAlgebra.lmul!(Dwork, D)        #D[:,:] = Dtemp*D
		LinearAlgebra.lmul!(Ework, E)        #D[:,:] = Dtemp*D

		# now use the Dwork array to hold the
		# column norms of the newly scaled P
		# so that we can compute the mean
		col_norms!(Dwork.diag, P)
		mean_col_norm_P = mean(Dwork.diag)
		inf_norm_q      = norm(q, Inf)

		if mean_col_norm_P  != 0. && inf_norm_q != 0.

			inf_norm_q = limit_scaling(inf_norm_q, set)
			scale_cost = max(inf_norm_q, mean_col_norm_P)
			scale_cost = limit_scaling(scale_cost, set)
			ctmp = one(T) / scale_cost

			# scale the penalty terms and overall scaling
			scalarmul!(P, ctmp)
			q       .*= ctmp
			c        *= ctmp
		end

	end #end Ruiz scaling loop

	# for certain cones we can only use a
	# a single scalar value.  In these cases
	# compute an adjustment to the overall scaling
	# so that the aggregate scaling on the cone
	# in questions turns out to be component-wise eq
	if rectify_set_scalings!(E, Ework, ws.p.C)
		#only rescale if the above returns true,
		#i.e. some cone scalings were rectified
		scale_data!(P, A, q, b, I, Ework, one(T))
		LinearAlgebra.lmul!(Ework, E)
	end

	issymmetric(ws.p.P) || symmetrize_full!(ws.p.P)
	#scale set components
	scale_sets!(E, ws.p.C)

	#update the inverse scaling data, c and c_inv
	ws.sm.Dinv.diag .= one(T) ./ D.diag
	ws.sm.Einv.diag .= one(T) ./ E.diag

	#These are Base.RefValue type so that
	#scaling can remain an immutable
	ws.sm.c[]        = c
	ws.sm.cinv[]     = one(T) ./ c

	# scale the potentially warm started variables
	ws.vars.x[:] = ws.sm.Dinv * ws.vars.x
	ws.vars.μ[:] = ws.sm.Einv * ws.vars.μ
	ws.vars.s.data[:] = ws.sm.E * ws.vars.s.data

	ws.vars.μ  .*= c

	return nothing
end

function inv_sqrt!(A::Vector{T}) where{T <: Real}
	@fastmath A .= one(T) ./ sqrt.(A)
end

function rectify_set_scalings!(E::AbstractMatrix{T}, Ework::AbstractMatrix{T}, C::AbstractConvexSet{T}) where {T <: AbstractFloat}

	Ework.diag         .= one(T)

	# FIX ME: split the vector Ework.diag over the sets.
	# This should NOT be done here.  Instead, the
	# splitting should happen at the time the scaling
	# structure is allocated, so that E is a Diagonal
	# of a SplitVector type.
	Es  = SplitVector(E.diag, C)
	Ews = SplitVector(Ework.diag, C)
	scale_changed = rectify_scaling!(Es, Ews, C)
	return scale_changed
end

function scale_sets!(E::AbstractMatrix{T}, C::AbstractConvexSet{T}) where {T <: AbstractFloat}

	# FIX ME: split the vector Ework.diag over the sets.
	# This should NOT be done here.  Instead, the
	# splitting should happen at the time the scaling
	# structure is allocated, so that E is a Diagonal
	# of a SplitVector type.
	Es = SplitVector(E.diag, C)
	scale!(C, Es)
end

function scale_data!(P::AbstractMatrix{T}, A::AbstractMatrix{T}, q::AbstractVector{T}, b::AbstractVector{T}, Ds::Union{IdentityMatrix, AbstractMatrix{T}}, Es::Union{IdentityMatrix, AbstractMatrix{T}}, cs::T = one(T)) where {T <: AbstractFloat}

	lrmul!(Ds, P, Ds) # P[:,:] = Ds*P*Ds
	lrmul!(Es, A, Ds) # A[:,:] = Es*A*Ds
	q[:] = Ds * q
	b[:] = Es * b
	if cs != one(T)
		scalarmul!(P,cs)
		q .*= cs
	end
	return nothing
end

function reverse_scaling!(ws::COSMO.Workspace)

	cinv = ws.sm.cinv[] #from the Base.RefValue type
	ws.vars.x[:] = ws.sm.D * ws.vars.x
	ws.vars.s[:] = ws.sm.Einv * ws.vars.s
	ws.vars.μ[:] = ws.sm.E * ws.vars.μ
	ws.vars.μ  .*= cinv

	# reverse scaling for model data
	if ws.flags.REVERSE_SCALE_PROBLEM_DATA
		scale_data!(ws.p.P, ws.p.A, ws.p.q, ws.p.b,
			ws.sm.Dinv, ws.sm.Einv, cinv)
	end
	return nothing
end


text:
@testset "NLCOP" begin
    if !is_travis
        @testset "Exact" begin
            @testset "Basic: uniform matroid" begin
                m = 2
                lw = Float64[1, 2, 3, 4, 5]
                nlw = Float64[5, 4, 3, 2, 1]
                # Maximum value for the sqrt(sum of the two chosen nlw): sqrt(5 + 4) = 3.0.
                # Minimum value for the sqrt(sum of the two chosen nlw): sqrt(2 + 1) ≈ 1.73.

                li = UniformMatroidInstance(lw, m, Maximise())
                nli = NonlinearCombinatorialInstance(
                    li,
                    lw,
                    nlw,
                    SquareRoot,
                    1.0,
                    DefaultLinearFormulation(Gurobi.Optimizer),
                    true,
                    DefaultLinearFormulation(Gurobi.Optimizer),
                )
                s = solve(nli, ExactNonlinearSolver(Gurobi.Optimizer))

                @test s.instance == li
                @test Set(s.variables) == Set([4, 5])
            end

            @testset "Basic: spanning tree" begin
                graph = complete_graph(5)
                lw = Dict{Edge{Int}, Float64}(e => src(e) for e in edges(graph))
                nlw = Dict{Edge{Int}, Float64}(
                    e => ne(graph) - src(e) for e in edges(graph)
                )

                li = SpanningTreeInstance(graph, lw, Maximise())
                nli = NonlinearCombinatorialInstance(
                    li,
                    lw,
                    nlw,
                    SquareRoot,
                    1.0,
                    DefaultLinearFormulation(Gurobi.Optimizer),
                    true,
                    DefaultLinearFormulation(Gurobi.Optimizer),
                )
                s = solve(nli, ExactNonlinearSolver(Gurobi.Optimizer))

                obj_lin = sum(lw[e] for e in s.variables)
                obj_nl = sqrt(sum(nlw[e] for e in s.variables))

                # There is a unique way to represent the solution, as the LP
                # formulation ensures that 1 is the root of the tree. Thus,
                # check explicitly the solution.
                @test s.instance == li
                @test Set(s.variables) ==
                      Set([Edge(1, 2), Edge(2, 3), Edge(3, 4), Edge(4, 5)])
                @test obj_lin + obj_nl ≈ 15.47 atol = 1.0e-2
            end

            @testset "Basic: elementary path" begin
                graph = complete_digraph(5)
                lw = Dict{Edge{Int}, Float64}(
                    e => min(src(e), dst(e)) for e in edges(graph)
                )
                nlw = Dict{Edge{Int}, Float64}(
                    e => ne(graph) - min(src(e), dst(e)) for e in edges(graph)
                )

                li = ElementaryPathInstance(graph, lw, 1, 5, Maximise())
                nli = NonlinearCombinatorialInstance(
                    li,
                    lw,
                    nlw,
                    SquareRoot,
                    1.0,
                    DefaultLinearFormulation(Gurobi.Optimizer),
                    true,
                    DefaultLinearFormulation(Gurobi.Optimizer),
                )
                s = solve(nli, ExactNonlinearSolver(Gurobi.Optimizer))

                obj_lin = sum(lw[e] for e in s.variables)
                obj_nl = sqrt(sum(nlw[e] for e in s.variables))

                @test s.instance == li
                @test Set(s.variables) ==
                      Set([Edge(1, 2), Edge(2, 3), Edge(3, 4), Edge(4, 5)])
                @test obj_lin + obj_nl ≈ 18.36 atol = 1.0e-2
            end
        end
    end

    @testset "Approximate" begin
        @testset "Basic: uniform matroid, DP-based" begin
            m = 2
            lw = Float64[1, 2, 3, 4, 5]
            nlw = Float64[5, 4, 3, 2, 1]
            # Maximum value for the sqrt(sum of the two chosen nlw): sqrt(5 + 4) = 3.0.
            # Minimum value for the sqrt(sum of the two chosen nlw): sqrt(2 + 1) ≈ 1.73.

            li = UniformMatroidInstance(lw, m, Maximise())
            nli = NonlinearCombinatorialInstance(
                li,
                lw,
                nlw,
                SquareRoot,
                0.01,
                DynamicProgramming(),
                true,
                DefaultLinearFormulation(nothing),
            )
            s = solve(nli, ApproximateNonlinearSolver(DynamicProgramming()))

            @test s.instance == li
            @test Set(s.variables) == Set([4, 5])
        end

        if !is_travis
            @testset "Basic: uniform matroid, LP-based" begin
                m = 2
                lw = Float64[1, 2, 3, 4, 5]
                nlw = Float64[5, 4, 3, 2, 1]
                # Maximum value for the sqrt(sum of the two chosen nlw): sqrt(5 + 4) = 3.0.
                # Minimum value for the sqrt(sum of the two chosen nlw): sqrt(2 + 1) ≈ 1.73.

                li = UniformMatroidInstance(lw, m, Maximise())
                nli = NonlinearCombinatorialInstance(
                    li,
                    lw,
                    nlw,
                    SquareRoot,
                    0.01,
                    DefaultLinearFormulation(Cbc.Optimizer),
                    true,
                    DefaultLinearFormulation(Cbc.Optimizer),
                )
                s = solve(nli, ApproximateNonlinearSolver(DefaultLinearFormulation(Cbc.Optimizer)))

                @test s.instance == li
                @test Set(s.variables) == Set([4, 5])
            end
        end

        @testset "Basic: spanning tree, DP-based" begin
            graph = complete_graph(5)
            lw = Dict{Edge{Int}, Float64}(
                e => min(src(e), dst(e)) for e in edges(graph)
            )
            nlw = Dict{Edge{Int}, Float64}(
                e => ne(graph) - min(src(e), dst(e)) for e in edges(graph)
            )

            li = SpanningTreeInstance(graph, lw, Maximise())
            @test fastest_exact(li) !== nothing
            nli = NonlinearCombinatorialInstance(
                li,
                lw,
                nlw,
                SquareRoot,
                1.0,
                DynamicProgramming(),
                true,
                DefaultLinearFormulation(nothing),
            )
            @test fastest_exact(nli) === nothing
            s = solve(nli, ApproximateNonlinearSolver(DynamicProgramming()))

            obj_lin = sum(lw[e] for e in s.variables)
            obj_nl = sqrt(sum(nlw[e] for e in s.variables))

            @test s.instance == li
            @test obj_lin + obj_nl ≈ 15.47 atol = 1.0e-2
        end

        @testset "Basic: spanning tree, LP-based" begin
            graph = complete_graph(5)
            lw = Dict{Edge{Int}, Float64}(
                e => min(src(e), dst(e)) for e in edges(graph)
            )
            nlw = Dict{Edge{Int}, Float64}(
                e => ne(graph) - min(src(e), dst(e)) for e in edges(graph)
            )

            li = SpanningTreeInstance(graph, lw, Maximise())
            nli = NonlinearCombinatorialInstance(
                li,
                lw,
                nlw,
                SquareRoot,
                1.0,
                DefaultLinearFormulation(Cbc.Optimizer),
                true,
                DefaultLinearFormulation(Cbc.Optimizer),
            )
            s = solve(nli, ApproximateNonlinearSolver(DefaultLinearFormulation(Cbc.Optimizer)))

            obj_lin = sum(e in keys(lw) ? lw[e] : lw[reverse(e)] for e in s.variables)
            obj_nl = sqrt(sum(e in keys(nlw) ? nlw[e] : lw[reverse(e)] for e in s.variables))

            @test s.instance == li
            @test obj_lin + obj_nl ≈ 15.47 atol = 1.0e-2
        end

        @testset "Basic: elementary path, DP-based" begin
            graph = complete_digraph(5)
            lw = Dict{Edge{Int}, Float64}(e => src(e) for e in edges(graph))
            nlw = Dict{Edge{Int}, Float64}(
                e => ne(graph) - src(e) for e in edges(graph)
            )

            li = ElementaryPathInstance(graph, lw, 1, 5, Maximise())
            nli = NonlinearCombinatorialInstance(
                li,
                lw,
                nlw,
                SquareRoot,
                1.0,
                DynamicProgramming(),
                true,
                DefaultLinearFormulation(nothing),
            )
            s = solve(nli, ApproximateNonlinearSolver(DynamicProgramming()))

            obj_lin = sum(lw[e] for e in s.variables)
            obj_nl = sqrt(sum(nlw[e] for e in s.variables))

            @test s.instance == li
            @test Set(s.variables) ==
                  Set([Edge(1, 2), Edge(2, 3), Edge(3, 4), Edge(4, 5)])
            @test obj_lin + obj_nl ≈ 18.36 atol = 1.0e-2
        end

        if !is_travis
            @testset "Basic: elementary path, LP-based" begin
                graph = complete_digraph(5)
                lw = Dict{Edge{Int}, Float64}(e => src(e) for e in edges(graph))
                nlw = Dict{Edge{Int}, Float64}(
                    e => ne(graph) - src(e) for e in edges(graph)
                )

                li = ElementaryPathInstance(graph, lw, 1, 5, Maximise())
                nli = NonlinearCombinatorialInstance(
                    li,
                    lw,
                    nlw,
                    SquareRoot,
                    1.0,
                    DefaultLinearFormulation(Gurobi.Optimizer),
                    true,
                    DefaultLinearFormulation(Gurobi.Optimizer),
                )
                s = solve(nli, ApproximateNonlinearSolver(DefaultLinearFormulation(Gurobi.Optimizer)))

                obj_lin = sum(lw[e] for e in s.variables)
                obj_nl = sqrt(sum(nlw[e] for e in s.variables))

                @test s.instance == li
            end
        end
    end
end


text:
module DatesInGerman

using Dates

const MONTHS = ("januar", "februar", "märz", "april", "mai", "juni", "juli", "august", "september", "oktober", "november", "dezember")

function parsefrom(date::String; inwords::Bool=true)::Date
    inwords ? date |> fromwords : date |> fromnumbers
end

function fromwords(date::String)::Date
    parts = match(r"(\d+).\s*(\S+)\s*(\d{4})", date).captures
    Date(parse(Int, parts[3]),
           findfirst(m -> m == lowercase(parts[2]), MONTHS),
           parse(Int, parts[1]))
end

function fromnumbers(date::String)::Date
    parts = match(r"(\d+)\.(\d+)\.(\d+)", date).captures
    Date(parse(Int, (parts[3])),
             parse(Int, parts[2]),
             parse(Int, parts[1]))
end

end # module


text:
PowerLawPlateau =  RheoModelClass(
        # Model name
        name="burgerliquid",
        # Model parameters,
        p = (:Gᵩ, :G₀, :τ, :α),
        # Relaxation modulus
        G = quote
                Gᵩ + (G₀ - Gᵩ)./(1 + t/τ).^(α)
            end,
        # Network
        info= "                "
        )


text:
using Test
using DomainSets, HCubature, StaticArrays

using DomainIntegrals

include("test_measures.jl")
include("test_integrals.jl")
include("test_gauss.jl")

@testset "Measures" begin
    test_measures()
end

@testset "Quadrature rules" begin
    test_gauss()
end

@testset "Selected integrals" begin
    test_integrals()
end


text:
using Documenter, DiseaseOutbreak

makedocs(;
    modules=[DiseaseOutbreak],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/chelseatrotter/DiseaseOutbreak.jl/blob/{commit}{path}#L{line}",
    sitename="DiseaseOutbreak.jl",
    authors="Sen Research Group, Saunak Sen, Hyeonju Kim, Gregory Farage, Chelsea Trotter",
    assets=String[],
)

deploydocs(;
    repo="github.com/chelseatrotter/DiseaseOutbreak.jl",
)


text:
@testset "Curve struct" begin
@testset "Uniformly distributed values" begin
    n = 20
    data = rand(Uniform(0,10), n)
    curve = CalculatedABC.ABCcurve(data)
    @test curve.effort == 0:0.05:1.0
    expected_yield = [0,
                      0.10223370412290836,
                      0.19630374973228518,
                      0.2798337861487627,
                      0.35309289814499817,
                      0.4251124119444249,
                      0.49223359866263033,
                      0.5559340313195282,
                      0.6192318028485928,
                      0.6706485917196718,
                      0.7210686469145637,
                      0.7696104643520768,
                      0.812596101423477,
                      0.8466456273448749,
                      0.8777009051387628,
                      0.9076823618805581,
                      0.9368834040175215,
                      0.9647084964342312,
                      0.9824989229044542,
                      0.994333140487581,
                      1.0000000000000002]
    @test isapprox(curve.yield, expected_yield)
end

@testset "Exponentially distributed values" begin
    n = 20
    data = rand(Exponential(1), n)
    curve = CalculatedABC.ABCcurve(data)
    @test curve.effort == 0:0.05:1.0
    expected_yield = [0.0,
                      0.1559704817922999,
                      0.259107787045309,
                      0.36017085915674457,
                      0.4567531613173713,
                      0.530483819367366,
                      0.5905820392527759,
                      0.6503414525690904,
                      0.6912925133668514,
                      0.7318440937106174,
                      0.7704643555492486,
                      0.8031964684921609,
                      0.8353247928350085,
                      0.8670259991394021,
                      0.8976554662905372,
                      0.9212877248856802,
                      0.9431001183934831,
                      0.964704167187443,
                      0.9810003669812493,
                      0.9911110632975999,
                      1.0]
    @test isapprox(curve.yield, expected_yield)
end

@testset "SwissInhabitants data" begin
    data = readdlm("swissinhabitants")[:,1]
    curve = CalculatedABC.ABCcurve(data)
    @test curve.effort == 0:0.0003453038674033149:1.0
    expected_yield = readdlm("expected_yield")[:,1]
    @test isapprox(curve.yield, expected_yield)
end
end


text:
using .StanSample

include("Stan/model_summary.jl")


text:
function adjacency_list(g::AbstractSimpleWeightedGraph)
    N = nv(g)
    Vector{Int}[outneighbors(g, i) for i = 1:N]
end


text:
### A Pluto.jl notebook ###
# v0.12.10

using Markdown
using InteractiveUtils

# ╔═╡ 11cdad30-10e1-11eb-0761-4181b97481a8
using Statistics

# ╔═╡ c35e0efe-100e-11eb-0f2f-178bcf665b41
begin
	import Pkg
	Pkg.activate(mktempdir())
end

# ╔═╡ ce8d10ae-100e-11eb-24b1-fb5b8786d728
begin
	Pkg.add(["Images", "ImageIO", "ImageMagick", "TestImages"])
	using Images
end

# ╔═╡ 26a3bca0-1013-11eb-0127-d30ccb95a826
md"""Saving a URL into Variable """

# ╔═╡ 111fb070-1011-11eb-0d49-2992c403e821
usflag = "https://foreres.files.wordpress.com/2010/04/small-american-flag.jpg"

# ╔═╡ 488c2820-1013-11eb-3aeb-597f0f1bf852
md"""Saving a URL variable to defined Location"""

# ╔═╡ 10e51870-1011-11eb-1781-3984a2cf5e94
US_Flag = download(usflag, "D:\\MyGit\\Colors\\Images\\US_Flag.jpg")

# ╔═╡ 8bba0860-1013-11eb-11f5-6dbba1917f11
md"""Opening the Image from the Saved Location"""

# ╔═╡ d16cb100-100e-11eb-1fb8-bf916eb323f1
US_Flag01 = load(US_Flag)

# ╔═╡ 9f3a4030-1013-11eb-2644-ef11c2ca7ab5
md"""Opening an image directly from Local Drive"""

# ╔═╡ 09439860-1013-11eb-0f92-db4886864291
ank = load("D:\\MyGit\\Colors\\Images\\US_Flag.jpg")

# ╔═╡ d12f0bc0-100e-11eb-2a60-17ab6eedc7d8
size(US_Flag01)

# ╔═╡ d0f6bdb0-100e-11eb-0c1d-b9082d7c2489
typeof(US_Flag01)

# ╔═╡ dbd14f10-1014-11eb-0cfd-7389af14ac04
md"""Created this Function Below to bring down the size of a pic, dont know how it works"""

# ╔═╡ 9127ff3e-1014-11eb-1078-43ad4981916f
decimate(image, ratio=5) = image[1:ratio:end, 1:ratio:end]

# ╔═╡ d0881d60-100e-11eb-0886-1b0e3ac7b2ef
usflsmall = let
	original = Images.load("D:\\MyGit\\Colors\\Images\\US_Flag.jpg")
	decimate(original, 3)
end

# ╔═╡ 1435cd42-1015-11eb-2fee-4122e629b624
md"""US_Flag01 (150,150 ) size brought down to 50,50"""

# ╔═╡ c8c39a10-1012-11eb-321c-c5bd5e28f744
size(usflsmall)

# ╔═╡ fd9e9170-1014-11eb-29b4-77028101ab8f
usflsmall[20:30,10:20,1]

# ╔═╡ fd6977b0-1014-11eb-3e6c-adf2e39e42b4
usflsmall[25,25]

# ╔═╡ fd360ba0-1014-11eb-3a59-b7599a434aec
md"""
#### Exercise 2.1
👉 Write a function **`mean_colors`** that accepts an object called `image`. It should calculate the mean (average) amounts of red, green and blue in the image and return a tuple `(r, g, b)` of those means.
"""

# ╔═╡ 12176060-10e1-11eb-2035-936ab2404587
#[ mean(float.(c.(img))) for c = [red,green,blue], img = [usflsmall] ]

# ╔═╡ 11e44270-10e1-11eb-02d3-d10adc187074
#mean(float.(red.(usflsmall)))

# ╔═╡ 11b60680-10e1-11eb-33bc-6f416cadb5dd
function mean_colors(image)
	missing = [ mean(float.(c.(img))) for c = [red,green,blue], img = [image] ]
	return 	missing
end

# ╔═╡ 11a0f7e0-10e1-11eb-0306-914b8bedeb72
mean_colors(usflsmall)

# ╔═╡ 9dcaa230-10e5-11eb-2b81-31318eda5ab0
typeof(usflsmall)

# ╔═╡ 118a62a0-10e1-11eb-362f-6f7ea7877251
size(usflsmall)

# ╔═╡ 1173cd60-10e1-11eb-0bed-196ed41f75fb
US_Flag01

# ╔═╡ 115b3c50-10e1-11eb-33d4-8b33fcb585bf
display(usflsmall

# ╔═╡ 1142ab40-10e1-11eb-2eec-39cfe8b06c73

# ╔═╡ 11353dc2-10e1-11eb-1ab8-67c8d170a5fe

# ╔═╡ 110a0f0e-10e1-11eb-319f-af3611722baf

# ╔═╡ c81fb8ee-1012-11eb-00a8-75af624da43b

# ╔═╡ Cell order:
# ╠═c35e0efe-100e-11eb-0f2f-178bcf665b41
# ╠═ce8d10ae-100e-11eb-24b1-fb5b8786d728
# ╠═26a3bca0-1013-11eb-0127-d30ccb95a826
# ╠═111fb070-1011-11eb-0d49-2992c403e821
# ╟─488c2820-1013-11eb-3aeb-597f0f1bf852
# ╠═10e51870-1011-11eb-1781-3984a2cf5e94
# ╟─8bba0860-1013-11eb-11f5-6dbba1917f11
# ╠═d16cb100-100e-11eb-1fb8-bf916eb323f1
# ╟─9f3a4030-1013-11eb-2644-ef11c2ca7ab5
# ╠═09439860-1013-11eb-0f92-db4886864291
# ╠═d12f0bc0-100e-11eb-2a60-17ab6eedc7d8
# ╠═d0f6bdb0-100e-11eb-0c1d-b9082d7c2489
# ╟─dbd14f10-1014-11eb-0cfd-7389af14ac04
# ╟─9127ff3e-1014-11eb-1078-43ad4981916f
# ╠═d0881d60-100e-11eb-0886-1b0e3ac7b2ef
# ╟─1435cd42-1015-11eb-2fee-4122e629b624
# ╠═c8c39a10-1012-11eb-321c-c5bd5e28f744
# ╠═fd9e9170-1014-11eb-29b4-77028101ab8f
# ╠═fd6977b0-1014-11eb-3e6c-adf2e39e42b4
# ╟─fd360ba0-1014-11eb-3a59-b7599a434aec
# ╠═12176060-10e1-11eb-2035-936ab2404587
# ╠═11e44270-10e1-11eb-02d3-d10adc187074
# ╠═11cdad30-10e1-11eb-0761-4181b97481a8
# ╠═11b60680-10e1-11eb-33bc-6f416cadb5dd
# ╠═11a0f7e0-10e1-11eb-0306-914b8bedeb72
# ╠═9dcaa230-10e5-11eb-2b81-31318eda5ab0
# ╠═118a62a0-10e1-11eb-362f-6f7ea7877251
# ╠═1173cd60-10e1-11eb-0bed-196ed41f75fb
# ╠═115b3c50-10e1-11eb-33d4-8b33fcb585bf
# ╠═1142ab40-10e1-11eb-2eec-39cfe8b06c73
# ╠═11353dc2-10e1-11eb-1ab8-67c8d170a5fe
# ╠═110a0f0e-10e1-11eb-319f-af3611722baf
# ╠═c81fb8ee-1012-11eb-00a8-75af624da43b


text:
######################################################################
# pairings.jl: Pairings on elliptic curves
######################################################################

export

######################################################################
# Miller's algorithm
######################################################################

function tangent(P::EllipticPoint)
end

function line(P::EllipticPoint, Q::EllipticPoint)
end

function miller(P::EllipticPoint, l::Integer)
	isneutral(l * P) || throw(ArgumentError("Given point in 'miller' is not a torsion point"))
	K = base_ring(P.curve)
	A, (x, y) = PolynomialRing(K, ["X", "Y"])
	E = P.curve
	O = neutral(E)

	f = 1
	T = P

	for b in bin(l)[2:end]
		D = tangent(P)
		V = line(2 * T, 0)
		f = f^2 * D // V
		if b == '1'
			...
		end
	end
end

######################################################################
# The Weil pairing
######################################################################

######################################################################
# The Tate pairing
######################################################################


text:
### A Pluto.jl notebook ###
# v0.17.2

# using Markdown
# using InteractiveUtils

# ╔═╡ ab7fc0f6-d065-4c03-aeb5-64a49214399e
const default_exprlist = [Symbol("@md_str"), Symbol("@benchmark"), :(using Markdown), :PLUTO_MANIFEST_TOML_CONTENTS, :PLUTO_PROJECT_TOML_CONTENTS]

# ╔═╡ 2ac5f3c6-96f8-4f24-8022-1f0470672383
#=╠═╡ notebook_exclusive
"""
    hasexpr(expr, exprlist)

Simple expression match; will return `true` if one of the exprs listed in `exprlist` can
be found inside `expr`. Checking for macrocalls can be done by including the `Symbol` of the macro among the expression list

```julia
hasexpr(:(@benchmark 3+2), Symbol("@benchmark")) == true
```
"""
function hasexpr(ex, exprlist::AbstractVector)
  result = false
  MacroTools.postwalk(ex) do y
    if y in exprlist
      result = true
    end
    return y
  end
  return result
end
  ╠═╡ notebook_exclusive =#

# ╔═╡ 1663bd40-4a0c-42cd-b716-4583d0a1bb66
# Fallback for when a single element is provided
hasexpr(ex,exprlist) = hasexpr(ex,[exprlist])

# ╔═╡ b42062a6-5bca-4c47-835c-caaeee1ceac6
"""
	include_mapexpr([exprlist])::Function
Returns a function that can be used as first argument (`mapexpr`) of [`include`](@ref) calls
to avoid including some of the expressions inside a Pluto notebook.

When applied to an expression `ex`, the output function traverse `ex` recursively and
returns `nothing` if `ex ∈ exprlist` and `ex` otherwise.

The expression walk is performed with [`MacroTools.postwalk`](@ref) so for example the name
of macros can be used to filter out nested expr which contains the macro call.

When called without arguments, the value of `exprlist` defaults to a vector containing some
filtering entries for Pluto notebooks:

	default_exprlist = $default_exprlist

```julia
include_mapexpr()(:(md"asd")) == nothing
include_mapexpr()(:(3 + 2)) == :(3+2)
include_mapexpr([3])(:(3 + 2)) == nothing
```
"""
include_mapexpr(exprlist=default_exprlist) = ex -> hasexpr(ex,exprlist) ? nothing : ex

# ╔═╡ b0f0df45-ed36-454a-9a9d-215f8c23836a
export include_mapexpr

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.0-rc2"
manifest_format = "2.0"

[deps]
"""

# ╔═╡ Cell order:
# ╠═ab7fc0f6-d065-4c03-aeb5-64a49214399e
# ╠═b0f0df45-ed36-454a-9a9d-215f8c23836a
# ╠═2ac5f3c6-96f8-4f24-8022-1f0470672383
# ╠═1663bd40-4a0c-42cd-b716-4583d0a1bb66
# ╠═b42062a6-5bca-4c47-835c-caaeee1ceac6
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
using Flux.Tracker: TrackedArray, TrackedReal
"""
	pairwisel2(x,y)

	pairwise distances between `x` and `y` using Euclidean Distance.

"""
pairwisel2(x,y) = -2 .* x' * y .+ sum(x.^2, dims=1)' .+ sum(y.^2, dims=1)

"""
	scaled_pairwisel2(x::M, y::M, σ::M) where {M<:Matrix}
	when x,y, and σ are matrices, than ``o_{i,j} = \\frac{x_i - y_j}{σ_i}`` where `σ` are treated as a column vectors.

	if σ is:
		- Number means a σ shared by all `x`
		- Vector means a σ is diagonal shared by all `x`
		- Transpose means each column of `x` has its own scalar σ
		- Matrix means each column of `x` has its own scalar diagonal σ

"""
scaled_pairwisel2(x, y, σ::T) where {T<: Union{Real,TrackedReal}} = pairwisel2(x, y) ./ σ^2
scaled_pairwisel2(x, y, σ::T) where {T<:Union{AbstractVector, TrackedArray{T, N, A} where {T, N, A<: AbstractVector}}} = pairwisel2(x ./σ, y./σ)
# scaled_pairwisel2(x, y, σ::T) where {T<:Union{Transpose, TrackedArray{T, N, A} where {T, N, A<: Transpose}}} = pairwisel2(x, y) ./ (σ'.^2)
scaled_pairwisel2(x::AbstractMatrix, y::AbstractMatrix, σ::AbstractMatrix) = (size(σ,1) > 1) ? _scaled_pairwisel2(x, y, σ) : pairwisel2(x, y) ./ (σ'.^2)

function _scaled_pairwisel2(x, y, σ)
	# @assert ((size(y, 1) == size(σ, 1) ) && (size(x, dims = 2) == size(σ, dims = 2)) && (size(x, 1) == size(y, 1)))
	@assert ((size(y, 1) == size(σ, 1) ) && (size(x, 2) == size(σ, 2)) && (size(x, 1) == size(y, 1)))
	o = zeros(size(x,2), size(y,2))
	@inbounds for i in 1:size(y,2)
		for j in 1:size(x,2)
			for k in 1:size(x,1)
				o[j,i] += ((x[k,j] - y[k,i])/σ[k,j])^2
			end
		end
	end
	o
end

function _scaled_pairwisel2_m_back(Δ, x, y, σ)
	@assert ((size(y, 1) == size(σ, 1) ) && (size(x, 2) == size(σ, 2)) && (size(x, 1) == size(y, 1)))
	∇x = zero(x)
	∇y = zero(y)
	∇σ = zero(σ)
	@inbounds for i in 1:size(y,2)
		for j in 1:size(x,2)
			for k in 1:size(x,1)
				δ = x[k,j] - y[k,i]
				∇x[k,j] += 2*Δ[j,i]*δ/σ[k,j]^2
				∇y[k,i] -= 2*Δ[j,i]*δ/σ[k,j]^2
				∇σ[k,j] -= 2*Δ[j,i]*δ^2/σ[k,j]^3
			end
		end
	end
	(∇x, ∇y, ∇σ)
end

for x in [AbstractMatrix,TrackedMatrix], y in  [AbstractMatrix,TrackedMatrix], z in [AbstractMatrix,TrackedMatrix]
 x == y == z == AbstractMatrix && continue
 eval(:(scaled_pairwisel2(x::TX, y::TY, σ::TZ) where {TX<:$x, TY<:$y, TZ<:$z} = (size(σ,1) == 1 ) ? pairwisel2(x, y) ./ (σ'.^2) : Flux.Tracker.track(scaled_pairwisel2, x, y, σ)))
end
Flux.Tracker.@grad function scaled_pairwisel2(x, y, σ)
  return(_scaled_pairwisel2(Flux.data(x), Flux.data(y), Flux.data(σ)), Δ -> _scaled_pairwisel2_m_back(Flux.data(Δ), Flux.data(x), Flux.data(y), Flux.data(σ)))
end

"""
	crosspdf_normal(x,c,σ2::T)

	probability density of Normal Distribution of samples in `x` (each column is
	one sample) with respect to a series of Normal Distributions defined
	by centers in `c` (each columns is one center) and standard deviation σ
	if σ is:

		- Number means a σ shared by all centers
		- Vector means a σ is diagonal shared by all centers
		- Transpose means each center has its own scalar σ
		- Matrix means each center has its own scalar diagonal σ

	This should be compatible with Flux
"""
crosspdf_normal(x, c, σ) = exp.(crosslog_normal(x, c, σ))

crosslog_normal(x, c, σ ::T) where {T<:Real} = - 0.5 .* scaled_pairwisel2(c, x, σ) .- size(x,1)*log(2π*σ^2)/2
# crosslog_normal(x, c, σ ::T) where {T<:Union{Transpose, TrackedArray{T, N, A} where {T, N, A<: Transpose}}} = - 0.5 .* scaled_pairwisel2(c, x, σ) .- size(x,1)*log(2π)/2 .- size(x,1)*log.(σ')
crosslog_normal(x, c, σ ::T) where {T<:AbstractVector} = - 0.5 .* scaled_pairwisel2(c, x, σ) .- size(x,1)*log(2π)/2 .- sum(log.(σ))
crosslog_normal(x, c, σ ::T) where {T<:AbstractMatrix} = - 0.5 .* scaled_pairwisel2(c, x, σ) .- size(x,1)*log(2π)/2 .- (1+ size(x,1) - size(σ,1)).*sum(log.(σ), dims=1)'

"""
		kldiv(μ,σ2)

		kl-divergence of a Gaussian min mean `μ` and diagonal variance `σ^2`
		to N(0,I)
"""
kldiv(μ,σ2) = - mean(sum((@.log(σ2) - μ^2 - σ2), dims=1))

"""
		log_normal(x,μ,σ2 = I)

		log-likelihood of x to the Normal with centre at mu
"""
log_normal(x) = - sum(x.^2, dims=1) / 2 .- size(x,1)*log(2π) / 2
log_normal(x,μ) = log_normal(x .- μ)
log_normal(x,μ, σ2::T) where {T<:Number} = - sum((@. ((x - μ)^2)/σ2), dims=1)/2 .- size(x,1)*log(σ2*2π)/2

log_bernoulli(x::AbstractMatrix,θ::AbstractVector) = log.(θ)' * x
log_bernoulli(x::AbstractMatrix,θ::AbstractMatrix) = sum(x .* log.(θ),dims=1)

max_stars_repo_name: pevnak/FluxExtensions

max_issues_repo_name: pevnak/FluxExtensions

max_forks_repo_name: pevnak/FluxExtensions


text:
abstract type GraphicElement3D <: GraphicElement end

#--- 3D POINTS -----------------------------------------
#-------------------------------------------------------
struct Vec3
    x::Real
    y::Real
    z::Real
end

Vec3(t::Tuple{<:Real,<:Real,<:Real}) = Vec3(t...)

string(v::Vec3) = "($(v.x),$(v.y),$(v.z))"

"""
    Point3D(x::Real,y::Real,z::Real; label="",pen=Pen())
    Point3D(P; label="",pen=Pen())

A graphics primitive representing a three-dimensional point.

`P` may be a 3-tuple of real numbers or a `Vec3`

# Examples
```julia-repl
julia> Point3D(0,4,5;pen="DarkGreen")
```
"""
struct Point3D <: GraphicElement3D
    P::Vec3
    label::AbstractString
    pen::Pen
end

const _DEFAULT_POINT3D_KWARGS =
    OrderedDict(
        :label=>"",
        :pen=>Pen()
    )

Point3D(v::Vec3;kwargs...) =
    Point3D(v,updatedvals(_DEFAULT_POINT3D_KWARGS,
                            process_pen_kwargs(kwargs))...)

Point3D(x::Real,y::Real,z::Real;kwargs...) =
                            Point3D(Vec3(x,y,z);kwargs...)

Point(x::Real,y::Real,z::Real;kwargs...) =
                            Point3D(x,y,z;kwargs...)

function AsyString(point::Point3D)
    label = point.label == "" ? "" : "L=$(enclosequote(point.label))"
    AsyString("""
    dot($(filterjoin(label,point.P,point.pen)));
    """)
end

function Base.show(P::Point3D)
    kwargs = kwargstring(P,_DEFAULT_POINT3D_KWARGS)
    print(io,"Point3D($(P.P.x),$(P.P.y),$(P.P.z)$kwargs)")
end

#--- 3D PATHS ------------------------------------------
#-------------------------------------------------------

"""
    Path3D(points;label="",pen=Pen(),arrow="",spline=false)

A graphics primitive representing a two-dimensional path

`points` may be an Array of `Vec3`s or an Array of 3-tuples.
Alternatively, iterables of coordinates may be supplied
separately as `x` and `y`

# Examples
```julia-repl
julia> Path3D([(0,0),(1,0),(1,1)];pen="MidnightBlue")
```
"""
struct Path3D <: GraphicElement3D
    points::Array{Vec3,1}
    label::AbstractString
    pen::Pen
    arrow::Arrow
    spline::Bool
end

const _DEFAULT_PATH3D_KWARGS =
    OrderedDict(
        :label => "",
        :pen => Pen(),
        :arrow => NoArrow(),
        :spline => false
    )

Path3D(points::Array{Vec3,1};kwargs...) =
    Path3D(points,updatedvals(_DEFAULT_PATH3D_KWARGS,
                                    process_pen_kwargs(kwargs))...)

Path3D(points::Array{<:Tuple{<:Real,<:Real,<:Real},1};kwargs...) =
                                Path3D(map(Vec3,points);kwargs...)

Path(points::Array{<:Tuple{<:Real,<:Real,<:Real},1};kwargs...) =
                                Path3D(map(Vec3,points);kwargs...)

Path3D(coords::Array{<:Real,2};kwargs...) =
            Path3D([Vec3(coords[i,:]...) for i=1:size(coords,1)];kwargs...)

function AsyString(P::Path3D)
    label = P.label == "" ? "" : "L=$(enclosequote(P.label))"
    pathname = "p{IDENTIFIER}"
    coords = vcat([[p.x p.y p.z] for p in P.points]...)
    spline = P.spline ? ".." : "--"
    arrow = P.arrow.name == "None" ? "" : "arrow=$(P.arrow)"
    if length(P.points) > 10
        AsyString("""
        file pathdata = input("path{IDENTIFIER}.csv");
        real[][] A = pathdata.csv().dimension(0,3);
        close(pathdata);
        guide3 $pathname;
        for(int i=0; i<A.length; ++i){
            $pathname = $pathname $spline (A[i][0],A[i][1],A[i][2]);
        }
        draw($(filterjoin(label,pathname,P.pen,arrow)));
        """,
        Dict("path{IDENTIFIER}.csv"=>coords))
    else
        AsyString("""
        path3 $pathname = $(join([string(p) for p in P.points],spline));
        draw($(filterjoin(label,pathname,P.pen,arrow)));
        """)
    end
end

function Base.show(io::IO,P::Path3D)
    l = length(P.points)
    s = l ≠ 1 ? "s" : ""
    kwargs = kwargstring(P,_DEFAULT_PATH3D_KWARGS)
    print(io,"""
    Path3D(<$l point$s>$kwargs)
    """)
end

#--- Surfaces ------------------------------------------
#-------------------------------------------------------

"""
    Surface(x::Array{<:Real},
            y::Array{<:Real},
            z::Array{<:Real,2};
            options)
    Surface(z::Array{<:Real},2)

A graphics primitive representing a surface in three dimensions
`x` and `y` may be one- or two-dimensional arrays

The surface passes through the points
    [x[i,j],y[i,j],z[i,j] for i=1:size(z,1),j=1:size(z,2)]

The options are
- `colors`: A vector of color names, for coloring
- `spline`: whether to draw a smooth or piecewise smooth surface
- `surfacepen`: a pen for drawing the surface
- `meshpen`: a pen for drawing the grid lines on the surface
- `clip`: either `false` or a boolean array of the same dimensions
          as `x`, `y`, and `z`, specifying patches to exclude
"""
struct Surface <: GraphicElement3D
    x::Array{<:Real,2}
    y::Array{<:Real,2}
    z::Array{<:Real,2}
    colors::Union{Vector,Array{Pen,2}}
    spline::Bool
    surfacepen::Pen
    meshpen::Pen
    clip # either false or an Array{Bool,2} of the same size as z
end

const _DEFAULT_SURFACE_KWARGS =
    OrderedDict(
         :colors => ["MidnightBlue",
                     "LightSeaGreen",
                     "Yellow",
                     "Tomato"],
         :spline => true,
         :surfacepen => Pen(),
         :meshpen => Pen(),
         :clip => false)

function Surface(x::Array{<:Real,2},
                 y::Array{<:Real,2},
                 z::Array{<:Real,2};kwargs...)
    Surface(x,y,z,updatedvals(_DEFAULT_SURFACE_KWARGS,kwargs)...)
end

function Surface(z::Array{<:Real,2};kwargs...)
    x = [i-1 for i=1:size(z,1),j=1:size(z,2)]
    y = [j-1 for i=1:size(z,1),j=1:size(z,2)]
    Surface(x,y,z;kwargs...)
end

function Surface(x::Union{Array{<:Real,1},AbstractRange},
                 y::Union{Array{<:Real,1},AbstractRange},
                 z::Array{<:Real,2};kwargs...)
    xtwo = [xi for xi=x,yi=y]
    ytwo = [yi for xi=x,yi=y]
    Surface(xtwo,ytwo,z;kwargs...)
end

function AsyString(S::Surface)
    r(P::Pen) = Float64(P.color.color.r)
    g(P::Pen) = Float64(P.color.color.g)
    b(P::Pen) = Float64(P.color.color.b)
    o(P::Pen) = Float64(P.opacity)
    asyData = Dict()
    for w in (:x,:y,:z)
        asyData[string(w)*"{IDENTIFIER}.csv"] = getfield(S,w)
    end
    D = Dict(k => getfield(S,k) for k in keys(_DEFAULT_SURFACE_KWARGS))
    spline = D[:spline] ? ",usplinetype=Spline,vsplinetype=Spline" : ""
    surfacepen = isdefault(D[:surfacepen]) ? "" : "surfacepen=$(D[:surfacepen])"
    meshpen = isdefault(D[:meshpen]) ? "" : "meshpen=$(D[:meshpen])"

    surfacepenmod = """
    p[i] += $(D[:surfacepen]);
    """
    if isdefault(D[:surfacepen])
        surfacepenmod = ""
    end

    safepaste(s) = length(s) > 1000 ? error("Misspecified colors") : s

    if isa(D[:colors],Array{Pen,2})
        merge!(asyData,Dict(
            "red{IDENTIFIER}.csv"=>r.(D[:colors]),
            "green{IDENTIFIER}.csv"=>g.(D[:colors]),
            "blue{IDENTIFIER}.csv"=>b.(D[:colors]),
            "alpha{IDENTIFIER}.csv"=>o.(D[:colors])))
        colorasy = """
        int m = $(size(D[:colors],1));
        int n = $(size(D[:colors],2));
        file reds = input("red{IDENTIFIER}.csv");
        real[][] redvalues = reds.csv().dimension(m,n);
        close(reds);
        file blues = input("blue{IDENTIFIER}.csv");
        real[][] bluevalues = blues.csv().dimension(m,n);
        close(blues);
        file greens = input("green{IDENTIFIER}.csv");
        real[][] greenvalues = greens.csv().dimension(m,n);
        close(greens);
        file alphas = input("alpha{IDENTIFIER}.csv");
        real[][] alphavalues = alphas.csv().dimension(m,n);
        close(alphas);

        pen[][] pixels = new pen[m][n];
        for(int i=0;i<m;++i){
          for(int j=0;j<n;++j){
              pixels[i][j] = (rgb(redvalues[i][j],
                                  greenvalues[i][j],
                                  bluevalues[i][j]) +
                                  opacity(alphavalues[i][j]));
          }
        }

        pen[][] patchpens = new pen[(m-1)*(n-1)][4];

        int splineoffset = $(D[:spline] ? 1 : 0);

        for(int i=0; i<m-1; ++i){
          for(int j=0; j<n-1; ++j){
            patchpens[(n-1)*i+j][0] = pixels[i][j];
            patchpens[(n-1)*i+j][1] = pixels[i+splineoffset][j];
            patchpens[(n-1)*i+j][2] = pixels[i+splineoffset][j+splineoffset];
            patchpens[(n-1)*i+j][3] = pixels[i][j+splineoffset];
          }
        }

        s.colors(patchpens);

        """
    else
        surfacepenmodloop = surfacepenmod == "" ? "" :
            """
            for(int i=0; i<p.length-1; ++i) {
                p[i] += $surfacepenmod;
            }
            """
        colorasy = """
        pen[] p = {$(safepaste((join(D[:colors],','))))};

        $surfacepenmodloop

        s.colors(palette(s.map(zpart),Gradient(100 ... p)));
        """
    end

    if D[:clip] ≠ false
        if (!isa(D[:clip],AbstractArray{<:Any,2})
            || size(D[:clip]) ≠ size(S.z))
            error("clip misspecified")
        end
        clipfile = """
        file condfile = input("cond{IDENTIFIER}.csv");
        int[][] cond = condfile.csv().dimension(m,n);
        close(condfile);
        """
        clipfunction = """
        bool cond(pair t) {
            int i = round(t.x);
            int j = round(t.y);
            return cond[i][j] == 0 ? false : true;
        }
        """
        cliparg = ",cond=cond"
        asyData["cond{IDENTIFIER}.csv"] =
                    convert(Array{Int64,2},D[:clip])
    else
        clipfile = ""
        clipfunction = ""
        cliparg = ""
    end

    AsyString("""
    int m = $(size(S.z,1));
    int n = $(size(S.z,2));

    file xdata = input("x{IDENTIFIER}.csv");
    real[][] x = xdata.csv().dimension(m,n);
    close(xdata);

    file ydata = input("y{IDENTIFIER}.csv");
    real[][] y = ydata.csv().dimension(m,n);
    close(ydata);

    file zdata = input("z{IDENTIFIER}.csv");
    real[][] z = zdata.csv().dimension(m,n);
    close(zdata);

    $clipfile

    triple f(pair t) {
      int i=round(t.x);
      int j=round(t.y);
      return (x[i][j],y[i][j],z[i][j]);
    }

    $clipfunction

    surface s = surface(f,(0,0),(m-1,n-1),
                            nu=m-1,nv=n-1$spline$cliparg);

    $colorasy

    draw($(filterjoin("s",surfacepen,meshpen,
                "render(merge=true)","light=nolight")));""",
    asyData)
end

function Base.show(io::IO,S::Surface)
    kwargs = kwargstring(S,_DEFAULT_SURFACE_KWARGS)
    print(io,"Surface(<$(size(S.z,1))×$(size(S.z,2))>$kwargs)")
end

#--- 3D POLYGONS -----------------------------------------
#---------------------------------------------------------

"""
    Polygon3D(points;pen=Pen(),
                     fillpen=Pen(color="white"),
                     spline=false)

A graphics primitive representing a three-dimensional polygon

`points` may be an Array of `Vec3`s or an Array of 3-tuples.

# Examples
```julia-repl
julia> Polygon3D([(0,0,0),(1,0,0),(1,1,0)];pen="MidnightBlue")
```
"""
struct Polygon3D <: GraphicElement3D
    points::Array{Vec3,1}
    pen::Pen
    fillpen::Pen
    spline::Bool
end

const _DEFAULT_POLYGON3D_KWARGS =
    OrderedDict(
        :pen => Pen(),
        :fillpen => Pen(color="white"),
        :spline => false
    )

Polygon3D(points::Array{Vec3,1};kwargs...) =
    Polygon3D(points,updatedvals(_DEFAULT_POLYGON3D_KWARGS,kwargs)...)

Polygon(points::Array{Vec3,1};kwargs...) = Polygon3D(points;kwargs...)

Polygon3D(points::Array{<:Tuple{<:Real,<:Real,<:Real},1};kwargs...) =
                                Polygon3D(map(Vec3,points);kwargs...)

Polygon3D(coords::Array{<:Real,2};kwargs...) =
            Polygon3D([Vec3(coords[i,:]...) for i=1:size(coords,1)];kwargs...)

Polygon3D(points::Array{<:RealOrComplex,1};kwargs...) =
                        Polygon3D(map(Vec3,points);kwargs...)

function Polygon(P::Path3D;kwargs...)
    lastindex = endof(P.points) - (P.points[1] == P.points[end] ? 1 : 0)
    Polygon(P.points[1:lastindex];kwargs...)
end

function AsyString(P::Polygon3D)
    pathname = "p{IDENTIFIER}"
    coords = vcat([[p.x p.y p.z] for p in P.points]...)
    spline = P.spline ? ".." : "--"
    surface = "surface($pathname)"
    if length(P.points) > 10
        AsyString("""
        file pathdata = input("path{IDENTIFIER}.csv");
        real[][] A = pathdata.csv().dimension(0,3);
        close(pathdata);
        guide3 $pathname;
        for(int i=0; i<A.length; ++i){
            $pathname = $pathname $spline (A[i][0],A[i][1],A[i][2]);
        }
        $pathname = $pathname $spline cycle;
        draw($(filterjoin(surface,P.fillpen,P.pen)));
        """,
        Dict("path{IDENTIFIER}.csv"=>coords))
    else
        AsyString("""
        path3 $pathname = $(join([string(p) for p in P.points],spline))--cycle;
        draw($(filterjoin(surface,P.fillpen,P.pen)));
        """)
    end
end

function Base.show(io::IO,P::Polygon3D)
    l = length(P.points)
    s = l ≠ 1 ? "s" : ""
    kwargs = kwargstring(P,_DEFAULT_POLYGON3D_KWARGS)
    print(io,"""
    Polygon3D(<$l point$s>$kwargs)
    """)
end

#--- 3D LABELS ----------------------------------------
#------------------------------------------------------

struct Label3D <: GraphicElement3D
    s::AbstractString
    location::Vec3
    pen::Pen
end

const _DEFAULT_LABEL3D_KWARGS =
    OrderedDict(
        :pen => Pen()
    )

Label3D(s::AbstractString,v::Vec3;kwargs...) =
    Label3D(s,v,updatedvals(_DEFAULT_LABEL3D_KWARGS,
                                process_pen_kwargs(kwargs))...)
Label(s::AbstractString,v::Vec3;kwargs...) = Label3D(s,v;kwargs...)
Label3D(s::AbstractString,
        P::Tuple{<:Real,<:Real,<:Real};kwargs...) =
                                        Label3D(s,Vec3(P);kwargs...)
Label(s::AbstractString,P::Tuple{<:Real,<:Real,<:Real};kwargs...) =
                                             Label3D(s,P;kwargs...)
Label3D(s,x::Real,y::Real;kwargs...) =
                Label3D(s,Vec3(x,y);kwargs...)
Label(s::AbstractString,x::Real,y::Real,z::Real;kwargs...) = Label3D(s,x,y;kwargs...)

function AsyString(L::Union{Label2D,Label3D})
    pen = isdefault(L.pen) ? "" : ",p=$(string(L.pen))"
    rotation = (isa(L,Label3D) || L.rotation == 0.0) ? "" : "rotate($(L.rotation))*"
    AsyString("label($rotation$(enclosequote(L.s)),$(string(L.location))$pen);")
end

function Base.show(io::IO,L::Label3D)
    kwargs = kwargstring(L,_DEFAULT_LABEL3D_KWARGS)
    P = L.location
    print(io,"Label3D($(enclosequote(L.s)),($(P.x),$(P.y),$(P.z))$kwargs)")
end

#--- 3D Plots ------------------------------------------
#-------------------------------------------------------

"""
    Plot3D(elements::Array{<:GraphicElement2D,1},
           options::Array{Any,1})

A container for a list of graphics primitives to draw,
together with drawing options

# Examples
```julia-repl
julia> Plot(Path([0 0 0; 1 1 1]),Path([0 0 0; 0 0 1]))
```
"""
struct Plot3D <: Plot
    elements::Array{<:GraphicElement,1}
    options::Array{Any,1}
end

+(G::GraphicElement, H::GraphicElement) = Plot([G,H])
+(P::Plot,G::GraphicElement) = P + Plot(G)
+(G::GraphicElement,P) = P + G

Plot3D() = Plot3D(GraphicElement[],[])

function Plot3D(elements::Array{<:GraphicElement};kwargs...)
    Plot3D(elements,collect(kwargs))
end

Plot3D(element::GraphicElement;kwargs...) = Plot3D([element];kwargs...)

function Plot(elements...;kwargs...)
    if length(elements) == 0
        Plot(GraphicElement[])
    end
    e = elements[1]
    if length(elements) > 1 || isa(e,GraphicElement)
        return Plot(collect(elements);kwargs...)
    end
    if isa(e,Vector{<:GraphicElement})
        if all(isa(g,Union{GraphicElement2D,RawString}) for g in e)
            return Plot2D(e;kwargs...)
        else
            return Plot3D(e;kwargs...)
        end
    else
        error("Argument to Plot should be a Vector{<:GraphicElement}")
    end
end

const _DEFAULT_PLOT3D_KWARGS =
    OrderedDict(
         :axes => true,
         :pdf => false,
         :render => 4,
         :xlabel => "",
         :ylabel => "",
         :zlabel => "",
         :packages => ["mathpazo"],
         :xmin => "-infinity",
         :xmax => "infinity",
         :ymin => "-infinity",
         :ymax => "infinity",
         :zmin => "-infinity",
         :zmax => "infinity",
         :xticks => "NoTicks3",
         :yticks => "NoTicks3",
         :zticks => "NoTicks3",
         :arrow => Arrow3(),
         :camera => nothing,
         :up => nothing,
         :projection => "perspective",
         :bgcolor => NamedColor("white"),
         :width => _DEFAULT_WIDTH,
         :ignoreaspect => false)

function AsyString(P::Plot3D)
    global _DEFAULT_PLOT3D_KWARGS
    D = copy(_DEFAULT_PLOT3D_KWARGS)
    for k in keys(Dict(P.options))
        if !(k in keys(D))
            error("Unknown Plot3D option $k")
        end
    end
    merge!(D,Dict(P.options))
    for s in [:xmin,:xmax,:ymin,:ymax,:zmin,:zmax]
        if isa(D[s],Real)
            D[s] = string(D[s])
        end
    end
    if isa(D[:bgcolor],AbstractString)
        D[:bgcolor] = NamedColor(D[:bgcolor])
    end
    if D[:axes]
        axesstring = """
        xaxis3(  L   = $(enclosequote(D[:xlabel])),
               arrow = $(D[:arrow]),
               xmin  = $(D[:xmin]),
               xmax  = $(D[:xmax]),
               ticks = $(D[:xticks]));
        yaxis3(  L   = $(enclosequote(D[:ylabel])),
               arrow = $(D[:arrow]),
               ymin  = $(D[:ymin]),
               ymax  = $(D[:ymax]),
               ticks = $(D[:yticks]));
        zaxis3(  L   = $(enclosequote(D[:zlabel])),
               arrow = $(D[:arrow]),
               zmin  = $(D[:zmin]),
               zmax  = $(D[:zmax]),
               ticks = $(D[:zticks]));
        """
    else
        axesstring = ""
    end

    if occursin("(",D[:projection])
        currentprojection = D[:projection]
    elseif D[:camera] != nothing
        a,b,c = D[:camera]
        up = isa(D[:up],Nothing) ? "" : "up=$(D[:up])"
        currentprojection = "$(D[:projection])($(filterjoin(a,b,c,up)));"
    else
        currentprojection = "$(D[:projection])(M.x+(M.x-m.x),M.y+0.25*(M.y-m.y),M.z+0.5*(M.z-m.z));"
    end

    bgcolor = "currentlight.background = $(string(D[:bgcolor]));"

    ignoreaspect = D[:ignoreaspect] ? ",IgnoreAspect" : ""
    pdfstring = """
    settings.outformat = "pdf";
    settings.render = $(D[:render]);
    """
    pdf = D[:pdf] ? pdfstring : ""
    packages = join("usepackage($(enclosequote(s)));"
                                    for s in D[:packages])
    shipout = D[:pdf] ? "shipout(bbox(3.0,invisible));" : ""

    asystrings = map(AsyString,P.elements)
    drawingcommands = join(replace(s.str,"{IDENTIFIER}"=>string(j))
                                for (j,s) in enumerate(asystrings))

    data = [A.data for A in asystrings]
    merged_data = merge([Dict(replace(k,"{IDENTIFIER}"=>j) => v
                for (k,v) in D) for (j,D) in enumerate(data)]...)

    AsyString("""
    import graph3;
    import palette;
    import x11colors;

    $packages

    $pdf

    size3($(D[:width])$ignoreaspect);

    $bgcolor

    $drawingcommands

    $axesstring

    triple M = currentpicture.userMax();
    triple m = currentpicture.userMin();
    currentprojection = $currentprojection
    $shipout
    """,
    merged_data)
end

function +(P::Plot3D,Q::Plot3D)
    global _DEFAULT_PLOT2D_KWARGS
    elements = vcat(P.elements,Q.elements)
    options = copy(_DEFAULT_PLOT3D_KWARGS)
    merge!(options,Dict(P.options))
    merge!(options,Dict(Q.options))
    return Plot3D(elements,collect(options))
end

#--- OUTPUT --------------------------------------------
#-------------------------------------------------------

function plot(z::Array{<:Real,2};kwargs...)
    global _DEFAULT_SURFACE_KWARGS
    surfkwargs, plotkwargs = splitkwargs(kwargs,_DEFAULT_SURFACE_KWARGS)
    Plot3D(Surface(z;surfkwargs...);ignoreaspect=true,plotkwargs...)
end

function plot(x,y,z::Array{<:Real,2};kwargs...)
    global _DEFAULT_SURFACE_KWARGS
    surfkwargs, plotkwargs = splitkwargs(kwargs,_DEFAULT_SURFACE_KWARGS)
    Plot3D(Surface(x,y,z;surfkwargs...);ignoreaspect=true,plotkwargs...)
end

function plot(f::Function,
              a::Real,
              b::Real,
              c::Real,
              d::Real;
              m::Integer=100,
              n::Integer=100,
              kwargs...)
    x = range(a,stop=b,length=m)
    y = range(c,stop=d,length=n)
    plot(x,y,[f(p,q) for p=x,q=y]; kwargs...)
end

plot(f::Function,x::Tuple{<:Real,<:Real},y::Tuple{<:Real,<:Real};kwargs...) =
    plot(f,x...,y...;kwargs...)

plot(xs,ys,f::Function; kwargs...) = plot(xs,ys,[f(x,y) for x in xs, y in ys]; kwargs...)

function interactive(P::Plot3D;kwargs...)
    check_asy_present()
    Q = deepcopy(P)
    filter!(p -> first(p) ≠ :pdf, Q.options)
    push!(Q.options,(:pdf,false))
    for (k,v) in kwargs
        filter!(p -> first(p) ≠ k, Q.options)
        push!(Q.options,(k,v))
    end
    directory = mktempdir()
    filename = "$directory/myplot.asy"
    save(filename,Q;runasy=false)
    originaldir = pwd()
    cd(directory)
    run(`asy -V myplot.asy`)
    cd(originaldir)
end

interactive(S::GraphicElement3D;kwargs...) =
        interactive(Plot3D([S],kwargs))

interactive(z::Array{<:Real,2};kwargs...) = interactive(plot(z;kwargs...))


text:
import Distributions
using Random
using SpecialFunctions: besseli

struct CPNormal{T,S} <: Distributions.ContinuousUnivariateDistribution
    μ::T
    σ::S
end

Base.minimum(::CPNormal) = -Inf
Base.maximum(::CPNormal) = Inf

const log2π = log(2*π)
function Distributions.logpdf(dist::CPNormal, x::Real)
    μ,σ = dist.μ, dist.σ
    s,c = sincos(x-μ)
    dθ = atan(s, c)
    return  -(abs2(dθ/σ) + log2π)/2 - log(σ)
end

function Base.rand(rng::AbstractRNG, d::CPNormal)
    return d.μ + d.σ*randn(rng)
end

Base.rand(rng::AbstractRNG, ::Type{Float64}, d::CPNormal) = rand(rng, d)